LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine dlatm5 ( integer  PRTYPE,
integer  M,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( ldb, * )  B,
integer  LDB,
double precision, dimension( ldc, * )  C,
integer  LDC,
double precision, dimension( ldd, * )  D,
integer  LDD,
double precision, dimension( lde, * )  E,
integer  LDE,
double precision, dimension( ldf, * )  F,
integer  LDF,
double precision, dimension( ldr, * )  R,
integer  LDR,
double precision, dimension( ldl, * )  L,
integer  LDL,
double precision  ALPHA,
integer  QBLCKA,
integer  QBLCKB 
)

DLATM5

Purpose: 
 DLATM5 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 futher 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
          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.
Date
June 2016
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 270 of file dlatm5.f.

270 *
271 * -- LAPACK computational routine (version 3.6.1) --
272 * -- LAPACK is a software package provided by Univ. of Tennessee, --
273 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
274 * June 2016
275 *
276 * .. Scalar Arguments ..
277  INTEGER lda, ldb, ldc, ldd, lde, ldf, ldl, ldr, m, n,
278  $ prtype, qblcka, qblckb
279  DOUBLE PRECISION alpha
280 * ..
281 * .. Array Arguments ..
282  DOUBLE PRECISION 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  DOUBLE PRECISION one, zero, twenty, half, two
291  parameter ( one = 1.0d+0, zero = 0.0d+0, twenty = 2.0d+1,
292  $ half = 0.5d+0, two = 2.0d+0 )
293 * ..
294 * .. Local Scalars ..
295  INTEGER i, j, k
296  DOUBLE PRECISION imeps, reeps
297 * ..
298 * .. Intrinsic Functions ..
299  INTRINSIC dble, mod, sin
300 * ..
301 * .. External Subroutines ..
302  EXTERNAL dgemm
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( dble( 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( dble( i ) ) )*two
349  d( i, j ) = ( half-sin( dble( 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( dble( i+j ) ) )*two
361  e( i, j ) = ( half-sin( dble( 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( dble( i*j ) ) )*twenty
372  l( i, j ) = ( half-sin( dble( 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( dble( i*j ) ) )*twenty
396  d( i, j ) = ( half-sin( dble( 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( dble( i+j ) ) )*twenty
403  e( i, j ) = ( half-sin( dble( 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( dble( j / i ) ) )*twenty
410  l( i, j ) = ( half-sin( dble( 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( dble( i*j ) ) )*alpha / twenty
420  l( i, j ) = ( half-sin( dble( 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 dgemm( 'N', 'N', m, n, m, one, a, lda, r, ldr, zero, c, ldc )
495  CALL dgemm( 'N', 'N', m, n, n, -one, l, ldl, b, ldb, one, c, ldc )
496  CALL dgemm( 'N', 'N', m, n, m, one, d, ldd, r, ldr, zero, f, ldf )
497  CALL dgemm( 'N', 'N', m, n, n, -one, l, ldl, e, lde, one, f, ldf )
498 *
499 * End of DLATM5
500 *
lde
logical function lde(RI, RJ, LR)
Definition: dblat2.f:2945
dgemm
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:189

Here is the call graph for this function:

Here is the caller graph for this function: