strsm.c
Go to the documentation of this file.
00001 /* strsm.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 strsm_(char *side, char *uplo, char *transa, char *diag, 
00017         integer *m, integer *n, real *alpha, real *a, integer *lda, real *b, 
00018         integer *ldb)
00019 {
00020     /* System generated locals */
00021     integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3;
00022 
00023     /* Local variables */
00024     integer i__, j, k, info;
00025     real temp;
00026     logical lside;
00027     extern logical lsame_(char *, char *);
00028     integer nrowa;
00029     logical upper;
00030     extern /* Subroutine */ int xerbla_(char *, integer *);
00031     logical nounit;
00032 
00033 /*     .. Scalar Arguments .. */
00034 /*     .. */
00035 /*     .. Array Arguments .. */
00036 /*     .. */
00037 
00038 /*  Purpose */
00039 /*  ======= */
00040 
00041 /*  STRSM  solves one of the matrix equations */
00042 
00043 /*     op( A )*X = alpha*B,   or   X*op( A ) = alpha*B, */
00044 
00045 /*  where alpha is a scalar, X and B are m by n matrices, A is a unit, or */
00046 /*  non-unit,  upper or lower triangular matrix  and  op( A )  is one  of */
00047 
00048 /*     op( A ) = A   or   op( A ) = A'. */
00049 
00050 /*  The matrix X is overwritten on B. */
00051 
00052 /*  Arguments */
00053 /*  ========== */
00054 
00055 /*  SIDE   - CHARACTER*1. */
00056 /*           On entry, SIDE specifies whether op( A ) appears on the left */
00057 /*           or right of X as follows: */
00058 
00059 /*              SIDE = 'L' or 'l'   op( A )*X = alpha*B. */
00060 
00061 /*              SIDE = 'R' or 'r'   X*op( A ) = alpha*B. */
00062 
00063 /*           Unchanged on exit. */
00064 
00065 /*  UPLO   - CHARACTER*1. */
00066 /*           On entry, UPLO specifies whether the matrix A is an upper or */
00067 /*           lower triangular matrix as follows: */
00068 
00069 /*              UPLO = 'U' or 'u'   A is an upper triangular matrix. */
00070 
00071 /*              UPLO = 'L' or 'l'   A is a lower triangular matrix. */
00072 
00073 /*           Unchanged on exit. */
00074 
00075 /*  TRANSA - CHARACTER*1. */
00076 /*           On entry, TRANSA specifies the form of op( A ) to be used in */
00077 /*           the matrix multiplication as follows: */
00078 
00079 /*              TRANSA = 'N' or 'n'   op( A ) = A. */
00080 
00081 /*              TRANSA = 'T' or 't'   op( A ) = A'. */
00082 
00083 /*              TRANSA = 'C' or 'c'   op( A ) = A'. */
00084 
00085 /*           Unchanged on exit. */
00086 
00087 /*  DIAG   - CHARACTER*1. */
00088 /*           On entry, DIAG specifies whether or not A is unit triangular */
00089 /*           as follows: */
00090 
00091 /*              DIAG = 'U' or 'u'   A is assumed to be unit triangular. */
00092 
00093 /*              DIAG = 'N' or 'n'   A is not assumed to be unit */
00094 /*                                  triangular. */
00095 
00096 /*           Unchanged on exit. */
00097 
00098 /*  M      - INTEGER. */
00099 /*           On entry, M specifies the number of rows of B. M must be at */
00100 /*           least zero. */
00101 /*           Unchanged on exit. */
00102 
00103 /*  N      - INTEGER. */
00104 /*           On entry, N specifies the number of columns of B.  N must be */
00105 /*           at least zero. */
00106 /*           Unchanged on exit. */
00107 
00108 /*  ALPHA  - REAL            . */
00109 /*           On entry,  ALPHA specifies the scalar  alpha. When  alpha is */
00110 /*           zero then  A is not referenced and  B need not be set before */
00111 /*           entry. */
00112 /*           Unchanged on exit. */
00113 
00114 /*  A      - REAL             array of DIMENSION ( LDA, k ), where k is m */
00115 /*           when  SIDE = 'L' or 'l'  and is  n  when  SIDE = 'R' or 'r'. */
00116 /*           Before entry  with  UPLO = 'U' or 'u',  the  leading  k by k */
00117 /*           upper triangular part of the array  A must contain the upper */
00118 /*           triangular matrix  and the strictly lower triangular part of */
00119 /*           A is not referenced. */
00120 /*           Before entry  with  UPLO = 'L' or 'l',  the  leading  k by k */
00121 /*           lower triangular part of the array  A must contain the lower */
00122 /*           triangular matrix  and the strictly upper triangular part of */
00123 /*           A is not referenced. */
00124 /*           Note that when  DIAG = 'U' or 'u',  the diagonal elements of */
00125 /*           A  are not referenced either,  but are assumed to be  unity. */
00126 /*           Unchanged on exit. */
00127 
00128 /*  LDA    - INTEGER. */
00129 /*           On entry, LDA specifies the first dimension of A as declared */
00130 /*           in the calling (sub) program.  When  SIDE = 'L' or 'l'  then */
00131 /*           LDA  must be at least  max( 1, m ),  when  SIDE = 'R' or 'r' */
00132 /*           then LDA must be at least max( 1, n ). */
00133 /*           Unchanged on exit. */
00134 
00135 /*  B      - REAL             array of DIMENSION ( LDB, n ). */
00136 /*           Before entry,  the leading  m by n part of the array  B must */
00137 /*           contain  the  right-hand  side  matrix  B,  and  on exit  is */
00138 /*           overwritten by the solution matrix  X. */
00139 
00140 /*  LDB    - INTEGER. */
00141 /*           On entry, LDB specifies the first dimension of B as declared */
00142 /*           in  the  calling  (sub)  program.   LDB  must  be  at  least */
00143 /*           max( 1, m ). */
00144 /*           Unchanged on exit. */
00145 
00146 
00147 /*  Level 3 Blas routine. */
00148 
00149 
00150 /*  -- Written on 8-February-1989. */
00151 /*     Jack Dongarra, Argonne National Laboratory. */
00152 /*     Iain Duff, AERE Harwell. */
00153 /*     Jeremy Du Croz, Numerical Algorithms Group Ltd. */
00154 /*     Sven Hammarling, Numerical Algorithms Group Ltd. */
00155 
00156 
00157 /*     .. External Functions .. */
00158 /*     .. */
00159 /*     .. External Subroutines .. */
00160 /*     .. */
00161 /*     .. Intrinsic Functions .. */
00162 /*     .. */
00163 /*     .. Local Scalars .. */
00164 /*     .. */
00165 /*     .. Parameters .. */
00166 /*     .. */
00167 
00168 /*     Test the input parameters. */
00169 
00170     /* Parameter adjustments */
00171     a_dim1 = *lda;
00172     a_offset = 1 + a_dim1;
00173     a -= a_offset;
00174     b_dim1 = *ldb;
00175     b_offset = 1 + b_dim1;
00176     b -= b_offset;
00177 
00178     /* Function Body */
00179     lside = lsame_(side, "L");
00180     if (lside) {
00181         nrowa = *m;
00182     } else {
00183         nrowa = *n;
00184     }
00185     nounit = lsame_(diag, "N");
00186     upper = lsame_(uplo, "U");
00187 
00188     info = 0;
00189     if (! lside && ! lsame_(side, "R")) {
00190         info = 1;
00191     } else if (! upper && ! lsame_(uplo, "L")) {
00192         info = 2;
00193     } else if (! lsame_(transa, "N") && ! lsame_(transa, 
00194              "T") && ! lsame_(transa, "C")) {
00195         info = 3;
00196     } else if (! lsame_(diag, "U") && ! lsame_(diag, 
00197             "N")) {
00198         info = 4;
00199     } else if (*m < 0) {
00200         info = 5;
00201     } else if (*n < 0) {
00202         info = 6;
00203     } else if (*lda < max(1,nrowa)) {
00204         info = 9;
00205     } else if (*ldb < max(1,*m)) {
00206         info = 11;
00207     }
00208     if (info != 0) {
00209         xerbla_("STRSM ", &info);
00210         return 0;
00211     }
00212 
00213 /*     Quick return if possible. */
00214 
00215     if (*m == 0 || *n == 0) {
00216         return 0;
00217     }
00218 
00219 /*     And when  alpha.eq.zero. */
00220 
00221     if (*alpha == 0.f) {
00222         i__1 = *n;
00223         for (j = 1; j <= i__1; ++j) {
00224             i__2 = *m;
00225             for (i__ = 1; i__ <= i__2; ++i__) {
00226                 b[i__ + j * b_dim1] = 0.f;
00227 /* L10: */
00228             }
00229 /* L20: */
00230         }
00231         return 0;
00232     }
00233 
00234 /*     Start the operations. */
00235 
00236     if (lside) {
00237         if (lsame_(transa, "N")) {
00238 
00239 /*           Form  B := alpha*inv( A )*B. */
00240 
00241             if (upper) {
00242                 i__1 = *n;
00243                 for (j = 1; j <= i__1; ++j) {
00244                     if (*alpha != 1.f) {
00245                         i__2 = *m;
00246                         for (i__ = 1; i__ <= i__2; ++i__) {
00247                             b[i__ + j * b_dim1] = *alpha * b[i__ + j * b_dim1]
00248                                     ;
00249 /* L30: */
00250                         }
00251                     }
00252                     for (k = *m; k >= 1; --k) {
00253                         if (b[k + j * b_dim1] != 0.f) {
00254                             if (nounit) {
00255                                 b[k + j * b_dim1] /= a[k + k * a_dim1];
00256                             }
00257                             i__2 = k - 1;
00258                             for (i__ = 1; i__ <= i__2; ++i__) {
00259                                 b[i__ + j * b_dim1] -= b[k + j * b_dim1] * a[
00260                                         i__ + k * a_dim1];
00261 /* L40: */
00262                             }
00263                         }
00264 /* L50: */
00265                     }
00266 /* L60: */
00267                 }
00268             } else {
00269                 i__1 = *n;
00270                 for (j = 1; j <= i__1; ++j) {
00271                     if (*alpha != 1.f) {
00272                         i__2 = *m;
00273                         for (i__ = 1; i__ <= i__2; ++i__) {
00274                             b[i__ + j * b_dim1] = *alpha * b[i__ + j * b_dim1]
00275                                     ;
00276 /* L70: */
00277                         }
00278                     }
00279                     i__2 = *m;
00280                     for (k = 1; k <= i__2; ++k) {
00281                         if (b[k + j * b_dim1] != 0.f) {
00282                             if (nounit) {
00283                                 b[k + j * b_dim1] /= a[k + k * a_dim1];
00284                             }
00285                             i__3 = *m;
00286                             for (i__ = k + 1; i__ <= i__3; ++i__) {
00287                                 b[i__ + j * b_dim1] -= b[k + j * b_dim1] * a[
00288                                         i__ + k * a_dim1];
00289 /* L80: */
00290                             }
00291                         }
00292 /* L90: */
00293                     }
00294 /* L100: */
00295                 }
00296             }
00297         } else {
00298 
00299 /*           Form  B := alpha*inv( A' )*B. */
00300 
00301             if (upper) {
00302                 i__1 = *n;
00303                 for (j = 1; j <= i__1; ++j) {
00304                     i__2 = *m;
00305                     for (i__ = 1; i__ <= i__2; ++i__) {
00306                         temp = *alpha * b[i__ + j * b_dim1];
00307                         i__3 = i__ - 1;
00308                         for (k = 1; k <= i__3; ++k) {
00309                             temp -= a[k + i__ * a_dim1] * b[k + j * b_dim1];
00310 /* L110: */
00311                         }
00312                         if (nounit) {
00313                             temp /= a[i__ + i__ * a_dim1];
00314                         }
00315                         b[i__ + j * b_dim1] = temp;
00316 /* L120: */
00317                     }
00318 /* L130: */
00319                 }
00320             } else {
00321                 i__1 = *n;
00322                 for (j = 1; j <= i__1; ++j) {
00323                     for (i__ = *m; i__ >= 1; --i__) {
00324                         temp = *alpha * b[i__ + j * b_dim1];
00325                         i__2 = *m;
00326                         for (k = i__ + 1; k <= i__2; ++k) {
00327                             temp -= a[k + i__ * a_dim1] * b[k + j * b_dim1];
00328 /* L140: */
00329                         }
00330                         if (nounit) {
00331                             temp /= a[i__ + i__ * a_dim1];
00332                         }
00333                         b[i__ + j * b_dim1] = temp;
00334 /* L150: */
00335                     }
00336 /* L160: */
00337                 }
00338             }
00339         }
00340     } else {
00341         if (lsame_(transa, "N")) {
00342 
00343 /*           Form  B := alpha*B*inv( A ). */
00344 
00345             if (upper) {
00346                 i__1 = *n;
00347                 for (j = 1; j <= i__1; ++j) {
00348                     if (*alpha != 1.f) {
00349                         i__2 = *m;
00350                         for (i__ = 1; i__ <= i__2; ++i__) {
00351                             b[i__ + j * b_dim1] = *alpha * b[i__ + j * b_dim1]
00352                                     ;
00353 /* L170: */
00354                         }
00355                     }
00356                     i__2 = j - 1;
00357                     for (k = 1; k <= i__2; ++k) {
00358                         if (a[k + j * a_dim1] != 0.f) {
00359                             i__3 = *m;
00360                             for (i__ = 1; i__ <= i__3; ++i__) {
00361                                 b[i__ + j * b_dim1] -= a[k + j * a_dim1] * b[
00362                                         i__ + k * b_dim1];
00363 /* L180: */
00364                             }
00365                         }
00366 /* L190: */
00367                     }
00368                     if (nounit) {
00369                         temp = 1.f / a[j + j * a_dim1];
00370                         i__2 = *m;
00371                         for (i__ = 1; i__ <= i__2; ++i__) {
00372                             b[i__ + j * b_dim1] = temp * b[i__ + j * b_dim1];
00373 /* L200: */
00374                         }
00375                     }
00376 /* L210: */
00377                 }
00378             } else {
00379                 for (j = *n; j >= 1; --j) {
00380                     if (*alpha != 1.f) {
00381                         i__1 = *m;
00382                         for (i__ = 1; i__ <= i__1; ++i__) {
00383                             b[i__ + j * b_dim1] = *alpha * b[i__ + j * b_dim1]
00384                                     ;
00385 /* L220: */
00386                         }
00387                     }
00388                     i__1 = *n;
00389                     for (k = j + 1; k <= i__1; ++k) {
00390                         if (a[k + j * a_dim1] != 0.f) {
00391                             i__2 = *m;
00392                             for (i__ = 1; i__ <= i__2; ++i__) {
00393                                 b[i__ + j * b_dim1] -= a[k + j * a_dim1] * b[
00394                                         i__ + k * b_dim1];
00395 /* L230: */
00396                             }
00397                         }
00398 /* L240: */
00399                     }
00400                     if (nounit) {
00401                         temp = 1.f / a[j + j * a_dim1];
00402                         i__1 = *m;
00403                         for (i__ = 1; i__ <= i__1; ++i__) {
00404                             b[i__ + j * b_dim1] = temp * b[i__ + j * b_dim1];
00405 /* L250: */
00406                         }
00407                     }
00408 /* L260: */
00409                 }
00410             }
00411         } else {
00412 
00413 /*           Form  B := alpha*B*inv( A' ). */
00414 
00415             if (upper) {
00416                 for (k = *n; k >= 1; --k) {
00417                     if (nounit) {
00418                         temp = 1.f / a[k + k * a_dim1];
00419                         i__1 = *m;
00420                         for (i__ = 1; i__ <= i__1; ++i__) {
00421                             b[i__ + k * b_dim1] = temp * b[i__ + k * b_dim1];
00422 /* L270: */
00423                         }
00424                     }
00425                     i__1 = k - 1;
00426                     for (j = 1; j <= i__1; ++j) {
00427                         if (a[j + k * a_dim1] != 0.f) {
00428                             temp = a[j + k * a_dim1];
00429                             i__2 = *m;
00430                             for (i__ = 1; i__ <= i__2; ++i__) {
00431                                 b[i__ + j * b_dim1] -= temp * b[i__ + k * 
00432                                         b_dim1];
00433 /* L280: */
00434                             }
00435                         }
00436 /* L290: */
00437                     }
00438                     if (*alpha != 1.f) {
00439                         i__1 = *m;
00440                         for (i__ = 1; i__ <= i__1; ++i__) {
00441                             b[i__ + k * b_dim1] = *alpha * b[i__ + k * b_dim1]
00442                                     ;
00443 /* L300: */
00444                         }
00445                     }
00446 /* L310: */
00447                 }
00448             } else {
00449                 i__1 = *n;
00450                 for (k = 1; k <= i__1; ++k) {
00451                     if (nounit) {
00452                         temp = 1.f / a[k + k * a_dim1];
00453                         i__2 = *m;
00454                         for (i__ = 1; i__ <= i__2; ++i__) {
00455                             b[i__ + k * b_dim1] = temp * b[i__ + k * b_dim1];
00456 /* L320: */
00457                         }
00458                     }
00459                     i__2 = *n;
00460                     for (j = k + 1; j <= i__2; ++j) {
00461                         if (a[j + k * a_dim1] != 0.f) {
00462                             temp = a[j + k * a_dim1];
00463                             i__3 = *m;
00464                             for (i__ = 1; i__ <= i__3; ++i__) {
00465                                 b[i__ + j * b_dim1] -= temp * b[i__ + k * 
00466                                         b_dim1];
00467 /* L330: */
00468                             }
00469                         }
00470 /* L340: */
00471                     }
00472                     if (*alpha != 1.f) {
00473                         i__2 = *m;
00474                         for (i__ = 1; i__ <= i__2; ++i__) {
00475                             b[i__ + k * b_dim1] = *alpha * b[i__ + k * b_dim1]
00476                                     ;
00477 /* L350: */
00478                         }
00479                     }
00480 /* L360: */
00481                 }
00482             }
00483         }
00484     }
00485 
00486     return 0;
00487 
00488 /*     End of STRSM . */
00489 
00490 } /* strsm_ */

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