1 SUBROUTINE ZLATM5( PRTYPE, M, N, A, LDA, B, LDB, C, LDC, D, LDD,
2 $ E, LDE, F, LDF, R, LDR, L, LDL, ALPHA, QBLCKA,
3 $ QBLCKB )
4 *
5 * -- LAPACK test routine (version 3.1) --
6 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
7 * November 2006
8 *
9 * .. Scalar Arguments ..
10 INTEGER LDA, LDB, LDC, LDD, LDE, LDF, LDL, LDR, M, N,
11 $ PRTYPE, QBLCKA, QBLCKB
12 DOUBLE PRECISION ALPHA
13 * ..
14 * .. Array Arguments ..
15 COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * ),
16 $ D( LDD, * ), E( LDE, * ), F( LDF, * ),
17 $ L( LDL, * ), R( LDR, * )
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 * ZLATM5 generates matrices involved in the Generalized Sylvester
24 * equation:
25 *
26 * A * R - L * B = C
27 * D * R - L * E = F
28 *
29 * They also satisfy (the diagonalization condition)
30 *
31 * [ I -L ] ( [ A -C ], [ D -F ] ) [ I R ] = ( [ A ], [ D ] )
32 * [ I ] ( [ B ] [ E ] ) [ I ] ( [ B ] [ E ] )
33 *
34 *
35 * Arguments
36 * =========
37 *
38 * PRTYPE (input) INTEGER
39 * "Points" to a certian type of the matrices to generate
40 * (see futher details).
41 *
42 * M (input) INTEGER
43 * Specifies the order of A and D and the number of rows in
44 * C, F, R and L.
45 *
46 * N (input) INTEGER
47 * Specifies the order of B and E and the number of columns in
48 * C, F, R and L.
49 *
50 * A (output) COMPLEX*16 array, dimension (LDA, M).
51 * On exit A M-by-M is initialized according to PRTYPE.
52 *
53 * LDA (input) INTEGER
54 * The leading dimension of A.
55 *
56 * B (output) COMPLEX*16 array, dimension (LDB, N).
57 * On exit B N-by-N is initialized according to PRTYPE.
58 *
59 * LDB (input) INTEGER
60 * The leading dimension of B.
61 *
62 * C (output) COMPLEX*16 array, dimension (LDC, N).
63 * On exit C M-by-N is initialized according to PRTYPE.
64 *
65 * LDC (input) INTEGER
66 * The leading dimension of C.
67 *
68 * D (output) COMPLEX*16 array, dimension (LDD, M).
69 * On exit D M-by-M is initialized according to PRTYPE.
70 *
71 * LDD (input) INTEGER
72 * The leading dimension of D.
73 *
74 * E (output) COMPLEX*16 array, dimension (LDE, N).
75 * On exit E N-by-N is initialized according to PRTYPE.
76 *
77 * LDE (input) INTEGER
78 * The leading dimension of E.
79 *
80 * F (output) COMPLEX*16 array, dimension (LDF, N).
81 * On exit F M-by-N is initialized according to PRTYPE.
82 *
83 * LDF (input) INTEGER
84 * The leading dimension of F.
85 *
86 * R (output) COMPLEX*16 array, dimension (LDR, N).
87 * On exit R M-by-N is initialized according to PRTYPE.
88 *
89 * LDR (input) INTEGER
90 * The leading dimension of R.
91 *
92 * L (output) COMPLEX*16 array, dimension (LDL, N).
93 * On exit L M-by-N is initialized according to PRTYPE.
94 *
95 * LDL (input) INTEGER
96 * The leading dimension of L.
97 *
98 * ALPHA (input) DOUBLE PRECISION
99 * Parameter used in generating PRTYPE = 1 and 5 matrices.
100 *
101 * QBLCKA (input) INTEGER
102 * When PRTYPE = 3, specifies the distance between 2-by-2
103 * blocks on the diagonal in A. Otherwise, QBLCKA is not
104 * referenced. QBLCKA > 1.
105 *
106 * QBLCKB (input) INTEGER
107 * When PRTYPE = 3, specifies the distance between 2-by-2
108 * blocks on the diagonal in B. Otherwise, QBLCKB is not
109 * referenced. QBLCKB > 1.
110 *
111 *
112 * Further Details
113 * ===============
114 *
115 * PRTYPE = 1: A and B are Jordan blocks, D and E are identity matrices
116 *
117 * A : if (i == j) then A(i, j) = 1.0
118 * if (j == i + 1) then A(i, j) = -1.0
119 * else A(i, j) = 0.0, i, j = 1...M
120 *
121 * B : if (i == j) then B(i, j) = 1.0 - ALPHA
122 * if (j == i + 1) then B(i, j) = 1.0
123 * else B(i, j) = 0.0, i, j = 1...N
124 *
125 * D : if (i == j) then D(i, j) = 1.0
126 * else D(i, j) = 0.0, i, j = 1...M
127 *
128 * E : if (i == j) then E(i, j) = 1.0
129 * else E(i, j) = 0.0, i, j = 1...N
130 *
131 * L = R are chosen from [-10...10],
132 * which specifies the right hand sides (C, F).
133 *
134 * PRTYPE = 2 or 3: Triangular and/or quasi- triangular.
135 *
136 * A : if (i <= j) then A(i, j) = [-1...1]
137 * else A(i, j) = 0.0, i, j = 1...M
138 *
139 * if (PRTYPE = 3) then
140 * A(k + 1, k + 1) = A(k, k)
141 * A(k + 1, k) = [-1...1]
142 * sign(A(k, k + 1) = -(sin(A(k + 1, k))
143 * k = 1, M - 1, QBLCKA
144 *
145 * B : if (i <= j) then B(i, j) = [-1...1]
146 * else B(i, j) = 0.0, i, j = 1...N
147 *
148 * if (PRTYPE = 3) then
149 * B(k + 1, k + 1) = B(k, k)
150 * B(k + 1, k) = [-1...1]
151 * sign(B(k, k + 1) = -(sign(B(k + 1, k))
152 * k = 1, N - 1, QBLCKB
153 *
154 * D : if (i <= j) then D(i, j) = [-1...1].
155 * else D(i, j) = 0.0, i, j = 1...M
156 *
157 *
158 * E : if (i <= j) then D(i, j) = [-1...1]
159 * else E(i, j) = 0.0, i, j = 1...N
160 *
161 * L, R are chosen from [-10...10],
162 * which specifies the right hand sides (C, F).
163 *
164 * PRTYPE = 4 Full
165 * A(i, j) = [-10...10]
166 * D(i, j) = [-1...1] i,j = 1...M
167 * B(i, j) = [-10...10]
168 * E(i, j) = [-1...1] i,j = 1...N
169 * R(i, j) = [-10...10]
170 * L(i, j) = [-1...1] i = 1..M ,j = 1...N
171 *
172 * L, R specifies the right hand sides (C, F).
173 *
174 * PRTYPE = 5 special case common and/or close eigs.
175 *
176 * =====================================================================
177 *
178 * .. Parameters ..
179 COMPLEX*16 ONE, TWO, ZERO, HALF, TWENTY
180 PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ),
181 $ TWO = ( 2.0D+0, 0.0D+0 ),
182 $ ZERO = ( 0.0D+0, 0.0D+0 ),
183 $ HALF = ( 0.5D+0, 0.0D+0 ),
184 $ TWENTY = ( 2.0D+1, 0.0D+0 ) )
185 * ..
186 * .. Local Scalars ..
187 INTEGER I, J, K
188 COMPLEX*16 IMEPS, REEPS
189 * ..
190 * .. Intrinsic Functions ..
191 INTRINSIC DCMPLX, MOD, SIN
192 * ..
193 * .. External Subroutines ..
194 EXTERNAL ZGEMM
195 * ..
196 * .. Executable Statements ..
197 *
198 IF( PRTYPE.EQ.1 ) THEN
199 DO 20 I = 1, M
200 DO 10 J = 1, M
201 IF( I.EQ.J ) THEN
202 A( I, J ) = ONE
203 D( I, J ) = ONE
204 ELSE IF( I.EQ.J-1 ) THEN
205 A( I, J ) = -ONE
206 D( I, J ) = ZERO
207 ELSE
208 A( I, J ) = ZERO
209 D( I, J ) = ZERO
210 END IF
211 10 CONTINUE
212 20 CONTINUE
213 *
214 DO 40 I = 1, N
215 DO 30 J = 1, N
216 IF( I.EQ.J ) THEN
217 B( I, J ) = ONE - ALPHA
218 E( I, J ) = ONE
219 ELSE IF( I.EQ.J-1 ) THEN
220 B( I, J ) = ONE
221 E( I, J ) = ZERO
222 ELSE
223 B( I, J ) = ZERO
224 E( I, J ) = ZERO
225 END IF
226 30 CONTINUE
227 40 CONTINUE
228 *
229 DO 60 I = 1, M
230 DO 50 J = 1, N
231 R( I, J ) = ( HALF-SIN( DCMPLX( I / J ) ) )*TWENTY
232 L( I, J ) = R( I, J )
233 50 CONTINUE
234 60 CONTINUE
235 *
236 ELSE IF( PRTYPE.EQ.2 .OR. PRTYPE.EQ.3 ) THEN
237 DO 80 I = 1, M
238 DO 70 J = 1, M
239 IF( I.LE.J ) THEN
240 A( I, J ) = ( HALF-SIN( DCMPLX( I ) ) )*TWO
241 D( I, J ) = ( HALF-SIN( DCMPLX( I*J ) ) )*TWO
242 ELSE
243 A( I, J ) = ZERO
244 D( I, J ) = ZERO
245 END IF
246 70 CONTINUE
247 80 CONTINUE
248 *
249 DO 100 I = 1, N
250 DO 90 J = 1, N
251 IF( I.LE.J ) THEN
252 B( I, J ) = ( HALF-SIN( DCMPLX( I+J ) ) )*TWO
253 E( I, J ) = ( HALF-SIN( DCMPLX( J ) ) )*TWO
254 ELSE
255 B( I, J ) = ZERO
256 E( I, J ) = ZERO
257 END IF
258 90 CONTINUE
259 100 CONTINUE
260 *
261 DO 120 I = 1, M
262 DO 110 J = 1, N
263 R( I, J ) = ( HALF-SIN( DCMPLX( I*J ) ) )*TWENTY
264 L( I, J ) = ( HALF-SIN( DCMPLX( I+J ) ) )*TWENTY
265 110 CONTINUE
266 120 CONTINUE
267 *
268 IF( PRTYPE.EQ.3 ) THEN
269 IF( QBLCKA.LE.1 )
270 $ QBLCKA = 2
271 DO 130 K = 1, M - 1, QBLCKA
272 A( K+1, K+1 ) = A( K, K )
273 A( K+1, K ) = -SIN( A( K, K+1 ) )
274 130 CONTINUE
275 *
276 IF( QBLCKB.LE.1 )
277 $ QBLCKB = 2
278 DO 140 K = 1, N - 1, QBLCKB
279 B( K+1, K+1 ) = B( K, K )
280 B( K+1, K ) = -SIN( B( K, K+1 ) )
281 140 CONTINUE
282 END IF
283 *
284 ELSE IF( PRTYPE.EQ.4 ) THEN
285 DO 160 I = 1, M
286 DO 150 J = 1, M
287 A( I, J ) = ( HALF-SIN( DCMPLX( I*J ) ) )*TWENTY
288 D( I, J ) = ( HALF-SIN( DCMPLX( I+J ) ) )*TWO
289 150 CONTINUE
290 160 CONTINUE
291 *
292 DO 180 I = 1, N
293 DO 170 J = 1, N
294 B( I, J ) = ( HALF-SIN( DCMPLX( I+J ) ) )*TWENTY
295 E( I, J ) = ( HALF-SIN( DCMPLX( I*J ) ) )*TWO
296 170 CONTINUE
297 180 CONTINUE
298 *
299 DO 200 I = 1, M
300 DO 190 J = 1, N
301 R( I, J ) = ( HALF-SIN( DCMPLX( J / I ) ) )*TWENTY
302 L( I, J ) = ( HALF-SIN( DCMPLX( I*J ) ) )*TWO
303 190 CONTINUE
304 200 CONTINUE
305 *
306 ELSE IF( PRTYPE.GE.5 ) THEN
307 REEPS = HALF*TWO*TWENTY / ALPHA
308 IMEPS = ( HALF-TWO ) / ALPHA
309 DO 220 I = 1, M
310 DO 210 J = 1, N
311 R( I, J ) = ( HALF-SIN( DCMPLX( I*J ) ) )*ALPHA / TWENTY
312 L( I, J ) = ( HALF-SIN( DCMPLX( I+J ) ) )*ALPHA / TWENTY
313 210 CONTINUE
314 220 CONTINUE
315 *
316 DO 230 I = 1, M
317 D( I, I ) = ONE
318 230 CONTINUE
319 *
320 DO 240 I = 1, M
321 IF( I.LE.4 ) THEN
322 A( I, I ) = ONE
323 IF( I.GT.2 )
324 $ A( I, I ) = ONE + REEPS
325 IF( MOD( I, 2 ).NE.0 .AND. I.LT.M ) THEN
326 A( I, I+1 ) = IMEPS
327 ELSE IF( I.GT.1 ) THEN
328 A( I, I-1 ) = -IMEPS
329 END IF
330 ELSE IF( I.LE.8 ) THEN
331 IF( I.LE.6 ) THEN
332 A( I, I ) = REEPS
333 ELSE
334 A( I, I ) = -REEPS
335 END IF
336 IF( MOD( I, 2 ).NE.0 .AND. I.LT.M ) THEN
337 A( I, I+1 ) = ONE
338 ELSE IF( I.GT.1 ) THEN
339 A( I, I-1 ) = -ONE
340 END IF
341 ELSE
342 A( I, I ) = ONE
343 IF( MOD( I, 2 ).NE.0 .AND. I.LT.M ) THEN
344 A( I, I+1 ) = IMEPS*2
345 ELSE IF( I.GT.1 ) THEN
346 A( I, I-1 ) = -IMEPS*2
347 END IF
348 END IF
349 240 CONTINUE
350 *
351 DO 250 I = 1, N
352 E( I, I ) = ONE
353 IF( I.LE.4 ) THEN
354 B( I, I ) = -ONE
355 IF( I.GT.2 )
356 $ B( I, I ) = ONE - REEPS
357 IF( MOD( I, 2 ).NE.0 .AND. I.LT.N ) THEN
358 B( I, I+1 ) = IMEPS
359 ELSE IF( I.GT.1 ) THEN
360 B( I, I-1 ) = -IMEPS
361 END IF
362 ELSE IF( I.LE.8 ) THEN
363 IF( I.LE.6 ) THEN
364 B( I, I ) = REEPS
365 ELSE
366 B( I, I ) = -REEPS
367 END IF
368 IF( MOD( I, 2 ).NE.0 .AND. I.LT.N ) THEN
369 B( I, I+1 ) = ONE + IMEPS
370 ELSE IF( I.GT.1 ) THEN
371 B( I, I-1 ) = -ONE - IMEPS
372 END IF
373 ELSE
374 B( I, I ) = ONE - REEPS
375 IF( MOD( I, 2 ).NE.0 .AND. I.LT.N ) THEN
376 B( I, I+1 ) = IMEPS*2
377 ELSE IF( I.GT.1 ) THEN
378 B( I, I-1 ) = -IMEPS*2
379 END IF
380 END IF
381 250 CONTINUE
382 END IF
383 *
384 * Compute rhs (C, F)
385 *
386 CALL ZGEMM( 'N', 'N', M, N, M, ONE, A, LDA, R, LDR, ZERO, C, LDC )
387 CALL ZGEMM( 'N', 'N', M, N, N, -ONE, L, LDL, B, LDB, ONE, C, LDC )
388 CALL ZGEMM( 'N', 'N', M, N, M, ONE, D, LDD, R, LDR, ZERO, F, LDF )
389 CALL ZGEMM( 'N', 'N', M, N, N, -ONE, L, LDL, E, LDE, ONE, F, LDF )
390 *
391 * End of ZLATM5
392 *
393 END
2 $ E, LDE, F, LDF, R, LDR, L, LDL, ALPHA, QBLCKA,
3 $ QBLCKB )
4 *
5 * -- LAPACK test routine (version 3.1) --
6 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
7 * November 2006
8 *
9 * .. Scalar Arguments ..
10 INTEGER LDA, LDB, LDC, LDD, LDE, LDF, LDL, LDR, M, N,
11 $ PRTYPE, QBLCKA, QBLCKB
12 DOUBLE PRECISION ALPHA
13 * ..
14 * .. Array Arguments ..
15 COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * ),
16 $ D( LDD, * ), E( LDE, * ), F( LDF, * ),
17 $ L( LDL, * ), R( LDR, * )
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 * ZLATM5 generates matrices involved in the Generalized Sylvester
24 * equation:
25 *
26 * A * R - L * B = C
27 * D * R - L * E = F
28 *
29 * They also satisfy (the diagonalization condition)
30 *
31 * [ I -L ] ( [ A -C ], [ D -F ] ) [ I R ] = ( [ A ], [ D ] )
32 * [ I ] ( [ B ] [ E ] ) [ I ] ( [ B ] [ E ] )
33 *
34 *
35 * Arguments
36 * =========
37 *
38 * PRTYPE (input) INTEGER
39 * "Points" to a certian type of the matrices to generate
40 * (see futher details).
41 *
42 * M (input) INTEGER
43 * Specifies the order of A and D and the number of rows in
44 * C, F, R and L.
45 *
46 * N (input) INTEGER
47 * Specifies the order of B and E and the number of columns in
48 * C, F, R and L.
49 *
50 * A (output) COMPLEX*16 array, dimension (LDA, M).
51 * On exit A M-by-M is initialized according to PRTYPE.
52 *
53 * LDA (input) INTEGER
54 * The leading dimension of A.
55 *
56 * B (output) COMPLEX*16 array, dimension (LDB, N).
57 * On exit B N-by-N is initialized according to PRTYPE.
58 *
59 * LDB (input) INTEGER
60 * The leading dimension of B.
61 *
62 * C (output) COMPLEX*16 array, dimension (LDC, N).
63 * On exit C M-by-N is initialized according to PRTYPE.
64 *
65 * LDC (input) INTEGER
66 * The leading dimension of C.
67 *
68 * D (output) COMPLEX*16 array, dimension (LDD, M).
69 * On exit D M-by-M is initialized according to PRTYPE.
70 *
71 * LDD (input) INTEGER
72 * The leading dimension of D.
73 *
74 * E (output) COMPLEX*16 array, dimension (LDE, N).
75 * On exit E N-by-N is initialized according to PRTYPE.
76 *
77 * LDE (input) INTEGER
78 * The leading dimension of E.
79 *
80 * F (output) COMPLEX*16 array, dimension (LDF, N).
81 * On exit F M-by-N is initialized according to PRTYPE.
82 *
83 * LDF (input) INTEGER
84 * The leading dimension of F.
85 *
86 * R (output) COMPLEX*16 array, dimension (LDR, N).
87 * On exit R M-by-N is initialized according to PRTYPE.
88 *
89 * LDR (input) INTEGER
90 * The leading dimension of R.
91 *
92 * L (output) COMPLEX*16 array, dimension (LDL, N).
93 * On exit L M-by-N is initialized according to PRTYPE.
94 *
95 * LDL (input) INTEGER
96 * The leading dimension of L.
97 *
98 * ALPHA (input) DOUBLE PRECISION
99 * Parameter used in generating PRTYPE = 1 and 5 matrices.
100 *
101 * QBLCKA (input) INTEGER
102 * When PRTYPE = 3, specifies the distance between 2-by-2
103 * blocks on the diagonal in A. Otherwise, QBLCKA is not
104 * referenced. QBLCKA > 1.
105 *
106 * QBLCKB (input) INTEGER
107 * When PRTYPE = 3, specifies the distance between 2-by-2
108 * blocks on the diagonal in B. Otherwise, QBLCKB is not
109 * referenced. QBLCKB > 1.
110 *
111 *
112 * Further Details
113 * ===============
114 *
115 * PRTYPE = 1: A and B are Jordan blocks, D and E are identity matrices
116 *
117 * A : if (i == j) then A(i, j) = 1.0
118 * if (j == i + 1) then A(i, j) = -1.0
119 * else A(i, j) = 0.0, i, j = 1...M
120 *
121 * B : if (i == j) then B(i, j) = 1.0 - ALPHA
122 * if (j == i + 1) then B(i, j) = 1.0
123 * else B(i, j) = 0.0, i, j = 1...N
124 *
125 * D : if (i == j) then D(i, j) = 1.0
126 * else D(i, j) = 0.0, i, j = 1...M
127 *
128 * E : if (i == j) then E(i, j) = 1.0
129 * else E(i, j) = 0.0, i, j = 1...N
130 *
131 * L = R are chosen from [-10...10],
132 * which specifies the right hand sides (C, F).
133 *
134 * PRTYPE = 2 or 3: Triangular and/or quasi- triangular.
135 *
136 * A : if (i <= j) then A(i, j) = [-1...1]
137 * else A(i, j) = 0.0, i, j = 1...M
138 *
139 * if (PRTYPE = 3) then
140 * A(k + 1, k + 1) = A(k, k)
141 * A(k + 1, k) = [-1...1]
142 * sign(A(k, k + 1) = -(sin(A(k + 1, k))
143 * k = 1, M - 1, QBLCKA
144 *
145 * B : if (i <= j) then B(i, j) = [-1...1]
146 * else B(i, j) = 0.0, i, j = 1...N
147 *
148 * if (PRTYPE = 3) then
149 * B(k + 1, k + 1) = B(k, k)
150 * B(k + 1, k) = [-1...1]
151 * sign(B(k, k + 1) = -(sign(B(k + 1, k))
152 * k = 1, N - 1, QBLCKB
153 *
154 * D : if (i <= j) then D(i, j) = [-1...1].
155 * else D(i, j) = 0.0, i, j = 1...M
156 *
157 *
158 * E : if (i <= j) then D(i, j) = [-1...1]
159 * else E(i, j) = 0.0, i, j = 1...N
160 *
161 * L, R are chosen from [-10...10],
162 * which specifies the right hand sides (C, F).
163 *
164 * PRTYPE = 4 Full
165 * A(i, j) = [-10...10]
166 * D(i, j) = [-1...1] i,j = 1...M
167 * B(i, j) = [-10...10]
168 * E(i, j) = [-1...1] i,j = 1...N
169 * R(i, j) = [-10...10]
170 * L(i, j) = [-1...1] i = 1..M ,j = 1...N
171 *
172 * L, R specifies the right hand sides (C, F).
173 *
174 * PRTYPE = 5 special case common and/or close eigs.
175 *
176 * =====================================================================
177 *
178 * .. Parameters ..
179 COMPLEX*16 ONE, TWO, ZERO, HALF, TWENTY
180 PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ),
181 $ TWO = ( 2.0D+0, 0.0D+0 ),
182 $ ZERO = ( 0.0D+0, 0.0D+0 ),
183 $ HALF = ( 0.5D+0, 0.0D+0 ),
184 $ TWENTY = ( 2.0D+1, 0.0D+0 ) )
185 * ..
186 * .. Local Scalars ..
187 INTEGER I, J, K
188 COMPLEX*16 IMEPS, REEPS
189 * ..
190 * .. Intrinsic Functions ..
191 INTRINSIC DCMPLX, MOD, SIN
192 * ..
193 * .. External Subroutines ..
194 EXTERNAL ZGEMM
195 * ..
196 * .. Executable Statements ..
197 *
198 IF( PRTYPE.EQ.1 ) THEN
199 DO 20 I = 1, M
200 DO 10 J = 1, M
201 IF( I.EQ.J ) THEN
202 A( I, J ) = ONE
203 D( I, J ) = ONE
204 ELSE IF( I.EQ.J-1 ) THEN
205 A( I, J ) = -ONE
206 D( I, J ) = ZERO
207 ELSE
208 A( I, J ) = ZERO
209 D( I, J ) = ZERO
210 END IF
211 10 CONTINUE
212 20 CONTINUE
213 *
214 DO 40 I = 1, N
215 DO 30 J = 1, N
216 IF( I.EQ.J ) THEN
217 B( I, J ) = ONE - ALPHA
218 E( I, J ) = ONE
219 ELSE IF( I.EQ.J-1 ) THEN
220 B( I, J ) = ONE
221 E( I, J ) = ZERO
222 ELSE
223 B( I, J ) = ZERO
224 E( I, J ) = ZERO
225 END IF
226 30 CONTINUE
227 40 CONTINUE
228 *
229 DO 60 I = 1, M
230 DO 50 J = 1, N
231 R( I, J ) = ( HALF-SIN( DCMPLX( I / J ) ) )*TWENTY
232 L( I, J ) = R( I, J )
233 50 CONTINUE
234 60 CONTINUE
235 *
236 ELSE IF( PRTYPE.EQ.2 .OR. PRTYPE.EQ.3 ) THEN
237 DO 80 I = 1, M
238 DO 70 J = 1, M
239 IF( I.LE.J ) THEN
240 A( I, J ) = ( HALF-SIN( DCMPLX( I ) ) )*TWO
241 D( I, J ) = ( HALF-SIN( DCMPLX( I*J ) ) )*TWO
242 ELSE
243 A( I, J ) = ZERO
244 D( I, J ) = ZERO
245 END IF
246 70 CONTINUE
247 80 CONTINUE
248 *
249 DO 100 I = 1, N
250 DO 90 J = 1, N
251 IF( I.LE.J ) THEN
252 B( I, J ) = ( HALF-SIN( DCMPLX( I+J ) ) )*TWO
253 E( I, J ) = ( HALF-SIN( DCMPLX( J ) ) )*TWO
254 ELSE
255 B( I, J ) = ZERO
256 E( I, J ) = ZERO
257 END IF
258 90 CONTINUE
259 100 CONTINUE
260 *
261 DO 120 I = 1, M
262 DO 110 J = 1, N
263 R( I, J ) = ( HALF-SIN( DCMPLX( I*J ) ) )*TWENTY
264 L( I, J ) = ( HALF-SIN( DCMPLX( I+J ) ) )*TWENTY
265 110 CONTINUE
266 120 CONTINUE
267 *
268 IF( PRTYPE.EQ.3 ) THEN
269 IF( QBLCKA.LE.1 )
270 $ QBLCKA = 2
271 DO 130 K = 1, M - 1, QBLCKA
272 A( K+1, K+1 ) = A( K, K )
273 A( K+1, K ) = -SIN( A( K, K+1 ) )
274 130 CONTINUE
275 *
276 IF( QBLCKB.LE.1 )
277 $ QBLCKB = 2
278 DO 140 K = 1, N - 1, QBLCKB
279 B( K+1, K+1 ) = B( K, K )
280 B( K+1, K ) = -SIN( B( K, K+1 ) )
281 140 CONTINUE
282 END IF
283 *
284 ELSE IF( PRTYPE.EQ.4 ) THEN
285 DO 160 I = 1, M
286 DO 150 J = 1, M
287 A( I, J ) = ( HALF-SIN( DCMPLX( I*J ) ) )*TWENTY
288 D( I, J ) = ( HALF-SIN( DCMPLX( I+J ) ) )*TWO
289 150 CONTINUE
290 160 CONTINUE
291 *
292 DO 180 I = 1, N
293 DO 170 J = 1, N
294 B( I, J ) = ( HALF-SIN( DCMPLX( I+J ) ) )*TWENTY
295 E( I, J ) = ( HALF-SIN( DCMPLX( I*J ) ) )*TWO
296 170 CONTINUE
297 180 CONTINUE
298 *
299 DO 200 I = 1, M
300 DO 190 J = 1, N
301 R( I, J ) = ( HALF-SIN( DCMPLX( J / I ) ) )*TWENTY
302 L( I, J ) = ( HALF-SIN( DCMPLX( I*J ) ) )*TWO
303 190 CONTINUE
304 200 CONTINUE
305 *
306 ELSE IF( PRTYPE.GE.5 ) THEN
307 REEPS = HALF*TWO*TWENTY / ALPHA
308 IMEPS = ( HALF-TWO ) / ALPHA
309 DO 220 I = 1, M
310 DO 210 J = 1, N
311 R( I, J ) = ( HALF-SIN( DCMPLX( I*J ) ) )*ALPHA / TWENTY
312 L( I, J ) = ( HALF-SIN( DCMPLX( I+J ) ) )*ALPHA / TWENTY
313 210 CONTINUE
314 220 CONTINUE
315 *
316 DO 230 I = 1, M
317 D( I, I ) = ONE
318 230 CONTINUE
319 *
320 DO 240 I = 1, M
321 IF( I.LE.4 ) THEN
322 A( I, I ) = ONE
323 IF( I.GT.2 )
324 $ A( I, I ) = ONE + REEPS
325 IF( MOD( I, 2 ).NE.0 .AND. I.LT.M ) THEN
326 A( I, I+1 ) = IMEPS
327 ELSE IF( I.GT.1 ) THEN
328 A( I, I-1 ) = -IMEPS
329 END IF
330 ELSE IF( I.LE.8 ) THEN
331 IF( I.LE.6 ) THEN
332 A( I, I ) = REEPS
333 ELSE
334 A( I, I ) = -REEPS
335 END IF
336 IF( MOD( I, 2 ).NE.0 .AND. I.LT.M ) THEN
337 A( I, I+1 ) = ONE
338 ELSE IF( I.GT.1 ) THEN
339 A( I, I-1 ) = -ONE
340 END IF
341 ELSE
342 A( I, I ) = ONE
343 IF( MOD( I, 2 ).NE.0 .AND. I.LT.M ) THEN
344 A( I, I+1 ) = IMEPS*2
345 ELSE IF( I.GT.1 ) THEN
346 A( I, I-1 ) = -IMEPS*2
347 END IF
348 END IF
349 240 CONTINUE
350 *
351 DO 250 I = 1, N
352 E( I, I ) = ONE
353 IF( I.LE.4 ) THEN
354 B( I, I ) = -ONE
355 IF( I.GT.2 )
356 $ B( I, I ) = ONE - REEPS
357 IF( MOD( I, 2 ).NE.0 .AND. I.LT.N ) THEN
358 B( I, I+1 ) = IMEPS
359 ELSE IF( I.GT.1 ) THEN
360 B( I, I-1 ) = -IMEPS
361 END IF
362 ELSE IF( I.LE.8 ) THEN
363 IF( I.LE.6 ) THEN
364 B( I, I ) = REEPS
365 ELSE
366 B( I, I ) = -REEPS
367 END IF
368 IF( MOD( I, 2 ).NE.0 .AND. I.LT.N ) THEN
369 B( I, I+1 ) = ONE + IMEPS
370 ELSE IF( I.GT.1 ) THEN
371 B( I, I-1 ) = -ONE - IMEPS
372 END IF
373 ELSE
374 B( I, I ) = ONE - REEPS
375 IF( MOD( I, 2 ).NE.0 .AND. I.LT.N ) THEN
376 B( I, I+1 ) = IMEPS*2
377 ELSE IF( I.GT.1 ) THEN
378 B( I, I-1 ) = -IMEPS*2
379 END IF
380 END IF
381 250 CONTINUE
382 END IF
383 *
384 * Compute rhs (C, F)
385 *
386 CALL ZGEMM( 'N', 'N', M, N, M, ONE, A, LDA, R, LDR, ZERO, C, LDC )
387 CALL ZGEMM( 'N', 'N', M, N, N, -ONE, L, LDL, B, LDB, ONE, C, LDC )
388 CALL ZGEMM( 'N', 'N', M, N, M, ONE, D, LDD, R, LDR, ZERO, F, LDF )
389 CALL ZGEMM( 'N', 'N', M, N, N, -ONE, L, LDL, E, LDE, ONE, F, LDF )
390 *
391 * End of ZLATM5
392 *
393 END