spbstf.c
Go to the documentation of this file.
00001 /* spbstf.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 real c_b9 = -1.f;
00020 
00021 /* Subroutine */ int spbstf_(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, m, km;
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 /*  SPBSTF 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 SSBGST. */
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) REAL 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_("SPBSTF", &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.f) {
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             r__1 = 1.f / ajj;
00212             sscal_(&km, &r__1, &ab[*kd + 1 - km + j * ab_dim1], &c__1);
00213             ssyr_("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.f) {
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                 r__1 = 1.f / ajj;
00240                 sscal_(&km, &r__1, &ab[*kd + (j + 1) * ab_dim1], &kld);
00241                 ssyr_("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.f) {
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             r__1 = 1.f / ajj;
00269             sscal_(&km, &r__1, &ab[km + 1 + (j - km) * ab_dim1], &kld);
00270             ssyr_("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.f) {
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                 r__1 = 1.f / ajj;
00297                 sscal_(&km, &r__1, &ab[j * ab_dim1 + 2], &c__1);
00298                 ssyr_("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 SPBSTF */
00311 
00312 } /* spbstf_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:56:12