c Triangular Solve dtrsl() c ---------------- c solves systems of the form c c t * x = b c or c trans(t) * x = b c c where t is a triangular matrix of order n. here trans(t) c denotes the transpose of the matrix t. c c on entry c c t double precision(ldt,n) c t contains the matrix of the system. the zero c elements of the matrix are not referenced, and c the corresponding elements of the array can be c used to store other information. c c ldt integer c ldt is the leading dimension of the array t. c c n integer c n is the order of the system. c c b double precision(n). c b contains the right hand side of the system. c c job integer c job specifies what kind of system is to be solved. c if job is c c 00 solve t*x=b, t lower triangular, c 01 solve t*x=b, t upper triangular, c 10 solve trans(t)*x=b, t lower triangular, c 11 solve trans(t)*x=b, t upper triangular. c c on return c c b b contains the solution, if info .eq. 0. c otherwise b is unaltered. c c info integer c info contains zero if the system is nonsingular. c otherwise info contains the index of c the first zero diagonal element of t. c c linpack. this version dated 08/14/78 . c g. w. stewart, university of maryland, argonne national lab. c c subroutines and functions c c blas: daxpy,ddot c fortran mod c subroutine dtrsl(t,ldt,n,b,job,info) integer ldt,n,job,info double precision t(ldt,1),b(1) c c internal variables c double precision ddot,temp integer case,j,jj c c begin block permitting ...exits to 150 c c check for zero diagonal elements. c do 10 info = 1, n if (t(info,info) .eq. 0.0d0) go to 150 c ......exit 10 continue info = 0 c c determine the task and go to it. c case = 1 if (mod(job,10) .ne. 0) case = 2 if (mod(job,100)/10 .ne. 0) case = case + 2 go to (20,50,80,110), case c C Case 1 (job = 00): c solve t*x=b for t lower triangular c 20 continue b(1) = b(1)/t(1,1) if (n .ge. 2) then do 30 j = 2, n temp = -b(j-1) call daxpy(n-j+1,temp,t(j,j-1),1,b(j),1) b(j) = b(j)/t(j,j) 30 continue endif go to 140 c C Case 2 (job = 01): c solve t*x=b for t upper triangular. c 50 continue b(n) = b(n)/t(n,n) if (n .ge. 2) then do 60 jj = 2, n j = n - jj + 1 temp = -b(j+1) call daxpy(j,temp,t(1,j+1),1,b(1),1) b(j) = b(j)/t(j,j) 60 continue endif go to 140 c C Case 3 (job = 10): c solve trans(t)*x=b for t lower triangular. c 80 continue b(n) = b(n)/t(n,n) if (n .ge. 2) then do 90 jj = 2, n j = n - jj + 1 b(j) = b(j) - ddot(jj-1,t(j+1,j),1,b(j+1),1) b(j) = b(j)/t(j,j) 90 continue endif go to 140 c C Case 4 (job = 11): c solve trans(t)*x=b for t upper triangular. c 110 continue b(1) = b(1)/t(1,1) if (n .ge. 2) then do 120 j = 2, n b(j) = b(j) - ddot(j-1,t(1,j),1,b(1),1) b(j) = b(j)/t(j,j) 120 continue endif C 140 continue c EXIT: 150 continue return end