00001 /* spbtf2.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 spbtf2_(char *uplo, integer *n, integer *kd, real *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 ssyr_(char *, integer *, real *, real *, 00036 integer *, real *, integer *); 00037 extern logical lsame_(char *, char *); 00038 extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *); 00039 logical upper; 00040 extern /* Subroutine */ int xerbla_(char *, integer *); 00041 00042 00043 /* -- LAPACK routine (version 3.2) -- */ 00044 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00045 /* November 2006 */ 00046 00047 /* .. Scalar Arguments .. */ 00048 /* .. */ 00049 /* .. Array Arguments .. */ 00050 /* .. */ 00051 00052 /* Purpose */ 00053 /* ======= */ 00054 00055 /* SPBTF2 computes the Cholesky factorization of a real symmetric */ 00056 /* positive definite band matrix A. */ 00057 00058 /* The factorization has the form */ 00059 /* A = U' * U , if UPLO = 'U', or */ 00060 /* A = L * L', if UPLO = 'L', */ 00061 /* where U is an upper triangular matrix, U' is the transpose of U, and */ 00062 /* L is lower triangular. */ 00063 00064 /* This is the unblocked version of the algorithm, calling Level 2 BLAS. */ 00065 00066 /* Arguments */ 00067 /* ========= */ 00068 00069 /* UPLO (input) CHARACTER*1 */ 00070 /* Specifies whether the upper or lower triangular part of the */ 00071 /* symmetric matrix A is stored: */ 00072 /* = 'U': Upper triangular */ 00073 /* = 'L': Lower triangular */ 00074 00075 /* N (input) INTEGER */ 00076 /* The order of the matrix A. N >= 0. */ 00077 00078 /* KD (input) INTEGER */ 00079 /* The number of super-diagonals of the matrix A if UPLO = 'U', */ 00080 /* or the number of sub-diagonals if UPLO = 'L'. KD >= 0. */ 00081 00082 /* AB (input/output) REAL array, dimension (LDAB,N) */ 00083 /* On entry, the upper or lower triangle of the symmetric band */ 00084 /* matrix A, stored in the first KD+1 rows of the array. The */ 00085 /* j-th column of A is stored in the j-th column of the array AB */ 00086 /* as follows: */ 00087 /* if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j; */ 00088 /* if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd). */ 00089 00090 /* On exit, if INFO = 0, the triangular factor U or L from the */ 00091 /* Cholesky factorization A = U'*U or A = L*L' of the band */ 00092 /* matrix A, in the same storage format as A. */ 00093 00094 /* LDAB (input) INTEGER */ 00095 /* The leading dimension of the array AB. LDAB >= KD+1. */ 00096 00097 /* INFO (output) INTEGER */ 00098 /* = 0: successful exit */ 00099 /* < 0: if INFO = -k, the k-th argument had an illegal value */ 00100 /* > 0: if INFO = k, the leading minor of order k is not */ 00101 /* positive definite, and the factorization could not be */ 00102 /* completed. */ 00103 00104 /* Further Details */ 00105 /* =============== */ 00106 00107 /* The band storage scheme is illustrated by the following example, when */ 00108 /* N = 6, KD = 2, and UPLO = 'U': */ 00109 00110 /* On entry: On exit: */ 00111 00112 /* * * a13 a24 a35 a46 * * u13 u24 u35 u46 */ 00113 /* * a12 a23 a34 a45 a56 * u12 u23 u34 u45 u56 */ 00114 /* a11 a22 a33 a44 a55 a66 u11 u22 u33 u44 u55 u66 */ 00115 00116 /* Similarly, if UPLO = 'L' the format of A is as follows: */ 00117 00118 /* On entry: On exit: */ 00119 00120 /* a11 a22 a33 a44 a55 a66 l11 l22 l33 l44 l55 l66 */ 00121 /* a21 a32 a43 a54 a65 * l21 l32 l43 l54 l65 * */ 00122 /* a31 a42 a53 a64 * * l31 l42 l53 l64 * * */ 00123 00124 /* Array elements marked * are not used by the routine. */ 00125 00126 /* ===================================================================== */ 00127 00128 /* .. Parameters .. */ 00129 /* .. */ 00130 /* .. Local Scalars .. */ 00131 /* .. */ 00132 /* .. External Functions .. */ 00133 /* .. */ 00134 /* .. External Subroutines .. */ 00135 /* .. */ 00136 /* .. Intrinsic Functions .. */ 00137 /* .. */ 00138 /* .. Executable Statements .. */ 00139 00140 /* Test the input parameters. */ 00141 00142 /* Parameter adjustments */ 00143 ab_dim1 = *ldab; 00144 ab_offset = 1 + ab_dim1; 00145 ab -= ab_offset; 00146 00147 /* Function Body */ 00148 *info = 0; 00149 upper = lsame_(uplo, "U"); 00150 if (! upper && ! lsame_(uplo, "L")) { 00151 *info = -1; 00152 } else if (*n < 0) { 00153 *info = -2; 00154 } else if (*kd < 0) { 00155 *info = -3; 00156 } else if (*ldab < *kd + 1) { 00157 *info = -5; 00158 } 00159 if (*info != 0) { 00160 i__1 = -(*info); 00161 xerbla_("SPBTF2", &i__1); 00162 return 0; 00163 } 00164 00165 /* Quick return if possible */ 00166 00167 if (*n == 0) { 00168 return 0; 00169 } 00170 00171 /* Computing MAX */ 00172 i__1 = 1, i__2 = *ldab - 1; 00173 kld = max(i__1,i__2); 00174 00175 if (upper) { 00176 00177 /* Compute the Cholesky factorization A = U'*U. */ 00178 00179 i__1 = *n; 00180 for (j = 1; j <= i__1; ++j) { 00181 00182 /* Compute U(J,J) and test for non-positive-definiteness. */ 00183 00184 ajj = ab[*kd + 1 + j * ab_dim1]; 00185 if (ajj <= 0.f) { 00186 goto L30; 00187 } 00188 ajj = sqrt(ajj); 00189 ab[*kd + 1 + j * ab_dim1] = ajj; 00190 00191 /* Compute elements J+1:J+KN of row J and update the */ 00192 /* trailing submatrix within the band. */ 00193 00194 /* Computing MIN */ 00195 i__2 = *kd, i__3 = *n - j; 00196 kn = min(i__2,i__3); 00197 if (kn > 0) { 00198 r__1 = 1.f / ajj; 00199 sscal_(&kn, &r__1, &ab[*kd + (j + 1) * ab_dim1], &kld); 00200 ssyr_("Upper", &kn, &c_b8, &ab[*kd + (j + 1) * ab_dim1], &kld, 00201 &ab[*kd + 1 + (j + 1) * ab_dim1], &kld); 00202 } 00203 /* L10: */ 00204 } 00205 } else { 00206 00207 /* Compute the Cholesky factorization A = L*L'. */ 00208 00209 i__1 = *n; 00210 for (j = 1; j <= i__1; ++j) { 00211 00212 /* Compute L(J,J) and test for non-positive-definiteness. */ 00213 00214 ajj = ab[j * ab_dim1 + 1]; 00215 if (ajj <= 0.f) { 00216 goto L30; 00217 } 00218 ajj = sqrt(ajj); 00219 ab[j * ab_dim1 + 1] = ajj; 00220 00221 /* Compute elements J+1:J+KN of column J and update the */ 00222 /* trailing submatrix within the band. */ 00223 00224 /* Computing MIN */ 00225 i__2 = *kd, i__3 = *n - j; 00226 kn = min(i__2,i__3); 00227 if (kn > 0) { 00228 r__1 = 1.f / ajj; 00229 sscal_(&kn, &r__1, &ab[j * ab_dim1 + 2], &c__1); 00230 ssyr_("Lower", &kn, &c_b8, &ab[j * ab_dim1 + 2], &c__1, &ab[( 00231 j + 1) * ab_dim1 + 1], &kld); 00232 } 00233 /* L20: */ 00234 } 00235 } 00236 return 0; 00237 00238 L30: 00239 *info = j; 00240 return 0; 00241 00242 /* End of SPBTF2 */ 00243 00244 } /* spbtf2_ */