00001 /* dpbstf.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 integer c__1 = 1; 00019 static doublereal c_b9 = -1.; 00020 00021 /* Subroutine */ int dpbstf_(char *uplo, integer *n, integer *kd, doublereal * 00022 ab, integer *ldab, integer *info) 00023 { 00024 /* System generated locals */ 00025 integer ab_dim1, ab_offset, i__1, i__2, i__3; 00026 doublereal d__1; 00027 00028 /* Builtin functions */ 00029 double sqrt(doublereal); 00030 00031 /* Local variables */ 00032 integer j, m, km; 00033 doublereal ajj; 00034 integer kld; 00035 extern /* Subroutine */ int dsyr_(char *, integer *, doublereal *, 00036 doublereal *, integer *, doublereal *, integer *), dscal_( 00037 integer *, doublereal *, doublereal *, integer *); 00038 extern logical lsame_(char *, char *); 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 /* DPBSTF computes a split Cholesky factorization of a real */ 00056 /* symmetric positive definite band matrix A. */ 00057 00058 /* This routine is designed to be used in conjunction with DSBGST. */ 00059 00060 /* The factorization has the form A = S**T*S where S is a band matrix */ 00061 /* of the same bandwidth as A and the following structure: */ 00062 00063 /* S = ( U ) */ 00064 /* ( M L ) */ 00065 00066 /* where U is upper triangular of order m = (n+kd)/2, and L is lower */ 00067 /* triangular of order n-m. */ 00068 00069 /* Arguments */ 00070 /* ========= */ 00071 00072 /* UPLO (input) CHARACTER*1 */ 00073 /* = 'U': Upper triangle of A is stored; */ 00074 /* = 'L': Lower triangle of A is stored. */ 00075 00076 /* N (input) INTEGER */ 00077 /* The order of the matrix A. N >= 0. */ 00078 00079 /* KD (input) INTEGER */ 00080 /* The number of superdiagonals of the matrix A if UPLO = 'U', */ 00081 /* or the number of subdiagonals if UPLO = 'L'. KD >= 0. */ 00082 00083 /* AB (input/output) DOUBLE PRECISION array, dimension (LDAB,N) */ 00084 /* On entry, the upper or lower triangle of the symmetric 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 factor S from the split Cholesky */ 00092 /* factorization A = S**T*S. See Further Details. */ 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 = -i, the i-th argument had an illegal value */ 00100 /* > 0: if INFO = i, the factorization could not be completed, */ 00101 /* because the updated element a(i,i) was negative; the */ 00102 /* matrix A is not positive definite. */ 00103 00104 /* Further Details */ 00105 /* =============== */ 00106 00107 /* The band storage scheme is illustrated by the following example, when */ 00108 /* N = 7, KD = 2: */ 00109 00110 /* S = ( s11 s12 s13 ) */ 00111 /* ( s22 s23 s24 ) */ 00112 /* ( s33 s34 ) */ 00113 /* ( s44 ) */ 00114 /* ( s53 s54 s55 ) */ 00115 /* ( s64 s65 s66 ) */ 00116 /* ( s75 s76 s77 ) */ 00117 00118 /* If UPLO = 'U', the array AB holds: */ 00119 00120 /* on entry: on exit: */ 00121 00122 /* * * a13 a24 a35 a46 a57 * * s13 s24 s53 s64 s75 */ 00123 /* * a12 a23 a34 a45 a56 a67 * s12 s23 s34 s54 s65 s76 */ 00124 /* a11 a22 a33 a44 a55 a66 a77 s11 s22 s33 s44 s55 s66 s77 */ 00125 00126 /* If UPLO = 'L', the array AB holds: */ 00127 00128 /* on entry: on exit: */ 00129 00130 /* a11 a22 a33 a44 a55 a66 a77 s11 s22 s33 s44 s55 s66 s77 */ 00131 /* a21 a32 a43 a54 a65 a76 * s12 s23 s34 s54 s65 s76 * */ 00132 /* a31 a42 a53 a64 a64 * * s13 s24 s53 s64 s75 * * */ 00133 00134 /* Array elements marked * are not used by the routine. */ 00135 00136 /* ===================================================================== */ 00137 00138 /* .. Parameters .. */ 00139 /* .. */ 00140 /* .. Local Scalars .. */ 00141 /* .. */ 00142 /* .. External Functions .. */ 00143 /* .. */ 00144 /* .. External Subroutines .. */ 00145 /* .. */ 00146 /* .. Intrinsic Functions .. */ 00147 /* .. */ 00148 /* .. Executable Statements .. */ 00149 00150 /* Test the input parameters. */ 00151 00152 /* Parameter adjustments */ 00153 ab_dim1 = *ldab; 00154 ab_offset = 1 + ab_dim1; 00155 ab -= ab_offset; 00156 00157 /* Function Body */ 00158 *info = 0; 00159 upper = lsame_(uplo, "U"); 00160 if (! upper && ! lsame_(uplo, "L")) { 00161 *info = -1; 00162 } else if (*n < 0) { 00163 *info = -2; 00164 } else if (*kd < 0) { 00165 *info = -3; 00166 } else if (*ldab < *kd + 1) { 00167 *info = -5; 00168 } 00169 if (*info != 0) { 00170 i__1 = -(*info); 00171 xerbla_("DPBSTF", &i__1); 00172 return 0; 00173 } 00174 00175 /* Quick return if possible */ 00176 00177 if (*n == 0) { 00178 return 0; 00179 } 00180 00181 /* Computing MAX */ 00182 i__1 = 1, i__2 = *ldab - 1; 00183 kld = max(i__1,i__2); 00184 00185 /* Set the splitting point m. */ 00186 00187 m = (*n + *kd) / 2; 00188 00189 if (upper) { 00190 00191 /* Factorize A(m+1:n,m+1:n) as L**T*L, and update A(1:m,1:m). */ 00192 00193 i__1 = m + 1; 00194 for (j = *n; j >= i__1; --j) { 00195 00196 /* Compute s(j,j) and test for non-positive-definiteness. */ 00197 00198 ajj = ab[*kd + 1 + j * ab_dim1]; 00199 if (ajj <= 0.) { 00200 goto L50; 00201 } 00202 ajj = sqrt(ajj); 00203 ab[*kd + 1 + j * ab_dim1] = ajj; 00204 /* Computing MIN */ 00205 i__2 = j - 1; 00206 km = min(i__2,*kd); 00207 00208 /* Compute elements j-km:j-1 of the j-th column and update the */ 00209 /* the leading submatrix within the band. */ 00210 00211 d__1 = 1. / ajj; 00212 dscal_(&km, &d__1, &ab[*kd + 1 - km + j * ab_dim1], &c__1); 00213 dsyr_("Upper", &km, &c_b9, &ab[*kd + 1 - km + j * ab_dim1], &c__1, 00214 &ab[*kd + 1 + (j - km) * ab_dim1], &kld); 00215 /* L10: */ 00216 } 00217 00218 /* Factorize the updated submatrix A(1:m,1:m) as U**T*U. */ 00219 00220 i__1 = m; 00221 for (j = 1; j <= i__1; ++j) { 00222 00223 /* Compute s(j,j) and test for non-positive-definiteness. */ 00224 00225 ajj = ab[*kd + 1 + j * ab_dim1]; 00226 if (ajj <= 0.) { 00227 goto L50; 00228 } 00229 ajj = sqrt(ajj); 00230 ab[*kd + 1 + j * ab_dim1] = ajj; 00231 /* Computing MIN */ 00232 i__2 = *kd, i__3 = m - j; 00233 km = min(i__2,i__3); 00234 00235 /* Compute elements j+1:j+km of the j-th row and update the */ 00236 /* trailing submatrix within the band. */ 00237 00238 if (km > 0) { 00239 d__1 = 1. / ajj; 00240 dscal_(&km, &d__1, &ab[*kd + (j + 1) * ab_dim1], &kld); 00241 dsyr_("Upper", &km, &c_b9, &ab[*kd + (j + 1) * ab_dim1], &kld, 00242 &ab[*kd + 1 + (j + 1) * ab_dim1], &kld); 00243 } 00244 /* L20: */ 00245 } 00246 } else { 00247 00248 /* Factorize A(m+1:n,m+1:n) as L**T*L, and update A(1:m,1:m). */ 00249 00250 i__1 = m + 1; 00251 for (j = *n; j >= i__1; --j) { 00252 00253 /* Compute s(j,j) and test for non-positive-definiteness. */ 00254 00255 ajj = ab[j * ab_dim1 + 1]; 00256 if (ajj <= 0.) { 00257 goto L50; 00258 } 00259 ajj = sqrt(ajj); 00260 ab[j * ab_dim1 + 1] = ajj; 00261 /* Computing MIN */ 00262 i__2 = j - 1; 00263 km = min(i__2,*kd); 00264 00265 /* Compute elements j-km:j-1 of the j-th row and update the */ 00266 /* trailing submatrix within the band. */ 00267 00268 d__1 = 1. / ajj; 00269 dscal_(&km, &d__1, &ab[km + 1 + (j - km) * ab_dim1], &kld); 00270 dsyr_("Lower", &km, &c_b9, &ab[km + 1 + (j - km) * ab_dim1], &kld, 00271 &ab[(j - km) * ab_dim1 + 1], &kld); 00272 /* L30: */ 00273 } 00274 00275 /* Factorize the updated submatrix A(1:m,1:m) as U**T*U. */ 00276 00277 i__1 = m; 00278 for (j = 1; j <= i__1; ++j) { 00279 00280 /* Compute s(j,j) and test for non-positive-definiteness. */ 00281 00282 ajj = ab[j * ab_dim1 + 1]; 00283 if (ajj <= 0.) { 00284 goto L50; 00285 } 00286 ajj = sqrt(ajj); 00287 ab[j * ab_dim1 + 1] = ajj; 00288 /* Computing MIN */ 00289 i__2 = *kd, i__3 = m - j; 00290 km = min(i__2,i__3); 00291 00292 /* Compute elements j+1:j+km of the j-th column and update the */ 00293 /* trailing submatrix within the band. */ 00294 00295 if (km > 0) { 00296 d__1 = 1. / ajj; 00297 dscal_(&km, &d__1, &ab[j * ab_dim1 + 2], &c__1); 00298 dsyr_("Lower", &km, &c_b9, &ab[j * ab_dim1 + 2], &c__1, &ab[( 00299 j + 1) * ab_dim1 + 1], &kld); 00300 } 00301 /* L40: */ 00302 } 00303 } 00304 return 0; 00305 00306 L50: 00307 *info = j; 00308 return 0; 00309 00310 /* End of DPBSTF */ 00311 00312 } /* dpbstf_ */