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


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