clatm6.c
Go to the documentation of this file.
00001 /* clatm6.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 integer c__1 = 1;
00019 static integer c__4 = 4;
00020 static integer c__8 = 8;
00021 static integer c__24 = 24;
00022 
00023 /* Subroutine */ int clatm6_(integer *type__, integer *n, complex *a, integer 
00024         *lda, complex *b, complex *x, integer *ldx, complex *y, integer *ldy, 
00025         complex *alpha, complex *beta, complex *wx, complex *wy, real *s, 
00026         real *dif)
00027 {
00028     /* System generated locals */
00029     integer a_dim1, a_offset, b_dim1, b_offset, x_dim1, x_offset, y_dim1, 
00030             y_offset, i__1, i__2, i__3;
00031     real r__1, r__2;
00032     complex q__1, q__2, q__3, q__4;
00033 
00034     /* Builtin functions */
00035     void r_cnjg(complex *, complex *);
00036     double c_abs(complex *), sqrt(doublereal);
00037 
00038     /* Local variables */
00039     integer i__, j;
00040     complex z__[64]     /* was [8][8] */;
00041     integer info;
00042     complex work[26];
00043     extern /* Subroutine */ int clakf2_(integer *, integer *, complex *, 
00044             integer *, complex *, complex *, complex *, complex *, integer *);
00045     real rwork[50];
00046     extern /* Subroutine */ int cgesvd_(char *, char *, integer *, integer *, 
00047             complex *, integer *, real *, complex *, integer *, complex *, 
00048             integer *, complex *, integer *, real *, integer *), clacpy_(char *, integer *, integer *, complex *, integer 
00049             *, complex *, integer *);
00050 
00051 
00052 /*  -- LAPACK test routine (version 3.1) -- */
00053 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00054 /*     November 2006 */
00055 
00056 /*     .. Scalar Arguments .. */
00057 /*     .. */
00058 /*     .. Array Arguments .. */
00059 /*     .. */
00060 
00061 /*  Purpose */
00062 /*  ======= */
00063 
00064 /*  CLATM6 generates test matrices for the generalized eigenvalue */
00065 /*  problem, their corresponding right and left eigenvector matrices, */
00066 /*  and also reciprocal condition numbers for all eigenvalues and */
00067 /*  the reciprocal condition numbers of eigenvectors corresponding to */
00068 /*  the 1th and 5th eigenvalues. */
00069 
00070 /*  Test Matrices */
00071 /*  ============= */
00072 
00073 /*  Two kinds of test matrix pairs */
00074 /*           (A, B) = inverse(YH) * (Da, Db) * inverse(X) */
00075 /*  are used in the tests: */
00076 
00077 /*  Type 1: */
00078 /*     Da = 1+a   0    0    0    0    Db = 1   0   0   0   0 */
00079 /*           0   2+a   0    0    0         0   1   0   0   0 */
00080 /*           0    0   3+a   0    0         0   0   1   0   0 */
00081 /*           0    0    0   4+a   0         0   0   0   1   0 */
00082 /*           0    0    0    0   5+a ,      0   0   0   0   1 */
00083 /*  and Type 2: */
00084 /*     Da = 1+i   0    0       0       0    Db = 1   0   0   0   0 */
00085 /*           0   1-i   0       0       0         0   1   0   0   0 */
00086 /*           0    0    1       0       0         0   0   1   0   0 */
00087 /*           0    0    0 (1+a)+(1+b)i  0         0   0   0   1   0 */
00088 /*           0    0    0       0 (1+a)-(1+b)i,   0   0   0   0   1 . */
00089 
00090 /*  In both cases the same inverse(YH) and inverse(X) are used to compute */
00091 /*  (A, B), giving the exact eigenvectors to (A,B) as (YH, X): */
00092 
00093 /*  YH:  =  1    0   -y    y   -y    X =  1   0  -x  -x   x */
00094 /*          0    1   -y    y   -y         0   1   x  -x  -x */
00095 /*          0    0    1    0    0         0   0   1   0   0 */
00096 /*          0    0    0    1    0         0   0   0   1   0 */
00097 /*          0    0    0    0    1,        0   0   0   0   1 , where */
00098 
00099 /*  a, b, x and y will have all values independently of each other. */
00100 
00101 /*  Arguments */
00102 /*  ========= */
00103 
00104 /*  TYPE    (input) INTEGER */
00105 /*          Specifies the problem type (see futher details). */
00106 
00107 /*  N       (input) INTEGER */
00108 /*          Size of the matrices A and B. */
00109 
00110 /*  A       (output) COMPLEX array, dimension (LDA, N). */
00111 /*          On exit A N-by-N is initialized according to TYPE. */
00112 
00113 /*  LDA     (input) INTEGER */
00114 /*          The leading dimension of A and of B. */
00115 
00116 /*  B       (output) COMPLEX array, dimension (LDA, N). */
00117 /*          On exit B N-by-N is initialized according to TYPE. */
00118 
00119 /*  X       (output) COMPLEX array, dimension (LDX, N). */
00120 /*          On exit X is the N-by-N matrix of right eigenvectors. */
00121 
00122 /*  LDX     (input) INTEGER */
00123 /*          The leading dimension of X. */
00124 
00125 /*  Y       (output) COMPLEX array, dimension (LDY, N). */
00126 /*          On exit Y is the N-by-N matrix of left eigenvectors. */
00127 
00128 /*  LDY     (input) INTEGER */
00129 /*          The leading dimension of Y. */
00130 
00131 /*  ALPHA   (input) COMPLEX */
00132 /*  BETA    (input) COMPLEX */
00133 /*          Weighting constants for matrix A. */
00134 
00135 /*  WX      (input) COMPLEX */
00136 /*          Constant for right eigenvector matrix. */
00137 
00138 /*  WY      (input) COMPLEX */
00139 /*          Constant for left eigenvector matrix. */
00140 
00141 /*  S       (output) REAL array, dimension (N) */
00142 /*          S(i) is the reciprocal condition number for eigenvalue i. */
00143 
00144 /*  DIF     (output) REAL array, dimension (N) */
00145 /*          DIF(i) is the reciprocal condition number for eigenvector i. */
00146 
00147 /*  ===================================================================== */
00148 
00149 /*     .. Parameters .. */
00150 /*     .. */
00151 /*     .. Local Scalars .. */
00152 /*     .. */
00153 /*     .. Local Arrays .. */
00154 /*     .. */
00155 /*     .. Intrinsic Functions .. */
00156 /*     .. */
00157 /*     .. External Subroutines .. */
00158 /*     .. */
00159 /*     .. Executable Statements .. */
00160 
00161 /*     Generate test problem ... */
00162 /*     (Da, Db) ... */
00163 
00164     /* Parameter adjustments */
00165     b_dim1 = *lda;
00166     b_offset = 1 + b_dim1;
00167     b -= b_offset;
00168     a_dim1 = *lda;
00169     a_offset = 1 + a_dim1;
00170     a -= a_offset;
00171     x_dim1 = *ldx;
00172     x_offset = 1 + x_dim1;
00173     x -= x_offset;
00174     y_dim1 = *ldy;
00175     y_offset = 1 + y_dim1;
00176     y -= y_offset;
00177     --s;
00178     --dif;
00179 
00180     /* Function Body */
00181     i__1 = *n;
00182     for (i__ = 1; i__ <= i__1; ++i__) {
00183         i__2 = *n;
00184         for (j = 1; j <= i__2; ++j) {
00185 
00186             if (i__ == j) {
00187                 i__3 = i__ + i__ * a_dim1;
00188                 q__2.r = (real) i__, q__2.i = 0.f;
00189                 q__1.r = q__2.r + alpha->r, q__1.i = q__2.i + alpha->i;
00190                 a[i__3].r = q__1.r, a[i__3].i = q__1.i;
00191                 i__3 = i__ + i__ * b_dim1;
00192                 b[i__3].r = 1.f, b[i__3].i = 0.f;
00193             } else {
00194                 i__3 = i__ + j * a_dim1;
00195                 a[i__3].r = 0.f, a[i__3].i = 0.f;
00196                 i__3 = i__ + j * b_dim1;
00197                 b[i__3].r = 0.f, b[i__3].i = 0.f;
00198             }
00199 
00200 /* L10: */
00201         }
00202 /* L20: */
00203     }
00204     if (*type__ == 2) {
00205         i__1 = a_dim1 + 1;
00206         a[i__1].r = 1.f, a[i__1].i = 1.f;
00207         i__1 = (a_dim1 << 1) + 2;
00208         r_cnjg(&q__1, &a[a_dim1 + 1]);
00209         a[i__1].r = q__1.r, a[i__1].i = q__1.i;
00210         i__1 = a_dim1 * 3 + 3;
00211         a[i__1].r = 1.f, a[i__1].i = 0.f;
00212         i__1 = (a_dim1 << 2) + 4;
00213         q__2.r = alpha->r + 1.f, q__2.i = alpha->i + 0.f;
00214         r__1 = q__2.r;
00215         q__3.r = beta->r + 1.f, q__3.i = beta->i + 0.f;
00216         r__2 = q__3.r;
00217         q__1.r = r__1, q__1.i = r__2;
00218         a[i__1].r = q__1.r, a[i__1].i = q__1.i;
00219         i__1 = a_dim1 * 5 + 5;
00220         r_cnjg(&q__1, &a[(a_dim1 << 2) + 4]);
00221         a[i__1].r = q__1.r, a[i__1].i = q__1.i;
00222     }
00223 
00224 /*     Form X and Y */
00225 
00226     clacpy_("F", n, n, &b[b_offset], lda, &y[y_offset], ldy);
00227     i__1 = y_dim1 + 3;
00228     r_cnjg(&q__2, wy);
00229     q__1.r = -q__2.r, q__1.i = -q__2.i;
00230     y[i__1].r = q__1.r, y[i__1].i = q__1.i;
00231     i__1 = y_dim1 + 4;
00232     r_cnjg(&q__1, wy);
00233     y[i__1].r = q__1.r, y[i__1].i = q__1.i;
00234     i__1 = y_dim1 + 5;
00235     r_cnjg(&q__2, wy);
00236     q__1.r = -q__2.r, q__1.i = -q__2.i;
00237     y[i__1].r = q__1.r, y[i__1].i = q__1.i;
00238     i__1 = (y_dim1 << 1) + 3;
00239     r_cnjg(&q__2, wy);
00240     q__1.r = -q__2.r, q__1.i = -q__2.i;
00241     y[i__1].r = q__1.r, y[i__1].i = q__1.i;
00242     i__1 = (y_dim1 << 1) + 4;
00243     r_cnjg(&q__1, wy);
00244     y[i__1].r = q__1.r, y[i__1].i = q__1.i;
00245     i__1 = (y_dim1 << 1) + 5;
00246     r_cnjg(&q__2, wy);
00247     q__1.r = -q__2.r, q__1.i = -q__2.i;
00248     y[i__1].r = q__1.r, y[i__1].i = q__1.i;
00249 
00250     clacpy_("F", n, n, &b[b_offset], lda, &x[x_offset], ldx);
00251     i__1 = x_dim1 * 3 + 1;
00252     q__1.r = -wx->r, q__1.i = -wx->i;
00253     x[i__1].r = q__1.r, x[i__1].i = q__1.i;
00254     i__1 = (x_dim1 << 2) + 1;
00255     q__1.r = -wx->r, q__1.i = -wx->i;
00256     x[i__1].r = q__1.r, x[i__1].i = q__1.i;
00257     i__1 = x_dim1 * 5 + 1;
00258     x[i__1].r = wx->r, x[i__1].i = wx->i;
00259     i__1 = x_dim1 * 3 + 2;
00260     x[i__1].r = wx->r, x[i__1].i = wx->i;
00261     i__1 = (x_dim1 << 2) + 2;
00262     q__1.r = -wx->r, q__1.i = -wx->i;
00263     x[i__1].r = q__1.r, x[i__1].i = q__1.i;
00264     i__1 = x_dim1 * 5 + 2;
00265     q__1.r = -wx->r, q__1.i = -wx->i;
00266     x[i__1].r = q__1.r, x[i__1].i = q__1.i;
00267 
00268 /*     Form (A, B) */
00269 
00270     i__1 = b_dim1 * 3 + 1;
00271     q__1.r = wx->r + wy->r, q__1.i = wx->i + wy->i;
00272     b[i__1].r = q__1.r, b[i__1].i = q__1.i;
00273     i__1 = b_dim1 * 3 + 2;
00274     q__2.r = -wx->r, q__2.i = -wx->i;
00275     q__1.r = q__2.r + wy->r, q__1.i = q__2.i + wy->i;
00276     b[i__1].r = q__1.r, b[i__1].i = q__1.i;
00277     i__1 = (b_dim1 << 2) + 1;
00278     q__1.r = wx->r - wy->r, q__1.i = wx->i - wy->i;
00279     b[i__1].r = q__1.r, b[i__1].i = q__1.i;
00280     i__1 = (b_dim1 << 2) + 2;
00281     q__1.r = wx->r - wy->r, q__1.i = wx->i - wy->i;
00282     b[i__1].r = q__1.r, b[i__1].i = q__1.i;
00283     i__1 = b_dim1 * 5 + 1;
00284     q__2.r = -wx->r, q__2.i = -wx->i;
00285     q__1.r = q__2.r + wy->r, q__1.i = q__2.i + wy->i;
00286     b[i__1].r = q__1.r, b[i__1].i = q__1.i;
00287     i__1 = b_dim1 * 5 + 2;
00288     q__1.r = wx->r + wy->r, q__1.i = wx->i + wy->i;
00289     b[i__1].r = q__1.r, b[i__1].i = q__1.i;
00290     i__1 = a_dim1 * 3 + 1;
00291     i__2 = a_dim1 + 1;
00292     q__2.r = wx->r * a[i__2].r - wx->i * a[i__2].i, q__2.i = wx->r * a[i__2]
00293             .i + wx->i * a[i__2].r;
00294     i__3 = a_dim1 * 3 + 3;
00295     q__3.r = wy->r * a[i__3].r - wy->i * a[i__3].i, q__3.i = wy->r * a[i__3]
00296             .i + wy->i * a[i__3].r;
00297     q__1.r = q__2.r + q__3.r, q__1.i = q__2.i + q__3.i;
00298     a[i__1].r = q__1.r, a[i__1].i = q__1.i;
00299     i__1 = a_dim1 * 3 + 2;
00300     q__3.r = -wx->r, q__3.i = -wx->i;
00301     i__2 = (a_dim1 << 1) + 2;
00302     q__2.r = q__3.r * a[i__2].r - q__3.i * a[i__2].i, q__2.i = q__3.r * a[
00303             i__2].i + q__3.i * a[i__2].r;
00304     i__3 = a_dim1 * 3 + 3;
00305     q__4.r = wy->r * a[i__3].r - wy->i * a[i__3].i, q__4.i = wy->r * a[i__3]
00306             .i + wy->i * a[i__3].r;
00307     q__1.r = q__2.r + q__4.r, q__1.i = q__2.i + q__4.i;
00308     a[i__1].r = q__1.r, a[i__1].i = q__1.i;
00309     i__1 = (a_dim1 << 2) + 1;
00310     i__2 = a_dim1 + 1;
00311     q__2.r = wx->r * a[i__2].r - wx->i * a[i__2].i, q__2.i = wx->r * a[i__2]
00312             .i + wx->i * a[i__2].r;
00313     i__3 = (a_dim1 << 2) + 4;
00314     q__3.r = wy->r * a[i__3].r - wy->i * a[i__3].i, q__3.i = wy->r * a[i__3]
00315             .i + wy->i * a[i__3].r;
00316     q__1.r = q__2.r - q__3.r, q__1.i = q__2.i - q__3.i;
00317     a[i__1].r = q__1.r, a[i__1].i = q__1.i;
00318     i__1 = (a_dim1 << 2) + 2;
00319     i__2 = (a_dim1 << 1) + 2;
00320     q__2.r = wx->r * a[i__2].r - wx->i * a[i__2].i, q__2.i = wx->r * a[i__2]
00321             .i + wx->i * a[i__2].r;
00322     i__3 = (a_dim1 << 2) + 4;
00323     q__3.r = wy->r * a[i__3].r - wy->i * a[i__3].i, q__3.i = wy->r * a[i__3]
00324             .i + wy->i * a[i__3].r;
00325     q__1.r = q__2.r - q__3.r, q__1.i = q__2.i - q__3.i;
00326     a[i__1].r = q__1.r, a[i__1].i = q__1.i;
00327     i__1 = a_dim1 * 5 + 1;
00328     q__3.r = -wx->r, q__3.i = -wx->i;
00329     i__2 = a_dim1 + 1;
00330     q__2.r = q__3.r * a[i__2].r - q__3.i * a[i__2].i, q__2.i = q__3.r * a[
00331             i__2].i + q__3.i * a[i__2].r;
00332     i__3 = a_dim1 * 5 + 5;
00333     q__4.r = wy->r * a[i__3].r - wy->i * a[i__3].i, q__4.i = wy->r * a[i__3]
00334             .i + wy->i * a[i__3].r;
00335     q__1.r = q__2.r + q__4.r, q__1.i = q__2.i + q__4.i;
00336     a[i__1].r = q__1.r, a[i__1].i = q__1.i;
00337     i__1 = a_dim1 * 5 + 2;
00338     i__2 = (a_dim1 << 1) + 2;
00339     q__2.r = wx->r * a[i__2].r - wx->i * a[i__2].i, q__2.i = wx->r * a[i__2]
00340             .i + wx->i * a[i__2].r;
00341     i__3 = a_dim1 * 5 + 5;
00342     q__3.r = wy->r * a[i__3].r - wy->i * a[i__3].i, q__3.i = wy->r * a[i__3]
00343             .i + wy->i * a[i__3].r;
00344     q__1.r = q__2.r + q__3.r, q__1.i = q__2.i + q__3.i;
00345     a[i__1].r = q__1.r, a[i__1].i = q__1.i;
00346 
00347 /*     Compute condition numbers */
00348 
00349     s[1] = 1.f / sqrt((c_abs(wy) * 3.f * c_abs(wy) + 1.f) / (c_abs(&a[a_dim1 
00350             + 1]) * c_abs(&a[a_dim1 + 1]) + 1.f));
00351     s[2] = 1.f / sqrt((c_abs(wy) * 3.f * c_abs(wy) + 1.f) / (c_abs(&a[(a_dim1 
00352             << 1) + 2]) * c_abs(&a[(a_dim1 << 1) + 2]) + 1.f));
00353     s[3] = 1.f / sqrt((c_abs(wx) * 2.f * c_abs(wx) + 1.f) / (c_abs(&a[a_dim1 *
00354              3 + 3]) * c_abs(&a[a_dim1 * 3 + 3]) + 1.f));
00355     s[4] = 1.f / sqrt((c_abs(wx) * 2.f * c_abs(wx) + 1.f) / (c_abs(&a[(a_dim1 
00356             << 2) + 4]) * c_abs(&a[(a_dim1 << 2) + 4]) + 1.f));
00357     s[5] = 1.f / sqrt((c_abs(wx) * 2.f * c_abs(wx) + 1.f) / (c_abs(&a[a_dim1 *
00358              5 + 5]) * c_abs(&a[a_dim1 * 5 + 5]) + 1.f));
00359 
00360     clakf2_(&c__1, &c__4, &a[a_offset], lda, &a[(a_dim1 << 1) + 2], &b[
00361             b_offset], &b[(b_dim1 << 1) + 2], z__, &c__8);
00362     cgesvd_("N", "N", &c__8, &c__8, z__, &c__8, rwork, work, &c__1, &work[1], 
00363             &c__1, &work[2], &c__24, &rwork[8], &info);
00364     dif[1] = rwork[7];
00365 
00366     clakf2_(&c__4, &c__1, &a[a_offset], lda, &a[a_dim1 * 5 + 5], &b[b_offset], 
00367              &b[b_dim1 * 5 + 5], z__, &c__8);
00368     cgesvd_("N", "N", &c__8, &c__8, z__, &c__8, rwork, work, &c__1, &work[1], 
00369             &c__1, &work[2], &c__24, &rwork[8], &info);
00370     dif[5] = rwork[7];
00371 
00372     return 0;
00373 
00374 /*     End of CLATM6 */
00375 
00376 } /* clatm6_ */

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