00001 /* cpbtf2.f -- translated by f2c (version 20061008). 00002 You must link the resulting object file with libf2c: 00003 on Microsoft Windows system, link with libf2c.lib; 00004 on Linux or Unix systems, link with .../path/to/libf2c.a -lm 00005 or, if you install libf2c.a in a standard place, with -lf2c -lm 00006 -- in that order, at the end of the command line, as in 00007 cc *.o -lf2c -lm 00008 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., 00009 00010 http://www.netlib.org/f2c/libf2c.zip 00011 */ 00012 00013 #include "f2c.h" 00014 #include "blaswrap.h" 00015 00016 /* Table of constant values */ 00017 00018 static real c_b8 = -1.f; 00019 static integer c__1 = 1; 00020 00021 /* Subroutine */ int cpbtf2_(char *uplo, integer *n, integer *kd, complex *ab, 00022 integer *ldab, integer *info) 00023 { 00024 /* System generated locals */ 00025 integer ab_dim1, ab_offset, i__1, i__2, i__3; 00026 real r__1; 00027 00028 /* Builtin functions */ 00029 double sqrt(doublereal); 00030 00031 /* Local variables */ 00032 integer j, kn; 00033 real ajj; 00034 integer kld; 00035 extern /* Subroutine */ int cher_(char *, integer *, real *, complex *, 00036 integer *, complex *, integer *); 00037 extern logical lsame_(char *, char *); 00038 logical upper; 00039 extern /* Subroutine */ int clacgv_(integer *, complex *, integer *), 00040 csscal_(integer *, real *, complex *, integer *), xerbla_(char *, 00041 integer *); 00042 00043 00044 /* -- LAPACK routine (version 3.2) -- */ 00045 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00046 /* November 2006 */ 00047 00048 /* .. Scalar Arguments .. */ 00049 /* .. */ 00050 /* .. Array Arguments .. */ 00051 /* .. */ 00052 00053 /* Purpose */ 00054 /* ======= */ 00055 00056 /* CPBTF2 computes the Cholesky factorization of a complex Hermitian */ 00057 /* positive definite band matrix A. */ 00058 00059 /* The factorization has the form */ 00060 /* A = U' * U , if UPLO = 'U', or */ 00061 /* A = L * L', if UPLO = 'L', */ 00062 /* where U is an upper triangular matrix, U' is the conjugate transpose */ 00063 /* of U, and L is lower triangular. */ 00064 00065 /* This is the unblocked version of the algorithm, calling Level 2 BLAS. */ 00066 00067 /* Arguments */ 00068 /* ========= */ 00069 00070 /* UPLO (input) CHARACTER*1 */ 00071 /* Specifies whether the upper or lower triangular part of the */ 00072 /* Hermitian matrix A is stored: */ 00073 /* = 'U': Upper triangular */ 00074 /* = 'L': Lower triangular */ 00075 00076 /* N (input) INTEGER */ 00077 /* The order of the matrix A. N >= 0. */ 00078 00079 /* KD (input) INTEGER */ 00080 /* The number of super-diagonals of the matrix A if UPLO = 'U', */ 00081 /* or the number of sub-diagonals if UPLO = 'L'. KD >= 0. */ 00082 00083 /* AB (input/output) COMPLEX array, dimension (LDAB,N) */ 00084 /* On entry, the upper or lower triangle of the Hermitian band */ 00085 /* matrix A, stored in the first KD+1 rows of the array. The */ 00086 /* j-th column of A is stored in the j-th column of the array AB */ 00087 /* as follows: */ 00088 /* if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j; */ 00089 /* if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd). */ 00090 00091 /* On exit, if INFO = 0, the triangular factor U or L from the */ 00092 /* Cholesky factorization A = U'*U or A = L*L' of the band */ 00093 /* matrix A, in the same storage format as A. */ 00094 00095 /* LDAB (input) INTEGER */ 00096 /* The leading dimension of the array AB. LDAB >= KD+1. */ 00097 00098 /* INFO (output) INTEGER */ 00099 /* = 0: successful exit */ 00100 /* < 0: if INFO = -k, the k-th argument had an illegal value */ 00101 /* > 0: if INFO = k, the leading minor of order k is not */ 00102 /* positive definite, and the factorization could not be */ 00103 /* completed. */ 00104 00105 /* Further Details */ 00106 /* =============== */ 00107 00108 /* The band storage scheme is illustrated by the following example, when */ 00109 /* N = 6, KD = 2, and UPLO = 'U': */ 00110 00111 /* On entry: On exit: */ 00112 00113 /* * * a13 a24 a35 a46 * * u13 u24 u35 u46 */ 00114 /* * a12 a23 a34 a45 a56 * u12 u23 u34 u45 u56 */ 00115 /* a11 a22 a33 a44 a55 a66 u11 u22 u33 u44 u55 u66 */ 00116 00117 /* Similarly, if UPLO = 'L' the format of A is as follows: */ 00118 00119 /* On entry: On exit: */ 00120 00121 /* a11 a22 a33 a44 a55 a66 l11 l22 l33 l44 l55 l66 */ 00122 /* a21 a32 a43 a54 a65 * l21 l32 l43 l54 l65 * */ 00123 /* a31 a42 a53 a64 * * l31 l42 l53 l64 * * */ 00124 00125 /* Array elements marked * are not used by the routine. */ 00126 00127 /* ===================================================================== */ 00128 00129 /* .. Parameters .. */ 00130 /* .. */ 00131 /* .. Local Scalars .. */ 00132 /* .. */ 00133 /* .. External Functions .. */ 00134 /* .. */ 00135 /* .. External Subroutines .. */ 00136 /* .. */ 00137 /* .. Intrinsic Functions .. */ 00138 /* .. */ 00139 /* .. Executable Statements .. */ 00140 00141 /* Test the input parameters. */ 00142 00143 /* Parameter adjustments */ 00144 ab_dim1 = *ldab; 00145 ab_offset = 1 + ab_dim1; 00146 ab -= ab_offset; 00147 00148 /* Function Body */ 00149 *info = 0; 00150 upper = lsame_(uplo, "U"); 00151 if (! upper && ! lsame_(uplo, "L")) { 00152 *info = -1; 00153 } else if (*n < 0) { 00154 *info = -2; 00155 } else if (*kd < 0) { 00156 *info = -3; 00157 } else if (*ldab < *kd + 1) { 00158 *info = -5; 00159 } 00160 if (*info != 0) { 00161 i__1 = -(*info); 00162 xerbla_("CPBTF2", &i__1); 00163 return 0; 00164 } 00165 00166 /* Quick return if possible */ 00167 00168 if (*n == 0) { 00169 return 0; 00170 } 00171 00172 /* Computing MAX */ 00173 i__1 = 1, i__2 = *ldab - 1; 00174 kld = max(i__1,i__2); 00175 00176 if (upper) { 00177 00178 /* Compute the Cholesky factorization A = U'*U. */ 00179 00180 i__1 = *n; 00181 for (j = 1; j <= i__1; ++j) { 00182 00183 /* Compute U(J,J) and test for non-positive-definiteness. */ 00184 00185 i__2 = *kd + 1 + j * ab_dim1; 00186 ajj = ab[i__2].r; 00187 if (ajj <= 0.f) { 00188 i__2 = *kd + 1 + j * ab_dim1; 00189 ab[i__2].r = ajj, ab[i__2].i = 0.f; 00190 goto L30; 00191 } 00192 ajj = sqrt(ajj); 00193 i__2 = *kd + 1 + j * ab_dim1; 00194 ab[i__2].r = ajj, ab[i__2].i = 0.f; 00195 00196 /* Compute elements J+1:J+KN of row J and update the */ 00197 /* trailing submatrix within the band. */ 00198 00199 /* Computing MIN */ 00200 i__2 = *kd, i__3 = *n - j; 00201 kn = min(i__2,i__3); 00202 if (kn > 0) { 00203 r__1 = 1.f / ajj; 00204 csscal_(&kn, &r__1, &ab[*kd + (j + 1) * ab_dim1], &kld); 00205 clacgv_(&kn, &ab[*kd + (j + 1) * ab_dim1], &kld); 00206 cher_("Upper", &kn, &c_b8, &ab[*kd + (j + 1) * ab_dim1], &kld, 00207 &ab[*kd + 1 + (j + 1) * ab_dim1], &kld); 00208 clacgv_(&kn, &ab[*kd + (j + 1) * ab_dim1], &kld); 00209 } 00210 /* L10: */ 00211 } 00212 } else { 00213 00214 /* Compute the Cholesky factorization A = L*L'. */ 00215 00216 i__1 = *n; 00217 for (j = 1; j <= i__1; ++j) { 00218 00219 /* Compute L(J,J) and test for non-positive-definiteness. */ 00220 00221 i__2 = j * ab_dim1 + 1; 00222 ajj = ab[i__2].r; 00223 if (ajj <= 0.f) { 00224 i__2 = j * ab_dim1 + 1; 00225 ab[i__2].r = ajj, ab[i__2].i = 0.f; 00226 goto L30; 00227 } 00228 ajj = sqrt(ajj); 00229 i__2 = j * ab_dim1 + 1; 00230 ab[i__2].r = ajj, ab[i__2].i = 0.f; 00231 00232 /* Compute elements J+1:J+KN of column J and update the */ 00233 /* trailing submatrix within the band. */ 00234 00235 /* Computing MIN */ 00236 i__2 = *kd, i__3 = *n - j; 00237 kn = min(i__2,i__3); 00238 if (kn > 0) { 00239 r__1 = 1.f / ajj; 00240 csscal_(&kn, &r__1, &ab[j * ab_dim1 + 2], &c__1); 00241 cher_("Lower", &kn, &c_b8, &ab[j * ab_dim1 + 2], &c__1, &ab[( 00242 j + 1) * ab_dim1 + 1], &kld); 00243 } 00244 /* L20: */ 00245 } 00246 } 00247 return 0; 00248 00249 L30: 00250 *info = j; 00251 return 0; 00252 00253 /* End of CPBTF2 */ 00254 00255 } /* cpbtf2_ */