![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DSYTRF 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DSYTRF + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsytrf.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsytrf.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsytrf.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DSYTRF( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO ) 00022 * 00023 * .. Scalar Arguments .. 00024 * CHARACTER UPLO 00025 * INTEGER INFO, LDA, LWORK, N 00026 * .. 00027 * .. Array Arguments .. 00028 * INTEGER IPIV( * ) 00029 * DOUBLE PRECISION A( LDA, * ), WORK( * ) 00030 * .. 00031 * 00032 * 00033 *> \par Purpose: 00034 * ============= 00035 *> 00036 *> \verbatim 00037 *> 00038 *> DSYTRF computes the factorization of a real symmetric matrix A using 00039 *> the Bunch-Kaufman diagonal pivoting method. The form of the 00040 *> factorization is 00041 *> 00042 *> A = U*D*U**T or A = L*D*L**T 00043 *> 00044 *> where U (or L) is a product of permutation and unit upper (lower) 00045 *> triangular matrices, and D is symmetric and block diagonal with 00046 *> 1-by-1 and 2-by-2 diagonal blocks. 00047 *> 00048 *> This is the blocked version of the algorithm, calling Level 3 BLAS. 00049 *> \endverbatim 00050 * 00051 * Arguments: 00052 * ========== 00053 * 00054 *> \param[in] UPLO 00055 *> \verbatim 00056 *> UPLO is CHARACTER*1 00057 *> = 'U': Upper triangle of A is stored; 00058 *> = 'L': Lower triangle of A is stored. 00059 *> \endverbatim 00060 *> 00061 *> \param[in] N 00062 *> \verbatim 00063 *> N is INTEGER 00064 *> The order of the matrix A. N >= 0. 00065 *> \endverbatim 00066 *> 00067 *> \param[in,out] A 00068 *> \verbatim 00069 *> A is DOUBLE PRECISION array, dimension (LDA,N) 00070 *> On entry, the symmetric matrix A. If UPLO = 'U', the leading 00071 *> N-by-N upper triangular part of A contains the upper 00072 *> triangular part of the matrix A, and the strictly lower 00073 *> triangular part of A is not referenced. If UPLO = 'L', the 00074 *> leading N-by-N lower triangular part of A contains the lower 00075 *> triangular part of the matrix A, and the strictly upper 00076 *> triangular part of A is not referenced. 00077 *> 00078 *> On exit, the block diagonal matrix D and the multipliers used 00079 *> to obtain the factor U or L (see below for further details). 00080 *> \endverbatim 00081 *> 00082 *> \param[in] LDA 00083 *> \verbatim 00084 *> LDA is INTEGER 00085 *> The leading dimension of the array A. LDA >= max(1,N). 00086 *> \endverbatim 00087 *> 00088 *> \param[out] IPIV 00089 *> \verbatim 00090 *> IPIV is INTEGER array, dimension (N) 00091 *> Details of the interchanges and the block structure of D. 00092 *> If IPIV(k) > 0, then rows and columns k and IPIV(k) were 00093 *> interchanged and D(k,k) is a 1-by-1 diagonal block. 00094 *> If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and 00095 *> columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k) 00096 *> is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) = 00097 *> IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were 00098 *> interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block. 00099 *> \endverbatim 00100 *> 00101 *> \param[out] WORK 00102 *> \verbatim 00103 *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) 00104 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00105 *> \endverbatim 00106 *> 00107 *> \param[in] LWORK 00108 *> \verbatim 00109 *> LWORK is INTEGER 00110 *> The length of WORK. LWORK >=1. For best performance 00111 *> LWORK >= N*NB, where NB is the block size returned by ILAENV. 00112 *> 00113 *> If LWORK = -1, then a workspace query is assumed; the routine 00114 *> only calculates the optimal size of the WORK array, returns 00115 *> this value as the first entry of the WORK array, and no error 00116 *> message related to LWORK is issued by XERBLA. 00117 *> \endverbatim 00118 *> 00119 *> \param[out] INFO 00120 *> \verbatim 00121 *> INFO is INTEGER 00122 *> = 0: successful exit 00123 *> < 0: if INFO = -i, the i-th argument had an illegal value 00124 *> > 0: if INFO = i, D(i,i) is exactly zero. The factorization 00125 *> has been completed, but the block diagonal matrix D is 00126 *> exactly singular, and division by zero will occur if it 00127 *> is used to solve a system of equations. 00128 *> \endverbatim 00129 * 00130 * Authors: 00131 * ======== 00132 * 00133 *> \author Univ. of Tennessee 00134 *> \author Univ. of California Berkeley 00135 *> \author Univ. of Colorado Denver 00136 *> \author NAG Ltd. 00137 * 00138 *> \date November 2011 00139 * 00140 *> \ingroup doubleSYcomputational 00141 * 00142 *> \par Further Details: 00143 * ===================== 00144 *> 00145 *> \verbatim 00146 *> 00147 *> If UPLO = 'U', then A = U*D*U**T, where 00148 *> U = P(n)*U(n)* ... *P(k)U(k)* ..., 00149 *> i.e., U is a product of terms P(k)*U(k), where k decreases from n to 00150 *> 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1 00151 *> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as 00152 *> defined by IPIV(k), and U(k) is a unit upper triangular matrix, such 00153 *> that if the diagonal block D(k) is of order s (s = 1 or 2), then 00154 *> 00155 *> ( I v 0 ) k-s 00156 *> U(k) = ( 0 I 0 ) s 00157 *> ( 0 0 I ) n-k 00158 *> k-s s n-k 00159 *> 00160 *> If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k). 00161 *> If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k), 00162 *> and A(k,k), and v overwrites A(1:k-2,k-1:k). 00163 *> 00164 *> If UPLO = 'L', then A = L*D*L**T, where 00165 *> L = P(1)*L(1)* ... *P(k)*L(k)* ..., 00166 *> i.e., L is a product of terms P(k)*L(k), where k increases from 1 to 00167 *> n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1 00168 *> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as 00169 *> defined by IPIV(k), and L(k) is a unit lower triangular matrix, such 00170 *> that if the diagonal block D(k) is of order s (s = 1 or 2), then 00171 *> 00172 *> ( I 0 0 ) k-1 00173 *> L(k) = ( 0 I 0 ) s 00174 *> ( 0 v I ) n-k-s+1 00175 *> k-1 s n-k-s+1 00176 *> 00177 *> If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k). 00178 *> If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k), 00179 *> and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1). 00180 *> \endverbatim 00181 *> 00182 * ===================================================================== 00183 SUBROUTINE DSYTRF( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO ) 00184 * 00185 * -- LAPACK computational routine (version 3.4.0) -- 00186 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00187 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00188 * November 2011 00189 * 00190 * .. Scalar Arguments .. 00191 CHARACTER UPLO 00192 INTEGER INFO, LDA, LWORK, N 00193 * .. 00194 * .. Array Arguments .. 00195 INTEGER IPIV( * ) 00196 DOUBLE PRECISION A( LDA, * ), WORK( * ) 00197 * .. 00198 * 00199 * ===================================================================== 00200 * 00201 * .. Local Scalars .. 00202 LOGICAL LQUERY, UPPER 00203 INTEGER IINFO, IWS, J, K, KB, LDWORK, LWKOPT, NB, NBMIN 00204 * .. 00205 * .. External Functions .. 00206 LOGICAL LSAME 00207 INTEGER ILAENV 00208 EXTERNAL LSAME, ILAENV 00209 * .. 00210 * .. External Subroutines .. 00211 EXTERNAL DLASYF, DSYTF2, XERBLA 00212 * .. 00213 * .. Intrinsic Functions .. 00214 INTRINSIC MAX 00215 * .. 00216 * .. Executable Statements .. 00217 * 00218 * Test the input parameters. 00219 * 00220 INFO = 0 00221 UPPER = LSAME( UPLO, 'U' ) 00222 LQUERY = ( LWORK.EQ.-1 ) 00223 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00224 INFO = -1 00225 ELSE IF( N.LT.0 ) THEN 00226 INFO = -2 00227 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00228 INFO = -4 00229 ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN 00230 INFO = -7 00231 END IF 00232 * 00233 IF( INFO.EQ.0 ) THEN 00234 * 00235 * Determine the block size 00236 * 00237 NB = ILAENV( 1, 'DSYTRF', UPLO, N, -1, -1, -1 ) 00238 LWKOPT = N*NB 00239 WORK( 1 ) = LWKOPT 00240 END IF 00241 * 00242 IF( INFO.NE.0 ) THEN 00243 CALL XERBLA( 'DSYTRF', -INFO ) 00244 RETURN 00245 ELSE IF( LQUERY ) THEN 00246 RETURN 00247 END IF 00248 * 00249 NBMIN = 2 00250 LDWORK = N 00251 IF( NB.GT.1 .AND. NB.LT.N ) THEN 00252 IWS = LDWORK*NB 00253 IF( LWORK.LT.IWS ) THEN 00254 NB = MAX( LWORK / LDWORK, 1 ) 00255 NBMIN = MAX( 2, ILAENV( 2, 'DSYTRF', UPLO, N, -1, -1, -1 ) ) 00256 END IF 00257 ELSE 00258 IWS = 1 00259 END IF 00260 IF( NB.LT.NBMIN ) 00261 $ NB = N 00262 * 00263 IF( UPPER ) THEN 00264 * 00265 * Factorize A as U*D*U**T using the upper triangle of A 00266 * 00267 * K is the main loop index, decreasing from N to 1 in steps of 00268 * KB, where KB is the number of columns factorized by DLASYF; 00269 * KB is either NB or NB-1, or K for the last block 00270 * 00271 K = N 00272 10 CONTINUE 00273 * 00274 * If K < 1, exit from loop 00275 * 00276 IF( K.LT.1 ) 00277 $ GO TO 40 00278 * 00279 IF( K.GT.NB ) THEN 00280 * 00281 * Factorize columns k-kb+1:k of A and use blocked code to 00282 * update columns 1:k-kb 00283 * 00284 CALL DLASYF( UPLO, K, NB, KB, A, LDA, IPIV, WORK, LDWORK, 00285 $ IINFO ) 00286 ELSE 00287 * 00288 * Use unblocked code to factorize columns 1:k of A 00289 * 00290 CALL DSYTF2( UPLO, K, A, LDA, IPIV, IINFO ) 00291 KB = K 00292 END IF 00293 * 00294 * Set INFO on the first occurrence of a zero pivot 00295 * 00296 IF( INFO.EQ.0 .AND. IINFO.GT.0 ) 00297 $ INFO = IINFO 00298 * 00299 * Decrease K and return to the start of the main loop 00300 * 00301 K = K - KB 00302 GO TO 10 00303 * 00304 ELSE 00305 * 00306 * Factorize A as L*D*L**T using the lower triangle of A 00307 * 00308 * K is the main loop index, increasing from 1 to N in steps of 00309 * KB, where KB is the number of columns factorized by DLASYF; 00310 * KB is either NB or NB-1, or N-K+1 for the last block 00311 * 00312 K = 1 00313 20 CONTINUE 00314 * 00315 * If K > N, exit from loop 00316 * 00317 IF( K.GT.N ) 00318 $ GO TO 40 00319 * 00320 IF( K.LE.N-NB ) THEN 00321 * 00322 * Factorize columns k:k+kb-1 of A and use blocked code to 00323 * update columns k+kb:n 00324 * 00325 CALL DLASYF( UPLO, N-K+1, NB, KB, A( K, K ), LDA, IPIV( K ), 00326 $ WORK, LDWORK, IINFO ) 00327 ELSE 00328 * 00329 * Use unblocked code to factorize columns k:n of A 00330 * 00331 CALL DSYTF2( UPLO, N-K+1, A( K, K ), LDA, IPIV( K ), IINFO ) 00332 KB = N - K + 1 00333 END IF 00334 * 00335 * Set INFO on the first occurrence of a zero pivot 00336 * 00337 IF( INFO.EQ.0 .AND. IINFO.GT.0 ) 00338 $ INFO = IINFO + K - 1 00339 * 00340 * Adjust IPIV 00341 * 00342 DO 30 J = K, K + KB - 1 00343 IF( IPIV( J ).GT.0 ) THEN 00344 IPIV( J ) = IPIV( J ) + K - 1 00345 ELSE 00346 IPIV( J ) = IPIV( J ) - K + 1 00347 END IF 00348 30 CONTINUE 00349 * 00350 * Increase K and return to the start of the main loop 00351 * 00352 K = K + KB 00353 GO TO 20 00354 * 00355 END IF 00356 * 00357 40 CONTINUE 00358 WORK( 1 ) = LWKOPT 00359 RETURN 00360 * 00361 * End of DSYTRF 00362 * 00363 END