swiftnav: cpbstf.c Source File

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