zlahilb.c
Go to the documentation of this file.
00001 /* zlahilb.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__2 = 2;
00019 static doublecomplex c_b6 = {0.,0.};
00020 
00021 /* Subroutine */ int zlahilb_(integer *n, integer *nrhs, doublecomplex *a, 
00022         integer *lda, doublecomplex *x, integer *ldx, doublecomplex *b, 
00023         integer *ldb, doublereal *work, integer *info, char *path)
00024 {
00025     /* Initialized data */
00026 
00027     static doublecomplex d1[8] = { {-1.,0.},{0.,1.},{-1.,-1.},{0.,-1.},{1.,0.}
00028             ,{-1.,1.},{1.,1.},{1.,-1.} };
00029     static doublecomplex d2[8] = { {-1.,0.},{0.,-1.},{-1.,1.},{0.,1.},{1.,0.},
00030             {-1.,-1.},{1.,-1.},{1.,1.} };
00031     static doublecomplex invd1[8] = { {-1.,0.},{0.,-1.},{-.5,.5},{0.,1.},{1.,
00032             0.},{-.5,-.5},{.5,-.5},{.5,.5} };
00033     static doublecomplex invd2[8] = { {-1.,0.},{0.,1.},{-.5,-.5},{0.,-1.},{1.,
00034             0.},{-.5,.5},{.5,.5},{.5,-.5} };
00035 
00036     /* System generated locals */
00037     integer a_dim1, a_offset, x_dim1, x_offset, b_dim1, b_offset, i__1, i__2, 
00038             i__3, i__4, i__5;
00039     doublereal d__1;
00040     doublecomplex z__1, z__2;
00041 
00042     /* Builtin functions */
00043     /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);
00044 
00045     /* Local variables */
00046     integer i__, j, m, r__;
00047     char c2[2];
00048     integer ti, tm;
00049     doublecomplex tmp;
00050     extern /* Subroutine */ int xerbla_(char *, integer *);
00051     extern logical lsamen_(integer *, char *, char *);
00052     extern /* Subroutine */ int zlaset_(char *, integer *, integer *, 
00053             doublecomplex *, doublecomplex *, doublecomplex *, integer *);
00054 
00055 
00056 /*  -- LAPACK auxiliary test routine (version 3.0) -- */
00057 /*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
00058 /*     Courant Institute, Argonne National Lab, and Rice University */
00059 /*     28 August, 2006 */
00060 
00061 /*     David Vu <dtv@cs.berkeley.edu> */
00062 /*     Yozo Hida <yozo@cs.berkeley.edu> */
00063 /*     Jason Riedy <ejr@cs.berkeley.edu> */
00064 /*     D. Halligan <dhalligan@berkeley.edu> */
00065 
00066 /*     .. Scalar Arguments .. */
00067 /*     .. Array Arguments .. */
00068 /*     .. */
00069 
00070 /*  Purpose */
00071 /*  ======= */
00072 
00073 /*  ZLAHILB generates an N by N scaled Hilbert matrix in A along with */
00074 /*  NRHS right-hand sides in B and solutions in X such that A*X=B. */
00075 
00076 /*  The Hilbert matrix is scaled by M = LCM(1, 2, ..., 2*N-1) so that all */
00077 /*  entries are integers.  The right-hand sides are the first NRHS */
00078 /*  columns of M * the identity matrix, and the solutions are the */
00079 /*  first NRHS columns of the inverse Hilbert matrix. */
00080 
00081 /*  The condition number of the Hilbert matrix grows exponentially with */
00082 /*  its size, roughly as O(e ** (3.5*N)).  Additionally, the inverse */
00083 /*  Hilbert matrices beyond a relatively small dimension cannot be */
00084 /*  generated exactly without extra precision.  Precision is exhausted */
00085 /*  when the largest entry in the inverse Hilbert matrix is greater than */
00086 /*  2 to the power of the number of bits in the fraction of the data type */
00087 /*  used plus one, which is 24 for single precision. */
00088 
00089 /*  In single, the generated solution is exact for N <= 6 and has */
00090 /*  small componentwise error for 7 <= N <= 11. */
00091 
00092 /*  Arguments */
00093 /*  ========= */
00094 
00095 /*  N       (input) INTEGER */
00096 /*          The dimension of the matrix A. */
00097 
00098 /*  NRHS    (input) NRHS */
00099 /*          The requested number of right-hand sides. */
00100 
00101 /*  A       (output) COMPLEX array, dimension (LDA, N) */
00102 /*          The generated scaled Hilbert matrix. */
00103 
00104 /*  LDA     (input) INTEGER */
00105 /*          The leading dimension of the array A.  LDA >= N. */
00106 
00107 /*  X       (output) COMPLEX array, dimension (LDX, NRHS) */
00108 /*          The generated exact solutions.  Currently, the first NRHS */
00109 /*          columns of the inverse Hilbert matrix. */
00110 
00111 /*  LDX     (input) INTEGER */
00112 /*          The leading dimension of the array X.  LDX >= N. */
00113 
00114 /*  B       (output) REAL array, dimension (LDB, NRHS) */
00115 /*          The generated right-hand sides.  Currently, the first NRHS */
00116 /*          columns of LCM(1, 2, ..., 2*N-1) * the identity matrix. */
00117 
00118 /*  LDB     (input) INTEGER */
00119 /*          The leading dimension of the array B.  LDB >= N. */
00120 
00121 /*  WORK    (workspace) REAL array, dimension (N) */
00122 
00123 
00124 /*  INFO    (output) INTEGER */
00125 /*          = 0: successful exit */
00126 /*          = 1: N is too large; the data is still generated but may not */
00127 /*               be not exact. */
00128 /*          < 0: if INFO = -i, the i-th argument had an illegal value */
00129 
00130 /*  ===================================================================== */
00131 /*     .. Local Scalars .. */
00132 /*     .. Parameters .. */
00133 /*     NMAX_EXACT   the largest dimension where the generated data is */
00134 /*                  exact. */
00135 /*     NMAX_APPROX  the largest dimension where the generated data has */
00136 /*                  a small componentwise relative error. */
00137 /*     ??? complex uses how many bits ??? */
00138 /*     d's are generated from random permuation of those eight elements. */
00139     /* Parameter adjustments */
00140     --work;
00141     a_dim1 = *lda;
00142     a_offset = 1 + a_dim1;
00143     a -= a_offset;
00144     x_dim1 = *ldx;
00145     x_offset = 1 + x_dim1;
00146     x -= x_offset;
00147     b_dim1 = *ldb;
00148     b_offset = 1 + b_dim1;
00149     b -= b_offset;
00150 
00151     /* Function Body */
00152 /*     .. */
00153 /*     .. External Functions */
00154 /*     .. */
00155 /*     .. Executable Statements .. */
00156     s_copy(c2, path + 1, (ftnlen)2, (ftnlen)2);
00157 
00158 /*     Test the input arguments */
00159 
00160     *info = 0;
00161     if (*n < 0 || *n > 11) {
00162         *info = -1;
00163     } else if (*nrhs < 0) {
00164         *info = -2;
00165     } else if (*lda < *n) {
00166         *info = -4;
00167     } else if (*ldx < *n) {
00168         *info = -6;
00169     } else if (*ldb < *n) {
00170         *info = -8;
00171     }
00172     if (*info < 0) {
00173         i__1 = -(*info);
00174         xerbla_("ZLAHILB", &i__1);
00175         return 0;
00176     }
00177     if (*n > 6) {
00178         *info = 1;
00179     }
00180 /*     Compute M = the LCM of the integers [1, 2*N-1].  The largest */
00181 /*     reasonable N is small enough that integers suffice (up to N = 11). */
00182     m = 1;
00183     i__1 = (*n << 1) - 1;
00184     for (i__ = 2; i__ <= i__1; ++i__) {
00185         tm = m;
00186         ti = i__;
00187         r__ = tm % ti;
00188         while(r__ != 0) {
00189             tm = ti;
00190             ti = r__;
00191             r__ = tm % ti;
00192         }
00193         m = m / ti * i__;
00194     }
00195 /*     Generate the scaled Hilbert matrix in A */
00196 /*     If we are testing SY routines, take D1_i = D2_i, else, D1_i = D2_i* */
00197     if (lsamen_(&c__2, c2, "SY")) {
00198         i__1 = *n;
00199         for (j = 1; j <= i__1; ++j) {
00200             i__2 = *n;
00201             for (i__ = 1; i__ <= i__2; ++i__) {
00202                 i__3 = i__ + j * a_dim1;
00203                 i__4 = j % 8;
00204                 d__1 = (doublereal) m / (i__ + j - 1);
00205                 z__2.r = d__1 * d1[i__4].r, z__2.i = d__1 * d1[i__4].i;
00206                 i__5 = i__ % 8;
00207                 z__1.r = z__2.r * d1[i__5].r - z__2.i * d1[i__5].i, z__1.i = 
00208                         z__2.r * d1[i__5].i + z__2.i * d1[i__5].r;
00209                 a[i__3].r = z__1.r, a[i__3].i = z__1.i;
00210             }
00211         }
00212     } else {
00213         i__1 = *n;
00214         for (j = 1; j <= i__1; ++j) {
00215             i__2 = *n;
00216             for (i__ = 1; i__ <= i__2; ++i__) {
00217                 i__3 = i__ + j * a_dim1;
00218                 i__4 = j % 8;
00219                 d__1 = (doublereal) m / (i__ + j - 1);
00220                 z__2.r = d__1 * d1[i__4].r, z__2.i = d__1 * d1[i__4].i;
00221                 i__5 = i__ % 8;
00222                 z__1.r = z__2.r * d2[i__5].r - z__2.i * d2[i__5].i, z__1.i = 
00223                         z__2.r * d2[i__5].i + z__2.i * d2[i__5].r;
00224                 a[i__3].r = z__1.r, a[i__3].i = z__1.i;
00225             }
00226         }
00227     }
00228 /*     Generate matrix B as simply the first NRHS columns of M * the */
00229 /*     identity. */
00230     d__1 = (doublereal) m;
00231     tmp.r = d__1, tmp.i = 0.;
00232     zlaset_("Full", n, nrhs, &c_b6, &tmp, &b[b_offset], ldb);
00233 /*     Generate the true solutions in X.  Because B = the first NRHS */
00234 /*     columns of M*I, the true solutions are just the first NRHS columns */
00235 /*     of the inverse Hilbert matrix. */
00236     work[1] = (doublereal) (*n);
00237     i__1 = *n;
00238     for (j = 2; j <= i__1; ++j) {
00239         work[j] = work[j - 1] / (j - 1) * (j - 1 - *n) / (j - 1) * (*n + j - 
00240                 1);
00241     }
00242 /*     If we are testing SY routines, take D1_i = D2_i, else, D1_i = D2_i* */
00243     if (lsamen_(&c__2, c2, "SY")) {
00244         i__1 = *nrhs;
00245         for (j = 1; j <= i__1; ++j) {
00246             i__2 = *n;
00247             for (i__ = 1; i__ <= i__2; ++i__) {
00248                 i__3 = i__ + j * x_dim1;
00249                 i__4 = j % 8;
00250                 d__1 = work[i__] * work[j] / (i__ + j - 1);
00251                 z__2.r = d__1 * invd1[i__4].r, z__2.i = d__1 * invd1[i__4].i;
00252                 i__5 = i__ % 8;
00253                 z__1.r = z__2.r * invd1[i__5].r - z__2.i * invd1[i__5].i, 
00254                         z__1.i = z__2.r * invd1[i__5].i + z__2.i * invd1[i__5]
00255                         .r;
00256                 x[i__3].r = z__1.r, x[i__3].i = z__1.i;
00257             }
00258         }
00259     } else {
00260         i__1 = *nrhs;
00261         for (j = 1; j <= i__1; ++j) {
00262             i__2 = *n;
00263             for (i__ = 1; i__ <= i__2; ++i__) {
00264                 i__3 = i__ + j * x_dim1;
00265                 i__4 = j % 8;
00266                 d__1 = work[i__] * work[j] / (i__ + j - 1);
00267                 z__2.r = d__1 * invd2[i__4].r, z__2.i = d__1 * invd2[i__4].i;
00268                 i__5 = i__ % 8;
00269                 z__1.r = z__2.r * invd1[i__5].r - z__2.i * invd1[i__5].i, 
00270                         z__1.i = z__2.r * invd1[i__5].i + z__2.i * invd1[i__5]
00271                         .r;
00272                 x[i__3].r = z__1.r, x[i__3].i = z__1.i;
00273             }
00274         }
00275     }
00276     return 0;
00277 } /* zlahilb_ */


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