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

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