slatms.c
Go to the documentation of this file.
00001 /* slatms.f -- translated by f2c (version 20061008).
00002    You must link the resulting object file with libf2c:
00003         on Microsoft Windows system, link with libf2c.lib;
00004         on Linux or Unix systems, link with .../path/to/libf2c.a -lm
00005         or, if you install libf2c.a in a standard place, with -lf2c -lm
00006         -- in that order, at the end of the command line, as in
00007                 cc *.o -lf2c -lm
00008         Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
00009 
00010                 http://www.netlib.org/f2c/libf2c.zip
00011 */
00012 
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015 
00016 /* Table of constant values */
00017 
00018 static integer c__1 = 1;
00019 static real c_b22 = 0.f;
00020 static logical c_true = TRUE_;
00021 static logical c_false = FALSE_;
00022 
00023 /* Subroutine */ int slatms_(integer *m, integer *n, char *dist, integer *
00024         iseed, char *sym, real *d__, integer *mode, real *cond, real *dmax__, 
00025         integer *kl, integer *ku, char *pack, real *a, integer *lda, real *
00026         work, integer *info)
00027 {
00028     /* System generated locals */
00029     integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5, i__6;
00030     real r__1, r__2, r__3;
00031     logical L__1;
00032 
00033     /* Builtin functions */
00034     double cos(doublereal), sin(doublereal);
00035 
00036     /* Local variables */
00037     real c__;
00038     integer i__, j, k;
00039     real s;
00040     integer ic, jc, nc, il, ir, jr, mr, ir1, ir2, jch, llb, jkl, jku, uub, 
00041             ilda, icol;
00042     real temp;
00043     integer irow, isym;
00044     real alpha, angle;
00045     integer ipack, ioffg;
00046     extern logical lsame_(char *, char *);
00047     integer iinfo;
00048     extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);
00049     integer idist, mnmin, iskew;
00050     real extra, dummy;
00051     extern /* Subroutine */ int scopy_(integer *, real *, integer *, real *, 
00052             integer *), slatm1_(integer *, real *, integer *, integer *, 
00053             integer *, real *, integer *, integer *);
00054     integer iendch, ipackg;
00055     extern /* Subroutine */ int slagge_(integer *, integer *, integer *, 
00056             integer *, real *, real *, integer *, integer *, real *, integer *
00057 );
00058     integer minlda;
00059     extern /* Subroutine */ int xerbla_(char *, integer *);
00060     extern doublereal slarnd_(integer *, integer *);
00061     logical iltemp, givens;
00062     integer ioffst, irsign;
00063     extern /* Subroutine */ int slartg_(real *, real *, real *, real *, real *
00064 ), slaset_(char *, integer *, integer *, real *, real *, real *, 
00065             integer *), slagsy_(integer *, integer *, real *, real *, 
00066             integer *, integer *, real *, integer *), slarot_(logical *, 
00067             logical *, logical *, integer *, real *, real *, real *, integer *
00068 , real *, real *);
00069     logical ilextr, topdwn;
00070     integer isympk;
00071 
00072 
00073 /*  -- LAPACK test routine (version 3.1) -- */
00074 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00075 /*     November 2006 */
00076 
00077 /*     .. Scalar Arguments .. */
00078 /*     .. */
00079 /*     .. Array Arguments .. */
00080 /*     .. */
00081 
00082 /*  Purpose */
00083 /*  ======= */
00084 
00085 /*     SLATMS generates random matrices with specified singular values */
00086 /*     (or symmetric/hermitian with specified eigenvalues) */
00087 /*     for testing LAPACK programs. */
00088 
00089 /*     SLATMS operates by applying the following sequence of */
00090 /*     operations: */
00091 
00092 /*       Set the diagonal to D, where D may be input or */
00093 /*          computed according to MODE, COND, DMAX, and SYM */
00094 /*          as described below. */
00095 
00096 /*       Generate a matrix with the appropriate band structure, by one */
00097 /*          of two methods: */
00098 
00099 /*       Method A: */
00100 /*           Generate a dense M x N matrix by multiplying D on the left */
00101 /*               and the right by random unitary matrices, then: */
00102 
00103 /*           Reduce the bandwidth according to KL and KU, using */
00104 /*           Householder transformations. */
00105 
00106 /*       Method B: */
00107 /*           Convert the bandwidth-0 (i.e., diagonal) matrix to a */
00108 /*               bandwidth-1 matrix using Givens rotations, "chasing" */
00109 /*               out-of-band elements back, much as in QR; then */
00110 /*               convert the bandwidth-1 to a bandwidth-2 matrix, etc. */
00111 /*               Note that for reasonably small bandwidths (relative to */
00112 /*               M and N) this requires less storage, as a dense matrix */
00113 /*               is not generated.  Also, for symmetric matrices, only */
00114 /*               one triangle is generated. */
00115 
00116 /*       Method A is chosen if the bandwidth is a large fraction of the */
00117 /*           order of the matrix, and LDA is at least M (so a dense */
00118 /*           matrix can be stored.)  Method B is chosen if the bandwidth */
00119 /*           is small (< 1/2 N for symmetric, < .3 N+M for */
00120 /*           non-symmetric), or LDA is less than M and not less than the */
00121 /*           bandwidth. */
00122 
00123 /*       Pack the matrix if desired. Options specified by PACK are: */
00124 /*          no packing */
00125 /*          zero out upper half (if symmetric) */
00126 /*          zero out lower half (if symmetric) */
00127 /*          store the upper half columnwise (if symmetric or upper */
00128 /*                triangular) */
00129 /*          store the lower half columnwise (if symmetric or lower */
00130 /*                triangular) */
00131 /*          store the lower triangle in banded format (if symmetric */
00132 /*                or lower triangular) */
00133 /*          store the upper triangle in banded format (if symmetric */
00134 /*                or upper triangular) */
00135 /*          store the entire matrix in banded format */
00136 /*       If Method B is chosen, and band format is specified, then the */
00137 /*          matrix will be generated in the band format, so no repacking */
00138 /*          will be necessary. */
00139 
00140 /*  Arguments */
00141 /*  ========= */
00142 
00143 /*  M      - INTEGER */
00144 /*           The number of rows of A. Not modified. */
00145 
00146 /*  N      - INTEGER */
00147 /*           The number of columns of A. Not modified. */
00148 
00149 /*  DIST   - CHARACTER*1 */
00150 /*           On entry, DIST specifies the type of distribution to be used */
00151 /*           to generate the random eigen-/singular values. */
00152 /*           'U' => UNIFORM( 0, 1 )  ( 'U' for uniform ) */
00153 /*           'S' => UNIFORM( -1, 1 ) ( 'S' for symmetric ) */
00154 /*           'N' => NORMAL( 0, 1 )   ( 'N' for normal ) */
00155 /*           Not modified. */
00156 
00157 /*  ISEED  - INTEGER array, dimension ( 4 ) */
00158 /*           On entry ISEED specifies the seed of the random number */
00159 /*           generator. They should lie between 0 and 4095 inclusive, */
00160 /*           and ISEED(4) should be odd. The random number generator */
00161 /*           uses a linear congruential sequence limited to small */
00162 /*           integers, and so should produce machine independent */
00163 /*           random numbers. The values of ISEED are changed on */
00164 /*           exit, and can be used in the next call to SLATMS */
00165 /*           to continue the same random number sequence. */
00166 /*           Changed on exit. */
00167 
00168 /*  SYM    - CHARACTER*1 */
00169 /*           If SYM='S' or 'H', the generated matrix is symmetric, with */
00170 /*             eigenvalues specified by D, COND, MODE, and DMAX; they */
00171 /*             may be positive, negative, or zero. */
00172 /*           If SYM='P', the generated matrix is symmetric, with */
00173 /*             eigenvalues (= singular values) specified by D, COND, */
00174 /*             MODE, and DMAX; they will not be negative. */
00175 /*           If SYM='N', the generated matrix is nonsymmetric, with */
00176 /*             singular values specified by D, COND, MODE, and DMAX; */
00177 /*             they will not be negative. */
00178 /*           Not modified. */
00179 
00180 /*  D      - REAL array, dimension ( MIN( M , N ) ) */
00181 /*           This array is used to specify the singular values or */
00182 /*           eigenvalues of A (see SYM, above.)  If MODE=0, then D is */
00183 /*           assumed to contain the singular/eigenvalues, otherwise */
00184 /*           they will be computed according to MODE, COND, and DMAX, */
00185 /*           and placed in D. */
00186 /*           Modified if MODE is nonzero. */
00187 
00188 /*  MODE   - INTEGER */
00189 /*           On entry this describes how the singular/eigenvalues are to */
00190 /*           be specified: */
00191 /*           MODE = 0 means use D as input */
00192 /*           MODE = 1 sets D(1)=1 and D(2:N)=1.0/COND */
00193 /*           MODE = 2 sets D(1:N-1)=1 and D(N)=1.0/COND */
00194 /*           MODE = 3 sets D(I)=COND**(-(I-1)/(N-1)) */
00195 /*           MODE = 4 sets D(i)=1 - (i-1)/(N-1)*(1 - 1/COND) */
00196 /*           MODE = 5 sets D to random numbers in the range */
00197 /*                    ( 1/COND , 1 ) such that their logarithms */
00198 /*                    are uniformly distributed. */
00199 /*           MODE = 6 set D to random numbers from same distribution */
00200 /*                    as the rest of the matrix. */
00201 /*           MODE < 0 has the same meaning as ABS(MODE), except that */
00202 /*              the order of the elements of D is reversed. */
00203 /*           Thus if MODE is positive, D has entries ranging from */
00204 /*              1 to 1/COND, if negative, from 1/COND to 1, */
00205 /*           If SYM='S' or 'H', and MODE is neither 0, 6, nor -6, then */
00206 /*              the elements of D will also be multiplied by a random */
00207 /*              sign (i.e., +1 or -1.) */
00208 /*           Not modified. */
00209 
00210 /*  COND   - REAL */
00211 /*           On entry, this is used as described under MODE above. */
00212 /*           If used, it must be >= 1. Not modified. */
00213 
00214 /*  DMAX   - REAL */
00215 /*           If MODE is neither -6, 0 nor 6, the contents of D, as */
00216 /*           computed according to MODE and COND, will be scaled by */
00217 /*           DMAX / max(abs(D(i))); thus, the maximum absolute eigen- or */
00218 /*           singular value (which is to say the norm) will be abs(DMAX). */
00219 /*           Note that DMAX need not be positive: if DMAX is negative */
00220 /*           (or zero), D will be scaled by a negative number (or zero). */
00221 /*           Not modified. */
00222 
00223 /*  KL     - INTEGER */
00224 /*           This specifies the lower bandwidth of the  matrix. For */
00225 /*           example, KL=0 implies upper triangular, KL=1 implies upper */
00226 /*           Hessenberg, and KL being at least M-1 means that the matrix */
00227 /*           has full lower bandwidth.  KL must equal KU if the matrix */
00228 /*           is symmetric. */
00229 /*           Not modified. */
00230 
00231 /*  KU     - INTEGER */
00232 /*           This specifies the upper bandwidth of the  matrix. For */
00233 /*           example, KU=0 implies lower triangular, KU=1 implies lower */
00234 /*           Hessenberg, and KU being at least N-1 means that the matrix */
00235 /*           has full upper bandwidth.  KL must equal KU if the matrix */
00236 /*           is symmetric. */
00237 /*           Not modified. */
00238 
00239 /*  PACK   - CHARACTER*1 */
00240 /*           This specifies packing of matrix as follows: */
00241 /*           'N' => no packing */
00242 /*           'U' => zero out all subdiagonal entries (if symmetric) */
00243 /*           'L' => zero out all superdiagonal entries (if symmetric) */
00244 /*           'C' => store the upper triangle columnwise */
00245 /*                  (only if the matrix is symmetric or upper triangular) */
00246 /*           'R' => store the lower triangle columnwise */
00247 /*                  (only if the matrix is symmetric or lower triangular) */
00248 /*           'B' => store the lower triangle in band storage scheme */
00249 /*                  (only if matrix symmetric or lower triangular) */
00250 /*           'Q' => store the upper triangle in band storage scheme */
00251 /*                  (only if matrix symmetric or upper triangular) */
00252 /*           'Z' => store the entire matrix in band storage scheme */
00253 /*                      (pivoting can be provided for by using this */
00254 /*                      option to store A in the trailing rows of */
00255 /*                      the allocated storage) */
00256 
00257 /*           Using these options, the various LAPACK packed and banded */
00258 /*           storage schemes can be obtained: */
00259 /*           GB               - use 'Z' */
00260 /*           PB, SB or TB     - use 'B' or 'Q' */
00261 /*           PP, SP or TP     - use 'C' or 'R' */
00262 
00263 /*           If two calls to SLATMS differ only in the PACK parameter, */
00264 /*           they will generate mathematically equivalent matrices. */
00265 /*           Not modified. */
00266 
00267 /*  A      - REAL array, dimension ( LDA, N ) */
00268 /*           On exit A is the desired test matrix.  A is first generated */
00269 /*           in full (unpacked) form, and then packed, if so specified */
00270 /*           by PACK.  Thus, the first M elements of the first N */
00271 /*           columns will always be modified.  If PACK specifies a */
00272 /*           packed or banded storage scheme, all LDA elements of the */
00273 /*           first N columns will be modified; the elements of the */
00274 /*           array which do not correspond to elements of the generated */
00275 /*           matrix are set to zero. */
00276 /*           Modified. */
00277 
00278 /*  LDA    - INTEGER */
00279 /*           LDA specifies the first dimension of A as declared in the */
00280 /*           calling program.  If PACK='N', 'U', 'L', 'C', or 'R', then */
00281 /*           LDA must be at least M.  If PACK='B' or 'Q', then LDA must */
00282 /*           be at least MIN( KL, M-1) (which is equal to MIN(KU,N-1)). */
00283 /*           If PACK='Z', LDA must be large enough to hold the packed */
00284 /*           array: MIN( KU, N-1) + MIN( KL, M-1) + 1. */
00285 /*           Not modified. */
00286 
00287 /*  WORK   - REAL array, dimension ( 3*MAX( N , M ) ) */
00288 /*           Workspace. */
00289 /*           Modified. */
00290 
00291 /*  INFO   - INTEGER */
00292 /*           Error code.  On exit, INFO will be set to one of the */
00293 /*           following values: */
00294 /*             0 => normal return */
00295 /*            -1 => M negative or unequal to N and SYM='S', 'H', or 'P' */
00296 /*            -2 => N negative */
00297 /*            -3 => DIST illegal string */
00298 /*            -5 => SYM illegal string */
00299 /*            -7 => MODE not in range -6 to 6 */
00300 /*            -8 => COND less than 1.0, and MODE neither -6, 0 nor 6 */
00301 /*           -10 => KL negative */
00302 /*           -11 => KU negative, or SYM='S' or 'H' and KU not equal to KL */
00303 /*           -12 => PACK illegal string, or PACK='U' or 'L', and SYM='N'; */
00304 /*                  or PACK='C' or 'Q' and SYM='N' and KL is not zero; */
00305 /*                  or PACK='R' or 'B' and SYM='N' and KU is not zero; */
00306 /*                  or PACK='U', 'L', 'C', 'R', 'B', or 'Q', and M is not */
00307 /*                  N. */
00308 /*           -14 => LDA is less than M, or PACK='Z' and LDA is less than */
00309 /*                  MIN(KU,N-1) + MIN(KL,M-1) + 1. */
00310 /*            1  => Error return from SLATM1 */
00311 /*            2  => Cannot scale to DMAX (max. sing. value is 0) */
00312 /*            3  => Error return from SLAGGE or SLAGSY */
00313 
00314 /*  ===================================================================== */
00315 
00316 /*     .. Parameters .. */
00317 /*     .. */
00318 /*     .. Local Scalars .. */
00319 /*     .. */
00320 /*     .. External Functions .. */
00321 /*     .. */
00322 /*     .. External Subroutines .. */
00323 /*     .. */
00324 /*     .. Intrinsic Functions .. */
00325 /*     .. */
00326 /*     .. Executable Statements .. */
00327 
00328 /*     1)      Decode and Test the input parameters. */
00329 /*             Initialize flags & seed. */
00330 
00331     /* Parameter adjustments */
00332     --iseed;
00333     --d__;
00334     a_dim1 = *lda;
00335     a_offset = 1 + a_dim1;
00336     a -= a_offset;
00337     --work;
00338 
00339     /* Function Body */
00340     *info = 0;
00341 
00342 /*     Quick return if possible */
00343 
00344     if (*m == 0 || *n == 0) {
00345         return 0;
00346     }
00347 
00348 /*     Decode DIST */
00349 
00350     if (lsame_(dist, "U")) {
00351         idist = 1;
00352     } else if (lsame_(dist, "S")) {
00353         idist = 2;
00354     } else if (lsame_(dist, "N")) {
00355         idist = 3;
00356     } else {
00357         idist = -1;
00358     }
00359 
00360 /*     Decode SYM */
00361 
00362     if (lsame_(sym, "N")) {
00363         isym = 1;
00364         irsign = 0;
00365     } else if (lsame_(sym, "P")) {
00366         isym = 2;
00367         irsign = 0;
00368     } else if (lsame_(sym, "S")) {
00369         isym = 2;
00370         irsign = 1;
00371     } else if (lsame_(sym, "H")) {
00372         isym = 2;
00373         irsign = 1;
00374     } else {
00375         isym = -1;
00376     }
00377 
00378 /*     Decode PACK */
00379 
00380     isympk = 0;
00381     if (lsame_(pack, "N")) {
00382         ipack = 0;
00383     } else if (lsame_(pack, "U")) {
00384         ipack = 1;
00385         isympk = 1;
00386     } else if (lsame_(pack, "L")) {
00387         ipack = 2;
00388         isympk = 1;
00389     } else if (lsame_(pack, "C")) {
00390         ipack = 3;
00391         isympk = 2;
00392     } else if (lsame_(pack, "R")) {
00393         ipack = 4;
00394         isympk = 3;
00395     } else if (lsame_(pack, "B")) {
00396         ipack = 5;
00397         isympk = 3;
00398     } else if (lsame_(pack, "Q")) {
00399         ipack = 6;
00400         isympk = 2;
00401     } else if (lsame_(pack, "Z")) {
00402         ipack = 7;
00403     } else {
00404         ipack = -1;
00405     }
00406 
00407 /*     Set certain internal parameters */
00408 
00409     mnmin = min(*m,*n);
00410 /* Computing MIN */
00411     i__1 = *kl, i__2 = *m - 1;
00412     llb = min(i__1,i__2);
00413 /* Computing MIN */
00414     i__1 = *ku, i__2 = *n - 1;
00415     uub = min(i__1,i__2);
00416 /* Computing MIN */
00417     i__1 = *m, i__2 = *n + llb;
00418     mr = min(i__1,i__2);
00419 /* Computing MIN */
00420     i__1 = *n, i__2 = *m + uub;
00421     nc = min(i__1,i__2);
00422 
00423     if (ipack == 5 || ipack == 6) {
00424         minlda = uub + 1;
00425     } else if (ipack == 7) {
00426         minlda = llb + uub + 1;
00427     } else {
00428         minlda = *m;
00429     }
00430 
00431 /*     Use Givens rotation method if bandwidth small enough, */
00432 /*     or if LDA is too small to store the matrix unpacked. */
00433 
00434     givens = FALSE_;
00435     if (isym == 1) {
00436 /* Computing MAX */
00437         i__1 = 1, i__2 = mr + nc;
00438         if ((real) (llb + uub) < (real) max(i__1,i__2) * .3f) {
00439             givens = TRUE_;
00440         }
00441     } else {
00442         if (llb << 1 < *m) {
00443             givens = TRUE_;
00444         }
00445     }
00446     if (*lda < *m && *lda >= minlda) {
00447         givens = TRUE_;
00448     }
00449 
00450 /*     Set INFO if an error */
00451 
00452     if (*m < 0) {
00453         *info = -1;
00454     } else if (*m != *n && isym != 1) {
00455         *info = -1;
00456     } else if (*n < 0) {
00457         *info = -2;
00458     } else if (idist == -1) {
00459         *info = -3;
00460     } else if (isym == -1) {
00461         *info = -5;
00462     } else if (abs(*mode) > 6) {
00463         *info = -7;
00464     } else if (*mode != 0 && abs(*mode) != 6 && *cond < 1.f) {
00465         *info = -8;
00466     } else if (*kl < 0) {
00467         *info = -10;
00468     } else if (*ku < 0 || isym != 1 && *kl != *ku) {
00469         *info = -11;
00470     } else if (ipack == -1 || isympk == 1 && isym == 1 || isympk == 2 && isym 
00471             == 1 && *kl > 0 || isympk == 3 && isym == 1 && *ku > 0 || isympk 
00472             != 0 && *m != *n) {
00473         *info = -12;
00474     } else if (*lda < max(1,minlda)) {
00475         *info = -14;
00476     }
00477 
00478     if (*info != 0) {
00479         i__1 = -(*info);
00480         xerbla_("SLATMS", &i__1);
00481         return 0;
00482     }
00483 
00484 /*     Initialize random number generator */
00485 
00486     for (i__ = 1; i__ <= 4; ++i__) {
00487         iseed[i__] = (i__1 = iseed[i__], abs(i__1)) % 4096;
00488 /* L10: */
00489     }
00490 
00491     if (iseed[4] % 2 != 1) {
00492         ++iseed[4];
00493     }
00494 
00495 /*     2)      Set up D  if indicated. */
00496 
00497 /*             Compute D according to COND and MODE */
00498 
00499     slatm1_(mode, cond, &irsign, &idist, &iseed[1], &d__[1], &mnmin, &iinfo);
00500     if (iinfo != 0) {
00501         *info = 1;
00502         return 0;
00503     }
00504 
00505 /*     Choose Top-Down if D is (apparently) increasing, */
00506 /*     Bottom-Up if D is (apparently) decreasing. */
00507 
00508     if (dabs(d__[1]) <= (r__1 = d__[mnmin], dabs(r__1))) {
00509         topdwn = TRUE_;
00510     } else {
00511         topdwn = FALSE_;
00512     }
00513 
00514     if (*mode != 0 && abs(*mode) != 6) {
00515 
00516 /*        Scale by DMAX */
00517 
00518         temp = dabs(d__[1]);
00519         i__1 = mnmin;
00520         for (i__ = 2; i__ <= i__1; ++i__) {
00521 /* Computing MAX */
00522             r__2 = temp, r__3 = (r__1 = d__[i__], dabs(r__1));
00523             temp = dmax(r__2,r__3);
00524 /* L20: */
00525         }
00526 
00527         if (temp > 0.f) {
00528             alpha = *dmax__ / temp;
00529         } else {
00530             *info = 2;
00531             return 0;
00532         }
00533 
00534         sscal_(&mnmin, &alpha, &d__[1], &c__1);
00535 
00536     }
00537 
00538 /*     3)      Generate Banded Matrix using Givens rotations. */
00539 /*             Also the special case of UUB=LLB=0 */
00540 
00541 /*               Compute Addressing constants to cover all */
00542 /*               storage formats.  Whether GE, SY, GB, or SB, */
00543 /*               upper or lower triangle or both, */
00544 /*               the (i,j)-th element is in */
00545 /*               A( i - ISKEW*j + IOFFST, j ) */
00546 
00547     if (ipack > 4) {
00548         ilda = *lda - 1;
00549         iskew = 1;
00550         if (ipack > 5) {
00551             ioffst = uub + 1;
00552         } else {
00553             ioffst = 1;
00554         }
00555     } else {
00556         ilda = *lda;
00557         iskew = 0;
00558         ioffst = 0;
00559     }
00560 
00561 /*     IPACKG is the format that the matrix is generated in. If this is */
00562 /*     different from IPACK, then the matrix must be repacked at the */
00563 /*     end.  It also signals how to compute the norm, for scaling. */
00564 
00565     ipackg = 0;
00566     slaset_("Full", lda, n, &c_b22, &c_b22, &a[a_offset], lda);
00567 
00568 /*     Diagonal Matrix -- We are done, unless it */
00569 /*     is to be stored SP/PP/TP (PACK='R' or 'C') */
00570 
00571     if (llb == 0 && uub == 0) {
00572         i__1 = ilda + 1;
00573         scopy_(&mnmin, &d__[1], &c__1, &a[1 - iskew + ioffst + a_dim1], &i__1)
00574                 ;
00575         if (ipack <= 2 || ipack >= 5) {
00576             ipackg = ipack;
00577         }
00578 
00579     } else if (givens) {
00580 
00581 /*        Check whether to use Givens rotations, */
00582 /*        Householder transformations, or nothing. */
00583 
00584         if (isym == 1) {
00585 
00586 /*           Non-symmetric -- A = U D V */
00587 
00588             if (ipack > 4) {
00589                 ipackg = ipack;
00590             } else {
00591                 ipackg = 0;
00592             }
00593 
00594             i__1 = ilda + 1;
00595             scopy_(&mnmin, &d__[1], &c__1, &a[1 - iskew + ioffst + a_dim1], &
00596                     i__1);
00597 
00598             if (topdwn) {
00599                 jkl = 0;
00600                 i__1 = uub;
00601                 for (jku = 1; jku <= i__1; ++jku) {
00602 
00603 /*                 Transform from bandwidth JKL, JKU-1 to JKL, JKU */
00604 
00605 /*                 Last row actually rotated is M */
00606 /*                 Last column actually rotated is MIN( M+JKU, N ) */
00607 
00608 /* Computing MIN */
00609                     i__3 = *m + jku;
00610                     i__2 = min(i__3,*n) + jkl - 1;
00611                     for (jr = 1; jr <= i__2; ++jr) {
00612                         extra = 0.f;
00613                         angle = slarnd_(&c__1, &iseed[1]) * 
00614                                 6.2831853071795864769252867663f;
00615                         c__ = cos(angle);
00616                         s = sin(angle);
00617 /* Computing MAX */
00618                         i__3 = 1, i__4 = jr - jkl;
00619                         icol = max(i__3,i__4);
00620                         if (jr < *m) {
00621 /* Computing MIN */
00622                             i__3 = *n, i__4 = jr + jku;
00623                             il = min(i__3,i__4) + 1 - icol;
00624                             L__1 = jr > jkl;
00625                             slarot_(&c_true, &L__1, &c_false, &il, &c__, &s, &
00626                                     a[jr - iskew * icol + ioffst + icol * 
00627                                     a_dim1], &ilda, &extra, &dummy);
00628                         }
00629 
00630 /*                    Chase "EXTRA" back up */
00631 
00632                         ir = jr;
00633                         ic = icol;
00634                         i__3 = -jkl - jku;
00635                         for (jch = jr - jkl; i__3 < 0 ? jch >= 1 : jch <= 1; 
00636                                 jch += i__3) {
00637                             if (ir < *m) {
00638                                 slartg_(&a[ir + 1 - iskew * (ic + 1) + ioffst 
00639                                         + (ic + 1) * a_dim1], &extra, &c__, &
00640                                         s, &dummy);
00641                             }
00642 /* Computing MAX */
00643                             i__4 = 1, i__5 = jch - jku;
00644                             irow = max(i__4,i__5);
00645                             il = ir + 2 - irow;
00646                             temp = 0.f;
00647                             iltemp = jch > jku;
00648                             r__1 = -s;
00649                             slarot_(&c_false, &iltemp, &c_true, &il, &c__, &
00650                                     r__1, &a[irow - iskew * ic + ioffst + ic *
00651                                      a_dim1], &ilda, &temp, &extra);
00652                             if (iltemp) {
00653                                 slartg_(&a[irow + 1 - iskew * (ic + 1) + 
00654                                         ioffst + (ic + 1) * a_dim1], &temp, &
00655                                         c__, &s, &dummy);
00656 /* Computing MAX */
00657                                 i__4 = 1, i__5 = jch - jku - jkl;
00658                                 icol = max(i__4,i__5);
00659                                 il = ic + 2 - icol;
00660                                 extra = 0.f;
00661                                 L__1 = jch > jku + jkl;
00662                                 r__1 = -s;
00663                                 slarot_(&c_true, &L__1, &c_true, &il, &c__, &
00664                                         r__1, &a[irow - iskew * icol + ioffst 
00665                                         + icol * a_dim1], &ilda, &extra, &
00666                                         temp);
00667                                 ic = icol;
00668                                 ir = irow;
00669                             }
00670 /* L30: */
00671                         }
00672 /* L40: */
00673                     }
00674 /* L50: */
00675                 }
00676 
00677                 jku = uub;
00678                 i__1 = llb;
00679                 for (jkl = 1; jkl <= i__1; ++jkl) {
00680 
00681 /*                 Transform from bandwidth JKL-1, JKU to JKL, JKU */
00682 
00683 /* Computing MIN */
00684                     i__3 = *n + jkl;
00685                     i__2 = min(i__3,*m) + jku - 1;
00686                     for (jc = 1; jc <= i__2; ++jc) {
00687                         extra = 0.f;
00688                         angle = slarnd_(&c__1, &iseed[1]) * 
00689                                 6.2831853071795864769252867663f;
00690                         c__ = cos(angle);
00691                         s = sin(angle);
00692 /* Computing MAX */
00693                         i__3 = 1, i__4 = jc - jku;
00694                         irow = max(i__3,i__4);
00695                         if (jc < *n) {
00696 /* Computing MIN */
00697                             i__3 = *m, i__4 = jc + jkl;
00698                             il = min(i__3,i__4) + 1 - irow;
00699                             L__1 = jc > jku;
00700                             slarot_(&c_false, &L__1, &c_false, &il, &c__, &s, 
00701                                     &a[irow - iskew * jc + ioffst + jc * 
00702                                     a_dim1], &ilda, &extra, &dummy);
00703                         }
00704 
00705 /*                    Chase "EXTRA" back up */
00706 
00707                         ic = jc;
00708                         ir = irow;
00709                         i__3 = -jkl - jku;
00710                         for (jch = jc - jku; i__3 < 0 ? jch >= 1 : jch <= 1; 
00711                                 jch += i__3) {
00712                             if (ic < *n) {
00713                                 slartg_(&a[ir + 1 - iskew * (ic + 1) + ioffst 
00714                                         + (ic + 1) * a_dim1], &extra, &c__, &
00715                                         s, &dummy);
00716                             }
00717 /* Computing MAX */
00718                             i__4 = 1, i__5 = jch - jkl;
00719                             icol = max(i__4,i__5);
00720                             il = ic + 2 - icol;
00721                             temp = 0.f;
00722                             iltemp = jch > jkl;
00723                             r__1 = -s;
00724                             slarot_(&c_true, &iltemp, &c_true, &il, &c__, &
00725                                     r__1, &a[ir - iskew * icol + ioffst + 
00726                                     icol * a_dim1], &ilda, &temp, &extra);
00727                             if (iltemp) {
00728                                 slartg_(&a[ir + 1 - iskew * (icol + 1) + 
00729                                         ioffst + (icol + 1) * a_dim1], &temp, 
00730                                         &c__, &s, &dummy);
00731 /* Computing MAX */
00732                                 i__4 = 1, i__5 = jch - jkl - jku;
00733                                 irow = max(i__4,i__5);
00734                                 il = ir + 2 - irow;
00735                                 extra = 0.f;
00736                                 L__1 = jch > jkl + jku;
00737                                 r__1 = -s;
00738                                 slarot_(&c_false, &L__1, &c_true, &il, &c__, &
00739                                         r__1, &a[irow - iskew * icol + ioffst 
00740                                         + icol * a_dim1], &ilda, &extra, &
00741                                         temp);
00742                                 ic = icol;
00743                                 ir = irow;
00744                             }
00745 /* L60: */
00746                         }
00747 /* L70: */
00748                     }
00749 /* L80: */
00750                 }
00751 
00752             } else {
00753 
00754 /*              Bottom-Up -- Start at the bottom right. */
00755 
00756                 jkl = 0;
00757                 i__1 = uub;
00758                 for (jku = 1; jku <= i__1; ++jku) {
00759 
00760 /*                 Transform from bandwidth JKL, JKU-1 to JKL, JKU */
00761 
00762 /*                 First row actually rotated is M */
00763 /*                 First column actually rotated is MIN( M+JKU, N ) */
00764 
00765 /* Computing MIN */
00766                     i__2 = *m, i__3 = *n + jkl;
00767                     iendch = min(i__2,i__3) - 1;
00768 /* Computing MIN */
00769                     i__2 = *m + jku;
00770                     i__3 = 1 - jkl;
00771                     for (jc = min(i__2,*n) - 1; jc >= i__3; --jc) {
00772                         extra = 0.f;
00773                         angle = slarnd_(&c__1, &iseed[1]) * 
00774                                 6.2831853071795864769252867663f;
00775                         c__ = cos(angle);
00776                         s = sin(angle);
00777 /* Computing MAX */
00778                         i__2 = 1, i__4 = jc - jku + 1;
00779                         irow = max(i__2,i__4);
00780                         if (jc > 0) {
00781 /* Computing MIN */
00782                             i__2 = *m, i__4 = jc + jkl + 1;
00783                             il = min(i__2,i__4) + 1 - irow;
00784                             L__1 = jc + jkl < *m;
00785                             slarot_(&c_false, &c_false, &L__1, &il, &c__, &s, 
00786                                     &a[irow - iskew * jc + ioffst + jc * 
00787                                     a_dim1], &ilda, &dummy, &extra);
00788                         }
00789 
00790 /*                    Chase "EXTRA" back down */
00791 
00792                         ic = jc;
00793                         i__2 = iendch;
00794                         i__4 = jkl + jku;
00795                         for (jch = jc + jkl; i__4 < 0 ? jch >= i__2 : jch <= 
00796                                 i__2; jch += i__4) {
00797                             ilextr = ic > 0;
00798                             if (ilextr) {
00799                                 slartg_(&a[jch - iskew * ic + ioffst + ic * 
00800                                         a_dim1], &extra, &c__, &s, &dummy);
00801                             }
00802                             ic = max(1,ic);
00803 /* Computing MIN */
00804                             i__5 = *n - 1, i__6 = jch + jku;
00805                             icol = min(i__5,i__6);
00806                             iltemp = jch + jku < *n;
00807                             temp = 0.f;
00808                             i__5 = icol + 2 - ic;
00809                             slarot_(&c_true, &ilextr, &iltemp, &i__5, &c__, &
00810                                     s, &a[jch - iskew * ic + ioffst + ic * 
00811                                     a_dim1], &ilda, &extra, &temp);
00812                             if (iltemp) {
00813                                 slartg_(&a[jch - iskew * icol + ioffst + icol 
00814                                         * a_dim1], &temp, &c__, &s, &dummy);
00815 /* Computing MIN */
00816                                 i__5 = iendch, i__6 = jch + jkl + jku;
00817                                 il = min(i__5,i__6) + 2 - jch;
00818                                 extra = 0.f;
00819                                 L__1 = jch + jkl + jku <= iendch;
00820                                 slarot_(&c_false, &c_true, &L__1, &il, &c__, &
00821                                         s, &a[jch - iskew * icol + ioffst + 
00822                                         icol * a_dim1], &ilda, &temp, &extra);
00823                                 ic = icol;
00824                             }
00825 /* L90: */
00826                         }
00827 /* L100: */
00828                     }
00829 /* L110: */
00830                 }
00831 
00832                 jku = uub;
00833                 i__1 = llb;
00834                 for (jkl = 1; jkl <= i__1; ++jkl) {
00835 
00836 /*                 Transform from bandwidth JKL-1, JKU to JKL, JKU */
00837 
00838 /*                 First row actually rotated is MIN( N+JKL, M ) */
00839 /*                 First column actually rotated is N */
00840 
00841 /* Computing MIN */
00842                     i__3 = *n, i__4 = *m + jku;
00843                     iendch = min(i__3,i__4) - 1;
00844 /* Computing MIN */
00845                     i__3 = *n + jkl;
00846                     i__4 = 1 - jku;
00847                     for (jr = min(i__3,*m) - 1; jr >= i__4; --jr) {
00848                         extra = 0.f;
00849                         angle = slarnd_(&c__1, &iseed[1]) * 
00850                                 6.2831853071795864769252867663f;
00851                         c__ = cos(angle);
00852                         s = sin(angle);
00853 /* Computing MAX */
00854                         i__3 = 1, i__2 = jr - jkl + 1;
00855                         icol = max(i__3,i__2);
00856                         if (jr > 0) {
00857 /* Computing MIN */
00858                             i__3 = *n, i__2 = jr + jku + 1;
00859                             il = min(i__3,i__2) + 1 - icol;
00860                             L__1 = jr + jku < *n;
00861                             slarot_(&c_true, &c_false, &L__1, &il, &c__, &s, &
00862                                     a[jr - iskew * icol + ioffst + icol * 
00863                                     a_dim1], &ilda, &dummy, &extra);
00864                         }
00865 
00866 /*                    Chase "EXTRA" back down */
00867 
00868                         ir = jr;
00869                         i__3 = iendch;
00870                         i__2 = jkl + jku;
00871                         for (jch = jr + jku; i__2 < 0 ? jch >= i__3 : jch <= 
00872                                 i__3; jch += i__2) {
00873                             ilextr = ir > 0;
00874                             if (ilextr) {
00875                                 slartg_(&a[ir - iskew * jch + ioffst + jch * 
00876                                         a_dim1], &extra, &c__, &s, &dummy);
00877                             }
00878                             ir = max(1,ir);
00879 /* Computing MIN */
00880                             i__5 = *m - 1, i__6 = jch + jkl;
00881                             irow = min(i__5,i__6);
00882                             iltemp = jch + jkl < *m;
00883                             temp = 0.f;
00884                             i__5 = irow + 2 - ir;
00885                             slarot_(&c_false, &ilextr, &iltemp, &i__5, &c__, &
00886                                     s, &a[ir - iskew * jch + ioffst + jch * 
00887                                     a_dim1], &ilda, &extra, &temp);
00888                             if (iltemp) {
00889                                 slartg_(&a[irow - iskew * jch + ioffst + jch *
00890                                          a_dim1], &temp, &c__, &s, &dummy);
00891 /* Computing MIN */
00892                                 i__5 = iendch, i__6 = jch + jkl + jku;
00893                                 il = min(i__5,i__6) + 2 - jch;
00894                                 extra = 0.f;
00895                                 L__1 = jch + jkl + jku <= iendch;
00896                                 slarot_(&c_true, &c_true, &L__1, &il, &c__, &
00897                                         s, &a[irow - iskew * jch + ioffst + 
00898                                         jch * a_dim1], &ilda, &temp, &extra);
00899                                 ir = irow;
00900                             }
00901 /* L120: */
00902                         }
00903 /* L130: */
00904                     }
00905 /* L140: */
00906                 }
00907             }
00908 
00909         } else {
00910 
00911 /*           Symmetric -- A = U D U' */
00912 
00913             ipackg = ipack;
00914             ioffg = ioffst;
00915 
00916             if (topdwn) {
00917 
00918 /*              Top-Down -- Generate Upper triangle only */
00919 
00920                 if (ipack >= 5) {
00921                     ipackg = 6;
00922                     ioffg = uub + 1;
00923                 } else {
00924                     ipackg = 1;
00925                 }
00926                 i__1 = ilda + 1;
00927                 scopy_(&mnmin, &d__[1], &c__1, &a[1 - iskew + ioffg + a_dim1], 
00928                          &i__1);
00929 
00930                 i__1 = uub;
00931                 for (k = 1; k <= i__1; ++k) {
00932                     i__4 = *n - 1;
00933                     for (jc = 1; jc <= i__4; ++jc) {
00934 /* Computing MAX */
00935                         i__2 = 1, i__3 = jc - k;
00936                         irow = max(i__2,i__3);
00937 /* Computing MIN */
00938                         i__2 = jc + 1, i__3 = k + 2;
00939                         il = min(i__2,i__3);
00940                         extra = 0.f;
00941                         temp = a[jc - iskew * (jc + 1) + ioffg + (jc + 1) * 
00942                                 a_dim1];
00943                         angle = slarnd_(&c__1, &iseed[1]) * 
00944                                 6.2831853071795864769252867663f;
00945                         c__ = cos(angle);
00946                         s = sin(angle);
00947                         L__1 = jc > k;
00948                         slarot_(&c_false, &L__1, &c_true, &il, &c__, &s, &a[
00949                                 irow - iskew * jc + ioffg + jc * a_dim1], &
00950                                 ilda, &extra, &temp);
00951 /* Computing MIN */
00952                         i__3 = k, i__5 = *n - jc;
00953                         i__2 = min(i__3,i__5) + 1;
00954                         slarot_(&c_true, &c_true, &c_false, &i__2, &c__, &s, &
00955                                 a[(1 - iskew) * jc + ioffg + jc * a_dim1], &
00956                                 ilda, &temp, &dummy);
00957 
00958 /*                    Chase EXTRA back up the matrix */
00959 
00960                         icol = jc;
00961                         i__2 = -k;
00962                         for (jch = jc - k; i__2 < 0 ? jch >= 1 : jch <= 1; 
00963                                 jch += i__2) {
00964                             slartg_(&a[jch + 1 - iskew * (icol + 1) + ioffg + 
00965                                     (icol + 1) * a_dim1], &extra, &c__, &s, &
00966                                     dummy);
00967                             temp = a[jch - iskew * (jch + 1) + ioffg + (jch + 
00968                                     1) * a_dim1];
00969                             i__3 = k + 2;
00970                             r__1 = -s;
00971                             slarot_(&c_true, &c_true, &c_true, &i__3, &c__, &
00972                                     r__1, &a[(1 - iskew) * jch + ioffg + jch *
00973                                      a_dim1], &ilda, &temp, &extra);
00974 /* Computing MAX */
00975                             i__3 = 1, i__5 = jch - k;
00976                             irow = max(i__3,i__5);
00977 /* Computing MIN */
00978                             i__3 = jch + 1, i__5 = k + 2;
00979                             il = min(i__3,i__5);
00980                             extra = 0.f;
00981                             L__1 = jch > k;
00982                             r__1 = -s;
00983                             slarot_(&c_false, &L__1, &c_true, &il, &c__, &
00984                                     r__1, &a[irow - iskew * jch + ioffg + jch 
00985                                     * a_dim1], &ilda, &extra, &temp);
00986                             icol = jch;
00987 /* L150: */
00988                         }
00989 /* L160: */
00990                     }
00991 /* L170: */
00992                 }
00993 
00994 /*              If we need lower triangle, copy from upper. Note that */
00995 /*              the order of copying is chosen to work for 'q' -> 'b' */
00996 
00997                 if (ipack != ipackg && ipack != 3) {
00998                     i__1 = *n;
00999                     for (jc = 1; jc <= i__1; ++jc) {
01000                         irow = ioffst - iskew * jc;
01001 /* Computing MIN */
01002                         i__2 = *n, i__3 = jc + uub;
01003                         i__4 = min(i__2,i__3);
01004                         for (jr = jc; jr <= i__4; ++jr) {
01005                             a[jr + irow + jc * a_dim1] = a[jc - iskew * jr + 
01006                                     ioffg + jr * a_dim1];
01007 /* L180: */
01008                         }
01009 /* L190: */
01010                     }
01011                     if (ipack == 5) {
01012                         i__1 = *n;
01013                         for (jc = *n - uub + 1; jc <= i__1; ++jc) {
01014                             i__4 = uub + 1;
01015                             for (jr = *n + 2 - jc; jr <= i__4; ++jr) {
01016                                 a[jr + jc * a_dim1] = 0.f;
01017 /* L200: */
01018                             }
01019 /* L210: */
01020                         }
01021                     }
01022                     if (ipackg == 6) {
01023                         ipackg = ipack;
01024                     } else {
01025                         ipackg = 0;
01026                     }
01027                 }
01028             } else {
01029 
01030 /*              Bottom-Up -- Generate Lower triangle only */
01031 
01032                 if (ipack >= 5) {
01033                     ipackg = 5;
01034                     if (ipack == 6) {
01035                         ioffg = 1;
01036                     }
01037                 } else {
01038                     ipackg = 2;
01039                 }
01040                 i__1 = ilda + 1;
01041                 scopy_(&mnmin, &d__[1], &c__1, &a[1 - iskew + ioffg + a_dim1], 
01042                          &i__1);
01043 
01044                 i__1 = uub;
01045                 for (k = 1; k <= i__1; ++k) {
01046                     for (jc = *n - 1; jc >= 1; --jc) {
01047 /* Computing MIN */
01048                         i__4 = *n + 1 - jc, i__2 = k + 2;
01049                         il = min(i__4,i__2);
01050                         extra = 0.f;
01051                         temp = a[(1 - iskew) * jc + 1 + ioffg + jc * a_dim1];
01052                         angle = slarnd_(&c__1, &iseed[1]) * 
01053                                 6.2831853071795864769252867663f;
01054                         c__ = cos(angle);
01055                         s = -sin(angle);
01056                         L__1 = *n - jc > k;
01057                         slarot_(&c_false, &c_true, &L__1, &il, &c__, &s, &a[(
01058                                 1 - iskew) * jc + ioffg + jc * a_dim1], &ilda, 
01059                                  &temp, &extra);
01060 /* Computing MAX */
01061                         i__4 = 1, i__2 = jc - k + 1;
01062                         icol = max(i__4,i__2);
01063                         i__4 = jc + 2 - icol;
01064                         slarot_(&c_true, &c_false, &c_true, &i__4, &c__, &s, &
01065                                 a[jc - iskew * icol + ioffg + icol * a_dim1], 
01066                                 &ilda, &dummy, &temp);
01067 
01068 /*                    Chase EXTRA back down the matrix */
01069 
01070                         icol = jc;
01071                         i__4 = *n - 1;
01072                         i__2 = k;
01073                         for (jch = jc + k; i__2 < 0 ? jch >= i__4 : jch <= 
01074                                 i__4; jch += i__2) {
01075                             slartg_(&a[jch - iskew * icol + ioffg + icol * 
01076                                     a_dim1], &extra, &c__, &s, &dummy);
01077                             temp = a[(1 - iskew) * jch + 1 + ioffg + jch * 
01078                                     a_dim1];
01079                             i__3 = k + 2;
01080                             slarot_(&c_true, &c_true, &c_true, &i__3, &c__, &
01081                                     s, &a[jch - iskew * icol + ioffg + icol * 
01082                                     a_dim1], &ilda, &extra, &temp);
01083 /* Computing MIN */
01084                             i__3 = *n + 1 - jch, i__5 = k + 2;
01085                             il = min(i__3,i__5);
01086                             extra = 0.f;
01087                             L__1 = *n - jch > k;
01088                             slarot_(&c_false, &c_true, &L__1, &il, &c__, &s, &
01089                                     a[(1 - iskew) * jch + ioffg + jch * 
01090                                     a_dim1], &ilda, &temp, &extra);
01091                             icol = jch;
01092 /* L220: */
01093                         }
01094 /* L230: */
01095                     }
01096 /* L240: */
01097                 }
01098 
01099 /*              If we need upper triangle, copy from lower. Note that */
01100 /*              the order of copying is chosen to work for 'b' -> 'q' */
01101 
01102                 if (ipack != ipackg && ipack != 4) {
01103                     for (jc = *n; jc >= 1; --jc) {
01104                         irow = ioffst - iskew * jc;
01105 /* Computing MAX */
01106                         i__2 = 1, i__4 = jc - uub;
01107                         i__1 = max(i__2,i__4);
01108                         for (jr = jc; jr >= i__1; --jr) {
01109                             a[jr + irow + jc * a_dim1] = a[jc - iskew * jr + 
01110                                     ioffg + jr * a_dim1];
01111 /* L250: */
01112                         }
01113 /* L260: */
01114                     }
01115                     if (ipack == 6) {
01116                         i__1 = uub;
01117                         for (jc = 1; jc <= i__1; ++jc) {
01118                             i__2 = uub + 1 - jc;
01119                             for (jr = 1; jr <= i__2; ++jr) {
01120                                 a[jr + jc * a_dim1] = 0.f;
01121 /* L270: */
01122                             }
01123 /* L280: */
01124                         }
01125                     }
01126                     if (ipackg == 5) {
01127                         ipackg = ipack;
01128                     } else {
01129                         ipackg = 0;
01130                     }
01131                 }
01132             }
01133         }
01134 
01135     } else {
01136 
01137 /*        4)      Generate Banded Matrix by first */
01138 /*                Rotating by random Unitary matrices, */
01139 /*                then reducing the bandwidth using Householder */
01140 /*                transformations. */
01141 
01142 /*                Note: we should get here only if LDA .ge. N */
01143 
01144         if (isym == 1) {
01145 
01146 /*           Non-symmetric -- A = U D V */
01147 
01148             slagge_(&mr, &nc, &llb, &uub, &d__[1], &a[a_offset], lda, &iseed[
01149                     1], &work[1], &iinfo);
01150         } else {
01151 
01152 /*           Symmetric -- A = U D U' */
01153 
01154             slagsy_(m, &llb, &d__[1], &a[a_offset], lda, &iseed[1], &work[1], 
01155                     &iinfo);
01156 
01157         }
01158         if (iinfo != 0) {
01159             *info = 3;
01160             return 0;
01161         }
01162     }
01163 
01164 /*     5)      Pack the matrix */
01165 
01166     if (ipack != ipackg) {
01167         if (ipack == 1) {
01168 
01169 /*           'U' -- Upper triangular, not packed */
01170 
01171             i__1 = *m;
01172             for (j = 1; j <= i__1; ++j) {
01173                 i__2 = *m;
01174                 for (i__ = j + 1; i__ <= i__2; ++i__) {
01175                     a[i__ + j * a_dim1] = 0.f;
01176 /* L290: */
01177                 }
01178 /* L300: */
01179             }
01180 
01181         } else if (ipack == 2) {
01182 
01183 /*           'L' -- Lower triangular, not packed */
01184 
01185             i__1 = *m;
01186             for (j = 2; j <= i__1; ++j) {
01187                 i__2 = j - 1;
01188                 for (i__ = 1; i__ <= i__2; ++i__) {
01189                     a[i__ + j * a_dim1] = 0.f;
01190 /* L310: */
01191                 }
01192 /* L320: */
01193             }
01194 
01195         } else if (ipack == 3) {
01196 
01197 /*           'C' -- Upper triangle packed Columnwise. */
01198 
01199             icol = 1;
01200             irow = 0;
01201             i__1 = *m;
01202             for (j = 1; j <= i__1; ++j) {
01203                 i__2 = j;
01204                 for (i__ = 1; i__ <= i__2; ++i__) {
01205                     ++irow;
01206                     if (irow > *lda) {
01207                         irow = 1;
01208                         ++icol;
01209                     }
01210                     a[irow + icol * a_dim1] = a[i__ + j * a_dim1];
01211 /* L330: */
01212                 }
01213 /* L340: */
01214             }
01215 
01216         } else if (ipack == 4) {
01217 
01218 /*           'R' -- Lower triangle packed Columnwise. */
01219 
01220             icol = 1;
01221             irow = 0;
01222             i__1 = *m;
01223             for (j = 1; j <= i__1; ++j) {
01224                 i__2 = *m;
01225                 for (i__ = j; i__ <= i__2; ++i__) {
01226                     ++irow;
01227                     if (irow > *lda) {
01228                         irow = 1;
01229                         ++icol;
01230                     }
01231                     a[irow + icol * a_dim1] = a[i__ + j * a_dim1];
01232 /* L350: */
01233                 }
01234 /* L360: */
01235             }
01236 
01237         } else if (ipack >= 5) {
01238 
01239 /*           'B' -- The lower triangle is packed as a band matrix. */
01240 /*           'Q' -- The upper triangle is packed as a band matrix. */
01241 /*           'Z' -- The whole matrix is packed as a band matrix. */
01242 
01243             if (ipack == 5) {
01244                 uub = 0;
01245             }
01246             if (ipack == 6) {
01247                 llb = 0;
01248             }
01249 
01250             i__1 = uub;
01251             for (j = 1; j <= i__1; ++j) {
01252 /* Computing MIN */
01253                 i__2 = j + llb;
01254                 for (i__ = min(i__2,*m); i__ >= 1; --i__) {
01255                     a[i__ - j + uub + 1 + j * a_dim1] = a[i__ + j * a_dim1];
01256 /* L370: */
01257                 }
01258 /* L380: */
01259             }
01260 
01261             i__1 = *n;
01262             for (j = uub + 2; j <= i__1; ++j) {
01263 /* Computing MIN */
01264                 i__4 = j + llb;
01265                 i__2 = min(i__4,*m);
01266                 for (i__ = j - uub; i__ <= i__2; ++i__) {
01267                     a[i__ - j + uub + 1 + j * a_dim1] = a[i__ + j * a_dim1];
01268 /* L390: */
01269                 }
01270 /* L400: */
01271             }
01272         }
01273 
01274 /*        If packed, zero out extraneous elements. */
01275 
01276 /*        Symmetric/Triangular Packed -- */
01277 /*        zero out everything after A(IROW,ICOL) */
01278 
01279         if (ipack == 3 || ipack == 4) {
01280             i__1 = *m;
01281             for (jc = icol; jc <= i__1; ++jc) {
01282                 i__2 = *lda;
01283                 for (jr = irow + 1; jr <= i__2; ++jr) {
01284                     a[jr + jc * a_dim1] = 0.f;
01285 /* L410: */
01286                 }
01287                 irow = 0;
01288 /* L420: */
01289             }
01290 
01291         } else if (ipack >= 5) {
01292 
01293 /*           Packed Band -- */
01294 /*              1st row is now in A( UUB+2-j, j), zero above it */
01295 /*              m-th row is now in A( M+UUB-j,j), zero below it */
01296 /*              last non-zero diagonal is now in A( UUB+LLB+1,j ), */
01297 /*                 zero below it, too. */
01298 
01299             ir1 = uub + llb + 2;
01300             ir2 = uub + *m + 2;
01301             i__1 = *n;
01302             for (jc = 1; jc <= i__1; ++jc) {
01303                 i__2 = uub + 1 - jc;
01304                 for (jr = 1; jr <= i__2; ++jr) {
01305                     a[jr + jc * a_dim1] = 0.f;
01306 /* L430: */
01307                 }
01308 /* Computing MAX */
01309 /* Computing MIN */
01310                 i__3 = ir1, i__5 = ir2 - jc;
01311                 i__2 = 1, i__4 = min(i__3,i__5);
01312                 i__6 = *lda;
01313                 for (jr = max(i__2,i__4); jr <= i__6; ++jr) {
01314                     a[jr + jc * a_dim1] = 0.f;
01315 /* L440: */
01316                 }
01317 /* L450: */
01318             }
01319         }
01320     }
01321 
01322     return 0;
01323 
01324 /*     End of SLATMS */
01325 
01326 } /* slatms_ */


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