Subroutine Dgefa(A,Lda,N,Ipvt,Info) Integer Lda,N,Ipvt(*),Info Double Precision A(Lda,*) C C Dgefa Factors A Double Precision Matrix By Gaussian Elimination. C C Dgefa Is Usually Called By Dgeco, But It Can Be Called C Directly With A Saving In Time If Rcond Is Not Needed. C (Time For Dgeco) = (1 + 9/N)*(Time For Dgefa) . C C On Entry C C A Double Precision(Lda, N) C The Matrix To Be Factored. C C Lda Integer C The Leading Dimension Of The Array A . C C N Integer C The Order Of The Matrix A . C C On Return C C A An Upper Triangular Matrix And The Multipliers C Which Were Used To Obtain It. C The Factorization Can Be Written A = L*U Where C L Is A Product Of Permutation And Unit Lower C Triangular Matrices And U Is Upper Triangular. C C Ipvt Integer(N) C An Integer Vector Of Pivot Indices. C C Info Integer C = 0 Normal Value. C = K If U(K,K) .eq. 0.0 . This Is Not An Error C Condition For This Subroutine, But It Does C Indicate That Dgesl Or Dgedi Will Divide By Zero C If Called. Use Rcond In Dgeco For A Reliable C Indication Of Singularity. C C Linpack. This Version Dated 08/14/78 . C Cleve Moler, University Of New Mexico, Argonne National Lab. C C Subroutines And Functions C C Blas Daxpy,Dscal,Idamax C C Internal Variables C Double Precision T Integer Idamax,J,K,Kp1,L,Nm1 C C C Gaussian Elimination With Partial Pivoting C Info = 0 Nm1 = N - 1 If (Nm1 .lt. 1) Go To 70 Do 60 K = 1, Nm1 Kp1 = K + 1 C C Find L = Pivot Index C L = Idamax(N-K+1,A(K,K),1) + K - 1 Ipvt(K) = L C C Zero Pivot Implies This Column Already Triangularized C If (A(L,K) .eq. 0.0d0) Go To 40 C C Interchange If Necessary C If (L .eq. K) Go To 10 T = A(L,K) A(L,K) = A(K,K) A(K,K) = T 10 Continue C C Compute Multipliers C T = -1.0d0/A(K,K) Call Dscal(N-K,T,A(K+1,K),1) C C Row Elimination With Column Indexing C Do 30 J = Kp1, N T = A(L,J) If (L .eq. K) Go To 20 A(L,J) = A(K,J) A(K,J) = T 20 Continue Call Daxpy(N-K,T,A(K+1,K),1,A(K+1,J),1) 30 Continue Go To 50 40 Continue Info = K 50 Continue 60 Continue 70 Continue Ipvt(N) = N If (A(N,N) .eq. 0.0d0) Info = N Return End Subroutine Dgesl(A,Lda,N,Ipvt,B,Job) Integer Lda,N,Ipvt(*),Job Double Precision A(Lda,*),B(*) C C Dgesl Solves The Double Precision System C A * X = B Or Trans(A) * X = B C Using The Factors Computed By Dgeco Or Dgefa. C C On Entry C C A Double Precision(Lda, N) C The Output From Dgeco Or Dgefa. C C Lda Integer C The Leading Dimension Of The Array A . C C N Integer C The Order Of The Matrix A . C C Ipvt Integer(N) C The Pivot Vector From Dgeco Or Dgefa. C C B Double Precision(N) C The Right Hand Side Vector. C C Job Integer C = 0 To Solve A*X = B , C = Nonzero To Solve Trans(A)*X = B Where C Trans(A) Is The Transpose. C C On Return C C B The Solution Vector X . C C Error Condition C 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 Dgeco Has Set Rcond .gt. 0.0 C Or Dgefa Has Set Info .eq. 0 . C C To Compute Inverse(A) * C Where C Is A Matrix C With P Columns C Call Dgeco(A,Lda,N,Ipvt,Rcond,Z) C If (Rcond Is Too Small) Go To ... C Do 10 J = 1, P C Call Dgesl(A,Lda,N,Ipvt,C(1,J),0) C 10 Continue C C Linpack. This Version Dated 08/14/78 . C Cleve Moler, University Of New Mexico, Argonne National Lab. C C Subroutines And Functions C C Blas Daxpy,Ddot C C Internal Variables C Double Precision Ddot,T Integer K,Kb,L,Nm1 C Nm1 = N - 1 If (Job .ne. 0) Go To 50 C C Job = 0 , Solve A * X = B C First Solve L*Y = B C If (Nm1 .lt. 1) Go To 30 Do 20 K = 1, Nm1 L = Ipvt(K) T = B(L) If (L .eq. K) Go To 10 B(L) = B(K) B(K) = T 10 Continue Call Daxpy(N-K,T,A(K+1,K),1,B(K+1),1) 20 Continue 30 Continue C C Now Solve U*X = Y C Do 40 Kb = 1, N K = N + 1 - Kb B(K) = B(K)/A(K,K) T = -B(K) Call Daxpy(K-1,T,A(1,K),1,B(1),1) 40 Continue Go To 100 50 Continue C C Job = Nonzero, Solve Trans(A) * X = B C First Solve Trans(U)*Y = B C Do 60 K = 1, N T = Ddot(K-1,A(1,K),1,B(1),1) B(K) = (B(K) - T)/A(K,K) 60 Continue C C Now Solve Trans(L)*X = Y C If (Nm1 .lt. 1) Go To 90 Do 80 Kb = 1, Nm1 K = N - Kb B(K) = B(K) + Ddot(N-K,A(K+1,K),1,B(K+1),1) L = Ipvt(K) If (L .eq. K) Go To 70 T = B(L) B(L) = B(K) B(K) = T 70 Continue 80 Continue 90 Continue 100 Continue Return End Subroutine Daxpy ( N, Alpha, X, Incx, Y, Incy ) Implicit Double Precision (A-H,O-Z) Integer N, Incx, Incy Dimension X( * ), Y( * ) C Daxpy Performs The Operation C C Y : = Alpha*X + Y C C C Modified Nag Fortran 77 Version Of The Blas Routine Daxpy . C C -- Written On 3-September-1982. C Sven Hammarling, Nag Central Office. Integer I , Ix , Iy Parameter ( Zero = 0.0d+0 ) If( N .lt. 1 )Return If( Alpha .eq. Zero )Return If( ( Incx .eq. Incy ) .and. ( Incx .gt. 0 ) )Then If ( Alpha .eq. 1.0d+0 ) Then Do 5, Ix = 1, 1 + ( N - 1 )*Incx, Incx Y( Ix ) = X( Ix ) + Y( Ix ) 5 Continue Else Do 10, Ix = 1, 1 + ( N - 1 )*Incx, Incx Y( Ix ) = Alpha*X( Ix ) + Y( Ix ) 10 Continue Endif Else If( Incy .ge. 0 )Then Iy = 1 Else Iy = 1 - ( N - 1 )*Incy End If If( Incx .gt. 0 )Then Do 20, Ix = 1, 1 + ( N - 1 )*Incx, Incx Y( Iy ) = Alpha*X( Ix ) + Y( Iy ) Iy = Iy + Incy 20 Continue Else Ix = 1 - ( N - 1 )*Incx Do 30, I = 1, N Y( Iy ) = Alpha*X( Ix ) + Y( Iy ) Ix = Ix + Incx Iy = Iy + Incy 30 Continue End If End If Return * End Of Daxpy . End Double Precision Function Ddot ( N, X, Incx, Y, Incy ) Implicit Double Precision (A-H,O-Z) Integer N, Incx, Incy Dimension X( * ), Y( * ) C Ddot Returns The Value C C Ddot = X'Y C C C Modified Nag Fortran 77 Version Of The Blas Routine Ddot . C C -- Written On 21-September-1982. C Sven Hammarling, Nag Central Office. Integer I , Ix , Iy Parameter ( Zero = 0.0d+0 ) Sum = Zero If( N .ge. 1 )Then If( ( Incx .eq. Incy ) .and. ( Incx .gt. 0 ) )Then Do 10, Ix = 1, 1 + ( N - 1 )*Incx, Incx Sum = Sum + X( Ix )*Y( Ix ) 10 Continue Else If( Incy .ge. 0 )Then Iy = 1 Else Iy = 1 - ( N - 1 )*Incy End If If( Incx .gt. 0 )Then Do 20, Ix = 1, 1 + ( N - 1 )*Incx, Incx Sum = Sum + X( Ix )*Y( Iy ) Iy = Iy + Incy 20 Continue Else Ix = 1 - ( N - 1 )*Incx Do 30, I = 1, N Sum = Sum + X( Ix )*Y( Iy ) Ix = Ix + Incx Iy = Iy + Incy 30 Continue End If End If End If Ddot = Sum Return * End Of Ddot . End Subroutine Dscal ( N, Alpha, X, Incx ) Implicit Double Precision (A-H,O-Z) Integer N, Incx Dimension X( * ) C Dscal Performs The Operation C C X : = Alpha*X C C C Modified Nag Fortran 77 Version Of The Blas Routine Dscal . C C -- Written On 26-November-1982. C Sven Hammarling, Nag Central Office. Integer Ix Parameter ( One = 1.0d+0, Zero = 0.0d+0 ) If( N .ge. 1 )Then If( Alpha .eq. Zero )Then Do 10, Ix = 1, 1 + ( N - 1 )*Incx, Incx X( Ix ) = Zero 10 Continue Else If( Alpha .eq. ( -One ) )Then Do 20, Ix = 1, 1 + ( N - 1 )*Incx, Incx X( Ix ) = -X( Ix ) 20 Continue Else If( Alpha .ne. One )Then Do 30, Ix = 1, 1 + ( N - 1 )*Incx, Incx X( Ix ) = Alpha*X( Ix ) 30 Continue End If End If Return * End Of Dscal . End Integer Function Idamax( N, X, Incx ) Implicit Double Precision (A-H,O-Z) Integer N, Incx Dimension X( * ) C Idamax Returns The Smallest Value Of I Such That C C Abs( X( I ) ) = Max( Abs( X( J ) ) ) C J C C Nag Fortran 77 Version Of The Blas Routine Idamax. C Nag Fortran 77 O( N ) Basic Linear Algebra Routine. C C -- Written On 31-May-1983. C Sven Hammarling, Nag Central Office. Intrinsic Abs Integer I , Imax , Ix If( N .lt. 1 )Then Idamax = 0 Return End If Imax = 1 If( N .gt. 1 )Then Xmax = Abs( X( 1 ) ) Ix = 1 Do 10, I = 2, N Ix = Ix + Incx If( Xmax .lt. Abs( X( Ix ) ) )Then Xmax = Abs( X( Ix ) ) Imax = I End If 10 Continue End If Idamax = Imax Return * End Of Idamax. End c =================================================================================== subroutine dcopy ( n, x, incx, y, incy ) implicit double precision (a-h,o-z) integer n, incx, incy dimension x( * ), y( * ) c dcopy performs the operation c c y := x c c nag fortran 77 version of the blas routine dcopy . c nag fortran 77 o( n ) basic linear algebra routine. c c -- written on 26-november-1982. c sven hammarling, nag central office. integer i , ix , iy if( n.lt.1 )return if( ( incx.eq.incy ).and.( incy.gt.0 ) )then do 10, iy = 1, 1 + ( n - 1 )*incy, incy y( iy ) = x( iy ) 10 continue else if( incx.ge.0 )then ix = 1 else ix = 1 - ( n - 1 )*incx end if if( incy.gt.0 )then do 20, iy = 1, 1 + ( n - 1 )*incy, incy y( iy ) = x( ix ) ix = ix + incx 20 continue else iy = 1 - ( n - 1 )*incy do 30, i = 1, n y( iy ) = x( ix ) iy = iy + incy ix = ix + incx 30 continue end if end if return * end of dcopy . end