1*> \brief \b SLATM5
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 SLATM5( 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*       REAL               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*> SLATM5 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 certian type of the matrices to generate
52*>          (see futher 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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*> \date November 2011
196*
197*> \ingroup real_matgen
198*
199*> \par Further Details:
200*  =====================
201*>
202*> \verbatim
203*>
204*>  PRTYPE = 1: A and B are Jordan blocks, D and E are identity matrices
205*>
206*>             A : if (i == j) then A(i, j) = 1.0
207*>                 if (j == i + 1) then A(i, j) = -1.0
208*>                 else A(i, j) = 0.0,            i, j = 1...M
209*>
210*>             B : if (i == j) then B(i, j) = 1.0 - ALPHA
211*>                 if (j == i + 1) then B(i, j) = 1.0
212*>                 else B(i, j) = 0.0,            i, j = 1...N
213*>
214*>             D : if (i == j) then D(i, j) = 1.0
215*>                 else D(i, j) = 0.0,            i, j = 1...M
216*>
217*>             E : if (i == j) then E(i, j) = 1.0
218*>                 else E(i, j) = 0.0,            i, j = 1...N
219*>
220*>             L =  R are chosen from [-10...10],
221*>                  which specifies the right hand sides (C, F).
222*>
223*>  PRTYPE = 2 or 3: Triangular and/or quasi- triangular.
224*>
225*>             A : if (i <= j) then A(i, j) = [-1...1]
226*>                 else A(i, j) = 0.0,             i, j = 1...M
227*>
228*>                 if (PRTYPE = 3) then
229*>                    A(k + 1, k + 1) = A(k, k)
230*>                    A(k + 1, k) = [-1...1]
231*>                    sign(A(k, k + 1) = -(sin(A(k + 1, k))
232*>                        k = 1, M - 1, QBLCKA
233*>
234*>             B : if (i <= j) then B(i, j) = [-1...1]
235*>                 else B(i, j) = 0.0,            i, j = 1...N
236*>
237*>                 if (PRTYPE = 3) then
238*>                    B(k + 1, k + 1) = B(k, k)
239*>                    B(k + 1, k) = [-1...1]
240*>                    sign(B(k, k + 1) = -(sign(B(k + 1, k))
241*>                        k = 1, N - 1, QBLCKB
242*>
243*>             D : if (i <= j) then D(i, j) = [-1...1].
244*>                 else D(i, j) = 0.0,            i, j = 1...M
245*>
246*>
247*>             E : if (i <= j) then D(i, j) = [-1...1]
248*>                 else E(i, j) = 0.0,            i, j = 1...N
249*>
250*>                 L, R are chosen from [-10...10],
251*>                 which specifies the right hand sides (C, F).
252*>
253*>  PRTYPE = 4 Full
254*>             A(i, j) = [-10...10]
255*>             D(i, j) = [-1...1]    i,j = 1...M
256*>             B(i, j) = [-10...10]
257*>             E(i, j) = [-1...1]    i,j = 1...N
258*>             R(i, j) = [-10...10]
259*>             L(i, j) = [-1...1]    i = 1..M ,j = 1...N
260*>
261*>             L, R specifies the right hand sides (C, F).
262*>
263*>  PRTYPE = 5 special case common and/or close eigs.
264*> \endverbatim
265*>
266*  =====================================================================
267      SUBROUTINE SLATM5( PRTYPE, M, N, A, LDA, B, LDB, C, LDC, D, LDD,
268     $                   E, LDE, F, LDF, R, LDR, L, LDL, ALPHA, QBLCKA,
269     $                   QBLCKB )
270*
271*  -- LAPACK computational routine (version 3.4.0) --
272*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
273*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
274*     November 2011
275*
276*     .. Scalar Arguments ..
277      INTEGER            LDA, LDB, LDC, LDD, LDE, LDF, LDL, LDR, M, N,
278     $                   PRTYPE, QBLCKA, QBLCKB
279      REAL               ALPHA
280*     ..
281*     .. Array Arguments ..
282      REAL               A( LDA, * ), B( LDB, * ), C( LDC, * ),
283     $                   D( LDD, * ), E( LDE, * ), F( LDF, * ),
284     $                   L( LDL, * ), R( LDR, * )
285*     ..
286*
287*  =====================================================================
288*
289*     .. Parameters ..
290      REAL               ONE, ZERO, TWENTY, HALF, TWO
291      PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0, TWENTY = 2.0E+1,
292     $                   HALF = 0.5E+0, TWO = 2.0E+0 )
293*     ..
294*     .. Local Scalars ..
295      INTEGER            I, J, K
296      REAL               IMEPS, REEPS
297*     ..
298*     .. Intrinsic Functions ..
299      INTRINSIC          MOD, REAL, SIN
300*     ..
301*     .. External Subroutines ..
302      EXTERNAL           SGEMM
303*     ..
304*     .. Executable Statements ..
305*
306      IF( PRTYPE.EQ.1 ) THEN
307         DO 20 I = 1, M
308            DO 10 J = 1, M
309               IF( I.EQ.J ) THEN
310                  A( I, J ) = ONE
311                  D( I, J ) = ONE
312               ELSE IF( I.EQ.J-1 ) THEN
313                  A( I, J ) = -ONE
314                  D( I, J ) = ZERO
315               ELSE
316                  A( I, J ) = ZERO
317                  D( I, J ) = ZERO
318               END IF
319   10       CONTINUE
320   20    CONTINUE
321*
322         DO 40 I = 1, N
323            DO 30 J = 1, N
324               IF( I.EQ.J ) THEN
325                  B( I, J ) = ONE - ALPHA
326                  E( I, J ) = ONE
327               ELSE IF( I.EQ.J-1 ) THEN
328                  B( I, J ) = ONE
329                  E( I, J ) = ZERO
330               ELSE
331                  B( I, J ) = ZERO
332                  E( I, J ) = ZERO
333               END IF
334   30       CONTINUE
335   40    CONTINUE
336*
337         DO 60 I = 1, M
338            DO 50 J = 1, N
339               R( I, J ) = ( HALF-SIN( REAL( I / J ) ) )*TWENTY
340               L( I, J ) = R( I, J )
341   50       CONTINUE
342   60    CONTINUE
343*
344      ELSE IF( PRTYPE.EQ.2 .OR. PRTYPE.EQ.3 ) THEN
345         DO 80 I = 1, M
346            DO 70 J = 1, M
347               IF( I.LE.J ) THEN
348                  A( I, J ) = ( HALF-SIN( REAL( I ) ) )*TWO
349                  D( I, J ) = ( HALF-SIN( REAL( I*J ) ) )*TWO
350               ELSE
351                  A( I, J ) = ZERO
352                  D( I, J ) = ZERO
353               END IF
354   70       CONTINUE
355   80    CONTINUE
356*
357         DO 100 I = 1, N
358            DO 90 J = 1, N
359               IF( I.LE.J ) THEN
360                  B( I, J ) = ( HALF-SIN( REAL( I+J ) ) )*TWO
361                  E( I, J ) = ( HALF-SIN( REAL( J ) ) )*TWO
362               ELSE
363                  B( I, J ) = ZERO
364                  E( I, J ) = ZERO
365               END IF
366   90       CONTINUE
367  100    CONTINUE
368*
369         DO 120 I = 1, M
370            DO 110 J = 1, N
371               R( I, J ) = ( HALF-SIN( REAL( I*J ) ) )*TWENTY
372               L( I, J ) = ( HALF-SIN( REAL( I+J ) ) )*TWENTY
373  110       CONTINUE
374  120    CONTINUE
375*
376         IF( PRTYPE.EQ.3 ) THEN
377            IF( QBLCKA.LE.1 )
378     $         QBLCKA = 2
379            DO 130 K = 1, M - 1, QBLCKA
380               A( K+1, K+1 ) = A( K, K )
381               A( K+1, K ) = -SIN( A( K, K+1 ) )
382  130       CONTINUE
383*
384            IF( QBLCKB.LE.1 )
385     $         QBLCKB = 2
386            DO 140 K = 1, N - 1, QBLCKB
387               B( K+1, K+1 ) = B( K, K )
388               B( K+1, K ) = -SIN( B( K, K+1 ) )
389  140       CONTINUE
390         END IF
391*
392      ELSE IF( PRTYPE.EQ.4 ) THEN
393         DO 160 I = 1, M
394            DO 150 J = 1, M
395               A( I, J ) = ( HALF-SIN( REAL( I*J ) ) )*TWENTY
396               D( I, J ) = ( HALF-SIN( REAL( I+J ) ) )*TWO
397  150       CONTINUE
398  160    CONTINUE
399*
400         DO 180 I = 1, N
401            DO 170 J = 1, N
402               B( I, J ) = ( HALF-SIN( REAL( I+J ) ) )*TWENTY
403               E( I, J ) = ( HALF-SIN( REAL( I*J ) ) )*TWO
404  170       CONTINUE
405  180    CONTINUE
406*
407         DO 200 I = 1, M
408            DO 190 J = 1, N
409               R( I, J ) = ( HALF-SIN( REAL( J / I ) ) )*TWENTY
410               L( I, J ) = ( HALF-SIN( REAL( I*J ) ) )*TWO
411  190       CONTINUE
412  200    CONTINUE
413*
414      ELSE IF( PRTYPE.GE.5 ) THEN
415         REEPS = HALF*TWO*TWENTY / ALPHA
416         IMEPS = ( HALF-TWO ) / ALPHA
417         DO 220 I = 1, M
418            DO 210 J = 1, N
419               R( I, J ) = ( HALF-SIN( REAL( I*J ) ) )*ALPHA / TWENTY
420               L( I, J ) = ( HALF-SIN( REAL( I+J ) ) )*ALPHA / TWENTY
421  210       CONTINUE
422  220    CONTINUE
423*
424         DO 230 I = 1, M
425            D( I, I ) = ONE
426  230    CONTINUE
427*
428         DO 240 I = 1, M
429            IF( I.LE.4 ) THEN
430               A( I, I ) = ONE
431               IF( I.GT.2 )
432     $            A( I, I ) = ONE + REEPS
433               IF( MOD( I, 2 ).NE.0 .AND. I.LT.M ) THEN
434                  A( I, I+1 ) = IMEPS
435               ELSE IF( I.GT.1 ) THEN
436                  A( I, I-1 ) = -IMEPS
437               END IF
438            ELSE IF( I.LE.8 ) THEN
439               IF( I.LE.6 ) THEN
440                  A( I, I ) = REEPS
441               ELSE
442                  A( I, I ) = -REEPS
443               END IF
444               IF( MOD( I, 2 ).NE.0 .AND. I.LT.M ) THEN
445                  A( I, I+1 ) = ONE
446               ELSE IF( I.GT.1 ) THEN
447                  A( I, I-1 ) = -ONE
448               END IF
449            ELSE
450               A( I, I ) = ONE
451               IF( MOD( I, 2 ).NE.0 .AND. I.LT.M ) THEN
452                  A( I, I+1 ) = IMEPS*2
453               ELSE IF( I.GT.1 ) THEN
454                  A( I, I-1 ) = -IMEPS*2
455               END IF
456            END IF
457  240    CONTINUE
458*
459         DO 250 I = 1, N
460            E( I, I ) = ONE
461            IF( I.LE.4 ) THEN
462               B( I, I ) = -ONE
463               IF( I.GT.2 )
464     $            B( I, I ) = ONE - REEPS
465               IF( MOD( I, 2 ).NE.0 .AND. I.LT.N ) THEN
466                  B( I, I+1 ) = IMEPS
467               ELSE IF( I.GT.1 ) THEN
468                  B( I, I-1 ) = -IMEPS
469               END IF
470            ELSE IF( I.LE.8 ) THEN
471               IF( I.LE.6 ) THEN
472                  B( I, I ) = REEPS
473               ELSE
474                  B( I, I ) = -REEPS
475               END IF
476               IF( MOD( I, 2 ).NE.0 .AND. I.LT.N ) THEN
477                  B( I, I+1 ) = ONE + IMEPS
478               ELSE IF( I.GT.1 ) THEN
479                  B( I, I-1 ) = -ONE - IMEPS
480               END IF
481            ELSE
482               B( I, I ) = ONE - REEPS
483               IF( MOD( I, 2 ).NE.0 .AND. I.LT.N ) THEN
484                  B( I, I+1 ) = IMEPS*2
485               ELSE IF( I.GT.1 ) THEN
486                  B( I, I-1 ) = -IMEPS*2
487               END IF
488            END IF
489  250    CONTINUE
490      END IF
491*
492*     Compute rhs (C, F)
493*
494      CALL SGEMM( 'N', 'N', M, N, M, ONE, A, LDA, R, LDR, ZERO, C, LDC )
495      CALL SGEMM( 'N', 'N', M, N, N, -ONE, L, LDL, B, LDB, ONE, C, LDC )
496      CALL SGEMM( 'N', 'N', M, N, M, ONE, D, LDD, R, LDR, ZERO, F, LDF )
497      CALL SGEMM( 'N', 'N', M, N, N, -ONE, L, LDL, E, LDE, ONE, F, LDF )
498*
499*     End of SLATM5
500*
501      END
502