swiftnav: csytrs.c Source File

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