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


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