![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SSYTRI 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SSYTRI + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssytri.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssytri.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssytri.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SSYTRI( UPLO, N, A, LDA, IPIV, WORK, INFO ) 00022 * 00023 * .. Scalar Arguments .. 00024 * CHARACTER UPLO 00025 * INTEGER INFO, LDA, N 00026 * .. 00027 * .. Array Arguments .. 00028 * INTEGER IPIV( * ) 00029 * REAL A( LDA, * ), WORK( * ) 00030 * .. 00031 * 00032 * 00033 *> \par Purpose: 00034 * ============= 00035 *> 00036 *> \verbatim 00037 *> 00038 *> SSYTRI computes the inverse of a real symmetric indefinite matrix 00039 *> A using the factorization A = U*D*U**T or A = L*D*L**T computed by 00040 *> SSYTRF. 00041 *> \endverbatim 00042 * 00043 * Arguments: 00044 * ========== 00045 * 00046 *> \param[in] UPLO 00047 *> \verbatim 00048 *> UPLO is CHARACTER*1 00049 *> Specifies whether the details of the factorization are stored 00050 *> as an upper or lower triangular matrix. 00051 *> = 'U': Upper triangular, form is A = U*D*U**T; 00052 *> = 'L': Lower triangular, form is A = L*D*L**T. 00053 *> \endverbatim 00054 *> 00055 *> \param[in] N 00056 *> \verbatim 00057 *> N is INTEGER 00058 *> The order of the matrix A. N >= 0. 00059 *> \endverbatim 00060 *> 00061 *> \param[in,out] A 00062 *> \verbatim 00063 *> A is REAL array, dimension (LDA,N) 00064 *> On entry, the block diagonal matrix D and the multipliers 00065 *> used to obtain the factor U or L as computed by SSYTRF. 00066 *> 00067 *> On exit, if INFO = 0, the (symmetric) inverse of the original 00068 *> matrix. If UPLO = 'U', the upper triangular part of the 00069 *> inverse is formed and the part of A below the diagonal is not 00070 *> referenced; if UPLO = 'L' the lower triangular part of the 00071 *> inverse is formed and the part of A above the diagonal is 00072 *> not referenced. 00073 *> \endverbatim 00074 *> 00075 *> \param[in] LDA 00076 *> \verbatim 00077 *> LDA is INTEGER 00078 *> The leading dimension of the array A. LDA >= max(1,N). 00079 *> \endverbatim 00080 *> 00081 *> \param[in] IPIV 00082 *> \verbatim 00083 *> IPIV is INTEGER array, dimension (N) 00084 *> Details of the interchanges and the block structure of D 00085 *> as determined by SSYTRF. 00086 *> \endverbatim 00087 *> 00088 *> \param[out] WORK 00089 *> \verbatim 00090 *> WORK is REAL array, dimension (N) 00091 *> \endverbatim 00092 *> 00093 *> \param[out] INFO 00094 *> \verbatim 00095 *> INFO is INTEGER 00096 *> = 0: successful exit 00097 *> < 0: if INFO = -i, the i-th argument had an illegal value 00098 *> > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its 00099 *> inverse could not be computed. 00100 *> \endverbatim 00101 * 00102 * Authors: 00103 * ======== 00104 * 00105 *> \author Univ. of Tennessee 00106 *> \author Univ. of California Berkeley 00107 *> \author Univ. of Colorado Denver 00108 *> \author NAG Ltd. 00109 * 00110 *> \date November 2011 00111 * 00112 *> \ingroup realSYcomputational 00113 * 00114 * ===================================================================== 00115 SUBROUTINE SSYTRI( UPLO, N, A, LDA, IPIV, WORK, INFO ) 00116 * 00117 * -- LAPACK computational routine (version 3.4.0) -- 00118 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00119 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00120 * November 2011 00121 * 00122 * .. Scalar Arguments .. 00123 CHARACTER UPLO 00124 INTEGER INFO, LDA, N 00125 * .. 00126 * .. Array Arguments .. 00127 INTEGER IPIV( * ) 00128 REAL A( LDA, * ), WORK( * ) 00129 * .. 00130 * 00131 * ===================================================================== 00132 * 00133 * .. Parameters .. 00134 REAL ONE, ZERO 00135 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 00136 * .. 00137 * .. Local Scalars .. 00138 LOGICAL UPPER 00139 INTEGER K, KP, KSTEP 00140 REAL AK, AKKP1, AKP1, D, T, TEMP 00141 * .. 00142 * .. External Functions .. 00143 LOGICAL LSAME 00144 REAL SDOT 00145 EXTERNAL LSAME, SDOT 00146 * .. 00147 * .. External Subroutines .. 00148 EXTERNAL SCOPY, SSWAP, SSYMV, XERBLA 00149 * .. 00150 * .. Intrinsic Functions .. 00151 INTRINSIC ABS, MAX 00152 * .. 00153 * .. Executable Statements .. 00154 * 00155 * Test the input parameters. 00156 * 00157 INFO = 0 00158 UPPER = LSAME( UPLO, 'U' ) 00159 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00160 INFO = -1 00161 ELSE IF( N.LT.0 ) THEN 00162 INFO = -2 00163 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00164 INFO = -4 00165 END IF 00166 IF( INFO.NE.0 ) THEN 00167 CALL XERBLA( 'SSYTRI', -INFO ) 00168 RETURN 00169 END IF 00170 * 00171 * Quick return if possible 00172 * 00173 IF( N.EQ.0 ) 00174 $ RETURN 00175 * 00176 * Check that the diagonal matrix D is nonsingular. 00177 * 00178 IF( UPPER ) THEN 00179 * 00180 * Upper triangular storage: examine D from bottom to top 00181 * 00182 DO 10 INFO = N, 1, -1 00183 IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO ) 00184 $ RETURN 00185 10 CONTINUE 00186 ELSE 00187 * 00188 * Lower triangular storage: examine D from top to bottom. 00189 * 00190 DO 20 INFO = 1, N 00191 IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO ) 00192 $ RETURN 00193 20 CONTINUE 00194 END IF 00195 INFO = 0 00196 * 00197 IF( UPPER ) THEN 00198 * 00199 * Compute inv(A) from the factorization A = U*D*U**T. 00200 * 00201 * K is the main loop index, increasing from 1 to N in steps of 00202 * 1 or 2, depending on the size of the diagonal blocks. 00203 * 00204 K = 1 00205 30 CONTINUE 00206 * 00207 * If K > N, exit from loop. 00208 * 00209 IF( K.GT.N ) 00210 $ GO TO 40 00211 * 00212 IF( IPIV( K ).GT.0 ) THEN 00213 * 00214 * 1 x 1 diagonal block 00215 * 00216 * Invert the diagonal block. 00217 * 00218 A( K, K ) = ONE / A( K, K ) 00219 * 00220 * Compute column K of the inverse. 00221 * 00222 IF( K.GT.1 ) THEN 00223 CALL SCOPY( K-1, A( 1, K ), 1, WORK, 1 ) 00224 CALL SSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO, 00225 $ A( 1, K ), 1 ) 00226 A( K, K ) = A( K, K ) - SDOT( K-1, WORK, 1, A( 1, K ), 00227 $ 1 ) 00228 END IF 00229 KSTEP = 1 00230 ELSE 00231 * 00232 * 2 x 2 diagonal block 00233 * 00234 * Invert the diagonal block. 00235 * 00236 T = ABS( A( K, K+1 ) ) 00237 AK = A( K, K ) / T 00238 AKP1 = A( K+1, K+1 ) / T 00239 AKKP1 = A( K, K+1 ) / T 00240 D = T*( AK*AKP1-ONE ) 00241 A( K, K ) = AKP1 / D 00242 A( K+1, K+1 ) = AK / D 00243 A( K, K+1 ) = -AKKP1 / D 00244 * 00245 * Compute columns K and K+1 of the inverse. 00246 * 00247 IF( K.GT.1 ) THEN 00248 CALL SCOPY( K-1, A( 1, K ), 1, WORK, 1 ) 00249 CALL SSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO, 00250 $ A( 1, K ), 1 ) 00251 A( K, K ) = A( K, K ) - SDOT( K-1, WORK, 1, A( 1, K ), 00252 $ 1 ) 00253 A( K, K+1 ) = A( K, K+1 ) - 00254 $ SDOT( K-1, A( 1, K ), 1, A( 1, K+1 ), 1 ) 00255 CALL SCOPY( K-1, A( 1, K+1 ), 1, WORK, 1 ) 00256 CALL SSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO, 00257 $ A( 1, K+1 ), 1 ) 00258 A( K+1, K+1 ) = A( K+1, K+1 ) - 00259 $ SDOT( K-1, WORK, 1, A( 1, K+1 ), 1 ) 00260 END IF 00261 KSTEP = 2 00262 END IF 00263 * 00264 KP = ABS( IPIV( K ) ) 00265 IF( KP.NE.K ) THEN 00266 * 00267 * Interchange rows and columns K and KP in the leading 00268 * submatrix A(1:k+1,1:k+1) 00269 * 00270 CALL SSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 ) 00271 CALL SSWAP( K-KP-1, A( KP+1, K ), 1, A( KP, KP+1 ), LDA ) 00272 TEMP = A( K, K ) 00273 A( K, K ) = A( KP, KP ) 00274 A( KP, KP ) = TEMP 00275 IF( KSTEP.EQ.2 ) THEN 00276 TEMP = A( K, K+1 ) 00277 A( K, K+1 ) = A( KP, K+1 ) 00278 A( KP, K+1 ) = TEMP 00279 END IF 00280 END IF 00281 * 00282 K = K + KSTEP 00283 GO TO 30 00284 40 CONTINUE 00285 * 00286 ELSE 00287 * 00288 * Compute inv(A) from the factorization A = L*D*L**T. 00289 * 00290 * K is the main loop index, increasing from 1 to N in steps of 00291 * 1 or 2, depending on the size of the diagonal blocks. 00292 * 00293 K = N 00294 50 CONTINUE 00295 * 00296 * If K < 1, exit from loop. 00297 * 00298 IF( K.LT.1 ) 00299 $ GO TO 60 00300 * 00301 IF( IPIV( K ).GT.0 ) THEN 00302 * 00303 * 1 x 1 diagonal block 00304 * 00305 * Invert the diagonal block. 00306 * 00307 A( K, K ) = ONE / A( K, K ) 00308 * 00309 * Compute column K of the inverse. 00310 * 00311 IF( K.LT.N ) THEN 00312 CALL SCOPY( N-K, A( K+1, K ), 1, WORK, 1 ) 00313 CALL SSYMV( UPLO, N-K, -ONE, A( K+1, K+1 ), LDA, WORK, 1, 00314 $ ZERO, A( K+1, K ), 1 ) 00315 A( K, K ) = A( K, K ) - SDOT( N-K, WORK, 1, A( K+1, K ), 00316 $ 1 ) 00317 END IF 00318 KSTEP = 1 00319 ELSE 00320 * 00321 * 2 x 2 diagonal block 00322 * 00323 * Invert the diagonal block. 00324 * 00325 T = ABS( A( K, K-1 ) ) 00326 AK = A( K-1, K-1 ) / T 00327 AKP1 = A( K, K ) / T 00328 AKKP1 = A( K, K-1 ) / T 00329 D = T*( AK*AKP1-ONE ) 00330 A( K-1, K-1 ) = AKP1 / D 00331 A( K, K ) = AK / D 00332 A( K, K-1 ) = -AKKP1 / D 00333 * 00334 * Compute columns K-1 and K of the inverse. 00335 * 00336 IF( K.LT.N ) THEN 00337 CALL SCOPY( N-K, A( K+1, K ), 1, WORK, 1 ) 00338 CALL SSYMV( UPLO, N-K, -ONE, A( K+1, K+1 ), LDA, WORK, 1, 00339 $ ZERO, A( K+1, K ), 1 ) 00340 A( K, K ) = A( K, K ) - SDOT( N-K, WORK, 1, A( K+1, K ), 00341 $ 1 ) 00342 A( K, K-1 ) = A( K, K-1 ) - 00343 $ SDOT( N-K, A( K+1, K ), 1, A( K+1, K-1 ), 00344 $ 1 ) 00345 CALL SCOPY( N-K, A( K+1, K-1 ), 1, WORK, 1 ) 00346 CALL SSYMV( UPLO, N-K, -ONE, A( K+1, K+1 ), LDA, WORK, 1, 00347 $ ZERO, A( K+1, K-1 ), 1 ) 00348 A( K-1, K-1 ) = A( K-1, K-1 ) - 00349 $ SDOT( N-K, WORK, 1, A( K+1, K-1 ), 1 ) 00350 END IF 00351 KSTEP = 2 00352 END IF 00353 * 00354 KP = ABS( IPIV( K ) ) 00355 IF( KP.NE.K ) THEN 00356 * 00357 * Interchange rows and columns K and KP in the trailing 00358 * submatrix A(k-1:n,k-1:n) 00359 * 00360 IF( KP.LT.N ) 00361 $ CALL SSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 ) 00362 CALL SSWAP( KP-K-1, A( K+1, K ), 1, A( KP, K+1 ), LDA ) 00363 TEMP = A( K, K ) 00364 A( K, K ) = A( KP, KP ) 00365 A( KP, KP ) = TEMP 00366 IF( KSTEP.EQ.2 ) THEN 00367 TEMP = A( K, K-1 ) 00368 A( K, K-1 ) = A( KP, K-1 ) 00369 A( KP, K-1 ) = TEMP 00370 END IF 00371 END IF 00372 * 00373 K = K - KSTEP 00374 GO TO 50 00375 60 CONTINUE 00376 END IF 00377 * 00378 RETURN 00379 * 00380 * End of SSYTRI 00381 * 00382 END