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*> \date December 2016
171*
172*> \ingroup single_eig
173*
174*  =====================================================================
175      SUBROUTINE SLATM4( ITYPE, N, NZ1, NZ2, ISIGN, AMAGN, RCOND,
176     $                   TRIANG, IDIST, ISEED, A, LDA )
177*
178*  -- LAPACK test routine (version 3.7.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*     December 2016
182*
183*     .. Scalar Arguments ..
184      INTEGER            IDIST, ISIGN, ITYPE, LDA, N, NZ1, NZ2
185      REAL               AMAGN, RCOND, TRIANG
186*     ..
187*     .. Array Arguments ..
188      INTEGER            ISEED( 4 )
189      REAL               A( LDA, * )
190*     ..
191*
192*  =====================================================================
193*
194*     .. Parameters ..
195      REAL               ZERO, ONE, TWO
196      PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0 )
197      REAL               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      REAL               ALPHA, CL, CR, SAFMIN, SL, SR, SV1, SV2, TEMP
204*     ..
205*     .. External Functions ..
206      REAL               SLAMCH, SLARAN, SLARND
207      EXTERNAL           SLAMCH, SLARAN, SLARND
208*     ..
209*     .. External Subroutines ..
210      EXTERNAL           SLASET
211*     ..
212*     .. Intrinsic Functions ..
213      INTRINSIC          ABS, EXP, LOG, MAX, MIN, MOD, REAL, SQRT
214*     ..
215*     .. Executable Statements ..
216*
217      IF( N.LE.0 )
218     $   RETURN
219      CALL SLASET( '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 ) = REAL( 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 / REAL( KLEN-1 ) )
309            DO 150 I = 2, KLEN
310               A( NZ1+I, NZ1+I ) = ALPHA**REAL( 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 ) / REAL( KLEN-1 )
321            DO 170 I = 2, KLEN
322               A( NZ1+I, NZ1+I ) = REAL( 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*SLARAN( 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 ) = SLARND( 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*REAL( A( JD, JD ) )
349  230    CONTINUE
350         DO 240 JD = ISDB, ISDE
351            A( JD+1, JD ) = AMAGN*REAL( 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( REAL( A( JD, JD ) ).NE.ZERO ) THEN
360                  IF( SLARAN( ISEED ).GT.HALF )
361     $               A( JD, JD ) = -A( JD, JD )
362               END IF
363  250       CONTINUE
364            DO 260 JD = ISDB, ISDE
365               IF( REAL( A( JD+1, JD ) ).NE.ZERO ) THEN
366                  IF( SLARAN( 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 = SLAMCH( 'S' )
392            DO 290 JD = KBEG, KEND - 1, 2
393               IF( SLARAN( ISEED ).GT.HALF ) THEN
394*
395*                 Rotation on left.
396*
397                  CL = TWO*SLARAN( ISEED ) - ONE
398                  SL = TWO*SLARAN( 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*SLARAN( ISEED ) - ONE
406                  SR = TWO*SLARAN( 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*SLARND( 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*SLARND( IDIST, ISEED )
441  310       CONTINUE
442  320    CONTINUE
443      END IF
444*
445      RETURN
446*
447*     End of SLATM4
448*
449      END
450