ztrsm.c
Go to the documentation of this file.
00001 /* ztrsm.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 
00020 /* Subroutine */ int ztrsm_(char *side, char *uplo, char *transa, char *diag, 
00021         integer *m, integer *n, doublecomplex *alpha, doublecomplex *a, 
00022         integer *lda, doublecomplex *b, integer *ldb)
00023 {
00024     /* System generated locals */
00025     integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3, i__4, i__5, 
00026             i__6, i__7;
00027     doublecomplex z__1, z__2, z__3;
00028 
00029     /* Builtin functions */
00030     void z_div(doublecomplex *, doublecomplex *, doublecomplex *), d_cnjg(
00031             doublecomplex *, doublecomplex *);
00032 
00033     /* Local variables */
00034     integer i__, j, k, info;
00035     doublecomplex temp;
00036     logical lside;
00037     extern logical lsame_(char *, char *);
00038     integer nrowa;
00039     logical upper;
00040     extern /* Subroutine */ int xerbla_(char *, integer *);
00041     logical noconj, nounit;
00042 
00043 /*     .. Scalar Arguments .. */
00044 /*     .. */
00045 /*     .. Array Arguments .. */
00046 /*     .. */
00047 
00048 /*  Purpose */
00049 /*  ======= */
00050 
00051 /*  ZTRSM  solves one of the matrix equations */
00052 
00053 /*     op( A )*X = alpha*B,   or   X*op( A ) = alpha*B, */
00054 
00055 /*  where alpha is a scalar, X and B are m by n matrices, A is a unit, or */
00056 /*  non-unit,  upper or lower triangular matrix  and  op( A )  is one  of */
00057 
00058 /*     op( A ) = A   or   op( A ) = A'   or   op( A ) = conjg( A' ). */
00059 
00060 /*  The matrix X is overwritten on B. */
00061 
00062 /*  Arguments */
00063 /*  ========== */
00064 
00065 /*  SIDE   - CHARACTER*1. */
00066 /*           On entry, SIDE specifies whether op( A ) appears on the left */
00067 /*           or right of X as follows: */
00068 
00069 /*              SIDE = 'L' or 'l'   op( A )*X = alpha*B. */
00070 
00071 /*              SIDE = 'R' or 'r'   X*op( A ) = alpha*B. */
00072 
00073 /*           Unchanged on exit. */
00074 
00075 /*  UPLO   - CHARACTER*1. */
00076 /*           On entry, UPLO specifies whether the matrix A is an upper or */
00077 /*           lower triangular matrix as follows: */
00078 
00079 /*              UPLO = 'U' or 'u'   A is an upper triangular matrix. */
00080 
00081 /*              UPLO = 'L' or 'l'   A is a lower triangular matrix. */
00082 
00083 /*           Unchanged on exit. */
00084 
00085 /*  TRANSA - CHARACTER*1. */
00086 /*           On entry, TRANSA specifies the form of op( A ) to be used in */
00087 /*           the matrix multiplication as follows: */
00088 
00089 /*              TRANSA = 'N' or 'n'   op( A ) = A. */
00090 
00091 /*              TRANSA = 'T' or 't'   op( A ) = A'. */
00092 
00093 /*              TRANSA = 'C' or 'c'   op( A ) = conjg( A' ). */
00094 
00095 /*           Unchanged on exit. */
00096 
00097 /*  DIAG   - CHARACTER*1. */
00098 /*           On entry, DIAG specifies whether or not A is unit triangular */
00099 /*           as follows: */
00100 
00101 /*              DIAG = 'U' or 'u'   A is assumed to be unit triangular. */
00102 
00103 /*              DIAG = 'N' or 'n'   A is not assumed to be unit */
00104 /*                                  triangular. */
00105 
00106 /*           Unchanged on exit. */
00107 
00108 /*  M      - INTEGER. */
00109 /*           On entry, M specifies the number of rows of B. M must be at */
00110 /*           least zero. */
00111 /*           Unchanged on exit. */
00112 
00113 /*  N      - INTEGER. */
00114 /*           On entry, N specifies the number of columns of B.  N must be */
00115 /*           at least zero. */
00116 /*           Unchanged on exit. */
00117 
00118 /*  ALPHA  - COMPLEX*16      . */
00119 /*           On entry,  ALPHA specifies the scalar  alpha. When  alpha is */
00120 /*           zero then  A is not referenced and  B need not be set before */
00121 /*           entry. */
00122 /*           Unchanged on exit. */
00123 
00124 /*  A      - COMPLEX*16       array of DIMENSION ( LDA, k ), where k is m */
00125 /*           when  SIDE = 'L' or 'l'  and is  n  when  SIDE = 'R' or 'r'. */
00126 /*           Before entry  with  UPLO = 'U' or 'u',  the  leading  k by k */
00127 /*           upper triangular part of the array  A must contain the upper */
00128 /*           triangular matrix  and the strictly lower triangular part of */
00129 /*           A is not referenced. */
00130 /*           Before entry  with  UPLO = 'L' or 'l',  the  leading  k by k */
00131 /*           lower triangular part of the array  A must contain the lower */
00132 /*           triangular matrix  and the strictly upper triangular part of */
00133 /*           A is not referenced. */
00134 /*           Note that when  DIAG = 'U' or 'u',  the diagonal elements of */
00135 /*           A  are not referenced either,  but are assumed to be  unity. */
00136 /*           Unchanged on exit. */
00137 
00138 /*  LDA    - INTEGER. */
00139 /*           On entry, LDA specifies the first dimension of A as declared */
00140 /*           in the calling (sub) program.  When  SIDE = 'L' or 'l'  then */
00141 /*           LDA  must be at least  max( 1, m ),  when  SIDE = 'R' or 'r' */
00142 /*           then LDA must be at least max( 1, n ). */
00143 /*           Unchanged on exit. */
00144 
00145 /*  B      - COMPLEX*16       array of DIMENSION ( LDB, n ). */
00146 /*           Before entry,  the leading  m by n part of the array  B must */
00147 /*           contain  the  right-hand  side  matrix  B,  and  on exit  is */
00148 /*           overwritten by the solution matrix  X. */
00149 
00150 /*  LDB    - INTEGER. */
00151 /*           On entry, LDB specifies the first dimension of B as declared */
00152 /*           in  the  calling  (sub)  program.   LDB  must  be  at  least */
00153 /*           max( 1, m ). */
00154 /*           Unchanged on exit. */
00155 
00156 
00157 /*  Level 3 Blas routine. */
00158 
00159 /*  -- Written on 8-February-1989. */
00160 /*     Jack Dongarra, Argonne National Laboratory. */
00161 /*     Iain Duff, AERE Harwell. */
00162 /*     Jeremy Du Croz, Numerical Algorithms Group Ltd. */
00163 /*     Sven Hammarling, Numerical Algorithms Group Ltd. */
00164 
00165 
00166 /*     .. External Functions .. */
00167 /*     .. */
00168 /*     .. External Subroutines .. */
00169 /*     .. */
00170 /*     .. Intrinsic Functions .. */
00171 /*     .. */
00172 /*     .. Local Scalars .. */
00173 /*     .. */
00174 /*     .. Parameters .. */
00175 /*     .. */
00176 
00177 /*     Test the input parameters. */
00178 
00179     /* Parameter adjustments */
00180     a_dim1 = *lda;
00181     a_offset = 1 + a_dim1;
00182     a -= a_offset;
00183     b_dim1 = *ldb;
00184     b_offset = 1 + b_dim1;
00185     b -= b_offset;
00186 
00187     /* Function Body */
00188     lside = lsame_(side, "L");
00189     if (lside) {
00190         nrowa = *m;
00191     } else {
00192         nrowa = *n;
00193     }
00194     noconj = lsame_(transa, "T");
00195     nounit = lsame_(diag, "N");
00196     upper = lsame_(uplo, "U");
00197 
00198     info = 0;
00199     if (! lside && ! lsame_(side, "R")) {
00200         info = 1;
00201     } else if (! upper && ! lsame_(uplo, "L")) {
00202         info = 2;
00203     } else if (! lsame_(transa, "N") && ! lsame_(transa, 
00204              "T") && ! lsame_(transa, "C")) {
00205         info = 3;
00206     } else if (! lsame_(diag, "U") && ! lsame_(diag, 
00207             "N")) {
00208         info = 4;
00209     } else if (*m < 0) {
00210         info = 5;
00211     } else if (*n < 0) {
00212         info = 6;
00213     } else if (*lda < max(1,nrowa)) {
00214         info = 9;
00215     } else if (*ldb < max(1,*m)) {
00216         info = 11;
00217     }
00218     if (info != 0) {
00219         xerbla_("ZTRSM ", &info);
00220         return 0;
00221     }
00222 
00223 /*     Quick return if possible. */
00224 
00225     if (*m == 0 || *n == 0) {
00226         return 0;
00227     }
00228 
00229 /*     And when  alpha.eq.zero. */
00230 
00231     if (alpha->r == 0. && alpha->i == 0.) {
00232         i__1 = *n;
00233         for (j = 1; j <= i__1; ++j) {
00234             i__2 = *m;
00235             for (i__ = 1; i__ <= i__2; ++i__) {
00236                 i__3 = i__ + j * b_dim1;
00237                 b[i__3].r = 0., b[i__3].i = 0.;
00238 /* L10: */
00239             }
00240 /* L20: */
00241         }
00242         return 0;
00243     }
00244 
00245 /*     Start the operations. */
00246 
00247     if (lside) {
00248         if (lsame_(transa, "N")) {
00249 
00250 /*           Form  B := alpha*inv( A )*B. */
00251 
00252             if (upper) {
00253                 i__1 = *n;
00254                 for (j = 1; j <= i__1; ++j) {
00255                     if (alpha->r != 1. || alpha->i != 0.) {
00256                         i__2 = *m;
00257                         for (i__ = 1; i__ <= i__2; ++i__) {
00258                             i__3 = i__ + j * b_dim1;
00259                             i__4 = i__ + j * b_dim1;
00260                             z__1.r = alpha->r * b[i__4].r - alpha->i * b[i__4]
00261                                     .i, z__1.i = alpha->r * b[i__4].i + 
00262                                     alpha->i * b[i__4].r;
00263                             b[i__3].r = z__1.r, b[i__3].i = z__1.i;
00264 /* L30: */
00265                         }
00266                     }
00267                     for (k = *m; k >= 1; --k) {
00268                         i__2 = k + j * b_dim1;
00269                         if (b[i__2].r != 0. || b[i__2].i != 0.) {
00270                             if (nounit) {
00271                                 i__2 = k + j * b_dim1;
00272                                 z_div(&z__1, &b[k + j * b_dim1], &a[k + k * 
00273                                         a_dim1]);
00274                                 b[i__2].r = z__1.r, b[i__2].i = z__1.i;
00275                             }
00276                             i__2 = k - 1;
00277                             for (i__ = 1; i__ <= i__2; ++i__) {
00278                                 i__3 = i__ + j * b_dim1;
00279                                 i__4 = i__ + j * b_dim1;
00280                                 i__5 = k + j * b_dim1;
00281                                 i__6 = i__ + k * a_dim1;
00282                                 z__2.r = b[i__5].r * a[i__6].r - b[i__5].i * 
00283                                         a[i__6].i, z__2.i = b[i__5].r * a[
00284                                         i__6].i + b[i__5].i * a[i__6].r;
00285                                 z__1.r = b[i__4].r - z__2.r, z__1.i = b[i__4]
00286                                         .i - z__2.i;
00287                                 b[i__3].r = z__1.r, b[i__3].i = z__1.i;
00288 /* L40: */
00289                             }
00290                         }
00291 /* L50: */
00292                     }
00293 /* L60: */
00294                 }
00295             } else {
00296                 i__1 = *n;
00297                 for (j = 1; j <= i__1; ++j) {
00298                     if (alpha->r != 1. || alpha->i != 0.) {
00299                         i__2 = *m;
00300                         for (i__ = 1; i__ <= i__2; ++i__) {
00301                             i__3 = i__ + j * b_dim1;
00302                             i__4 = i__ + j * b_dim1;
00303                             z__1.r = alpha->r * b[i__4].r - alpha->i * b[i__4]
00304                                     .i, z__1.i = alpha->r * b[i__4].i + 
00305                                     alpha->i * b[i__4].r;
00306                             b[i__3].r = z__1.r, b[i__3].i = z__1.i;
00307 /* L70: */
00308                         }
00309                     }
00310                     i__2 = *m;
00311                     for (k = 1; k <= i__2; ++k) {
00312                         i__3 = k + j * b_dim1;
00313                         if (b[i__3].r != 0. || b[i__3].i != 0.) {
00314                             if (nounit) {
00315                                 i__3 = k + j * b_dim1;
00316                                 z_div(&z__1, &b[k + j * b_dim1], &a[k + k * 
00317                                         a_dim1]);
00318                                 b[i__3].r = z__1.r, b[i__3].i = z__1.i;
00319                             }
00320                             i__3 = *m;
00321                             for (i__ = k + 1; i__ <= i__3; ++i__) {
00322                                 i__4 = i__ + j * b_dim1;
00323                                 i__5 = i__ + j * b_dim1;
00324                                 i__6 = k + j * b_dim1;
00325                                 i__7 = i__ + k * a_dim1;
00326                                 z__2.r = b[i__6].r * a[i__7].r - b[i__6].i * 
00327                                         a[i__7].i, z__2.i = b[i__6].r * a[
00328                                         i__7].i + b[i__6].i * a[i__7].r;
00329                                 z__1.r = b[i__5].r - z__2.r, z__1.i = b[i__5]
00330                                         .i - z__2.i;
00331                                 b[i__4].r = z__1.r, b[i__4].i = z__1.i;
00332 /* L80: */
00333                             }
00334                         }
00335 /* L90: */
00336                     }
00337 /* L100: */
00338                 }
00339             }
00340         } else {
00341 
00342 /*           Form  B := alpha*inv( A' )*B */
00343 /*           or    B := alpha*inv( conjg( A' ) )*B. */
00344 
00345             if (upper) {
00346                 i__1 = *n;
00347                 for (j = 1; j <= i__1; ++j) {
00348                     i__2 = *m;
00349                     for (i__ = 1; i__ <= i__2; ++i__) {
00350                         i__3 = i__ + j * b_dim1;
00351                         z__1.r = alpha->r * b[i__3].r - alpha->i * b[i__3].i, 
00352                                 z__1.i = alpha->r * b[i__3].i + alpha->i * b[
00353                                 i__3].r;
00354                         temp.r = z__1.r, temp.i = z__1.i;
00355                         if (noconj) {
00356                             i__3 = i__ - 1;
00357                             for (k = 1; k <= i__3; ++k) {
00358                                 i__4 = k + i__ * a_dim1;
00359                                 i__5 = k + j * b_dim1;
00360                                 z__2.r = a[i__4].r * b[i__5].r - a[i__4].i * 
00361                                         b[i__5].i, z__2.i = a[i__4].r * b[
00362                                         i__5].i + a[i__4].i * b[i__5].r;
00363                                 z__1.r = temp.r - z__2.r, z__1.i = temp.i - 
00364                                         z__2.i;
00365                                 temp.r = z__1.r, temp.i = z__1.i;
00366 /* L110: */
00367                             }
00368                             if (nounit) {
00369                                 z_div(&z__1, &temp, &a[i__ + i__ * a_dim1]);
00370                                 temp.r = z__1.r, temp.i = z__1.i;
00371                             }
00372                         } else {
00373                             i__3 = i__ - 1;
00374                             for (k = 1; k <= i__3; ++k) {
00375                                 d_cnjg(&z__3, &a[k + i__ * a_dim1]);
00376                                 i__4 = k + j * b_dim1;
00377                                 z__2.r = z__3.r * b[i__4].r - z__3.i * b[i__4]
00378                                         .i, z__2.i = z__3.r * b[i__4].i + 
00379                                         z__3.i * b[i__4].r;
00380                                 z__1.r = temp.r - z__2.r, z__1.i = temp.i - 
00381                                         z__2.i;
00382                                 temp.r = z__1.r, temp.i = z__1.i;
00383 /* L120: */
00384                             }
00385                             if (nounit) {
00386                                 d_cnjg(&z__2, &a[i__ + i__ * a_dim1]);
00387                                 z_div(&z__1, &temp, &z__2);
00388                                 temp.r = z__1.r, temp.i = z__1.i;
00389                             }
00390                         }
00391                         i__3 = i__ + j * b_dim1;
00392                         b[i__3].r = temp.r, b[i__3].i = temp.i;
00393 /* L130: */
00394                     }
00395 /* L140: */
00396                 }
00397             } else {
00398                 i__1 = *n;
00399                 for (j = 1; j <= i__1; ++j) {
00400                     for (i__ = *m; i__ >= 1; --i__) {
00401                         i__2 = i__ + j * b_dim1;
00402                         z__1.r = alpha->r * b[i__2].r - alpha->i * b[i__2].i, 
00403                                 z__1.i = alpha->r * b[i__2].i + alpha->i * b[
00404                                 i__2].r;
00405                         temp.r = z__1.r, temp.i = z__1.i;
00406                         if (noconj) {
00407                             i__2 = *m;
00408                             for (k = i__ + 1; k <= i__2; ++k) {
00409                                 i__3 = k + i__ * a_dim1;
00410                                 i__4 = k + j * b_dim1;
00411                                 z__2.r = a[i__3].r * b[i__4].r - a[i__3].i * 
00412                                         b[i__4].i, z__2.i = a[i__3].r * b[
00413                                         i__4].i + a[i__3].i * b[i__4].r;
00414                                 z__1.r = temp.r - z__2.r, z__1.i = temp.i - 
00415                                         z__2.i;
00416                                 temp.r = z__1.r, temp.i = z__1.i;
00417 /* L150: */
00418                             }
00419                             if (nounit) {
00420                                 z_div(&z__1, &temp, &a[i__ + i__ * a_dim1]);
00421                                 temp.r = z__1.r, temp.i = z__1.i;
00422                             }
00423                         } else {
00424                             i__2 = *m;
00425                             for (k = i__ + 1; k <= i__2; ++k) {
00426                                 d_cnjg(&z__3, &a[k + i__ * a_dim1]);
00427                                 i__3 = k + j * b_dim1;
00428                                 z__2.r = z__3.r * b[i__3].r - z__3.i * b[i__3]
00429                                         .i, z__2.i = z__3.r * b[i__3].i + 
00430                                         z__3.i * b[i__3].r;
00431                                 z__1.r = temp.r - z__2.r, z__1.i = temp.i - 
00432                                         z__2.i;
00433                                 temp.r = z__1.r, temp.i = z__1.i;
00434 /* L160: */
00435                             }
00436                             if (nounit) {
00437                                 d_cnjg(&z__2, &a[i__ + i__ * a_dim1]);
00438                                 z_div(&z__1, &temp, &z__2);
00439                                 temp.r = z__1.r, temp.i = z__1.i;
00440                             }
00441                         }
00442                         i__2 = i__ + j * b_dim1;
00443                         b[i__2].r = temp.r, b[i__2].i = temp.i;
00444 /* L170: */
00445                     }
00446 /* L180: */
00447                 }
00448             }
00449         }
00450     } else {
00451         if (lsame_(transa, "N")) {
00452 
00453 /*           Form  B := alpha*B*inv( A ). */
00454 
00455             if (upper) {
00456                 i__1 = *n;
00457                 for (j = 1; j <= i__1; ++j) {
00458                     if (alpha->r != 1. || alpha->i != 0.) {
00459                         i__2 = *m;
00460                         for (i__ = 1; i__ <= i__2; ++i__) {
00461                             i__3 = i__ + j * b_dim1;
00462                             i__4 = i__ + j * b_dim1;
00463                             z__1.r = alpha->r * b[i__4].r - alpha->i * b[i__4]
00464                                     .i, z__1.i = alpha->r * b[i__4].i + 
00465                                     alpha->i * b[i__4].r;
00466                             b[i__3].r = z__1.r, b[i__3].i = z__1.i;
00467 /* L190: */
00468                         }
00469                     }
00470                     i__2 = j - 1;
00471                     for (k = 1; k <= i__2; ++k) {
00472                         i__3 = k + j * a_dim1;
00473                         if (a[i__3].r != 0. || a[i__3].i != 0.) {
00474                             i__3 = *m;
00475                             for (i__ = 1; i__ <= i__3; ++i__) {
00476                                 i__4 = i__ + j * b_dim1;
00477                                 i__5 = i__ + j * b_dim1;
00478                                 i__6 = k + j * a_dim1;
00479                                 i__7 = i__ + k * b_dim1;
00480                                 z__2.r = a[i__6].r * b[i__7].r - a[i__6].i * 
00481                                         b[i__7].i, z__2.i = a[i__6].r * b[
00482                                         i__7].i + a[i__6].i * b[i__7].r;
00483                                 z__1.r = b[i__5].r - z__2.r, z__1.i = b[i__5]
00484                                         .i - z__2.i;
00485                                 b[i__4].r = z__1.r, b[i__4].i = z__1.i;
00486 /* L200: */
00487                             }
00488                         }
00489 /* L210: */
00490                     }
00491                     if (nounit) {
00492                         z_div(&z__1, &c_b1, &a[j + j * a_dim1]);
00493                         temp.r = z__1.r, temp.i = z__1.i;
00494                         i__2 = *m;
00495                         for (i__ = 1; i__ <= i__2; ++i__) {
00496                             i__3 = i__ + j * b_dim1;
00497                             i__4 = i__ + j * b_dim1;
00498                             z__1.r = temp.r * b[i__4].r - temp.i * b[i__4].i, 
00499                                     z__1.i = temp.r * b[i__4].i + temp.i * b[
00500                                     i__4].r;
00501                             b[i__3].r = z__1.r, b[i__3].i = z__1.i;
00502 /* L220: */
00503                         }
00504                     }
00505 /* L230: */
00506                 }
00507             } else {
00508                 for (j = *n; j >= 1; --j) {
00509                     if (alpha->r != 1. || alpha->i != 0.) {
00510                         i__1 = *m;
00511                         for (i__ = 1; i__ <= i__1; ++i__) {
00512                             i__2 = i__ + j * b_dim1;
00513                             i__3 = i__ + j * b_dim1;
00514                             z__1.r = alpha->r * b[i__3].r - alpha->i * b[i__3]
00515                                     .i, z__1.i = alpha->r * b[i__3].i + 
00516                                     alpha->i * b[i__3].r;
00517                             b[i__2].r = z__1.r, b[i__2].i = z__1.i;
00518 /* L240: */
00519                         }
00520                     }
00521                     i__1 = *n;
00522                     for (k = j + 1; k <= i__1; ++k) {
00523                         i__2 = k + j * a_dim1;
00524                         if (a[i__2].r != 0. || a[i__2].i != 0.) {
00525                             i__2 = *m;
00526                             for (i__ = 1; i__ <= i__2; ++i__) {
00527                                 i__3 = i__ + j * b_dim1;
00528                                 i__4 = i__ + j * b_dim1;
00529                                 i__5 = k + j * a_dim1;
00530                                 i__6 = i__ + k * b_dim1;
00531                                 z__2.r = a[i__5].r * b[i__6].r - a[i__5].i * 
00532                                         b[i__6].i, z__2.i = a[i__5].r * b[
00533                                         i__6].i + a[i__5].i * b[i__6].r;
00534                                 z__1.r = b[i__4].r - z__2.r, z__1.i = b[i__4]
00535                                         .i - z__2.i;
00536                                 b[i__3].r = z__1.r, b[i__3].i = z__1.i;
00537 /* L250: */
00538                             }
00539                         }
00540 /* L260: */
00541                     }
00542                     if (nounit) {
00543                         z_div(&z__1, &c_b1, &a[j + j * a_dim1]);
00544                         temp.r = z__1.r, temp.i = z__1.i;
00545                         i__1 = *m;
00546                         for (i__ = 1; i__ <= i__1; ++i__) {
00547                             i__2 = i__ + j * b_dim1;
00548                             i__3 = i__ + j * b_dim1;
00549                             z__1.r = temp.r * b[i__3].r - temp.i * b[i__3].i, 
00550                                     z__1.i = temp.r * b[i__3].i + temp.i * b[
00551                                     i__3].r;
00552                             b[i__2].r = z__1.r, b[i__2].i = z__1.i;
00553 /* L270: */
00554                         }
00555                     }
00556 /* L280: */
00557                 }
00558             }
00559         } else {
00560 
00561 /*           Form  B := alpha*B*inv( A' ) */
00562 /*           or    B := alpha*B*inv( conjg( A' ) ). */
00563 
00564             if (upper) {
00565                 for (k = *n; k >= 1; --k) {
00566                     if (nounit) {
00567                         if (noconj) {
00568                             z_div(&z__1, &c_b1, &a[k + k * a_dim1]);
00569                             temp.r = z__1.r, temp.i = z__1.i;
00570                         } else {
00571                             d_cnjg(&z__2, &a[k + k * a_dim1]);
00572                             z_div(&z__1, &c_b1, &z__2);
00573                             temp.r = z__1.r, temp.i = z__1.i;
00574                         }
00575                         i__1 = *m;
00576                         for (i__ = 1; i__ <= i__1; ++i__) {
00577                             i__2 = i__ + k * b_dim1;
00578                             i__3 = i__ + k * b_dim1;
00579                             z__1.r = temp.r * b[i__3].r - temp.i * b[i__3].i, 
00580                                     z__1.i = temp.r * b[i__3].i + temp.i * b[
00581                                     i__3].r;
00582                             b[i__2].r = z__1.r, b[i__2].i = z__1.i;
00583 /* L290: */
00584                         }
00585                     }
00586                     i__1 = k - 1;
00587                     for (j = 1; j <= i__1; ++j) {
00588                         i__2 = j + k * a_dim1;
00589                         if (a[i__2].r != 0. || a[i__2].i != 0.) {
00590                             if (noconj) {
00591                                 i__2 = j + k * a_dim1;
00592                                 temp.r = a[i__2].r, temp.i = a[i__2].i;
00593                             } else {
00594                                 d_cnjg(&z__1, &a[j + k * a_dim1]);
00595                                 temp.r = z__1.r, temp.i = z__1.i;
00596                             }
00597                             i__2 = *m;
00598                             for (i__ = 1; i__ <= i__2; ++i__) {
00599                                 i__3 = i__ + j * b_dim1;
00600                                 i__4 = i__ + j * b_dim1;
00601                                 i__5 = i__ + k * b_dim1;
00602                                 z__2.r = temp.r * b[i__5].r - temp.i * b[i__5]
00603                                         .i, z__2.i = temp.r * b[i__5].i + 
00604                                         temp.i * b[i__5].r;
00605                                 z__1.r = b[i__4].r - z__2.r, z__1.i = b[i__4]
00606                                         .i - z__2.i;
00607                                 b[i__3].r = z__1.r, b[i__3].i = z__1.i;
00608 /* L300: */
00609                             }
00610                         }
00611 /* L310: */
00612                     }
00613                     if (alpha->r != 1. || alpha->i != 0.) {
00614                         i__1 = *m;
00615                         for (i__ = 1; i__ <= i__1; ++i__) {
00616                             i__2 = i__ + k * b_dim1;
00617                             i__3 = i__ + k * b_dim1;
00618                             z__1.r = alpha->r * b[i__3].r - alpha->i * b[i__3]
00619                                     .i, z__1.i = alpha->r * b[i__3].i + 
00620                                     alpha->i * b[i__3].r;
00621                             b[i__2].r = z__1.r, b[i__2].i = z__1.i;
00622 /* L320: */
00623                         }
00624                     }
00625 /* L330: */
00626                 }
00627             } else {
00628                 i__1 = *n;
00629                 for (k = 1; k <= i__1; ++k) {
00630                     if (nounit) {
00631                         if (noconj) {
00632                             z_div(&z__1, &c_b1, &a[k + k * a_dim1]);
00633                             temp.r = z__1.r, temp.i = z__1.i;
00634                         } else {
00635                             d_cnjg(&z__2, &a[k + k * a_dim1]);
00636                             z_div(&z__1, &c_b1, &z__2);
00637                             temp.r = z__1.r, temp.i = z__1.i;
00638                         }
00639                         i__2 = *m;
00640                         for (i__ = 1; i__ <= i__2; ++i__) {
00641                             i__3 = i__ + k * b_dim1;
00642                             i__4 = i__ + k * b_dim1;
00643                             z__1.r = temp.r * b[i__4].r - temp.i * b[i__4].i, 
00644                                     z__1.i = temp.r * b[i__4].i + temp.i * b[
00645                                     i__4].r;
00646                             b[i__3].r = z__1.r, b[i__3].i = z__1.i;
00647 /* L340: */
00648                         }
00649                     }
00650                     i__2 = *n;
00651                     for (j = k + 1; j <= i__2; ++j) {
00652                         i__3 = j + k * a_dim1;
00653                         if (a[i__3].r != 0. || a[i__3].i != 0.) {
00654                             if (noconj) {
00655                                 i__3 = j + k * a_dim1;
00656                                 temp.r = a[i__3].r, temp.i = a[i__3].i;
00657                             } else {
00658                                 d_cnjg(&z__1, &a[j + k * a_dim1]);
00659                                 temp.r = z__1.r, temp.i = z__1.i;
00660                             }
00661                             i__3 = *m;
00662                             for (i__ = 1; i__ <= i__3; ++i__) {
00663                                 i__4 = i__ + j * b_dim1;
00664                                 i__5 = i__ + j * b_dim1;
00665                                 i__6 = i__ + k * b_dim1;
00666                                 z__2.r = temp.r * b[i__6].r - temp.i * b[i__6]
00667                                         .i, z__2.i = temp.r * b[i__6].i + 
00668                                         temp.i * b[i__6].r;
00669                                 z__1.r = b[i__5].r - z__2.r, z__1.i = b[i__5]
00670                                         .i - z__2.i;
00671                                 b[i__4].r = z__1.r, b[i__4].i = z__1.i;
00672 /* L350: */
00673                             }
00674                         }
00675 /* L360: */
00676                     }
00677                     if (alpha->r != 1. || alpha->i != 0.) {
00678                         i__2 = *m;
00679                         for (i__ = 1; i__ <= i__2; ++i__) {
00680                             i__3 = i__ + k * b_dim1;
00681                             i__4 = i__ + k * b_dim1;
00682                             z__1.r = alpha->r * b[i__4].r - alpha->i * b[i__4]
00683                                     .i, z__1.i = alpha->r * b[i__4].i + 
00684                                     alpha->i * b[i__4].r;
00685                             b[i__3].r = z__1.r, b[i__3].i = z__1.i;
00686 /* L370: */
00687                         }
00688                     }
00689 /* L380: */
00690                 }
00691             }
00692         }
00693     }
00694 
00695     return 0;
00696 
00697 /*     End of ZTRSM . */
00698 
00699 } /* ztrsm_ */


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