stfsm.c
Go to the documentation of this file.
00001 /* stfsm.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 real c_b23 = -1.f;
00019 static real c_b27 = 1.f;
00020 
00021 /* Subroutine */ int stfsm_(char *transr, char *side, char *uplo, char *trans, 
00022          char *diag, integer *m, integer *n, real *alpha, real *a, real *b, 
00023         integer *ldb)
00024 {
00025     /* System generated locals */
00026     integer b_dim1, b_offset, i__1, i__2;
00027 
00028     /* Local variables */
00029     integer i__, j, k, m1, m2, n1, n2, info;
00030     logical normaltransr, lside;
00031     extern logical lsame_(char *, char *);
00032     extern /* Subroutine */ int sgemm_(char *, char *, integer *, integer *, 
00033             integer *, real *, real *, integer *, real *, integer *, real *, 
00034             real *, integer *);
00035     logical lower;
00036     extern /* Subroutine */ int strsm_(char *, char *, char *, char *, 
00037             integer *, integer *, real *, real *, integer *, real *, integer *
00038 ), xerbla_(char *, integer *);
00039     logical misodd, nisodd, notrans;
00040 
00041 
00042 /*  -- LAPACK routine (version 3.2.1)                                    -- */
00043 
00044 /*  -- Contributed by Fred Gustavson of the IBM Watson Research Center -- */
00045 /*  -- April 2009                                                      -- */
00046 
00047 /*  -- LAPACK is a software package provided by Univ. of Tennessee,    -- */
00048 /*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
00049 
00050 /*     .. */
00051 /*     .. Scalar Arguments .. */
00052 /*     .. */
00053 /*     .. Array Arguments .. */
00054 /*     .. */
00055 
00056 /*  Purpose */
00057 /*  ======= */
00058 
00059 /*  Level 3 BLAS like routine for A in RFP Format. */
00060 
00061 /*  STFSM  solves the matrix equation */
00062 
00063 /*     op( A )*X = alpha*B  or  X*op( A ) = alpha*B */
00064 
00065 /*  where alpha is a scalar, X and B are m by n matrices, A is a unit, or */
00066 /*  non-unit,  upper or lower triangular matrix  and  op( A )  is one  of */
00067 
00068 /*     op( A ) = A   or   op( A ) = A'. */
00069 
00070 /*  A is in Rectangular Full Packed (RFP) Format. */
00071 
00072 /*  The matrix X is overwritten on B. */
00073 
00074 /*  Arguments */
00075 /*  ========== */
00076 
00077 /*  TRANSR - (input) CHARACTER */
00078 /*          = 'N':  The Normal Form of RFP A is stored; */
00079 /*          = 'T':  The Transpose Form of RFP A is stored. */
00080 
00081 /*  SIDE   - (input) CHARACTER */
00082 /*           On entry, SIDE specifies whether op( A ) appears on the left */
00083 /*           or right of X as follows: */
00084 
00085 /*              SIDE = 'L' or 'l'   op( A )*X = alpha*B. */
00086 
00087 /*              SIDE = 'R' or 'r'   X*op( A ) = alpha*B. */
00088 
00089 /*           Unchanged on exit. */
00090 
00091 /*  UPLO   - (input) CHARACTER */
00092 /*           On entry, UPLO specifies whether the RFP matrix A came from */
00093 /*           an upper or lower triangular matrix as follows: */
00094 /*           UPLO = 'U' or 'u' RFP A came from an upper triangular matrix */
00095 /*           UPLO = 'L' or 'l' RFP A came from a  lower triangular matrix */
00096 
00097 /*           Unchanged on exit. */
00098 
00099 /*  TRANS  - (input) CHARACTER */
00100 /*           On entry, TRANS  specifies the form of op( A ) to be used */
00101 /*           in the matrix multiplication as follows: */
00102 
00103 /*              TRANS  = 'N' or 'n'   op( A ) = A. */
00104 
00105 /*              TRANS  = 'T' or 't'   op( A ) = A'. */
00106 
00107 /*           Unchanged on exit. */
00108 
00109 /*  DIAG   - (input) CHARACTER */
00110 /*           On entry, DIAG specifies whether or not RFP A is unit */
00111 /*           triangular as follows: */
00112 
00113 /*              DIAG = 'U' or 'u'   A is assumed to be unit triangular. */
00114 
00115 /*              DIAG = 'N' or 'n'   A is not assumed to be unit */
00116 /*                                  triangular. */
00117 
00118 /*           Unchanged on exit. */
00119 
00120 /*  M      - (input) INTEGER. */
00121 /*           On entry, M specifies the number of rows of B. M must be at */
00122 /*           least zero. */
00123 /*           Unchanged on exit. */
00124 
00125 /*  N      - (input) INTEGER. */
00126 /*           On entry, N specifies the number of columns of B.  N must be */
00127 /*           at least zero. */
00128 /*           Unchanged on exit. */
00129 
00130 /*  ALPHA  - (input) REAL. */
00131 /*           On entry,  ALPHA specifies the scalar  alpha. When  alpha is */
00132 /*           zero then  A is not referenced and  B need not be set before */
00133 /*           entry. */
00134 /*           Unchanged on exit. */
00135 
00136 /*  A      - (input) REAL array, dimension (NT); */
00137 /*           NT = N*(N+1)/2. On entry, the matrix A in RFP Format. */
00138 /*           RFP Format is described by TRANSR, UPLO and N as follows: */
00139 /*           If TRANSR='N' then RFP A is (0:N,0:K-1) when N is even; */
00140 /*           K=N/2. RFP A is (0:N-1,0:K) when N is odd; K=N/2. If */
00141 /*           TRANSR = 'T' then RFP is the transpose of RFP A as */
00142 /*           defined when TRANSR = 'N'. The contents of RFP A are defined */
00143 /*           by UPLO as follows: If UPLO = 'U' the RFP A contains the NT */
00144 /*           elements of upper packed A either in normal or */
00145 /*           transpose Format. If UPLO = 'L' the RFP A contains */
00146 /*           the NT elements of lower packed A either in normal or */
00147 /*           transpose Format. The LDA of RFP A is (N+1)/2 when */
00148 /*           TRANSR = 'T'. When TRANSR is 'N' the LDA is N+1 when N is */
00149 /*           even and is N when is odd. */
00150 /*           See the Note below for more details. Unchanged on exit. */
00151 
00152 /*  B      - (input/ouptut) REAL array,  DIMENSION (LDB,N) */
00153 /*           Before entry,  the leading  m by n part of the array  B must */
00154 /*           contain  the  right-hand  side  matrix  B,  and  on exit  is */
00155 /*           overwritten by the solution matrix  X. */
00156 
00157 /*  LDB    - (input) INTEGER. */
00158 /*           On entry, LDB specifies the first dimension of B as declared */
00159 /*           in  the  calling  (sub)  program.   LDB  must  be  at  least */
00160 /*           max( 1, m ). */
00161 /*           Unchanged on exit. */
00162 
00163 /*  Further Details */
00164 /*  =============== */
00165 
00166 /*  We first consider Rectangular Full Packed (RFP) Format when N is */
00167 /*  even. We give an example where N = 6. */
00168 
00169 /*      AP is Upper             AP is Lower */
00170 
00171 /*   00 01 02 03 04 05       00 */
00172 /*      11 12 13 14 15       10 11 */
00173 /*         22 23 24 25       20 21 22 */
00174 /*            33 34 35       30 31 32 33 */
00175 /*               44 45       40 41 42 43 44 */
00176 /*                  55       50 51 52 53 54 55 */
00177 
00178 
00179 /*  Let TRANSR = 'N'. RFP holds AP as follows: */
00180 /*  For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last */
00181 /*  three columns of AP upper. The lower triangle A(4:6,0:2) consists of */
00182 /*  the transpose of the first three columns of AP upper. */
00183 /*  For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first */
00184 /*  three columns of AP lower. The upper triangle A(0:2,0:2) consists of */
00185 /*  the transpose of the last three columns of AP lower. */
00186 /*  This covers the case N even and TRANSR = 'N'. */
00187 
00188 /*         RFP A                   RFP A */
00189 
00190 /*        03 04 05                33 43 53 */
00191 /*        13 14 15                00 44 54 */
00192 /*        23 24 25                10 11 55 */
00193 /*        33 34 35                20 21 22 */
00194 /*        00 44 45                30 31 32 */
00195 /*        01 11 55                40 41 42 */
00196 /*        02 12 22                50 51 52 */
00197 
00198 /*  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the */
00199 /*  transpose of RFP A above. One therefore gets: */
00200 
00201 
00202 /*           RFP A                   RFP A */
00203 
00204 /*     03 13 23 33 00 01 02    33 00 10 20 30 40 50 */
00205 /*     04 14 24 34 44 11 12    43 44 11 21 31 41 51 */
00206 /*     05 15 25 35 45 55 22    53 54 55 22 32 42 52 */
00207 
00208 
00209 /*  We first consider Rectangular Full Packed (RFP) Format when N is */
00210 /*  odd. We give an example where N = 5. */
00211 
00212 /*     AP is Upper                 AP is Lower */
00213 
00214 /*   00 01 02 03 04              00 */
00215 /*      11 12 13 14              10 11 */
00216 /*         22 23 24              20 21 22 */
00217 /*            33 34              30 31 32 33 */
00218 /*               44              40 41 42 43 44 */
00219 
00220 
00221 /*  Let TRANSR = 'N'. RFP holds AP as follows: */
00222 /*  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last */
00223 /*  three columns of AP upper. The lower triangle A(3:4,0:1) consists of */
00224 /*  the transpose of the first two columns of AP upper. */
00225 /*  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first */
00226 /*  three columns of AP lower. The upper triangle A(0:1,1:2) consists of */
00227 /*  the transpose of the last two columns of AP lower. */
00228 /*  This covers the case N odd and TRANSR = 'N'. */
00229 
00230 /*         RFP A                   RFP A */
00231 
00232 /*        02 03 04                00 33 43 */
00233 /*        12 13 14                10 11 44 */
00234 /*        22 23 24                20 21 22 */
00235 /*        00 33 34                30 31 32 */
00236 /*        01 11 44                40 41 42 */
00237 
00238 /*  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the */
00239 /*  transpose of RFP A above. One therefore gets: */
00240 
00241 /*           RFP A                   RFP A */
00242 
00243 /*     02 12 22 00 01             00 10 20 30 40 50 */
00244 /*     03 13 23 33 11             33 11 21 31 41 51 */
00245 /*     04 14 24 34 44             43 44 22 32 42 52 */
00246 
00247 /*  Reference */
00248 /*  ========= */
00249 
00250 /*  ===================================================================== */
00251 
00252 /*     .. */
00253 /*     .. Parameters .. */
00254 /*     .. */
00255 /*     .. Local Scalars .. */
00256 /*     .. */
00257 /*     .. External Functions .. */
00258 /*     .. */
00259 /*     .. External Subroutines .. */
00260 /*     .. */
00261 /*     .. Intrinsic Functions .. */
00262 /*     .. */
00263 /*     .. Executable Statements .. */
00264 
00265 /*     Test the input parameters. */
00266 
00267     /* Parameter adjustments */
00268     b_dim1 = *ldb - 1 - 0 + 1;
00269     b_offset = 0 + b_dim1 * 0;
00270     b -= b_offset;
00271 
00272     /* Function Body */
00273     info = 0;
00274     normaltransr = lsame_(transr, "N");
00275     lside = lsame_(side, "L");
00276     lower = lsame_(uplo, "L");
00277     notrans = lsame_(trans, "N");
00278     if (! normaltransr && ! lsame_(transr, "T")) {
00279         info = -1;
00280     } else if (! lside && ! lsame_(side, "R")) {
00281         info = -2;
00282     } else if (! lower && ! lsame_(uplo, "U")) {
00283         info = -3;
00284     } else if (! notrans && ! lsame_(trans, "T")) {
00285         info = -4;
00286     } else if (! lsame_(diag, "N") && ! lsame_(diag, 
00287             "U")) {
00288         info = -5;
00289     } else if (*m < 0) {
00290         info = -6;
00291     } else if (*n < 0) {
00292         info = -7;
00293     } else if (*ldb < max(1,*m)) {
00294         info = -11;
00295     }
00296     if (info != 0) {
00297         i__1 = -info;
00298         xerbla_("STFSM ", &i__1);
00299         return 0;
00300     }
00301 
00302 /*     Quick return when ( (N.EQ.0).OR.(M.EQ.0) ) */
00303 
00304     if (*m == 0 || *n == 0) {
00305         return 0;
00306     }
00307 
00308 /*     Quick return when ALPHA.EQ.(0D+0) */
00309 
00310     if (*alpha == 0.f) {
00311         i__1 = *n - 1;
00312         for (j = 0; j <= i__1; ++j) {
00313             i__2 = *m - 1;
00314             for (i__ = 0; i__ <= i__2; ++i__) {
00315                 b[i__ + j * b_dim1] = 0.f;
00316 /* L10: */
00317             }
00318 /* L20: */
00319         }
00320         return 0;
00321     }
00322 
00323     if (lside) {
00324 
00325 /*        SIDE = 'L' */
00326 
00327 /*        A is M-by-M. */
00328 /*        If M is odd, set NISODD = .TRUE., and M1 and M2. */
00329 /*        If M is even, NISODD = .FALSE., and M. */
00330 
00331         if (*m % 2 == 0) {
00332             misodd = FALSE_;
00333             k = *m / 2;
00334         } else {
00335             misodd = TRUE_;
00336             if (lower) {
00337                 m2 = *m / 2;
00338                 m1 = *m - m2;
00339             } else {
00340                 m1 = *m / 2;
00341                 m2 = *m - m1;
00342             }
00343         }
00344 
00345         if (misodd) {
00346 
00347 /*           SIDE = 'L' and N is odd */
00348 
00349             if (normaltransr) {
00350 
00351 /*              SIDE = 'L', N is odd, and TRANSR = 'N' */
00352 
00353                 if (lower) {
00354 
00355 /*                 SIDE  ='L', N is odd, TRANSR = 'N', and UPLO = 'L' */
00356 
00357                     if (notrans) {
00358 
00359 /*                    SIDE  ='L', N is odd, TRANSR = 'N', UPLO = 'L', and */
00360 /*                    TRANS = 'N' */
00361 
00362                         if (*m == 1) {
00363                             strsm_("L", "L", "N", diag, &m1, n, alpha, a, m, &
00364                                     b[b_offset], ldb);
00365                         } else {
00366                             strsm_("L", "L", "N", diag, &m1, n, alpha, a, m, &
00367                                     b[b_offset], ldb);
00368                             sgemm_("N", "N", &m2, n, &m1, &c_b23, &a[m1], m, &
00369                                     b[b_offset], ldb, alpha, &b[m1], ldb);
00370                             strsm_("L", "U", "T", diag, &m2, n, &c_b27, &a[*m]
00371 , m, &b[m1], ldb);
00372                         }
00373 
00374                     } else {
00375 
00376 /*                    SIDE  ='L', N is odd, TRANSR = 'N', UPLO = 'L', and */
00377 /*                    TRANS = 'T' */
00378 
00379                         if (*m == 1) {
00380                             strsm_("L", "L", "T", diag, &m1, n, alpha, a, m, &
00381                                     b[b_offset], ldb);
00382                         } else {
00383                             strsm_("L", "U", "N", diag, &m2, n, alpha, &a[*m], 
00384                                      m, &b[m1], ldb);
00385                             sgemm_("T", "N", &m1, n, &m2, &c_b23, &a[m1], m, &
00386                                     b[m1], ldb, alpha, &b[b_offset], ldb);
00387                             strsm_("L", "L", "T", diag, &m1, n, &c_b27, a, m, 
00388                                     &b[b_offset], ldb);
00389                         }
00390 
00391                     }
00392 
00393                 } else {
00394 
00395 /*                 SIDE  ='L', N is odd, TRANSR = 'N', and UPLO = 'U' */
00396 
00397                     if (! notrans) {
00398 
00399 /*                    SIDE  ='L', N is odd, TRANSR = 'N', UPLO = 'U', and */
00400 /*                    TRANS = 'N' */
00401 
00402                         strsm_("L", "L", "N", diag, &m1, n, alpha, &a[m2], m, 
00403                                 &b[b_offset], ldb);
00404                         sgemm_("T", "N", &m2, n, &m1, &c_b23, a, m, &b[
00405                                 b_offset], ldb, alpha, &b[m1], ldb);
00406                         strsm_("L", "U", "T", diag, &m2, n, &c_b27, &a[m1], m, 
00407                                  &b[m1], ldb);
00408 
00409                     } else {
00410 
00411 /*                    SIDE  ='L', N is odd, TRANSR = 'N', UPLO = 'U', and */
00412 /*                    TRANS = 'T' */
00413 
00414                         strsm_("L", "U", "N", diag, &m2, n, alpha, &a[m1], m, 
00415                                 &b[m1], ldb);
00416                         sgemm_("N", "N", &m1, n, &m2, &c_b23, a, m, &b[m1], 
00417                                 ldb, alpha, &b[b_offset], ldb);
00418                         strsm_("L", "L", "T", diag, &m1, n, &c_b27, &a[m2], m, 
00419                                  &b[b_offset], ldb);
00420 
00421                     }
00422 
00423                 }
00424 
00425             } else {
00426 
00427 /*              SIDE = 'L', N is odd, and TRANSR = 'T' */
00428 
00429                 if (lower) {
00430 
00431 /*                 SIDE  ='L', N is odd, TRANSR = 'T', and UPLO = 'L' */
00432 
00433                     if (notrans) {
00434 
00435 /*                    SIDE  ='L', N is odd, TRANSR = 'T', UPLO = 'L', and */
00436 /*                    TRANS = 'N' */
00437 
00438                         if (*m == 1) {
00439                             strsm_("L", "U", "T", diag, &m1, n, alpha, a, &m1, 
00440                                      &b[b_offset], ldb);
00441                         } else {
00442                             strsm_("L", "U", "T", diag, &m1, n, alpha, a, &m1, 
00443                                      &b[b_offset], ldb);
00444                             sgemm_("T", "N", &m2, n, &m1, &c_b23, &a[m1 * m1], 
00445                                      &m1, &b[b_offset], ldb, alpha, &b[m1], 
00446                                     ldb);
00447                             strsm_("L", "L", "N", diag, &m2, n, &c_b27, &a[1], 
00448                                      &m1, &b[m1], ldb);
00449                         }
00450 
00451                     } else {
00452 
00453 /*                    SIDE  ='L', N is odd, TRANSR = 'T', UPLO = 'L', and */
00454 /*                    TRANS = 'T' */
00455 
00456                         if (*m == 1) {
00457                             strsm_("L", "U", "N", diag, &m1, n, alpha, a, &m1, 
00458                                      &b[b_offset], ldb);
00459                         } else {
00460                             strsm_("L", "L", "T", diag, &m2, n, alpha, &a[1], 
00461                                     &m1, &b[m1], ldb);
00462                             sgemm_("N", "N", &m1, n, &m2, &c_b23, &a[m1 * m1], 
00463                                      &m1, &b[m1], ldb, alpha, &b[b_offset], 
00464                                     ldb);
00465                             strsm_("L", "U", "N", diag, &m1, n, &c_b27, a, &
00466                                     m1, &b[b_offset], ldb);
00467                         }
00468 
00469                     }
00470 
00471                 } else {
00472 
00473 /*                 SIDE  ='L', N is odd, TRANSR = 'T', and UPLO = 'U' */
00474 
00475                     if (! notrans) {
00476 
00477 /*                    SIDE  ='L', N is odd, TRANSR = 'T', UPLO = 'U', and */
00478 /*                    TRANS = 'N' */
00479 
00480                         strsm_("L", "U", "T", diag, &m1, n, alpha, &a[m2 * m2]
00481 , &m2, &b[b_offset], ldb);
00482                         sgemm_("N", "N", &m2, n, &m1, &c_b23, a, &m2, &b[
00483                                 b_offset], ldb, alpha, &b[m1], ldb);
00484                         strsm_("L", "L", "N", diag, &m2, n, &c_b27, &a[m1 * 
00485                                 m2], &m2, &b[m1], ldb);
00486 
00487                     } else {
00488 
00489 /*                    SIDE  ='L', N is odd, TRANSR = 'T', UPLO = 'U', and */
00490 /*                    TRANS = 'T' */
00491 
00492                         strsm_("L", "L", "T", diag, &m2, n, alpha, &a[m1 * m2]
00493 , &m2, &b[m1], ldb);
00494                         sgemm_("T", "N", &m1, n, &m2, &c_b23, a, &m2, &b[m1], 
00495                                 ldb, alpha, &b[b_offset], ldb);
00496                         strsm_("L", "U", "N", diag, &m1, n, &c_b27, &a[m2 * 
00497                                 m2], &m2, &b[b_offset], ldb);
00498 
00499                     }
00500 
00501                 }
00502 
00503             }
00504 
00505         } else {
00506 
00507 /*           SIDE = 'L' and N is even */
00508 
00509             if (normaltransr) {
00510 
00511 /*              SIDE = 'L', N is even, and TRANSR = 'N' */
00512 
00513                 if (lower) {
00514 
00515 /*                 SIDE  ='L', N is even, TRANSR = 'N', and UPLO = 'L' */
00516 
00517                     if (notrans) {
00518 
00519 /*                    SIDE  ='L', N is even, TRANSR = 'N', UPLO = 'L', */
00520 /*                    and TRANS = 'N' */
00521 
00522                         i__1 = *m + 1;
00523                         strsm_("L", "L", "N", diag, &k, n, alpha, &a[1], &
00524                                 i__1, &b[b_offset], ldb);
00525                         i__1 = *m + 1;
00526                         sgemm_("N", "N", &k, n, &k, &c_b23, &a[k + 1], &i__1, 
00527                                 &b[b_offset], ldb, alpha, &b[k], ldb);
00528                         i__1 = *m + 1;
00529                         strsm_("L", "U", "T", diag, &k, n, &c_b27, a, &i__1, &
00530                                 b[k], ldb);
00531 
00532                     } else {
00533 
00534 /*                    SIDE  ='L', N is even, TRANSR = 'N', UPLO = 'L', */
00535 /*                    and TRANS = 'T' */
00536 
00537                         i__1 = *m + 1;
00538                         strsm_("L", "U", "N", diag, &k, n, alpha, a, &i__1, &
00539                                 b[k], ldb);
00540                         i__1 = *m + 1;
00541                         sgemm_("T", "N", &k, n, &k, &c_b23, &a[k + 1], &i__1, 
00542                                 &b[k], ldb, alpha, &b[b_offset], ldb);
00543                         i__1 = *m + 1;
00544                         strsm_("L", "L", "T", diag, &k, n, &c_b27, &a[1], &
00545                                 i__1, &b[b_offset], ldb);
00546 
00547                     }
00548 
00549                 } else {
00550 
00551 /*                 SIDE  ='L', N is even, TRANSR = 'N', and UPLO = 'U' */
00552 
00553                     if (! notrans) {
00554 
00555 /*                    SIDE  ='L', N is even, TRANSR = 'N', UPLO = 'U', */
00556 /*                    and TRANS = 'N' */
00557 
00558                         i__1 = *m + 1;
00559                         strsm_("L", "L", "N", diag, &k, n, alpha, &a[k + 1], &
00560                                 i__1, &b[b_offset], ldb);
00561                         i__1 = *m + 1;
00562                         sgemm_("T", "N", &k, n, &k, &c_b23, a, &i__1, &b[
00563                                 b_offset], ldb, alpha, &b[k], ldb);
00564                         i__1 = *m + 1;
00565                         strsm_("L", "U", "T", diag, &k, n, &c_b27, &a[k], &
00566                                 i__1, &b[k], ldb);
00567 
00568                     } else {
00569 
00570 /*                    SIDE  ='L', N is even, TRANSR = 'N', UPLO = 'U', */
00571 /*                    and TRANS = 'T' */
00572                         i__1 = *m + 1;
00573                         strsm_("L", "U", "N", diag, &k, n, alpha, &a[k], &
00574                                 i__1, &b[k], ldb);
00575                         i__1 = *m + 1;
00576                         sgemm_("N", "N", &k, n, &k, &c_b23, a, &i__1, &b[k], 
00577                                 ldb, alpha, &b[b_offset], ldb);
00578                         i__1 = *m + 1;
00579                         strsm_("L", "L", "T", diag, &k, n, &c_b27, &a[k + 1], 
00580                                 &i__1, &b[b_offset], ldb);
00581 
00582                     }
00583 
00584                 }
00585 
00586             } else {
00587 
00588 /*              SIDE = 'L', N is even, and TRANSR = 'T' */
00589 
00590                 if (lower) {
00591 
00592 /*                 SIDE  ='L', N is even, TRANSR = 'T', and UPLO = 'L' */
00593 
00594                     if (notrans) {
00595 
00596 /*                    SIDE  ='L', N is even, TRANSR = 'T', UPLO = 'L', */
00597 /*                    and TRANS = 'N' */
00598 
00599                         strsm_("L", "U", "T", diag, &k, n, alpha, &a[k], &k, &
00600                                 b[b_offset], ldb);
00601                         sgemm_("T", "N", &k, n, &k, &c_b23, &a[k * (k + 1)], &
00602                                 k, &b[b_offset], ldb, alpha, &b[k], ldb);
00603                         strsm_("L", "L", "N", diag, &k, n, &c_b27, a, &k, &b[
00604                                 k], ldb);
00605 
00606                     } else {
00607 
00608 /*                    SIDE  ='L', N is even, TRANSR = 'T', UPLO = 'L', */
00609 /*                    and TRANS = 'T' */
00610 
00611                         strsm_("L", "L", "T", diag, &k, n, alpha, a, &k, &b[k]
00612 , ldb);
00613                         sgemm_("N", "N", &k, n, &k, &c_b23, &a[k * (k + 1)], &
00614                                 k, &b[k], ldb, alpha, &b[b_offset], ldb);
00615                         strsm_("L", "U", "N", diag, &k, n, &c_b27, &a[k], &k, 
00616                                 &b[b_offset], ldb);
00617 
00618                     }
00619 
00620                 } else {
00621 
00622 /*                 SIDE  ='L', N is even, TRANSR = 'T', and UPLO = 'U' */
00623 
00624                     if (! notrans) {
00625 
00626 /*                    SIDE  ='L', N is even, TRANSR = 'T', UPLO = 'U', */
00627 /*                    and TRANS = 'N' */
00628 
00629                         strsm_("L", "U", "T", diag, &k, n, alpha, &a[k * (k + 
00630                                 1)], &k, &b[b_offset], ldb);
00631                         sgemm_("N", "N", &k, n, &k, &c_b23, a, &k, &b[
00632                                 b_offset], ldb, alpha, &b[k], ldb);
00633                         strsm_("L", "L", "N", diag, &k, n, &c_b27, &a[k * k], 
00634                                 &k, &b[k], ldb);
00635 
00636                     } else {
00637 
00638 /*                    SIDE  ='L', N is even, TRANSR = 'T', UPLO = 'U', */
00639 /*                    and TRANS = 'T' */
00640 
00641                         strsm_("L", "L", "T", diag, &k, n, alpha, &a[k * k], &
00642                                 k, &b[k], ldb);
00643                         sgemm_("T", "N", &k, n, &k, &c_b23, a, &k, &b[k], ldb, 
00644                                  alpha, &b[b_offset], ldb);
00645                         strsm_("L", "U", "N", diag, &k, n, &c_b27, &a[k * (k 
00646                                 + 1)], &k, &b[b_offset], ldb);
00647 
00648                     }
00649 
00650                 }
00651 
00652             }
00653 
00654         }
00655 
00656     } else {
00657 
00658 /*        SIDE = 'R' */
00659 
00660 /*        A is N-by-N. */
00661 /*        If N is odd, set NISODD = .TRUE., and N1 and N2. */
00662 /*        If N is even, NISODD = .FALSE., and K. */
00663 
00664         if (*n % 2 == 0) {
00665             nisodd = FALSE_;
00666             k = *n / 2;
00667         } else {
00668             nisodd = TRUE_;
00669             if (lower) {
00670                 n2 = *n / 2;
00671                 n1 = *n - n2;
00672             } else {
00673                 n1 = *n / 2;
00674                 n2 = *n - n1;
00675             }
00676         }
00677 
00678         if (nisodd) {
00679 
00680 /*           SIDE = 'R' and N is odd */
00681 
00682             if (normaltransr) {
00683 
00684 /*              SIDE = 'R', N is odd, and TRANSR = 'N' */
00685 
00686                 if (lower) {
00687 
00688 /*                 SIDE  ='R', N is odd, TRANSR = 'N', and UPLO = 'L' */
00689 
00690                     if (notrans) {
00691 
00692 /*                    SIDE  ='R', N is odd, TRANSR = 'N', UPLO = 'L', and */
00693 /*                    TRANS = 'N' */
00694 
00695                         strsm_("R", "U", "T", diag, m, &n2, alpha, &a[*n], n, 
00696                                 &b[n1 * b_dim1], ldb);
00697                         sgemm_("N", "N", m, &n1, &n2, &c_b23, &b[n1 * b_dim1], 
00698                                  ldb, &a[n1], n, alpha, b, ldb);
00699                         strsm_("R", "L", "N", diag, m, &n1, &c_b27, a, n, b, 
00700                                 ldb);
00701 
00702                     } else {
00703 
00704 /*                    SIDE  ='R', N is odd, TRANSR = 'N', UPLO = 'L', and */
00705 /*                    TRANS = 'T' */
00706 
00707                         strsm_("R", "L", "T", diag, m, &n1, alpha, a, n, b, 
00708                                 ldb);
00709                         sgemm_("N", "T", m, &n2, &n1, &c_b23, b, ldb, &a[n1], 
00710                                 n, alpha, &b[n1 * b_dim1], ldb);
00711                         strsm_("R", "U", "N", diag, m, &n2, &c_b27, &a[*n], n, 
00712                                  &b[n1 * b_dim1], ldb);
00713 
00714                     }
00715 
00716                 } else {
00717 
00718 /*                 SIDE  ='R', N is odd, TRANSR = 'N', and UPLO = 'U' */
00719 
00720                     if (notrans) {
00721 
00722 /*                    SIDE  ='R', N is odd, TRANSR = 'N', UPLO = 'U', and */
00723 /*                    TRANS = 'N' */
00724 
00725                         strsm_("R", "L", "T", diag, m, &n1, alpha, &a[n2], n, 
00726                                 b, ldb);
00727                         sgemm_("N", "N", m, &n2, &n1, &c_b23, b, ldb, a, n, 
00728                                 alpha, &b[n1 * b_dim1], ldb);
00729                         strsm_("R", "U", "N", diag, m, &n2, &c_b27, &a[n1], n, 
00730                                  &b[n1 * b_dim1], ldb);
00731 
00732                     } else {
00733 
00734 /*                    SIDE  ='R', N is odd, TRANSR = 'N', UPLO = 'U', and */
00735 /*                    TRANS = 'T' */
00736 
00737                         strsm_("R", "U", "T", diag, m, &n2, alpha, &a[n1], n, 
00738                                 &b[n1 * b_dim1], ldb);
00739                         sgemm_("N", "T", m, &n1, &n2, &c_b23, &b[n1 * b_dim1], 
00740                                  ldb, a, n, alpha, b, ldb);
00741                         strsm_("R", "L", "N", diag, m, &n1, &c_b27, &a[n2], n, 
00742                                  b, ldb);
00743 
00744                     }
00745 
00746                 }
00747 
00748             } else {
00749 
00750 /*              SIDE = 'R', N is odd, and TRANSR = 'T' */
00751 
00752                 if (lower) {
00753 
00754 /*                 SIDE  ='R', N is odd, TRANSR = 'T', and UPLO = 'L' */
00755 
00756                     if (notrans) {
00757 
00758 /*                    SIDE  ='R', N is odd, TRANSR = 'T', UPLO = 'L', and */
00759 /*                    TRANS = 'N' */
00760 
00761                         strsm_("R", "L", "N", diag, m, &n2, alpha, &a[1], &n1, 
00762                                  &b[n1 * b_dim1], ldb);
00763                         sgemm_("N", "T", m, &n1, &n2, &c_b23, &b[n1 * b_dim1], 
00764                                  ldb, &a[n1 * n1], &n1, alpha, b, ldb);
00765                         strsm_("R", "U", "T", diag, m, &n1, &c_b27, a, &n1, b, 
00766                                  ldb);
00767 
00768                     } else {
00769 
00770 /*                    SIDE  ='R', N is odd, TRANSR = 'T', UPLO = 'L', and */
00771 /*                    TRANS = 'T' */
00772 
00773                         strsm_("R", "U", "N", diag, m, &n1, alpha, a, &n1, b, 
00774                                 ldb);
00775                         sgemm_("N", "N", m, &n2, &n1, &c_b23, b, ldb, &a[n1 * 
00776                                 n1], &n1, alpha, &b[n1 * b_dim1], ldb);
00777                         strsm_("R", "L", "T", diag, m, &n2, &c_b27, &a[1], &
00778                                 n1, &b[n1 * b_dim1], ldb);
00779 
00780                     }
00781 
00782                 } else {
00783 
00784 /*                 SIDE  ='R', N is odd, TRANSR = 'T', and UPLO = 'U' */
00785 
00786                     if (notrans) {
00787 
00788 /*                    SIDE  ='R', N is odd, TRANSR = 'T', UPLO = 'U', and */
00789 /*                    TRANS = 'N' */
00790 
00791                         strsm_("R", "U", "N", diag, m, &n1, alpha, &a[n2 * n2]
00792 , &n2, b, ldb);
00793                         sgemm_("N", "T", m, &n2, &n1, &c_b23, b, ldb, a, &n2, 
00794                                 alpha, &b[n1 * b_dim1], ldb);
00795                         strsm_("R", "L", "T", diag, m, &n2, &c_b27, &a[n1 * 
00796                                 n2], &n2, &b[n1 * b_dim1], ldb);
00797 
00798                     } else {
00799 
00800 /*                    SIDE  ='R', N is odd, TRANSR = 'T', UPLO = 'U', and */
00801 /*                    TRANS = 'T' */
00802 
00803                         strsm_("R", "L", "N", diag, m, &n2, alpha, &a[n1 * n2]
00804 , &n2, &b[n1 * b_dim1], ldb);
00805                         sgemm_("N", "N", m, &n1, &n2, &c_b23, &b[n1 * b_dim1], 
00806                                  ldb, a, &n2, alpha, b, ldb);
00807                         strsm_("R", "U", "T", diag, m, &n1, &c_b27, &a[n2 * 
00808                                 n2], &n2, b, ldb);
00809 
00810                     }
00811 
00812                 }
00813 
00814             }
00815 
00816         } else {
00817 
00818 /*           SIDE = 'R' and N is even */
00819 
00820             if (normaltransr) {
00821 
00822 /*              SIDE = 'R', N is even, and TRANSR = 'N' */
00823 
00824                 if (lower) {
00825 
00826 /*                 SIDE  ='R', N is even, TRANSR = 'N', and UPLO = 'L' */
00827 
00828                     if (notrans) {
00829 
00830 /*                    SIDE  ='R', N is even, TRANSR = 'N', UPLO = 'L', */
00831 /*                    and TRANS = 'N' */
00832 
00833                         i__1 = *n + 1;
00834                         strsm_("R", "U", "T", diag, m, &k, alpha, a, &i__1, &
00835                                 b[k * b_dim1], ldb);
00836                         i__1 = *n + 1;
00837                         sgemm_("N", "N", m, &k, &k, &c_b23, &b[k * b_dim1], 
00838                                 ldb, &a[k + 1], &i__1, alpha, b, ldb);
00839                         i__1 = *n + 1;
00840                         strsm_("R", "L", "N", diag, m, &k, &c_b27, &a[1], &
00841                                 i__1, b, ldb);
00842 
00843                     } else {
00844 
00845 /*                    SIDE  ='R', N is even, TRANSR = 'N', UPLO = 'L', */
00846 /*                    and TRANS = 'T' */
00847 
00848                         i__1 = *n + 1;
00849                         strsm_("R", "L", "T", diag, m, &k, alpha, &a[1], &
00850                                 i__1, b, ldb);
00851                         i__1 = *n + 1;
00852                         sgemm_("N", "T", m, &k, &k, &c_b23, b, ldb, &a[k + 1], 
00853                                  &i__1, alpha, &b[k * b_dim1], ldb);
00854                         i__1 = *n + 1;
00855                         strsm_("R", "U", "N", diag, m, &k, &c_b27, a, &i__1, &
00856                                 b[k * b_dim1], ldb);
00857 
00858                     }
00859 
00860                 } else {
00861 
00862 /*                 SIDE  ='R', N is even, TRANSR = 'N', and UPLO = 'U' */
00863 
00864                     if (notrans) {
00865 
00866 /*                    SIDE  ='R', N is even, TRANSR = 'N', UPLO = 'U', */
00867 /*                    and TRANS = 'N' */
00868 
00869                         i__1 = *n + 1;
00870                         strsm_("R", "L", "T", diag, m, &k, alpha, &a[k + 1], &
00871                                 i__1, b, ldb);
00872                         i__1 = *n + 1;
00873                         sgemm_("N", "N", m, &k, &k, &c_b23, b, ldb, a, &i__1, 
00874                                 alpha, &b[k * b_dim1], ldb);
00875                         i__1 = *n + 1;
00876                         strsm_("R", "U", "N", diag, m, &k, &c_b27, &a[k], &
00877                                 i__1, &b[k * b_dim1], ldb);
00878 
00879                     } else {
00880 
00881 /*                    SIDE  ='R', N is even, TRANSR = 'N', UPLO = 'U', */
00882 /*                    and TRANS = 'T' */
00883 
00884                         i__1 = *n + 1;
00885                         strsm_("R", "U", "T", diag, m, &k, alpha, &a[k], &
00886                                 i__1, &b[k * b_dim1], ldb);
00887                         i__1 = *n + 1;
00888                         sgemm_("N", "T", m, &k, &k, &c_b23, &b[k * b_dim1], 
00889                                 ldb, a, &i__1, alpha, b, ldb);
00890                         i__1 = *n + 1;
00891                         strsm_("R", "L", "N", diag, m, &k, &c_b27, &a[k + 1], 
00892                                 &i__1, b, ldb);
00893 
00894                     }
00895 
00896                 }
00897 
00898             } else {
00899 
00900 /*              SIDE = 'R', N is even, and TRANSR = 'T' */
00901 
00902                 if (lower) {
00903 
00904 /*                 SIDE  ='R', N is even, TRANSR = 'T', and UPLO = 'L' */
00905 
00906                     if (notrans) {
00907 
00908 /*                    SIDE  ='R', N is even, TRANSR = 'T', UPLO = 'L', */
00909 /*                    and TRANS = 'N' */
00910 
00911                         strsm_("R", "L", "N", diag, m, &k, alpha, a, &k, &b[k 
00912                                 * b_dim1], ldb);
00913                         sgemm_("N", "T", m, &k, &k, &c_b23, &b[k * b_dim1], 
00914                                 ldb, &a[(k + 1) * k], &k, alpha, b, ldb);
00915                         strsm_("R", "U", "T", diag, m, &k, &c_b27, &a[k], &k, 
00916                                 b, ldb);
00917 
00918                     } else {
00919 
00920 /*                    SIDE  ='R', N is even, TRANSR = 'T', UPLO = 'L', */
00921 /*                    and TRANS = 'T' */
00922 
00923                         strsm_("R", "U", "N", diag, m, &k, alpha, &a[k], &k, 
00924                                 b, ldb);
00925                         sgemm_("N", "N", m, &k, &k, &c_b23, b, ldb, &a[(k + 1)
00926                                  * k], &k, alpha, &b[k * b_dim1], ldb);
00927                         strsm_("R", "L", "T", diag, m, &k, &c_b27, a, &k, &b[
00928                                 k * b_dim1], ldb);
00929 
00930                     }
00931 
00932                 } else {
00933 
00934 /*                 SIDE  ='R', N is even, TRANSR = 'T', and UPLO = 'U' */
00935 
00936                     if (notrans) {
00937 
00938 /*                    SIDE  ='R', N is even, TRANSR = 'T', UPLO = 'U', */
00939 /*                    and TRANS = 'N' */
00940 
00941                         strsm_("R", "U", "N", diag, m, &k, alpha, &a[(k + 1) *
00942                                  k], &k, b, ldb);
00943                         sgemm_("N", "T", m, &k, &k, &c_b23, b, ldb, a, &k, 
00944                                 alpha, &b[k * b_dim1], ldb);
00945                         strsm_("R", "L", "T", diag, m, &k, &c_b27, &a[k * k], 
00946                                 &k, &b[k * b_dim1], ldb);
00947 
00948                     } else {
00949 
00950 /*                    SIDE  ='R', N is even, TRANSR = 'T', UPLO = 'U', */
00951 /*                    and TRANS = 'T' */
00952 
00953                         strsm_("R", "L", "N", diag, m, &k, alpha, &a[k * k], &
00954                                 k, &b[k * b_dim1], ldb);
00955                         sgemm_("N", "N", m, &k, &k, &c_b23, &b[k * b_dim1], 
00956                                 ldb, a, &k, alpha, b, ldb);
00957                         strsm_("R", "U", "T", diag, m, &k, &c_b27, &a[(k + 1) 
00958                                 * k], &k, b, ldb);
00959 
00960                     }
00961 
00962                 }
00963 
00964             }
00965 
00966         }
00967     }
00968 
00969     return 0;
00970 
00971 /*     End of STFSM */
00972 
00973 } /* stfsm_ */

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