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


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