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

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