swiftnav: dlatm6.c Source File

Go to the documentation of this file.
00001 /* dlatm6.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__12 = 12;
00021 static integer c__8 = 8;
00022 static integer c__40 = 40;
00023 static integer c__2 = 2;
00024 static integer c__3 = 3;
00025 static integer c__60 = 60;
00026 
00027 /* Subroutine */ int dlatm6_(integer *type__, integer *n, doublereal *a, 
00028         integer *lda, doublereal *b, doublereal *x, integer *ldx, doublereal *
00029         y, integer *ldy, doublereal *alpha, doublereal *beta, doublereal *wx, 
00030         doublereal *wy, doublereal *s, doublereal *dif)
00031 {
00032     /* System generated locals */
00033     integer a_dim1, a_offset, b_dim1, b_offset, x_dim1, x_offset, y_dim1, 
00034             y_offset, i__1, i__2;
00035 
00036     /* Builtin functions */
00037     double sqrt(doublereal);
00038 
00039     /* Local variables */
00040     integer i__, j;
00041     doublereal z__[144] /* was [12][12] */;
00042     integer info;
00043     doublereal work[100];
00044     extern /* Subroutine */ int dlakf2_(integer *, integer *, doublereal *, 
00045             integer *, doublereal *, doublereal *, doublereal *, doublereal *, 
00046              integer *), dgesvd_(char *, char *, integer *, integer *, 
00047             doublereal *, integer *, doublereal *, doublereal *, integer *, 
00048             doublereal *, integer *, doublereal *, integer *, integer *), dlacpy_(char *, integer *, integer *, doublereal 
00049             *, integer *, doublereal *, 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 /*  DLATM6 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 
00075 /*        (A, B) = inverse(YH) * (Da, Db) * inverse(X) */
00076 
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 , and */
00085 
00086 /*  Type 2: */
00087 /*     Da =  1   -1    0    0    0    Db = 1   0   0   0   0 */
00088 /*           1    1    0    0    0         0   1   0   0   0 */
00089 /*           0    0    1    0    0         0   0   1   0   0 */
00090 /*           0    0    0   1+a  1+b        0   0   0   1   0 */
00091 /*           0    0    0  -1-b  1+a ,      0   0   0   0   1 . */
00092 
00093 /*  In both cases the same inverse(YH) and inverse(X) are used to compute */
00094 /*  (A, B), giving the exact eigenvectors to (A,B) as (YH, X): */
00095 
00096 /*  YH:  =  1    0   -y    y   -y    X =  1   0  -x  -x   x */
00097 /*          0    1   -y    y   -y         0   1   x  -x  -x */
00098 /*          0    0    1    0    0         0   0   1   0   0 */
00099 /*          0    0    0    1    0         0   0   0   1   0 */
00100 /*          0    0    0    0    1,        0   0   0   0   1 , */
00101 
00102 /* where a, b, x and y will have all values independently of each other. */
00103 
00104 /*  Arguments */
00105 /*  ========= */
00106 
00107 /*  TYPE    (input) INTEGER */
00108 /*          Specifies the problem type (see futher details). */
00109 
00110 /*  N       (input) INTEGER */
00111 /*          Size of the matrices A and B. */
00112 
00113 /*  A       (output) DOUBLE PRECISION array, dimension (LDA, N). */
00114 /*          On exit A N-by-N is initialized according to TYPE. */
00115 
00116 /*  LDA     (input) INTEGER */
00117 /*          The leading dimension of A and of B. */
00118 
00119 /*  B       (output) DOUBLE PRECISION array, dimension (LDA, N). */
00120 /*          On exit B N-by-N is initialized according to TYPE. */
00121 
00122 /*  X       (output) DOUBLE PRECISION array, dimension (LDX, N). */
00123 /*          On exit X is the N-by-N matrix of right eigenvectors. */
00124 
00125 /*  LDX     (input) INTEGER */
00126 /*          The leading dimension of X. */
00127 
00128 /*  Y       (output) DOUBLE PRECISION array, dimension (LDY, N). */
00129 /*          On exit Y is the N-by-N matrix of left eigenvectors. */
00130 
00131 /*  LDY     (input) INTEGER */
00132 /*          The leading dimension of Y. */
00133 
00134 /*  ALPHA   (input) DOUBLE PRECISION */
00135 /*  BETA    (input) DOUBLE PRECISION */
00136 /*          Weighting constants for matrix A. */
00137 
00138 /*  WX      (input) DOUBLE PRECISION */
00139 /*          Constant for right eigenvector matrix. */
00140 
00141 /*  WY      (input) DOUBLE PRECISION */
00142 /*          Constant for left eigenvector matrix. */
00143 
00144 /*  S       (output) DOUBLE PRECISION array, dimension (N) */
00145 /*          S(i) is the reciprocal condition number for eigenvalue i. */
00146 
00147 /*  DIF     (output) DOUBLE PRECISION array, dimension (N) */
00148 /*          DIF(i) is the reciprocal condition number for eigenvector i. */
00149 
00150 /*  ===================================================================== */
00151 
00152 /*     .. Parameters .. */
00153 /*     .. */
00154 /*     .. Local Scalars .. */
00155 /*     .. */
00156 /*     .. Local Arrays .. */
00157 /*     .. */
00158 /*     .. Intrinsic Functions .. */
00159 /*     .. */
00160 /*     .. External Subroutines .. */
00161 /*     .. */
00162 /*     .. Executable Statements .. */
00163 
00164 /*     Generate test problem ... */
00165 /*     (Da, Db) ... */
00166 
00167     /* Parameter adjustments */
00168     b_dim1 = *lda;
00169     b_offset = 1 + b_dim1;
00170     b -= b_offset;
00171     a_dim1 = *lda;
00172     a_offset = 1 + a_dim1;
00173     a -= a_offset;
00174     x_dim1 = *ldx;
00175     x_offset = 1 + x_dim1;
00176     x -= x_offset;
00177     y_dim1 = *ldy;
00178     y_offset = 1 + y_dim1;
00179     y -= y_offset;
00180     --s;
00181     --dif;
00182 
00183     /* Function Body */
00184     i__1 = *n;
00185     for (i__ = 1; i__ <= i__1; ++i__) {
00186         i__2 = *n;
00187         for (j = 1; j <= i__2; ++j) {
00188 
00189             if (i__ == j) {
00190                 a[i__ + i__ * a_dim1] = (doublereal) i__ + *alpha;
00191                 b[i__ + i__ * b_dim1] = 1.;
00192             } else {
00193                 a[i__ + j * a_dim1] = 0.;
00194                 b[i__ + j * b_dim1] = 0.;
00195             }
00196 
00197 /* L10: */
00198         }
00199 /* L20: */
00200     }
00201 
00202 /*     Form X and Y */
00203 
00204     dlacpy_("F", n, n, &b[b_offset], lda, &y[y_offset], ldy);
00205     y[y_dim1 + 3] = -(*wy);
00206     y[y_dim1 + 4] = *wy;
00207     y[y_dim1 + 5] = -(*wy);
00208     y[(y_dim1 << 1) + 3] = -(*wy);
00209     y[(y_dim1 << 1) + 4] = *wy;
00210     y[(y_dim1 << 1) + 5] = -(*wy);
00211 
00212     dlacpy_("F", n, n, &b[b_offset], lda, &x[x_offset], ldx);
00213     x[x_dim1 * 3 + 1] = -(*wx);
00214     x[(x_dim1 << 2) + 1] = -(*wx);
00215     x[x_dim1 * 5 + 1] = *wx;
00216     x[x_dim1 * 3 + 2] = *wx;
00217     x[(x_dim1 << 2) + 2] = -(*wx);
00218     x[x_dim1 * 5 + 2] = -(*wx);
00219 
00220 /*     Form (A, B) */
00221 
00222     b[b_dim1 * 3 + 1] = *wx + *wy;
00223     b[b_dim1 * 3 + 2] = -(*wx) + *wy;
00224     b[(b_dim1 << 2) + 1] = *wx - *wy;
00225     b[(b_dim1 << 2) + 2] = *wx - *wy;
00226     b[b_dim1 * 5 + 1] = -(*wx) + *wy;
00227     b[b_dim1 * 5 + 2] = *wx + *wy;
00228     if (*type__ == 1) {
00229         a[a_dim1 * 3 + 1] = *wx * a[a_dim1 + 1] + *wy * a[a_dim1 * 3 + 3];
00230         a[a_dim1 * 3 + 2] = -(*wx) * a[(a_dim1 << 1) + 2] + *wy * a[a_dim1 * 
00231                 3 + 3];
00232         a[(a_dim1 << 2) + 1] = *wx * a[a_dim1 + 1] - *wy * a[(a_dim1 << 2) + 
00233                 4];
00234         a[(a_dim1 << 2) + 2] = *wx * a[(a_dim1 << 1) + 2] - *wy * a[(a_dim1 <<
00235                  2) + 4];
00236         a[a_dim1 * 5 + 1] = -(*wx) * a[a_dim1 + 1] + *wy * a[a_dim1 * 5 + 5];
00237         a[a_dim1 * 5 + 2] = *wx * a[(a_dim1 << 1) + 2] + *wy * a[a_dim1 * 5 + 
00238                 5];
00239     } else if (*type__ == 2) {
00240         a[a_dim1 * 3 + 1] = *wx * 2. + *wy;
00241         a[a_dim1 * 3 + 2] = *wy;
00242         a[(a_dim1 << 2) + 1] = -(*wy) * (*alpha + 2. + *beta);
00243         a[(a_dim1 << 2) + 2] = *wx * 2. - *wy * (*alpha + 2. + *beta);
00244         a[a_dim1 * 5 + 1] = *wx * -2. + *wy * (*alpha - *beta);
00245         a[a_dim1 * 5 + 2] = *wy * (*alpha - *beta);
00246         a[a_dim1 + 1] = 1.;
00247         a[(a_dim1 << 1) + 1] = -1.;
00248         a[a_dim1 + 2] = 1.;
00249         a[(a_dim1 << 1) + 2] = a[a_dim1 + 1];
00250         a[a_dim1 * 3 + 3] = 1.;
00251         a[(a_dim1 << 2) + 4] = *alpha + 1.;
00252         a[a_dim1 * 5 + 4] = *beta + 1.;
00253         a[(a_dim1 << 2) + 5] = -a[a_dim1 * 5 + 4];
00254         a[a_dim1 * 5 + 5] = a[(a_dim1 << 2) + 4];
00255     }
00256 
00257 /*     Compute condition numbers */
00258 
00259     if (*type__ == 1) {
00260 
00261         s[1] = 1. / sqrt((*wy * 3. * *wy + 1.) / (a[a_dim1 + 1] * a[a_dim1 + 
00262                 1] + 1.));
00263         s[2] = 1. / sqrt((*wy * 3. * *wy + 1.) / (a[(a_dim1 << 1) + 2] * a[(
00264                 a_dim1 << 1) + 2] + 1.));
00265         s[3] = 1. / sqrt((*wx * 2. * *wx + 1.) / (a[a_dim1 * 3 + 3] * a[
00266                 a_dim1 * 3 + 3] + 1.));
00267         s[4] = 1. / sqrt((*wx * 2. * *wx + 1.) / (a[(a_dim1 << 2) + 4] * a[(
00268                 a_dim1 << 2) + 4] + 1.));
00269         s[5] = 1. / sqrt((*wx * 2. * *wx + 1.) / (a[a_dim1 * 5 + 5] * a[
00270                 a_dim1 * 5 + 5] + 1.));
00271 
00272         dlakf2_(&c__1, &c__4, &a[a_offset], lda, &a[(a_dim1 << 1) + 2], &b[
00273                 b_offset], &b[(b_dim1 << 1) + 2], z__, &c__12);
00274         dgesvd_("N", "N", &c__8, &c__8, z__, &c__12, work, &work[8], &c__1, &
00275                 work[9], &c__1, &work[10], &c__40, &info);
00276         dif[1] = work[7];
00277 
00278         dlakf2_(&c__4, &c__1, &a[a_offset], lda, &a[a_dim1 * 5 + 5], &b[
00279                 b_offset], &b[b_dim1 * 5 + 5], z__, &c__12);
00280         dgesvd_("N", "N", &c__8, &c__8, z__, &c__12, work, &work[8], &c__1, &
00281                 work[9], &c__1, &work[10], &c__40, &info);
00282         dif[5] = work[7];
00283 
00284     } else if (*type__ == 2) {
00285 
00286         s[1] = 1. / sqrt(*wy * *wy + .33333333333333331);
00287         s[2] = s[1];
00288         s[3] = 1. / sqrt(*wx * *wx + .5);
00289         s[4] = 1. / sqrt((*wx * 2. * *wx + 1.) / ((*alpha + 1.) * (*alpha + 
00290                 1.) + 1. + (*beta + 1.) * (*beta + 1.)));
00291         s[5] = s[4];
00292 
00293         dlakf2_(&c__2, &c__3, &a[a_offset], lda, &a[a_dim1 * 3 + 3], &b[
00294                 b_offset], &b[b_dim1 * 3 + 3], z__, &c__12);
00295         dgesvd_("N", "N", &c__12, &c__12, z__, &c__12, work, &work[12], &c__1, 
00296                  &work[13], &c__1, &work[14], &c__60, &info);
00297         dif[1] = work[11];
00298 
00299         dlakf2_(&c__3, &c__2, &a[a_offset], lda, &a[(a_dim1 << 2) + 4], &b[
00300                 b_offset], &b[(b_dim1 << 2) + 4], z__, &c__12);
00301         dgesvd_("N", "N", &c__12, &c__12, z__, &c__12, work, &work[12], &c__1, 
00302                  &work[13], &c__1, &work[14], &c__60, &info);
00303         dif[5] = work[11];
00304 
00305     }
00306 
00307     return 0;
00308 
00309 /*     End of DLATM6 */
00310 
00311 } /* dlatm6_ */