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


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