zlatm5.c
Go to the documentation of this file.
00001 /* zlatm5.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 doublecomplex c_b1 = {1.,0.};
00019 static doublecomplex c_b3 = {0.,0.};
00020 static doublecomplex c_b5 = {20.,0.};
00021 
00022 /* Subroutine */ int zlatm5_(integer *prtype, integer *m, integer *n, 
00023         doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb, 
00024         doublecomplex *c__, integer *ldc, doublecomplex *d__, integer *ldd, 
00025         doublecomplex *e, integer *lde, doublecomplex *f, integer *ldf, 
00026         doublecomplex *r__, integer *ldr, doublecomplex *l, integer *ldl, 
00027         doublereal *alpha, integer *qblcka, integer *qblckb)
00028 {
00029     /* System generated locals */
00030     integer a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, d_dim1, 
00031             d_offset, e_dim1, e_offset, f_dim1, f_offset, l_dim1, l_offset, 
00032             r_dim1, r_offset, i__1, i__2, i__3, i__4;
00033     doublereal d__1;
00034     doublecomplex z__1, z__2, z__3, z__4, z__5;
00035 
00036     /* Builtin functions */
00037     void z_sin(doublecomplex *, doublecomplex *), z_div(doublecomplex *, 
00038             doublecomplex *, doublecomplex *);
00039 
00040     /* Local variables */
00041     integer i__, j, k;
00042     doublecomplex imeps, reeps;
00043     extern /* Subroutine */ int zgemm_(char *, char *, integer *, integer *, 
00044             integer *, doublecomplex *, doublecomplex *, integer *, 
00045             doublecomplex *, integer *, doublecomplex *, doublecomplex *, 
00046             integer *);
00047 
00048 
00049 /*  -- LAPACK test routine (version 3.1) -- */
00050 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00051 /*     November 2006 */
00052 
00053 /*     .. Scalar Arguments .. */
00054 /*     .. */
00055 /*     .. Array Arguments .. */
00056 /*     .. */
00057 
00058 /*  Purpose */
00059 /*  ======= */
00060 
00061 /*  ZLATM5 generates matrices involved in the Generalized Sylvester */
00062 /*  equation: */
00063 
00064 /*      A * R - L * B = C */
00065 /*      D * R - L * E = F */
00066 
00067 /*  They also satisfy (the diagonalization condition) */
00068 
00069 /*   [ I -L ] ( [ A  -C ], [ D -F ] ) [ I  R ] = ( [ A    ], [ D    ] ) */
00070 /*   [    I ] ( [     B ]  [    E ] ) [    I ]   ( [    B ]  [    E ] ) */
00071 
00072 
00073 /*  Arguments */
00074 /*  ========= */
00075 
00076 /*  PRTYPE  (input) INTEGER */
00077 /*          "Points" to a certian type of the matrices to generate */
00078 /*          (see futher details). */
00079 
00080 /*  M       (input) INTEGER */
00081 /*          Specifies the order of A and D and the number of rows in */
00082 /*          C, F,  R and L. */
00083 
00084 /*  N       (input) INTEGER */
00085 /*          Specifies the order of B and E and the number of columns in */
00086 /*          C, F, R and L. */
00087 
00088 /*  A       (output) COMPLEX*16 array, dimension (LDA, M). */
00089 /*          On exit A M-by-M is initialized according to PRTYPE. */
00090 
00091 /*  LDA     (input) INTEGER */
00092 /*          The leading dimension of A. */
00093 
00094 /*  B       (output) COMPLEX*16 array, dimension (LDB, N). */
00095 /*          On exit B N-by-N is initialized according to PRTYPE. */
00096 
00097 /*  LDB     (input) INTEGER */
00098 /*          The leading dimension of B. */
00099 
00100 /*  C       (output) COMPLEX*16 array, dimension (LDC, N). */
00101 /*          On exit C M-by-N is initialized according to PRTYPE. */
00102 
00103 /*  LDC     (input) INTEGER */
00104 /*          The leading dimension of C. */
00105 
00106 /*  D       (output) COMPLEX*16 array, dimension (LDD, M). */
00107 /*          On exit D M-by-M is initialized according to PRTYPE. */
00108 
00109 /*  LDD     (input) INTEGER */
00110 /*          The leading dimension of D. */
00111 
00112 /*  E       (output) COMPLEX*16 array, dimension (LDE, N). */
00113 /*          On exit E N-by-N is initialized according to PRTYPE. */
00114 
00115 /*  LDE     (input) INTEGER */
00116 /*          The leading dimension of E. */
00117 
00118 /*  F       (output) COMPLEX*16 array, dimension (LDF, N). */
00119 /*          On exit F M-by-N is initialized according to PRTYPE. */
00120 
00121 /*  LDF     (input) INTEGER */
00122 /*          The leading dimension of F. */
00123 
00124 /*  R       (output) COMPLEX*16 array, dimension (LDR, N). */
00125 /*          On exit R M-by-N is initialized according to PRTYPE. */
00126 
00127 /*  LDR     (input) INTEGER */
00128 /*          The leading dimension of R. */
00129 
00130 /*  L       (output) COMPLEX*16 array, dimension (LDL, N). */
00131 /*          On exit L M-by-N is initialized according to PRTYPE. */
00132 
00133 /*  LDL     (input) INTEGER */
00134 /*          The leading dimension of L. */
00135 
00136 /*  ALPHA   (input) DOUBLE PRECISION */
00137 /*          Parameter used in generating PRTYPE = 1 and 5 matrices. */
00138 
00139 /*  QBLCKA  (input) INTEGER */
00140 /*          When PRTYPE = 3, specifies the distance between 2-by-2 */
00141 /*          blocks on the diagonal in A. Otherwise, QBLCKA is not */
00142 /*          referenced. QBLCKA > 1. */
00143 
00144 /*  QBLCKB  (input) INTEGER */
00145 /*          When PRTYPE = 3, specifies the distance between 2-by-2 */
00146 /*          blocks on the diagonal in B. Otherwise, QBLCKB is not */
00147 /*          referenced. QBLCKB > 1. */
00148 
00149 
00150 /*  Further Details */
00151 /*  =============== */
00152 
00153 /*  PRTYPE = 1: A and B are Jordan blocks, D and E are identity matrices */
00154 
00155 /*             A : if (i == j) then A(i, j) = 1.0 */
00156 /*                 if (j == i + 1) then A(i, j) = -1.0 */
00157 /*                 else A(i, j) = 0.0,            i, j = 1...M */
00158 
00159 /*             B : if (i == j) then B(i, j) = 1.0 - ALPHA */
00160 /*                 if (j == i + 1) then B(i, j) = 1.0 */
00161 /*                 else B(i, j) = 0.0,            i, j = 1...N */
00162 
00163 /*             D : if (i == j) then D(i, j) = 1.0 */
00164 /*                 else D(i, j) = 0.0,            i, j = 1...M */
00165 
00166 /*             E : if (i == j) then E(i, j) = 1.0 */
00167 /*                 else E(i, j) = 0.0,            i, j = 1...N */
00168 
00169 /*             L =  R are chosen from [-10...10], */
00170 /*                  which specifies the right hand sides (C, F). */
00171 
00172 /*  PRTYPE = 2 or 3: Triangular and/or quasi- triangular. */
00173 
00174 /*             A : if (i <= j) then A(i, j) = [-1...1] */
00175 /*                 else A(i, j) = 0.0,             i, j = 1...M */
00176 
00177 /*                 if (PRTYPE = 3) then */
00178 /*                    A(k + 1, k + 1) = A(k, k) */
00179 /*                    A(k + 1, k) = [-1...1] */
00180 /*                    sign(A(k, k + 1) = -(sin(A(k + 1, k)) */
00181 /*                        k = 1, M - 1, QBLCKA */
00182 
00183 /*             B : if (i <= j) then B(i, j) = [-1...1] */
00184 /*                 else B(i, j) = 0.0,            i, j = 1...N */
00185 
00186 /*                 if (PRTYPE = 3) then */
00187 /*                    B(k + 1, k + 1) = B(k, k) */
00188 /*                    B(k + 1, k) = [-1...1] */
00189 /*                    sign(B(k, k + 1) = -(sign(B(k + 1, k)) */
00190 /*                        k = 1, N - 1, QBLCKB */
00191 
00192 /*             D : if (i <= j) then D(i, j) = [-1...1]. */
00193 /*                 else D(i, j) = 0.0,            i, j = 1...M */
00194 
00195 
00196 /*             E : if (i <= j) then D(i, j) = [-1...1] */
00197 /*                 else E(i, j) = 0.0,            i, j = 1...N */
00198 
00199 /*                 L, R are chosen from [-10...10], */
00200 /*                 which specifies the right hand sides (C, F). */
00201 
00202 /*  PRTYPE = 4 Full */
00203 /*             A(i, j) = [-10...10] */
00204 /*             D(i, j) = [-1...1]    i,j = 1...M */
00205 /*             B(i, j) = [-10...10] */
00206 /*             E(i, j) = [-1...1]    i,j = 1...N */
00207 /*             R(i, j) = [-10...10] */
00208 /*             L(i, j) = [-1...1]    i = 1..M ,j = 1...N */
00209 
00210 /*             L, R specifies the right hand sides (C, F). */
00211 
00212 /*  PRTYPE = 5 special case common and/or close eigs. */
00213 
00214 /*  ===================================================================== */
00215 
00216 /*     .. Parameters .. */
00217 /*     .. */
00218 /*     .. Local Scalars .. */
00219 /*     .. */
00220 /*     .. Intrinsic Functions .. */
00221 /*     .. */
00222 /*     .. External Subroutines .. */
00223 /*     .. */
00224 /*     .. Executable Statements .. */
00225 
00226     /* Parameter adjustments */
00227     a_dim1 = *lda;
00228     a_offset = 1 + a_dim1;
00229     a -= a_offset;
00230     b_dim1 = *ldb;
00231     b_offset = 1 + b_dim1;
00232     b -= b_offset;
00233     c_dim1 = *ldc;
00234     c_offset = 1 + c_dim1;
00235     c__ -= c_offset;
00236     d_dim1 = *ldd;
00237     d_offset = 1 + d_dim1;
00238     d__ -= d_offset;
00239     e_dim1 = *lde;
00240     e_offset = 1 + e_dim1;
00241     e -= e_offset;
00242     f_dim1 = *ldf;
00243     f_offset = 1 + f_dim1;
00244     f -= f_offset;
00245     r_dim1 = *ldr;
00246     r_offset = 1 + r_dim1;
00247     r__ -= r_offset;
00248     l_dim1 = *ldl;
00249     l_offset = 1 + l_dim1;
00250     l -= l_offset;
00251 
00252     /* Function Body */
00253     if (*prtype == 1) {
00254         i__1 = *m;
00255         for (i__ = 1; i__ <= i__1; ++i__) {
00256             i__2 = *m;
00257             for (j = 1; j <= i__2; ++j) {
00258                 if (i__ == j) {
00259                     i__3 = i__ + j * a_dim1;
00260                     a[i__3].r = 1., a[i__3].i = 0.;
00261                     i__3 = i__ + j * d_dim1;
00262                     d__[i__3].r = 1., d__[i__3].i = 0.;
00263                 } else if (i__ == j - 1) {
00264                     i__3 = i__ + j * a_dim1;
00265                     z__1.r = -1., z__1.i = -0.;
00266                     a[i__3].r = z__1.r, a[i__3].i = z__1.i;
00267                     i__3 = i__ + j * d_dim1;
00268                     d__[i__3].r = 0., d__[i__3].i = 0.;
00269                 } else {
00270                     i__3 = i__ + j * a_dim1;
00271                     a[i__3].r = 0., a[i__3].i = 0.;
00272                     i__3 = i__ + j * d_dim1;
00273                     d__[i__3].r = 0., d__[i__3].i = 0.;
00274                 }
00275 /* L10: */
00276             }
00277 /* L20: */
00278         }
00279 
00280         i__1 = *n;
00281         for (i__ = 1; i__ <= i__1; ++i__) {
00282             i__2 = *n;
00283             for (j = 1; j <= i__2; ++j) {
00284                 if (i__ == j) {
00285                     i__3 = i__ + j * b_dim1;
00286                     z__1.r = 1. - *alpha, z__1.i = 0.;
00287                     b[i__3].r = z__1.r, b[i__3].i = z__1.i;
00288                     i__3 = i__ + j * e_dim1;
00289                     e[i__3].r = 1., e[i__3].i = 0.;
00290                 } else if (i__ == j - 1) {
00291                     i__3 = i__ + j * b_dim1;
00292                     b[i__3].r = 1., b[i__3].i = 0.;
00293                     i__3 = i__ + j * e_dim1;
00294                     e[i__3].r = 0., e[i__3].i = 0.;
00295                 } else {
00296                     i__3 = i__ + j * b_dim1;
00297                     b[i__3].r = 0., b[i__3].i = 0.;
00298                     i__3 = i__ + j * e_dim1;
00299                     e[i__3].r = 0., e[i__3].i = 0.;
00300                 }
00301 /* L30: */
00302             }
00303 /* L40: */
00304         }
00305 
00306         i__1 = *m;
00307         for (i__ = 1; i__ <= i__1; ++i__) {
00308             i__2 = *n;
00309             for (j = 1; j <= i__2; ++j) {
00310                 i__3 = i__ + j * r_dim1;
00311                 i__4 = i__ / j;
00312                 z__4.r = (doublereal) i__4, z__4.i = 0.;
00313                 z_sin(&z__3, &z__4);
00314                 z__2.r = .5 - z__3.r, z__2.i = 0. - z__3.i;
00315                 z__1.r = z__2.r * 20. - z__2.i * 0., z__1.i = z__2.r * 0. + 
00316                         z__2.i * 20.;
00317                 r__[i__3].r = z__1.r, r__[i__3].i = z__1.i;
00318                 i__3 = i__ + j * l_dim1;
00319                 i__4 = i__ + j * r_dim1;
00320                 l[i__3].r = r__[i__4].r, l[i__3].i = r__[i__4].i;
00321 /* L50: */
00322             }
00323 /* L60: */
00324         }
00325 
00326     } else if (*prtype == 2 || *prtype == 3) {
00327         i__1 = *m;
00328         for (i__ = 1; i__ <= i__1; ++i__) {
00329             i__2 = *m;
00330             for (j = 1; j <= i__2; ++j) {
00331                 if (i__ <= j) {
00332                     i__3 = i__ + j * a_dim1;
00333                     z__4.r = (doublereal) i__, z__4.i = 0.;
00334                     z_sin(&z__3, &z__4);
00335                     z__2.r = .5 - z__3.r, z__2.i = 0. - z__3.i;
00336                     z__1.r = z__2.r * 2. - z__2.i * 0., z__1.i = z__2.r * 0. 
00337                             + z__2.i * 2.;
00338                     a[i__3].r = z__1.r, a[i__3].i = z__1.i;
00339                     i__3 = i__ + j * d_dim1;
00340                     i__4 = i__ * j;
00341                     z__4.r = (doublereal) i__4, z__4.i = 0.;
00342                     z_sin(&z__3, &z__4);
00343                     z__2.r = .5 - z__3.r, z__2.i = 0. - z__3.i;
00344                     z__1.r = z__2.r * 2. - z__2.i * 0., z__1.i = z__2.r * 0. 
00345                             + z__2.i * 2.;
00346                     d__[i__3].r = z__1.r, d__[i__3].i = z__1.i;
00347                 } else {
00348                     i__3 = i__ + j * a_dim1;
00349                     a[i__3].r = 0., a[i__3].i = 0.;
00350                     i__3 = i__ + j * d_dim1;
00351                     d__[i__3].r = 0., d__[i__3].i = 0.;
00352                 }
00353 /* L70: */
00354             }
00355 /* L80: */
00356         }
00357 
00358         i__1 = *n;
00359         for (i__ = 1; i__ <= i__1; ++i__) {
00360             i__2 = *n;
00361             for (j = 1; j <= i__2; ++j) {
00362                 if (i__ <= j) {
00363                     i__3 = i__ + j * b_dim1;
00364                     i__4 = i__ + j;
00365                     z__4.r = (doublereal) i__4, z__4.i = 0.;
00366                     z_sin(&z__3, &z__4);
00367                     z__2.r = .5 - z__3.r, z__2.i = 0. - z__3.i;
00368                     z__1.r = z__2.r * 2. - z__2.i * 0., z__1.i = z__2.r * 0. 
00369                             + z__2.i * 2.;
00370                     b[i__3].r = z__1.r, b[i__3].i = z__1.i;
00371                     i__3 = i__ + j * e_dim1;
00372                     z__4.r = (doublereal) j, z__4.i = 0.;
00373                     z_sin(&z__3, &z__4);
00374                     z__2.r = .5 - z__3.r, z__2.i = 0. - z__3.i;
00375                     z__1.r = z__2.r * 2. - z__2.i * 0., z__1.i = z__2.r * 0. 
00376                             + z__2.i * 2.;
00377                     e[i__3].r = z__1.r, e[i__3].i = z__1.i;
00378                 } else {
00379                     i__3 = i__ + j * b_dim1;
00380                     b[i__3].r = 0., b[i__3].i = 0.;
00381                     i__3 = i__ + j * e_dim1;
00382                     e[i__3].r = 0., e[i__3].i = 0.;
00383                 }
00384 /* L90: */
00385             }
00386 /* L100: */
00387         }
00388 
00389         i__1 = *m;
00390         for (i__ = 1; i__ <= i__1; ++i__) {
00391             i__2 = *n;
00392             for (j = 1; j <= i__2; ++j) {
00393                 i__3 = i__ + j * r_dim1;
00394                 i__4 = i__ * j;
00395                 z__4.r = (doublereal) i__4, z__4.i = 0.;
00396                 z_sin(&z__3, &z__4);
00397                 z__2.r = .5 - z__3.r, z__2.i = 0. - z__3.i;
00398                 z__1.r = z__2.r * 20. - z__2.i * 0., z__1.i = z__2.r * 0. + 
00399                         z__2.i * 20.;
00400                 r__[i__3].r = z__1.r, r__[i__3].i = z__1.i;
00401                 i__3 = i__ + j * l_dim1;
00402                 i__4 = i__ + j;
00403                 z__4.r = (doublereal) i__4, z__4.i = 0.;
00404                 z_sin(&z__3, &z__4);
00405                 z__2.r = .5 - z__3.r, z__2.i = 0. - z__3.i;
00406                 z__1.r = z__2.r * 20. - z__2.i * 0., z__1.i = z__2.r * 0. + 
00407                         z__2.i * 20.;
00408                 l[i__3].r = z__1.r, l[i__3].i = z__1.i;
00409 /* L110: */
00410             }
00411 /* L120: */
00412         }
00413 
00414         if (*prtype == 3) {
00415             if (*qblcka <= 1) {
00416                 *qblcka = 2;
00417             }
00418             i__1 = *m - 1;
00419             i__2 = *qblcka;
00420             for (k = 1; i__2 < 0 ? k >= i__1 : k <= i__1; k += i__2) {
00421                 i__3 = k + 1 + (k + 1) * a_dim1;
00422                 i__4 = k + k * a_dim1;
00423                 a[i__3].r = a[i__4].r, a[i__3].i = a[i__4].i;
00424                 i__3 = k + 1 + k * a_dim1;
00425                 z_sin(&z__2, &a[k + (k + 1) * a_dim1]);
00426                 z__1.r = -z__2.r, z__1.i = -z__2.i;
00427                 a[i__3].r = z__1.r, a[i__3].i = z__1.i;
00428 /* L130: */
00429             }
00430 
00431             if (*qblckb <= 1) {
00432                 *qblckb = 2;
00433             }
00434             i__2 = *n - 1;
00435             i__1 = *qblckb;
00436             for (k = 1; i__1 < 0 ? k >= i__2 : k <= i__2; k += i__1) {
00437                 i__3 = k + 1 + (k + 1) * b_dim1;
00438                 i__4 = k + k * b_dim1;
00439                 b[i__3].r = b[i__4].r, b[i__3].i = b[i__4].i;
00440                 i__3 = k + 1 + k * b_dim1;
00441                 z_sin(&z__2, &b[k + (k + 1) * b_dim1]);
00442                 z__1.r = -z__2.r, z__1.i = -z__2.i;
00443                 b[i__3].r = z__1.r, b[i__3].i = z__1.i;
00444 /* L140: */
00445             }
00446         }
00447 
00448     } else if (*prtype == 4) {
00449         i__1 = *m;
00450         for (i__ = 1; i__ <= i__1; ++i__) {
00451             i__2 = *m;
00452             for (j = 1; j <= i__2; ++j) {
00453                 i__3 = i__ + j * a_dim1;
00454                 i__4 = i__ * j;
00455                 z__4.r = (doublereal) i__4, z__4.i = 0.;
00456                 z_sin(&z__3, &z__4);
00457                 z__2.r = .5 - z__3.r, z__2.i = 0. - z__3.i;
00458                 z__1.r = z__2.r * 20. - z__2.i * 0., z__1.i = z__2.r * 0. + 
00459                         z__2.i * 20.;
00460                 a[i__3].r = z__1.r, a[i__3].i = z__1.i;
00461                 i__3 = i__ + j * d_dim1;
00462                 i__4 = i__ + j;
00463                 z__4.r = (doublereal) i__4, z__4.i = 0.;
00464                 z_sin(&z__3, &z__4);
00465                 z__2.r = .5 - z__3.r, z__2.i = 0. - z__3.i;
00466                 z__1.r = z__2.r * 2. - z__2.i * 0., z__1.i = z__2.r * 0. + 
00467                         z__2.i * 2.;
00468                 d__[i__3].r = z__1.r, d__[i__3].i = z__1.i;
00469 /* L150: */
00470             }
00471 /* L160: */
00472         }
00473 
00474         i__1 = *n;
00475         for (i__ = 1; i__ <= i__1; ++i__) {
00476             i__2 = *n;
00477             for (j = 1; j <= i__2; ++j) {
00478                 i__3 = i__ + j * b_dim1;
00479                 i__4 = i__ + j;
00480                 z__4.r = (doublereal) i__4, z__4.i = 0.;
00481                 z_sin(&z__3, &z__4);
00482                 z__2.r = .5 - z__3.r, z__2.i = 0. - z__3.i;
00483                 z__1.r = z__2.r * 20. - z__2.i * 0., z__1.i = z__2.r * 0. + 
00484                         z__2.i * 20.;
00485                 b[i__3].r = z__1.r, b[i__3].i = z__1.i;
00486                 i__3 = i__ + j * e_dim1;
00487                 i__4 = i__ * j;
00488                 z__4.r = (doublereal) i__4, z__4.i = 0.;
00489                 z_sin(&z__3, &z__4);
00490                 z__2.r = .5 - z__3.r, z__2.i = 0. - z__3.i;
00491                 z__1.r = z__2.r * 2. - z__2.i * 0., z__1.i = z__2.r * 0. + 
00492                         z__2.i * 2.;
00493                 e[i__3].r = z__1.r, e[i__3].i = z__1.i;
00494 /* L170: */
00495             }
00496 /* L180: */
00497         }
00498 
00499         i__1 = *m;
00500         for (i__ = 1; i__ <= i__1; ++i__) {
00501             i__2 = *n;
00502             for (j = 1; j <= i__2; ++j) {
00503                 i__3 = i__ + j * r_dim1;
00504                 i__4 = j / i__;
00505                 z__4.r = (doublereal) i__4, z__4.i = 0.;
00506                 z_sin(&z__3, &z__4);
00507                 z__2.r = .5 - z__3.r, z__2.i = 0. - z__3.i;
00508                 z__1.r = z__2.r * 20. - z__2.i * 0., z__1.i = z__2.r * 0. + 
00509                         z__2.i * 20.;
00510                 r__[i__3].r = z__1.r, r__[i__3].i = z__1.i;
00511                 i__3 = i__ + j * l_dim1;
00512                 i__4 = i__ * j;
00513                 z__4.r = (doublereal) i__4, z__4.i = 0.;
00514                 z_sin(&z__3, &z__4);
00515                 z__2.r = .5 - z__3.r, z__2.i = 0. - z__3.i;
00516                 z__1.r = z__2.r * 2. - z__2.i * 0., z__1.i = z__2.r * 0. + 
00517                         z__2.i * 2.;
00518                 l[i__3].r = z__1.r, l[i__3].i = z__1.i;
00519 /* L190: */
00520             }
00521 /* L200: */
00522         }
00523 
00524     } else if (*prtype >= 5) {
00525         z__3.r = 1., z__3.i = 0.;
00526         z__2.r = z__3.r * 20. - z__3.i * 0., z__2.i = z__3.r * 0. + z__3.i * 
00527                 20.;
00528         z__1.r = z__2.r / *alpha, z__1.i = z__2.i / *alpha;
00529         reeps.r = z__1.r, reeps.i = z__1.i;
00530         z__2.r = -1.5, z__2.i = 0.;
00531         z__1.r = z__2.r / *alpha, z__1.i = z__2.i / *alpha;
00532         imeps.r = z__1.r, imeps.i = z__1.i;
00533         i__1 = *m;
00534         for (i__ = 1; i__ <= i__1; ++i__) {
00535             i__2 = *n;
00536             for (j = 1; j <= i__2; ++j) {
00537                 i__3 = i__ + j * r_dim1;
00538                 i__4 = i__ * j;
00539                 z__5.r = (doublereal) i__4, z__5.i = 0.;
00540                 z_sin(&z__4, &z__5);
00541                 z__3.r = .5 - z__4.r, z__3.i = 0. - z__4.i;
00542                 z__2.r = *alpha * z__3.r, z__2.i = *alpha * z__3.i;
00543                 z_div(&z__1, &z__2, &c_b5);
00544                 r__[i__3].r = z__1.r, r__[i__3].i = z__1.i;
00545                 i__3 = i__ + j * l_dim1;
00546                 i__4 = i__ + j;
00547                 z__5.r = (doublereal) i__4, z__5.i = 0.;
00548                 z_sin(&z__4, &z__5);
00549                 z__3.r = .5 - z__4.r, z__3.i = 0. - z__4.i;
00550                 z__2.r = *alpha * z__3.r, z__2.i = *alpha * z__3.i;
00551                 z_div(&z__1, &z__2, &c_b5);
00552                 l[i__3].r = z__1.r, l[i__3].i = z__1.i;
00553 /* L210: */
00554             }
00555 /* L220: */
00556         }
00557 
00558         i__1 = *m;
00559         for (i__ = 1; i__ <= i__1; ++i__) {
00560             i__2 = i__ + i__ * d_dim1;
00561             d__[i__2].r = 1., d__[i__2].i = 0.;
00562 /* L230: */
00563         }
00564 
00565         i__1 = *m;
00566         for (i__ = 1; i__ <= i__1; ++i__) {
00567             if (i__ <= 4) {
00568                 i__2 = i__ + i__ * a_dim1;
00569                 a[i__2].r = 1., a[i__2].i = 0.;
00570                 if (i__ > 2) {
00571                     i__2 = i__ + i__ * a_dim1;
00572                     z__1.r = reeps.r + 1., z__1.i = reeps.i + 0.;
00573                     a[i__2].r = z__1.r, a[i__2].i = z__1.i;
00574                 }
00575                 if (i__ % 2 != 0 && i__ < *m) {
00576                     i__2 = i__ + (i__ + 1) * a_dim1;
00577                     a[i__2].r = imeps.r, a[i__2].i = imeps.i;
00578                 } else if (i__ > 1) {
00579                     i__2 = i__ + (i__ - 1) * a_dim1;
00580                     z__1.r = -imeps.r, z__1.i = -imeps.i;
00581                     a[i__2].r = z__1.r, a[i__2].i = z__1.i;
00582                 }
00583             } else if (i__ <= 8) {
00584                 if (i__ <= 6) {
00585                     i__2 = i__ + i__ * a_dim1;
00586                     a[i__2].r = reeps.r, a[i__2].i = reeps.i;
00587                 } else {
00588                     i__2 = i__ + i__ * a_dim1;
00589                     z__1.r = -reeps.r, z__1.i = -reeps.i;
00590                     a[i__2].r = z__1.r, a[i__2].i = z__1.i;
00591                 }
00592                 if (i__ % 2 != 0 && i__ < *m) {
00593                     i__2 = i__ + (i__ + 1) * a_dim1;
00594                     a[i__2].r = 1., a[i__2].i = 0.;
00595                 } else if (i__ > 1) {
00596                     i__2 = i__ + (i__ - 1) * a_dim1;
00597                     z__1.r = -1., z__1.i = -0.;
00598                     a[i__2].r = z__1.r, a[i__2].i = z__1.i;
00599                 }
00600             } else {
00601                 i__2 = i__ + i__ * a_dim1;
00602                 a[i__2].r = 1., a[i__2].i = 0.;
00603                 if (i__ % 2 != 0 && i__ < *m) {
00604                     i__2 = i__ + (i__ + 1) * a_dim1;
00605                     d__1 = 2.;
00606                     z__1.r = d__1 * imeps.r, z__1.i = d__1 * imeps.i;
00607                     a[i__2].r = z__1.r, a[i__2].i = z__1.i;
00608                 } else if (i__ > 1) {
00609                     i__2 = i__ + (i__ - 1) * a_dim1;
00610                     z__2.r = -imeps.r, z__2.i = -imeps.i;
00611                     d__1 = 2.;
00612                     z__1.r = d__1 * z__2.r, z__1.i = d__1 * z__2.i;
00613                     a[i__2].r = z__1.r, a[i__2].i = z__1.i;
00614                 }
00615             }
00616 /* L240: */
00617         }
00618 
00619         i__1 = *n;
00620         for (i__ = 1; i__ <= i__1; ++i__) {
00621             i__2 = i__ + i__ * e_dim1;
00622             e[i__2].r = 1., e[i__2].i = 0.;
00623             if (i__ <= 4) {
00624                 i__2 = i__ + i__ * b_dim1;
00625                 z__1.r = -1., z__1.i = -0.;
00626                 b[i__2].r = z__1.r, b[i__2].i = z__1.i;
00627                 if (i__ > 2) {
00628                     i__2 = i__ + i__ * b_dim1;
00629                     z__1.r = 1. - reeps.r, z__1.i = 0. - reeps.i;
00630                     b[i__2].r = z__1.r, b[i__2].i = z__1.i;
00631                 }
00632                 if (i__ % 2 != 0 && i__ < *n) {
00633                     i__2 = i__ + (i__ + 1) * b_dim1;
00634                     b[i__2].r = imeps.r, b[i__2].i = imeps.i;
00635                 } else if (i__ > 1) {
00636                     i__2 = i__ + (i__ - 1) * b_dim1;
00637                     z__1.r = -imeps.r, z__1.i = -imeps.i;
00638                     b[i__2].r = z__1.r, b[i__2].i = z__1.i;
00639                 }
00640             } else if (i__ <= 8) {
00641                 if (i__ <= 6) {
00642                     i__2 = i__ + i__ * b_dim1;
00643                     b[i__2].r = reeps.r, b[i__2].i = reeps.i;
00644                 } else {
00645                     i__2 = i__ + i__ * b_dim1;
00646                     z__1.r = -reeps.r, z__1.i = -reeps.i;
00647                     b[i__2].r = z__1.r, b[i__2].i = z__1.i;
00648                 }
00649                 if (i__ % 2 != 0 && i__ < *n) {
00650                     i__2 = i__ + (i__ + 1) * b_dim1;
00651                     z__1.r = imeps.r + 1., z__1.i = imeps.i + 0.;
00652                     b[i__2].r = z__1.r, b[i__2].i = z__1.i;
00653                 } else if (i__ > 1) {
00654                     i__2 = i__ + (i__ - 1) * b_dim1;
00655                     z__2.r = -1., z__2.i = -0.;
00656                     z__1.r = z__2.r - imeps.r, z__1.i = z__2.i - imeps.i;
00657                     b[i__2].r = z__1.r, b[i__2].i = z__1.i;
00658                 }
00659             } else {
00660                 i__2 = i__ + i__ * b_dim1;
00661                 z__1.r = 1. - reeps.r, z__1.i = 0. - reeps.i;
00662                 b[i__2].r = z__1.r, b[i__2].i = z__1.i;
00663                 if (i__ % 2 != 0 && i__ < *n) {
00664                     i__2 = i__ + (i__ + 1) * b_dim1;
00665                     d__1 = 2.;
00666                     z__1.r = d__1 * imeps.r, z__1.i = d__1 * imeps.i;
00667                     b[i__2].r = z__1.r, b[i__2].i = z__1.i;
00668                 } else if (i__ > 1) {
00669                     i__2 = i__ + (i__ - 1) * b_dim1;
00670                     z__2.r = -imeps.r, z__2.i = -imeps.i;
00671                     d__1 = 2.;
00672                     z__1.r = d__1 * z__2.r, z__1.i = d__1 * z__2.i;
00673                     b[i__2].r = z__1.r, b[i__2].i = z__1.i;
00674                 }
00675             }
00676 /* L250: */
00677         }
00678     }
00679 
00680 /*     Compute rhs (C, F) */
00681 
00682     zgemm_("N", "N", m, n, m, &c_b1, &a[a_offset], lda, &r__[r_offset], ldr, &
00683             c_b3, &c__[c_offset], ldc);
00684     z__1.r = -1., z__1.i = -0.;
00685     zgemm_("N", "N", m, n, n, &z__1, &l[l_offset], ldl, &b[b_offset], ldb, &
00686             c_b1, &c__[c_offset], ldc);
00687     zgemm_("N", "N", m, n, m, &c_b1, &d__[d_offset], ldd, &r__[r_offset], ldr, 
00688              &c_b3, &f[f_offset], ldf);
00689     z__1.r = -1., z__1.i = -0.;
00690     zgemm_("N", "N", m, n, n, &z__1, &l[l_offset], ldl, &e[e_offset], lde, &
00691             c_b1, &f[f_offset], ldf);
00692 
00693 /*     End of ZLATM5 */
00694 
00695     return 0;
00696 } /* zlatm5_ */


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