LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ slatm5()

subroutine slatm5 ( integer  PRTYPE,
integer  M,
integer  N,
real, dimension( lda, * )  A,
integer  LDA,
real, dimension( ldb, * )  B,
integer  LDB,
real, dimension( ldc, * )  C,
integer  LDC,
real, dimension( ldd, * )  D,
integer  LDD,
real, dimension( lde, * )  E,
integer  LDE,
real, dimension( ldf, * )  F,
integer  LDF,
real, dimension( ldr, * )  R,
integer  LDR,
real, dimension( ldl, * )  L,
integer  LDL,
real  ALPHA,
integer  QBLCKA,
integer  QBLCKB 
)

SLATM5

Purpose:
 SLATM5 generates matrices involved in the Generalized Sylvester
 equation:

     A * R - L * B = C
     D * R - L * E = F

 They also satisfy (the diagonalization condition)

  [ I -L ] ( [ A  -C ], [ D -F ] ) [ I  R ] = ( [ A    ], [ D    ] )
  [    I ] ( [     B ]  [    E ] ) [    I ]   ( [    B ]  [    E ] )
Parameters
[in]PRTYPE
          PRTYPE is INTEGER
          "Points" to a certain type of the matrices to generate
          (see further details).
[in]M
          M is INTEGER
          Specifies the order of A and D and the number of rows in
          C, F,  R and L.
[in]N
          N is INTEGER
          Specifies the order of B and E and the number of columns in
          C, F, R and L.
[out]A
          A is REAL array, dimension (LDA, M).
          On exit A M-by-M is initialized according to PRTYPE.
[in]LDA
          LDA is INTEGER
          The leading dimension of A.
[out]B
          B is REAL array, dimension (LDB, N).
          On exit B N-by-N is initialized according to PRTYPE.
[in]LDB
          LDB is INTEGER
          The leading dimension of B.
[out]C
          C is REAL array, dimension (LDC, N).
          On exit C M-by-N is initialized according to PRTYPE.
[in]LDC
          LDC is INTEGER
          The leading dimension of C.
[out]D
          D is REAL array, dimension (LDD, M).
          On exit D M-by-M is initialized according to PRTYPE.
[in]LDD
          LDD is INTEGER
          The leading dimension of D.
[out]E
          E is REAL array, dimension (LDE, N).
          On exit E N-by-N is initialized according to PRTYPE.
[in]LDE
          LDE is INTEGER
          The leading dimension of E.
[out]F
          F is REAL array, dimension (LDF, N).
          On exit F M-by-N is initialized according to PRTYPE.
[in]LDF
          LDF is INTEGER
          The leading dimension of F.
[out]R
          R is REAL array, dimension (LDR, N).
          On exit R M-by-N is initialized according to PRTYPE.
[in]LDR
          LDR is INTEGER
          The leading dimension of R.
[out]L
          L is REAL array, dimension (LDL, N).
          On exit L M-by-N is initialized according to PRTYPE.
[in]LDL
          LDL is INTEGER
          The leading dimension of L.
[in]ALPHA
          ALPHA is REAL
          Parameter used in generating PRTYPE = 1 and 5 matrices.
[in]QBLCKA
          QBLCKA is INTEGER
          When PRTYPE = 3, specifies the distance between 2-by-2
          blocks on the diagonal in A. Otherwise, QBLCKA is not
          referenced. QBLCKA > 1.
[in]QBLCKB
          QBLCKB is INTEGER
          When PRTYPE = 3, specifies the distance between 2-by-2
          blocks on the diagonal in B. Otherwise, QBLCKB is not
          referenced. QBLCKB > 1.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
  PRTYPE = 1: A and B are Jordan blocks, D and E are identity matrices

             A : if (i == j) then A(i, j) = 1.0
                 if (j == i + 1) then A(i, j) = -1.0
                 else A(i, j) = 0.0,            i, j = 1...M

             B : if (i == j) then B(i, j) = 1.0 - ALPHA
                 if (j == i + 1) then B(i, j) = 1.0
                 else B(i, j) = 0.0,            i, j = 1...N

             D : if (i == j) then D(i, j) = 1.0
                 else D(i, j) = 0.0,            i, j = 1...M

             E : if (i == j) then E(i, j) = 1.0
                 else E(i, j) = 0.0,            i, j = 1...N

             L =  R are chosen from [-10...10],
                  which specifies the right hand sides (C, F).

  PRTYPE = 2 or 3: Triangular and/or quasi- triangular.

             A : if (i <= j) then A(i, j) = [-1...1]
                 else A(i, j) = 0.0,             i, j = 1...M

                 if (PRTYPE = 3) then
                    A(k + 1, k + 1) = A(k, k)
                    A(k + 1, k) = [-1...1]
                    sign(A(k, k + 1) = -(sin(A(k + 1, k))
                        k = 1, M - 1, QBLCKA

             B : if (i <= j) then B(i, j) = [-1...1]
                 else B(i, j) = 0.0,            i, j = 1...N

                 if (PRTYPE = 3) then
                    B(k + 1, k + 1) = B(k, k)
                    B(k + 1, k) = [-1...1]
                    sign(B(k, k + 1) = -(sign(B(k + 1, k))
                        k = 1, N - 1, QBLCKB

             D : if (i <= j) then D(i, j) = [-1...1].
                 else D(i, j) = 0.0,            i, j = 1...M


             E : if (i <= j) then D(i, j) = [-1...1]
                 else E(i, j) = 0.0,            i, j = 1...N

                 L, R are chosen from [-10...10],
                 which specifies the right hand sides (C, F).

  PRTYPE = 4 Full
             A(i, j) = [-10...10]
             D(i, j) = [-1...1]    i,j = 1...M
             B(i, j) = [-10...10]
             E(i, j) = [-1...1]    i,j = 1...N
             R(i, j) = [-10...10]
             L(i, j) = [-1...1]    i = 1..M ,j = 1...N

             L, R specifies the right hand sides (C, F).

  PRTYPE = 5 special case common and/or close eigs.

Definition at line 265 of file slatm5.f.

268 *
269 * -- LAPACK computational routine --
270 * -- LAPACK is a software package provided by Univ. of Tennessee, --
271 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
272 *
273 * .. Scalar Arguments ..
274  INTEGER LDA, LDB, LDC, LDD, LDE, LDF, LDL, LDR, M, N,
275  $ PRTYPE, QBLCKA, QBLCKB
276  REAL ALPHA
277 * ..
278 * .. Array Arguments ..
279  REAL A( LDA, * ), B( LDB, * ), C( LDC, * ),
280  $ D( LDD, * ), E( LDE, * ), F( LDF, * ),
281  $ L( LDL, * ), R( LDR, * )
282 * ..
283 *
284 * =====================================================================
285 *
286 * .. Parameters ..
287  REAL ONE, ZERO, TWENTY, HALF, TWO
288  parameter( one = 1.0e+0, zero = 0.0e+0, twenty = 2.0e+1,
289  $ half = 0.5e+0, two = 2.0e+0 )
290 * ..
291 * .. Local Scalars ..
292  INTEGER I, J, K
293  REAL IMEPS, REEPS
294 * ..
295 * .. Intrinsic Functions ..
296  INTRINSIC mod, real, sin
297 * ..
298 * .. External Subroutines ..
299  EXTERNAL sgemm
300 * ..
301 * .. Executable Statements ..
302 *
303  IF( prtype.EQ.1 ) THEN
304  DO 20 i = 1, m
305  DO 10 j = 1, m
306  IF( i.EQ.j ) THEN
307  a( i, j ) = one
308  d( i, j ) = one
309  ELSE IF( i.EQ.j-1 ) THEN
310  a( i, j ) = -one
311  d( i, j ) = zero
312  ELSE
313  a( i, j ) = zero
314  d( i, j ) = zero
315  END IF
316  10 CONTINUE
317  20 CONTINUE
318 *
319  DO 40 i = 1, n
320  DO 30 j = 1, n
321  IF( i.EQ.j ) THEN
322  b( i, j ) = one - alpha
323  e( i, j ) = one
324  ELSE IF( i.EQ.j-1 ) THEN
325  b( i, j ) = one
326  e( i, j ) = zero
327  ELSE
328  b( i, j ) = zero
329  e( i, j ) = zero
330  END IF
331  30 CONTINUE
332  40 CONTINUE
333 *
334  DO 60 i = 1, m
335  DO 50 j = 1, n
336  r( i, j ) = ( half-sin( real( i / j ) ) )*twenty
337  l( i, j ) = r( i, j )
338  50 CONTINUE
339  60 CONTINUE
340 *
341  ELSE IF( prtype.EQ.2 .OR. prtype.EQ.3 ) THEN
342  DO 80 i = 1, m
343  DO 70 j = 1, m
344  IF( i.LE.j ) THEN
345  a( i, j ) = ( half-sin( real( i ) ) )*two
346  d( i, j ) = ( half-sin( real( i*j ) ) )*two
347  ELSE
348  a( i, j ) = zero
349  d( i, j ) = zero
350  END IF
351  70 CONTINUE
352  80 CONTINUE
353 *
354  DO 100 i = 1, n
355  DO 90 j = 1, n
356  IF( i.LE.j ) THEN
357  b( i, j ) = ( half-sin( real( i+j ) ) )*two
358  e( i, j ) = ( half-sin( real( j ) ) )*two
359  ELSE
360  b( i, j ) = zero
361  e( i, j ) = zero
362  END IF
363  90 CONTINUE
364  100 CONTINUE
365 *
366  DO 120 i = 1, m
367  DO 110 j = 1, n
368  r( i, j ) = ( half-sin( real( i*j ) ) )*twenty
369  l( i, j ) = ( half-sin( real( i+j ) ) )*twenty
370  110 CONTINUE
371  120 CONTINUE
372 *
373  IF( prtype.EQ.3 ) THEN
374  IF( qblcka.LE.1 )
375  $ qblcka = 2
376  DO 130 k = 1, m - 1, qblcka
377  a( k+1, k+1 ) = a( k, k )
378  a( k+1, k ) = -sin( a( k, k+1 ) )
379  130 CONTINUE
380 *
381  IF( qblckb.LE.1 )
382  $ qblckb = 2
383  DO 140 k = 1, n - 1, qblckb
384  b( k+1, k+1 ) = b( k, k )
385  b( k+1, k ) = -sin( b( k, k+1 ) )
386  140 CONTINUE
387  END IF
388 *
389  ELSE IF( prtype.EQ.4 ) THEN
390  DO 160 i = 1, m
391  DO 150 j = 1, m
392  a( i, j ) = ( half-sin( real( i*j ) ) )*twenty
393  d( i, j ) = ( half-sin( real( i+j ) ) )*two
394  150 CONTINUE
395  160 CONTINUE
396 *
397  DO 180 i = 1, n
398  DO 170 j = 1, n
399  b( i, j ) = ( half-sin( real( i+j ) ) )*twenty
400  e( i, j ) = ( half-sin( real( i*j ) ) )*two
401  170 CONTINUE
402  180 CONTINUE
403 *
404  DO 200 i = 1, m
405  DO 190 j = 1, n
406  r( i, j ) = ( half-sin( real( j / i ) ) )*twenty
407  l( i, j ) = ( half-sin( real( i*j ) ) )*two
408  190 CONTINUE
409  200 CONTINUE
410 *
411  ELSE IF( prtype.GE.5 ) THEN
412  reeps = half*two*twenty / alpha
413  imeps = ( half-two ) / alpha
414  DO 220 i = 1, m
415  DO 210 j = 1, n
416  r( i, j ) = ( half-sin( real( i*j ) ) )*alpha / twenty
417  l( i, j ) = ( half-sin( real( i+j ) ) )*alpha / twenty
418  210 CONTINUE
419  220 CONTINUE
420 *
421  DO 230 i = 1, m
422  d( i, i ) = one
423  230 CONTINUE
424 *
425  DO 240 i = 1, m
426  IF( i.LE.4 ) THEN
427  a( i, i ) = one
428  IF( i.GT.2 )
429  $ a( i, i ) = one + reeps
430  IF( mod( i, 2 ).NE.0 .AND. i.LT.m ) THEN
431  a( i, i+1 ) = imeps
432  ELSE IF( i.GT.1 ) THEN
433  a( i, i-1 ) = -imeps
434  END IF
435  ELSE IF( i.LE.8 ) THEN
436  IF( i.LE.6 ) THEN
437  a( i, i ) = reeps
438  ELSE
439  a( i, i ) = -reeps
440  END IF
441  IF( mod( i, 2 ).NE.0 .AND. i.LT.m ) THEN
442  a( i, i+1 ) = one
443  ELSE IF( i.GT.1 ) THEN
444  a( i, i-1 ) = -one
445  END IF
446  ELSE
447  a( i, i ) = one
448  IF( mod( i, 2 ).NE.0 .AND. i.LT.m ) THEN
449  a( i, i+1 ) = imeps*2
450  ELSE IF( i.GT.1 ) THEN
451  a( i, i-1 ) = -imeps*2
452  END IF
453  END IF
454  240 CONTINUE
455 *
456  DO 250 i = 1, n
457  e( i, i ) = one
458  IF( i.LE.4 ) THEN
459  b( i, i ) = -one
460  IF( i.GT.2 )
461  $ b( i, i ) = one - reeps
462  IF( mod( i, 2 ).NE.0 .AND. i.LT.n ) THEN
463  b( i, i+1 ) = imeps
464  ELSE IF( i.GT.1 ) THEN
465  b( i, i-1 ) = -imeps
466  END IF
467  ELSE IF( i.LE.8 ) THEN
468  IF( i.LE.6 ) THEN
469  b( i, i ) = reeps
470  ELSE
471  b( i, i ) = -reeps
472  END IF
473  IF( mod( i, 2 ).NE.0 .AND. i.LT.n ) THEN
474  b( i, i+1 ) = one + imeps
475  ELSE IF( i.GT.1 ) THEN
476  b( i, i-1 ) = -one - imeps
477  END IF
478  ELSE
479  b( i, i ) = one - reeps
480  IF( mod( i, 2 ).NE.0 .AND. i.LT.n ) THEN
481  b( i, i+1 ) = imeps*2
482  ELSE IF( i.GT.1 ) THEN
483  b( i, i-1 ) = -imeps*2
484  END IF
485  END IF
486  250 CONTINUE
487  END IF
488 *
489 * Compute rhs (C, F)
490 *
491  CALL sgemm( 'N', 'N', m, n, m, one, a, lda, r, ldr, zero, c, ldc )
492  CALL sgemm( 'N', 'N', m, n, n, -one, l, ldl, b, ldb, one, c, ldc )
493  CALL sgemm( 'N', 'N', m, n, m, one, d, ldd, r, ldr, zero, f, ldf )
494  CALL sgemm( 'N', 'N', m, n, n, -one, l, ldl, e, lde, one, f, ldf )
495 *
496 * End of SLATM5
497 *
logical function lde(RI, RJ, LR)
Definition: dblat2.f:2942
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:187
Here is the call graph for this function:
Here is the caller graph for this function: