LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dlatme.f
Go to the documentation of this file.
1 *> \brief \b DLATME
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 * Definition:
9 * ===========
10 *
11 * SUBROUTINE DLATME( N, DIST, ISEED, D, MODE, COND, DMAX, EI,
12 * RSIGN,
13 * UPPER, SIM, DS, MODES, CONDS, KL, KU, ANORM,
14 * A,
15 * LDA, WORK, INFO )
16 *
17 * .. Scalar Arguments ..
18 * CHARACTER DIST, RSIGN, SIM, UPPER
19 * INTEGER INFO, KL, KU, LDA, MODE, MODES, N
20 * DOUBLE PRECISION ANORM, COND, CONDS, DMAX
21 * ..
22 * .. Array Arguments ..
23 * CHARACTER EI( * )
24 * INTEGER ISEED( 4 )
25 * DOUBLE PRECISION A( LDA, * ), D( * ), DS( * ), WORK( * )
26 * ..
27 *
28 *
29 *> \par Purpose:
30 * =============
31 *>
32 *> \verbatim
33 *>
34 *> DLATME generates random non-symmetric square matrices with
35 *> specified eigenvalues for testing LAPACK programs.
36 *>
37 *> DLATME operates by applying the following sequence of
38 *> operations:
39 *>
40 *> 1. Set the diagonal to D, where D may be input or
41 *> computed according to MODE, COND, DMAX, and RSIGN
42 *> as described below.
43 *>
44 *> 2. If complex conjugate pairs are desired (MODE=0 and EI(1)='R',
45 *> or MODE=5), certain pairs of adjacent elements of D are
46 *> interpreted as the real and complex parts of a complex
47 *> conjugate pair; A thus becomes block diagonal, with 1x1
48 *> and 2x2 blocks.
49 *>
50 *> 3. If UPPER='T', the upper triangle of A is set to random values
51 *> out of distribution DIST.
52 *>
53 *> 4. If SIM='T', A is multiplied on the left by a random matrix
54 *> X, whose singular values are specified by DS, MODES, and
55 *> CONDS, and on the right by X inverse.
56 *>
57 *> 5. If KL < N-1, the lower bandwidth is reduced to KL using
58 *> Householder transformations. If KU < N-1, the upper
59 *> bandwidth is reduced to KU.
60 *>
61 *> 6. If ANORM is not negative, the matrix is scaled to have
62 *> maximum-element-norm ANORM.
63 *>
64 *> (Note: since the matrix cannot be reduced beyond Hessenberg form,
65 *> no packing options are available.)
66 *> \endverbatim
67 *
68 * Arguments:
69 * ==========
70 *
71 *> \param[in] N
72 *> \verbatim
73 *> N is INTEGER
74 *> The number of columns (or rows) of A. Not modified.
75 *> \endverbatim
76 *>
77 *> \param[in] DIST
78 *> \verbatim
79 *> DIST is CHARACTER*1
80 *> On entry, DIST specifies the type of distribution to be used
81 *> to generate the random eigen-/singular values, and for the
82 *> upper triangle (see UPPER).
83 *> 'U' => UNIFORM( 0, 1 ) ( 'U' for uniform )
84 *> 'S' => UNIFORM( -1, 1 ) ( 'S' for symmetric )
85 *> 'N' => NORMAL( 0, 1 ) ( 'N' for normal )
86 *> Not modified.
87 *> \endverbatim
88 *>
89 *> \param[in,out] ISEED
90 *> \verbatim
91 *> ISEED is INTEGER array, dimension ( 4 )
92 *> On entry ISEED specifies the seed of the random number
93 *> generator. They should lie between 0 and 4095 inclusive,
94 *> and ISEED(4) should be odd. The random number generator
95 *> uses a linear congruential sequence limited to small
96 *> integers, and so should produce machine independent
97 *> random numbers. The values of ISEED are changed on
98 *> exit, and can be used in the next call to DLATME
99 *> to continue the same random number sequence.
100 *> Changed on exit.
101 *> \endverbatim
102 *>
103 *> \param[in,out] D
104 *> \verbatim
105 *> D is DOUBLE PRECISION array, dimension ( N )
106 *> This array is used to specify the eigenvalues of A. If
107 *> MODE=0, then D is assumed to contain the eigenvalues (but
108 *> see the description of EI), otherwise they will be
109 *> computed according to MODE, COND, DMAX, and RSIGN and
110 *> placed in D.
111 *> Modified if MODE is nonzero.
112 *> \endverbatim
113 *>
114 *> \param[in] MODE
115 *> \verbatim
116 *> MODE is INTEGER
117 *> On entry this describes how the eigenvalues are to
118 *> be specified:
119 *> MODE = 0 means use D (with EI) as input
120 *> MODE = 1 sets D(1)=1 and D(2:N)=1.0/COND
121 *> MODE = 2 sets D(1:N-1)=1 and D(N)=1.0/COND
122 *> MODE = 3 sets D(I)=COND**(-(I-1)/(N-1))
123 *> MODE = 4 sets D(i)=1 - (i-1)/(N-1)*(1 - 1/COND)
124 *> MODE = 5 sets D to random numbers in the range
125 *> ( 1/COND , 1 ) such that their logarithms
126 *> are uniformly distributed. Each odd-even pair
127 *> of elements will be either used as two real
128 *> eigenvalues or as the real and imaginary part
129 *> of a complex conjugate pair of eigenvalues;
130 *> the choice of which is done is random, with
131 *> 50-50 probability, for each pair.
132 *> MODE = 6 set D to random numbers from same distribution
133 *> as the rest of the matrix.
134 *> MODE < 0 has the same meaning as ABS(MODE), except that
135 *> the order of the elements of D is reversed.
136 *> Thus if MODE is between 1 and 4, D has entries ranging
137 *> from 1 to 1/COND, if between -1 and -4, D has entries
138 *> ranging from 1/COND to 1,
139 *> Not modified.
140 *> \endverbatim
141 *>
142 *> \param[in] COND
143 *> \verbatim
144 *> COND is DOUBLE PRECISION
145 *> On entry, this is used as described under MODE above.
146 *> If used, it must be >= 1. Not modified.
147 *> \endverbatim
148 *>
149 *> \param[in] DMAX
150 *> \verbatim
151 *> DMAX is DOUBLE PRECISION
152 *> If MODE is neither -6, 0 nor 6, the contents of D, as
153 *> computed according to MODE and COND, will be scaled by
154 *> DMAX / max(abs(D(i))). Note that DMAX need not be
155 *> positive: if DMAX is negative (or zero), D will be
156 *> scaled by a negative number (or zero).
157 *> Not modified.
158 *> \endverbatim
159 *>
160 *> \param[in] EI
161 *> \verbatim
162 *> EI is CHARACTER*1 array, dimension ( N )
163 *> If MODE is 0, and EI(1) is not ' ' (space character),
164 *> this array specifies which elements of D (on input) are
165 *> real eigenvalues and which are the real and imaginary parts
166 *> of a complex conjugate pair of eigenvalues. The elements
167 *> of EI may then only have the values 'R' and 'I'. If
168 *> EI(j)='R' and EI(j+1)='I', then the j-th eigenvalue is
169 *> CMPLX( D(j) , D(j+1) ), and the (j+1)-th is the complex
170 *> conjugate thereof. If EI(j)=EI(j+1)='R', then the j-th
171 *> eigenvalue is D(j) (i.e., real). EI(1) may not be 'I',
172 *> nor may two adjacent elements of EI both have the value 'I'.
173 *> If MODE is not 0, then EI is ignored. If MODE is 0 and
174 *> EI(1)=' ', then the eigenvalues will all be real.
175 *> Not modified.
176 *> \endverbatim
177 *>
178 *> \param[in] RSIGN
179 *> \verbatim
180 *> RSIGN is CHARACTER*1
181 *> If MODE is not 0, 6, or -6, and RSIGN='T', then the
182 *> elements of D, as computed according to MODE and COND, will
183 *> be multiplied by a random sign (+1 or -1). If RSIGN='F',
184 *> they will not be. RSIGN may only have the values 'T' or
185 *> 'F'.
186 *> Not modified.
187 *> \endverbatim
188 *>
189 *> \param[in] UPPER
190 *> \verbatim
191 *> UPPER is CHARACTER*1
192 *> If UPPER='T', then the elements of A above the diagonal
193 *> (and above the 2x2 diagonal blocks, if A has complex
194 *> eigenvalues) will be set to random numbers out of DIST.
195 *> If UPPER='F', they will not. UPPER may only have the
196 *> values 'T' or 'F'.
197 *> Not modified.
198 *> \endverbatim
199 *>
200 *> \param[in] SIM
201 *> \verbatim
202 *> SIM is CHARACTER*1
203 *> If SIM='T', then A will be operated on by a "similarity
204 *> transform", i.e., multiplied on the left by a matrix X and
205 *> on the right by X inverse. X = U S V, where U and V are
206 *> random unitary matrices and S is a (diagonal) matrix of
207 *> singular values specified by DS, MODES, and CONDS. If
208 *> SIM='F', then A will not be transformed.
209 *> Not modified.
210 *> \endverbatim
211 *>
212 *> \param[in,out] DS
213 *> \verbatim
214 *> DS is DOUBLE PRECISION array, dimension ( N )
215 *> This array is used to specify the singular values of X,
216 *> in the same way that D specifies the eigenvalues of A.
217 *> If MODE=0, the DS contains the singular values, which
218 *> may not be zero.
219 *> Modified if MODE is nonzero.
220 *> \endverbatim
221 *>
222 *> \param[in] MODES
223 *> \verbatim
224 *> MODES is INTEGER
225 *> \endverbatim
226 *>
227 *> \param[in] CONDS
228 *> \verbatim
229 *> CONDS is DOUBLE PRECISION
230 *> Same as MODE and COND, but for specifying the diagonal
231 *> of S. MODES=-6 and +6 are not allowed (since they would
232 *> result in randomly ill-conditioned eigenvalues.)
233 *> \endverbatim
234 *>
235 *> \param[in] KL
236 *> \verbatim
237 *> KL is INTEGER
238 *> This specifies the lower bandwidth of the matrix. KL=1
239 *> specifies upper Hessenberg form. If KL is at least N-1,
240 *> then A will have full lower bandwidth. KL must be at
241 *> least 1.
242 *> Not modified.
243 *> \endverbatim
244 *>
245 *> \param[in] KU
246 *> \verbatim
247 *> KU is INTEGER
248 *> This specifies the upper bandwidth of the matrix. KU=1
249 *> specifies lower Hessenberg form. If KU is at least N-1,
250 *> then A will have full upper bandwidth; if KU and KL
251 *> are both at least N-1, then A will be dense. Only one of
252 *> KU and KL may be less than N-1. KU must be at least 1.
253 *> Not modified.
254 *> \endverbatim
255 *>
256 *> \param[in] ANORM
257 *> \verbatim
258 *> ANORM is DOUBLE PRECISION
259 *> If ANORM is not negative, then A will be scaled by a non-
260 *> negative real number to make the maximum-element-norm of A
261 *> to be ANORM.
262 *> Not modified.
263 *> \endverbatim
264 *>
265 *> \param[out] A
266 *> \verbatim
267 *> A is DOUBLE PRECISION array, dimension ( LDA, N )
268 *> On exit A is the desired test matrix.
269 *> Modified.
270 *> \endverbatim
271 *>
272 *> \param[in] LDA
273 *> \verbatim
274 *> LDA is INTEGER
275 *> LDA specifies the first dimension of A as declared in the
276 *> calling program. LDA must be at least N.
277 *> Not modified.
278 *> \endverbatim
279 *>
280 *> \param[out] WORK
281 *> \verbatim
282 *> WORK is DOUBLE PRECISION array, dimension ( 3*N )
283 *> Workspace.
284 *> Modified.
285 *> \endverbatim
286 *>
287 *> \param[out] INFO
288 *> \verbatim
289 *> INFO is INTEGER
290 *> Error code. On exit, INFO will be set to one of the
291 *> following values:
292 *> 0 => normal return
293 *> -1 => N negative
294 *> -2 => DIST illegal string
295 *> -5 => MODE not in range -6 to 6
296 *> -6 => COND less than 1.0, and MODE neither -6, 0 nor 6
297 *> -8 => EI(1) is not ' ' or 'R', EI(j) is not 'R' or 'I', or
298 *> two adjacent elements of EI are 'I'.
299 *> -9 => RSIGN is not 'T' or 'F'
300 *> -10 => UPPER is not 'T' or 'F'
301 *> -11 => SIM is not 'T' or 'F'
302 *> -12 => MODES=0 and DS has a zero singular value.
303 *> -13 => MODES is not in the range -5 to 5.
304 *> -14 => MODES is nonzero and CONDS is less than 1.
305 *> -15 => KL is less than 1.
306 *> -16 => KU is less than 1, or KL and KU are both less than
307 *> N-1.
308 *> -19 => LDA is less than N.
309 *> 1 => Error return from DLATM1 (computing D)
310 *> 2 => Cannot scale to DMAX (max. eigenvalue is 0)
311 *> 3 => Error return from DLATM1 (computing DS)
312 *> 4 => Error return from DLARGE
313 *> 5 => Zero singular value from DLATM1.
314 *> \endverbatim
315 *
316 * Authors:
317 * ========
318 *
319 *> \author Univ. of Tennessee
320 *> \author Univ. of California Berkeley
321 *> \author Univ. of Colorado Denver
322 *> \author NAG Ltd.
323 *
324 *> \date November 2011
325 *
326 *> \ingroup double_matgen
327 *
328 * =====================================================================
329  SUBROUTINE dlatme( N, DIST, ISEED, D, MODE, COND, DMAX, EI,
330  $ rsign,
331  $ upper, sim, ds, modes, conds, kl, ku, anorm,
332  $ a,
333  $ lda, work, info )
334 *
335 * -- LAPACK computational routine (version 3.4.0) --
336 * -- LAPACK is a software package provided by Univ. of Tennessee, --
337 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
338 * November 2011
339 *
340 * .. Scalar Arguments ..
341  CHARACTER dist, rsign, sim, upper
342  INTEGER info, kl, ku, lda, mode, modes, n
343  DOUBLE PRECISION anorm, cond, conds, dmax
344 * ..
345 * .. Array Arguments ..
346  CHARACTER ei( * )
347  INTEGER iseed( 4 )
348  DOUBLE PRECISION a( lda, * ), d( * ), ds( * ), work( * )
349 * ..
350 *
351 * =====================================================================
352 *
353 * .. Parameters ..
354  DOUBLE PRECISION zero
355  parameter( zero = 0.0d0 )
356  DOUBLE PRECISION one
357  parameter( one = 1.0d0 )
358  DOUBLE PRECISION half
359  parameter( half = 1.0d0 / 2.0d0 )
360 * ..
361 * .. Local Scalars ..
362  LOGICAL badei, bads, useei
363  INTEGER i, ic, icols, idist, iinfo, ir, irows, irsign,
364  $ isim, iupper, j, jc, jcr, jr
365  DOUBLE PRECISION alpha, tau, temp, xnorms
366 * ..
367 * .. Local Arrays ..
368  DOUBLE PRECISION tempa( 1 )
369 * ..
370 * .. External Functions ..
371  LOGICAL lsame
372  DOUBLE PRECISION dlange, dlaran
373  EXTERNAL lsame, dlange, dlaran
374 * ..
375 * .. External Subroutines ..
376  EXTERNAL dcopy, dgemv, dger, dlarfg, dlarge, dlarnv,
378 * ..
379 * .. Intrinsic Functions ..
380  INTRINSIC abs, max, mod
381 * ..
382 * .. Executable Statements ..
383 *
384 * 1) Decode and Test the input parameters.
385 * Initialize flags & seed.
386 *
387  info = 0
388 *
389 * Quick return if possible
390 *
391  IF( n.EQ.0 )
392  $ return
393 *
394 * Decode DIST
395 *
396  IF( lsame( dist, 'U' ) ) THEN
397  idist = 1
398  ELSE IF( lsame( dist, 'S' ) ) THEN
399  idist = 2
400  ELSE IF( lsame( dist, 'N' ) ) THEN
401  idist = 3
402  ELSE
403  idist = -1
404  END IF
405 *
406 * Check EI
407 *
408  useei = .true.
409  badei = .false.
410  IF( lsame( ei( 1 ), ' ' ) .OR. mode.NE.0 ) THEN
411  useei = .false.
412  ELSE
413  IF( lsame( ei( 1 ), 'R' ) ) THEN
414  DO 10 j = 2, n
415  IF( lsame( ei( j ), 'I' ) ) THEN
416  IF( lsame( ei( j-1 ), 'I' ) )
417  $ badei = .true.
418  ELSE
419  IF( .NOT.lsame( ei( j ), 'R' ) )
420  $ badei = .true.
421  END IF
422  10 continue
423  ELSE
424  badei = .true.
425  END IF
426  END IF
427 *
428 * Decode RSIGN
429 *
430  IF( lsame( rsign, 'T' ) ) THEN
431  irsign = 1
432  ELSE IF( lsame( rsign, 'F' ) ) THEN
433  irsign = 0
434  ELSE
435  irsign = -1
436  END IF
437 *
438 * Decode UPPER
439 *
440  IF( lsame( upper, 'T' ) ) THEN
441  iupper = 1
442  ELSE IF( lsame( upper, 'F' ) ) THEN
443  iupper = 0
444  ELSE
445  iupper = -1
446  END IF
447 *
448 * Decode SIM
449 *
450  IF( lsame( sim, 'T' ) ) THEN
451  isim = 1
452  ELSE IF( lsame( sim, 'F' ) ) THEN
453  isim = 0
454  ELSE
455  isim = -1
456  END IF
457 *
458 * Check DS, if MODES=0 and ISIM=1
459 *
460  bads = .false.
461  IF( modes.EQ.0 .AND. isim.EQ.1 ) THEN
462  DO 20 j = 1, n
463  IF( ds( j ).EQ.zero )
464  $ bads = .true.
465  20 continue
466  END IF
467 *
468 * Set INFO if an error
469 *
470  IF( n.LT.0 ) THEN
471  info = -1
472  ELSE IF( idist.EQ.-1 ) THEN
473  info = -2
474  ELSE IF( abs( mode ).GT.6 ) THEN
475  info = -5
476  ELSE IF( ( mode.NE.0 .AND. abs( mode ).NE.6 ) .AND. cond.LT.one )
477  $ THEN
478  info = -6
479  ELSE IF( badei ) THEN
480  info = -8
481  ELSE IF( irsign.EQ.-1 ) THEN
482  info = -9
483  ELSE IF( iupper.EQ.-1 ) THEN
484  info = -10
485  ELSE IF( isim.EQ.-1 ) THEN
486  info = -11
487  ELSE IF( bads ) THEN
488  info = -12
489  ELSE IF( isim.EQ.1 .AND. abs( modes ).GT.5 ) THEN
490  info = -13
491  ELSE IF( isim.EQ.1 .AND. modes.NE.0 .AND. conds.LT.one ) THEN
492  info = -14
493  ELSE IF( kl.LT.1 ) THEN
494  info = -15
495  ELSE IF( ku.LT.1 .OR. ( ku.LT.n-1 .AND. kl.LT.n-1 ) ) THEN
496  info = -16
497  ELSE IF( lda.LT.max( 1, n ) ) THEN
498  info = -19
499  END IF
500 *
501  IF( info.NE.0 ) THEN
502  CALL xerbla( 'DLATME', -info )
503  return
504  END IF
505 *
506 * Initialize random number generator
507 *
508  DO 30 i = 1, 4
509  iseed( i ) = mod( abs( iseed( i ) ), 4096 )
510  30 continue
511 *
512  IF( mod( iseed( 4 ), 2 ).NE.1 )
513  $ iseed( 4 ) = iseed( 4 ) + 1
514 *
515 * 2) Set up diagonal of A
516 *
517 * Compute D according to COND and MODE
518 *
519  CALL dlatm1( mode, cond, irsign, idist, iseed, d, n, iinfo )
520  IF( iinfo.NE.0 ) THEN
521  info = 1
522  return
523  END IF
524  IF( mode.NE.0 .AND. abs( mode ).NE.6 ) THEN
525 *
526 * Scale by DMAX
527 *
528  temp = abs( d( 1 ) )
529  DO 40 i = 2, n
530  temp = max( temp, abs( d( i ) ) )
531  40 continue
532 *
533  IF( temp.GT.zero ) THEN
534  alpha = dmax / temp
535  ELSE IF( dmax.NE.zero ) THEN
536  info = 2
537  return
538  ELSE
539  alpha = zero
540  END IF
541 *
542  CALL dscal( n, alpha, d, 1 )
543 *
544  END IF
545 *
546  CALL dlaset( 'Full', n, n, zero, zero, a, lda )
547  CALL dcopy( n, d, 1, a, lda+1 )
548 *
549 * Set up complex conjugate pairs
550 *
551  IF( mode.EQ.0 ) THEN
552  IF( useei ) THEN
553  DO 50 j = 2, n
554  IF( lsame( ei( j ), 'I' ) ) THEN
555  a( j-1, j ) = a( j, j )
556  a( j, j-1 ) = -a( j, j )
557  a( j, j ) = a( j-1, j-1 )
558  END IF
559  50 continue
560  END IF
561 *
562  ELSE IF( abs( mode ).EQ.5 ) THEN
563 *
564  DO 60 j = 2, n, 2
565  IF( dlaran( iseed ).GT.half ) THEN
566  a( j-1, j ) = a( j, j )
567  a( j, j-1 ) = -a( j, j )
568  a( j, j ) = a( j-1, j-1 )
569  END IF
570  60 continue
571  END IF
572 *
573 * 3) If UPPER='T', set upper triangle of A to random numbers.
574 * (but don't modify the corners of 2x2 blocks.)
575 *
576  IF( iupper.NE.0 ) THEN
577  DO 70 jc = 2, n
578  IF( a( jc-1, jc ).NE.zero ) THEN
579  jr = jc - 2
580  ELSE
581  jr = jc - 1
582  END IF
583  CALL dlarnv( idist, iseed, jr, a( 1, jc ) )
584  70 continue
585  END IF
586 *
587 * 4) If SIM='T', apply similarity transformation.
588 *
589 * -1
590 * Transform is X A X , where X = U S V, thus
591 *
592 * it is U S V A V' (1/S) U'
593 *
594  IF( isim.NE.0 ) THEN
595 *
596 * Compute S (singular values of the eigenvector matrix)
597 * according to CONDS and MODES
598 *
599  CALL dlatm1( modes, conds, 0, 0, iseed, ds, n, iinfo )
600  IF( iinfo.NE.0 ) THEN
601  info = 3
602  return
603  END IF
604 *
605 * Multiply by V and V'
606 *
607  CALL dlarge( n, a, lda, iseed, work, iinfo )
608  IF( iinfo.NE.0 ) THEN
609  info = 4
610  return
611  END IF
612 *
613 * Multiply by S and (1/S)
614 *
615  DO 80 j = 1, n
616  CALL dscal( n, ds( j ), a( j, 1 ), lda )
617  IF( ds( j ).NE.zero ) THEN
618  CALL dscal( n, one / ds( j ), a( 1, j ), 1 )
619  ELSE
620  info = 5
621  return
622  END IF
623  80 continue
624 *
625 * Multiply by U and U'
626 *
627  CALL dlarge( n, a, lda, iseed, work, iinfo )
628  IF( iinfo.NE.0 ) THEN
629  info = 4
630  return
631  END IF
632  END IF
633 *
634 * 5) Reduce the bandwidth.
635 *
636  IF( kl.LT.n-1 ) THEN
637 *
638 * Reduce bandwidth -- kill column
639 *
640  DO 90 jcr = kl + 1, n - 1
641  ic = jcr - kl
642  irows = n + 1 - jcr
643  icols = n + kl - jcr
644 *
645  CALL dcopy( irows, a( jcr, ic ), 1, work, 1 )
646  xnorms = work( 1 )
647  CALL dlarfg( irows, xnorms, work( 2 ), 1, tau )
648  work( 1 ) = one
649 *
650  CALL dgemv( 'T', irows, icols, one, a( jcr, ic+1 ), lda,
651  $ work, 1, zero, work( irows+1 ), 1 )
652  CALL dger( irows, icols, -tau, work, 1, work( irows+1 ), 1,
653  $ a( jcr, ic+1 ), lda )
654 *
655  CALL dgemv( 'N', n, irows, one, a( 1, jcr ), lda, work, 1,
656  $ zero, work( irows+1 ), 1 )
657  CALL dger( n, irows, -tau, work( irows+1 ), 1, work, 1,
658  $ a( 1, jcr ), lda )
659 *
660  a( jcr, ic ) = xnorms
661  CALL dlaset( 'Full', irows-1, 1, zero, zero, a( jcr+1, ic ),
662  $ lda )
663  90 continue
664  ELSE IF( ku.LT.n-1 ) THEN
665 *
666 * Reduce upper bandwidth -- kill a row at a time.
667 *
668  DO 100 jcr = ku + 1, n - 1
669  ir = jcr - ku
670  irows = n + ku - jcr
671  icols = n + 1 - jcr
672 *
673  CALL dcopy( icols, a( ir, jcr ), lda, work, 1 )
674  xnorms = work( 1 )
675  CALL dlarfg( icols, xnorms, work( 2 ), 1, tau )
676  work( 1 ) = one
677 *
678  CALL dgemv( 'N', irows, icols, one, a( ir+1, jcr ), lda,
679  $ work, 1, zero, work( icols+1 ), 1 )
680  CALL dger( irows, icols, -tau, work( icols+1 ), 1, work, 1,
681  $ a( ir+1, jcr ), lda )
682 *
683  CALL dgemv( 'C', icols, n, one, a( jcr, 1 ), lda, work, 1,
684  $ zero, work( icols+1 ), 1 )
685  CALL dger( icols, n, -tau, work, 1, work( icols+1 ), 1,
686  $ a( jcr, 1 ), lda )
687 *
688  a( ir, jcr ) = xnorms
689  CALL dlaset( 'Full', 1, icols-1, zero, zero, a( ir, jcr+1 ),
690  $ lda )
691  100 continue
692  END IF
693 *
694 * Scale the matrix to have norm ANORM
695 *
696  IF( anorm.GE.zero ) THEN
697  temp = dlange( 'M', n, n, a, lda, tempa )
698  IF( temp.GT.zero ) THEN
699  alpha = anorm / temp
700  DO 110 j = 1, n
701  CALL dscal( n, alpha, a( 1, j ), 1 )
702  110 continue
703  END IF
704  END IF
705 *
706  return
707 *
708 * End of DLATME
709 *
710  END