LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ zlatm5()

subroutine zlatm5 ( integer  prtype,
integer  m,
integer  n,
complex*16, dimension( lda, * )  a,
integer  lda,
complex*16, dimension( ldb, * )  b,
integer  ldb,
complex*16, dimension( ldc, * )  c,
integer  ldc,
complex*16, dimension( ldd, * )  d,
integer  ldd,
complex*16, dimension( lde, * )  e,
integer  lde,
complex*16, dimension( ldf, * )  f,
integer  ldf,
complex*16, dimension( ldr, * )  r,
integer  ldr,
complex*16, dimension( ldl, * )  l,
integer  ldl,
double precision  alpha,
integer  qblcka,
integer  qblckb 
)

ZLATM5

Purpose:
 ZLATM5 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 COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 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.
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 zlatm5.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 DOUBLE PRECISION ALPHA
277* ..
278* .. Array Arguments ..
279 COMPLEX*16 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 COMPLEX*16 ONE, TWO, ZERO, HALF, TWENTY
288 parameter( one = ( 1.0d+0, 0.0d+0 ),
289 $ two = ( 2.0d+0, 0.0d+0 ),
290 $ zero = ( 0.0d+0, 0.0d+0 ),
291 $ half = ( 0.5d+0, 0.0d+0 ),
292 $ twenty = ( 2.0d+1, 0.0d+0 ) )
293* ..
294* .. Local Scalars ..
295 INTEGER I, J, K
296 COMPLEX*16 IMEPS, REEPS
297* ..
298* .. Intrinsic Functions ..
299 INTRINSIC dcmplx, mod, sin
300* ..
301* .. External Subroutines ..
302 EXTERNAL zgemm
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( dcmplx( 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( dcmplx( i ) ) )*two
349 d( i, j ) = ( half-sin( dcmplx( 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( dcmplx( i+j ) ) )*two
361 e( i, j ) = ( half-sin( dcmplx( 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( dcmplx( i*j ) ) )*twenty
372 l( i, j ) = ( half-sin( dcmplx( 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( dcmplx( i*j ) ) )*twenty
396 d( i, j ) = ( half-sin( dcmplx( 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( dcmplx( i+j ) ) )*twenty
403 e( i, j ) = ( half-sin( dcmplx( 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( dcmplx( j / i ) ) )*twenty
410 l( i, j ) = ( half-sin( dcmplx( 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( dcmplx( i*j ) ) )*alpha / twenty
420 l( i, j ) = ( half-sin( dcmplx( 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 zgemm( 'N', 'N', m, n, m, one, a, lda, r, ldr, zero, c, ldc )
495 CALL zgemm( 'N', 'N', m, n, n, -one, l, ldl, b, ldb, one, c, ldc )
496 CALL zgemm( 'N', 'N', m, n, m, one, d, ldd, r, ldr, zero, f, ldf )
497 CALL zgemm( 'N', 'N', m, n, n, -one, l, ldl, e, lde, one, f, ldf )
498*
499* End of ZLATM5
500*
lde
logical function lde(ri, rj, lr)
Definition dblat2.f:2970
zgemm
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
Definition zgemm.f:188
Here is the call graph for this function:
Here is the caller graph for this function: