00001 SUBROUTINE SSYTRI( UPLO, N, A, LDA, IPIV, WORK, INFO )
00002
00003
00004
00005
00006
00007
00008
00009 CHARACTER UPLO
00010 INTEGER INFO, LDA, N
00011
00012
00013 INTEGER IPIV( * )
00014 REAL A( LDA, * ), WORK( * )
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 REAL ONE, ZERO
00066 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 )
00067
00068
00069 LOGICAL UPPER
00070 INTEGER K, KP, KSTEP
00071 REAL AK, AKKP1, AKP1, D, T, TEMP
00072
00073
00074 LOGICAL LSAME
00075 REAL SDOT
00076 EXTERNAL LSAME, SDOT
00077
00078
00079 EXTERNAL SCOPY, SSWAP, SSYMV, XERBLA
00080
00081
00082 INTRINSIC ABS, MAX
00083
00084
00085
00086
00087
00088 INFO = 0
00089 UPPER = LSAME( UPLO, 'U' )
00090 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00091 INFO = -1
00092 ELSE IF( N.LT.0 ) THEN
00093 INFO = -2
00094 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00095 INFO = -4
00096 END IF
00097 IF( INFO.NE.0 ) THEN
00098 CALL XERBLA( 'SSYTRI', -INFO )
00099 RETURN
00100 END IF
00101
00102
00103
00104 IF( N.EQ.0 )
00105 $ RETURN
00106
00107
00108
00109 IF( UPPER ) THEN
00110
00111
00112
00113 DO 10 INFO = N, 1, -1
00114 IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
00115 $ RETURN
00116 10 CONTINUE
00117 ELSE
00118
00119
00120
00121 DO 20 INFO = 1, N
00122 IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
00123 $ RETURN
00124 20 CONTINUE
00125 END IF
00126 INFO = 0
00127
00128 IF( UPPER ) THEN
00129
00130
00131
00132
00133
00134
00135 K = 1
00136 30 CONTINUE
00137
00138
00139
00140 IF( K.GT.N )
00141 $ GO TO 40
00142
00143 IF( IPIV( K ).GT.0 ) THEN
00144
00145
00146
00147
00148
00149 A( K, K ) = ONE / A( K, K )
00150
00151
00152
00153 IF( K.GT.1 ) THEN
00154 CALL SCOPY( K-1, A( 1, K ), 1, WORK, 1 )
00155 CALL SSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO,
00156 $ A( 1, K ), 1 )
00157 A( K, K ) = A( K, K ) - SDOT( K-1, WORK, 1, A( 1, K ),
00158 $ 1 )
00159 END IF
00160 KSTEP = 1
00161 ELSE
00162
00163
00164
00165
00166
00167 T = ABS( A( K, K+1 ) )
00168 AK = A( K, K ) / T
00169 AKP1 = A( K+1, K+1 ) / T
00170 AKKP1 = A( K, K+1 ) / T
00171 D = T*( AK*AKP1-ONE )
00172 A( K, K ) = AKP1 / D
00173 A( K+1, K+1 ) = AK / D
00174 A( K, K+1 ) = -AKKP1 / D
00175
00176
00177
00178 IF( K.GT.1 ) THEN
00179 CALL SCOPY( K-1, A( 1, K ), 1, WORK, 1 )
00180 CALL SSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO,
00181 $ A( 1, K ), 1 )
00182 A( K, K ) = A( K, K ) - SDOT( K-1, WORK, 1, A( 1, K ),
00183 $ 1 )
00184 A( K, K+1 ) = A( K, K+1 ) -
00185 $ SDOT( K-1, A( 1, K ), 1, A( 1, K+1 ), 1 )
00186 CALL SCOPY( K-1, A( 1, K+1 ), 1, WORK, 1 )
00187 CALL SSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO,
00188 $ A( 1, K+1 ), 1 )
00189 A( K+1, K+1 ) = A( K+1, K+1 ) -
00190 $ SDOT( K-1, WORK, 1, A( 1, K+1 ), 1 )
00191 END IF
00192 KSTEP = 2
00193 END IF
00194
00195 KP = ABS( IPIV( K ) )
00196 IF( KP.NE.K ) THEN
00197
00198
00199
00200
00201 CALL SSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
00202 CALL SSWAP( K-KP-1, A( KP+1, K ), 1, A( KP, KP+1 ), LDA )
00203 TEMP = A( K, K )
00204 A( K, K ) = A( KP, KP )
00205 A( KP, KP ) = TEMP
00206 IF( KSTEP.EQ.2 ) THEN
00207 TEMP = A( K, K+1 )
00208 A( K, K+1 ) = A( KP, K+1 )
00209 A( KP, K+1 ) = TEMP
00210 END IF
00211 END IF
00212
00213 K = K + KSTEP
00214 GO TO 30
00215 40 CONTINUE
00216
00217 ELSE
00218
00219
00220
00221
00222
00223
00224 K = N
00225 50 CONTINUE
00226
00227
00228
00229 IF( K.LT.1 )
00230 $ GO TO 60
00231
00232 IF( IPIV( K ).GT.0 ) THEN
00233
00234
00235
00236
00237
00238 A( K, K ) = ONE / A( K, K )
00239
00240
00241
00242 IF( K.LT.N ) THEN
00243 CALL SCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
00244 CALL SSYMV( UPLO, N-K, -ONE, A( K+1, K+1 ), LDA, WORK, 1,
00245 $ ZERO, A( K+1, K ), 1 )
00246 A( K, K ) = A( K, K ) - SDOT( N-K, WORK, 1, A( K+1, K ),
00247 $ 1 )
00248 END IF
00249 KSTEP = 1
00250 ELSE
00251
00252
00253
00254
00255
00256 T = ABS( A( K, K-1 ) )
00257 AK = A( K-1, K-1 ) / T
00258 AKP1 = A( K, K ) / T
00259 AKKP1 = A( K, K-1 ) / T
00260 D = T*( AK*AKP1-ONE )
00261 A( K-1, K-1 ) = AKP1 / D
00262 A( K, K ) = AK / D
00263 A( K, K-1 ) = -AKKP1 / D
00264
00265
00266
00267 IF( K.LT.N ) THEN
00268 CALL SCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
00269 CALL SSYMV( UPLO, N-K, -ONE, A( K+1, K+1 ), LDA, WORK, 1,
00270 $ ZERO, A( K+1, K ), 1 )
00271 A( K, K ) = A( K, K ) - SDOT( N-K, WORK, 1, A( K+1, K ),
00272 $ 1 )
00273 A( K, K-1 ) = A( K, K-1 ) -
00274 $ SDOT( N-K, A( K+1, K ), 1, A( K+1, K-1 ),
00275 $ 1 )
00276 CALL SCOPY( N-K, A( K+1, K-1 ), 1, WORK, 1 )
00277 CALL SSYMV( UPLO, N-K, -ONE, A( K+1, K+1 ), LDA, WORK, 1,
00278 $ ZERO, A( K+1, K-1 ), 1 )
00279 A( K-1, K-1 ) = A( K-1, K-1 ) -
00280 $ SDOT( N-K, WORK, 1, A( K+1, K-1 ), 1 )
00281 END IF
00282 KSTEP = 2
00283 END IF
00284
00285 KP = ABS( IPIV( K ) )
00286 IF( KP.NE.K ) THEN
00287
00288
00289
00290
00291 IF( KP.LT.N )
00292 $ CALL SSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
00293 CALL SSWAP( KP-K-1, A( K+1, K ), 1, A( KP, K+1 ), LDA )
00294 TEMP = A( K, K )
00295 A( K, K ) = A( KP, KP )
00296 A( KP, KP ) = TEMP
00297 IF( KSTEP.EQ.2 ) THEN
00298 TEMP = A( K, K-1 )
00299 A( K, K-1 ) = A( KP, K-1 )
00300 A( KP, K-1 ) = TEMP
00301 END IF
00302 END IF
00303
00304 K = K - KSTEP
00305 GO TO 50
00306 60 CONTINUE
00307 END IF
00308
00309 RETURN
00310
00311
00312
00313 END