zlattr.c
Go to the documentation of this file.
00001 /* zlattr.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__5 = 5;
00019 static integer c__2 = 2;
00020 static integer c__1 = 1;
00021 static integer c__4 = 4;
00022 static doublereal c_b92 = 2.;
00023 static integer c_n1 = -1;
00024 
00025 /* Subroutine */ int zlattr_(integer *imat, char *uplo, char *trans, char *
00026         diag, integer *iseed, integer *n, doublecomplex *a, integer *lda, 
00027         doublecomplex *b, doublecomplex *work, doublereal *rwork, integer *
00028         info)
00029 {
00030     /* System generated locals */
00031     integer a_dim1, a_offset, i__1, i__2, i__3, i__4;
00032     doublereal d__1, d__2;
00033     doublecomplex z__1, z__2;
00034 
00035     /* Builtin functions */
00036     /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);
00037     void z_div(doublecomplex *, doublecomplex *, doublecomplex *);
00038     double pow_dd(doublereal *, doublereal *), sqrt(doublereal);
00039     void d_cnjg(doublecomplex *, doublecomplex *);
00040     double z_abs(doublecomplex *);
00041 
00042     /* Local variables */
00043     doublereal c__;
00044     integer i__, j;
00045     doublecomplex s;
00046     doublereal x, y, z__;
00047     doublecomplex ra, rb;
00048     integer kl, ku, iy;
00049     doublereal ulp, sfac;
00050     integer mode;
00051     char path[3], dist[1];
00052     doublereal unfl, rexp;
00053     char type__[1];
00054     doublereal texp;
00055     extern /* Subroutine */ int zrot_(integer *, doublecomplex *, integer *, 
00056             doublecomplex *, integer *, doublereal *, doublecomplex *);
00057     doublecomplex star1, plus1, plus2;
00058     doublereal bscal;
00059     extern logical lsame_(char *, char *);
00060     doublereal tscal, anorm, bnorm, tleft;
00061     logical upper;
00062     extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *, 
00063             doublecomplex *, integer *), zrotg_(doublecomplex *, 
00064             doublecomplex *, doublereal *, doublecomplex *), zswap_(integer *, 
00065              doublecomplex *, integer *, doublecomplex *, integer *), zlatb4_(
00066             char *, integer *, integer *, integer *, char *, integer *, 
00067             integer *, doublereal *, integer *, doublereal *, char *), dlabad_(doublereal *, doublereal *);
00068     extern doublereal dlamch_(char *), dlarnd_(integer *, integer *);
00069     extern /* Subroutine */ int zdscal_(integer *, doublereal *, 
00070             doublecomplex *, integer *);
00071     doublereal bignum, cndnum;
00072     extern /* Subroutine */ int dlarnv_(integer *, integer *, integer *, 
00073             doublereal *);
00074     extern integer izamax_(integer *, doublecomplex *, integer *);
00075     extern /* Double Complex */ VOID zlarnd_(doublecomplex *, integer *, 
00076             integer *);
00077     integer jcount;
00078     extern /* Subroutine */ int zlatms_(integer *, integer *, char *, integer 
00079             *, char *, doublereal *, integer *, doublereal *, doublereal *, 
00080             integer *, integer *, char *, doublecomplex *, integer *, 
00081             doublecomplex *, integer *);
00082     doublereal smlnum;
00083     extern /* Subroutine */ int zlarnv_(integer *, integer *, integer *, 
00084             doublecomplex *);
00085 
00086 
00087 /*  -- LAPACK test routine (version 3.1) -- */
00088 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00089 /*     November 2006 */
00090 
00091 /*     .. Scalar Arguments .. */
00092 /*     .. */
00093 /*     .. Array Arguments .. */
00094 /*     .. */
00095 
00096 /*  Purpose */
00097 /*  ======= */
00098 
00099 /*  ZLATTR generates a triangular test matrix in 2-dimensional storage. */
00100 /*  IMAT and UPLO uniquely specify the properties of the test matrix, */
00101 /*  which is returned in the array A. */
00102 
00103 /*  Arguments */
00104 /*  ========= */
00105 
00106 /*  IMAT    (input) INTEGER */
00107 /*          An integer key describing which matrix to generate for this */
00108 /*          path. */
00109 
00110 /*  UPLO    (input) CHARACTER*1 */
00111 /*          Specifies whether the matrix A will be upper or lower */
00112 /*          triangular. */
00113 /*          = 'U':  Upper triangular */
00114 /*          = 'L':  Lower triangular */
00115 
00116 /*  TRANS   (input) CHARACTER*1 */
00117 /*          Specifies whether the matrix or its transpose will be used. */
00118 /*          = 'N':  No transpose */
00119 /*          = 'T':  Transpose */
00120 /*          = 'C':  Conjugate transpose */
00121 
00122 /*  DIAG    (output) CHARACTER*1 */
00123 /*          Specifies whether or not the matrix A is unit triangular. */
00124 /*          = 'N':  Non-unit triangular */
00125 /*          = 'U':  Unit triangular */
00126 
00127 /*  ISEED   (input/output) INTEGER array, dimension (4) */
00128 /*          The seed vector for the random number generator (used in */
00129 /*          ZLATMS).  Modified on exit. */
00130 
00131 /*  N       (input) INTEGER */
00132 /*          The order of the matrix to be generated. */
00133 
00134 /*  A       (output) COMPLEX*16 array, dimension (LDA,N) */
00135 /*          The triangular matrix A.  If UPLO = 'U', the leading N x N */
00136 /*          upper triangular part of the array A contains the upper */
00137 /*          triangular matrix, and the strictly lower triangular part of */
00138 /*          A is not referenced.  If UPLO = 'L', the leading N x N lower */
00139 /*          triangular part of the array A contains the lower triangular */
00140 /*          matrix and the strictly upper triangular part of A is not */
00141 /*          referenced. */
00142 
00143 /*  LDA     (input) INTEGER */
00144 /*          The leading dimension of the array A.  LDA >= max(1,N). */
00145 
00146 /*  B       (output) COMPLEX*16 array, dimension (N) */
00147 /*          The right hand side vector, if IMAT > 10. */
00148 
00149 /*  WORK    (workspace) COMPLEX*16 array, dimension (2*N) */
00150 
00151 /*  RWORK   (workspace) DOUBLE PRECISION array, dimension (N) */
00152 
00153 /*  INFO    (output) INTEGER */
00154 /*          = 0:  successful exit */
00155 /*          < 0:  if INFO = -i, the i-th argument had an illegal value */
00156 
00157 /*  ===================================================================== */
00158 
00159 /*     .. Parameters .. */
00160 /*     .. */
00161 /*     .. Local Scalars .. */
00162 /*     .. */
00163 /*     .. External Functions .. */
00164 /*     .. */
00165 /*     .. External Subroutines .. */
00166 /*     .. */
00167 /*     .. Intrinsic Functions .. */
00168 /*     .. */
00169 /*     .. Executable Statements .. */
00170 
00171     /* Parameter adjustments */
00172     --iseed;
00173     a_dim1 = *lda;
00174     a_offset = 1 + a_dim1;
00175     a -= a_offset;
00176     --b;
00177     --work;
00178     --rwork;
00179 
00180     /* Function Body */
00181     s_copy(path, "Zomplex precision", (ftnlen)1, (ftnlen)17);
00182     s_copy(path + 1, "TR", (ftnlen)2, (ftnlen)2);
00183     unfl = dlamch_("Safe minimum");
00184     ulp = dlamch_("Epsilon") * dlamch_("Base");
00185     smlnum = unfl;
00186     bignum = (1. - ulp) / smlnum;
00187     dlabad_(&smlnum, &bignum);
00188     if (*imat >= 7 && *imat <= 10 || *imat == 18) {
00189         *(unsigned char *)diag = 'U';
00190     } else {
00191         *(unsigned char *)diag = 'N';
00192     }
00193     *info = 0;
00194 
00195 /*     Quick return if N.LE.0. */
00196 
00197     if (*n <= 0) {
00198         return 0;
00199     }
00200 
00201 /*     Call ZLATB4 to set parameters for CLATMS. */
00202 
00203     upper = lsame_(uplo, "U");
00204     if (upper) {
00205         zlatb4_(path, imat, n, n, type__, &kl, &ku, &anorm, &mode, &cndnum, 
00206                 dist);
00207     } else {
00208         i__1 = -(*imat);
00209         zlatb4_(path, &i__1, n, n, type__, &kl, &ku, &anorm, &mode, &cndnum, 
00210                 dist);
00211     }
00212 
00213 /*     IMAT <= 6:  Non-unit triangular matrix */
00214 
00215     if (*imat <= 6) {
00216         zlatms_(n, n, dist, &iseed[1], type__, &rwork[1], &mode, &cndnum, &
00217                 anorm, &kl, &ku, "No packing", &a[a_offset], lda, &work[1], 
00218                 info);
00219 
00220 /*     IMAT > 6:  Unit triangular matrix */
00221 /*     The diagonal is deliberately set to something other than 1. */
00222 
00223 /*     IMAT = 7:  Matrix is the identity */
00224 
00225     } else if (*imat == 7) {
00226         if (upper) {
00227             i__1 = *n;
00228             for (j = 1; j <= i__1; ++j) {
00229                 i__2 = j - 1;
00230                 for (i__ = 1; i__ <= i__2; ++i__) {
00231                     i__3 = i__ + j * a_dim1;
00232                     a[i__3].r = 0., a[i__3].i = 0.;
00233 /* L10: */
00234                 }
00235                 i__2 = j + j * a_dim1;
00236                 a[i__2].r = (doublereal) j, a[i__2].i = 0.;
00237 /* L20: */
00238             }
00239         } else {
00240             i__1 = *n;
00241             for (j = 1; j <= i__1; ++j) {
00242                 i__2 = j + j * a_dim1;
00243                 a[i__2].r = (doublereal) j, a[i__2].i = 0.;
00244                 i__2 = *n;
00245                 for (i__ = j + 1; i__ <= i__2; ++i__) {
00246                     i__3 = i__ + j * a_dim1;
00247                     a[i__3].r = 0., a[i__3].i = 0.;
00248 /* L30: */
00249                 }
00250 /* L40: */
00251             }
00252         }
00253 
00254 /*     IMAT > 7:  Non-trivial unit triangular matrix */
00255 
00256 /*     Generate a unit triangular matrix T with condition CNDNUM by */
00257 /*     forming a triangular matrix with known singular values and */
00258 /*     filling in the zero entries with Givens rotations. */
00259 
00260     } else if (*imat <= 10) {
00261         if (upper) {
00262             i__1 = *n;
00263             for (j = 1; j <= i__1; ++j) {
00264                 i__2 = j - 1;
00265                 for (i__ = 1; i__ <= i__2; ++i__) {
00266                     i__3 = i__ + j * a_dim1;
00267                     a[i__3].r = 0., a[i__3].i = 0.;
00268 /* L50: */
00269                 }
00270                 i__2 = j + j * a_dim1;
00271                 a[i__2].r = (doublereal) j, a[i__2].i = 0.;
00272 /* L60: */
00273             }
00274         } else {
00275             i__1 = *n;
00276             for (j = 1; j <= i__1; ++j) {
00277                 i__2 = j + j * a_dim1;
00278                 a[i__2].r = (doublereal) j, a[i__2].i = 0.;
00279                 i__2 = *n;
00280                 for (i__ = j + 1; i__ <= i__2; ++i__) {
00281                     i__3 = i__ + j * a_dim1;
00282                     a[i__3].r = 0., a[i__3].i = 0.;
00283 /* L70: */
00284                 }
00285 /* L80: */
00286             }
00287         }
00288 
00289 /*        Since the trace of a unit triangular matrix is 1, the product */
00290 /*        of its singular values must be 1.  Let s = sqrt(CNDNUM), */
00291 /*        x = sqrt(s) - 1/sqrt(s), y = sqrt(2/(n-2))*x, and z = x**2. */
00292 /*        The following triangular matrix has singular values s, 1, 1, */
00293 /*        ..., 1, 1/s: */
00294 
00295 /*        1  y  y  y  ...  y  y  z */
00296 /*           1  0  0  ...  0  0  y */
00297 /*              1  0  ...  0  0  y */
00298 /*                 .  ...  .  .  . */
00299 /*                     .   .  .  . */
00300 /*                         1  0  y */
00301 /*                            1  y */
00302 /*                               1 */
00303 
00304 /*        To fill in the zeros, we first multiply by a matrix with small */
00305 /*        condition number of the form */
00306 
00307 /*        1  0  0  0  0  ... */
00308 /*           1  +  *  0  0  ... */
00309 /*              1  +  0  0  0 */
00310 /*                 1  +  *  0  0 */
00311 /*                    1  +  0  0 */
00312 /*                       ... */
00313 /*                          1  +  0 */
00314 /*                             1  0 */
00315 /*                                1 */
00316 
00317 /*        Each element marked with a '*' is formed by taking the product */
00318 /*        of the adjacent elements marked with '+'.  The '*'s can be */
00319 /*        chosen freely, and the '+'s are chosen so that the inverse of */
00320 /*        T will have elements of the same magnitude as T.  If the *'s in */
00321 /*        both T and inv(T) have small magnitude, T is well conditioned. */
00322 /*        The two offdiagonals of T are stored in WORK. */
00323 
00324 /*        The product of these two matrices has the form */
00325 
00326 /*        1  y  y  y  y  y  .  y  y  z */
00327 /*           1  +  *  0  0  .  0  0  y */
00328 /*              1  +  0  0  .  0  0  y */
00329 /*                 1  +  *  .  .  .  . */
00330 /*                    1  +  .  .  .  . */
00331 /*                       .  .  .  .  . */
00332 /*                          .  .  .  . */
00333 /*                             1  +  y */
00334 /*                                1  y */
00335 /*                                   1 */
00336 
00337 /*        Now we multiply by Givens rotations, using the fact that */
00338 
00339 /*              [  c   s ] [  1   w ] [ -c  -s ] =  [  1  -w ] */
00340 /*              [ -s   c ] [  0   1 ] [  s  -c ]    [  0   1 ] */
00341 /*        and */
00342 /*              [ -c  -s ] [  1   0 ] [  c   s ] =  [  1   0 ] */
00343 /*              [  s  -c ] [  w   1 ] [ -s   c ]    [ -w   1 ] */
00344 
00345 /*        where c = w / sqrt(w**2+4) and s = 2 / sqrt(w**2+4). */
00346 
00347         zlarnd_(&z__2, &c__5, &iseed[1]);
00348         z__1.r = z__2.r * .25, z__1.i = z__2.i * .25;
00349         star1.r = z__1.r, star1.i = z__1.i;
00350         sfac = .5;
00351         zlarnd_(&z__2, &c__5, &iseed[1]);
00352         z__1.r = sfac * z__2.r, z__1.i = sfac * z__2.i;
00353         plus1.r = z__1.r, plus1.i = z__1.i;
00354         i__1 = *n;
00355         for (j = 1; j <= i__1; j += 2) {
00356             z_div(&z__1, &star1, &plus1);
00357             plus2.r = z__1.r, plus2.i = z__1.i;
00358             i__2 = j;
00359             work[i__2].r = plus1.r, work[i__2].i = plus1.i;
00360             i__2 = *n + j;
00361             work[i__2].r = star1.r, work[i__2].i = star1.i;
00362             if (j + 1 <= *n) {
00363                 i__2 = j + 1;
00364                 work[i__2].r = plus2.r, work[i__2].i = plus2.i;
00365                 i__2 = *n + j + 1;
00366                 work[i__2].r = 0., work[i__2].i = 0.;
00367                 z_div(&z__1, &star1, &plus2);
00368                 plus1.r = z__1.r, plus1.i = z__1.i;
00369                 rexp = dlarnd_(&c__2, &iseed[1]);
00370                 if (rexp < 0.) {
00371                     d__2 = 1. - rexp;
00372                     d__1 = -pow_dd(&sfac, &d__2);
00373                     zlarnd_(&z__2, &c__5, &iseed[1]);
00374                     z__1.r = d__1 * z__2.r, z__1.i = d__1 * z__2.i;
00375                     star1.r = z__1.r, star1.i = z__1.i;
00376                 } else {
00377                     d__2 = rexp + 1.;
00378                     d__1 = pow_dd(&sfac, &d__2);
00379                     zlarnd_(&z__2, &c__5, &iseed[1]);
00380                     z__1.r = d__1 * z__2.r, z__1.i = d__1 * z__2.i;
00381                     star1.r = z__1.r, star1.i = z__1.i;
00382                 }
00383             }
00384 /* L90: */
00385         }
00386 
00387         x = sqrt(cndnum) - 1 / sqrt(cndnum);
00388         if (*n > 2) {
00389             y = sqrt(2. / (*n - 2)) * x;
00390         } else {
00391             y = 0.;
00392         }
00393         z__ = x * x;
00394 
00395         if (upper) {
00396             if (*n > 3) {
00397                 i__1 = *n - 3;
00398                 i__2 = *lda + 1;
00399                 zcopy_(&i__1, &work[1], &c__1, &a[a_dim1 * 3 + 2], &i__2);
00400                 if (*n > 4) {
00401                     i__1 = *n - 4;
00402                     i__2 = *lda + 1;
00403                     zcopy_(&i__1, &work[*n + 1], &c__1, &a[(a_dim1 << 2) + 2], 
00404                              &i__2);
00405                 }
00406             }
00407             i__1 = *n - 1;
00408             for (j = 2; j <= i__1; ++j) {
00409                 i__2 = j * a_dim1 + 1;
00410                 a[i__2].r = y, a[i__2].i = 0.;
00411                 i__2 = j + *n * a_dim1;
00412                 a[i__2].r = y, a[i__2].i = 0.;
00413 /* L100: */
00414             }
00415             i__1 = *n * a_dim1 + 1;
00416             a[i__1].r = z__, a[i__1].i = 0.;
00417         } else {
00418             if (*n > 3) {
00419                 i__1 = *n - 3;
00420                 i__2 = *lda + 1;
00421                 zcopy_(&i__1, &work[1], &c__1, &a[(a_dim1 << 1) + 3], &i__2);
00422                 if (*n > 4) {
00423                     i__1 = *n - 4;
00424                     i__2 = *lda + 1;
00425                     zcopy_(&i__1, &work[*n + 1], &c__1, &a[(a_dim1 << 1) + 4], 
00426                              &i__2);
00427                 }
00428             }
00429             i__1 = *n - 1;
00430             for (j = 2; j <= i__1; ++j) {
00431                 i__2 = j + a_dim1;
00432                 a[i__2].r = y, a[i__2].i = 0.;
00433                 i__2 = *n + j * a_dim1;
00434                 a[i__2].r = y, a[i__2].i = 0.;
00435 /* L110: */
00436             }
00437             i__1 = *n + a_dim1;
00438             a[i__1].r = z__, a[i__1].i = 0.;
00439         }
00440 
00441 /*        Fill in the zeros using Givens rotations. */
00442 
00443         if (upper) {
00444             i__1 = *n - 1;
00445             for (j = 1; j <= i__1; ++j) {
00446                 i__2 = j + (j + 1) * a_dim1;
00447                 ra.r = a[i__2].r, ra.i = a[i__2].i;
00448                 rb.r = 2., rb.i = 0.;
00449                 zrotg_(&ra, &rb, &c__, &s);
00450 
00451 /*              Multiply by [ c  s; -conjg(s)  c] on the left. */
00452 
00453                 if (*n > j + 1) {
00454                     i__2 = *n - j - 1;
00455                     zrot_(&i__2, &a[j + (j + 2) * a_dim1], lda, &a[j + 1 + (j 
00456                             + 2) * a_dim1], lda, &c__, &s);
00457                 }
00458 
00459 /*              Multiply by [-c -s;  conjg(s) -c] on the right. */
00460 
00461                 if (j > 1) {
00462                     i__2 = j - 1;
00463                     d__1 = -c__;
00464                     z__1.r = -s.r, z__1.i = -s.i;
00465                     zrot_(&i__2, &a[(j + 1) * a_dim1 + 1], &c__1, &a[j * 
00466                             a_dim1 + 1], &c__1, &d__1, &z__1);
00467                 }
00468 
00469 /*              Negate A(J,J+1). */
00470 
00471                 i__2 = j + (j + 1) * a_dim1;
00472                 i__3 = j + (j + 1) * a_dim1;
00473                 z__1.r = -a[i__3].r, z__1.i = -a[i__3].i;
00474                 a[i__2].r = z__1.r, a[i__2].i = z__1.i;
00475 /* L120: */
00476             }
00477         } else {
00478             i__1 = *n - 1;
00479             for (j = 1; j <= i__1; ++j) {
00480                 i__2 = j + 1 + j * a_dim1;
00481                 ra.r = a[i__2].r, ra.i = a[i__2].i;
00482                 rb.r = 2., rb.i = 0.;
00483                 zrotg_(&ra, &rb, &c__, &s);
00484                 d_cnjg(&z__1, &s);
00485                 s.r = z__1.r, s.i = z__1.i;
00486 
00487 /*              Multiply by [ c -s;  conjg(s) c] on the right. */
00488 
00489                 if (*n > j + 1) {
00490                     i__2 = *n - j - 1;
00491                     z__1.r = -s.r, z__1.i = -s.i;
00492                     zrot_(&i__2, &a[j + 2 + (j + 1) * a_dim1], &c__1, &a[j + 
00493                             2 + j * a_dim1], &c__1, &c__, &z__1);
00494                 }
00495 
00496 /*              Multiply by [-c  s; -conjg(s) -c] on the left. */
00497 
00498                 if (j > 1) {
00499                     i__2 = j - 1;
00500                     d__1 = -c__;
00501                     zrot_(&i__2, &a[j + a_dim1], lda, &a[j + 1 + a_dim1], lda, 
00502                              &d__1, &s);
00503                 }
00504 
00505 /*              Negate A(J+1,J). */
00506 
00507                 i__2 = j + 1 + j * a_dim1;
00508                 i__3 = j + 1 + j * a_dim1;
00509                 z__1.r = -a[i__3].r, z__1.i = -a[i__3].i;
00510                 a[i__2].r = z__1.r, a[i__2].i = z__1.i;
00511 /* L130: */
00512             }
00513         }
00514 
00515 /*     IMAT > 10:  Pathological test cases.  These triangular matrices */
00516 /*     are badly scaled or badly conditioned, so when used in solving a */
00517 /*     triangular system they may cause overflow in the solution vector. */
00518 
00519     } else if (*imat == 11) {
00520 
00521 /*        Type 11:  Generate a triangular matrix with elements between */
00522 /*        -1 and 1. Give the diagonal norm 2 to make it well-conditioned. */
00523 /*        Make the right hand side large so that it requires scaling. */
00524 
00525         if (upper) {
00526             i__1 = *n;
00527             for (j = 1; j <= i__1; ++j) {
00528                 i__2 = j - 1;
00529                 zlarnv_(&c__4, &iseed[1], &i__2, &a[j * a_dim1 + 1]);
00530                 i__2 = j + j * a_dim1;
00531                 zlarnd_(&z__2, &c__5, &iseed[1]);
00532                 z__1.r = z__2.r * 2., z__1.i = z__2.i * 2.;
00533                 a[i__2].r = z__1.r, a[i__2].i = z__1.i;
00534 /* L140: */
00535             }
00536         } else {
00537             i__1 = *n;
00538             for (j = 1; j <= i__1; ++j) {
00539                 if (j < *n) {
00540                     i__2 = *n - j;
00541                     zlarnv_(&c__4, &iseed[1], &i__2, &a[j + 1 + j * a_dim1]);
00542                 }
00543                 i__2 = j + j * a_dim1;
00544                 zlarnd_(&z__2, &c__5, &iseed[1]);
00545                 z__1.r = z__2.r * 2., z__1.i = z__2.i * 2.;
00546                 a[i__2].r = z__1.r, a[i__2].i = z__1.i;
00547 /* L150: */
00548             }
00549         }
00550 
00551 /*        Set the right hand side so that the largest value is BIGNUM. */
00552 
00553         zlarnv_(&c__2, &iseed[1], n, &b[1]);
00554         iy = izamax_(n, &b[1], &c__1);
00555         bnorm = z_abs(&b[iy]);
00556         bscal = bignum / max(1.,bnorm);
00557         zdscal_(n, &bscal, &b[1], &c__1);
00558 
00559     } else if (*imat == 12) {
00560 
00561 /*        Type 12:  Make the first diagonal element in the solve small to */
00562 /*        cause immediate overflow when dividing by T(j,j). */
00563 /*        In type 12, the offdiagonal elements are small (CNORM(j) < 1). */
00564 
00565         zlarnv_(&c__2, &iseed[1], n, &b[1]);
00566 /* Computing MAX */
00567         d__1 = 1., d__2 = (doublereal) (*n - 1);
00568         tscal = 1. / max(d__1,d__2);
00569         if (upper) {
00570             i__1 = *n;
00571             for (j = 1; j <= i__1; ++j) {
00572                 i__2 = j - 1;
00573                 zlarnv_(&c__4, &iseed[1], &i__2, &a[j * a_dim1 + 1]);
00574                 i__2 = j - 1;
00575                 zdscal_(&i__2, &tscal, &a[j * a_dim1 + 1], &c__1);
00576                 i__2 = j + j * a_dim1;
00577                 zlarnd_(&z__1, &c__5, &iseed[1]);
00578                 a[i__2].r = z__1.r, a[i__2].i = z__1.i;
00579 /* L160: */
00580             }
00581             i__1 = *n + *n * a_dim1;
00582             i__2 = *n + *n * a_dim1;
00583             z__1.r = smlnum * a[i__2].r, z__1.i = smlnum * a[i__2].i;
00584             a[i__1].r = z__1.r, a[i__1].i = z__1.i;
00585         } else {
00586             i__1 = *n;
00587             for (j = 1; j <= i__1; ++j) {
00588                 if (j < *n) {
00589                     i__2 = *n - j;
00590                     zlarnv_(&c__4, &iseed[1], &i__2, &a[j + 1 + j * a_dim1]);
00591                     i__2 = *n - j;
00592                     zdscal_(&i__2, &tscal, &a[j + 1 + j * a_dim1], &c__1);
00593                 }
00594                 i__2 = j + j * a_dim1;
00595                 zlarnd_(&z__1, &c__5, &iseed[1]);
00596                 a[i__2].r = z__1.r, a[i__2].i = z__1.i;
00597 /* L170: */
00598             }
00599             i__1 = a_dim1 + 1;
00600             i__2 = a_dim1 + 1;
00601             z__1.r = smlnum * a[i__2].r, z__1.i = smlnum * a[i__2].i;
00602             a[i__1].r = z__1.r, a[i__1].i = z__1.i;
00603         }
00604 
00605     } else if (*imat == 13) {
00606 
00607 /*        Type 13:  Make the first diagonal element in the solve small to */
00608 /*        cause immediate overflow when dividing by T(j,j). */
00609 /*        In type 13, the offdiagonal elements are O(1) (CNORM(j) > 1). */
00610 
00611         zlarnv_(&c__2, &iseed[1], n, &b[1]);
00612         if (upper) {
00613             i__1 = *n;
00614             for (j = 1; j <= i__1; ++j) {
00615                 i__2 = j - 1;
00616                 zlarnv_(&c__4, &iseed[1], &i__2, &a[j * a_dim1 + 1]);
00617                 i__2 = j + j * a_dim1;
00618                 zlarnd_(&z__1, &c__5, &iseed[1]);
00619                 a[i__2].r = z__1.r, a[i__2].i = z__1.i;
00620 /* L180: */
00621             }
00622             i__1 = *n + *n * a_dim1;
00623             i__2 = *n + *n * a_dim1;
00624             z__1.r = smlnum * a[i__2].r, z__1.i = smlnum * a[i__2].i;
00625             a[i__1].r = z__1.r, a[i__1].i = z__1.i;
00626         } else {
00627             i__1 = *n;
00628             for (j = 1; j <= i__1; ++j) {
00629                 if (j < *n) {
00630                     i__2 = *n - j;
00631                     zlarnv_(&c__4, &iseed[1], &i__2, &a[j + 1 + j * a_dim1]);
00632                 }
00633                 i__2 = j + j * a_dim1;
00634                 zlarnd_(&z__1, &c__5, &iseed[1]);
00635                 a[i__2].r = z__1.r, a[i__2].i = z__1.i;
00636 /* L190: */
00637             }
00638             i__1 = a_dim1 + 1;
00639             i__2 = a_dim1 + 1;
00640             z__1.r = smlnum * a[i__2].r, z__1.i = smlnum * a[i__2].i;
00641             a[i__1].r = z__1.r, a[i__1].i = z__1.i;
00642         }
00643 
00644     } else if (*imat == 14) {
00645 
00646 /*        Type 14:  T is diagonal with small numbers on the diagonal to */
00647 /*        make the growth factor underflow, but a small right hand side */
00648 /*        chosen so that the solution does not overflow. */
00649 
00650         if (upper) {
00651             jcount = 1;
00652             for (j = *n; j >= 1; --j) {
00653                 i__1 = j - 1;
00654                 for (i__ = 1; i__ <= i__1; ++i__) {
00655                     i__2 = i__ + j * a_dim1;
00656                     a[i__2].r = 0., a[i__2].i = 0.;
00657 /* L200: */
00658                 }
00659                 if (jcount <= 2) {
00660                     i__1 = j + j * a_dim1;
00661                     zlarnd_(&z__2, &c__5, &iseed[1]);
00662                     z__1.r = smlnum * z__2.r, z__1.i = smlnum * z__2.i;
00663                     a[i__1].r = z__1.r, a[i__1].i = z__1.i;
00664                 } else {
00665                     i__1 = j + j * a_dim1;
00666                     zlarnd_(&z__1, &c__5, &iseed[1]);
00667                     a[i__1].r = z__1.r, a[i__1].i = z__1.i;
00668                 }
00669                 ++jcount;
00670                 if (jcount > 4) {
00671                     jcount = 1;
00672                 }
00673 /* L210: */
00674             }
00675         } else {
00676             jcount = 1;
00677             i__1 = *n;
00678             for (j = 1; j <= i__1; ++j) {
00679                 i__2 = *n;
00680                 for (i__ = j + 1; i__ <= i__2; ++i__) {
00681                     i__3 = i__ + j * a_dim1;
00682                     a[i__3].r = 0., a[i__3].i = 0.;
00683 /* L220: */
00684                 }
00685                 if (jcount <= 2) {
00686                     i__2 = j + j * a_dim1;
00687                     zlarnd_(&z__2, &c__5, &iseed[1]);
00688                     z__1.r = smlnum * z__2.r, z__1.i = smlnum * z__2.i;
00689                     a[i__2].r = z__1.r, a[i__2].i = z__1.i;
00690                 } else {
00691                     i__2 = j + j * a_dim1;
00692                     zlarnd_(&z__1, &c__5, &iseed[1]);
00693                     a[i__2].r = z__1.r, a[i__2].i = z__1.i;
00694                 }
00695                 ++jcount;
00696                 if (jcount > 4) {
00697                     jcount = 1;
00698                 }
00699 /* L230: */
00700             }
00701         }
00702 
00703 /*        Set the right hand side alternately zero and small. */
00704 
00705         if (upper) {
00706             b[1].r = 0., b[1].i = 0.;
00707             for (i__ = *n; i__ >= 2; i__ += -2) {
00708                 i__1 = i__;
00709                 b[i__1].r = 0., b[i__1].i = 0.;
00710                 i__1 = i__ - 1;
00711                 zlarnd_(&z__2, &c__5, &iseed[1]);
00712                 z__1.r = smlnum * z__2.r, z__1.i = smlnum * z__2.i;
00713                 b[i__1].r = z__1.r, b[i__1].i = z__1.i;
00714 /* L240: */
00715             }
00716         } else {
00717             i__1 = *n;
00718             b[i__1].r = 0., b[i__1].i = 0.;
00719             i__1 = *n - 1;
00720             for (i__ = 1; i__ <= i__1; i__ += 2) {
00721                 i__2 = i__;
00722                 b[i__2].r = 0., b[i__2].i = 0.;
00723                 i__2 = i__ + 1;
00724                 zlarnd_(&z__2, &c__5, &iseed[1]);
00725                 z__1.r = smlnum * z__2.r, z__1.i = smlnum * z__2.i;
00726                 b[i__2].r = z__1.r, b[i__2].i = z__1.i;
00727 /* L250: */
00728             }
00729         }
00730 
00731     } else if (*imat == 15) {
00732 
00733 /*        Type 15:  Make the diagonal elements small to cause gradual */
00734 /*        overflow when dividing by T(j,j).  To control the amount of */
00735 /*        scaling needed, the matrix is bidiagonal. */
00736 
00737 /* Computing MAX */
00738         d__1 = 1., d__2 = (doublereal) (*n - 1);
00739         texp = 1. / max(d__1,d__2);
00740         tscal = pow_dd(&smlnum, &texp);
00741         zlarnv_(&c__4, &iseed[1], n, &b[1]);
00742         if (upper) {
00743             i__1 = *n;
00744             for (j = 1; j <= i__1; ++j) {
00745                 i__2 = j - 2;
00746                 for (i__ = 1; i__ <= i__2; ++i__) {
00747                     i__3 = i__ + j * a_dim1;
00748                     a[i__3].r = 0., a[i__3].i = 0.;
00749 /* L260: */
00750                 }
00751                 if (j > 1) {
00752                     i__2 = j - 1 + j * a_dim1;
00753                     a[i__2].r = -1., a[i__2].i = -1.;
00754                 }
00755                 i__2 = j + j * a_dim1;
00756                 zlarnd_(&z__2, &c__5, &iseed[1]);
00757                 z__1.r = tscal * z__2.r, z__1.i = tscal * z__2.i;
00758                 a[i__2].r = z__1.r, a[i__2].i = z__1.i;
00759 /* L270: */
00760             }
00761             i__1 = *n;
00762             b[i__1].r = 1., b[i__1].i = 1.;
00763         } else {
00764             i__1 = *n;
00765             for (j = 1; j <= i__1; ++j) {
00766                 i__2 = *n;
00767                 for (i__ = j + 2; i__ <= i__2; ++i__) {
00768                     i__3 = i__ + j * a_dim1;
00769                     a[i__3].r = 0., a[i__3].i = 0.;
00770 /* L280: */
00771                 }
00772                 if (j < *n) {
00773                     i__2 = j + 1 + j * a_dim1;
00774                     a[i__2].r = -1., a[i__2].i = -1.;
00775                 }
00776                 i__2 = j + j * a_dim1;
00777                 zlarnd_(&z__2, &c__5, &iseed[1]);
00778                 z__1.r = tscal * z__2.r, z__1.i = tscal * z__2.i;
00779                 a[i__2].r = z__1.r, a[i__2].i = z__1.i;
00780 /* L290: */
00781             }
00782             b[1].r = 1., b[1].i = 1.;
00783         }
00784 
00785     } else if (*imat == 16) {
00786 
00787 /*        Type 16:  One zero diagonal element. */
00788 
00789         iy = *n / 2 + 1;
00790         if (upper) {
00791             i__1 = *n;
00792             for (j = 1; j <= i__1; ++j) {
00793                 i__2 = j - 1;
00794                 zlarnv_(&c__4, &iseed[1], &i__2, &a[j * a_dim1 + 1]);
00795                 if (j != iy) {
00796                     i__2 = j + j * a_dim1;
00797                     zlarnd_(&z__2, &c__5, &iseed[1]);
00798                     z__1.r = z__2.r * 2., z__1.i = z__2.i * 2.;
00799                     a[i__2].r = z__1.r, a[i__2].i = z__1.i;
00800                 } else {
00801                     i__2 = j + j * a_dim1;
00802                     a[i__2].r = 0., a[i__2].i = 0.;
00803                 }
00804 /* L300: */
00805             }
00806         } else {
00807             i__1 = *n;
00808             for (j = 1; j <= i__1; ++j) {
00809                 if (j < *n) {
00810                     i__2 = *n - j;
00811                     zlarnv_(&c__4, &iseed[1], &i__2, &a[j + 1 + j * a_dim1]);
00812                 }
00813                 if (j != iy) {
00814                     i__2 = j + j * a_dim1;
00815                     zlarnd_(&z__2, &c__5, &iseed[1]);
00816                     z__1.r = z__2.r * 2., z__1.i = z__2.i * 2.;
00817                     a[i__2].r = z__1.r, a[i__2].i = z__1.i;
00818                 } else {
00819                     i__2 = j + j * a_dim1;
00820                     a[i__2].r = 0., a[i__2].i = 0.;
00821                 }
00822 /* L310: */
00823             }
00824         }
00825         zlarnv_(&c__2, &iseed[1], n, &b[1]);
00826         zdscal_(n, &c_b92, &b[1], &c__1);
00827 
00828     } else if (*imat == 17) {
00829 
00830 /*        Type 17:  Make the offdiagonal elements large to cause overflow */
00831 /*        when adding a column of T.  In the non-transposed case, the */
00832 /*        matrix is constructed to cause overflow when adding a column in */
00833 /*        every other step. */
00834 
00835         tscal = unfl / ulp;
00836         tscal = (1. - ulp) / tscal;
00837         i__1 = *n;
00838         for (j = 1; j <= i__1; ++j) {
00839             i__2 = *n;
00840             for (i__ = 1; i__ <= i__2; ++i__) {
00841                 i__3 = i__ + j * a_dim1;
00842                 a[i__3].r = 0., a[i__3].i = 0.;
00843 /* L320: */
00844             }
00845 /* L330: */
00846         }
00847         texp = 1.;
00848         if (upper) {
00849             for (j = *n; j >= 2; j += -2) {
00850                 i__1 = j * a_dim1 + 1;
00851                 d__1 = -tscal / (doublereal) (*n + 1);
00852                 a[i__1].r = d__1, a[i__1].i = 0.;
00853                 i__1 = j + j * a_dim1;
00854                 a[i__1].r = 1., a[i__1].i = 0.;
00855                 i__1 = j;
00856                 d__1 = texp * (1. - ulp);
00857                 b[i__1].r = d__1, b[i__1].i = 0.;
00858                 i__1 = (j - 1) * a_dim1 + 1;
00859                 d__1 = -(tscal / (doublereal) (*n + 1)) / (doublereal) (*n + 
00860                         2);
00861                 a[i__1].r = d__1, a[i__1].i = 0.;
00862                 i__1 = j - 1 + (j - 1) * a_dim1;
00863                 a[i__1].r = 1., a[i__1].i = 0.;
00864                 i__1 = j - 1;
00865                 d__1 = texp * (doublereal) (*n * *n + *n - 1);
00866                 b[i__1].r = d__1, b[i__1].i = 0.;
00867                 texp *= 2.;
00868 /* L340: */
00869             }
00870             d__1 = (doublereal) (*n + 1) / (doublereal) (*n + 2) * tscal;
00871             b[1].r = d__1, b[1].i = 0.;
00872         } else {
00873             i__1 = *n - 1;
00874             for (j = 1; j <= i__1; j += 2) {
00875                 i__2 = *n + j * a_dim1;
00876                 d__1 = -tscal / (doublereal) (*n + 1);
00877                 a[i__2].r = d__1, a[i__2].i = 0.;
00878                 i__2 = j + j * a_dim1;
00879                 a[i__2].r = 1., a[i__2].i = 0.;
00880                 i__2 = j;
00881                 d__1 = texp * (1. - ulp);
00882                 b[i__2].r = d__1, b[i__2].i = 0.;
00883                 i__2 = *n + (j + 1) * a_dim1;
00884                 d__1 = -(tscal / (doublereal) (*n + 1)) / (doublereal) (*n + 
00885                         2);
00886                 a[i__2].r = d__1, a[i__2].i = 0.;
00887                 i__2 = j + 1 + (j + 1) * a_dim1;
00888                 a[i__2].r = 1., a[i__2].i = 0.;
00889                 i__2 = j + 1;
00890                 d__1 = texp * (doublereal) (*n * *n + *n - 1);
00891                 b[i__2].r = d__1, b[i__2].i = 0.;
00892                 texp *= 2.;
00893 /* L350: */
00894             }
00895             i__1 = *n;
00896             d__1 = (doublereal) (*n + 1) / (doublereal) (*n + 2) * tscal;
00897             b[i__1].r = d__1, b[i__1].i = 0.;
00898         }
00899 
00900     } else if (*imat == 18) {
00901 
00902 /*        Type 18:  Generate a unit triangular matrix with elements */
00903 /*        between -1 and 1, and make the right hand side large so that it */
00904 /*        requires scaling. */
00905 
00906         if (upper) {
00907             i__1 = *n;
00908             for (j = 1; j <= i__1; ++j) {
00909                 i__2 = j - 1;
00910                 zlarnv_(&c__4, &iseed[1], &i__2, &a[j * a_dim1 + 1]);
00911                 i__2 = j + j * a_dim1;
00912                 a[i__2].r = 0., a[i__2].i = 0.;
00913 /* L360: */
00914             }
00915         } else {
00916             i__1 = *n;
00917             for (j = 1; j <= i__1; ++j) {
00918                 if (j < *n) {
00919                     i__2 = *n - j;
00920                     zlarnv_(&c__4, &iseed[1], &i__2, &a[j + 1 + j * a_dim1]);
00921                 }
00922                 i__2 = j + j * a_dim1;
00923                 a[i__2].r = 0., a[i__2].i = 0.;
00924 /* L370: */
00925             }
00926         }
00927 
00928 /*        Set the right hand side so that the largest value is BIGNUM. */
00929 
00930         zlarnv_(&c__2, &iseed[1], n, &b[1]);
00931         iy = izamax_(n, &b[1], &c__1);
00932         bnorm = z_abs(&b[iy]);
00933         bscal = bignum / max(1.,bnorm);
00934         zdscal_(n, &bscal, &b[1], &c__1);
00935 
00936     } else if (*imat == 19) {
00937 
00938 /*        Type 19:  Generate a triangular matrix with elements between */
00939 /*        BIGNUM/(n-1) and BIGNUM so that at least one of the column */
00940 /*        norms will exceed BIGNUM. */
00941 /*        1/3/91:  ZLATRS no longer can handle this case */
00942 
00943 /* Computing MAX */
00944         d__1 = 1., d__2 = (doublereal) (*n - 1);
00945         tleft = bignum / max(d__1,d__2);
00946 /* Computing MAX */
00947         d__1 = 1., d__2 = (doublereal) (*n);
00948         tscal = bignum * ((doublereal) (*n - 1) / max(d__1,d__2));
00949         if (upper) {
00950             i__1 = *n;
00951             for (j = 1; j <= i__1; ++j) {
00952                 zlarnv_(&c__5, &iseed[1], &j, &a[j * a_dim1 + 1]);
00953                 dlarnv_(&c__1, &iseed[1], &j, &rwork[1]);
00954                 i__2 = j;
00955                 for (i__ = 1; i__ <= i__2; ++i__) {
00956                     i__3 = i__ + j * a_dim1;
00957                     i__4 = i__ + j * a_dim1;
00958                     d__1 = tleft + rwork[i__] * tscal;
00959                     z__1.r = d__1 * a[i__4].r, z__1.i = d__1 * a[i__4].i;
00960                     a[i__3].r = z__1.r, a[i__3].i = z__1.i;
00961 /* L380: */
00962                 }
00963 /* L390: */
00964             }
00965         } else {
00966             i__1 = *n;
00967             for (j = 1; j <= i__1; ++j) {
00968                 i__2 = *n - j + 1;
00969                 zlarnv_(&c__5, &iseed[1], &i__2, &a[j + j * a_dim1]);
00970                 i__2 = *n - j + 1;
00971                 dlarnv_(&c__1, &iseed[1], &i__2, &rwork[1]);
00972                 i__2 = *n;
00973                 for (i__ = j; i__ <= i__2; ++i__) {
00974                     i__3 = i__ + j * a_dim1;
00975                     i__4 = i__ + j * a_dim1;
00976                     d__1 = tleft + rwork[i__ - j + 1] * tscal;
00977                     z__1.r = d__1 * a[i__4].r, z__1.i = d__1 * a[i__4].i;
00978                     a[i__3].r = z__1.r, a[i__3].i = z__1.i;
00979 /* L400: */
00980                 }
00981 /* L410: */
00982             }
00983         }
00984         zlarnv_(&c__2, &iseed[1], n, &b[1]);
00985         zdscal_(n, &c_b92, &b[1], &c__1);
00986     }
00987 
00988 /*     Flip the matrix if the transpose will be used. */
00989 
00990     if (! lsame_(trans, "N")) {
00991         if (upper) {
00992             i__1 = *n / 2;
00993             for (j = 1; j <= i__1; ++j) {
00994                 i__2 = *n - (j << 1) + 1;
00995                 zswap_(&i__2, &a[j + j * a_dim1], lda, &a[j + 1 + (*n - j + 1)
00996                          * a_dim1], &c_n1);
00997 /* L420: */
00998             }
00999         } else {
01000             i__1 = *n / 2;
01001             for (j = 1; j <= i__1; ++j) {
01002                 i__2 = *n - (j << 1) + 1;
01003                 i__3 = -(*lda);
01004                 zswap_(&i__2, &a[j + j * a_dim1], &c__1, &a[*n - j + 1 + (j + 
01005                         1) * a_dim1], &i__3);
01006 /* L430: */
01007             }
01008         }
01009     }
01010 
01011     return 0;
01012 
01013 /*     End of ZLATTR */
01014 
01015 } /* zlattr_ */


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