slatme.c
Go to the documentation of this file.
00001 /* slatme.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_b23 = 0.f;
00020 static integer c__0 = 0;
00021 static real c_b39 = 1.f;
00022 
00023 /* Subroutine */ int slatme_(integer *n, char *dist, integer *iseed, real *
00024         d__, integer *mode, real *cond, real *dmax__, char *ei, char *rsign, 
00025         char *upper, char *sim, real *ds, integer *modes, real *conds, 
00026         integer *kl, integer *ku, real *anorm, real *a, integer *lda, real *
00027         work, integer *info)
00028 {
00029     /* System generated locals */
00030     integer a_dim1, a_offset, i__1, i__2;
00031     real r__1, r__2, r__3;
00032 
00033     /* Local variables */
00034     integer i__, j, ic, jc, ir, jr, jcr;
00035     real tau;
00036     logical bads;
00037     extern /* Subroutine */ int sger_(integer *, integer *, real *, real *, 
00038             integer *, real *, integer *, real *, integer *);
00039     integer isim;
00040     real temp;
00041     logical badei;
00042     real alpha;
00043     extern logical lsame_(char *, char *);
00044     integer iinfo;
00045     extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);
00046     real tempa[1];
00047     integer icols;
00048     logical useei;
00049     integer idist;
00050     extern /* Subroutine */ int sgemv_(char *, integer *, integer *, real *, 
00051             real *, integer *, real *, integer *, real *, real *, integer *), scopy_(integer *, real *, integer *, real *, integer *);
00052     integer irows;
00053     extern /* Subroutine */ int slatm1_(integer *, real *, integer *, integer 
00054             *, integer *, real *, integer *, integer *);
00055     extern doublereal slange_(char *, integer *, integer *, real *, integer *, 
00056              real *);
00057     extern /* Subroutine */ int slarge_(integer *, real *, integer *, integer 
00058             *, real *, integer *), slarfg_(integer *, real *, real *, integer 
00059             *, real *), xerbla_(char *, integer *);
00060     extern doublereal slaran_(integer *);
00061     integer irsign;
00062     extern /* Subroutine */ int slaset_(char *, integer *, integer *, real *, 
00063             real *, real *, integer *);
00064     integer iupper;
00065     extern /* Subroutine */ int slarnv_(integer *, integer *, integer *, real 
00066             *);
00067     real xnorms;
00068 
00069 
00070 /*  -- LAPACK test routine (version 3.1) -- */
00071 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00072 /*     November 2006 */
00073 
00074 /*     .. Scalar Arguments .. */
00075 /*     .. */
00076 /*     .. Array Arguments .. */
00077 /*     .. */
00078 
00079 /*  Purpose */
00080 /*  ======= */
00081 
00082 /*     SLATME generates random non-symmetric square matrices with */
00083 /*     specified eigenvalues for testing LAPACK programs. */
00084 
00085 /*     SLATME operates by applying the following sequence of */
00086 /*     operations: */
00087 
00088 /*     1. Set the diagonal to D, where D may be input or */
00089 /*          computed according to MODE, COND, DMAX, and RSIGN */
00090 /*          as described below. */
00091 
00092 /*     2. If complex conjugate pairs are desired (MODE=0 and EI(1)='R', */
00093 /*          or MODE=5), certain pairs of adjacent elements of D are */
00094 /*          interpreted as the real and complex parts of a complex */
00095 /*          conjugate pair; A thus becomes block diagonal, with 1x1 */
00096 /*          and 2x2 blocks. */
00097 
00098 /*     3. If UPPER='T', the upper triangle of A is set to random values */
00099 /*          out of distribution DIST. */
00100 
00101 /*     4. If SIM='T', A is multiplied on the left by a random matrix */
00102 /*          X, whose singular values are specified by DS, MODES, and */
00103 /*          CONDS, and on the right by X inverse. */
00104 
00105 /*     5. If KL < N-1, the lower bandwidth is reduced to KL using */
00106 /*          Householder transformations.  If KU < N-1, the upper */
00107 /*          bandwidth is reduced to KU. */
00108 
00109 /*     6. If ANORM is not negative, the matrix is scaled to have */
00110 /*          maximum-element-norm ANORM. */
00111 
00112 /*     (Note: since the matrix cannot be reduced beyond Hessenberg form, */
00113 /*      no packing options are available.) */
00114 
00115 /*  Arguments */
00116 /*  ========= */
00117 
00118 /*  N      - INTEGER */
00119 /*           The number of columns (or rows) of A. Not modified. */
00120 
00121 /*  DIST   - CHARACTER*1 */
00122 /*           On entry, DIST specifies the type of distribution to be used */
00123 /*           to generate the random eigen-/singular values, and for the */
00124 /*           upper triangle (see UPPER). */
00125 /*           'U' => UNIFORM( 0, 1 )  ( 'U' for uniform ) */
00126 /*           'S' => UNIFORM( -1, 1 ) ( 'S' for symmetric ) */
00127 /*           'N' => NORMAL( 0, 1 )   ( 'N' for normal ) */
00128 /*           Not modified. */
00129 
00130 /*  ISEED  - INTEGER array, dimension ( 4 ) */
00131 /*           On entry ISEED specifies the seed of the random number */
00132 /*           generator. They should lie between 0 and 4095 inclusive, */
00133 /*           and ISEED(4) should be odd. The random number generator */
00134 /*           uses a linear congruential sequence limited to small */
00135 /*           integers, and so should produce machine independent */
00136 /*           random numbers. The values of ISEED are changed on */
00137 /*           exit, and can be used in the next call to SLATME */
00138 /*           to continue the same random number sequence. */
00139 /*           Changed on exit. */
00140 
00141 /*  D      - REAL array, dimension ( N ) */
00142 /*           This array is used to specify the eigenvalues of A.  If */
00143 /*           MODE=0, then D is assumed to contain the eigenvalues (but */
00144 /*           see the description of EI), otherwise they will be */
00145 /*           computed according to MODE, COND, DMAX, and RSIGN and */
00146 /*           placed in D. */
00147 /*           Modified if MODE is nonzero. */
00148 
00149 /*  MODE   - INTEGER */
00150 /*           On entry this describes how the eigenvalues are to */
00151 /*           be specified: */
00152 /*           MODE = 0 means use D (with EI) as input */
00153 /*           MODE = 1 sets D(1)=1 and D(2:N)=1.0/COND */
00154 /*           MODE = 2 sets D(1:N-1)=1 and D(N)=1.0/COND */
00155 /*           MODE = 3 sets D(I)=COND**(-(I-1)/(N-1)) */
00156 /*           MODE = 4 sets D(i)=1 - (i-1)/(N-1)*(1 - 1/COND) */
00157 /*           MODE = 5 sets D to random numbers in the range */
00158 /*                    ( 1/COND , 1 ) such that their logarithms */
00159 /*                    are uniformly distributed.  Each odd-even pair */
00160 /*                    of elements will be either used as two real */
00161 /*                    eigenvalues or as the real and imaginary part */
00162 /*                    of a complex conjugate pair of eigenvalues; */
00163 /*                    the choice of which is done is random, with */
00164 /*                    50-50 probability, for each pair. */
00165 /*           MODE = 6 set D to random numbers from same distribution */
00166 /*                    as the rest of the matrix. */
00167 /*           MODE < 0 has the same meaning as ABS(MODE), except that */
00168 /*              the order of the elements of D is reversed. */
00169 /*           Thus if MODE is between 1 and 4, D has entries ranging */
00170 /*              from 1 to 1/COND, if between -1 and -4, D has entries */
00171 /*              ranging from 1/COND to 1, */
00172 /*           Not modified. */
00173 
00174 /*  COND   - REAL */
00175 /*           On entry, this is used as described under MODE above. */
00176 /*           If used, it must be >= 1. Not modified. */
00177 
00178 /*  DMAX   - REAL */
00179 /*           If MODE is neither -6, 0 nor 6, the contents of D, as */
00180 /*           computed according to MODE and COND, will be scaled by */
00181 /*           DMAX / max(abs(D(i))).  Note that DMAX need not be */
00182 /*           positive: if DMAX is negative (or zero), D will be */
00183 /*           scaled by a negative number (or zero). */
00184 /*           Not modified. */
00185 
00186 /*  EI     - CHARACTER*1 array, dimension ( N ) */
00187 /*           If MODE is 0, and EI(1) is not ' ' (space character), */
00188 /*           this array specifies which elements of D (on input) are */
00189 /*           real eigenvalues and which are the real and imaginary parts */
00190 /*           of a complex conjugate pair of eigenvalues.  The elements */
00191 /*           of EI may then only have the values 'R' and 'I'.  If */
00192 /*           EI(j)='R' and EI(j+1)='I', then the j-th eigenvalue is */
00193 /*           CMPLX( D(j) , D(j+1) ), and the (j+1)-th is the complex */
00194 /*           conjugate thereof.  If EI(j)=EI(j+1)='R', then the j-th */
00195 /*           eigenvalue is D(j) (i.e., real).  EI(1) may not be 'I', */
00196 /*           nor may two adjacent elements of EI both have the value 'I'. */
00197 /*           If MODE is not 0, then EI is ignored.  If MODE is 0 and */
00198 /*           EI(1)=' ', then the eigenvalues will all be real. */
00199 /*           Not modified. */
00200 
00201 /*  RSIGN  - CHARACTER*1 */
00202 /*           If MODE is not 0, 6, or -6, and RSIGN='T', then the */
00203 /*           elements of D, as computed according to MODE and COND, will */
00204 /*           be multiplied by a random sign (+1 or -1).  If RSIGN='F', */
00205 /*           they will not be.  RSIGN may only have the values 'T' or */
00206 /*           'F'. */
00207 /*           Not modified. */
00208 
00209 /*  UPPER  - CHARACTER*1 */
00210 /*           If UPPER='T', then the elements of A above the diagonal */
00211 /*           (and above the 2x2 diagonal blocks, if A has complex */
00212 /*           eigenvalues) will be set to random numbers out of DIST. */
00213 /*           If UPPER='F', they will not.  UPPER may only have the */
00214 /*           values 'T' or 'F'. */
00215 /*           Not modified. */
00216 
00217 /*  SIM    - CHARACTER*1 */
00218 /*           If SIM='T', then A will be operated on by a "similarity */
00219 /*           transform", i.e., multiplied on the left by a matrix X and */
00220 /*           on the right by X inverse.  X = U S V, where U and V are */
00221 /*           random unitary matrices and S is a (diagonal) matrix of */
00222 /*           singular values specified by DS, MODES, and CONDS.  If */
00223 /*           SIM='F', then A will not be transformed. */
00224 /*           Not modified. */
00225 
00226 /*  DS     - REAL array, dimension ( N ) */
00227 /*           This array is used to specify the singular values of X, */
00228 /*           in the same way that D specifies the eigenvalues of A. */
00229 /*           If MODE=0, the DS contains the singular values, which */
00230 /*           may not be zero. */
00231 /*           Modified if MODE is nonzero. */
00232 
00233 /*  MODES  - INTEGER */
00234 /*  CONDS  - REAL */
00235 /*           Same as MODE and COND, but for specifying the diagonal */
00236 /*           of S.  MODES=-6 and +6 are not allowed (since they would */
00237 /*           result in randomly ill-conditioned eigenvalues.) */
00238 
00239 /*  KL     - INTEGER */
00240 /*           This specifies the lower bandwidth of the  matrix.  KL=1 */
00241 /*           specifies upper Hessenberg form.  If KL is at least N-1, */
00242 /*           then A will have full lower bandwidth.  KL must be at */
00243 /*           least 1. */
00244 /*           Not modified. */
00245 
00246 /*  KU     - INTEGER */
00247 /*           This specifies the upper bandwidth of the  matrix.  KU=1 */
00248 /*           specifies lower Hessenberg form.  If KU is at least N-1, */
00249 /*           then A will have full upper bandwidth; if KU and KL */
00250 /*           are both at least N-1, then A will be dense.  Only one of */
00251 /*           KU and KL may be less than N-1.  KU must be at least 1. */
00252 /*           Not modified. */
00253 
00254 /*  ANORM  - REAL */
00255 /*           If ANORM is not negative, then A will be scaled by a non- */
00256 /*           negative real number to make the maximum-element-norm of A */
00257 /*           to be ANORM. */
00258 /*           Not modified. */
00259 
00260 /*  A      - REAL array, dimension ( LDA, N ) */
00261 /*           On exit A is the desired test matrix. */
00262 /*           Modified. */
00263 
00264 /*  LDA    - INTEGER */
00265 /*           LDA specifies the first dimension of A as declared in the */
00266 /*           calling program.  LDA must be at least N. */
00267 /*           Not modified. */
00268 
00269 /*  WORK   - REAL array, dimension ( 3*N ) */
00270 /*           Workspace. */
00271 /*           Modified. */
00272 
00273 /*  INFO   - INTEGER */
00274 /*           Error code.  On exit, INFO will be set to one of the */
00275 /*           following values: */
00276 /*             0 => normal return */
00277 /*            -1 => N negative */
00278 /*            -2 => DIST illegal string */
00279 /*            -5 => MODE not in range -6 to 6 */
00280 /*            -6 => COND less than 1.0, and MODE neither -6, 0 nor 6 */
00281 /*            -8 => EI(1) is not ' ' or 'R', EI(j) is not 'R' or 'I', or */
00282 /*                  two adjacent elements of EI are 'I'. */
00283 /*            -9 => RSIGN is not 'T' or 'F' */
00284 /*           -10 => UPPER is not 'T' or 'F' */
00285 /*           -11 => SIM   is not 'T' or 'F' */
00286 /*           -12 => MODES=0 and DS has a zero singular value. */
00287 /*           -13 => MODES is not in the range -5 to 5. */
00288 /*           -14 => MODES is nonzero and CONDS is less than 1. */
00289 /*           -15 => KL is less than 1. */
00290 /*           -16 => KU is less than 1, or KL and KU are both less than */
00291 /*                  N-1. */
00292 /*           -19 => LDA is less than N. */
00293 /*            1  => Error return from SLATM1 (computing D) */
00294 /*            2  => Cannot scale to DMAX (max. eigenvalue is 0) */
00295 /*            3  => Error return from SLATM1 (computing DS) */
00296 /*            4  => Error return from SLARGE */
00297 /*            5  => Zero singular value from SLATM1. */
00298 
00299 /*  ===================================================================== */
00300 
00301 /*     .. Parameters .. */
00302 /*     .. */
00303 /*     .. Local Scalars .. */
00304 /*     .. */
00305 /*     .. Local Arrays .. */
00306 /*     .. */
00307 /*     .. External Functions .. */
00308 /*     .. */
00309 /*     .. External Subroutines .. */
00310 /*     .. */
00311 /*     .. Intrinsic Functions .. */
00312 /*     .. */
00313 /*     .. Executable Statements .. */
00314 
00315 /*     1)      Decode and Test the input parameters. */
00316 /*             Initialize flags & seed. */
00317 
00318     /* Parameter adjustments */
00319     --iseed;
00320     --d__;
00321     --ei;
00322     --ds;
00323     a_dim1 = *lda;
00324     a_offset = 1 + a_dim1;
00325     a -= a_offset;
00326     --work;
00327 
00328     /* Function Body */
00329     *info = 0;
00330 
00331 /*     Quick return if possible */
00332 
00333     if (*n == 0) {
00334         return 0;
00335     }
00336 
00337 /*     Decode DIST */
00338 
00339     if (lsame_(dist, "U")) {
00340         idist = 1;
00341     } else if (lsame_(dist, "S")) {
00342         idist = 2;
00343     } else if (lsame_(dist, "N")) {
00344         idist = 3;
00345     } else {
00346         idist = -1;
00347     }
00348 
00349 /*     Check EI */
00350 
00351     useei = TRUE_;
00352     badei = FALSE_;
00353     if (lsame_(ei + 1, " ") || *mode != 0) {
00354         useei = FALSE_;
00355     } else {
00356         if (lsame_(ei + 1, "R")) {
00357             i__1 = *n;
00358             for (j = 2; j <= i__1; ++j) {
00359                 if (lsame_(ei + j, "I")) {
00360                     if (lsame_(ei + (j - 1), "I")) {
00361                         badei = TRUE_;
00362                     }
00363                 } else {
00364                     if (! lsame_(ei + j, "R")) {
00365                         badei = TRUE_;
00366                     }
00367                 }
00368 /* L10: */
00369             }
00370         } else {
00371             badei = TRUE_;
00372         }
00373     }
00374 
00375 /*     Decode RSIGN */
00376 
00377     if (lsame_(rsign, "T")) {
00378         irsign = 1;
00379     } else if (lsame_(rsign, "F")) {
00380         irsign = 0;
00381     } else {
00382         irsign = -1;
00383     }
00384 
00385 /*     Decode UPPER */
00386 
00387     if (lsame_(upper, "T")) {
00388         iupper = 1;
00389     } else if (lsame_(upper, "F")) {
00390         iupper = 0;
00391     } else {
00392         iupper = -1;
00393     }
00394 
00395 /*     Decode SIM */
00396 
00397     if (lsame_(sim, "T")) {
00398         isim = 1;
00399     } else if (lsame_(sim, "F")) {
00400         isim = 0;
00401     } else {
00402         isim = -1;
00403     }
00404 
00405 /*     Check DS, if MODES=0 and ISIM=1 */
00406 
00407     bads = FALSE_;
00408     if (*modes == 0 && isim == 1) {
00409         i__1 = *n;
00410         for (j = 1; j <= i__1; ++j) {
00411             if (ds[j] == 0.f) {
00412                 bads = TRUE_;
00413             }
00414 /* L20: */
00415         }
00416     }
00417 
00418 /*     Set INFO if an error */
00419 
00420     if (*n < 0) {
00421         *info = -1;
00422     } else if (idist == -1) {
00423         *info = -2;
00424     } else if (abs(*mode) > 6) {
00425         *info = -5;
00426     } else if (*mode != 0 && abs(*mode) != 6 && *cond < 1.f) {
00427         *info = -6;
00428     } else if (badei) {
00429         *info = -8;
00430     } else if (irsign == -1) {
00431         *info = -9;
00432     } else if (iupper == -1) {
00433         *info = -10;
00434     } else if (isim == -1) {
00435         *info = -11;
00436     } else if (bads) {
00437         *info = -12;
00438     } else if (isim == 1 && abs(*modes) > 5) {
00439         *info = -13;
00440     } else if (isim == 1 && *modes != 0 && *conds < 1.f) {
00441         *info = -14;
00442     } else if (*kl < 1) {
00443         *info = -15;
00444     } else if (*ku < 1 || *ku < *n - 1 && *kl < *n - 1) {
00445         *info = -16;
00446     } else if (*lda < max(1,*n)) {
00447         *info = -19;
00448     }
00449 
00450     if (*info != 0) {
00451         i__1 = -(*info);
00452         xerbla_("SLATME", &i__1);
00453         return 0;
00454     }
00455 
00456 /*     Initialize random number generator */
00457 
00458     for (i__ = 1; i__ <= 4; ++i__) {
00459         iseed[i__] = (i__1 = iseed[i__], abs(i__1)) % 4096;
00460 /* L30: */
00461     }
00462 
00463     if (iseed[4] % 2 != 1) {
00464         ++iseed[4];
00465     }
00466 
00467 /*     2)      Set up diagonal of A */
00468 
00469 /*             Compute D according to COND and MODE */
00470 
00471     slatm1_(mode, cond, &irsign, &idist, &iseed[1], &d__[1], n, &iinfo);
00472     if (iinfo != 0) {
00473         *info = 1;
00474         return 0;
00475     }
00476     if (*mode != 0 && abs(*mode) != 6) {
00477 
00478 /*        Scale by DMAX */
00479 
00480         temp = dabs(d__[1]);
00481         i__1 = *n;
00482         for (i__ = 2; i__ <= i__1; ++i__) {
00483 /* Computing MAX */
00484             r__2 = temp, r__3 = (r__1 = d__[i__], dabs(r__1));
00485             temp = dmax(r__2,r__3);
00486 /* L40: */
00487         }
00488 
00489         if (temp > 0.f) {
00490             alpha = *dmax__ / temp;
00491         } else if (*dmax__ != 0.f) {
00492             *info = 2;
00493             return 0;
00494         } else {
00495             alpha = 0.f;
00496         }
00497 
00498         sscal_(n, &alpha, &d__[1], &c__1);
00499 
00500     }
00501 
00502     slaset_("Full", n, n, &c_b23, &c_b23, &a[a_offset], lda);
00503     i__1 = *lda + 1;
00504     scopy_(n, &d__[1], &c__1, &a[a_offset], &i__1);
00505 
00506 /*     Set up complex conjugate pairs */
00507 
00508     if (*mode == 0) {
00509         if (useei) {
00510             i__1 = *n;
00511             for (j = 2; j <= i__1; ++j) {
00512                 if (lsame_(ei + j, "I")) {
00513                     a[j - 1 + j * a_dim1] = a[j + j * a_dim1];
00514                     a[j + (j - 1) * a_dim1] = -a[j + j * a_dim1];
00515                     a[j + j * a_dim1] = a[j - 1 + (j - 1) * a_dim1];
00516                 }
00517 /* L50: */
00518             }
00519         }
00520 
00521     } else if (abs(*mode) == 5) {
00522 
00523         i__1 = *n;
00524         for (j = 2; j <= i__1; j += 2) {
00525             if (slaran_(&iseed[1]) > .5f) {
00526                 a[j - 1 + j * a_dim1] = a[j + j * a_dim1];
00527                 a[j + (j - 1) * a_dim1] = -a[j + j * a_dim1];
00528                 a[j + j * a_dim1] = a[j - 1 + (j - 1) * a_dim1];
00529             }
00530 /* L60: */
00531         }
00532     }
00533 
00534 /*     3)      If UPPER='T', set upper triangle of A to random numbers. */
00535 /*             (but don't modify the corners of 2x2 blocks.) */
00536 
00537     if (iupper != 0) {
00538         i__1 = *n;
00539         for (jc = 2; jc <= i__1; ++jc) {
00540             if (a[jc - 1 + jc * a_dim1] != 0.f) {
00541                 jr = jc - 2;
00542             } else {
00543                 jr = jc - 1;
00544             }
00545             slarnv_(&idist, &iseed[1], &jr, &a[jc * a_dim1 + 1]);
00546 /* L70: */
00547         }
00548     }
00549 
00550 /*     4)      If SIM='T', apply similarity transformation. */
00551 
00552 /*                                -1 */
00553 /*             Transform is  X A X  , where X = U S V, thus */
00554 
00555 /*             it is  U S V A V' (1/S) U' */
00556 
00557     if (isim != 0) {
00558 
00559 /*        Compute S (singular values of the eigenvector matrix) */
00560 /*        according to CONDS and MODES */
00561 
00562         slatm1_(modes, conds, &c__0, &c__0, &iseed[1], &ds[1], n, &iinfo);
00563         if (iinfo != 0) {
00564             *info = 3;
00565             return 0;
00566         }
00567 
00568 /*        Multiply by V and V' */
00569 
00570         slarge_(n, &a[a_offset], lda, &iseed[1], &work[1], &iinfo);
00571         if (iinfo != 0) {
00572             *info = 4;
00573             return 0;
00574         }
00575 
00576 /*        Multiply by S and (1/S) */
00577 
00578         i__1 = *n;
00579         for (j = 1; j <= i__1; ++j) {
00580             sscal_(n, &ds[j], &a[j + a_dim1], lda);
00581             if (ds[j] != 0.f) {
00582                 r__1 = 1.f / ds[j];
00583                 sscal_(n, &r__1, &a[j * a_dim1 + 1], &c__1);
00584             } else {
00585                 *info = 5;
00586                 return 0;
00587             }
00588 /* L80: */
00589         }
00590 
00591 /*        Multiply by U and U' */
00592 
00593         slarge_(n, &a[a_offset], lda, &iseed[1], &work[1], &iinfo);
00594         if (iinfo != 0) {
00595             *info = 4;
00596             return 0;
00597         }
00598     }
00599 
00600 /*     5)      Reduce the bandwidth. */
00601 
00602     if (*kl < *n - 1) {
00603 
00604 /*        Reduce bandwidth -- kill column */
00605 
00606         i__1 = *n - 1;
00607         for (jcr = *kl + 1; jcr <= i__1; ++jcr) {
00608             ic = jcr - *kl;
00609             irows = *n + 1 - jcr;
00610             icols = *n + *kl - jcr;
00611 
00612             scopy_(&irows, &a[jcr + ic * a_dim1], &c__1, &work[1], &c__1);
00613             xnorms = work[1];
00614             slarfg_(&irows, &xnorms, &work[2], &c__1, &tau);
00615             work[1] = 1.f;
00616 
00617             sgemv_("T", &irows, &icols, &c_b39, &a[jcr + (ic + 1) * a_dim1], 
00618                     lda, &work[1], &c__1, &c_b23, &work[irows + 1], &c__1);
00619             r__1 = -tau;
00620             sger_(&irows, &icols, &r__1, &work[1], &c__1, &work[irows + 1], &
00621                     c__1, &a[jcr + (ic + 1) * a_dim1], lda);
00622 
00623             sgemv_("N", n, &irows, &c_b39, &a[jcr * a_dim1 + 1], lda, &work[1]
00624 , &c__1, &c_b23, &work[irows + 1], &c__1);
00625             r__1 = -tau;
00626             sger_(n, &irows, &r__1, &work[irows + 1], &c__1, &work[1], &c__1, 
00627                     &a[jcr * a_dim1 + 1], lda);
00628 
00629             a[jcr + ic * a_dim1] = xnorms;
00630             i__2 = irows - 1;
00631             slaset_("Full", &i__2, &c__1, &c_b23, &c_b23, &a[jcr + 1 + ic * 
00632                     a_dim1], lda);
00633 /* L90: */
00634         }
00635     } else if (*ku < *n - 1) {
00636 
00637 /*        Reduce upper bandwidth -- kill a row at a time. */
00638 
00639         i__1 = *n - 1;
00640         for (jcr = *ku + 1; jcr <= i__1; ++jcr) {
00641             ir = jcr - *ku;
00642             irows = *n + *ku - jcr;
00643             icols = *n + 1 - jcr;
00644 
00645             scopy_(&icols, &a[ir + jcr * a_dim1], lda, &work[1], &c__1);
00646             xnorms = work[1];
00647             slarfg_(&icols, &xnorms, &work[2], &c__1, &tau);
00648             work[1] = 1.f;
00649 
00650             sgemv_("N", &irows, &icols, &c_b39, &a[ir + 1 + jcr * a_dim1], 
00651                     lda, &work[1], &c__1, &c_b23, &work[icols + 1], &c__1);
00652             r__1 = -tau;
00653             sger_(&irows, &icols, &r__1, &work[icols + 1], &c__1, &work[1], &
00654                     c__1, &a[ir + 1 + jcr * a_dim1], lda);
00655 
00656             sgemv_("C", &icols, n, &c_b39, &a[jcr + a_dim1], lda, &work[1], &
00657                     c__1, &c_b23, &work[icols + 1], &c__1);
00658             r__1 = -tau;
00659             sger_(&icols, n, &r__1, &work[1], &c__1, &work[icols + 1], &c__1, 
00660                     &a[jcr + a_dim1], lda);
00661 
00662             a[ir + jcr * a_dim1] = xnorms;
00663             i__2 = icols - 1;
00664             slaset_("Full", &c__1, &i__2, &c_b23, &c_b23, &a[ir + (jcr + 1) * 
00665                     a_dim1], lda);
00666 /* L100: */
00667         }
00668     }
00669 
00670 /*     Scale the matrix to have norm ANORM */
00671 
00672     if (*anorm >= 0.f) {
00673         temp = slange_("M", n, n, &a[a_offset], lda, tempa);
00674         if (temp > 0.f) {
00675             alpha = *anorm / temp;
00676             i__1 = *n;
00677             for (j = 1; j <= i__1; ++j) {
00678                 sscal_(n, &alpha, &a[j * a_dim1 + 1], &c__1);
00679 /* L110: */
00680             }
00681         }
00682     }
00683 
00684     return 0;
00685 
00686 /*     End of SLATME */
00687 
00688 } /* slatme_ */


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