swiftnav: slatm5.c Source File

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