clavsp.c
Go to the documentation of this file.
00001 /* clavsp.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 clavsp_(char *uplo, char *trans, char *diag, integer *n, 
00022         integer *nrhs, complex *a, integer *ipiv, complex *b, integer *ldb, 
00023         integer *info)
00024 {
00025     /* System generated locals */
00026     integer b_dim1, b_offset, i__1, i__2;
00027     complex q__1, q__2, q__3;
00028 
00029     /* Local variables */
00030     integer j, k;
00031     complex t1, t2, d11, d12, d21, d22;
00032     integer kc, kp;
00033     extern /* Subroutine */ int cscal_(integer *, complex *, complex *, 
00034             integer *);
00035     extern logical lsame_(char *, char *);
00036     extern /* Subroutine */ int cgemv_(char *, integer *, integer *, complex *
00037 , complex *, integer *, complex *, integer *, complex *, complex *
00038 , integer *), cgeru_(integer *, integer *, complex *, 
00039             complex *, integer *, complex *, integer *, complex *, integer *),
00040              cswap_(integer *, complex *, integer *, complex *, integer *), 
00041             xerbla_(char *, integer *);
00042     integer kcnext;
00043     logical nounit;
00044 
00045 
00046 /*  -- LAPACK auxiliary routine (version 3.1) -- */
00047 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00048 /*     November 2006 */
00049 
00050 /*     .. Scalar Arguments .. */
00051 /*     .. */
00052 /*     .. Array Arguments .. */
00053 /*     .. */
00054 
00055 /*  Purpose */
00056 /*  ======= */
00057 
00058 /*     CLAVSP  performs one of the matrix-vector operations */
00059 /*        x := A*x  or  x := A^T*x, */
00060 /*     where x is an N element vector and  A is one of the factors */
00061 /*     from the symmetric factorization computed by CSPTRF. */
00062 /*     CSPTRF produces a factorization of the form */
00063 /*          U * D * U^T     or     L * D * L^T, */
00064 /*     where U (or L) is a product of permutation and unit upper (lower) */
00065 /*     triangular matrices, U^T (or L^T) is the transpose of */
00066 /*     U (or L), and D is symmetric and block diagonal with 1 x 1 and */
00067 /*     2 x 2 diagonal blocks.  The multipliers for the transformations */
00068 /*     and the upper or lower triangular parts of the diagonal blocks */
00069 /*     are stored columnwise in packed format in the linear array A. */
00070 
00071 /*     If TRANS = 'N' or 'n', CLAVSP multiplies either by U or U * D */
00072 /*     (or L or L * D). */
00073 /*     If TRANS = 'C' or 'c', CLAVSP multiplies either by U^T or D * U^T */
00074 /*     (or L^T or D * L^T ). */
00075 
00076 /*  Arguments */
00077 /*  ========== */
00078 
00079 /*  UPLO   - CHARACTER*1 */
00080 /*           On entry, UPLO specifies whether the triangular matrix */
00081 /*           stored in A is upper or lower triangular. */
00082 /*              UPLO = 'U' or 'u'   The matrix is upper triangular. */
00083 /*              UPLO = 'L' or 'l'   The matrix is lower triangular. */
00084 /*           Unchanged on exit. */
00085 
00086 /*  TRANS  - CHARACTER*1 */
00087 /*           On entry, TRANS specifies the operation to be performed as */
00088 /*           follows: */
00089 /*              TRANS = 'N' or 'n'   x := A*x. */
00090 /*              TRANS = 'T' or 't'   x := A^T*x. */
00091 /*           Unchanged on exit. */
00092 
00093 /*  DIAG   - CHARACTER*1 */
00094 /*           On entry, DIAG specifies whether the diagonal blocks are */
00095 /*           assumed to be unit matrices, as follows: */
00096 /*              DIAG = 'U' or 'u'   Diagonal blocks are unit matrices. */
00097 /*              DIAG = 'N' or 'n'   Diagonal blocks are non-unit. */
00098 /*           Unchanged on exit. */
00099 
00100 /*  N      - INTEGER */
00101 /*           On entry, N specifies the order of the matrix A. */
00102 /*           N must be at least zero. */
00103 /*           Unchanged on exit. */
00104 
00105 /*  NRHS   - INTEGER */
00106 /*           On entry, NRHS specifies the number of right hand sides, */
00107 /*           i.e., the number of vectors x to be multiplied by A. */
00108 /*           NRHS must be at least zero. */
00109 /*           Unchanged on exit. */
00110 
00111 /*  A      - COMPLEX array, dimension( N*(N+1)/2 ) */
00112 /*           On entry, A contains a block diagonal matrix and the */
00113 /*           multipliers of the transformations used to obtain it, */
00114 /*           stored as a packed triangular matrix. */
00115 /*           Unchanged on exit. */
00116 
00117 /*  IPIV   - INTEGER array, dimension( N ) */
00118 /*           On entry, IPIV contains the vector of pivot indices as */
00119 /*           determined by CSPTRF. */
00120 /*           If IPIV( K ) = K, no interchange was done. */
00121 /*           If IPIV( K ) <> K but IPIV( K ) > 0, then row K was inter- */
00122 /*           changed with row IPIV( K ) and a 1 x 1 pivot block was used. */
00123 /*           If IPIV( K ) < 0 and UPLO = 'U', then row K-1 was exchanged */
00124 /*           with row | IPIV( K ) | and a 2 x 2 pivot block was used. */
00125 /*           If IPIV( K ) < 0 and UPLO = 'L', then row K+1 was exchanged */
00126 /*           with row | IPIV( K ) | and a 2 x 2 pivot block was used. */
00127 
00128 /*  B      - COMPLEX array, dimension( LDB, NRHS ) */
00129 /*           On entry, B contains NRHS vectors of length N. */
00130 /*           On exit, B is overwritten with the product A * B. */
00131 
00132 /*  LDB    - INTEGER */
00133 /*           On entry, LDB contains the leading dimension of B as */
00134 /*           declared in the calling program.  LDB must be at least */
00135 /*           max( 1, N ). */
00136 /*           Unchanged on exit. */
00137 
00138 /*  INFO   - INTEGER */
00139 /*           INFO is the error flag. */
00140 /*           On exit, a value of 0 indicates a successful exit. */
00141 /*           A negative value, say -K, indicates that the K-th argument */
00142 /*           has an illegal value. */
00143 
00144 /*  ===================================================================== */
00145 
00146 /*     .. Parameters .. */
00147 /*     .. */
00148 /*     .. Local Scalars .. */
00149 /*     .. */
00150 /*     .. External Functions .. */
00151 /*     .. */
00152 /*     .. External Subroutines .. */
00153 /*     .. */
00154 /*     .. Intrinsic Functions .. */
00155 /*     .. */
00156 /*     .. Executable Statements .. */
00157 
00158 /*     Test the input parameters. */
00159 
00160     /* Parameter adjustments */
00161     --a;
00162     --ipiv;
00163     b_dim1 = *ldb;
00164     b_offset = 1 + b_dim1;
00165     b -= b_offset;
00166 
00167     /* Function Body */
00168     *info = 0;
00169     if (! lsame_(uplo, "U") && ! lsame_(uplo, "L")) {
00170         *info = -1;
00171     } else if (! lsame_(trans, "N") && ! lsame_(trans, 
00172             "T")) {
00173         *info = -2;
00174     } else if (! lsame_(diag, "U") && ! lsame_(diag, 
00175             "N")) {
00176         *info = -3;
00177     } else if (*n < 0) {
00178         *info = -4;
00179     } else if (*ldb < max(1,*n)) {
00180         *info = -8;
00181     }
00182     if (*info != 0) {
00183         i__1 = -(*info);
00184         xerbla_("CLAVSP ", &i__1);
00185         return 0;
00186     }
00187 
00188 /*     Quick return if possible. */
00189 
00190     if (*n == 0) {
00191         return 0;
00192     }
00193 
00194     nounit = lsame_(diag, "N");
00195 /* ------------------------------------------ */
00196 
00197 /*     Compute  B := A * B  (No transpose) */
00198 
00199 /* ------------------------------------------ */
00200     if (lsame_(trans, "N")) {
00201 
00202 /*        Compute  B := U*B */
00203 /*        where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1)) */
00204 
00205         if (lsame_(uplo, "U")) {
00206 
00207 /*        Loop forward applying the transformations. */
00208 
00209             k = 1;
00210             kc = 1;
00211 L10:
00212             if (k > *n) {
00213                 goto L30;
00214             }
00215 
00216 /*           1 x 1 pivot block */
00217 
00218             if (ipiv[k] > 0) {
00219 
00220 /*              Multiply by the diagonal element if forming U * D. */
00221 
00222                 if (nounit) {
00223                     cscal_(nrhs, &a[kc + k - 1], &b[k + b_dim1], ldb);
00224                 }
00225 
00226 /*              Multiply by P(K) * inv(U(K))  if K > 1. */
00227 
00228                 if (k > 1) {
00229 
00230 /*                 Apply the transformation. */
00231 
00232                     i__1 = k - 1;
00233                     cgeru_(&i__1, nrhs, &c_b1, &a[kc], &c__1, &b[k + b_dim1], 
00234                             ldb, &b[b_dim1 + 1], ldb);
00235 
00236 /*                 Interchange if P(K) != I. */
00237 
00238                     kp = ipiv[k];
00239                     if (kp != k) {
00240                         cswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1], 
00241                                 ldb);
00242                     }
00243                 }
00244                 kc += k;
00245                 ++k;
00246             } else {
00247 
00248 /*              2 x 2 pivot block */
00249 
00250                 kcnext = kc + k;
00251 
00252 /*              Multiply by the diagonal block if forming U * D. */
00253 
00254                 if (nounit) {
00255                     i__1 = kcnext - 1;
00256                     d11.r = a[i__1].r, d11.i = a[i__1].i;
00257                     i__1 = kcnext + k;
00258                     d22.r = a[i__1].r, d22.i = a[i__1].i;
00259                     i__1 = kcnext + k - 1;
00260                     d12.r = a[i__1].r, d12.i = a[i__1].i;
00261                     d21.r = d12.r, d21.i = d12.i;
00262                     i__1 = *nrhs;
00263                     for (j = 1; j <= i__1; ++j) {
00264                         i__2 = k + j * b_dim1;
00265                         t1.r = b[i__2].r, t1.i = b[i__2].i;
00266                         i__2 = k + 1 + j * b_dim1;
00267                         t2.r = b[i__2].r, t2.i = b[i__2].i;
00268                         i__2 = k + j * b_dim1;
00269                         q__2.r = d11.r * t1.r - d11.i * t1.i, q__2.i = d11.r *
00270                                  t1.i + d11.i * t1.r;
00271                         q__3.r = d12.r * t2.r - d12.i * t2.i, q__3.i = d12.r *
00272                                  t2.i + d12.i * t2.r;
00273                         q__1.r = q__2.r + q__3.r, q__1.i = q__2.i + q__3.i;
00274                         b[i__2].r = q__1.r, b[i__2].i = q__1.i;
00275                         i__2 = k + 1 + j * b_dim1;
00276                         q__2.r = d21.r * t1.r - d21.i * t1.i, q__2.i = d21.r *
00277                                  t1.i + d21.i * t1.r;
00278                         q__3.r = d22.r * t2.r - d22.i * t2.i, q__3.i = d22.r *
00279                                  t2.i + d22.i * t2.r;
00280                         q__1.r = q__2.r + q__3.r, q__1.i = q__2.i + q__3.i;
00281                         b[i__2].r = q__1.r, b[i__2].i = q__1.i;
00282 /* L20: */
00283                     }
00284                 }
00285 
00286 /*              Multiply by  P(K) * inv(U(K))  if K > 1. */
00287 
00288                 if (k > 1) {
00289 
00290 /*                 Apply the transformations. */
00291 
00292                     i__1 = k - 1;
00293                     cgeru_(&i__1, nrhs, &c_b1, &a[kc], &c__1, &b[k + b_dim1], 
00294                             ldb, &b[b_dim1 + 1], ldb);
00295                     i__1 = k - 1;
00296                     cgeru_(&i__1, nrhs, &c_b1, &a[kcnext], &c__1, &b[k + 1 + 
00297                             b_dim1], ldb, &b[b_dim1 + 1], ldb);
00298 
00299 /*                 Interchange if P(K) != I. */
00300 
00301                     kp = (i__1 = ipiv[k], abs(i__1));
00302                     if (kp != k) {
00303                         cswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1], 
00304                                 ldb);
00305                     }
00306                 }
00307                 kc = kcnext + k + 1;
00308                 k += 2;
00309             }
00310             goto L10;
00311 L30:
00312 
00313 /*        Compute  B := L*B */
00314 /*        where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m)) . */
00315 
00316             ;
00317         } else {
00318 
00319 /*           Loop backward applying the transformations to B. */
00320 
00321             k = *n;
00322             kc = *n * (*n + 1) / 2 + 1;
00323 L40:
00324             if (k < 1) {
00325                 goto L60;
00326             }
00327             kc -= *n - k + 1;
00328 
00329 /*           Test the pivot index.  If greater than zero, a 1 x 1 */
00330 /*           pivot was used, otherwise a 2 x 2 pivot was used. */
00331 
00332             if (ipiv[k] > 0) {
00333 
00334 /*              1 x 1 pivot block: */
00335 
00336 /*              Multiply by the diagonal element if forming L * D. */
00337 
00338                 if (nounit) {
00339                     cscal_(nrhs, &a[kc], &b[k + b_dim1], ldb);
00340                 }
00341 
00342 /*              Multiply by  P(K) * inv(L(K))  if K < N. */
00343 
00344                 if (k != *n) {
00345                     kp = ipiv[k];
00346 
00347 /*                 Apply the transformation. */
00348 
00349                     i__1 = *n - k;
00350                     cgeru_(&i__1, nrhs, &c_b1, &a[kc + 1], &c__1, &b[k + 
00351                             b_dim1], ldb, &b[k + 1 + b_dim1], ldb);
00352 
00353 /*                 Interchange if a permutation was applied at the */
00354 /*                 K-th step of the factorization. */
00355 
00356                     if (kp != k) {
00357                         cswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1], 
00358                                 ldb);
00359                     }
00360                 }
00361                 --k;
00362 
00363             } else {
00364 
00365 /*              2 x 2 pivot block: */
00366 
00367                 kcnext = kc - (*n - k + 2);
00368 
00369 /*              Multiply by the diagonal block if forming L * D. */
00370 
00371                 if (nounit) {
00372                     i__1 = kcnext;
00373                     d11.r = a[i__1].r, d11.i = a[i__1].i;
00374                     i__1 = kc;
00375                     d22.r = a[i__1].r, d22.i = a[i__1].i;
00376                     i__1 = kcnext + 1;
00377                     d21.r = a[i__1].r, d21.i = a[i__1].i;
00378                     d12.r = d21.r, d12.i = d21.i;
00379                     i__1 = *nrhs;
00380                     for (j = 1; j <= i__1; ++j) {
00381                         i__2 = k - 1 + j * b_dim1;
00382                         t1.r = b[i__2].r, t1.i = b[i__2].i;
00383                         i__2 = k + j * b_dim1;
00384                         t2.r = b[i__2].r, t2.i = b[i__2].i;
00385                         i__2 = k - 1 + j * b_dim1;
00386                         q__2.r = d11.r * t1.r - d11.i * t1.i, q__2.i = d11.r *
00387                                  t1.i + d11.i * t1.r;
00388                         q__3.r = d12.r * t2.r - d12.i * t2.i, q__3.i = d12.r *
00389                                  t2.i + d12.i * t2.r;
00390                         q__1.r = q__2.r + q__3.r, q__1.i = q__2.i + q__3.i;
00391                         b[i__2].r = q__1.r, b[i__2].i = q__1.i;
00392                         i__2 = k + j * b_dim1;
00393                         q__2.r = d21.r * t1.r - d21.i * t1.i, q__2.i = d21.r *
00394                                  t1.i + d21.i * t1.r;
00395                         q__3.r = d22.r * t2.r - d22.i * t2.i, q__3.i = d22.r *
00396                                  t2.i + d22.i * t2.r;
00397                         q__1.r = q__2.r + q__3.r, q__1.i = q__2.i + q__3.i;
00398                         b[i__2].r = q__1.r, b[i__2].i = q__1.i;
00399 /* L50: */
00400                     }
00401                 }
00402 
00403 /*              Multiply by  P(K) * inv(L(K))  if K < N. */
00404 
00405                 if (k != *n) {
00406 
00407 /*                 Apply the transformation. */
00408 
00409                     i__1 = *n - k;
00410                     cgeru_(&i__1, nrhs, &c_b1, &a[kc + 1], &c__1, &b[k + 
00411                             b_dim1], ldb, &b[k + 1 + b_dim1], ldb);
00412                     i__1 = *n - k;
00413                     cgeru_(&i__1, nrhs, &c_b1, &a[kcnext + 2], &c__1, &b[k - 
00414                             1 + b_dim1], ldb, &b[k + 1 + b_dim1], ldb);
00415 
00416 /*                 Interchange if a permutation was applied at the */
00417 /*                 K-th step of the factorization. */
00418 
00419                     kp = (i__1 = ipiv[k], abs(i__1));
00420                     if (kp != k) {
00421                         cswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1], 
00422                                 ldb);
00423                     }
00424                 }
00425                 kc = kcnext;
00426                 k += -2;
00427             }
00428             goto L40;
00429 L60:
00430             ;
00431         }
00432 /* ------------------------------------------------- */
00433 
00434 /*     Compute  B := A^T * B  (transpose) */
00435 
00436 /* ------------------------------------------------- */
00437     } else {
00438 
00439 /*        Form  B := U^T*B */
00440 /*        where U  = P(m)*inv(U(m))* ... *P(1)*inv(U(1)) */
00441 /*        and   U^T = inv(U^T(1))*P(1)* ... *inv(U^T(m))*P(m) */
00442 
00443         if (lsame_(uplo, "U")) {
00444 
00445 /*           Loop backward applying the transformations. */
00446 
00447             k = *n;
00448             kc = *n * (*n + 1) / 2 + 1;
00449 L70:
00450             if (k < 1) {
00451                 goto L90;
00452             }
00453             kc -= k;
00454 
00455 /*           1 x 1 pivot block. */
00456 
00457             if (ipiv[k] > 0) {
00458                 if (k > 1) {
00459 
00460 /*                 Interchange if P(K) != I. */
00461 
00462                     kp = ipiv[k];
00463                     if (kp != k) {
00464                         cswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1], 
00465                                 ldb);
00466                     }
00467 
00468 /*                 Apply the transformation: */
00469 /*                    y := y - B' * conjg(x) */
00470 /*                 where x is a column of A and y is a row of B. */
00471 
00472                     i__1 = k - 1;
00473                     cgemv_("Transpose", &i__1, nrhs, &c_b1, &b[b_offset], ldb, 
00474                              &a[kc], &c__1, &c_b1, &b[k + b_dim1], ldb);
00475                 }
00476                 if (nounit) {
00477                     cscal_(nrhs, &a[kc + k - 1], &b[k + b_dim1], ldb);
00478                 }
00479                 --k;
00480 
00481 /*           2 x 2 pivot block. */
00482 
00483             } else {
00484                 kcnext = kc - (k - 1);
00485                 if (k > 2) {
00486 
00487 /*                 Interchange if P(K) != I. */
00488 
00489                     kp = (i__1 = ipiv[k], abs(i__1));
00490                     if (kp != k - 1) {
00491                         cswap_(nrhs, &b[k - 1 + b_dim1], ldb, &b[kp + b_dim1], 
00492                                  ldb);
00493                     }
00494 
00495 /*                 Apply the transformations. */
00496 
00497                     i__1 = k - 2;
00498                     cgemv_("Transpose", &i__1, nrhs, &c_b1, &b[b_offset], ldb, 
00499                              &a[kc], &c__1, &c_b1, &b[k + b_dim1], ldb);
00500 
00501                     i__1 = k - 2;
00502                     cgemv_("Transpose", &i__1, nrhs, &c_b1, &b[b_offset], ldb, 
00503                              &a[kcnext], &c__1, &c_b1, &b[k - 1 + b_dim1], 
00504                             ldb);
00505                 }
00506 
00507 /*              Multiply by the diagonal block if non-unit. */
00508 
00509                 if (nounit) {
00510                     i__1 = kc - 1;
00511                     d11.r = a[i__1].r, d11.i = a[i__1].i;
00512                     i__1 = kc + k - 1;
00513                     d22.r = a[i__1].r, d22.i = a[i__1].i;
00514                     i__1 = kc + k - 2;
00515                     d12.r = a[i__1].r, d12.i = a[i__1].i;
00516                     d21.r = d12.r, d21.i = d12.i;
00517                     i__1 = *nrhs;
00518                     for (j = 1; j <= i__1; ++j) {
00519                         i__2 = k - 1 + j * b_dim1;
00520                         t1.r = b[i__2].r, t1.i = b[i__2].i;
00521                         i__2 = k + j * b_dim1;
00522                         t2.r = b[i__2].r, t2.i = b[i__2].i;
00523                         i__2 = k - 1 + j * b_dim1;
00524                         q__2.r = d11.r * t1.r - d11.i * t1.i, q__2.i = d11.r *
00525                                  t1.i + d11.i * t1.r;
00526                         q__3.r = d12.r * t2.r - d12.i * t2.i, q__3.i = d12.r *
00527                                  t2.i + d12.i * t2.r;
00528                         q__1.r = q__2.r + q__3.r, q__1.i = q__2.i + q__3.i;
00529                         b[i__2].r = q__1.r, b[i__2].i = q__1.i;
00530                         i__2 = k + j * b_dim1;
00531                         q__2.r = d21.r * t1.r - d21.i * t1.i, q__2.i = d21.r *
00532                                  t1.i + d21.i * t1.r;
00533                         q__3.r = d22.r * t2.r - d22.i * t2.i, q__3.i = d22.r *
00534                                  t2.i + d22.i * t2.r;
00535                         q__1.r = q__2.r + q__3.r, q__1.i = q__2.i + q__3.i;
00536                         b[i__2].r = q__1.r, b[i__2].i = q__1.i;
00537 /* L80: */
00538                     }
00539                 }
00540                 kc = kcnext;
00541                 k += -2;
00542             }
00543             goto L70;
00544 L90:
00545 
00546 /*        Form  B := L^T*B */
00547 /*        where L  = P(1)*inv(L(1))* ... *P(m)*inv(L(m)) */
00548 /*        and   L^T = inv(L(m))*P(m)* ... *inv(L(1))*P(1) */
00549 
00550             ;
00551         } else {
00552 
00553 /*           Loop forward applying the L-transformations. */
00554 
00555             k = 1;
00556             kc = 1;
00557 L100:
00558             if (k > *n) {
00559                 goto L120;
00560             }
00561 
00562 /*           1 x 1 pivot block */
00563 
00564             if (ipiv[k] > 0) {
00565                 if (k < *n) {
00566 
00567 /*                 Interchange if P(K) != I. */
00568 
00569                     kp = ipiv[k];
00570                     if (kp != k) {
00571                         cswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1], 
00572                                 ldb);
00573                     }
00574 
00575 /*                 Apply the transformation */
00576 
00577                     i__1 = *n - k;
00578                     cgemv_("Transpose", &i__1, nrhs, &c_b1, &b[k + 1 + b_dim1]
00579 , ldb, &a[kc + 1], &c__1, &c_b1, &b[k + b_dim1], 
00580                             ldb);
00581                 }
00582                 if (nounit) {
00583                     cscal_(nrhs, &a[kc], &b[k + b_dim1], ldb);
00584                 }
00585                 kc = kc + *n - k + 1;
00586                 ++k;
00587 
00588 /*           2 x 2 pivot block. */
00589 
00590             } else {
00591                 kcnext = kc + *n - k + 1;
00592                 if (k < *n - 1) {
00593 
00594 /*              Interchange if P(K) != I. */
00595 
00596                     kp = (i__1 = ipiv[k], abs(i__1));
00597                     if (kp != k + 1) {
00598                         cswap_(nrhs, &b[k + 1 + b_dim1], ldb, &b[kp + b_dim1], 
00599                                  ldb);
00600                     }
00601 
00602 /*                 Apply the transformation */
00603 
00604                     i__1 = *n - k - 1;
00605                     cgemv_("Transpose", &i__1, nrhs, &c_b1, &b[k + 2 + b_dim1]
00606 , ldb, &a[kcnext + 1], &c__1, &c_b1, &b[k + 1 + 
00607                             b_dim1], ldb);
00608 
00609                     i__1 = *n - k - 1;
00610                     cgemv_("Transpose", &i__1, nrhs, &c_b1, &b[k + 2 + b_dim1]
00611 , ldb, &a[kc + 2], &c__1, &c_b1, &b[k + b_dim1], 
00612                             ldb);
00613                 }
00614 
00615 /*              Multiply by the diagonal block if non-unit. */
00616 
00617                 if (nounit) {
00618                     i__1 = kc;
00619                     d11.r = a[i__1].r, d11.i = a[i__1].i;
00620                     i__1 = kcnext;
00621                     d22.r = a[i__1].r, d22.i = a[i__1].i;
00622                     i__1 = kc + 1;
00623                     d21.r = a[i__1].r, d21.i = a[i__1].i;
00624                     d12.r = d21.r, d12.i = d21.i;
00625                     i__1 = *nrhs;
00626                     for (j = 1; j <= i__1; ++j) {
00627                         i__2 = k + j * b_dim1;
00628                         t1.r = b[i__2].r, t1.i = b[i__2].i;
00629                         i__2 = k + 1 + j * b_dim1;
00630                         t2.r = b[i__2].r, t2.i = b[i__2].i;
00631                         i__2 = k + j * b_dim1;
00632                         q__2.r = d11.r * t1.r - d11.i * t1.i, q__2.i = d11.r *
00633                                  t1.i + d11.i * t1.r;
00634                         q__3.r = d12.r * t2.r - d12.i * t2.i, q__3.i = d12.r *
00635                                  t2.i + d12.i * t2.r;
00636                         q__1.r = q__2.r + q__3.r, q__1.i = q__2.i + q__3.i;
00637                         b[i__2].r = q__1.r, b[i__2].i = q__1.i;
00638                         i__2 = k + 1 + j * b_dim1;
00639                         q__2.r = d21.r * t1.r - d21.i * t1.i, q__2.i = d21.r *
00640                                  t1.i + d21.i * t1.r;
00641                         q__3.r = d22.r * t2.r - d22.i * t2.i, q__3.i = d22.r *
00642                                  t2.i + d22.i * t2.r;
00643                         q__1.r = q__2.r + q__3.r, q__1.i = q__2.i + q__3.i;
00644                         b[i__2].r = q__1.r, b[i__2].i = q__1.i;
00645 /* L110: */
00646                     }
00647                 }
00648                 kc = kcnext + (*n - k);
00649                 k += 2;
00650             }
00651             goto L100;
00652 L120:
00653             ;
00654         }
00655 
00656     }
00657     return 0;
00658 
00659 /*     End of CLAVSP */
00660 
00661 } /* clavsp_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:55:32