dlatm5.f

Go to the documentation of this file.

*> \brief \b DLATM5
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at 
*            http://www.netlib.org/lapack/explore-html/ 
*
*  Definition:
*  ===========
*
*       SUBROUTINE DLATM5( PRTYPE, M, N, A, LDA, B, LDB, C, LDC, D, LDD,
*                          E, LDE, F, LDF, R, LDR, L, LDL, ALPHA, QBLCKA,
*                          QBLCKB )
* 
*       .. Scalar Arguments ..
*       INTEGER            LDA, LDB, LDC, LDD, LDE, LDF, LDL, LDR, M, N,
*      $                   PRTYPE, QBLCKA, QBLCKB
*       DOUBLE PRECISION   ALPHA
*       ..
*       .. Array Arguments ..
*       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), C( LDC, * ),
*      $                   D( LDD, * ), E( LDE, * ), F( LDF, * ),
*      $                   L( LDL, * ), R( LDR, * )
*       ..
*  
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> 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 ] )
*>
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] PRTYPE
*> \verbatim
*>          PRTYPE is INTEGER
*>          "Points" to a certain type of the matrices to generate
*>          (see futher details).
*> \endverbatim
*>
*> \param[in] M
*> \verbatim
*>          M is INTEGER
*>          Specifies the order of A and D and the number of rows in
*>          C, F,  R and L.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          Specifies the order of B and E and the number of columns in
*>          C, F, R and L.
*> \endverbatim
*>
*> \param[out] A
*> \verbatim
*>          A is DOUBLE PRECISION array, dimension (LDA, M).
*>          On exit A M-by-M is initialized according to PRTYPE.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER
*>          The leading dimension of A.
*> \endverbatim
*>
*> \param[out] B
*> \verbatim
*>          B is DOUBLE PRECISION array, dimension (LDB, N).
*>          On exit B N-by-N is initialized according to PRTYPE.
*> \endverbatim
*>
*> \param[in] LDB
*> \verbatim
*>          LDB is INTEGER
*>          The leading dimension of B.
*> \endverbatim
*>
*> \param[out] C
*> \verbatim
*>          C is DOUBLE PRECISION array, dimension (LDC, N).
*>          On exit C M-by-N is initialized according to PRTYPE.
*> \endverbatim
*>
*> \param[in] LDC
*> \verbatim
*>          LDC is INTEGER
*>          The leading dimension of C.
*> \endverbatim
*>
*> \param[out] D
*> \verbatim
*>          D is DOUBLE PRECISION array, dimension (LDD, M).
*>          On exit D M-by-M is initialized according to PRTYPE.
*> \endverbatim
*>
*> \param[in] LDD
*> \verbatim
*>          LDD is INTEGER
*>          The leading dimension of D.
*> \endverbatim
*>
*> \param[out] E
*> \verbatim
*>          E is DOUBLE PRECISION array, dimension (LDE, N).
*>          On exit E N-by-N is initialized according to PRTYPE.
*> \endverbatim
*>
*> \param[in] LDE
*> \verbatim
*>          LDE is INTEGER
*>          The leading dimension of E.
*> \endverbatim
*>
*> \param[out] F
*> \verbatim
*>          F is DOUBLE PRECISION array, dimension (LDF, N).
*>          On exit F M-by-N is initialized according to PRTYPE.
*> \endverbatim
*>
*> \param[in] LDF
*> \verbatim
*>          LDF is INTEGER
*>          The leading dimension of F.
*> \endverbatim
*>
*> \param[out] R
*> \verbatim
*>          R is DOUBLE PRECISION array, dimension (LDR, N).
*>          On exit R M-by-N is initialized according to PRTYPE.
*> \endverbatim
*>
*> \param[in] LDR
*> \verbatim
*>          LDR is INTEGER
*>          The leading dimension of R.
*> \endverbatim
*>
*> \param[out] L
*> \verbatim
*>          L is DOUBLE PRECISION array, dimension (LDL, N).
*>          On exit L M-by-N is initialized according to PRTYPE.
*> \endverbatim
*>
*> \param[in] LDL
*> \verbatim
*>          LDL is INTEGER
*>          The leading dimension of L.
*> \endverbatim
*>
*> \param[in] ALPHA
*> \verbatim
*>          ALPHA is DOUBLE PRECISION
*>          Parameter used in generating PRTYPE = 1 and 5 matrices.
*> \endverbatim
*>
*> \param[in] QBLCKA
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[in] QBLCKB
*> \verbatim
*>          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.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee 
*> \author Univ. of California Berkeley 
*> \author Univ. of Colorado Denver 
*> \author NAG Ltd. 
*
*> \date June 2016
*
*> \ingroup double_matgen
*
*> \par Further Details:
*  =====================
*>
*> \verbatim
*>
*>  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.
*> \endverbatim
*>
*  =====================================================================
      SUBROUTINE dlatm5( PRTYPE, M, N, A, LDA, B, LDB, C, LDC, D, LDD,
     $                   e, lde, f, ldf, r, ldr, l, ldl, alpha, qblcka,
     $                   qblckb )
*
*  -- LAPACK computational routine (version 3.6.1) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     June 2016
*
*     .. Scalar Arguments ..
      INTEGER            LDA, LDB, LDC, LDD, LDE, LDF, LDL, LDR, M, N,
     $                   prtype, qblcka, qblckb
      DOUBLE PRECISION   ALPHA
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   A( lda, * ), B( ldb, * ), C( ldc, * ),
     $                   d( ldd, * ), e( lde, * ), f( ldf, * ),
     $                   l( ldl, * ), r( ldr, * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ONE, ZERO, TWENTY, HALF, TWO
      parameter                ( one = 1.0d+0, zero = 0.0d+0, twenty = 2.0d+1,
     $                   half = 0.5d+0, two = 2.0d+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            I, J, K
      DOUBLE PRECISION   IMEPS, REEPS
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          dble, mod, sin
*     ..
*     .. External Subroutines ..
      EXTERNAL           dgemm
*     ..
*     .. Executable Statements ..
*
      IF( prtype.EQ.1 ) THEN
         DO 20 i = 1, m
            DO 10 j = 1, m
               IF( i.EQ.j ) THEN
                  a( i, j ) = one
                  d( i, j ) = one
               ELSE IF( i.EQ.j-1 ) THEN
                  a( i, j ) = -one
                  d( i, j ) = zero
               ELSE
                  a( i, j ) = zero
                  d( i, j ) = zero
               END IF
   10       CONTINUE
   20    CONTINUE
*
         DO 40 i = 1, n
            DO 30 j = 1, n
               IF( i.EQ.j ) THEN
                  b( i, j ) = one - alpha
                  e( i, j ) = one
               ELSE IF( i.EQ.j-1 ) THEN
                  b( i, j ) = one
                  e( i, j ) = zero
               ELSE
                  b( i, j ) = zero
                  e( i, j ) = zero
               END IF
   30       CONTINUE
   40    CONTINUE
*
         DO 60 i = 1, m
            DO 50 j = 1, n
               r( i, j ) = ( half-sin( dble( i / j ) ) )*twenty
               l( i, j ) = r( i, j )
   50       CONTINUE
   60    CONTINUE
*
      ELSE IF( prtype.EQ.2 .OR. prtype.EQ.3 ) THEN
         DO 80 i = 1, m
            DO 70 j = 1, m
               IF( i.LE.j ) THEN
                  a( i, j ) = ( half-sin( dble( i ) ) )*two
                  d( i, j ) = ( half-sin( dble( i*j ) ) )*two
               ELSE
                  a( i, j ) = zero
                  d( i, j ) = zero
               END IF
   70       CONTINUE
   80    CONTINUE
*
         DO 100 i = 1, n
            DO 90 j = 1, n
               IF( i.LE.j ) THEN
                  b( i, j ) = ( half-sin( dble( i+j ) ) )*two
                  e( i, j ) = ( half-sin( dble( j ) ) )*two
               ELSE
                  b( i, j ) = zero
                  e( i, j ) = zero
               END IF
   90       CONTINUE
  100    CONTINUE
*
         DO 120 i = 1, m
            DO 110 j = 1, n
               r( i, j ) = ( half-sin( dble( i*j ) ) )*twenty
               l( i, j ) = ( half-sin( dble( i+j ) ) )*twenty
  110       CONTINUE
  120    CONTINUE
*
         IF( prtype.EQ.3 ) THEN
            IF( qblcka.LE.1 )
     $         qblcka = 2
            DO 130 k = 1, m - 1, qblcka
               a( k+1, k+1 ) = a( k, k )
               a( k+1, k ) = -sin( a( k, k+1 ) )
  130       CONTINUE
*
            IF( qblckb.LE.1 )
     $         qblckb = 2
            DO 140 k = 1, n - 1, qblckb
               b( k+1, k+1 ) = b( k, k )
               b( k+1, k ) = -sin( b( k, k+1 ) )
  140       CONTINUE
         END IF
*
      ELSE IF( prtype.EQ.4 ) THEN
         DO 160 i = 1, m
            DO 150 j = 1, m
               a( i, j ) = ( half-sin( dble( i*j ) ) )*twenty
               d( i, j ) = ( half-sin( dble( i+j ) ) )*two
  150       CONTINUE
  160    CONTINUE
*
         DO 180 i = 1, n
            DO 170 j = 1, n
               b( i, j ) = ( half-sin( dble( i+j ) ) )*twenty
               e( i, j ) = ( half-sin( dble( i*j ) ) )*two
  170       CONTINUE
  180    CONTINUE
*
         DO 200 i = 1, m
            DO 190 j = 1, n
               r( i, j ) = ( half-sin( dble( j / i ) ) )*twenty
               l( i, j ) = ( half-sin( dble( i*j ) ) )*two
  190       CONTINUE
  200    CONTINUE
*
      ELSE IF( prtype.GE.5 ) THEN
         reeps = half*two*twenty / alpha
         imeps = ( half-two ) / alpha
         DO 220 i = 1, m
            DO 210 j = 1, n
               r( i, j ) = ( half-sin( dble( i*j ) ) )*alpha / twenty
               l( i, j ) = ( half-sin( dble( i+j ) ) )*alpha / twenty
  210       CONTINUE
  220    CONTINUE
*
         DO 230 i = 1, m
            d( i, i ) = one
  230    CONTINUE
*
         DO 240 i = 1, m
            IF( i.LE.4 ) THEN
               a( i, i ) = one
               IF( i.GT.2 )
     $            a( i, i ) = one + reeps
               IF( mod( i, 2 ).NE.0 .AND. i.LT.m ) THEN
                  a( i, i+1 ) = imeps
               ELSE IF( i.GT.1 ) THEN
                  a( i, i-1 ) = -imeps
               END IF
            ELSE IF( i.LE.8 ) THEN
               IF( i.LE.6 ) THEN
                  a( i, i ) = reeps
               ELSE
                  a( i, i ) = -reeps
               END IF
               IF( mod( i, 2 ).NE.0 .AND. i.LT.m ) THEN
                  a( i, i+1 ) = one
               ELSE IF( i.GT.1 ) THEN
                  a( i, i-1 ) = -one
               END IF
            ELSE
               a( i, i ) = one
               IF( mod( i, 2 ).NE.0 .AND. i.LT.m ) THEN
                  a( i, i+1 ) = imeps*2
               ELSE IF( i.GT.1 ) THEN
                  a( i, i-1 ) = -imeps*2
               END IF
            END IF
  240    CONTINUE
*
         DO 250 i = 1, n
            e( i, i ) = one
            IF( i.LE.4 ) THEN
               b( i, i ) = -one
               IF( i.GT.2 )
     $            b( i, i ) = one - reeps
               IF( mod( i, 2 ).NE.0 .AND. i.LT.n ) THEN
                  b( i, i+1 ) = imeps
               ELSE IF( i.GT.1 ) THEN
                  b( i, i-1 ) = -imeps
               END IF
            ELSE IF( i.LE.8 ) THEN
               IF( i.LE.6 ) THEN
                  b( i, i ) = reeps
               ELSE
                  b( i, i ) = -reeps
               END IF
               IF( mod( i, 2 ).NE.0 .AND. i.LT.n ) THEN
                  b( i, i+1 ) = one + imeps
               ELSE IF( i.GT.1 ) THEN
                  b( i, i-1 ) = -one - imeps
               END IF
            ELSE
               b( i, i ) = one - reeps
               IF( mod( i, 2 ).NE.0 .AND. i.LT.n ) THEN
                  b( i, i+1 ) = imeps*2
               ELSE IF( i.GT.1 ) THEN
                  b( i, i-1 ) = -imeps*2
               END IF
            END IF
  250    CONTINUE
      END IF
*
*     Compute rhs (C, F)
*
      CALL dgemm( 'N', 'N', m, n, m, one, a, lda, r, ldr, zero, c, ldc )
      CALL dgemm( 'N', 'N', m, n, n, -one, l, ldl, b, ldb, one, c, ldc )
      CALL dgemm( 'N', 'N', m, n, m, one, d, ldd, r, ldr, zero, f, ldf )
      CALL dgemm( 'N', 'N', m, n, n, -one, l, ldl, e, lde, one, f, ldf )
*
*     End of DLATM5
*
      END