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


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