SUBROUTINE SGESL( A, LDA, N, IPVT, B, JOB ) 2,4 c Solves the real system c A * X = B or transpose(A) * X = B c using the factors computed by SGECO or SGEFA. c Revision date: 8/1/82 c Author: Moler, C. B. (U. of New Mexico) c On entry c A REAL(LDA, N) c the output from SGECO or SGEFA. c LDA INTEGER c the leading dimension of the array A c N INTEGER c the order of the matrix A c IPVT INTEGER(N) c the pivot vector from SGECO or SGEFA. c B REAL(N) c the right hand side vector. c JOB INTEGER c = 0 to solve A*X = B , c = nonzero to solve transpose(A)*X = B c On return c B the solution vector X c Error condition c A division by zero will occur if the input factor contains a c zero on the diagonal. Technically, this indicates singularity, c but it is often caused by improper arguments or improper c setting of LDA. It will not occur if the subroutines are c called correctly and if SGECO has set RCOND .GT. 0.0 c or SGEFA has set INFO .EQ. 0 . c To compute inverse(a) * c where c is a matrix c with p columns c call sgeco(a,lda,n,ipvt,rcond,z) c if (rcond is too small) go to ... c do 10 j = 1, p c call sgesl(a,lda,n,ipvt,c(1,j),0) c 10 continue c --------------------------------------------------------------------- c .. Scalar Arguments .. INTEGER JOB, LDA, N c .. c .. Array Arguments .. INTEGER IPVT( * ) REAL A( LDA, * ), B( * ) c .. c .. Local Scalars .. INTEGER K, KB, L, NM1 REAL T c .. c .. External Functions .. REAL SDOT EXTERNAL SDOT c .. c .. External Subroutines .. EXTERNAL SAXPY c .. NM1 = N - 1 IF( JOB.EQ.0 ) THEN c ** solve A * X = B c ** first solve L*Y = B DO 10 K = 1, NM1 L = IPVT( K ) T = B( L ) IF( L.NE.K ) THEN B( L ) = B( K ) B( K ) = T END IF CALL SAXPY( N - K, T, A( K+1, K ), 1, B( K+1 ), 1 ) 10 CONTINUE c ** now solve U*X = Y DO 20 KB = 1, N K = N + 1 - KB B( K ) = B( K ) / A( K, K ) T = - B( K ) CALL SAXPY( K-1, T, A( 1, K ), 1, B(1), 1 ) 20 CONTINUE ELSE c ** solve trans(A) * X = B c ** first solve trans(U)*Y = B DO 30 K = 1, N T = SDOT( K - 1, A( 1,K ), 1, B( 1 ), 1 ) B( K ) = ( B( K ) - T ) / A( K, K ) 30 CONTINUE c ** now solve trans(l)*x = y DO 40 KB = 1, NM1 K = N - KB B( K ) = B( K ) + SDOT( N - K, A( K+1, K ), 1, B( K+1 ), 1) L = IPVT( K ) IF( L.NE.K ) THEN T = B( L ) B( L ) = B( K ) B( K ) = T END IF 40 CONTINUE END IF END