LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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 *> \date November 2011
171 *
172 *> \ingroup double_eig
173 *
174 * =====================================================================
175  SUBROUTINE dlatm4( ITYPE, N, NZ1, NZ2, ISIGN, AMAGN, RCOND,
176  $ triang, idist, iseed, a, lda )
177 *
178 * -- LAPACK test routine (version 3.4.0) --
179 * -- LAPACK is a software package provided by Univ. of Tennessee, --
180 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
181 * November 2011
182 *
183 * .. Scalar Arguments ..
184  INTEGER idist, isign, itype, lda, n, nz1, nz2
185  DOUBLE PRECISION amagn, rcond, triang
186 * ..
187 * .. Array Arguments ..
188  INTEGER iseed( 4 )
189  DOUBLE PRECISION a( lda, * )
190 * ..
191 *
192 * =====================================================================
193 *
194 * .. Parameters ..
195  DOUBLE PRECISION zero, one, two
196  parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
197  DOUBLE PRECISION half
198  parameter( half = one / two )
199 * ..
200 * .. Local Scalars ..
201  INTEGER i, ioff, isdb, isde, jc, jd, jr, k, kbeg, kend,
202  $ klen
203  DOUBLE PRECISION alpha, cl, cr, safmin, sl, sr, sv1, sv2, temp
204 * ..
205 * .. External Functions ..
206  DOUBLE PRECISION dlamch, dlaran, dlarnd
207  EXTERNAL dlamch, dlaran, dlarnd
208 * ..
209 * .. External Subroutines ..
210  EXTERNAL dlaset
211 * ..
212 * .. Intrinsic Functions ..
213  INTRINSIC abs, dble, exp, log, max, min, mod, sqrt
214 * ..
215 * .. Executable Statements ..
216 *
217  IF( n.LE.0 )
218  $ return
219  CALL dlaset( 'Full', n, n, zero, zero, a, lda )
220 *
221 * Insure a correct ISEED
222 *
223  IF( mod( iseed( 4 ), 2 ).NE.1 )
224  $ iseed( 4 ) = iseed( 4 ) + 1
225 *
226 * Compute diagonal and subdiagonal according to ITYPE, NZ1, NZ2,
227 * and RCOND
228 *
229  IF( itype.NE.0 ) THEN
230  IF( abs( itype ).GE.4 ) THEN
231  kbeg = max( 1, min( n, nz1+1 ) )
232  kend = max( kbeg, min( n, n-nz2 ) )
233  klen = kend + 1 - kbeg
234  ELSE
235  kbeg = 1
236  kend = n
237  klen = n
238  END IF
239  isdb = 1
240  isde = 0
241  go to( 10, 30, 50, 80, 100, 120, 140, 160,
242  $ 180, 200 )abs( itype )
243 *
244 * abs(ITYPE) = 1: Identity
245 *
246  10 continue
247  DO 20 jd = 1, n
248  a( jd, jd ) = one
249  20 continue
250  go to 220
251 *
252 * abs(ITYPE) = 2: Transposed Jordan block
253 *
254  30 continue
255  DO 40 jd = 1, n - 1
256  a( jd+1, jd ) = one
257  40 continue
258  isdb = 1
259  isde = n - 1
260  go to 220
261 *
262 * abs(ITYPE) = 3: Transposed Jordan block, followed by the
263 * identity.
264 *
265  50 continue
266  k = ( n-1 ) / 2
267  DO 60 jd = 1, k
268  a( jd+1, jd ) = one
269  60 continue
270  isdb = 1
271  isde = k
272  DO 70 jd = k + 2, 2*k + 1
273  a( jd, jd ) = one
274  70 continue
275  go to 220
276 *
277 * abs(ITYPE) = 4: 1,...,k
278 *
279  80 continue
280  DO 90 jd = kbeg, kend
281  a( jd, jd ) = dble( jd-nz1 )
282  90 continue
283  go to 220
284 *
285 * abs(ITYPE) = 5: One large D value:
286 *
287  100 continue
288  DO 110 jd = kbeg + 1, kend
289  a( jd, jd ) = rcond
290  110 continue
291  a( kbeg, kbeg ) = one
292  go to 220
293 *
294 * abs(ITYPE) = 6: One small D value:
295 *
296  120 continue
297  DO 130 jd = kbeg, kend - 1
298  a( jd, jd ) = one
299  130 continue
300  a( kend, kend ) = rcond
301  go to 220
302 *
303 * abs(ITYPE) = 7: Exponentially distributed D values:
304 *
305  140 continue
306  a( kbeg, kbeg ) = one
307  IF( klen.GT.1 ) THEN
308  alpha = rcond**( one / dble( klen-1 ) )
309  DO 150 i = 2, klen
310  a( nz1+i, nz1+i ) = alpha**dble( i-1 )
311  150 continue
312  END IF
313  go to 220
314 *
315 * abs(ITYPE) = 8: Arithmetically distributed D values:
316 *
317  160 continue
318  a( kbeg, kbeg ) = one
319  IF( klen.GT.1 ) THEN
320  alpha = ( one-rcond ) / dble( klen-1 )
321  DO 170 i = 2, klen
322  a( nz1+i, nz1+i ) = dble( klen-i )*alpha + rcond
323  170 continue
324  END IF
325  go to 220
326 *
327 * abs(ITYPE) = 9: Randomly distributed D values on ( RCOND, 1):
328 *
329  180 continue
330  alpha = log( rcond )
331  DO 190 jd = kbeg, kend
332  a( jd, jd ) = exp( alpha*dlaran( iseed ) )
333  190 continue
334  go to 220
335 *
336 * abs(ITYPE) = 10: Randomly distributed D values from DIST
337 *
338  200 continue
339  DO 210 jd = kbeg, kend
340  a( jd, jd ) = dlarnd( idist, iseed )
341  210 continue
342 *
343  220 continue
344 *
345 * Scale by AMAGN
346 *
347  DO 230 jd = kbeg, kend
348  a( jd, jd ) = amagn*dble( a( jd, jd ) )
349  230 continue
350  DO 240 jd = isdb, isde
351  a( jd+1, jd ) = amagn*dble( a( jd+1, jd ) )
352  240 continue
353 *
354 * If ISIGN = 1 or 2, assign random signs to diagonal and
355 * subdiagonal
356 *
357  IF( isign.GT.0 ) THEN
358  DO 250 jd = kbeg, kend
359  IF( dble( a( jd, jd ) ).NE.zero ) THEN
360  IF( dlaran( iseed ).GT.half )
361  $ a( jd, jd ) = -a( jd, jd )
362  END IF
363  250 continue
364  DO 260 jd = isdb, isde
365  IF( dble( a( jd+1, jd ) ).NE.zero ) THEN
366  IF( dlaran( iseed ).GT.half )
367  $ a( jd+1, jd ) = -a( jd+1, jd )
368  END IF
369  260 continue
370  END IF
371 *
372 * Reverse if ITYPE < 0
373 *
374  IF( itype.LT.0 ) THEN
375  DO 270 jd = kbeg, ( kbeg+kend-1 ) / 2
376  temp = a( jd, jd )
377  a( jd, jd ) = a( kbeg+kend-jd, kbeg+kend-jd )
378  a( kbeg+kend-jd, kbeg+kend-jd ) = temp
379  270 continue
380  DO 280 jd = 1, ( n-1 ) / 2
381  temp = a( jd+1, jd )
382  a( jd+1, jd ) = a( n+1-jd, n-jd )
383  a( n+1-jd, n-jd ) = temp
384  280 continue
385  END IF
386 *
387 * If ISIGN = 2, and no subdiagonals already, then apply
388 * random rotations to make 2x2 blocks.
389 *
390  IF( isign.EQ.2 .AND. itype.NE.2 .AND. itype.NE.3 ) THEN
391  safmin = dlamch( 'S' )
392  DO 290 jd = kbeg, kend - 1, 2
393  IF( dlaran( iseed ).GT.half ) THEN
394 *
395 * Rotation on left.
396 *
397  cl = two*dlaran( iseed ) - one
398  sl = two*dlaran( iseed ) - one
399  temp = one / max( safmin, sqrt( cl**2+sl**2 ) )
400  cl = cl*temp
401  sl = sl*temp
402 *
403 * Rotation on right.
404 *
405  cr = two*dlaran( iseed ) - one
406  sr = two*dlaran( iseed ) - one
407  temp = one / max( safmin, sqrt( cr**2+sr**2 ) )
408  cr = cr*temp
409  sr = sr*temp
410 *
411 * Apply
412 *
413  sv1 = a( jd, jd )
414  sv2 = a( jd+1, jd+1 )
415  a( jd, jd ) = cl*cr*sv1 + sl*sr*sv2
416  a( jd+1, jd ) = -sl*cr*sv1 + cl*sr*sv2
417  a( jd, jd+1 ) = -cl*sr*sv1 + sl*cr*sv2
418  a( jd+1, jd+1 ) = sl*sr*sv1 + cl*cr*sv2
419  END IF
420  290 continue
421  END IF
422 *
423  END IF
424 *
425 * Fill in upper triangle (except for 2x2 blocks)
426 *
427  IF( triang.NE.zero ) THEN
428  IF( isign.NE.2 .OR. itype.EQ.2 .OR. itype.EQ.3 ) THEN
429  ioff = 1
430  ELSE
431  ioff = 2
432  DO 300 jr = 1, n - 1
433  IF( a( jr+1, jr ).EQ.zero )
434  $ a( jr, jr+1 ) = triang*dlarnd( idist, iseed )
435  300 continue
436  END IF
437 *
438  DO 320 jc = 2, n
439  DO 310 jr = 1, jc - ioff
440  a( jr, jc ) = triang*dlarnd( idist, iseed )
441  310 continue
442  320 continue
443  END IF
444 *
445  return
446 *
447 * End of DLATM4
448 *
449  END