function  [tau,N,B,CP,Q,R,tau2,N2,H,MR,VR]=absorb(P);
% format: [tau,N,B,CP,Q,R,tau2,N2,H,MR,VR]=absorb(P); 
% input: P an absorbing Markov transition matrix in 'from rows
%        to column' form.  Rows are probability vectors and
%        must sum to 1.0.  At least one main diagonal element=1
% output: tau: Time before absorption from a transient state.
%        B: P(ending up in absorbing state from a transient state)
%        tau2: variance matrix for tau;N: fundamental matrix
%        N2: covariance matrix for N; 
%        CP: transition matrix in canonical form.
%        Q: Submatrix of transient states, R: from transient to 
%        absorbing state submatrix,
%        H: hij the probability process will ever go to transient 
%           state j starting in transient state i (not counting the
%           initial state (K & S, p. 61))
%        MR: expected number of times that a Markov process remains
%            in a transient state once entered (including entering step)
%        VR: Variance for MR (K & S, Theorem 3.5.6, p. 61)
% written by E. Gallagher, ECOS
% UMASS/Boston email: Eugene.Gallagher@umb.edu
% written 9/26/93, revised: 3/26/09
% refs: 
% Kemeny, J. G. and J. L. Snell.  1976.  Finite Markov chains.  
%        Springer-Verlag, New York, New York, U.S.A. 
% Roberts, F. 1976.  Discrete Mathematical Models with applications
%        to social, biological, and environmental problems. Prentice-Hall
%
[pdim,pdim]=size(P);
if sum([sum(P')-ones(1,pdim)].^2)>0.001
      error('Rows of P must sum to 1.0')
end
dP=diag(P);
ri=find(dP==1.0);
if isempty(ri)
       error('No absorbing states (reenter P or use ergodic.m)')
end
rdim=length(ri);
qi=find(dP~=1.0);qdim=length(qi);
I=eye(qdim,qdim);
Q=P(qi,qi);
N=inv(I-Q);   % the fundamental matrix
tau=sum(N')';
CP=P([ri' qi'],[ri' qi']);  % the canonical form of P
R=CP(rdim+1:pdim,1:rdim);
B=N*R;
if nargout>6          % only if variances requested 
     Ndg=diag(diag(N));
     N2=N*(2*Ndg-I)-N.^2;
     tau2=(2*N-I)*tau-tau.^2;
     H=(N-I)/Ndg;
     dQ=diag(Q);oneq=ones(qdim,1);
     MR=oneq./(oneq-dQ);
     VR=dQ./(oneq-dQ).^2;
end