|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE CPBTRF( 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 COMPLEX AB( LDAB, * ) 00014 * .. 00015 * 00016 * Purpose 00017 * ======= 00018 * 00019 * CPBTRF computes the Cholesky factorization of a complex Hermitian 00020 * positive definite band matrix A. 00021 * 00022 * The factorization has the form 00023 * A = U**H * U, if UPLO = 'U', or 00024 * A = L * L**H, 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) COMPLEX array, dimension (LDAB,N) 00042 * On entry, the upper or lower triangle of the Hermitian 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**H*U or A = L*L**H 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 REAL ONE, ZERO 00092 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 00093 COMPLEX CONE 00094 PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ) ) 00095 INTEGER NBMAX, LDWORK 00096 PARAMETER ( NBMAX = 32, LDWORK = NBMAX+1 ) 00097 * .. 00098 * .. Local Scalars .. 00099 INTEGER I, I2, I3, IB, II, J, JJ, NB 00100 * .. 00101 * .. Local Arrays .. 00102 COMPLEX WORK( LDWORK, NBMAX ) 00103 * .. 00104 * .. External Functions .. 00105 LOGICAL LSAME 00106 INTEGER ILAENV 00107 EXTERNAL LSAME, ILAENV 00108 * .. 00109 * .. External Subroutines .. 00110 EXTERNAL CGEMM, CHERK, CPBTF2, CPOTF2, CTRSM, XERBLA 00111 * .. 00112 * .. Intrinsic Functions .. 00113 INTRINSIC MIN 00114 * .. 00115 * .. Executable Statements .. 00116 * 00117 * Test the input parameters. 00118 * 00119 INFO = 0 00120 IF( ( .NOT.LSAME( UPLO, 'U' ) ) .AND. 00121 $ ( .NOT.LSAME( UPLO, 'L' ) ) ) THEN 00122 INFO = -1 00123 ELSE IF( N.LT.0 ) THEN 00124 INFO = -2 00125 ELSE IF( KD.LT.0 ) THEN 00126 INFO = -3 00127 ELSE IF( LDAB.LT.KD+1 ) THEN 00128 INFO = -5 00129 END IF 00130 IF( INFO.NE.0 ) THEN 00131 CALL XERBLA( 'CPBTRF', -INFO ) 00132 RETURN 00133 END IF 00134 * 00135 * Quick return if possible 00136 * 00137 IF( N.EQ.0 ) 00138 $ RETURN 00139 * 00140 * Determine the block size for this environment 00141 * 00142 NB = ILAENV( 1, 'CPBTRF', UPLO, N, KD, -1, -1 ) 00143 * 00144 * The block size must not exceed the semi-bandwidth KD, and must not 00145 * exceed the limit set by the size of the local array WORK. 00146 * 00147 NB = MIN( NB, NBMAX ) 00148 * 00149 IF( NB.LE.1 .OR. NB.GT.KD ) THEN 00150 * 00151 * Use unblocked code 00152 * 00153 CALL CPBTF2( UPLO, N, KD, AB, LDAB, INFO ) 00154 ELSE 00155 * 00156 * Use blocked code 00157 * 00158 IF( LSAME( UPLO, 'U' ) ) THEN 00159 * 00160 * Compute the Cholesky factorization of a Hermitian band 00161 * matrix, given the upper triangle of the matrix in band 00162 * storage. 00163 * 00164 * Zero the upper triangle of the work array. 00165 * 00166 DO 20 J = 1, NB 00167 DO 10 I = 1, J - 1 00168 WORK( I, J ) = ZERO 00169 10 CONTINUE 00170 20 CONTINUE 00171 * 00172 * Process the band matrix one diagonal block at a time. 00173 * 00174 DO 70 I = 1, N, NB 00175 IB = MIN( NB, N-I+1 ) 00176 * 00177 * Factorize the diagonal block 00178 * 00179 CALL CPOTF2( UPLO, IB, AB( KD+1, I ), LDAB-1, II ) 00180 IF( II.NE.0 ) THEN 00181 INFO = I + II - 1 00182 GO TO 150 00183 END IF 00184 IF( I+IB.LE.N ) THEN 00185 * 00186 * Update the relevant part of the trailing submatrix. 00187 * If A11 denotes the diagonal block which has just been 00188 * factorized, then we need to update the remaining 00189 * blocks in the diagram: 00190 * 00191 * A11 A12 A13 00192 * A22 A23 00193 * A33 00194 * 00195 * The numbers of rows and columns in the partitioning 00196 * are IB, I2, I3 respectively. The blocks A12, A22 and 00197 * A23 are empty if IB = KD. The upper triangle of A13 00198 * lies outside the band. 00199 * 00200 I2 = MIN( KD-IB, N-I-IB+1 ) 00201 I3 = MIN( IB, N-I-KD+1 ) 00202 * 00203 IF( I2.GT.0 ) THEN 00204 * 00205 * Update A12 00206 * 00207 CALL CTRSM( 'Left', 'Upper', 'Conjugate transpose', 00208 $ 'Non-unit', IB, I2, CONE, 00209 $ AB( KD+1, I ), LDAB-1, 00210 $ AB( KD+1-IB, I+IB ), LDAB-1 ) 00211 * 00212 * Update A22 00213 * 00214 CALL CHERK( 'Upper', 'Conjugate transpose', I2, IB, 00215 $ -ONE, AB( KD+1-IB, I+IB ), LDAB-1, ONE, 00216 $ AB( KD+1, I+IB ), LDAB-1 ) 00217 END IF 00218 * 00219 IF( I3.GT.0 ) THEN 00220 * 00221 * Copy the lower triangle of A13 into the work array. 00222 * 00223 DO 40 JJ = 1, I3 00224 DO 30 II = JJ, IB 00225 WORK( II, JJ ) = AB( II-JJ+1, JJ+I+KD-1 ) 00226 30 CONTINUE 00227 40 CONTINUE 00228 * 00229 * Update A13 (in the work array). 00230 * 00231 CALL CTRSM( 'Left', 'Upper', 'Conjugate transpose', 00232 $ 'Non-unit', IB, I3, CONE, 00233 $ AB( KD+1, I ), LDAB-1, WORK, LDWORK ) 00234 * 00235 * Update A23 00236 * 00237 IF( I2.GT.0 ) 00238 $ CALL CGEMM( 'Conjugate transpose', 00239 $ 'No transpose', I2, I3, IB, -CONE, 00240 $ AB( KD+1-IB, I+IB ), LDAB-1, WORK, 00241 $ LDWORK, CONE, AB( 1+IB, I+KD ), 00242 $ LDAB-1 ) 00243 * 00244 * Update A33 00245 * 00246 CALL CHERK( 'Upper', 'Conjugate transpose', I3, IB, 00247 $ -ONE, WORK, LDWORK, ONE, 00248 $ AB( KD+1, I+KD ), LDAB-1 ) 00249 * 00250 * Copy the lower triangle of A13 back into place. 00251 * 00252 DO 60 JJ = 1, I3 00253 DO 50 II = JJ, IB 00254 AB( II-JJ+1, JJ+I+KD-1 ) = WORK( II, JJ ) 00255 50 CONTINUE 00256 60 CONTINUE 00257 END IF 00258 END IF 00259 70 CONTINUE 00260 ELSE 00261 * 00262 * Compute the Cholesky factorization of a Hermitian band 00263 * matrix, given the lower triangle of the matrix in band 00264 * storage. 00265 * 00266 * Zero the lower triangle of the work array. 00267 * 00268 DO 90 J = 1, NB 00269 DO 80 I = J + 1, NB 00270 WORK( I, J ) = ZERO 00271 80 CONTINUE 00272 90 CONTINUE 00273 * 00274 * Process the band matrix one diagonal block at a time. 00275 * 00276 DO 140 I = 1, N, NB 00277 IB = MIN( NB, N-I+1 ) 00278 * 00279 * Factorize the diagonal block 00280 * 00281 CALL CPOTF2( UPLO, IB, AB( 1, I ), LDAB-1, II ) 00282 IF( II.NE.0 ) THEN 00283 INFO = I + II - 1 00284 GO TO 150 00285 END IF 00286 IF( I+IB.LE.N ) THEN 00287 * 00288 * Update the relevant part of the trailing submatrix. 00289 * If A11 denotes the diagonal block which has just been 00290 * factorized, then we need to update the remaining 00291 * blocks in the diagram: 00292 * 00293 * A11 00294 * A21 A22 00295 * A31 A32 A33 00296 * 00297 * The numbers of rows and columns in the partitioning 00298 * are IB, I2, I3 respectively. The blocks A21, A22 and 00299 * A32 are empty if IB = KD. The lower triangle of A31 00300 * lies outside the band. 00301 * 00302 I2 = MIN( KD-IB, N-I-IB+1 ) 00303 I3 = MIN( IB, N-I-KD+1 ) 00304 * 00305 IF( I2.GT.0 ) THEN 00306 * 00307 * Update A21 00308 * 00309 CALL CTRSM( 'Right', 'Lower', 00310 $ 'Conjugate transpose', 'Non-unit', I2, 00311 $ IB, CONE, AB( 1, I ), LDAB-1, 00312 $ AB( 1+IB, I ), LDAB-1 ) 00313 * 00314 * Update A22 00315 * 00316 CALL CHERK( 'Lower', 'No transpose', I2, IB, -ONE, 00317 $ AB( 1+IB, I ), LDAB-1, ONE, 00318 $ AB( 1, I+IB ), LDAB-1 ) 00319 END IF 00320 * 00321 IF( I3.GT.0 ) THEN 00322 * 00323 * Copy the upper triangle of A31 into the work array. 00324 * 00325 DO 110 JJ = 1, IB 00326 DO 100 II = 1, MIN( JJ, I3 ) 00327 WORK( II, JJ ) = AB( KD+1-JJ+II, JJ+I-1 ) 00328 100 CONTINUE 00329 110 CONTINUE 00330 * 00331 * Update A31 (in the work array). 00332 * 00333 CALL CTRSM( 'Right', 'Lower', 00334 $ 'Conjugate transpose', 'Non-unit', I3, 00335 $ IB, CONE, AB( 1, I ), LDAB-1, WORK, 00336 $ LDWORK ) 00337 * 00338 * Update A32 00339 * 00340 IF( I2.GT.0 ) 00341 $ CALL CGEMM( 'No transpose', 00342 $ 'Conjugate transpose', I3, I2, IB, 00343 $ -CONE, WORK, LDWORK, AB( 1+IB, I ), 00344 $ LDAB-1, CONE, AB( 1+KD-IB, I+IB ), 00345 $ LDAB-1 ) 00346 * 00347 * Update A33 00348 * 00349 CALL CHERK( 'Lower', 'No transpose', I3, IB, -ONE, 00350 $ WORK, LDWORK, ONE, AB( 1, I+KD ), 00351 $ LDAB-1 ) 00352 * 00353 * Copy the upper triangle of A31 back into place. 00354 * 00355 DO 130 JJ = 1, IB 00356 DO 120 II = 1, MIN( JJ, I3 ) 00357 AB( KD+1-JJ+II, JJ+I-1 ) = WORK( II, JJ ) 00358 120 CONTINUE 00359 130 CONTINUE 00360 END IF 00361 END IF 00362 140 CONTINUE 00363 END IF 00364 END IF 00365 RETURN 00366 * 00367 150 CONTINUE 00368 RETURN 00369 * 00370 * End of CPBTRF 00371 * 00372 END
1.7.3 
