dpbtf2.c
Go to the documentation of this file.
00001 /* dpbtf2.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 doublereal c_b8 = -1.;
00019 static integer c__1 = 1;
00020 
00021 /* Subroutine */ int dpbtf2_(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, kn;
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 /*  DPBTF2 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) DOUBLE PRECISION 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_("DPBTF2", &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.) {
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                 d__1 = 1. / ajj;
00199                 dscal_(&kn, &d__1, &ab[*kd + (j + 1) * ab_dim1], &kld);
00200                 dsyr_("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.) {
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                 d__1 = 1. / ajj;
00229                 dscal_(&kn, &d__1, &ab[j * ab_dim1 + 2], &c__1);
00230                 dsyr_("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 DPBTF2 */
00243 
00244 } /* dpbtf2_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:55:47