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


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