 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE ZPBTF2( UPLO, N, KD, AB, LDAB, INFO ) 00002 * 00003 * -- LAPACK routine (version 3.3.1) -- 00004 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00005 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00006 * -- April 2011 -- 00007 * 00008 * .. Scalar Arguments .. 00009 CHARACTER UPLO 00010 INTEGER INFO, KD, LDAB, N 00011 * .. 00012 * .. Array Arguments .. 00013 COMPLEX*16 AB( LDAB, * ) 00014 * .. 00015 * 00016 * Purpose 00017 * ======= 00018 * 00019 * ZPBTF2 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, U**H is the conjugate transpose 00026 * of U, and L is lower triangular. 00027 * 00028 * This is the unblocked version of the algorithm, calling Level 2 BLAS. 00029 * 00030 * Arguments 00031 * ========= 00032 * 00033 * UPLO (input) CHARACTER*1 00034 * Specifies whether the upper or lower triangular part of the 00035 * Hermitian matrix A is stored: 00036 * = 'U': Upper triangular 00037 * = 'L': Lower triangular 00038 * 00039 * N (input) INTEGER 00040 * The order of the matrix A. N >= 0. 00041 * 00042 * KD (input) INTEGER 00043 * The number of super-diagonals of the matrix A if UPLO = 'U', 00044 * or the number of sub-diagonals if UPLO = 'L'. KD >= 0. 00045 * 00046 * AB (input/output) COMPLEX*16 array, dimension (LDAB,N) 00047 * On entry, the upper or lower triangle of the Hermitian band 00048 * matrix A, stored in the first KD+1 rows of the array. The 00049 * j-th column of A is stored in the j-th column of the array AB 00050 * as follows: 00051 * if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j; 00052 * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd). 00053 * 00054 * On exit, if INFO = 0, the triangular factor U or L from the 00055 * Cholesky factorization A = U**H *U or A = L*L**H of the band 00056 * matrix A, in the same storage format as A. 00057 * 00058 * LDAB (input) INTEGER 00059 * The leading dimension of the array AB. LDAB >= KD+1. 00060 * 00061 * INFO (output) INTEGER 00062 * = 0: successful exit 00063 * < 0: if INFO = -k, the k-th argument had an illegal value 00064 * > 0: if INFO = k, the leading minor of order k is not 00065 * positive definite, and the factorization could not be 00066 * completed. 00067 * 00068 * Further Details 00069 * =============== 00070 * 00071 * The band storage scheme is illustrated by the following example, when 00072 * N = 6, KD = 2, and UPLO = 'U': 00073 * 00074 * On entry: On exit: 00075 * 00076 * * * a13 a24 a35 a46 * * u13 u24 u35 u46 00077 * * a12 a23 a34 a45 a56 * u12 u23 u34 u45 u56 00078 * a11 a22 a33 a44 a55 a66 u11 u22 u33 u44 u55 u66 00079 * 00080 * Similarly, if UPLO = 'L' the format of A is as follows: 00081 * 00082 * On entry: On exit: 00083 * 00084 * a11 a22 a33 a44 a55 a66 l11 l22 l33 l44 l55 l66 00085 * a21 a32 a43 a54 a65 * l21 l32 l43 l54 l65 * 00086 * a31 a42 a53 a64 * * l31 l42 l53 l64 * * 00087 * 00088 * Array elements marked * are not used by the routine. 00089 * 00090 * ===================================================================== 00091 * 00092 * .. Parameters .. 00093 DOUBLE PRECISION ONE, ZERO 00094 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 00095 * .. 00096 * .. Local Scalars .. 00097 LOGICAL UPPER 00098 INTEGER J, KLD, KN 00099 DOUBLE PRECISION AJJ 00100 * .. 00101 * .. External Functions .. 00102 LOGICAL LSAME 00103 EXTERNAL LSAME 00104 * .. 00105 * .. External Subroutines .. 00106 EXTERNAL XERBLA, ZDSCAL, ZHER, ZLACGV 00107 * .. 00108 * .. Intrinsic Functions .. 00109 INTRINSIC DBLE, MAX, MIN, SQRT 00110 * .. 00111 * .. Executable Statements .. 00112 * 00113 * Test the input parameters. 00114 * 00115 INFO = 0 00116 UPPER = LSAME( UPLO, 'U' ) 00117 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00118 INFO = -1 00119 ELSE IF( N.LT.0 ) THEN 00120 INFO = -2 00121 ELSE IF( KD.LT.0 ) THEN 00122 INFO = -3 00123 ELSE IF( LDAB.LT.KD+1 ) THEN 00124 INFO = -5 00125 END IF 00126 IF( INFO.NE.0 ) THEN 00127 CALL XERBLA( 'ZPBTF2', -INFO ) 00128 RETURN 00129 END IF 00130 * 00131 * Quick return if possible 00132 * 00133 IF( N.EQ.0 ) 00134 $ RETURN 00135 * 00136 KLD = MAX( 1, LDAB-1 ) 00137 * 00138 IF( UPPER ) THEN 00139 * 00140 * Compute the Cholesky factorization A = U**H * U. 00141 * 00142 DO 10 J = 1, N 00143 * 00144 * Compute U(J,J) and test for non-positive-definiteness. 00145 * 00146 AJJ = DBLE( AB( KD+1, J ) ) 00147 IF( AJJ.LE.ZERO ) THEN 00148 AB( KD+1, J ) = AJJ 00149 GO TO 30 00150 END IF 00151 AJJ = SQRT( AJJ ) 00152 AB( KD+1, J ) = AJJ 00153 * 00154 * Compute elements J+1:J+KN of row J and update the 00155 * trailing submatrix within the band. 00156 * 00157 KN = MIN( KD, N-J ) 00158 IF( KN.GT.0 ) THEN 00159 CALL ZDSCAL( KN, ONE / AJJ, AB( KD, J+1 ), KLD ) 00160 CALL ZLACGV( KN, AB( KD, J+1 ), KLD ) 00161 CALL ZHER( 'Upper', KN, -ONE, AB( KD, J+1 ), KLD, 00162 $ AB( KD+1, J+1 ), KLD ) 00163 CALL ZLACGV( KN, AB( KD, J+1 ), KLD ) 00164 END IF 00165 10 CONTINUE 00166 ELSE 00167 * 00168 * Compute the Cholesky factorization A = L*L**H. 00169 * 00170 DO 20 J = 1, N 00171 * 00172 * Compute L(J,J) and test for non-positive-definiteness. 00173 * 00174 AJJ = DBLE( AB( 1, J ) ) 00175 IF( AJJ.LE.ZERO ) THEN 00176 AB( 1, J ) = AJJ 00177 GO TO 30 00178 END IF 00179 AJJ = SQRT( AJJ ) 00180 AB( 1, J ) = AJJ 00181 * 00182 * Compute elements J+1:J+KN of column J and update the 00183 * trailing submatrix within the band. 00184 * 00185 KN = MIN( KD, N-J ) 00186 IF( KN.GT.0 ) THEN 00187 CALL ZDSCAL( KN, ONE / AJJ, AB( 2, J ), 1 ) 00188 CALL ZHER( 'Lower', KN, -ONE, AB( 2, J ), 1, 00189 $ AB( 1, J+1 ), KLD ) 00190 END IF 00191 20 CONTINUE 00192 END IF 00193 RETURN 00194 * 00195 30 CONTINUE 00196 INFO = J 00197 RETURN 00198 * 00199 * End of ZPBTF2 00200 * 00201 END
1.7.3 
