swiftnav: chptrs.c Source File

Go to the documentation of this file.
00001 /* chptrs.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 complex c_b1 = {1.f,0.f};
00019 static integer c__1 = 1;
00020 
00021 /* Subroutine */ int chptrs_(char *uplo, integer *n, integer *nrhs, complex *
00022         ap, integer *ipiv, complex *b, integer *ldb, integer *info)
00023 {
00024     /* System generated locals */
00025     integer b_dim1, b_offset, i__1, i__2;
00026     complex q__1, q__2, q__3;
00027 
00028     /* Builtin functions */
00029     void c_div(complex *, complex *, complex *), r_cnjg(complex *, complex *);
00030 
00031     /* Local variables */
00032     integer j, k;
00033     real s;
00034     complex ak, bk;
00035     integer kc, kp;
00036     complex akm1, bkm1, akm1k;
00037     extern logical lsame_(char *, char *);
00038     complex denom;
00039     extern /* Subroutine */ int cgemv_(char *, integer *, integer *, complex *
00040 , complex *, integer *, complex *, integer *, complex *, complex *
00041 , integer *), cgeru_(integer *, integer *, complex *, 
00042             complex *, integer *, complex *, integer *, complex *, integer *),
00043              cswap_(integer *, complex *, integer *, complex *, integer *);
00044     logical upper;
00045     extern /* Subroutine */ int clacgv_(integer *, complex *, integer *), 
00046             csscal_(integer *, real *, complex *, integer *), xerbla_(char *, 
00047             integer *);
00048 
00049 
00050 /*  -- LAPACK routine (version 3.2) -- */
00051 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00052 /*     November 2006 */
00053 
00054 /*     .. Scalar Arguments .. */
00055 /*     .. */
00056 /*     .. Array Arguments .. */
00057 /*     .. */
00058 
00059 /*  Purpose */
00060 /*  ======= */
00061 
00062 /*  CHPTRS solves a system of linear equations A*X = B with a complex */
00063 /*  Hermitian matrix A stored in packed format using the factorization */
00064 /*  A = U*D*U**H or A = L*D*L**H computed by CHPTRF. */
00065 
00066 /*  Arguments */
00067 /*  ========= */
00068 
00069 /*  UPLO    (input) CHARACTER*1 */
00070 /*          Specifies whether the details of the factorization are stored */
00071 /*          as an upper or lower triangular matrix. */
00072 /*          = 'U':  Upper triangular, form is A = U*D*U**H; */
00073 /*          = 'L':  Lower triangular, form is A = L*D*L**H. */
00074 
00075 /*  N       (input) INTEGER */
00076 /*          The order of the matrix A.  N >= 0. */
00077 
00078 /*  NRHS    (input) INTEGER */
00079 /*          The number of right hand sides, i.e., the number of columns */
00080 /*          of the matrix B.  NRHS >= 0. */
00081 
00082 /*  AP      (input) COMPLEX array, dimension (N*(N+1)/2) */
00083 /*          The block diagonal matrix D and the multipliers used to */
00084 /*          obtain the factor U or L as computed by CHPTRF, stored as a */
00085 /*          packed triangular matrix. */
00086 
00087 /*  IPIV    (input) INTEGER array, dimension (N) */
00088 /*          Details of the interchanges and the block structure of D */
00089 /*          as determined by CHPTRF. */
00090 
00091 /*  B       (input/output) COMPLEX array, dimension (LDB,NRHS) */
00092 /*          On entry, the right hand side matrix B. */
00093 /*          On exit, the solution matrix X. */
00094 
00095 /*  LDB     (input) INTEGER */
00096 /*          The leading dimension of the array B.  LDB >= max(1,N). */
00097 
00098 /*  INFO    (output) INTEGER */
00099 /*          = 0:  successful exit */
00100 /*          < 0: if INFO = -i, the i-th argument had an illegal value */
00101 
00102 /*  ===================================================================== */
00103 
00104 /*     .. Parameters .. */
00105 /*     .. */
00106 /*     .. Local Scalars .. */
00107 /*     .. */
00108 /*     .. External Functions .. */
00109 /*     .. */
00110 /*     .. External Subroutines .. */
00111 /*     .. */
00112 /*     .. Intrinsic Functions .. */
00113 /*     .. */
00114 /*     .. Executable Statements .. */
00115 
00116     /* Parameter adjustments */
00117     --ap;
00118     --ipiv;
00119     b_dim1 = *ldb;
00120     b_offset = 1 + b_dim1;
00121     b -= b_offset;
00122 
00123     /* Function Body */
00124     *info = 0;
00125     upper = lsame_(uplo, "U");
00126     if (! upper && ! lsame_(uplo, "L")) {
00127         *info = -1;
00128     } else if (*n < 0) {
00129         *info = -2;
00130     } else if (*nrhs < 0) {
00131         *info = -3;
00132     } else if (*ldb < max(1,*n)) {
00133         *info = -7;
00134     }
00135     if (*info != 0) {
00136         i__1 = -(*info);
00137         xerbla_("CHPTRS", &i__1);
00138         return 0;
00139     }
00140 
00141 /*     Quick return if possible */
00142 
00143     if (*n == 0 || *nrhs == 0) {
00144         return 0;
00145     }
00146 
00147     if (upper) {
00148 
00149 /*        Solve A*X = B, where A = U*D*U'. */
00150 
00151 /*        First solve U*D*X = B, overwriting B with X. */
00152 
00153 /*        K is the main loop index, decreasing from N to 1 in steps of */
00154 /*        1 or 2, depending on the size of the diagonal blocks. */
00155 
00156         k = *n;
00157         kc = *n * (*n + 1) / 2 + 1;
00158 L10:
00159 
00160 /*        If K < 1, exit from loop. */
00161 
00162         if (k < 1) {
00163             goto L30;
00164         }
00165 
00166         kc -= k;
00167         if (ipiv[k] > 0) {
00168 
00169 /*           1 x 1 diagonal block */
00170 
00171 /*           Interchange rows K and IPIV(K). */
00172 
00173             kp = ipiv[k];
00174             if (kp != k) {
00175                 cswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1], ldb);
00176             }
00177 
00178 /*           Multiply by inv(U(K)), where U(K) is the transformation */
00179 /*           stored in column K of A. */
00180 
00181             i__1 = k - 1;
00182             q__1.r = -1.f, q__1.i = -0.f;
00183             cgeru_(&i__1, nrhs, &q__1, &ap[kc], &c__1, &b[k + b_dim1], ldb, &
00184                     b[b_dim1 + 1], ldb);
00185 
00186 /*           Multiply by the inverse of the diagonal block. */
00187 
00188             i__1 = kc + k - 1;
00189             s = 1.f / ap[i__1].r;
00190             csscal_(nrhs, &s, &b[k + b_dim1], ldb);
00191             --k;
00192         } else {
00193 
00194 /*           2 x 2 diagonal block */
00195 
00196 /*           Interchange rows K-1 and -IPIV(K). */
00197 
00198             kp = -ipiv[k];
00199             if (kp != k - 1) {
00200                 cswap_(nrhs, &b[k - 1 + b_dim1], ldb, &b[kp + b_dim1], ldb);
00201             }
00202 
00203 /*           Multiply by inv(U(K)), where U(K) is the transformation */
00204 /*           stored in columns K-1 and K of A. */
00205 
00206             i__1 = k - 2;
00207             q__1.r = -1.f, q__1.i = -0.f;
00208             cgeru_(&i__1, nrhs, &q__1, &ap[kc], &c__1, &b[k + b_dim1], ldb, &
00209                     b[b_dim1 + 1], ldb);
00210             i__1 = k - 2;
00211             q__1.r = -1.f, q__1.i = -0.f;
00212             cgeru_(&i__1, nrhs, &q__1, &ap[kc - (k - 1)], &c__1, &b[k - 1 + 
00213                     b_dim1], ldb, &b[b_dim1 + 1], ldb);
00214 
00215 /*           Multiply by the inverse of the diagonal block. */
00216 
00217             i__1 = kc + k - 2;
00218             akm1k.r = ap[i__1].r, akm1k.i = ap[i__1].i;
00219             c_div(&q__1, &ap[kc - 1], &akm1k);
00220             akm1.r = q__1.r, akm1.i = q__1.i;
00221             r_cnjg(&q__2, &akm1k);
00222             c_div(&q__1, &ap[kc + k - 1], &q__2);
00223             ak.r = q__1.r, ak.i = q__1.i;
00224             q__2.r = akm1.r * ak.r - akm1.i * ak.i, q__2.i = akm1.r * ak.i + 
00225                     akm1.i * ak.r;
00226             q__1.r = q__2.r - 1.f, q__1.i = q__2.i - 0.f;
00227             denom.r = q__1.r, denom.i = q__1.i;
00228             i__1 = *nrhs;
00229             for (j = 1; j <= i__1; ++j) {
00230                 c_div(&q__1, &b[k - 1 + j * b_dim1], &akm1k);
00231                 bkm1.r = q__1.r, bkm1.i = q__1.i;
00232                 r_cnjg(&q__2, &akm1k);
00233                 c_div(&q__1, &b[k + j * b_dim1], &q__2);
00234                 bk.r = q__1.r, bk.i = q__1.i;
00235                 i__2 = k - 1 + j * b_dim1;
00236                 q__3.r = ak.r * bkm1.r - ak.i * bkm1.i, q__3.i = ak.r * 
00237                         bkm1.i + ak.i * bkm1.r;
00238                 q__2.r = q__3.r - bk.r, q__2.i = q__3.i - bk.i;
00239                 c_div(&q__1, &q__2, &denom);
00240                 b[i__2].r = q__1.r, b[i__2].i = q__1.i;
00241                 i__2 = k + j * b_dim1;
00242                 q__3.r = akm1.r * bk.r - akm1.i * bk.i, q__3.i = akm1.r * 
00243                         bk.i + akm1.i * bk.r;
00244                 q__2.r = q__3.r - bkm1.r, q__2.i = q__3.i - bkm1.i;
00245                 c_div(&q__1, &q__2, &denom);
00246                 b[i__2].r = q__1.r, b[i__2].i = q__1.i;
00247 /* L20: */
00248             }
00249             kc = kc - k + 1;
00250             k += -2;
00251         }
00252 
00253         goto L10;
00254 L30:
00255 
00256 /*        Next solve U'*X = B, overwriting B with X. */
00257 
00258 /*        K is the main loop index, increasing from 1 to N in steps of */
00259 /*        1 or 2, depending on the size of the diagonal blocks. */
00260 
00261         k = 1;
00262         kc = 1;
00263 L40:
00264 
00265 /*        If K > N, exit from loop. */
00266 
00267         if (k > *n) {
00268             goto L50;
00269         }
00270 
00271         if (ipiv[k] > 0) {
00272 
00273 /*           1 x 1 diagonal block */
00274 
00275 /*           Multiply by inv(U'(K)), where U(K) is the transformation */
00276 /*           stored in column K of A. */
00277 
00278             if (k > 1) {
00279                 clacgv_(nrhs, &b[k + b_dim1], ldb);
00280                 i__1 = k - 1;
00281                 q__1.r = -1.f, q__1.i = -0.f;
00282                 cgemv_("Conjugate transpose", &i__1, nrhs, &q__1, &b[b_offset]
00283 , ldb, &ap[kc], &c__1, &c_b1, &b[k + b_dim1], ldb);
00284                 clacgv_(nrhs, &b[k + b_dim1], ldb);
00285             }
00286 
00287 /*           Interchange rows K and IPIV(K). */
00288 
00289             kp = ipiv[k];
00290             if (kp != k) {
00291                 cswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1], ldb);
00292             }
00293             kc += k;
00294             ++k;
00295         } else {
00296 
00297 /*           2 x 2 diagonal block */
00298 
00299 /*           Multiply by inv(U'(K+1)), where U(K+1) is the transformation */
00300 /*           stored in columns K and K+1 of A. */
00301 
00302             if (k > 1) {
00303                 clacgv_(nrhs, &b[k + b_dim1], ldb);
00304                 i__1 = k - 1;
00305                 q__1.r = -1.f, q__1.i = -0.f;
00306                 cgemv_("Conjugate transpose", &i__1, nrhs, &q__1, &b[b_offset]
00307 , ldb, &ap[kc], &c__1, &c_b1, &b[k + b_dim1], ldb);
00308                 clacgv_(nrhs, &b[k + b_dim1], ldb);
00309 
00310                 clacgv_(nrhs, &b[k + 1 + b_dim1], ldb);
00311                 i__1 = k - 1;
00312                 q__1.r = -1.f, q__1.i = -0.f;
00313                 cgemv_("Conjugate transpose", &i__1, nrhs, &q__1, &b[b_offset]
00314 , ldb, &ap[kc + k], &c__1, &c_b1, &b[k + 1 + b_dim1], 
00315                         ldb);
00316                 clacgv_(nrhs, &b[k + 1 + b_dim1], ldb);
00317             }
00318 
00319 /*           Interchange rows K and -IPIV(K). */
00320 
00321             kp = -ipiv[k];
00322             if (kp != k) {
00323                 cswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1], ldb);
00324             }
00325             kc = kc + (k << 1) + 1;
00326             k += 2;
00327         }
00328 
00329         goto L40;
00330 L50:
00331 
00332         ;
00333     } else {
00334 
00335 /*        Solve A*X = B, where A = L*D*L'. */
00336 
00337 /*        First solve L*D*X = B, overwriting B with X. */
00338 
00339 /*        K is the main loop index, increasing from 1 to N in steps of */
00340 /*        1 or 2, depending on the size of the diagonal blocks. */
00341 
00342         k = 1;
00343         kc = 1;
00344 L60:
00345 
00346 /*        If K > N, exit from loop. */
00347 
00348         if (k > *n) {
00349             goto L80;
00350         }
00351 
00352         if (ipiv[k] > 0) {
00353 
00354 /*           1 x 1 diagonal block */
00355 
00356 /*           Interchange rows K and IPIV(K). */
00357 
00358             kp = ipiv[k];
00359             if (kp != k) {
00360                 cswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1], ldb);
00361             }
00362 
00363 /*           Multiply by inv(L(K)), where L(K) is the transformation */
00364 /*           stored in column K of A. */
00365 
00366             if (k < *n) {
00367                 i__1 = *n - k;
00368                 q__1.r = -1.f, q__1.i = -0.f;
00369                 cgeru_(&i__1, nrhs, &q__1, &ap[kc + 1], &c__1, &b[k + b_dim1], 
00370                          ldb, &b[k + 1 + b_dim1], ldb);
00371             }
00372 
00373 /*           Multiply by the inverse of the diagonal block. */
00374 
00375             i__1 = kc;
00376             s = 1.f / ap[i__1].r;
00377             csscal_(nrhs, &s, &b[k + b_dim1], ldb);
00378             kc = kc + *n - k + 1;
00379             ++k;
00380         } else {
00381 
00382 /*           2 x 2 diagonal block */
00383 
00384 /*           Interchange rows K+1 and -IPIV(K). */
00385 
00386             kp = -ipiv[k];
00387             if (kp != k + 1) {
00388                 cswap_(nrhs, &b[k + 1 + b_dim1], ldb, &b[kp + b_dim1], ldb);
00389             }
00390 
00391 /*           Multiply by inv(L(K)), where L(K) is the transformation */
00392 /*           stored in columns K and K+1 of A. */
00393 
00394             if (k < *n - 1) {
00395                 i__1 = *n - k - 1;
00396                 q__1.r = -1.f, q__1.i = -0.f;
00397                 cgeru_(&i__1, nrhs, &q__1, &ap[kc + 2], &c__1, &b[k + b_dim1], 
00398                          ldb, &b[k + 2 + b_dim1], ldb);
00399                 i__1 = *n - k - 1;
00400                 q__1.r = -1.f, q__1.i = -0.f;
00401                 cgeru_(&i__1, nrhs, &q__1, &ap[kc + *n - k + 2], &c__1, &b[k 
00402                         + 1 + b_dim1], ldb, &b[k + 2 + b_dim1], ldb);
00403             }
00404 
00405 /*           Multiply by the inverse of the diagonal block. */
00406 
00407             i__1 = kc + 1;
00408             akm1k.r = ap[i__1].r, akm1k.i = ap[i__1].i;
00409             r_cnjg(&q__2, &akm1k);
00410             c_div(&q__1, &ap[kc], &q__2);
00411             akm1.r = q__1.r, akm1.i = q__1.i;
00412             c_div(&q__1, &ap[kc + *n - k + 1], &akm1k);
00413             ak.r = q__1.r, ak.i = q__1.i;
00414             q__2.r = akm1.r * ak.r - akm1.i * ak.i, q__2.i = akm1.r * ak.i + 
00415                     akm1.i * ak.r;
00416             q__1.r = q__2.r - 1.f, q__1.i = q__2.i - 0.f;
00417             denom.r = q__1.r, denom.i = q__1.i;
00418             i__1 = *nrhs;
00419             for (j = 1; j <= i__1; ++j) {
00420                 r_cnjg(&q__2, &akm1k);
00421                 c_div(&q__1, &b[k + j * b_dim1], &q__2);
00422                 bkm1.r = q__1.r, bkm1.i = q__1.i;
00423                 c_div(&q__1, &b[k + 1 + j * b_dim1], &akm1k);
00424                 bk.r = q__1.r, bk.i = q__1.i;
00425                 i__2 = k + j * b_dim1;
00426                 q__3.r = ak.r * bkm1.r - ak.i * bkm1.i, q__3.i = ak.r * 
00427                         bkm1.i + ak.i * bkm1.r;
00428                 q__2.r = q__3.r - bk.r, q__2.i = q__3.i - bk.i;
00429                 c_div(&q__1, &q__2, &denom);
00430                 b[i__2].r = q__1.r, b[i__2].i = q__1.i;
00431                 i__2 = k + 1 + j * b_dim1;
00432                 q__3.r = akm1.r * bk.r - akm1.i * bk.i, q__3.i = akm1.r * 
00433                         bk.i + akm1.i * bk.r;
00434                 q__2.r = q__3.r - bkm1.r, q__2.i = q__3.i - bkm1.i;
00435                 c_div(&q__1, &q__2, &denom);
00436                 b[i__2].r = q__1.r, b[i__2].i = q__1.i;
00437 /* L70: */
00438             }
00439             kc = kc + (*n - k << 1) + 1;
00440             k += 2;
00441         }
00442 
00443         goto L60;
00444 L80:
00445 
00446 /*        Next solve L'*X = B, overwriting B with X. */
00447 
00448 /*        K is the main loop index, decreasing from N to 1 in steps of */
00449 /*        1 or 2, depending on the size of the diagonal blocks. */
00450 
00451         k = *n;
00452         kc = *n * (*n + 1) / 2 + 1;
00453 L90:
00454 
00455 /*        If K < 1, exit from loop. */
00456 
00457         if (k < 1) {
00458             goto L100;
00459         }
00460 
00461         kc -= *n - k + 1;
00462         if (ipiv[k] > 0) {
00463 
00464 /*           1 x 1 diagonal block */
00465 
00466 /*           Multiply by inv(L'(K)), where L(K) is the transformation */
00467 /*           stored in column K of A. */
00468 
00469             if (k < *n) {
00470                 clacgv_(nrhs, &b[k + b_dim1], ldb);
00471                 i__1 = *n - k;
00472                 q__1.r = -1.f, q__1.i = -0.f;
00473                 cgemv_("Conjugate transpose", &i__1, nrhs, &q__1, &b[k + 1 + 
00474                         b_dim1], ldb, &ap[kc + 1], &c__1, &c_b1, &b[k + 
00475                         b_dim1], ldb);
00476                 clacgv_(nrhs, &b[k + b_dim1], ldb);
00477             }
00478 
00479 /*           Interchange rows K and IPIV(K). */
00480 
00481             kp = ipiv[k];
00482             if (kp != k) {
00483                 cswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1], ldb);
00484             }
00485             --k;
00486         } else {
00487 
00488 /*           2 x 2 diagonal block */
00489 
00490 /*           Multiply by inv(L'(K-1)), where L(K-1) is the transformation */
00491 /*           stored in columns K-1 and K of A. */
00492 
00493             if (k < *n) {
00494                 clacgv_(nrhs, &b[k + b_dim1], ldb);
00495                 i__1 = *n - k;
00496                 q__1.r = -1.f, q__1.i = -0.f;
00497                 cgemv_("Conjugate transpose", &i__1, nrhs, &q__1, &b[k + 1 + 
00498                         b_dim1], ldb, &ap[kc + 1], &c__1, &c_b1, &b[k + 
00499                         b_dim1], ldb);
00500                 clacgv_(nrhs, &b[k + b_dim1], ldb);
00501 
00502                 clacgv_(nrhs, &b[k - 1 + b_dim1], ldb);
00503                 i__1 = *n - k;
00504                 q__1.r = -1.f, q__1.i = -0.f;
00505                 cgemv_("Conjugate transpose", &i__1, nrhs, &q__1, &b[k + 1 + 
00506                         b_dim1], ldb, &ap[kc - (*n - k)], &c__1, &c_b1, &b[k 
00507                         - 1 + b_dim1], ldb);
00508                 clacgv_(nrhs, &b[k - 1 + b_dim1], ldb);
00509             }
00510 
00511 /*           Interchange rows K and -IPIV(K). */
00512 
00513             kp = -ipiv[k];
00514             if (kp != k) {
00515                 cswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1], ldb);
00516             }
00517             kc -= *n - k + 2;
00518             k += -2;
00519         }
00520 
00521         goto L90;
00522 L100:
00523         ;
00524     }
00525 
00526     return 0;
00527 
00528 /*     End of CHPTRS */
00529 
00530 } /* chptrs_ */