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

◆ 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
!>           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.

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