zpbtf2.c
Go to the documentation of this file.
00001 /* zpbtf2.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 zpbtf2_(char *uplo, integer *n, integer *kd, 
00022         doublecomplex *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 zher_(char *, integer *, doublereal *, 
00036             doublecomplex *, integer *, doublecomplex *, integer *);
00037     extern logical lsame_(char *, char *);
00038     logical upper;
00039     extern /* Subroutine */ int xerbla_(char *, integer *), zdscal_(
00040             integer *, doublereal *, doublecomplex *, integer *), zlacgv_(
00041             integer *, doublecomplex *, 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 /*  ZPBTF2 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*16 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_("ZPBTF2", &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.) {
00188                 i__2 = *kd + 1 + j * ab_dim1;
00189                 ab[i__2].r = ajj, ab[i__2].i = 0.;
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.;
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                 d__1 = 1. / ajj;
00204                 zdscal_(&kn, &d__1, &ab[*kd + (j + 1) * ab_dim1], &kld);
00205                 zlacgv_(&kn, &ab[*kd + (j + 1) * ab_dim1], &kld);
00206                 zher_("Upper", &kn, &c_b8, &ab[*kd + (j + 1) * ab_dim1], &kld, 
00207                          &ab[*kd + 1 + (j + 1) * ab_dim1], &kld);
00208                 zlacgv_(&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.) {
00224                 i__2 = j * ab_dim1 + 1;
00225                 ab[i__2].r = ajj, ab[i__2].i = 0.;
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.;
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                 d__1 = 1. / ajj;
00240                 zdscal_(&kn, &d__1, &ab[j * ab_dim1 + 2], &c__1);
00241                 zher_("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 ZPBTF2 */
00254 
00255 } /* zpbtf2_ */


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