[t,ldt,n,det,job,info]=ctrdi(t,ldt,n,det,job,info);

function [t,ldt,n,det,job,info]=ctrdi(t,ldt,n,det,job,info);
%***BEGIN PROLOGUE  CTRDI
%***PURPOSE  Compute the determinant and inverse of a triangular matrix.
%***LIBRARY   SLATEC (LINPACK)
%***CATEGORY  D2C3, D3C3
%***TYPE      COMPLEX (STRDI-S, DTRDI-D, CTRDI-C)
%***KEYWORDS  DETERMINANT, INVERSE, LINEAR ALGEBRA, LINPACK,
%             TRIANGULAR MATRIX
%***AUTHOR  Moler, C. B., (U. of New Mexico)
%***DESCRIPTION
%
%     CTRDI computes the determinant and inverse of a complex
%     triangular matrix.
%
%     On Entry
%
%        T       COMPLEX(LDT,N)
%                T contains the triangular matrix.  The zero
%                elements of the matrix are not referenced, and
%                the corresponding elements of the array can be
%                used to store other information.
%
%        LDT     INTEGER
%                LDT is the leading dimension of the array T.
%
%        N       INTEGER
%                N is the order of the system.
%
%        JOB     INTEGER
%                = 010       no det, inverse of lower triangular.
%                = 011       no det, inverse of upper triangular.
%                = 100       det, no inverse.
%                = 110       det, inverse of lower triangular.
%                = 111       det, inverse of upper triangular.
%
%     On Return
%
%        T       inverse of original matrix if requested.
%                Otherwise unchanged.
%
%        DET     COMPLEX(2)
%                determinant of original matrix if requested.
%                Otherwise not referenced.
%                Determinant = DET(1) * 10.0**DET(2)
%                with  1.0 .LE. CABS1(DET(1)) .LT. 10.0
%                or  DET(1) .EQ. 0.0 .
%
%        INFO    INTEGER
%                INFO contains zero if the system is nonsingular
%                and the inverse is requested.
%                Otherwise INFO contains the index of
%                a zero diagonal element of T.
%
%
%***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
%                 Stewart, LINPACK Users' Guide, SIAM, 1979.
%***ROUTINES CALLED  CAXPY, CSCAL
%***REVISION HISTORY  (YYMMDD)
%   780814  DATE WRITTEN
%   890831  Modified array declarations.  (WRB)
%   890831  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900326  Removed duplicate information from DESCRIPTION section.
%           (WRB)
%   920501  Reformatted the REFERENCES section.  (WRB)
%***end PROLOGUE  CTRDI
persistent i j k kb km1 kp1 temp ten zdum ; 

t_shape=size(t);t=reshape([t(:).',zeros(1,ceil(numel(t)./prod([ldt])).*prod([ldt])-numel(t))],ldt,[]);
%
if isempty(temp), temp=0; end;
if isempty(ten), ten=0; end;
if isempty(i), i=0; end;
if isempty(j), j=0; end;
if isempty(k), k=0; end;
if isempty(kb), kb=0; end;
if isempty(km1), km1=0; end;
if isempty(kp1), kp1=0; end;
if isempty(zdum), zdum=0; end;
% cabs1= @(zdum)  abs(real(zdum)) + abs(aimag(zdum));real :: cabs1;
cabs1= @(zdum)  abs(real(zdum)) + abs(imag(zdum));
%***FIRST EXECUTABLE STATEMENT  CTRDI
%
%        COMPUTE DETERMINANT
%
if( fix(job./100)~=0 )
det(1) =complex(1.0e0,0.0e0);
det(2) =complex(0.0e0,0.0e0);
ten = 10.0e0;
for i = 1 : n;
det(1) = t(i,i).*det(1);
cabs1= @(zdum)  abs(real(zdum)) + abs(imag(zdum));
if( cabs1(det(1))==0.0e0 )
break;
end;
cabs1= @(zdum)  abs(real(zdum)) + abs(imag(zdum));
while( cabs1(det(1))<1.0e0 );
det(1) = complex(ten,0.0e0).*det(1);
det(2) = det(2) -complex(1.0e0,0.0e0);
end;
cabs1= @(zdum)  abs(real(zdum)) + abs(imag(zdum));
while( cabs1(det(1))>=ten );
det(1) = det(1)./complex(ten,0.0e0);
det(2) = det(2) +complex(1.0e0,0.0e0);
end;
end;
end;
%
%        COMPUTE INVERSE OF UPPER TRIANGULAR
%
if( rem(fix(job./10),10)~=0 )
if( rem(job,10)==0 )
%
%              COMPUTE INVERSE OF LOWER TRIANGULAR
%
for kb = 1 : n;
k = fix(n + 1 - kb);
info = fix(k);
cabs1= @(zdum)  abs(real(zdum)) + abs(imag(zdum));
if( cabs1(t(k,k))==0.0e0 )
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
return;
end;
t(k,k) =complex(1.0e0,0.0e0)./t(k,k);
temp = -t(k,k);
if( k~=n )
[dumvar1,temp,t(sub2ind(size(t),k+1,k):end)]=cscal(n-k,temp,t(sub2ind(size(t),k+1,k):end),1);
end;
km1 = fix(k - 1);
if( km1>=1 )
for j = 1 : km1;
temp = t(k,j);
t(k,j) =complex(0.0e0,0.0e0);
[dumvar1,temp,dumvar3,dumvar4,dumvar5]=caxpy(n-k+1,temp,t(sub2ind(size(t),k,k):end),1,t(sub2ind(size(t),k,j):end),1);      t(sub2ind(size(t),k,k):end)=dumvar3; t(sub2ind(size(t),k,j):end)=dumvar5; 
end; j = fix(km1+1);
end;
end; kb = fix(n+1);
info = 0;
else;
for k = 1 : n;
info = fix(k);
cabs1= @(zdum)  abs(real(zdum)) + abs(imag(zdum));
if( cabs1(t(k,k))==0.0e0 )
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
return;
end;
t(k,k) =complex(1.0e0,0.0e0)./t(k,k);
temp = -t(k,k);
[dumvar1,temp,t(sub2ind(size(t),1,k):end)]=cscal(k-1,temp,t(sub2ind(size(t),1,k):end),1);
kp1 = fix(k + 1);
if( n>=kp1 )
for j = kp1 : n;
temp = t(k,j);
t(k,j) =complex(0.0e0,0.0e0);
[k,temp,dumvar3,dumvar4,dumvar5]=caxpy(k,temp,t(sub2ind(size(t),1,k):end),1,t(sub2ind(size(t),1,j):end),1);      t(sub2ind(size(t),1,k):end)=dumvar3; t(sub2ind(size(t),1,j):end)=dumvar5; 
end; j = fix(n+1);
end;
end; k = fix(n+1);
info = 0;
end;
end;
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
end %subroutine ctrdi
%DECK CTRMM