LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dlatm4.f
Go to the documentation of this file.
1*> \brief \b DLATM4
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 DLATM4( ITYPE, N, NZ1, NZ2, ISIGN, AMAGN, RCOND,
12* TRIANG, IDIST, ISEED, A, LDA )
13*
14* .. Scalar Arguments ..
15* INTEGER IDIST, ISIGN, ITYPE, LDA, N, NZ1, NZ2
16* DOUBLE PRECISION AMAGN, RCOND, TRIANG
17* ..
18* .. Array Arguments ..
19* INTEGER ISEED( 4 )
20* DOUBLE PRECISION A( LDA, * )
21* ..
22*
23*
24*> \par Purpose:
25* =============
26*>
27*> \verbatim
28*>
29*> DLATM4 generates basic square matrices, which may later be
30*> multiplied by others in order to produce test matrices. It is
31*> intended mainly to be used to test the generalized eigenvalue
32*> routines.
33*>
34*> It first generates the diagonal and (possibly) subdiagonal,
35*> according to the value of ITYPE, NZ1, NZ2, ISIGN, AMAGN, and RCOND.
36*> It then fills in the upper triangle with random numbers, if TRIANG is
37*> non-zero.
38*> \endverbatim
39*
40* Arguments:
41* ==========
42*
43*> \param[in] ITYPE
44*> \verbatim
45*> ITYPE is INTEGER
46*> The "type" of matrix on the diagonal and sub-diagonal.
47*> If ITYPE < 0, then type abs(ITYPE) is generated and then
48*> swapped end for end (A(I,J) := A'(N-J,N-I).) See also
49*> the description of AMAGN and ISIGN.
50*>
51*> Special types:
52*> = 0: the zero matrix.
53*> = 1: the identity.
54*> = 2: a transposed Jordan block.
55*> = 3: If N is odd, then a k+1 x k+1 transposed Jordan block
56*> followed by a k x k identity block, where k=(N-1)/2.
57*> If N is even, then k=(N-2)/2, and a zero diagonal entry
58*> is tacked onto the end.
59*>
60*> Diagonal types. The diagonal consists of NZ1 zeros, then
61*> k=N-NZ1-NZ2 nonzeros. The subdiagonal is zero. ITYPE
62*> specifies the nonzero diagonal entries as follows:
63*> = 4: 1, ..., k
64*> = 5: 1, RCOND, ..., RCOND
65*> = 6: 1, ..., 1, RCOND
66*> = 7: 1, a, a^2, ..., a^(k-1)=RCOND
67*> = 8: 1, 1-d, 1-2*d, ..., 1-(k-1)*d=RCOND
68*> = 9: random numbers chosen from (RCOND,1)
69*> = 10: random numbers with distribution IDIST (see DLARND.)
70*> \endverbatim
71*>
72*> \param[in] N
73*> \verbatim
74*> N is INTEGER
75*> The order of the matrix.
76*> \endverbatim
77*>
78*> \param[in] NZ1
79*> \verbatim
80*> NZ1 is INTEGER
81*> If abs(ITYPE) > 3, then the first NZ1 diagonal entries will
82*> be zero.
83*> \endverbatim
84*>
85*> \param[in] NZ2
86*> \verbatim
87*> NZ2 is INTEGER
88*> If abs(ITYPE) > 3, then the last NZ2 diagonal entries will
89*> be zero.
90*> \endverbatim
91*>
92*> \param[in] ISIGN
93*> \verbatim
94*> ISIGN is INTEGER
95*> = 0: The sign of the diagonal and subdiagonal entries will
96*> be left unchanged.
97*> = 1: The diagonal and subdiagonal entries will have their
98*> sign changed at random.
99*> = 2: If ITYPE is 2 or 3, then the same as ISIGN=1.
100*> Otherwise, with probability 0.5, odd-even pairs of
101*> diagonal entries A(2*j-1,2*j-1), A(2*j,2*j) will be
102*> converted to a 2x2 block by pre- and post-multiplying
103*> by distinct random orthogonal rotations. The remaining
104*> diagonal entries will have their sign changed at random.
105*> \endverbatim
106*>
107*> \param[in] AMAGN
108*> \verbatim
109*> AMAGN is DOUBLE PRECISION
110*> The diagonal and subdiagonal entries will be multiplied by
111*> AMAGN.
112*> \endverbatim
113*>
114*> \param[in] RCOND
115*> \verbatim
116*> RCOND is DOUBLE PRECISION
117*> If abs(ITYPE) > 4, then the smallest diagonal entry will be
118*> entry will be RCOND. RCOND must be between 0 and 1.
119*> \endverbatim
120*>
121*> \param[in] TRIANG
122*> \verbatim
123*> TRIANG is DOUBLE PRECISION
124*> The entries above the diagonal will be random numbers with
125*> magnitude bounded by TRIANG (i.e., random numbers multiplied
126*> by TRIANG.)
127*> \endverbatim
128*>
129*> \param[in] IDIST
130*> \verbatim
131*> IDIST is INTEGER
132*> Specifies the type of distribution to be used to generate a
133*> random matrix.
134*> = 1: UNIFORM( 0, 1 )
135*> = 2: UNIFORM( -1, 1 )
136*> = 3: NORMAL ( 0, 1 )
137*> \endverbatim
138*>
139*> \param[in,out] ISEED
140*> \verbatim
141*> ISEED is INTEGER array, dimension (4)
142*> On entry ISEED specifies the seed of the random number
143*> generator. The values of ISEED are changed on exit, and can
144*> be used in the next call to DLATM4 to continue the same
145*> random number sequence.
146*> Note: ISEED(4) should be odd, for the random number generator
147*> used at present.
148*> \endverbatim
149*>
150*> \param[out] A
151*> \verbatim
152*> A is DOUBLE PRECISION array, dimension (LDA, N)
153*> Array to be computed.
154*> \endverbatim
155*>
156*> \param[in] LDA
157*> \verbatim
158*> LDA is INTEGER
159*> Leading dimension of A. Must be at least 1 and at least N.
160*> \endverbatim
161*
162* Authors:
163* ========
164*
165*> \author Univ. of Tennessee
166*> \author Univ. of California Berkeley
167*> \author Univ. of Colorado Denver
168*> \author NAG Ltd.
169*
170*> \ingroup double_eig
171*
172* =====================================================================
173 SUBROUTINE dlatm4( ITYPE, N, NZ1, NZ2, ISIGN, AMAGN, RCOND,
174 $ TRIANG, IDIST, ISEED, A, LDA )
175*
176* -- LAPACK test routine --
177* -- LAPACK is a software package provided by Univ. of Tennessee, --
178* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
179*
180* .. Scalar Arguments ..
181 INTEGER IDIST, ISIGN, ITYPE, LDA, N, NZ1, NZ2
182 DOUBLE PRECISION AMAGN, RCOND, TRIANG
183* ..
184* .. Array Arguments ..
185 INTEGER ISEED( 4 )
186 DOUBLE PRECISION A( LDA, * )
187* ..
188*
189* =====================================================================
190*
191* .. Parameters ..
192 DOUBLE PRECISION ZERO, ONE, TWO
193 parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
194 DOUBLE PRECISION HALF
195 parameter( half = one / two )
196* ..
197* .. Local Scalars ..
198 INTEGER I, IOFF, ISDB, ISDE, JC, JD, JR, K, KBEG, KEND,
199 $ klen
200 DOUBLE PRECISION ALPHA, CL, CR, SAFMIN, SL, SR, SV1, SV2, TEMP
201* ..
202* .. External Functions ..
203 DOUBLE PRECISION DLAMCH, DLARAN, DLARND
204 EXTERNAL dlamch, dlaran, dlarnd
205* ..
206* .. External Subroutines ..
207 EXTERNAL dlaset
208* ..
209* .. Intrinsic Functions ..
210 INTRINSIC abs, dble, exp, log, max, min, mod, sqrt
211* ..
212* .. Executable Statements ..
213*
214 IF( n.LE.0 )
215 $ RETURN
216 CALL dlaset( 'Full', n, n, zero, zero, a, lda )
217*
218* Insure a correct ISEED
219*
220 IF( mod( iseed( 4 ), 2 ).NE.1 )
221 $ iseed( 4 ) = iseed( 4 ) + 1
222*
223* Compute diagonal and subdiagonal according to ITYPE, NZ1, NZ2,
224* and RCOND
225*
226 IF( itype.NE.0 ) THEN
227 IF( abs( itype ).GE.4 ) THEN
228 kbeg = max( 1, min( n, nz1+1 ) )
229 kend = max( kbeg, min( n, n-nz2 ) )
230 klen = kend + 1 - kbeg
231 ELSE
232 kbeg = 1
233 kend = n
234 klen = n
235 END IF
236 isdb = 1
237 isde = 0
238 GO TO ( 10, 30, 50, 80, 100, 120, 140, 160,
239 $ 180, 200 )abs( itype )
240*
241* abs(ITYPE) = 1: Identity
242*
243 10 CONTINUE
244 DO 20 jd = 1, n
245 a( jd, jd ) = one
246 20 CONTINUE
247 GO TO 220
248*
249* abs(ITYPE) = 2: Transposed Jordan block
250*
251 30 CONTINUE
252 DO 40 jd = 1, n - 1
253 a( jd+1, jd ) = one
254 40 CONTINUE
255 isdb = 1
256 isde = n - 1
257 GO TO 220
258*
259* abs(ITYPE) = 3: Transposed Jordan block, followed by the
260* identity.
261*
262 50 CONTINUE
263 k = ( n-1 ) / 2
264 DO 60 jd = 1, k
265 a( jd+1, jd ) = one
266 60 CONTINUE
267 isdb = 1
268 isde = k
269 DO 70 jd = k + 2, 2*k + 1
270 a( jd, jd ) = one
271 70 CONTINUE
272 GO TO 220
273*
274* abs(ITYPE) = 4: 1,...,k
275*
276 80 CONTINUE
277 DO 90 jd = kbeg, kend
278 a( jd, jd ) = dble( jd-nz1 )
279 90 CONTINUE
280 GO TO 220
281*
282* abs(ITYPE) = 5: One large D value:
283*
284 100 CONTINUE
285 DO 110 jd = kbeg + 1, kend
286 a( jd, jd ) = rcond
287 110 CONTINUE
288 a( kbeg, kbeg ) = one
289 GO TO 220
290*
291* abs(ITYPE) = 6: One small D value:
292*
293 120 CONTINUE
294 DO 130 jd = kbeg, kend - 1
295 a( jd, jd ) = one
296 130 CONTINUE
297 a( kend, kend ) = rcond
298 GO TO 220
299*
300* abs(ITYPE) = 7: Exponentially distributed D values:
301*
302 140 CONTINUE
303 a( kbeg, kbeg ) = one
304 IF( klen.GT.1 ) THEN
305 alpha = rcond**( one / dble( klen-1 ) )
306 DO 150 i = 2, klen
307 a( nz1+i, nz1+i ) = alpha**dble( i-1 )
308 150 CONTINUE
309 END IF
310 GO TO 220
311*
312* abs(ITYPE) = 8: Arithmetically distributed D values:
313*
314 160 CONTINUE
315 a( kbeg, kbeg ) = one
316 IF( klen.GT.1 ) THEN
317 alpha = ( one-rcond ) / dble( klen-1 )
318 DO 170 i = 2, klen
319 a( nz1+i, nz1+i ) = dble( klen-i )*alpha + rcond
320 170 CONTINUE
321 END IF
322 GO TO 220
323*
324* abs(ITYPE) = 9: Randomly distributed D values on ( RCOND, 1):
325*
326 180 CONTINUE
327 alpha = log( rcond )
328 DO 190 jd = kbeg, kend
329 a( jd, jd ) = exp( alpha*dlaran( iseed ) )
330 190 CONTINUE
331 GO TO 220
332*
333* abs(ITYPE) = 10: Randomly distributed D values from DIST
334*
335 200 CONTINUE
336 DO 210 jd = kbeg, kend
337 a( jd, jd ) = dlarnd( idist, iseed )
338 210 CONTINUE
339*
340 220 CONTINUE
341*
342* Scale by AMAGN
343*
344 DO 230 jd = kbeg, kend
345 a( jd, jd ) = amagn*dble( a( jd, jd ) )
346 230 CONTINUE
347 DO 240 jd = isdb, isde
348 a( jd+1, jd ) = amagn*dble( a( jd+1, jd ) )
349 240 CONTINUE
350*
351* If ISIGN = 1 or 2, assign random signs to diagonal and
352* subdiagonal
353*
354 IF( isign.GT.0 ) THEN
355 DO 250 jd = kbeg, kend
356 IF( dble( a( jd, jd ) ).NE.zero ) THEN
357 IF( dlaran( iseed ).GT.half )
358 $ a( jd, jd ) = -a( jd, jd )
359 END IF
360 250 CONTINUE
361 DO 260 jd = isdb, isde
362 IF( dble( a( jd+1, jd ) ).NE.zero ) THEN
363 IF( dlaran( iseed ).GT.half )
364 $ a( jd+1, jd ) = -a( jd+1, jd )
365 END IF
366 260 CONTINUE
367 END IF
368*
369* Reverse if ITYPE < 0
370*
371 IF( itype.LT.0 ) THEN
372 DO 270 jd = kbeg, ( kbeg+kend-1 ) / 2
373 temp = a( jd, jd )
374 a( jd, jd ) = a( kbeg+kend-jd, kbeg+kend-jd )
375 a( kbeg+kend-jd, kbeg+kend-jd ) = temp
376 270 CONTINUE
377 DO 280 jd = 1, ( n-1 ) / 2
378 temp = a( jd+1, jd )
379 a( jd+1, jd ) = a( n+1-jd, n-jd )
380 a( n+1-jd, n-jd ) = temp
381 280 CONTINUE
382 END IF
383*
384* If ISIGN = 2, and no subdiagonals already, then apply
385* random rotations to make 2x2 blocks.
386*
387 IF( isign.EQ.2 .AND. itype.NE.2 .AND. itype.NE.3 ) THEN
388 safmin = dlamch( 'S' )
389 DO 290 jd = kbeg, kend - 1, 2
390 IF( dlaran( iseed ).GT.half ) THEN
391*
392* Rotation on left.
393*
394 cl = two*dlaran( iseed ) - one
395 sl = two*dlaran( iseed ) - one
396 temp = one / max( safmin, sqrt( cl**2+sl**2 ) )
397 cl = cl*temp
398 sl = sl*temp
399*
400* Rotation on right.
401*
402 cr = two*dlaran( iseed ) - one
403 sr = two*dlaran( iseed ) - one
404 temp = one / max( safmin, sqrt( cr**2+sr**2 ) )
405 cr = cr*temp
406 sr = sr*temp
407*
408* Apply
409*
410 sv1 = a( jd, jd )
411 sv2 = a( jd+1, jd+1 )
412 a( jd, jd ) = cl*cr*sv1 + sl*sr*sv2
413 a( jd+1, jd ) = -sl*cr*sv1 + cl*sr*sv2
414 a( jd, jd+1 ) = -cl*sr*sv1 + sl*cr*sv2
415 a( jd+1, jd+1 ) = sl*sr*sv1 + cl*cr*sv2
416 END IF
417 290 CONTINUE
418 END IF
419*
420 END IF
421*
422* Fill in upper triangle (except for 2x2 blocks)
423*
424 IF( triang.NE.zero ) THEN
425 IF( isign.NE.2 .OR. itype.EQ.2 .OR. itype.EQ.3 ) THEN
426 ioff = 1
427 ELSE
428 ioff = 2
429 DO 300 jr = 1, n - 1
430 IF( a( jr+1, jr ).EQ.zero )
431 $ a( jr, jr+1 ) = triang*dlarnd( idist, iseed )
432 300 CONTINUE
433 END IF
434*
435 DO 320 jc = 2, n
436 DO 310 jr = 1, jc - ioff
437 a( jr, jc ) = triang*dlarnd( idist, iseed )
438 310 CONTINUE
439 320 CONTINUE
440 END IF
441*
442 RETURN
443*
444* End of DLATM4
445*
446 END
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: dlaset.f:110
subroutine dlatm4(ITYPE, N, NZ1, NZ2, ISIGN, AMAGN, RCOND, TRIANG, IDIST, ISEED, A, LDA)
DLATM4
Definition: dlatm4.f:175