swiftnav: dlatm5.c Source File

Go to the documentation of this file.
00001 /* dlatm5.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 doublereal c_b29 = 1.;
00019 static doublereal c_b30 = 0.;
00020 static doublereal c_b33 = -1.;
00021 
00022 /* Subroutine */ int dlatm5_(integer *prtype, integer *m, integer *n, 
00023         doublereal *a, integer *lda, doublereal *b, integer *ldb, doublereal *
00024         c__, integer *ldc, doublereal *d__, integer *ldd, doublereal *e, 
00025         integer *lde, doublereal *f, integer *ldf, doublereal *r__, integer *
00026         ldr, doublereal *l, integer *ldl, doublereal *alpha, integer *qblcka, 
00027         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;
00033 
00034     /* Builtin functions */
00035     double sin(doublereal);
00036 
00037     /* Local variables */
00038     integer i__, j, k;
00039     extern /* Subroutine */ int dgemm_(char *, char *, integer *, integer *, 
00040             integer *, doublereal *, doublereal *, integer *, doublereal *, 
00041             integer *, doublereal *, doublereal *, integer *);
00042     doublereal imeps, reeps;
00043 
00044 
00045 /*  -- LAPACK test routine (version 3.1) -- */
00046 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00047 /*     November 2006 */
00048 
00049 /*     .. Scalar Arguments .. */
00050 /*     .. */
00051 /*     .. Array Arguments .. */
00052 /*     .. */
00053 
00054 /*  Purpose */
00055 /*  ======= */
00056 
00057 /*  DLATM5 generates matrices involved in the Generalized Sylvester */
00058 /*  equation: */
00059 
00060 /*      A * R - L * B = C */
00061 /*      D * R - L * E = F */
00062 
00063 /*  They also satisfy (the diagonalization condition) */
00064 
00065 /*   [ I -L ] ( [ A  -C ], [ D -F ] ) [ I  R ] = ( [ A    ], [ D    ] ) */
00066 /*   [    I ] ( [     B ]  [    E ] ) [    I ]   ( [    B ]  [    E ] ) */
00067 
00068 
00069 /*  Arguments */
00070 /*  ========= */
00071 
00072 /*  PRTYPE  (input) INTEGER */
00073 /*          "Points" to a certian type of the matrices to generate */
00074 /*          (see futher details). */
00075 
00076 /*  M       (input) INTEGER */
00077 /*          Specifies the order of A and D and the number of rows in */
00078 /*          C, F,  R and L. */
00079 
00080 /*  N       (input) INTEGER */
00081 /*          Specifies the order of B and E and the number of columns in */
00082 /*          C, F, R and L. */
00083 
00084 /*  A       (output) DOUBLE PRECISION array, dimension (LDA, M). */
00085 /*          On exit A M-by-M is initialized according to PRTYPE. */
00086 
00087 /*  LDA     (input) INTEGER */
00088 /*          The leading dimension of A. */
00089 
00090 /*  B       (output) DOUBLE PRECISION array, dimension (LDB, N). */
00091 /*          On exit B N-by-N is initialized according to PRTYPE. */
00092 
00093 /*  LDB     (input) INTEGER */
00094 /*          The leading dimension of B. */
00095 
00096 /*  C       (output) DOUBLE PRECISION array, dimension (LDC, N). */
00097 /*          On exit C M-by-N is initialized according to PRTYPE. */
00098 
00099 /*  LDC     (input) INTEGER */
00100 /*          The leading dimension of C. */
00101 
00102 /*  D       (output) DOUBLE PRECISION array, dimension (LDD, M). */
00103 /*          On exit D M-by-M is initialized according to PRTYPE. */
00104 
00105 /*  LDD     (input) INTEGER */
00106 /*          The leading dimension of D. */
00107 
00108 /*  E       (output) DOUBLE PRECISION array, dimension (LDE, N). */
00109 /*          On exit E N-by-N is initialized according to PRTYPE. */
00110 
00111 /*  LDE     (input) INTEGER */
00112 /*          The leading dimension of E. */
00113 
00114 /*  F       (output) DOUBLE PRECISION array, dimension (LDF, N). */
00115 /*          On exit F M-by-N is initialized according to PRTYPE. */
00116 
00117 /*  LDF     (input) INTEGER */
00118 /*          The leading dimension of F. */
00119 
00120 /*  R       (output) DOUBLE PRECISION array, dimension (LDR, N). */
00121 /*          On exit R M-by-N is initialized according to PRTYPE. */
00122 
00123 /*  LDR     (input) INTEGER */
00124 /*          The leading dimension of R. */
00125 
00126 /*  L       (output) DOUBLE PRECISION array, dimension (LDL, N). */
00127 /*          On exit L M-by-N is initialized according to PRTYPE. */
00128 
00129 /*  LDL     (input) INTEGER */
00130 /*          The leading dimension of L. */
00131 
00132 /*  ALPHA   (input) DOUBLE PRECISION */
00133 /*          Parameter used in generating PRTYPE = 1 and 5 matrices. */
00134 
00135 /*  QBLCKA  (input) INTEGER */
00136 /*          When PRTYPE = 3, specifies the distance between 2-by-2 */
00137 /*          blocks on the diagonal in A. Otherwise, QBLCKA is not */
00138 /*          referenced. QBLCKA > 1. */
00139 
00140 /*  QBLCKB  (input) INTEGER */
00141 /*          When PRTYPE = 3, specifies the distance between 2-by-2 */
00142 /*          blocks on the diagonal in B. Otherwise, QBLCKB is not */
00143 /*          referenced. QBLCKB > 1. */
00144 
00145 
00146 /*  Further Details */
00147 /*  =============== */
00148 
00149 /*  PRTYPE = 1: A and B are Jordan blocks, D and E are identity matrices */
00150 
00151 /*             A : if (i == j) then A(i, j) = 1.0 */
00152 /*                 if (j == i + 1) then A(i, j) = -1.0 */
00153 /*                 else A(i, j) = 0.0,            i, j = 1...M */
00154 
00155 /*             B : if (i == j) then B(i, j) = 1.0 - ALPHA */
00156 /*                 if (j == i + 1) then B(i, j) = 1.0 */
00157 /*                 else B(i, j) = 0.0,            i, j = 1...N */
00158 
00159 /*             D : if (i == j) then D(i, j) = 1.0 */
00160 /*                 else D(i, j) = 0.0,            i, j = 1...M */
00161 
00162 /*             E : if (i == j) then E(i, j) = 1.0 */
00163 /*                 else E(i, j) = 0.0,            i, j = 1...N */
00164 
00165 /*             L =  R are chosen from [-10...10], */
00166 /*                  which specifies the right hand sides (C, F). */
00167 
00168 /*  PRTYPE = 2 or 3: Triangular and/or quasi- triangular. */
00169 
00170 /*             A : if (i <= j) then A(i, j) = [-1...1] */
00171 /*                 else A(i, j) = 0.0,             i, j = 1...M */
00172 
00173 /*                 if (PRTYPE = 3) then */
00174 /*                    A(k + 1, k + 1) = A(k, k) */
00175 /*                    A(k + 1, k) = [-1...1] */
00176 /*                    sign(A(k, k + 1) = -(sin(A(k + 1, k)) */
00177 /*                        k = 1, M - 1, QBLCKA */
00178 
00179 /*             B : if (i <= j) then B(i, j) = [-1...1] */
00180 /*                 else B(i, j) = 0.0,            i, j = 1...N */
00181 
00182 /*                 if (PRTYPE = 3) then */
00183 /*                    B(k + 1, k + 1) = B(k, k) */
00184 /*                    B(k + 1, k) = [-1...1] */
00185 /*                    sign(B(k, k + 1) = -(sign(B(k + 1, k)) */
00186 /*                        k = 1, N - 1, QBLCKB */
00187 
00188 /*             D : if (i <= j) then D(i, j) = [-1...1]. */
00189 /*                 else D(i, j) = 0.0,            i, j = 1...M */
00190 
00191 
00192 /*             E : if (i <= j) then D(i, j) = [-1...1] */
00193 /*                 else E(i, j) = 0.0,            i, j = 1...N */
00194 
00195 /*                 L, R are chosen from [-10...10], */
00196 /*                 which specifies the right hand sides (C, F). */
00197 
00198 /*  PRTYPE = 4 Full */
00199 /*             A(i, j) = [-10...10] */
00200 /*             D(i, j) = [-1...1]    i,j = 1...M */
00201 /*             B(i, j) = [-10...10] */
00202 /*             E(i, j) = [-1...1]    i,j = 1...N */
00203 /*             R(i, j) = [-10...10] */
00204 /*             L(i, j) = [-1...1]    i = 1..M ,j = 1...N */
00205 
00206 /*             L, R specifies the right hand sides (C, F). */
00207 
00208 /*  PRTYPE = 5 special case common and/or close eigs. */
00209 
00210 /*  ===================================================================== */
00211 
00212 /*     .. Parameters .. */
00213 /*     .. */
00214 /*     .. Local Scalars .. */
00215 /*     .. */
00216 /*     .. Intrinsic Functions .. */
00217 /*     .. */
00218 /*     .. External Subroutines .. */
00219 /*     .. */
00220 /*     .. Executable Statements .. */
00221 
00222     /* Parameter adjustments */
00223     a_dim1 = *lda;
00224     a_offset = 1 + a_dim1;
00225     a -= a_offset;
00226     b_dim1 = *ldb;
00227     b_offset = 1 + b_dim1;
00228     b -= b_offset;
00229     c_dim1 = *ldc;
00230     c_offset = 1 + c_dim1;
00231     c__ -= c_offset;
00232     d_dim1 = *ldd;
00233     d_offset = 1 + d_dim1;
00234     d__ -= d_offset;
00235     e_dim1 = *lde;
00236     e_offset = 1 + e_dim1;
00237     e -= e_offset;
00238     f_dim1 = *ldf;
00239     f_offset = 1 + f_dim1;
00240     f -= f_offset;
00241     r_dim1 = *ldr;
00242     r_offset = 1 + r_dim1;
00243     r__ -= r_offset;
00244     l_dim1 = *ldl;
00245     l_offset = 1 + l_dim1;
00246     l -= l_offset;
00247 
00248     /* Function Body */
00249     if (*prtype == 1) {
00250         i__1 = *m;
00251         for (i__ = 1; i__ <= i__1; ++i__) {
00252             i__2 = *m;
00253             for (j = 1; j <= i__2; ++j) {
00254                 if (i__ == j) {
00255                     a[i__ + j * a_dim1] = 1.;
00256                     d__[i__ + j * d_dim1] = 1.;
00257                 } else if (i__ == j - 1) {
00258                     a[i__ + j * a_dim1] = -1.;
00259                     d__[i__ + j * d_dim1] = 0.;
00260                 } else {
00261                     a[i__ + j * a_dim1] = 0.;
00262                     d__[i__ + j * d_dim1] = 0.;
00263                 }
00264 /* L10: */
00265             }
00266 /* L20: */
00267         }
00268 
00269         i__1 = *n;
00270         for (i__ = 1; i__ <= i__1; ++i__) {
00271             i__2 = *n;
00272             for (j = 1; j <= i__2; ++j) {
00273                 if (i__ == j) {
00274                     b[i__ + j * b_dim1] = 1. - *alpha;
00275                     e[i__ + j * e_dim1] = 1.;
00276                 } else if (i__ == j - 1) {
00277                     b[i__ + j * b_dim1] = 1.;
00278                     e[i__ + j * e_dim1] = 0.;
00279                 } else {
00280                     b[i__ + j * b_dim1] = 0.;
00281                     e[i__ + j * e_dim1] = 0.;
00282                 }
00283 /* L30: */
00284             }
00285 /* L40: */
00286         }
00287 
00288         i__1 = *m;
00289         for (i__ = 1; i__ <= i__1; ++i__) {
00290             i__2 = *n;
00291             for (j = 1; j <= i__2; ++j) {
00292                 r__[i__ + j * r_dim1] = (.5 - sin((doublereal) (i__ / j))) * 
00293                         20.;
00294                 l[i__ + j * l_dim1] = r__[i__ + j * r_dim1];
00295 /* L50: */
00296             }
00297 /* L60: */
00298         }
00299 
00300     } else if (*prtype == 2 || *prtype == 3) {
00301         i__1 = *m;
00302         for (i__ = 1; i__ <= i__1; ++i__) {
00303             i__2 = *m;
00304             for (j = 1; j <= i__2; ++j) {
00305                 if (i__ <= j) {
00306                     a[i__ + j * a_dim1] = (.5 - sin((doublereal) i__)) * 2.;
00307                     d__[i__ + j * d_dim1] = (.5 - sin((doublereal) (i__ * j)))
00308                              * 2.;
00309                 } else {
00310                     a[i__ + j * a_dim1] = 0.;
00311                     d__[i__ + j * d_dim1] = 0.;
00312                 }
00313 /* L70: */
00314             }
00315 /* L80: */
00316         }
00317 
00318         i__1 = *n;
00319         for (i__ = 1; i__ <= i__1; ++i__) {
00320             i__2 = *n;
00321             for (j = 1; j <= i__2; ++j) {
00322                 if (i__ <= j) {
00323                     b[i__ + j * b_dim1] = (.5 - sin((doublereal) (i__ + j))) *
00324                              2.;
00325                     e[i__ + j * e_dim1] = (.5 - sin((doublereal) j)) * 2.;
00326                 } else {
00327                     b[i__ + j * b_dim1] = 0.;
00328                     e[i__ + j * e_dim1] = 0.;
00329                 }
00330 /* L90: */
00331             }
00332 /* L100: */
00333         }
00334 
00335         i__1 = *m;
00336         for (i__ = 1; i__ <= i__1; ++i__) {
00337             i__2 = *n;
00338             for (j = 1; j <= i__2; ++j) {
00339                 r__[i__ + j * r_dim1] = (.5 - sin((doublereal) (i__ * j))) * 
00340                         20.;
00341                 l[i__ + j * l_dim1] = (.5 - sin((doublereal) (i__ + j))) * 
00342                         20.;
00343 /* L110: */
00344             }
00345 /* L120: */
00346         }
00347 
00348         if (*prtype == 3) {
00349             if (*qblcka <= 1) {
00350                 *qblcka = 2;
00351             }
00352             i__1 = *m - 1;
00353             i__2 = *qblcka;
00354             for (k = 1; i__2 < 0 ? k >= i__1 : k <= i__1; k += i__2) {
00355                 a[k + 1 + (k + 1) * a_dim1] = a[k + k * a_dim1];
00356                 a[k + 1 + k * a_dim1] = -sin(a[k + (k + 1) * a_dim1]);
00357 /* L130: */
00358             }
00359 
00360             if (*qblckb <= 1) {
00361                 *qblckb = 2;
00362             }
00363             i__2 = *n - 1;
00364             i__1 = *qblckb;
00365             for (k = 1; i__1 < 0 ? k >= i__2 : k <= i__2; k += i__1) {
00366                 b[k + 1 + (k + 1) * b_dim1] = b[k + k * b_dim1];
00367                 b[k + 1 + k * b_dim1] = -sin(b[k + (k + 1) * b_dim1]);
00368 /* L140: */
00369             }
00370         }
00371 
00372     } else if (*prtype == 4) {
00373         i__1 = *m;
00374         for (i__ = 1; i__ <= i__1; ++i__) {
00375             i__2 = *m;
00376             for (j = 1; j <= i__2; ++j) {
00377                 a[i__ + j * a_dim1] = (.5 - sin((doublereal) (i__ * j))) * 
00378                         20.;
00379                 d__[i__ + j * d_dim1] = (.5 - sin((doublereal) (i__ + j))) * 
00380                         2.;
00381 /* L150: */
00382             }
00383 /* L160: */
00384         }
00385 
00386         i__1 = *n;
00387         for (i__ = 1; i__ <= i__1; ++i__) {
00388             i__2 = *n;
00389             for (j = 1; j <= i__2; ++j) {
00390                 b[i__ + j * b_dim1] = (.5 - sin((doublereal) (i__ + j))) * 
00391                         20.;
00392                 e[i__ + j * e_dim1] = (.5 - sin((doublereal) (i__ * j))) * 2.;
00393 /* L170: */
00394             }
00395 /* L180: */
00396         }
00397 
00398         i__1 = *m;
00399         for (i__ = 1; i__ <= i__1; ++i__) {
00400             i__2 = *n;
00401             for (j = 1; j <= i__2; ++j) {
00402                 r__[i__ + j * r_dim1] = (.5 - sin((doublereal) (j / i__))) * 
00403                         20.;
00404                 l[i__ + j * l_dim1] = (.5 - sin((doublereal) (i__ * j))) * 2.;
00405 /* L190: */
00406             }
00407 /* L200: */
00408         }
00409 
00410     } else if (*prtype >= 5) {
00411         reeps = 20. / *alpha;
00412         imeps = -1.5 / *alpha;
00413         i__1 = *m;
00414         for (i__ = 1; i__ <= i__1; ++i__) {
00415             i__2 = *n;
00416             for (j = 1; j <= i__2; ++j) {
00417                 r__[i__ + j * r_dim1] = (.5 - sin((doublereal) (i__ * j))) * *
00418                         alpha / 20.;
00419                 l[i__ + j * l_dim1] = (.5 - sin((doublereal) (i__ + j))) * *
00420                         alpha / 20.;
00421 /* L210: */
00422             }
00423 /* L220: */
00424         }
00425 
00426         i__1 = *m;
00427         for (i__ = 1; i__ <= i__1; ++i__) {
00428             d__[i__ + i__ * d_dim1] = 1.;
00429 /* L230: */
00430         }
00431 
00432         i__1 = *m;
00433         for (i__ = 1; i__ <= i__1; ++i__) {
00434             if (i__ <= 4) {
00435                 a[i__ + i__ * a_dim1] = 1.;
00436                 if (i__ > 2) {
00437                     a[i__ + i__ * a_dim1] = reeps + 1.;
00438                 }
00439                 if (i__ % 2 != 0 && i__ < *m) {
00440                     a[i__ + (i__ + 1) * a_dim1] = imeps;
00441                 } else if (i__ > 1) {
00442                     a[i__ + (i__ - 1) * a_dim1] = -imeps;
00443                 }
00444             } else if (i__ <= 8) {
00445                 if (i__ <= 6) {
00446                     a[i__ + i__ * a_dim1] = reeps;
00447                 } else {
00448                     a[i__ + i__ * a_dim1] = -reeps;
00449                 }
00450                 if (i__ % 2 != 0 && i__ < *m) {
00451                     a[i__ + (i__ + 1) * a_dim1] = 1.;
00452                 } else if (i__ > 1) {
00453                     a[i__ + (i__ - 1) * a_dim1] = -1.;
00454                 }
00455             } else {
00456                 a[i__ + i__ * a_dim1] = 1.;
00457                 if (i__ % 2 != 0 && i__ < *m) {
00458                     a[i__ + (i__ + 1) * a_dim1] = imeps * 2;
00459                 } else if (i__ > 1) {
00460                     a[i__ + (i__ - 1) * a_dim1] = -imeps * 2;
00461                 }
00462             }
00463 /* L240: */
00464         }
00465 
00466         i__1 = *n;
00467         for (i__ = 1; i__ <= i__1; ++i__) {
00468             e[i__ + i__ * e_dim1] = 1.;
00469             if (i__ <= 4) {
00470                 b[i__ + i__ * b_dim1] = -1.;
00471                 if (i__ > 2) {
00472                     b[i__ + i__ * b_dim1] = 1. - reeps;
00473                 }
00474                 if (i__ % 2 != 0 && i__ < *n) {
00475                     b[i__ + (i__ + 1) * b_dim1] = imeps;
00476                 } else if (i__ > 1) {
00477                     b[i__ + (i__ - 1) * b_dim1] = -imeps;
00478                 }
00479             } else if (i__ <= 8) {
00480                 if (i__ <= 6) {
00481                     b[i__ + i__ * b_dim1] = reeps;
00482                 } else {
00483                     b[i__ + i__ * b_dim1] = -reeps;
00484                 }
00485                 if (i__ % 2 != 0 && i__ < *n) {
00486                     b[i__ + (i__ + 1) * b_dim1] = imeps + 1.;
00487                 } else if (i__ > 1) {
00488                     b[i__ + (i__ - 1) * b_dim1] = -1. - imeps;
00489                 }
00490             } else {
00491                 b[i__ + i__ * b_dim1] = 1. - reeps;
00492                 if (i__ % 2 != 0 && i__ < *n) {
00493                     b[i__ + (i__ + 1) * b_dim1] = imeps * 2;
00494                 } else if (i__ > 1) {
00495                     b[i__ + (i__ - 1) * b_dim1] = -imeps * 2;
00496                 }
00497             }
00498 /* L250: */
00499         }
00500     }
00501 
00502 /*     Compute rhs (C, F) */
00503 
00504     dgemm_("N", "N", m, n, m, &c_b29, &a[a_offset], lda, &r__[r_offset], ldr, 
00505             &c_b30, &c__[c_offset], ldc);
00506     dgemm_("N", "N", m, n, n, &c_b33, &l[l_offset], ldl, &b[b_offset], ldb, &
00507             c_b29, &c__[c_offset], ldc);
00508     dgemm_("N", "N", m, n, m, &c_b29, &d__[d_offset], ldd, &r__[r_offset], 
00509             ldr, &c_b30, &f[f_offset], ldf);
00510     dgemm_("N", "N", m, n, n, &c_b33, &l[l_offset], ldl, &e[e_offset], lde, &
00511             c_b29, &f[f_offset], ldf);
00512 
00513 /*     End of DLATM5 */
00514 
00515     return 0;
00516 } /* dlatm5_ */