swiftnav: slavsy.c Source File

Go to the documentation of this file.
00001 /* slavsy.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 real c_b15 = 1.f;
00019 static integer c__1 = 1;
00020 
00021 /* Subroutine */ int slavsy_(char *uplo, char *trans, char *diag, integer *n, 
00022         integer *nrhs, real *a, integer *lda, integer *ipiv, real *b, integer 
00023         *ldb, integer *info)
00024 {
00025     /* System generated locals */
00026     integer a_dim1, a_offset, b_dim1, b_offset, i__1;
00027 
00028     /* Local variables */
00029     integer j, k;
00030     real t1, t2, d11, d12, d21, d22;
00031     integer kp;
00032     extern /* Subroutine */ int sger_(integer *, integer *, real *, real *, 
00033             integer *, real *, integer *, real *, integer *);
00034     extern logical lsame_(char *, char *);
00035     extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *), 
00036             sgemv_(char *, integer *, integer *, real *, real *, integer *, 
00037             real *, integer *, real *, real *, integer *), sswap_(
00038             integer *, real *, integer *, real *, integer *), xerbla_(char *, 
00039             integer *);
00040     logical nounit;
00041 
00042 
00043 /*  -- LAPACK auxiliary routine (version 3.1) -- */
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 /*  SLAVSY  performs one of the matrix-vector operations */
00056 /*     x := A*x  or  x := A'*x, */
00057 /*  where x is an N element vector and A is one of the factors */
00058 /*  from the block U*D*U' or L*D*L' factorization computed by SSYTRF. */
00059 
00060 /*  If TRANS = 'N', multiplies by U  or U * D  (or L  or L * D) */
00061 /*  If TRANS = 'T', multiplies by U' or D * U' (or L' or D * L') */
00062 /*  If TRANS = 'C', multiplies by U' or D * U' (or L' or D * L') */
00063 
00064 /*  Arguments */
00065 /*  ========= */
00066 
00067 /*  UPLO    (input) CHARACTER*1 */
00068 /*          Specifies whether the factor stored in A is upper or lower */
00069 /*          triangular. */
00070 /*          = 'U':  Upper triangular */
00071 /*          = 'L':  Lower triangular */
00072 
00073 /*  TRANS   (input) CHARACTER*1 */
00074 /*          Specifies the operation to be performed: */
00075 /*          = 'N':  x := A*x */
00076 /*          = 'T':  x := A'*x */
00077 /*          = 'C':  x := A'*x */
00078 
00079 /*  DIAG    (input) CHARACTER*1 */
00080 /*          Specifies whether or not the diagonal blocks are unit */
00081 /*          matrices.  If the diagonal blocks are assumed to be unit, */
00082 /*          then A = U or A = L, otherwise A = U*D or A = L*D. */
00083 /*          = 'U':  Diagonal blocks are assumed to be unit matrices. */
00084 /*          = 'N':  Diagonal blocks are assumed to be non-unit matrices. */
00085 
00086 /*  N       (input) INTEGER */
00087 /*          The number of rows and columns of the matrix A.  N >= 0. */
00088 
00089 /*  NRHS    (input) INTEGER */
00090 /*          The number of right hand sides, i.e., the number of vectors */
00091 /*          x to be multiplied by A.  NRHS >= 0. */
00092 
00093 /*  A       (input) REAL array, dimension (LDA,N) */
00094 /*          The block diagonal matrix D and the multipliers used to */
00095 /*          obtain the factor U or L as computed by SSYTRF. */
00096 
00097 /*  LDA     (input) INTEGER */
00098 /*          The leading dimension of the array A.  LDA >= max(1,N). */
00099 
00100 /*  IPIV    (input) INTEGER array, dimension (N) */
00101 /*          The pivot indices from SSYTRF. */
00102 
00103 /*  B       (input/output) REAL array, dimension (LDB,NRHS) */
00104 /*          On entry, B contains NRHS vectors of length N. */
00105 /*          On exit, B is overwritten with the product A * B. */
00106 
00107 /*  LDB     (input) INTEGER */
00108 /*          The leading dimension of the array B.  LDB >= max(1,N). */
00109 
00110 /*  INFO    (output) INTEGER */
00111 /*          = 0: successful exit */
00112 /*          < 0: if INFO = -k, the k-th argument had an illegal value */
00113 
00114 /*  ===================================================================== */
00115 
00116 /*     .. Parameters .. */
00117 /*     .. */
00118 /*     .. Local Scalars .. */
00119 /*     .. */
00120 /*     .. External Functions .. */
00121 /*     .. */
00122 /*     .. External Subroutines .. */
00123 /*     .. */
00124 /*     .. Intrinsic Functions .. */
00125 /*     .. */
00126 /*     .. Executable Statements .. */
00127 
00128 /*     Test the input parameters. */
00129 
00130     /* Parameter adjustments */
00131     a_dim1 = *lda;
00132     a_offset = 1 + a_dim1;
00133     a -= a_offset;
00134     --ipiv;
00135     b_dim1 = *ldb;
00136     b_offset = 1 + b_dim1;
00137     b -= b_offset;
00138 
00139     /* Function Body */
00140     *info = 0;
00141     if (! lsame_(uplo, "U") && ! lsame_(uplo, "L")) {
00142         *info = -1;
00143     } else if (! lsame_(trans, "N") && ! lsame_(trans, 
00144             "T") && ! lsame_(trans, "C")) {
00145         *info = -2;
00146     } else if (! lsame_(diag, "U") && ! lsame_(diag, 
00147             "N")) {
00148         *info = -3;
00149     } else if (*n < 0) {
00150         *info = -4;
00151     } else if (*lda < max(1,*n)) {
00152         *info = -6;
00153     } else if (*ldb < max(1,*n)) {
00154         *info = -9;
00155     }
00156     if (*info != 0) {
00157         i__1 = -(*info);
00158         xerbla_("SLAVSY ", &i__1);
00159         return 0;
00160     }
00161 
00162 /*     Quick return if possible. */
00163 
00164     if (*n == 0) {
00165         return 0;
00166     }
00167 
00168     nounit = lsame_(diag, "N");
00169 /* ------------------------------------------ */
00170 
00171 /*     Compute  B := A * B  (No transpose) */
00172 
00173 /* ------------------------------------------ */
00174     if (lsame_(trans, "N")) {
00175 
00176 /*        Compute  B := U*B */
00177 /*        where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1)) */
00178 
00179         if (lsame_(uplo, "U")) {
00180 
00181 /*        Loop forward applying the transformations. */
00182 
00183             k = 1;
00184 L10:
00185             if (k > *n) {
00186                 goto L30;
00187             }
00188             if (ipiv[k] > 0) {
00189 
00190 /*              1 x 1 pivot block */
00191 
00192 /*              Multiply by the diagonal element if forming U * D. */
00193 
00194                 if (nounit) {
00195                     sscal_(nrhs, &a[k + k * a_dim1], &b[k + b_dim1], ldb);
00196                 }
00197 
00198 /*              Multiply by  P(K) * inv(U(K))  if K > 1. */
00199 
00200                 if (k > 1) {
00201 
00202 /*                 Apply the transformation. */
00203 
00204                     i__1 = k - 1;
00205                     sger_(&i__1, nrhs, &c_b15, &a[k * a_dim1 + 1], &c__1, &b[
00206                             k + b_dim1], ldb, &b[b_dim1 + 1], ldb);
00207 
00208 /*                 Interchange if P(K) .ne. I. */
00209 
00210                     kp = ipiv[k];
00211                     if (kp != k) {
00212                         sswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1], 
00213                                 ldb);
00214                     }
00215                 }
00216                 ++k;
00217             } else {
00218 
00219 /*              2 x 2 pivot block */
00220 
00221 /*              Multiply by the diagonal block if forming U * D. */
00222 
00223                 if (nounit) {
00224                     d11 = a[k + k * a_dim1];
00225                     d22 = a[k + 1 + (k + 1) * a_dim1];
00226                     d12 = a[k + (k + 1) * a_dim1];
00227                     d21 = d12;
00228                     i__1 = *nrhs;
00229                     for (j = 1; j <= i__1; ++j) {
00230                         t1 = b[k + j * b_dim1];
00231                         t2 = b[k + 1 + j * b_dim1];
00232                         b[k + j * b_dim1] = d11 * t1 + d12 * t2;
00233                         b[k + 1 + j * b_dim1] = d21 * t1 + d22 * t2;
00234 /* L20: */
00235                     }
00236                 }
00237 
00238 /*              Multiply by  P(K) * inv(U(K))  if K > 1. */
00239 
00240                 if (k > 1) {
00241 
00242 /*                 Apply the transformations. */
00243 
00244                     i__1 = k - 1;
00245                     sger_(&i__1, nrhs, &c_b15, &a[k * a_dim1 + 1], &c__1, &b[
00246                             k + b_dim1], ldb, &b[b_dim1 + 1], ldb);
00247                     i__1 = k - 1;
00248                     sger_(&i__1, nrhs, &c_b15, &a[(k + 1) * a_dim1 + 1], &
00249                             c__1, &b[k + 1 + b_dim1], ldb, &b[b_dim1 + 1], 
00250                             ldb);
00251 
00252 /*                 Interchange if P(K) .ne. I. */
00253 
00254                     kp = (i__1 = ipiv[k], abs(i__1));
00255                     if (kp != k) {
00256                         sswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1], 
00257                                 ldb);
00258                     }
00259                 }
00260                 k += 2;
00261             }
00262             goto L10;
00263 L30:
00264 
00265 /*        Compute  B := L*B */
00266 /*        where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m)) . */
00267 
00268             ;
00269         } else {
00270 
00271 /*           Loop backward applying the transformations to B. */
00272 
00273             k = *n;
00274 L40:
00275             if (k < 1) {
00276                 goto L60;
00277             }
00278 
00279 /*           Test the pivot index.  If greater than zero, a 1 x 1 */
00280 /*           pivot was used, otherwise a 2 x 2 pivot was used. */
00281 
00282             if (ipiv[k] > 0) {
00283 
00284 /*              1 x 1 pivot block: */
00285 
00286 /*              Multiply by the diagonal element if forming L * D. */
00287 
00288                 if (nounit) {
00289                     sscal_(nrhs, &a[k + k * a_dim1], &b[k + b_dim1], ldb);
00290                 }
00291 
00292 /*              Multiply by  P(K) * inv(L(K))  if K < N. */
00293 
00294                 if (k != *n) {
00295                     kp = ipiv[k];
00296 
00297 /*                 Apply the transformation. */
00298 
00299                     i__1 = *n - k;
00300                     sger_(&i__1, nrhs, &c_b15, &a[k + 1 + k * a_dim1], &c__1, 
00301                             &b[k + b_dim1], ldb, &b[k + 1 + b_dim1], ldb);
00302 
00303 /*                 Interchange if a permutation was applied at the */
00304 /*                 K-th step of the factorization. */
00305 
00306                     if (kp != k) {
00307                         sswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1], 
00308                                 ldb);
00309                     }
00310                 }
00311                 --k;
00312 
00313             } else {
00314 
00315 /*              2 x 2 pivot block: */
00316 
00317 /*              Multiply by the diagonal block if forming L * D. */
00318 
00319                 if (nounit) {
00320                     d11 = a[k - 1 + (k - 1) * a_dim1];
00321                     d22 = a[k + k * a_dim1];
00322                     d21 = a[k + (k - 1) * a_dim1];
00323                     d12 = d21;
00324                     i__1 = *nrhs;
00325                     for (j = 1; j <= i__1; ++j) {
00326                         t1 = b[k - 1 + j * b_dim1];
00327                         t2 = b[k + j * b_dim1];
00328                         b[k - 1 + j * b_dim1] = d11 * t1 + d12 * t2;
00329                         b[k + j * b_dim1] = d21 * t1 + d22 * t2;
00330 /* L50: */
00331                     }
00332                 }
00333 
00334 /*              Multiply by  P(K) * inv(L(K))  if K < N. */
00335 
00336                 if (k != *n) {
00337 
00338 /*                 Apply the transformation. */
00339 
00340                     i__1 = *n - k;
00341                     sger_(&i__1, nrhs, &c_b15, &a[k + 1 + k * a_dim1], &c__1, 
00342                             &b[k + b_dim1], ldb, &b[k + 1 + b_dim1], ldb);
00343                     i__1 = *n - k;
00344                     sger_(&i__1, nrhs, &c_b15, &a[k + 1 + (k - 1) * a_dim1], &
00345                             c__1, &b[k - 1 + b_dim1], ldb, &b[k + 1 + b_dim1], 
00346                              ldb);
00347 
00348 /*                 Interchange if a permutation was applied at the */
00349 /*                 K-th step of the factorization. */
00350 
00351                     kp = (i__1 = ipiv[k], abs(i__1));
00352                     if (kp != k) {
00353                         sswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1], 
00354                                 ldb);
00355                     }
00356                 }
00357                 k += -2;
00358             }
00359             goto L40;
00360 L60:
00361             ;
00362         }
00363 /* ---------------------------------------- */
00364 
00365 /*     Compute  B := A' * B  (transpose) */
00366 
00367 /* ---------------------------------------- */
00368     } else {
00369 
00370 /*        Form  B := U'*B */
00371 /*        where U  = P(m)*inv(U(m))* ... *P(1)*inv(U(1)) */
00372 /*        and   U' = inv(U'(1))*P(1)* ... *inv(U'(m))*P(m) */
00373 
00374         if (lsame_(uplo, "U")) {
00375 
00376 /*           Loop backward applying the transformations. */
00377 
00378             k = *n;
00379 L70:
00380             if (k < 1) {
00381                 goto L90;
00382             }
00383 
00384 /*           1 x 1 pivot block. */
00385 
00386             if (ipiv[k] > 0) {
00387                 if (k > 1) {
00388 
00389 /*                 Interchange if P(K) .ne. I. */
00390 
00391                     kp = ipiv[k];
00392                     if (kp != k) {
00393                         sswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1], 
00394                                 ldb);
00395                     }
00396 
00397 /*                 Apply the transformation */
00398 
00399                     i__1 = k - 1;
00400                     sgemv_("Transpose", &i__1, nrhs, &c_b15, &b[b_offset], 
00401                             ldb, &a[k * a_dim1 + 1], &c__1, &c_b15, &b[k + 
00402                             b_dim1], ldb);
00403                 }
00404                 if (nounit) {
00405                     sscal_(nrhs, &a[k + k * a_dim1], &b[k + b_dim1], ldb);
00406                 }
00407                 --k;
00408 
00409 /*           2 x 2 pivot block. */
00410 
00411             } else {
00412                 if (k > 2) {
00413 
00414 /*                 Interchange if P(K) .ne. I. */
00415 
00416                     kp = (i__1 = ipiv[k], abs(i__1));
00417                     if (kp != k - 1) {
00418                         sswap_(nrhs, &b[k - 1 + b_dim1], ldb, &b[kp + b_dim1], 
00419                                  ldb);
00420                     }
00421 
00422 /*                 Apply the transformations */
00423 
00424                     i__1 = k - 2;
00425                     sgemv_("Transpose", &i__1, nrhs, &c_b15, &b[b_offset], 
00426                             ldb, &a[k * a_dim1 + 1], &c__1, &c_b15, &b[k + 
00427                             b_dim1], ldb);
00428                     i__1 = k - 2;
00429                     sgemv_("Transpose", &i__1, nrhs, &c_b15, &b[b_offset], 
00430                             ldb, &a[(k - 1) * a_dim1 + 1], &c__1, &c_b15, &b[
00431                             k - 1 + b_dim1], ldb);
00432                 }
00433 
00434 /*              Multiply by the diagonal block if non-unit. */
00435 
00436                 if (nounit) {
00437                     d11 = a[k - 1 + (k - 1) * a_dim1];
00438                     d22 = a[k + k * a_dim1];
00439                     d12 = a[k - 1 + k * a_dim1];
00440                     d21 = d12;
00441                     i__1 = *nrhs;
00442                     for (j = 1; j <= i__1; ++j) {
00443                         t1 = b[k - 1 + j * b_dim1];
00444                         t2 = b[k + j * b_dim1];
00445                         b[k - 1 + j * b_dim1] = d11 * t1 + d12 * t2;
00446                         b[k + j * b_dim1] = d21 * t1 + d22 * t2;
00447 /* L80: */
00448                     }
00449                 }
00450                 k += -2;
00451             }
00452             goto L70;
00453 L90:
00454 
00455 /*        Form  B := L'*B */
00456 /*        where L  = P(1)*inv(L(1))* ... *P(m)*inv(L(m)) */
00457 /*        and   L' = inv(L'(m))*P(m)* ... *inv(L'(1))*P(1) */
00458 
00459             ;
00460         } else {
00461 
00462 /*           Loop forward applying the L-transformations. */
00463 
00464             k = 1;
00465 L100:
00466             if (k > *n) {
00467                 goto L120;
00468             }
00469 
00470 /*           1 x 1 pivot block */
00471 
00472             if (ipiv[k] > 0) {
00473                 if (k < *n) {
00474 
00475 /*                 Interchange if P(K) .ne. I. */
00476 
00477                     kp = ipiv[k];
00478                     if (kp != k) {
00479                         sswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1], 
00480                                 ldb);
00481                     }
00482 
00483 /*                 Apply the transformation */
00484 
00485                     i__1 = *n - k;
00486                     sgemv_("Transpose", &i__1, nrhs, &c_b15, &b[k + 1 + 
00487                             b_dim1], ldb, &a[k + 1 + k * a_dim1], &c__1, &
00488                             c_b15, &b[k + b_dim1], ldb);
00489                 }
00490                 if (nounit) {
00491                     sscal_(nrhs, &a[k + k * a_dim1], &b[k + b_dim1], ldb);
00492                 }
00493                 ++k;
00494 
00495 /*           2 x 2 pivot block. */
00496 
00497             } else {
00498                 if (k < *n - 1) {
00499 
00500 /*              Interchange if P(K) .ne. I. */
00501 
00502                     kp = (i__1 = ipiv[k], abs(i__1));
00503                     if (kp != k + 1) {
00504                         sswap_(nrhs, &b[k + 1 + b_dim1], ldb, &b[kp + b_dim1], 
00505                                  ldb);
00506                     }
00507 
00508 /*                 Apply the transformation */
00509 
00510                     i__1 = *n - k - 1;
00511                     sgemv_("Transpose", &i__1, nrhs, &c_b15, &b[k + 2 + 
00512                             b_dim1], ldb, &a[k + 2 + (k + 1) * a_dim1], &c__1, 
00513                              &c_b15, &b[k + 1 + b_dim1], ldb);
00514                     i__1 = *n - k - 1;
00515                     sgemv_("Transpose", &i__1, nrhs, &c_b15, &b[k + 2 + 
00516                             b_dim1], ldb, &a[k + 2 + k * a_dim1], &c__1, &
00517                             c_b15, &b[k + b_dim1], ldb);
00518                 }
00519 
00520 /*              Multiply by the diagonal block if non-unit. */
00521 
00522                 if (nounit) {
00523                     d11 = a[k + k * a_dim1];
00524                     d22 = a[k + 1 + (k + 1) * a_dim1];
00525                     d21 = a[k + 1 + k * a_dim1];
00526                     d12 = d21;
00527                     i__1 = *nrhs;
00528                     for (j = 1; j <= i__1; ++j) {
00529                         t1 = b[k + j * b_dim1];
00530                         t2 = b[k + 1 + j * b_dim1];
00531                         b[k + j * b_dim1] = d11 * t1 + d12 * t2;
00532                         b[k + 1 + j * b_dim1] = d21 * t1 + d22 * t2;
00533 /* L110: */
00534                     }
00535                 }
00536                 k += 2;
00537             }
00538             goto L100;
00539 L120:
00540             ;
00541         }
00542 
00543     }
00544     return 0;
00545 
00546 /*     End of SLAVSY */
00547 
00548 } /* slavsy_ */