00001 SUBROUTINE DGETRF( M, N, A, LDA, IPIV, INFO )
00002 IMPLICIT NONE
00003
00004
00005
00006
00007
00008
00009 INTEGER INFO, LDA, M, N
00010
00011
00012 INTEGER IPIV( * )
00013 DOUBLE PRECISION A( LDA, * )
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096 DOUBLE PRECISION ONE, ZERO, NEGONE
00097 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
00098 PARAMETER ( NEGONE = -1.0D+0 )
00099
00100
00101 DOUBLE PRECISION SFMIN, TMP
00102 INTEGER I, J, JP, NSTEP, NTOPIV, NPIVED, KAHEAD
00103 INTEGER KSTART, IPIVSTART, JPIVSTART, KCOLS
00104
00105
00106 DOUBLE PRECISION DLAMCH
00107 INTEGER IDAMAX
00108 LOGICAL DISNAN
00109 EXTERNAL DLAMCH, IDAMAX, DISNAN
00110
00111
00112 EXTERNAL DTRSM, DSCAL, XERBLA, DLASWP
00113
00114
00115 INTRINSIC MAX, MIN, IAND
00116
00117
00118
00119
00120
00121 INFO = 0
00122 IF( M.LT.0 ) THEN
00123 INFO = -1
00124 ELSE IF( N.LT.0 ) THEN
00125 INFO = -2
00126 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00127 INFO = -4
00128 END IF
00129 IF( INFO.NE.0 ) THEN
00130 CALL XERBLA( 'DGETRF', -INFO )
00131 RETURN
00132 END IF
00133
00134
00135
00136 IF( M.EQ.0 .OR. N.EQ.0 )
00137 $ RETURN
00138
00139
00140
00141 SFMIN = DLAMCH( 'S' )
00142
00143 NSTEP = MIN( M, N )
00144 DO J = 1, NSTEP
00145 KAHEAD = IAND( J, -J )
00146 KSTART = J + 1 - KAHEAD
00147 KCOLS = MIN( KAHEAD, M-J )
00148
00149
00150
00151 JP = J - 1 + IDAMAX( M-J+1, A( J, J ), 1 )
00152 IPIV( J ) = JP
00153
00154
00155 IF (JP .NE. J) THEN
00156 TMP = A( J, J )
00157 A( J, J ) = A( JP, J )
00158 A( JP, J ) = TMP
00159 END IF
00160
00161
00162 NTOPIV = 1
00163 IPIVSTART = J
00164 JPIVSTART = J - NTOPIV
00165 DO WHILE ( NTOPIV .LT. KAHEAD )
00166 CALL DLASWP( NTOPIV, A( 1, JPIVSTART ), LDA, IPIVSTART, J,
00167 $ IPIV, 1 )
00168 IPIVSTART = IPIVSTART - NTOPIV;
00169 NTOPIV = NTOPIV * 2;
00170 JPIVSTART = JPIVSTART - NTOPIV;
00171 END DO
00172
00173
00174 CALL DLASWP( KCOLS, A( 1,J+1 ), LDA, KSTART, J, IPIV, 1 )
00175
00176
00177 IF( A( J, J ).NE.ZERO .AND. .NOT.DISNAN( A( J, J ) ) ) THEN
00178 IF( ABS(A( J, J )) .GE. SFMIN ) THEN
00179 CALL DSCAL( M-J, ONE / A( J, J ), A( J+1, J ), 1 )
00180 ELSE
00181 DO I = 1, M-J
00182 A( J+I, J ) = A( J+I, J ) / A( J, J )
00183 END DO
00184 END IF
00185 ELSE IF( A( J,J ) .EQ. ZERO .AND. INFO .EQ. 0 ) THEN
00186 INFO = J
00187 END IF
00188
00189
00190 CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Unit', KAHEAD,
00191 $ KCOLS, ONE, A( KSTART, KSTART ), LDA,
00192 $ A( KSTART, J+1 ), LDA )
00193
00194 CALL DGEMM( 'No transpose', 'No transpose', M-J,
00195 $ KCOLS, KAHEAD, NEGONE, A( J+1, KSTART ), LDA,
00196 $ A( KSTART, J+1 ), LDA, ONE, A( J+1, J+1 ), LDA )
00197 END DO
00198
00199
00200 NPIVED = IAND( NSTEP, -NSTEP )
00201 J = NSTEP - NPIVED
00202 DO WHILE ( J .GT. 0 )
00203 NTOPIV = IAND( J, -J )
00204 CALL DLASWP( NTOPIV, A( 1, J-NTOPIV+1 ), LDA, J+1, NSTEP,
00205 $ IPIV, 1 )
00206 J = J - NTOPIV
00207 END DO
00208
00209
00210 IF ( M .LT. N ) THEN
00211 CALL DLASWP( N-M, A( 1, M+KCOLS+1 ), LDA, 1, M, IPIV, 1 )
00212 CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Unit', M,
00213 $ N-M, ONE, A, LDA, A( 1,M+KCOLS+1 ), LDA )
00214 END IF
00215
00216 RETURN
00217
00218
00219
00220 END