1       SUBROUTINE DLATM4( ITYPE, N, NZ1, NZ2, ISIGN, AMAGN, RCOND,
  2      $                   TRIANG, IDIST, ISEED, A, LDA )
  3 *
  4 *  -- LAPACK auxiliary test routine (version 3.1) --
  5 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
  6 *     November 2006
  7 *
  8 *     .. Scalar Arguments ..
  9       INTEGER            IDIST, ISIGN, ITYPE, LDA, N, NZ1, NZ2
 10       DOUBLE PRECISION   AMAGN, RCOND, TRIANG
 11 *     ..
 12 *     .. Array Arguments ..
 13       INTEGER            ISEED( 4 )
 14       DOUBLE PRECISION   A( LDA, * )
 15 *     ..
 16 *
 17 *  Purpose
 18 *  =======
 19 *
 20 *  DLATM4 generates basic square matrices, which may later be
 21 *  multiplied by others in order to produce test matrices.  It is
 22 *  intended mainly to be used to test the generalized eigenvalue
 23 *  routines.
 24 *
 25 *  It first generates the diagonal and (possibly) subdiagonal,
 26 *  according to the value of ITYPE, NZ1, NZ2, ISIGN, AMAGN, and RCOND.
 27 *  It then fills in the upper triangle with random numbers, if TRIANG is
 28 *  non-zero.
 29 *
 30 *  Arguments
 31 *  =========
 32 *
 33 *  ITYPE   (input) INTEGER
 34 *          The "type" of matrix on the diagonal and sub-diagonal.
 35 *          If ITYPE < 0, then type abs(ITYPE) is generated and then
 36 *             swapped end for end (A(I,J) := A'(N-J,N-I).)  See also
 37 *             the description of AMAGN and ISIGN.
 38 *
 39 *          Special types:
 40 *          = 0:  the zero matrix.
 41 *          = 1:  the identity.
 42 *          = 2:  a transposed Jordan block.
 43 *          = 3:  If N is odd, then a k+1 x k+1 transposed Jordan block
 44 *                followed by a k x k identity block, where k=(N-1)/2.
 45 *                If N is even, then k=(N-2)/2, and a zero diagonal entry
 46 *                is tacked onto the end.
 47 *
 48 *          Diagonal types.  The diagonal consists of NZ1 zeros, then
 49 *             k=N-NZ1-NZ2 nonzeros.  The subdiagonal is zero.  ITYPE
 50 *             specifies the nonzero diagonal entries as follows:
 51 *          = 4:  1, ..., k
 52 *          = 5:  1, RCOND, ..., RCOND
 53 *          = 6:  1, ..., 1, RCOND
 54 *          = 7:  1, a, a^2, ..., a^(k-1)=RCOND
 55 *          = 8:  1, 1-d, 1-2*d, ..., 1-(k-1)*d=RCOND
 56 *          = 9:  random numbers chosen from (RCOND,1)
 57 *          = 10: random numbers with distribution IDIST (see DLARND.)
 58 *
 59 *  N       (input) INTEGER
 60 *          The order of the matrix.
 61 *
 62 *  NZ1     (input) INTEGER
 63 *          If abs(ITYPE) > 3, then the first NZ1 diagonal entries will
 64 *          be zero.
 65 *
 66 *  NZ2     (input) INTEGER
 67 *          If abs(ITYPE) > 3, then the last NZ2 diagonal entries will
 68 *          be zero.
 69 *
 70 *  ISIGN   (input) INTEGER
 71 *          = 0: The sign of the diagonal and subdiagonal entries will
 72 *               be left unchanged.
 73 *          = 1: The diagonal and subdiagonal entries will have their
 74 *               sign changed at random.
 75 *          = 2: If ITYPE is 2 or 3, then the same as ISIGN=1.
 76 *               Otherwise, with probability 0.5, odd-even pairs of
 77 *               diagonal entries A(2*j-1,2*j-1), A(2*j,2*j) will be
 78 *               converted to a 2x2 block by pre- and post-multiplying
 79 *               by distinct random orthogonal rotations.  The remaining
 80 *               diagonal entries will have their sign changed at random.
 81 *
 82 *  AMAGN   (input) DOUBLE PRECISION
 83 *          The diagonal and subdiagonal entries will be multiplied by
 84 *          AMAGN.
 85 *
 86 *  RCOND   (input) DOUBLE PRECISION
 87 *          If abs(ITYPE) > 4, then the smallest diagonal entry will be
 88 *          entry will be RCOND.  RCOND must be between 0 and 1.
 89 *
 90 *  TRIANG  (input) DOUBLE PRECISION
 91 *          The entries above the diagonal will be random numbers with
 92 *          magnitude bounded by TRIANG (i.e., random numbers multiplied
 93 *          by TRIANG.)
 94 *
 95 *  IDIST   (input) INTEGER
 96 *          Specifies the type of distribution to be used to generate a
 97 *          random matrix.
 98 *          = 1:  UNIFORM( 0, 1 )
 99 *          = 2:  UNIFORM( -1, 1 )
100 *          = 3:  NORMAL ( 0, 1 )
101 *
102 *  ISEED   (input/output) INTEGER array, dimension (4)
103 *          On entry ISEED specifies the seed of the random number
104 *          generator.  The values of ISEED are changed on exit, and can
105 *          be used in the next call to DLATM4 to continue the same
106 *          random number sequence.
107 *          Note: ISEED(4) should be odd, for the random number generator
108 *          used at present.
109 *
110 *  A       (output) DOUBLE PRECISION array, dimension (LDA, N)
111 *          Array to be computed.
112 *
113 *  LDA     (input) INTEGER
114 *          Leading dimension of A.  Must be at least 1 and at least N.
115 *
116 *  =====================================================================
117 *
118 *     .. Parameters ..
119       DOUBLE PRECISION   ZERO, ONE, TWO
120       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
121       DOUBLE PRECISION   HALF
122       PARAMETER          ( HALF = ONE / TWO )
123 *     ..
124 *     .. Local Scalars ..
125       INTEGER            I, IOFF, ISDB, ISDE, JC, JD, JR, K, KBEG, KEND,
126      $                   KLEN
127       DOUBLE PRECISION   ALPHA, CL, CR, SAFMIN, SL, SR, SV1, SV2, TEMP
128 *     ..
129 *     .. External Functions ..
130       DOUBLE PRECISION   DLAMCH, DLARAN, DLARND
131       EXTERNAL           DLAMCH, DLARAN, DLARND
132 *     ..
133 *     .. External Subroutines ..
134       EXTERNAL           DLASET
135 *     ..
136 *     .. Intrinsic Functions ..
137       INTRINSIC          ABSDBLEEXPLOGMAXMINMODSQRT
138 *     ..
139 *     .. Executable Statements ..
140 *
141       IF( N.LE.0 )
142      $   RETURN
143       CALL DLASET( 'Full', N, N, ZERO, ZERO, A, LDA )
144 *
145 *     Insure a correct ISEED
146 *
147       IFMOD( ISEED( 4 ), 2 ).NE.1 )
148      $   ISEED( 4 ) = ISEED( 4 ) + 1
149 *
150 *     Compute diagonal and subdiagonal according to ITYPE, NZ1, NZ2,
151 *     and RCOND
152 *
153       IF( ITYPE.NE.0 ) THEN
154          IFABS( ITYPE ).GE.4 ) THEN
155             KBEG = MAX1MIN( N, NZ1+1 ) )
156             KEND = MAX( KBEG, MIN( N, N-NZ2 ) )
157             KLEN = KEND + 1 - KBEG
158          ELSE
159             KBEG = 1
160             KEND = N
161             KLEN = N
162          END IF
163          ISDB = 1
164          ISDE = 0
165          GO TO ( 10305080100120140160,
166      $           180200 )ABS( ITYPE )
167 *
168 *        abs(ITYPE) = 1: Identity
169 *
170    10    CONTINUE
171          DO 20 JD = 1, N
172             A( JD, JD ) = ONE
173    20    CONTINUE
174          GO TO 220
175 *
176 *        abs(ITYPE) = 2: Transposed Jordan block
177 *
178    30    CONTINUE
179          DO 40 JD = 1, N - 1
180             A( JD+1, JD ) = ONE
181    40    CONTINUE
182          ISDB = 1
183          ISDE = N - 1
184          GO TO 220
185 *
186 *        abs(ITYPE) = 3: Transposed Jordan block, followed by the
187 *                        identity.
188 *
189    50    CONTINUE
190          K = ( N-1 ) / 2
191          DO 60 JD = 1, K
192             A( JD+1, JD ) = ONE
193    60    CONTINUE
194          ISDB = 1
195          ISDE = K
196          DO 70 JD = K + 22*+ 1
197             A( JD, JD ) = ONE
198    70    CONTINUE
199          GO TO 220
200 *
201 *        abs(ITYPE) = 4: 1,...,k
202 *
203    80    CONTINUE
204          DO 90 JD = KBEG, KEND
205             A( JD, JD ) = DBLE( JD-NZ1 )
206    90    CONTINUE
207          GO TO 220
208 *
209 *        abs(ITYPE) = 5: One large D value:
210 *
211   100    CONTINUE
212          DO 110 JD = KBEG + 1, KEND
213             A( JD, JD ) = RCOND
214   110    CONTINUE
215          A( KBEG, KBEG ) = ONE
216          GO TO 220
217 *
218 *        abs(ITYPE) = 6: One small D value:
219 *
220   120    CONTINUE
221          DO 130 JD = KBEG, KEND - 1
222             A( JD, JD ) = ONE
223   130    CONTINUE
224          A( KEND, KEND ) = RCOND
225          GO TO 220
226 *
227 *        abs(ITYPE) = 7: Exponentially distributed D values:
228 *
229   140    CONTINUE
230          A( KBEG, KBEG ) = ONE
231          IF( KLEN.GT.1 ) THEN
232             ALPHA = RCOND**( ONE / DBLE( KLEN-1 ) )
233             DO 150 I = 2, KLEN
234                A( NZ1+I, NZ1+I ) = ALPHA**DBLE( I-1 )
235   150       CONTINUE
236          END IF
237          GO TO 220
238 *
239 *        abs(ITYPE) = 8: Arithmetically distributed D values:
240 *
241   160    CONTINUE
242          A( KBEG, KBEG ) = ONE
243          IF( KLEN.GT.1 ) THEN
244             ALPHA = ( ONE-RCOND ) / DBLE( KLEN-1 )
245             DO 170 I = 2, KLEN
246                A( NZ1+I, NZ1+I ) = DBLE( KLEN-I )*ALPHA + RCOND
247   170       CONTINUE
248          END IF
249          GO TO 220
250 *
251 *        abs(ITYPE) = 9: Randomly distributed D values on ( RCOND, 1):
252 *
253   180    CONTINUE
254          ALPHA = LOG( RCOND )
255          DO 190 JD = KBEG, KEND
256             A( JD, JD ) = EXP( ALPHA*DLARAN( ISEED ) )
257   190    CONTINUE
258          GO TO 220
259 *
260 *        abs(ITYPE) = 10: Randomly distributed D values from DIST
261 *
262   200    CONTINUE
263          DO 210 JD = KBEG, KEND
264             A( JD, JD ) = DLARND( IDIST, ISEED )
265   210    CONTINUE
266 *
267   220    CONTINUE
268 *
269 *        Scale by AMAGN
270 *
271          DO 230 JD = KBEG, KEND
272             A( JD, JD ) = AMAGN*DBLE( A( JD, JD ) )
273   230    CONTINUE
274          DO 240 JD = ISDB, ISDE
275             A( JD+1, JD ) = AMAGN*DBLE( A( JD+1, JD ) )
276   240    CONTINUE
277 *
278 *        If ISIGN = 1 or 2, assign random signs to diagonal and
279 *        subdiagonal
280 *
281          IFISIGN.GT.0 ) THEN
282             DO 250 JD = KBEG, KEND
283                IFDBLE( A( JD, JD ) ).NE.ZERO ) THEN
284                   IF( DLARAN( ISEED ).GT.HALF )
285      $               A( JD, JD ) = -A( JD, JD )
286                END IF
287   250       CONTINUE
288             DO 260 JD = ISDB, ISDE
289                IFDBLE( A( JD+1, JD ) ).NE.ZERO ) THEN
290                   IF( DLARAN( ISEED ).GT.HALF )
291      $               A( JD+1, JD ) = -A( JD+1, JD )
292                END IF
293   260       CONTINUE
294          END IF
295 *
296 *        Reverse if ITYPE < 0
297 *
298          IF( ITYPE.LT.0 ) THEN
299             DO 270 JD = KBEG, ( KBEG+KEND-1 ) / 2
300                TEMP = A( JD, JD )
301                A( JD, JD ) = A( KBEG+KEND-JD, KBEG+KEND-JD )
302                A( KBEG+KEND-JD, KBEG+KEND-JD ) = TEMP
303   270       CONTINUE
304             DO 280 JD = 1, ( N-1 ) / 2
305                TEMP = A( JD+1, JD )
306                A( JD+1, JD ) = A( N+1-JD, N-JD )
307                A( N+1-JD, N-JD ) = TEMP
308   280       CONTINUE
309          END IF
310 *
311 *        If ISIGN = 2, and no subdiagonals already, then apply
312 *        random rotations to make 2x2 blocks.
313 *
314          IFISIGN.EQ.2 .AND. ITYPE.NE.2 .AND. ITYPE.NE.3 ) THEN
315             SAFMIN = DLAMCH( 'S' )
316             DO 290 JD = KBEG, KEND - 12
317                IF( DLARAN( ISEED ).GT.HALF ) THEN
318 *
319 *                 Rotation on left.
320 *
321                   CL = TWO*DLARAN( ISEED ) - ONE
322                   SL = TWO*DLARAN( ISEED ) - ONE
323                   TEMP = ONE / MAX( SAFMIN, SQRT( CL**2+SL**2 ) )
324                   CL = CL*TEMP
325                   SL = SL*TEMP
326 *
327 *                 Rotation on right.
328 *
329                   CR = TWO*DLARAN( ISEED ) - ONE
330                   SR = TWO*DLARAN( ISEED ) - ONE
331                   TEMP = ONE / MAX( SAFMIN, SQRT( CR**2+SR**2 ) )
332                   CR = CR*TEMP
333                   SR = SR*TEMP
334 *
335 *                 Apply
336 *
337                   SV1 = A( JD, JD )
338                   SV2 = A( JD+1, JD+1 )
339                   A( JD, JD ) = CL*CR*SV1 + SL*SR*SV2
340                   A( JD+1, JD ) = -SL*CR*SV1 + CL*SR*SV2
341                   A( JD, JD+1 ) = -CL*SR*SV1 + SL*CR*SV2
342                   A( JD+1, JD+1 ) = SL*SR*SV1 + CL*CR*SV2
343                END IF
344   290       CONTINUE
345          END IF
346 *
347       END IF
348 *
349 *     Fill in upper triangle (except for 2x2 blocks)
350 *
351       IF( TRIANG.NE.ZERO ) THEN
352          IFISIGN.NE.2 .OR. ITYPE.EQ.2 .OR. ITYPE.EQ.3 ) THEN
353             IOFF = 1
354          ELSE
355             IOFF = 2
356             DO 300 JR = 1, N - 1
357                IF( A( JR+1, JR ).EQ.ZERO )
358      $            A( JR, JR+1 ) = TRIANG*DLARND( IDIST, ISEED )
359   300       CONTINUE
360          END IF
361 *
362          DO 320 JC = 2, N
363             DO 310 JR = 1, JC - IOFF
364                A( JR, JC ) = TRIANG*DLARND( IDIST, ISEED )
365   310       CONTINUE
366   320    CONTINUE
367       END IF
368 *
369       RETURN
370 *
371 *     End of DLATM4
372 *
373       END