LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
slatm4.f
Go to the documentation of this file.
1 *> \brief \b SLATM4
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 SLATM4( 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 * REAL AMAGN, RCOND, TRIANG
17 * ..
18 * .. Array Arguments ..
19 * INTEGER ISEED( 4 )
20 * REAL A( LDA, * )
21 * ..
22 *
23 *
24 *> \par Purpose:
25 * =============
26 *>
27 *> \verbatim
28 *>
29 *> SLATM4 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 SLARND.)
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 REAL
110 *> The diagonal and subdiagonal entries will be multiplied by
111 *> AMAGN.
112 *> \endverbatim
113 *>
114 *> \param[in] RCOND
115 *> \verbatim
116 *> RCOND is REAL
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 REAL
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 SLATM4 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 REAL 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 single_eig
171 *
172 * =====================================================================
173  SUBROUTINE slatm4( 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  REAL AMAGN, RCOND, TRIANG
183 * ..
184 * .. Array Arguments ..
185  INTEGER ISEED( 4 )
186  REAL A( LDA, * )
187 * ..
188 *
189 * =====================================================================
190 *
191 * .. Parameters ..
192  REAL ZERO, ONE, TWO
193  parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
194  REAL 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  REAL ALPHA, CL, CR, SAFMIN, SL, SR, SV1, SV2, TEMP
201 * ..
202 * .. External Functions ..
203  REAL SLAMCH, SLARAN, SLARND
204  EXTERNAL slamch, slaran, slarnd
205 * ..
206 * .. External Subroutines ..
207  EXTERNAL slaset
208 * ..
209 * .. Intrinsic Functions ..
210  INTRINSIC abs, exp, log, max, min, mod, real, sqrt
211 * ..
212 * .. Executable Statements ..
213 *
214  IF( n.LE.0 )
215  $ RETURN
216  CALL slaset( '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 ) = real( 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 / real( klen-1 ) )
306  DO 150 i = 2, klen
307  a( nz1+i, nz1+i ) = alpha**real( 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 ) / real( klen-1 )
318  DO 170 i = 2, klen
319  a( nz1+i, nz1+i ) = real( 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*slaran( 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 ) = slarnd( 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*real( a( jd, jd ) )
346  230 CONTINUE
347  DO 240 jd = isdb, isde
348  a( jd+1, jd ) = amagn*real( 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( real( a( jd, jd ) ).NE.zero ) THEN
357  IF( slaran( iseed ).GT.half )
358  $ a( jd, jd ) = -a( jd, jd )
359  END IF
360  250 CONTINUE
361  DO 260 jd = isdb, isde
362  IF( real( a( jd+1, jd ) ).NE.zero ) THEN
363  IF( slaran( 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 = slamch( 'S' )
389  DO 290 jd = kbeg, kend - 1, 2
390  IF( slaran( iseed ).GT.half ) THEN
391 *
392 * Rotation on left.
393 *
394  cl = two*slaran( iseed ) - one
395  sl = two*slaran( 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*slaran( iseed ) - one
403  sr = two*slaran( 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*slarnd( 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*slarnd( idist, iseed )
438  310 CONTINUE
439  320 CONTINUE
440  END IF
441 *
442  RETURN
443 *
444 * End of SLATM4
445 *
446  END
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: slaset.f:110
subroutine slatm4(ITYPE, N, NZ1, NZ2, ISIGN, AMAGN, RCOND, TRIANG, IDIST, ISEED, A, LDA)
SLATM4
Definition: slatm4.f:175