00001 SUBROUTINE ZGETRF( 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 COMPLEX*16 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 COMPLEX*16 ONE, NEGONE
00097 DOUBLE PRECISION ZERO
00098 PARAMETER ( ONE = (1.0D+0, 0.0D+0) )
00099 PARAMETER ( NEGONE = (-1.0D+0, 0.0D+0) )
00100 PARAMETER ( ZERO = 0.0D+0 )
00101
00102
00103 DOUBLE PRECISION SFMIN, PIVMAG
00104 COMPLEX*16 TMP
00105 INTEGER I, J, JP, NSTEP, NTOPIV, NPIVED, KAHEAD
00106 INTEGER KSTART, IPIVSTART, JPIVSTART, KCOLS
00107
00108
00109 DOUBLE PRECISION DLAMCH
00110 INTEGER IZAMAX
00111 LOGICAL DISNAN
00112 EXTERNAL DLAMCH, IZAMAX, DISNAN
00113
00114
00115 EXTERNAL ZTRSM, ZSCAL, XERBLA, ZLASWP
00116
00117
00118 INTRINSIC MAX, MIN, IAND, ABS
00119
00120
00121
00122
00123
00124 INFO = 0
00125 IF( M.LT.0 ) THEN
00126 INFO = -1
00127 ELSE IF( N.LT.0 ) THEN
00128 INFO = -2
00129 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00130 INFO = -4
00131 END IF
00132 IF( INFO.NE.0 ) THEN
00133 CALL XERBLA( 'ZGETRF', -INFO )
00134 RETURN
00135 END IF
00136
00137
00138
00139 IF( M.EQ.0 .OR. N.EQ.0 )
00140 $ RETURN
00141
00142
00143
00144 SFMIN = DLAMCH( 'S' )
00145
00146 NSTEP = MIN( M, N )
00147 DO J = 1, NSTEP
00148 KAHEAD = IAND( J, -J )
00149 KSTART = J + 1 - KAHEAD
00150 KCOLS = MIN( KAHEAD, M-J )
00151
00152
00153
00154 JP = J - 1 + IZAMAX( M-J+1, A( J, J ), 1 )
00155 IPIV( J ) = JP
00156
00157
00158 IF (JP .NE. J) THEN
00159 TMP = A( J, J )
00160 A( J, J ) = A( JP, J )
00161 A( JP, J ) = TMP
00162 END IF
00163
00164
00165 NTOPIV = 1
00166 IPIVSTART = J
00167 JPIVSTART = J - NTOPIV
00168 DO WHILE ( NTOPIV .LT. KAHEAD )
00169 CALL ZLASWP( NTOPIV, A( 1, JPIVSTART ), LDA, IPIVSTART, J,
00170 $ IPIV, 1 )
00171 IPIVSTART = IPIVSTART - NTOPIV;
00172 NTOPIV = NTOPIV * 2;
00173 JPIVSTART = JPIVSTART - NTOPIV;
00174 END DO
00175
00176
00177 CALL ZLASWP( KCOLS, A( 1,J+1 ), LDA, KSTART, J, IPIV, 1 )
00178
00179
00180 PIVMAG = ABS( A( J, J ) )
00181 IF( PIVMAG.NE.ZERO .AND. .NOT.DISNAN( PIVMAG ) ) THEN
00182 IF( PIVMAG .GE. SFMIN ) THEN
00183 CALL ZSCAL( M-J, ONE / A( J, J ), A( J+1, J ), 1 )
00184 ELSE
00185 DO I = 1, M-J
00186 A( J+I, J ) = A( J+I, J ) / A( J, J )
00187 END DO
00188 END IF
00189 ELSE IF( PIVMAG .EQ. ZERO .AND. INFO .EQ. 0 ) THEN
00190 INFO = J
00191 END IF
00192
00193
00194 CALL ZTRSM( 'Left', 'Lower', 'No transpose', 'Unit', KAHEAD,
00195 $ KCOLS, ONE, A( KSTART, KSTART ), LDA,
00196 $ A( KSTART, J+1 ), LDA )
00197
00198 CALL ZGEMM( 'No transpose', 'No transpose', M-J,
00199 $ KCOLS, KAHEAD, NEGONE, A( J+1, KSTART ), LDA,
00200 $ A( KSTART, J+1 ), LDA, ONE, A( J+1, J+1 ), LDA )
00201 END DO
00202
00203
00204 NPIVED = IAND( NSTEP, -NSTEP )
00205 J = NSTEP - NPIVED
00206 DO WHILE ( J .GT. 0 )
00207 NTOPIV = IAND( J, -J )
00208 CALL ZLASWP( NTOPIV, A( 1, J-NTOPIV+1 ), LDA, J+1, NSTEP,
00209 $ IPIV, 1 )
00210 J = J - NTOPIV
00211 END DO
00212
00213
00214 IF ( M .LT. N ) THEN
00215 CALL ZLASWP( N-M, A( 1, M+KCOLS+1 ), LDA, 1, M, IPIV, 1 )
00216 CALL ZTRSM( 'Left', 'Lower', 'No transpose', 'Unit', M,
00217 $ N-M, ONE, A, LDA, A( 1,M+KCOLS+1 ), LDA )
00218 END IF
00219
00220 RETURN
00221
00222
00223
00224 END