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


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