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


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:56:42