|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE DPBTRF( UPLO, N, KD, AB, LDAB, INFO ) 00002 * 00003 * -- LAPACK routine (version 3.2) -- 00004 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00005 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00006 * November 2006 00007 * 00008 * .. Scalar Arguments .. 00009 CHARACTER UPLO 00010 INTEGER INFO, KD, LDAB, N 00011 * .. 00012 * .. Array Arguments .. 00013 DOUBLE PRECISION AB( LDAB, * ) 00014 * .. 00015 * 00016 * Purpose 00017 * ======= 00018 * 00019 * DPBTRF computes the Cholesky factorization of a real symmetric 00020 * positive definite band matrix A. 00021 * 00022 * The factorization has the form 00023 * A = U**T * U, if UPLO = 'U', or 00024 * A = L * L**T, if UPLO = 'L', 00025 * where U is an upper triangular matrix and L is lower triangular. 00026 * 00027 * Arguments 00028 * ========= 00029 * 00030 * UPLO (input) CHARACTER*1 00031 * = 'U': Upper triangle of A is stored; 00032 * = 'L': Lower triangle of A is stored. 00033 * 00034 * N (input) INTEGER 00035 * The order of the matrix A. N >= 0. 00036 * 00037 * KD (input) INTEGER 00038 * The number of superdiagonals of the matrix A if UPLO = 'U', 00039 * or the number of subdiagonals if UPLO = 'L'. KD >= 0. 00040 * 00041 * AB (input/output) DOUBLE PRECISION array, dimension (LDAB,N) 00042 * On entry, the upper or lower triangle of the symmetric band 00043 * matrix A, stored in the first KD+1 rows of the array. The 00044 * j-th column of A is stored in the j-th column of the array AB 00045 * as follows: 00046 * if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j; 00047 * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd). 00048 * 00049 * On exit, if INFO = 0, the triangular factor U or L from the 00050 * Cholesky factorization A = U**T*U or A = L*L**T of the band 00051 * matrix A, in the same storage format as A. 00052 * 00053 * LDAB (input) INTEGER 00054 * The leading dimension of the array AB. LDAB >= KD+1. 00055 * 00056 * INFO (output) INTEGER 00057 * = 0: successful exit 00058 * < 0: if INFO = -i, the i-th argument had an illegal value 00059 * > 0: if INFO = i, the leading minor of order i is not 00060 * positive definite, and the factorization could not be 00061 * completed. 00062 * 00063 * Further Details 00064 * =============== 00065 * 00066 * The band storage scheme is illustrated by the following example, when 00067 * N = 6, KD = 2, and UPLO = 'U': 00068 * 00069 * On entry: On exit: 00070 * 00071 * * * a13 a24 a35 a46 * * u13 u24 u35 u46 00072 * * a12 a23 a34 a45 a56 * u12 u23 u34 u45 u56 00073 * a11 a22 a33 a44 a55 a66 u11 u22 u33 u44 u55 u66 00074 * 00075 * Similarly, if UPLO = 'L' the format of A is as follows: 00076 * 00077 * On entry: On exit: 00078 * 00079 * a11 a22 a33 a44 a55 a66 l11 l22 l33 l44 l55 l66 00080 * a21 a32 a43 a54 a65 * l21 l32 l43 l54 l65 * 00081 * a31 a42 a53 a64 * * l31 l42 l53 l64 * * 00082 * 00083 * Array elements marked * are not used by the routine. 00084 * 00085 * Contributed by 00086 * Peter Mayes and Giuseppe Radicati, IBM ECSEC, Rome, March 23, 1989 00087 * 00088 * ===================================================================== 00089 * 00090 * .. Parameters .. 00091 DOUBLE PRECISION ONE, ZERO 00092 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 00093 INTEGER NBMAX, LDWORK 00094 PARAMETER ( NBMAX = 32, LDWORK = NBMAX+1 ) 00095 * .. 00096 * .. Local Scalars .. 00097 INTEGER I, I2, I3, IB, II, J, JJ, NB 00098 * .. 00099 * .. Local Arrays .. 00100 DOUBLE PRECISION WORK( LDWORK, NBMAX ) 00101 * .. 00102 * .. External Functions .. 00103 LOGICAL LSAME 00104 INTEGER ILAENV 00105 EXTERNAL LSAME, ILAENV 00106 * .. 00107 * .. External Subroutines .. 00108 EXTERNAL DGEMM, DPBTF2, DPOTF2, DSYRK, DTRSM, XERBLA 00109 * .. 00110 * .. Intrinsic Functions .. 00111 INTRINSIC MIN 00112 * .. 00113 * .. Executable Statements .. 00114 * 00115 * Test the input parameters. 00116 * 00117 INFO = 0 00118 IF( ( .NOT.LSAME( UPLO, 'U' ) ) .AND. 00119 $ ( .NOT.LSAME( UPLO, 'L' ) ) ) THEN 00120 INFO = -1 00121 ELSE IF( N.LT.0 ) THEN 00122 INFO = -2 00123 ELSE IF( KD.LT.0 ) THEN 00124 INFO = -3 00125 ELSE IF( LDAB.LT.KD+1 ) THEN 00126 INFO = -5 00127 END IF 00128 IF( INFO.NE.0 ) THEN 00129 CALL XERBLA( 'DPBTRF', -INFO ) 00130 RETURN 00131 END IF 00132 * 00133 * Quick return if possible 00134 * 00135 IF( N.EQ.0 ) 00136 $ RETURN 00137 * 00138 * Determine the block size for this environment 00139 * 00140 NB = ILAENV( 1, 'DPBTRF', UPLO, N, KD, -1, -1 ) 00141 * 00142 * The block size must not exceed the semi-bandwidth KD, and must not 00143 * exceed the limit set by the size of the local array WORK. 00144 * 00145 NB = MIN( NB, NBMAX ) 00146 * 00147 IF( NB.LE.1 .OR. NB.GT.KD ) THEN 00148 * 00149 * Use unblocked code 00150 * 00151 CALL DPBTF2( UPLO, N, KD, AB, LDAB, INFO ) 00152 ELSE 00153 * 00154 * Use blocked code 00155 * 00156 IF( LSAME( UPLO, 'U' ) ) THEN 00157 * 00158 * Compute the Cholesky factorization of a symmetric band 00159 * matrix, given the upper triangle of the matrix in band 00160 * storage. 00161 * 00162 * Zero the upper triangle of the work array. 00163 * 00164 DO 20 J = 1, NB 00165 DO 10 I = 1, J - 1 00166 WORK( I, J ) = ZERO 00167 10 CONTINUE 00168 20 CONTINUE 00169 * 00170 * Process the band matrix one diagonal block at a time. 00171 * 00172 DO 70 I = 1, N, NB 00173 IB = MIN( NB, N-I+1 ) 00174 * 00175 * Factorize the diagonal block 00176 * 00177 CALL DPOTF2( UPLO, IB, AB( KD+1, I ), LDAB-1, II ) 00178 IF( II.NE.0 ) THEN 00179 INFO = I + II - 1 00180 GO TO 150 00181 END IF 00182 IF( I+IB.LE.N ) THEN 00183 * 00184 * Update the relevant part of the trailing submatrix. 00185 * If A11 denotes the diagonal block which has just been 00186 * factorized, then we need to update the remaining 00187 * blocks in the diagram: 00188 * 00189 * A11 A12 A13 00190 * A22 A23 00191 * A33 00192 * 00193 * The numbers of rows and columns in the partitioning 00194 * are IB, I2, I3 respectively. The blocks A12, A22 and 00195 * A23 are empty if IB = KD. The upper triangle of A13 00196 * lies outside the band. 00197 * 00198 I2 = MIN( KD-IB, N-I-IB+1 ) 00199 I3 = MIN( IB, N-I-KD+1 ) 00200 * 00201 IF( I2.GT.0 ) THEN 00202 * 00203 * Update A12 00204 * 00205 CALL DTRSM( 'Left', 'Upper', 'Transpose', 00206 $ 'Non-unit', IB, I2, ONE, AB( KD+1, I ), 00207 $ LDAB-1, AB( KD+1-IB, I+IB ), LDAB-1 ) 00208 * 00209 * Update A22 00210 * 00211 CALL DSYRK( 'Upper', 'Transpose', I2, IB, -ONE, 00212 $ AB( KD+1-IB, I+IB ), LDAB-1, ONE, 00213 $ AB( KD+1, I+IB ), LDAB-1 ) 00214 END IF 00215 * 00216 IF( I3.GT.0 ) THEN 00217 * 00218 * Copy the lower triangle of A13 into the work array. 00219 * 00220 DO 40 JJ = 1, I3 00221 DO 30 II = JJ, IB 00222 WORK( II, JJ ) = AB( II-JJ+1, JJ+I+KD-1 ) 00223 30 CONTINUE 00224 40 CONTINUE 00225 * 00226 * Update A13 (in the work array). 00227 * 00228 CALL DTRSM( 'Left', 'Upper', 'Transpose', 00229 $ 'Non-unit', IB, I3, ONE, AB( KD+1, I ), 00230 $ LDAB-1, WORK, LDWORK ) 00231 * 00232 * Update A23 00233 * 00234 IF( I2.GT.0 ) 00235 $ CALL DGEMM( 'Transpose', 'No Transpose', I2, I3, 00236 $ IB, -ONE, AB( KD+1-IB, I+IB ), 00237 $ LDAB-1, WORK, LDWORK, ONE, 00238 $ AB( 1+IB, I+KD ), LDAB-1 ) 00239 * 00240 * Update A33 00241 * 00242 CALL DSYRK( 'Upper', 'Transpose', I3, IB, -ONE, 00243 $ WORK, LDWORK, ONE, AB( KD+1, I+KD ), 00244 $ LDAB-1 ) 00245 * 00246 * Copy the lower triangle of A13 back into place. 00247 * 00248 DO 60 JJ = 1, I3 00249 DO 50 II = JJ, IB 00250 AB( II-JJ+1, JJ+I+KD-1 ) = WORK( II, JJ ) 00251 50 CONTINUE 00252 60 CONTINUE 00253 END IF 00254 END IF 00255 70 CONTINUE 00256 ELSE 00257 * 00258 * Compute the Cholesky factorization of a symmetric band 00259 * matrix, given the lower triangle of the matrix in band 00260 * storage. 00261 * 00262 * Zero the lower triangle of the work array. 00263 * 00264 DO 90 J = 1, NB 00265 DO 80 I = J + 1, NB 00266 WORK( I, J ) = ZERO 00267 80 CONTINUE 00268 90 CONTINUE 00269 * 00270 * Process the band matrix one diagonal block at a time. 00271 * 00272 DO 140 I = 1, N, NB 00273 IB = MIN( NB, N-I+1 ) 00274 * 00275 * Factorize the diagonal block 00276 * 00277 CALL DPOTF2( UPLO, IB, AB( 1, I ), LDAB-1, II ) 00278 IF( II.NE.0 ) THEN 00279 INFO = I + II - 1 00280 GO TO 150 00281 END IF 00282 IF( I+IB.LE.N ) THEN 00283 * 00284 * Update the relevant part of the trailing submatrix. 00285 * If A11 denotes the diagonal block which has just been 00286 * factorized, then we need to update the remaining 00287 * blocks in the diagram: 00288 * 00289 * A11 00290 * A21 A22 00291 * A31 A32 A33 00292 * 00293 * The numbers of rows and columns in the partitioning 00294 * are IB, I2, I3 respectively. The blocks A21, A22 and 00295 * A32 are empty if IB = KD. The lower triangle of A31 00296 * lies outside the band. 00297 * 00298 I2 = MIN( KD-IB, N-I-IB+1 ) 00299 I3 = MIN( IB, N-I-KD+1 ) 00300 * 00301 IF( I2.GT.0 ) THEN 00302 * 00303 * Update A21 00304 * 00305 CALL DTRSM( 'Right', 'Lower', 'Transpose', 00306 $ 'Non-unit', I2, IB, ONE, AB( 1, I ), 00307 $ LDAB-1, AB( 1+IB, I ), LDAB-1 ) 00308 * 00309 * Update A22 00310 * 00311 CALL DSYRK( 'Lower', 'No Transpose', I2, IB, -ONE, 00312 $ AB( 1+IB, I ), LDAB-1, ONE, 00313 $ AB( 1, I+IB ), LDAB-1 ) 00314 END IF 00315 * 00316 IF( I3.GT.0 ) THEN 00317 * 00318 * Copy the upper triangle of A31 into the work array. 00319 * 00320 DO 110 JJ = 1, IB 00321 DO 100 II = 1, MIN( JJ, I3 ) 00322 WORK( II, JJ ) = AB( KD+1-JJ+II, JJ+I-1 ) 00323 100 CONTINUE 00324 110 CONTINUE 00325 * 00326 * Update A31 (in the work array). 00327 * 00328 CALL DTRSM( 'Right', 'Lower', 'Transpose', 00329 $ 'Non-unit', I3, IB, ONE, AB( 1, I ), 00330 $ LDAB-1, WORK, LDWORK ) 00331 * 00332 * Update A32 00333 * 00334 IF( I2.GT.0 ) 00335 $ CALL DGEMM( 'No transpose', 'Transpose', I3, I2, 00336 $ IB, -ONE, WORK, LDWORK, 00337 $ AB( 1+IB, I ), LDAB-1, ONE, 00338 $ AB( 1+KD-IB, I+IB ), LDAB-1 ) 00339 * 00340 * Update A33 00341 * 00342 CALL DSYRK( 'Lower', 'No Transpose', I3, IB, -ONE, 00343 $ WORK, LDWORK, ONE, AB( 1, I+KD ), 00344 $ LDAB-1 ) 00345 * 00346 * Copy the upper triangle of A31 back into place. 00347 * 00348 DO 130 JJ = 1, IB 00349 DO 120 II = 1, MIN( JJ, I3 ) 00350 AB( KD+1-JJ+II, JJ+I-1 ) = WORK( II, JJ ) 00351 120 CONTINUE 00352 130 CONTINUE 00353 END IF 00354 END IF 00355 140 CONTINUE 00356 END IF 00357 END IF 00358 RETURN 00359 * 00360 150 CONTINUE 00361 RETURN 00362 * 00363 * End of DPBTRF 00364 * 00365 END
1.7.3 
