LAPACK  3.4.2 LAPACK: Linear Algebra PACKage
slatm5.f
Go to the documentation of this file.
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