LAPACK 3.3.1
Linear Algebra PACKage

dlatme.f

Go to the documentation of this file.
00001       SUBROUTINE DLATME( N, DIST, ISEED, D, MODE, COND, DMAX, EI, 
00002      $  RSIGN, 
00003      $                   UPPER, SIM, DS, MODES, CONDS, KL, KU, ANORM, 
00004      $  A, 
00005      $                   LDA, WORK, INFO )
00006 *
00007 *  -- LAPACK test routine (version 3.1) --
00008 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00009 *     June 2010
00010 *
00011 *     .. Scalar Arguments ..
00012       CHARACTER          DIST, RSIGN, SIM, UPPER
00013       INTEGER            INFO, KL, KU, LDA, MODE, MODES, N
00014       DOUBLE PRECISION   ANORM, COND, CONDS, DMAX
00015 *     ..
00016 *     .. Array Arguments ..
00017       CHARACTER          EI( * )
00018       INTEGER            ISEED( 4 )
00019       DOUBLE PRECISION   A( LDA, * ), D( * ), DS( * ), WORK( * )
00020 *     ..
00021 *
00022 *  Purpose
00023 *  =======
00024 *
00025 *     DLATME generates random non-symmetric square matrices with
00026 *     specified eigenvalues for testing LAPACK programs.
00027 *
00028 *     DLATME operates by applying the following sequence of
00029 *     operations:
00030 *
00031 *     1. Set the diagonal to D, where D may be input or
00032 *          computed according to MODE, COND, DMAX, and RSIGN
00033 *          as described below.
00034 *
00035 *     2. If complex conjugate pairs are desired (MODE=0 and EI(1)='R',
00036 *          or MODE=5), certain pairs of adjacent elements of D are
00037 *          interpreted as the real and complex parts of a complex
00038 *          conjugate pair; A thus becomes block diagonal, with 1x1
00039 *          and 2x2 blocks.
00040 *
00041 *     3. If UPPER='T', the upper triangle of A is set to random values
00042 *          out of distribution DIST.
00043 *
00044 *     4. If SIM='T', A is multiplied on the left by a random matrix
00045 *          X, whose singular values are specified by DS, MODES, and
00046 *          CONDS, and on the right by X inverse.
00047 *
00048 *     5. If KL < N-1, the lower bandwidth is reduced to KL using
00049 *          Householder transformations.  If KU < N-1, the upper
00050 *          bandwidth is reduced to KU.
00051 *
00052 *     6. If ANORM is not negative, the matrix is scaled to have
00053 *          maximum-element-norm ANORM.
00054 *
00055 *     (Note: since the matrix cannot be reduced beyond Hessenberg form,
00056 *      no packing options are available.)
00057 *
00058 *  Arguments
00059 *  =========
00060 *
00061 *  N        (input) INTEGER
00062 *           The number of columns (or rows) of A. Not modified.
00063 *
00064 *  DIST     (input) CHARACTER*1
00065 *           On entry, DIST specifies the type of distribution to be used
00066 *           to generate the random eigen-/singular values, and for the
00067 *           upper triangle (see UPPER).
00068 *           'U' => UNIFORM( 0, 1 )  ( 'U' for uniform )
00069 *           'S' => UNIFORM( -1, 1 ) ( 'S' for symmetric )
00070 *           'N' => NORMAL( 0, 1 )   ( 'N' for normal )
00071 *           Not modified.
00072 *
00073 *  ISEED    (input/output) INTEGER array, dimension ( 4 )
00074 *           On entry ISEED specifies the seed of the random number
00075 *           generator. They should lie between 0 and 4095 inclusive,
00076 *           and ISEED(4) should be odd. The random number generator
00077 *           uses a linear congruential sequence limited to small
00078 *           integers, and so should produce machine independent
00079 *           random numbers. The values of ISEED are changed on
00080 *           exit, and can be used in the next call to DLATME
00081 *           to continue the same random number sequence.
00082 *           Changed on exit.
00083 *
00084 *  D        (input/output) DOUBLE PRECISION array, dimension ( N )
00085 *           This array is used to specify the eigenvalues of A.  If
00086 *           MODE=0, then D is assumed to contain the eigenvalues (but
00087 *           see the description of EI), otherwise they will be
00088 *           computed according to MODE, COND, DMAX, and RSIGN and
00089 *           placed in D.
00090 *           Modified if MODE is nonzero.
00091 *
00092 *  MODE     (input) INTEGER
00093 *           On entry this describes how the eigenvalues are to
00094 *           be specified:
00095 *           MODE = 0 means use D (with EI) as input
00096 *           MODE = 1 sets D(1)=1 and D(2:N)=1.0/COND
00097 *           MODE = 2 sets D(1:N-1)=1 and D(N)=1.0/COND
00098 *           MODE = 3 sets D(I)=COND**(-(I-1)/(N-1))
00099 *           MODE = 4 sets D(i)=1 - (i-1)/(N-1)*(1 - 1/COND)
00100 *           MODE = 5 sets D to random numbers in the range
00101 *                    ( 1/COND , 1 ) such that their logarithms
00102 *                    are uniformly distributed.  Each odd-even pair
00103 *                    of elements will be either used as two real
00104 *                    eigenvalues or as the real and imaginary part
00105 *                    of a complex conjugate pair of eigenvalues;
00106 *                    the choice of which is done is random, with
00107 *                    50-50 probability, for each pair.
00108 *           MODE = 6 set D to random numbers from same distribution
00109 *                    as the rest of the matrix.
00110 *           MODE < 0 has the same meaning as ABS(MODE), except that
00111 *              the order of the elements of D is reversed.
00112 *           Thus if MODE is between 1 and 4, D has entries ranging
00113 *              from 1 to 1/COND, if between -1 and -4, D has entries
00114 *              ranging from 1/COND to 1,
00115 *           Not modified.
00116 *
00117 *  COND     (input) DOUBLE PRECISION
00118 *           On entry, this is used as described under MODE above.
00119 *           If used, it must be >= 1. Not modified.
00120 *
00121 *  DMAX     (input) DOUBLE PRECISION
00122 *           If MODE is neither -6, 0 nor 6, the contents of D, as
00123 *           computed according to MODE and COND, will be scaled by
00124 *           DMAX / max(abs(D(i))).  Note that DMAX need not be
00125 *           positive: if DMAX is negative (or zero), D will be
00126 *           scaled by a negative number (or zero).
00127 *           Not modified.
00128 *
00129 *  EI       (input) CHARACTER*1 array, dimension ( N )
00130 *           If MODE is 0, and EI(1) is not ' ' (space character),
00131 *           this array specifies which elements of D (on input) are
00132 *           real eigenvalues and which are the real and imaginary parts
00133 *           of a complex conjugate pair of eigenvalues.  The elements
00134 *           of EI may then only have the values 'R' and 'I'.  If
00135 *           EI(j)='R' and EI(j+1)='I', then the j-th eigenvalue is
00136 *           CMPLX( D(j) , D(j+1) ), and the (j+1)-th is the complex
00137 *           conjugate thereof.  If EI(j)=EI(j+1)='R', then the j-th
00138 *           eigenvalue is D(j) (i.e., real).  EI(1) may not be 'I',
00139 *           nor may two adjacent elements of EI both have the value 'I'.
00140 *           If MODE is not 0, then EI is ignored.  If MODE is 0 and
00141 *           EI(1)=' ', then the eigenvalues will all be real.
00142 *           Not modified.
00143 *
00144 *  RSIGN    (input) CHARACTER*1
00145 *           If MODE is not 0, 6, or -6, and RSIGN='T', then the
00146 *           elements of D, as computed according to MODE and COND, will
00147 *           be multiplied by a random sign (+1 or -1).  If RSIGN='F',
00148 *           they will not be.  RSIGN may only have the values 'T' or
00149 *           'F'.
00150 *           Not modified.
00151 *
00152 *  UPPER    (input) CHARACTER*1
00153 *           If UPPER='T', then the elements of A above the diagonal
00154 *           (and above the 2x2 diagonal blocks, if A has complex
00155 *           eigenvalues) will be set to random numbers out of DIST.
00156 *           If UPPER='F', they will not.  UPPER may only have the
00157 *           values 'T' or 'F'.
00158 *           Not modified.
00159 *
00160 *  SIM      (input) CHARACTER*1
00161 *           If SIM='T', then A will be operated on by a "similarity
00162 *           transform", i.e., multiplied on the left by a matrix X and
00163 *           on the right by X inverse.  X = U S V, where U and V are
00164 *           random unitary matrices and S is a (diagonal) matrix of
00165 *           singular values specified by DS, MODES, and CONDS.  If
00166 *           SIM='F', then A will not be transformed.
00167 *           Not modified.
00168 *
00169 *  DS       (input/output) DOUBLE PRECISION array, dimension ( N )
00170 *           This array is used to specify the singular values of X,
00171 *           in the same way that D specifies the eigenvalues of A.
00172 *           If MODE=0, the DS contains the singular values, which
00173 *           may not be zero.
00174 *           Modified if MODE is nonzero.
00175 *
00176 *  MODES    (input) INTEGER
00177 *
00178 *  CONDS    (input) DOUBLE PRECISION
00179 *           Same as MODE and COND, but for specifying the diagonal
00180 *           of S.  MODES=-6 and +6 are not allowed (since they would
00181 *           result in randomly ill-conditioned eigenvalues.)
00182 *
00183 *  KL       (input) INTEGER
00184 *           This specifies the lower bandwidth of the  matrix.  KL=1
00185 *           specifies upper Hessenberg form.  If KL is at least N-1,
00186 *           then A will have full lower bandwidth.  KL must be at
00187 *           least 1.
00188 *           Not modified.
00189 *
00190 *  KU       (input) INTEGER
00191 *           This specifies the upper bandwidth of the  matrix.  KU=1
00192 *           specifies lower Hessenberg form.  If KU is at least N-1,
00193 *           then A will have full upper bandwidth; if KU and KL
00194 *           are both at least N-1, then A will be dense.  Only one of
00195 *           KU and KL may be less than N-1.  KU must be at least 1.
00196 *           Not modified.
00197 *
00198 *  ANORM    (input) DOUBLE PRECISION
00199 *           If ANORM is not negative, then A will be scaled by a non-
00200 *           negative real number to make the maximum-element-norm of A
00201 *           to be ANORM.
00202 *           Not modified.
00203 *
00204 *  A        (output) DOUBLE PRECISION array, dimension ( LDA, N )
00205 *           On exit A is the desired test matrix.
00206 *           Modified.
00207 *
00208 *  LDA      (input) INTEGER
00209 *           LDA specifies the first dimension of A as declared in the
00210 *           calling program.  LDA must be at least N.
00211 *           Not modified.
00212 *
00213 *  WORK     (workspace) DOUBLE PRECISION array, dimension ( 3*N )
00214 *           Workspace.
00215 *           Modified.
00216 *
00217 *  INFO     (output) INTEGER
00218 *           Error code.  On exit, INFO will be set to one of the
00219 *           following values:
00220 *             0 => normal return
00221 *            -1 => N negative
00222 *            -2 => DIST illegal string
00223 *            -5 => MODE not in range -6 to 6
00224 *            -6 => COND less than 1.0, and MODE neither -6, 0 nor 6
00225 *            -8 => EI(1) is not ' ' or 'R', EI(j) is not 'R' or 'I', or
00226 *                  two adjacent elements of EI are 'I'.
00227 *            -9 => RSIGN is not 'T' or 'F'
00228 *           -10 => UPPER is not 'T' or 'F'
00229 *           -11 => SIM   is not 'T' or 'F'
00230 *           -12 => MODES=0 and DS has a zero singular value.
00231 *           -13 => MODES is not in the range -5 to 5.
00232 *           -14 => MODES is nonzero and CONDS is less than 1.
00233 *           -15 => KL is less than 1.
00234 *           -16 => KU is less than 1, or KL and KU are both less than
00235 *                  N-1.
00236 *           -19 => LDA is less than N.
00237 *            1  => Error return from DLATM1 (computing D)
00238 *            2  => Cannot scale to DMAX (max. eigenvalue is 0)
00239 *            3  => Error return from DLATM1 (computing DS)
00240 *            4  => Error return from DLARGE
00241 *            5  => Zero singular value from DLATM1.
00242 *
00243 *  =====================================================================
00244 *
00245 *     .. Parameters ..
00246       DOUBLE PRECISION   ZERO
00247       PARAMETER          ( ZERO = 0.0D0 )
00248       DOUBLE PRECISION   ONE
00249       PARAMETER          ( ONE = 1.0D0 )
00250       DOUBLE PRECISION   HALF
00251       PARAMETER          ( HALF = 1.0D0 / 2.0D0 )
00252 *     ..
00253 *     .. Local Scalars ..
00254       LOGICAL            BADEI, BADS, USEEI
00255       INTEGER            I, IC, ICOLS, IDIST, IINFO, IR, IROWS, IRSIGN,
00256      $                   ISIM, IUPPER, J, JC, JCR, JR
00257       DOUBLE PRECISION   ALPHA, TAU, TEMP, XNORMS
00258 *     ..
00259 *     .. Local Arrays ..
00260       DOUBLE PRECISION   TEMPA( 1 )
00261 *     ..
00262 *     .. External Functions ..
00263       LOGICAL            LSAME
00264       DOUBLE PRECISION   DLANGE, DLARAN
00265       EXTERNAL           LSAME, DLANGE, DLARAN
00266 *     ..
00267 *     .. External Subroutines ..
00268       EXTERNAL           DCOPY, DGEMV, DGER, DLARFG, DLARGE, DLARNV,
00269      $                   DLASET, DLATM1, DSCAL, XERBLA
00270 *     ..
00271 *     .. Intrinsic Functions ..
00272       INTRINSIC          ABS, MAX, MOD
00273 *     ..
00274 *     .. Executable Statements ..
00275 *
00276 *     1)      Decode and Test the input parameters.
00277 *             Initialize flags & seed.
00278 *
00279       INFO = 0
00280 *
00281 *     Quick return if possible
00282 *
00283       IF( N.EQ.0 )
00284      $   RETURN
00285 *
00286 *     Decode DIST
00287 *
00288       IF( LSAME( DIST, 'U' ) ) THEN
00289          IDIST = 1
00290       ELSE IF( LSAME( DIST, 'S' ) ) THEN
00291          IDIST = 2
00292       ELSE IF( LSAME( DIST, 'N' ) ) THEN
00293          IDIST = 3
00294       ELSE
00295          IDIST = -1
00296       END IF
00297 *
00298 *     Check EI
00299 *
00300       USEEI = .TRUE.
00301       BADEI = .FALSE.
00302       IF( LSAME( EI( 1 ), ' ' ) .OR. MODE.NE.0 ) THEN
00303          USEEI = .FALSE.
00304       ELSE
00305          IF( LSAME( EI( 1 ), 'R' ) ) THEN
00306             DO 10 J = 2, N
00307                IF( LSAME( EI( J ), 'I' ) ) THEN
00308                   IF( LSAME( EI( J-1 ), 'I' ) )
00309      $               BADEI = .TRUE.
00310                ELSE
00311                   IF( .NOT.LSAME( EI( J ), 'R' ) )
00312      $               BADEI = .TRUE.
00313                END IF
00314    10       CONTINUE
00315          ELSE
00316             BADEI = .TRUE.
00317          END IF
00318       END IF
00319 *
00320 *     Decode RSIGN
00321 *
00322       IF( LSAME( RSIGN, 'T' ) ) THEN
00323          IRSIGN = 1
00324       ELSE IF( LSAME( RSIGN, 'F' ) ) THEN
00325          IRSIGN = 0
00326       ELSE
00327          IRSIGN = -1
00328       END IF
00329 *
00330 *     Decode UPPER
00331 *
00332       IF( LSAME( UPPER, 'T' ) ) THEN
00333          IUPPER = 1
00334       ELSE IF( LSAME( UPPER, 'F' ) ) THEN
00335          IUPPER = 0
00336       ELSE
00337          IUPPER = -1
00338       END IF
00339 *
00340 *     Decode SIM
00341 *
00342       IF( LSAME( SIM, 'T' ) ) THEN
00343          ISIM = 1
00344       ELSE IF( LSAME( SIM, 'F' ) ) THEN
00345          ISIM = 0
00346       ELSE
00347          ISIM = -1
00348       END IF
00349 *
00350 *     Check DS, if MODES=0 and ISIM=1
00351 *
00352       BADS = .FALSE.
00353       IF( MODES.EQ.0 .AND. ISIM.EQ.1 ) THEN
00354          DO 20 J = 1, N
00355             IF( DS( J ).EQ.ZERO )
00356      $         BADS = .TRUE.
00357    20    CONTINUE
00358       END IF
00359 *
00360 *     Set INFO if an error
00361 *
00362       IF( N.LT.0 ) THEN
00363          INFO = -1
00364       ELSE IF( IDIST.EQ.-1 ) THEN
00365          INFO = -2
00366       ELSE IF( ABS( MODE ).GT.6 ) THEN
00367          INFO = -5
00368       ELSE IF( ( MODE.NE.0 .AND. ABS( MODE ).NE.6 ) .AND. COND.LT.ONE )
00369      $          THEN
00370          INFO = -6
00371       ELSE IF( BADEI ) THEN
00372          INFO = -8
00373       ELSE IF( IRSIGN.EQ.-1 ) THEN
00374          INFO = -9
00375       ELSE IF( IUPPER.EQ.-1 ) THEN
00376          INFO = -10
00377       ELSE IF( ISIM.EQ.-1 ) THEN
00378          INFO = -11
00379       ELSE IF( BADS ) THEN
00380          INFO = -12
00381       ELSE IF( ISIM.EQ.1 .AND. ABS( MODES ).GT.5 ) THEN
00382          INFO = -13
00383       ELSE IF( ISIM.EQ.1 .AND. MODES.NE.0 .AND. CONDS.LT.ONE ) THEN
00384          INFO = -14
00385       ELSE IF( KL.LT.1 ) THEN
00386          INFO = -15
00387       ELSE IF( KU.LT.1 .OR. ( KU.LT.N-1 .AND. KL.LT.N-1 ) ) THEN
00388          INFO = -16
00389       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00390          INFO = -19
00391       END IF
00392 *
00393       IF( INFO.NE.0 ) THEN
00394          CALL XERBLA( 'DLATME', -INFO )
00395          RETURN
00396       END IF
00397 *
00398 *     Initialize random number generator
00399 *
00400       DO 30 I = 1, 4
00401          ISEED( I ) = MOD( ABS( ISEED( I ) ), 4096 )
00402    30 CONTINUE
00403 *
00404       IF( MOD( ISEED( 4 ), 2 ).NE.1 )
00405      $   ISEED( 4 ) = ISEED( 4 ) + 1
00406 *
00407 *     2)      Set up diagonal of A
00408 *
00409 *             Compute D according to COND and MODE
00410 *
00411       CALL DLATM1( MODE, COND, IRSIGN, IDIST, ISEED, D, N, IINFO )
00412       IF( IINFO.NE.0 ) THEN
00413          INFO = 1
00414          RETURN
00415       END IF
00416       IF( MODE.NE.0 .AND. ABS( MODE ).NE.6 ) THEN
00417 *
00418 *        Scale by DMAX
00419 *
00420          TEMP = ABS( D( 1 ) )
00421          DO 40 I = 2, N
00422             TEMP = MAX( TEMP, ABS( D( I ) ) )
00423    40    CONTINUE
00424 *
00425          IF( TEMP.GT.ZERO ) THEN
00426             ALPHA = DMAX / TEMP
00427          ELSE IF( DMAX.NE.ZERO ) THEN
00428             INFO = 2
00429             RETURN
00430          ELSE
00431             ALPHA = ZERO
00432          END IF
00433 *
00434          CALL DSCAL( N, ALPHA, D, 1 )
00435 *
00436       END IF
00437 *
00438       CALL DLASET( 'Full', N, N, ZERO, ZERO, A, LDA )
00439       CALL DCOPY( N, D, 1, A, LDA+1 )
00440 *
00441 *     Set up complex conjugate pairs
00442 *
00443       IF( MODE.EQ.0 ) THEN
00444          IF( USEEI ) THEN
00445             DO 50 J = 2, N
00446                IF( LSAME( EI( J ), 'I' ) ) THEN
00447                   A( J-1, J ) = A( J, J )
00448                   A( J, J-1 ) = -A( J, J )
00449                   A( J, J ) = A( J-1, J-1 )
00450                END IF
00451    50       CONTINUE
00452          END IF
00453 *
00454       ELSE IF( ABS( MODE ).EQ.5 ) THEN
00455 *
00456          DO 60 J = 2, N, 2
00457             IF( DLARAN( ISEED ).GT.HALF ) THEN
00458                A( J-1, J ) = A( J, J )
00459                A( J, J-1 ) = -A( J, J )
00460                A( J, J ) = A( J-1, J-1 )
00461             END IF
00462    60    CONTINUE
00463       END IF
00464 *
00465 *     3)      If UPPER='T', set upper triangle of A to random numbers.
00466 *             (but don't modify the corners of 2x2 blocks.)
00467 *
00468       IF( IUPPER.NE.0 ) THEN
00469          DO 70 JC = 2, N
00470             IF( A( JC-1, JC ).NE.ZERO ) THEN
00471                JR = JC - 2
00472             ELSE
00473                JR = JC - 1
00474             END IF
00475             CALL DLARNV( IDIST, ISEED, JR, A( 1, JC ) )
00476    70    CONTINUE
00477       END IF
00478 *
00479 *     4)      If SIM='T', apply similarity transformation.
00480 *
00481 *                                -1
00482 *             Transform is  X A X  , where X = U S V, thus
00483 *
00484 *             it is  U S V A V' (1/S) U'
00485 *
00486       IF( ISIM.NE.0 ) THEN
00487 *
00488 *        Compute S (singular values of the eigenvector matrix)
00489 *        according to CONDS and MODES
00490 *
00491          CALL DLATM1( MODES, CONDS, 0, 0, ISEED, DS, N, IINFO )
00492          IF( IINFO.NE.0 ) THEN
00493             INFO = 3
00494             RETURN
00495          END IF
00496 *
00497 *        Multiply by V and V'
00498 *
00499          CALL DLARGE( N, A, LDA, ISEED, WORK, IINFO )
00500          IF( IINFO.NE.0 ) THEN
00501             INFO = 4
00502             RETURN
00503          END IF
00504 *
00505 *        Multiply by S and (1/S)
00506 *
00507          DO 80 J = 1, N
00508             CALL DSCAL( N, DS( J ), A( J, 1 ), LDA )
00509             IF( DS( J ).NE.ZERO ) THEN
00510                CALL DSCAL( N, ONE / DS( J ), A( 1, J ), 1 )
00511             ELSE
00512                INFO = 5
00513                RETURN
00514             END IF
00515    80    CONTINUE
00516 *
00517 *        Multiply by U and U'
00518 *
00519          CALL DLARGE( N, A, LDA, ISEED, WORK, IINFO )
00520          IF( IINFO.NE.0 ) THEN
00521             INFO = 4
00522             RETURN
00523          END IF
00524       END IF
00525 *
00526 *     5)      Reduce the bandwidth.
00527 *
00528       IF( KL.LT.N-1 ) THEN
00529 *
00530 *        Reduce bandwidth -- kill column
00531 *
00532          DO 90 JCR = KL + 1, N - 1
00533             IC = JCR - KL
00534             IROWS = N + 1 - JCR
00535             ICOLS = N + KL - JCR
00536 *
00537             CALL DCOPY( IROWS, A( JCR, IC ), 1, WORK, 1 )
00538             XNORMS = WORK( 1 )
00539             CALL DLARFG( IROWS, XNORMS, WORK( 2 ), 1, TAU )
00540             WORK( 1 ) = ONE
00541 *
00542             CALL DGEMV( 'T', IROWS, ICOLS, ONE, A( JCR, IC+1 ), LDA,
00543      $                  WORK, 1, ZERO, WORK( IROWS+1 ), 1 )
00544             CALL DGER( IROWS, ICOLS, -TAU, WORK, 1, WORK( IROWS+1 ), 1,
00545      $                 A( JCR, IC+1 ), LDA )
00546 *
00547             CALL DGEMV( 'N', N, IROWS, ONE, A( 1, JCR ), LDA, WORK, 1,
00548      $                  ZERO, WORK( IROWS+1 ), 1 )
00549             CALL DGER( N, IROWS, -TAU, WORK( IROWS+1 ), 1, WORK, 1,
00550      $                 A( 1, JCR ), LDA )
00551 *
00552             A( JCR, IC ) = XNORMS
00553             CALL DLASET( 'Full', IROWS-1, 1, ZERO, ZERO, A( JCR+1, IC ),
00554      $                   LDA )
00555    90    CONTINUE
00556       ELSE IF( KU.LT.N-1 ) THEN
00557 *
00558 *        Reduce upper bandwidth -- kill a row at a time.
00559 *
00560          DO 100 JCR = KU + 1, N - 1
00561             IR = JCR - KU
00562             IROWS = N + KU - JCR
00563             ICOLS = N + 1 - JCR
00564 *
00565             CALL DCOPY( ICOLS, A( IR, JCR ), LDA, WORK, 1 )
00566             XNORMS = WORK( 1 )
00567             CALL DLARFG( ICOLS, XNORMS, WORK( 2 ), 1, TAU )
00568             WORK( 1 ) = ONE
00569 *
00570             CALL DGEMV( 'N', IROWS, ICOLS, ONE, A( IR+1, JCR ), LDA,
00571      $                  WORK, 1, ZERO, WORK( ICOLS+1 ), 1 )
00572             CALL DGER( IROWS, ICOLS, -TAU, WORK( ICOLS+1 ), 1, WORK, 1,
00573      $                 A( IR+1, JCR ), LDA )
00574 *
00575             CALL DGEMV( 'C', ICOLS, N, ONE, A( JCR, 1 ), LDA, WORK, 1,
00576      $                  ZERO, WORK( ICOLS+1 ), 1 )
00577             CALL DGER( ICOLS, N, -TAU, WORK, 1, WORK( ICOLS+1 ), 1,
00578      $                 A( JCR, 1 ), LDA )
00579 *
00580             A( IR, JCR ) = XNORMS
00581             CALL DLASET( 'Full', 1, ICOLS-1, ZERO, ZERO, A( IR, JCR+1 ),
00582      $                   LDA )
00583   100    CONTINUE
00584       END IF
00585 *
00586 *     Scale the matrix to have norm ANORM
00587 *
00588       IF( ANORM.GE.ZERO ) THEN
00589          TEMP = DLANGE( 'M', N, N, A, LDA, TEMPA )
00590          IF( TEMP.GT.ZERO ) THEN
00591             ALPHA = ANORM / TEMP
00592             DO 110 J = 1, N
00593                CALL DSCAL( N, ALPHA, A( 1, J ), 1 )
00594   110       CONTINUE
00595          END IF
00596       END IF
00597 *
00598       RETURN
00599 *
00600 *     End of DLATME
00601 *
00602       END
 All Files Functions