LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
clatm5.f
Go to the documentation of this file.
1*> \brief \b CLATM5
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 CLATM5( PRTYPE, M, N, A, LDA, B, LDB, C, LDC, D, LDD,
12* E, LDE, F, LDF, R, LDR, L, LDL, ALPHA, QBLCKA,
13* QBLCKB )
14*
15* .. Scalar Arguments ..
16* INTEGER LDA, LDB, LDC, LDD, LDE, LDF, LDL, LDR, M, N,
17* $ PRTYPE, QBLCKA, QBLCKB
18* REAL ALPHA
19* ..
20* .. Array Arguments ..
21* COMPLEX A( LDA, * ), B( LDB, * ), C( LDC, * ),
22* $ D( LDD, * ), E( LDE, * ), F( LDF, * ),
23* $ L( LDL, * ), R( LDR, * )
24* ..
25*
26*
27*> \par Purpose:
28* =============
29*>
30*> \verbatim
31*>
32*> CLATM5 generates matrices involved in the Generalized Sylvester
33*> equation:
34*>
35*> A * R - L * B = C
36*> D * R - L * E = F
37*>
38*> They also satisfy (the diagonalization condition)
39*>
40*> [ I -L ] ( [ A -C ], [ D -F ] ) [ I R ] = ( [ A ], [ D ] )
41*> [ I ] ( [ B ] [ E ] ) [ I ] ( [ B ] [ E ] )
42*>
43*> \endverbatim
44*
45* Arguments:
46* ==========
47*
48*> \param[in] PRTYPE
49*> \verbatim
50*> PRTYPE is INTEGER
51*> "Points" to a certain type of the matrices to generate
52*> (see further details).
53*> \endverbatim
54*>
55*> \param[in] M
56*> \verbatim
57*> M is INTEGER
58*> Specifies the order of A and D and the number of rows in
59*> C, F, R and L.
60*> \endverbatim
61*>
62*> \param[in] N
63*> \verbatim
64*> N is INTEGER
65*> Specifies the order of B and E and the number of columns in
66*> C, F, R and L.
67*> \endverbatim
68*>
69*> \param[out] A
70*> \verbatim
71*> A is COMPLEX array, dimension (LDA, M).
72*> On exit A M-by-M is initialized according to PRTYPE.
73*> \endverbatim
74*>
75*> \param[in] LDA
76*> \verbatim
77*> LDA is INTEGER
78*> The leading dimension of A.
79*> \endverbatim
80*>
81*> \param[out] B
82*> \verbatim
83*> B is COMPLEX array, dimension (LDB, N).
84*> On exit B N-by-N is initialized according to PRTYPE.
85*> \endverbatim
86*>
87*> \param[in] LDB
88*> \verbatim
89*> LDB is INTEGER
90*> The leading dimension of B.
91*> \endverbatim
92*>
93*> \param[out] C
94*> \verbatim
95*> C is COMPLEX array, dimension (LDC, N).
96*> On exit C M-by-N is initialized according to PRTYPE.
97*> \endverbatim
98*>
99*> \param[in] LDC
100*> \verbatim
101*> LDC is INTEGER
102*> The leading dimension of C.
103*> \endverbatim
104*>
105*> \param[out] D
106*> \verbatim
107*> D is COMPLEX array, dimension (LDD, M).
108*> On exit D M-by-M is initialized according to PRTYPE.
109*> \endverbatim
110*>
111*> \param[in] LDD
112*> \verbatim
113*> LDD is INTEGER
114*> The leading dimension of D.
115*> \endverbatim
116*>
117*> \param[out] E
118*> \verbatim
119*> E is COMPLEX array, dimension (LDE, N).
120*> On exit E N-by-N is initialized according to PRTYPE.
121*> \endverbatim
122*>
123*> \param[in] LDE
124*> \verbatim
125*> LDE is INTEGER
126*> The leading dimension of E.
127*> \endverbatim
128*>
129*> \param[out] F
130*> \verbatim
131*> F is COMPLEX array, dimension (LDF, N).
132*> On exit F M-by-N is initialized according to PRTYPE.
133*> \endverbatim
134*>
135*> \param[in] LDF
136*> \verbatim
137*> LDF is INTEGER
138*> The leading dimension of F.
139*> \endverbatim
140*>
141*> \param[out] R
142*> \verbatim
143*> R is COMPLEX array, dimension (LDR, N).
144*> On exit R M-by-N is initialized according to PRTYPE.
145*> \endverbatim
146*>
147*> \param[in] LDR
148*> \verbatim
149*> LDR is INTEGER
150*> The leading dimension of R.
151*> \endverbatim
152*>
153*> \param[out] L
154*> \verbatim
155*> L is COMPLEX array, dimension (LDL, N).
156*> On exit L M-by-N is initialized according to PRTYPE.
157*> \endverbatim
158*>
159*> \param[in] LDL
160*> \verbatim
161*> LDL is INTEGER
162*> The leading dimension of L.
163*> \endverbatim
164*>
165*> \param[in] ALPHA
166*> \verbatim
167*> ALPHA is REAL
168*> Parameter used in generating PRTYPE = 1 and 5 matrices.
169*> \endverbatim
170*>
171*> \param[in] QBLCKA
172*> \verbatim
173*> QBLCKA is INTEGER
174*> When PRTYPE = 3, specifies the distance between 2-by-2
175*> blocks on the diagonal in A. Otherwise, QBLCKA is not
176*> referenced. QBLCKA > 1.
177*> \endverbatim
178*>
179*> \param[in] QBLCKB
180*> \verbatim
181*> QBLCKB is INTEGER
182*> When PRTYPE = 3, specifies the distance between 2-by-2
183*> blocks on the diagonal in B. Otherwise, QBLCKB is not
184*> referenced. QBLCKB > 1.
185*> \endverbatim
186*
187* Authors:
188* ========
189*
190*> \author Univ. of Tennessee
191*> \author Univ. of California Berkeley
192*> \author Univ. of Colorado Denver
193*> \author NAG Ltd.
194*
195*> \ingroup complex_matgen
196*
197*> \par Further Details:
198* =====================
199*>
200*> \verbatim
201*>
202*> PRTYPE = 1: A and B are Jordan blocks, D and E are identity matrices
203*>
204*> A : if (i == j) then A(i, j) = 1.0
205*> if (j == i + 1) then A(i, j) = -1.0
206*> else A(i, j) = 0.0, i, j = 1...M
207*>
208*> B : if (i == j) then B(i, j) = 1.0 - ALPHA
209*> if (j == i + 1) then B(i, j) = 1.0
210*> else B(i, j) = 0.0, i, j = 1...N
211*>
212*> D : if (i == j) then D(i, j) = 1.0
213*> else D(i, j) = 0.0, i, j = 1...M
214*>
215*> E : if (i == j) then E(i, j) = 1.0
216*> else E(i, j) = 0.0, i, j = 1...N
217*>
218*> L = R are chosen from [-10...10],
219*> which specifies the right hand sides (C, F).
220*>
221*> PRTYPE = 2 or 3: Triangular and/or quasi- triangular.
222*>
223*> A : if (i <= j) then A(i, j) = [-1...1]
224*> else A(i, j) = 0.0, i, j = 1...M
225*>
226*> if (PRTYPE = 3) then
227*> A(k + 1, k + 1) = A(k, k)
228*> A(k + 1, k) = [-1...1]
229*> sign(A(k, k + 1) = -(sin(A(k + 1, k))
230*> k = 1, M - 1, QBLCKA
231*>
232*> B : if (i <= j) then B(i, j) = [-1...1]
233*> else B(i, j) = 0.0, i, j = 1...N
234*>
235*> if (PRTYPE = 3) then
236*> B(k + 1, k + 1) = B(k, k)
237*> B(k + 1, k) = [-1...1]
238*> sign(B(k, k + 1) = -(sign(B(k + 1, k))
239*> k = 1, N - 1, QBLCKB
240*>
241*> D : if (i <= j) then D(i, j) = [-1...1].
242*> else D(i, j) = 0.0, i, j = 1...M
243*>
244*>
245*> E : if (i <= j) then D(i, j) = [-1...1]
246*> else E(i, j) = 0.0, i, j = 1...N
247*>
248*> L, R are chosen from [-10...10],
249*> which specifies the right hand sides (C, F).
250*>
251*> PRTYPE = 4 Full
252*> A(i, j) = [-10...10]
253*> D(i, j) = [-1...1] i,j = 1...M
254*> B(i, j) = [-10...10]
255*> E(i, j) = [-1...1] i,j = 1...N
256*> R(i, j) = [-10...10]
257*> L(i, j) = [-1...1] i = 1..M ,j = 1...N
258*>
259*> L, R specifies the right hand sides (C, F).
260*>
261*> PRTYPE = 5 special case common and/or close eigs.
262*> \endverbatim
263*>
264* =====================================================================
265 SUBROUTINE clatm5( PRTYPE, M, N, A, LDA, B, LDB, C, LDC, D,
266 $ LDD,
267 $ E, LDE, F, LDF, R, LDR, L, LDL, ALPHA, QBLCKA,
268 $ QBLCKB )
269*
270* -- LAPACK computational routine --
271* -- LAPACK is a software package provided by Univ. of Tennessee, --
272* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
273*
274* .. Scalar Arguments ..
275 INTEGER LDA, LDB, LDC, LDD, LDE, LDF, LDL, LDR, M, N,
276 $ PRTYPE, QBLCKA, QBLCKB
277 REAL ALPHA
278* ..
279* .. Array Arguments ..
280 COMPLEX A( LDA, * ), B( LDB, * ), C( LDC, * ),
281 $ D( LDD, * ), E( LDE, * ), F( LDF, * ),
282 $ L( LDL, * ), R( LDR, * )
283* ..
284*
285* =====================================================================
286*
287* .. Parameters ..
288 COMPLEX ONE, TWO, ZERO, HALF, TWENTY
289 PARAMETER ( ONE = ( 1.0e+0, 0.0e+0 ),
290 $ two = ( 2.0e+0, 0.0e+0 ),
291 $ zero = ( 0.0e+0, 0.0e+0 ),
292 $ half = ( 0.5e+0, 0.0e+0 ),
293 $ twenty = ( 2.0e+1, 0.0e+0 ) )
294* ..
295* .. Local Scalars ..
296 INTEGER I, J, K
297 COMPLEX IMEPS, REEPS
298* ..
299* .. Intrinsic Functions ..
300 INTRINSIC cmplx, mod, sin
301* ..
302* .. External Subroutines ..
303 EXTERNAL cgemm
304* ..
305* .. Executable Statements ..
306*
307 IF( prtype.EQ.1 ) THEN
308 DO 20 i = 1, m
309 DO 10 j = 1, m
310 IF( i.EQ.j ) THEN
311 a( i, j ) = one
312 d( i, j ) = one
313 ELSE IF( i.EQ.j-1 ) THEN
314 a( i, j ) = -one
315 d( i, j ) = zero
316 ELSE
317 a( i, j ) = zero
318 d( i, j ) = zero
319 END IF
320 10 CONTINUE
321 20 CONTINUE
322*
323 DO 40 i = 1, n
324 DO 30 j = 1, n
325 IF( i.EQ.j ) THEN
326 b( i, j ) = one - alpha
327 e( i, j ) = one
328 ELSE IF( i.EQ.j-1 ) THEN
329 b( i, j ) = one
330 e( i, j ) = zero
331 ELSE
332 b( i, j ) = zero
333 e( i, j ) = zero
334 END IF
335 30 CONTINUE
336 40 CONTINUE
337*
338 DO 60 i = 1, m
339 DO 50 j = 1, n
340 r( i, j ) = ( half-sin( cmplx( i / j ) ) )*twenty
341 l( i, j ) = r( i, j )
342 50 CONTINUE
343 60 CONTINUE
344*
345 ELSE IF( prtype.EQ.2 .OR. prtype.EQ.3 ) THEN
346 DO 80 i = 1, m
347 DO 70 j = 1, m
348 IF( i.LE.j ) THEN
349 a( i, j ) = ( half-sin( cmplx( i ) ) )*two
350 d( i, j ) = ( half-sin( cmplx( i*j ) ) )*two
351 ELSE
352 a( i, j ) = zero
353 d( i, j ) = zero
354 END IF
355 70 CONTINUE
356 80 CONTINUE
357*
358 DO 100 i = 1, n
359 DO 90 j = 1, n
360 IF( i.LE.j ) THEN
361 b( i, j ) = ( half-sin( cmplx( i+j ) ) )*two
362 e( i, j ) = ( half-sin( cmplx( j ) ) )*two
363 ELSE
364 b( i, j ) = zero
365 e( i, j ) = zero
366 END IF
367 90 CONTINUE
368 100 CONTINUE
369*
370 DO 120 i = 1, m
371 DO 110 j = 1, n
372 r( i, j ) = ( half-sin( cmplx( i*j ) ) )*twenty
373 l( i, j ) = ( half-sin( cmplx( i+j ) ) )*twenty
374 110 CONTINUE
375 120 CONTINUE
376*
377 IF( prtype.EQ.3 ) THEN
378 IF( qblcka.LE.1 )
379 $ qblcka = 2
380 DO 130 k = 1, m - 1, qblcka
381 a( k+1, k+1 ) = a( k, k )
382 a( k+1, k ) = -sin( a( k, k+1 ) )
383 130 CONTINUE
384*
385 IF( qblckb.LE.1 )
386 $ qblckb = 2
387 DO 140 k = 1, n - 1, qblckb
388 b( k+1, k+1 ) = b( k, k )
389 b( k+1, k ) = -sin( b( k, k+1 ) )
390 140 CONTINUE
391 END IF
392*
393 ELSE IF( prtype.EQ.4 ) THEN
394 DO 160 i = 1, m
395 DO 150 j = 1, m
396 a( i, j ) = ( half-sin( cmplx( i*j ) ) )*twenty
397 d( i, j ) = ( half-sin( cmplx( i+j ) ) )*two
398 150 CONTINUE
399 160 CONTINUE
400*
401 DO 180 i = 1, n
402 DO 170 j = 1, n
403 b( i, j ) = ( half-sin( cmplx( i+j ) ) )*twenty
404 e( i, j ) = ( half-sin( cmplx( i*j ) ) )*two
405 170 CONTINUE
406 180 CONTINUE
407*
408 DO 200 i = 1, m
409 DO 190 j = 1, n
410 r( i, j ) = ( half-sin( cmplx( j / i ) ) )*twenty
411 l( i, j ) = ( half-sin( cmplx( i*j ) ) )*two
412 190 CONTINUE
413 200 CONTINUE
414*
415 ELSE IF( prtype.GE.5 ) THEN
416 reeps = half*two*twenty / alpha
417 imeps = ( half-two ) / alpha
418 DO 220 i = 1, m
419 DO 210 j = 1, n
420 r( i, j ) = ( half-sin( cmplx( i*j ) ) )*alpha / twenty
421 l( i, j ) = ( half-sin( cmplx( i+j ) ) )*alpha / twenty
422 210 CONTINUE
423 220 CONTINUE
424*
425 DO 230 i = 1, m
426 d( i, i ) = one
427 230 CONTINUE
428*
429 DO 240 i = 1, m
430 IF( i.LE.4 ) THEN
431 a( i, i ) = one
432 IF( i.GT.2 )
433 $ a( i, i ) = one + reeps
434 IF( mod( i, 2 ).NE.0 .AND. i.LT.m ) THEN
435 a( i, i+1 ) = imeps
436 ELSE IF( i.GT.1 ) THEN
437 a( i, i-1 ) = -imeps
438 END IF
439 ELSE IF( i.LE.8 ) THEN
440 IF( i.LE.6 ) THEN
441 a( i, i ) = reeps
442 ELSE
443 a( i, i ) = -reeps
444 END IF
445 IF( mod( i, 2 ).NE.0 .AND. i.LT.m ) THEN
446 a( i, i+1 ) = one
447 ELSE IF( i.GT.1 ) THEN
448 a( i, i-1 ) = -one
449 END IF
450 ELSE
451 a( i, i ) = one
452 IF( mod( i, 2 ).NE.0 .AND. i.LT.m ) THEN
453 a( i, i+1 ) = imeps*2
454 ELSE IF( i.GT.1 ) THEN
455 a( i, i-1 ) = -imeps*2
456 END IF
457 END IF
458 240 CONTINUE
459*
460 DO 250 i = 1, n
461 e( i, i ) = one
462 IF( i.LE.4 ) THEN
463 b( i, i ) = -one
464 IF( i.GT.2 )
465 $ b( i, i ) = one - reeps
466 IF( mod( i, 2 ).NE.0 .AND. i.LT.n ) THEN
467 b( i, i+1 ) = imeps
468 ELSE IF( i.GT.1 ) THEN
469 b( i, i-1 ) = -imeps
470 END IF
471 ELSE IF( i.LE.8 ) THEN
472 IF( i.LE.6 ) THEN
473 b( i, i ) = reeps
474 ELSE
475 b( i, i ) = -reeps
476 END IF
477 IF( mod( i, 2 ).NE.0 .AND. i.LT.n ) THEN
478 b( i, i+1 ) = one + imeps
479 ELSE IF( i.GT.1 ) THEN
480 b( i, i-1 ) = -one - imeps
481 END IF
482 ELSE
483 b( i, i ) = one - reeps
484 IF( mod( i, 2 ).NE.0 .AND. i.LT.n ) THEN
485 b( i, i+1 ) = imeps*2
486 ELSE IF( i.GT.1 ) THEN
487 b( i, i-1 ) = -imeps*2
488 END IF
489 END IF
490 250 CONTINUE
491 END IF
492*
493* Compute rhs (C, F)
494*
495 CALL cgemm( 'N', 'N', m, n, m, one, a, lda, r, ldr, zero, c, ldc )
496 CALL cgemm( 'N', 'N', m, n, n, -one, l, ldl, b, ldb, one, c, ldc )
497 CALL cgemm( 'N', 'N', m, n, m, one, d, ldd, r, ldr, zero, f, ldf )
498 CALL cgemm( 'N', 'N', m, n, n, -one, l, ldl, e, lde, one, f, ldf )
499*
500* End of CLATM5
501*
502 END
clatm5
subroutine clatm5(prtype, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, r, ldr, l, ldl, alpha, qblcka, qblckb)
CLATM5
Definition clatm5.f:269
cgemm
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
Definition cgemm.f:188