"Fossies" - the Fresh Open Source Software Archive

Member "lapack-3.9.1/TESTING/EIG/cbdt01.f" (25 Mar 2021, 8536 Bytes) of package /linux/misc/lapack-3.9.1.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) Fortran 77 source code syntax highlighting (style: standard) with prefixed line numbers. Alternatively you can here view or download the uninterpreted source code file. See also the latest Fossies "Diffs" side-by-side code changes report for "cbdt01.f": 3.9.0_vs_3.9.1.

    1 *> \brief \b CBDT01
    2 *
    3 *  =========== DOCUMENTATION ===========
    4 *
    5 * Online html documentation available at
    6 *            http://www.netlib.org/lapack/explore-html/
    7 *
    8 *  Definition:
    9 *  ===========
   10 *
   11 *       SUBROUTINE CBDT01( M, N, KD, A, LDA, Q, LDQ, D, E, PT, LDPT, WORK,
   12 *                          RWORK, RESID )
   13 *
   14 *       .. Scalar Arguments ..
   15 *       INTEGER            KD, LDA, LDPT, LDQ, M, N
   16 *       REAL               RESID
   17 *       ..
   18 *       .. Array Arguments ..
   19 *       REAL               D( * ), E( * ), RWORK( * )
   20 *       COMPLEX            A( LDA, * ), PT( LDPT, * ), Q( LDQ, * ),
   21 *      $                   WORK( * )
   22 *       ..
   23 *
   24 *
   25 *> \par Purpose:
   26 *  =============
   27 *>
   28 *> \verbatim
   29 *>
   30 *> CBDT01 reconstructs a general matrix A from its bidiagonal form
   31 *>    A = Q * B * P'
   32 *> where Q (m by min(m,n)) and P' (min(m,n) by n) are unitary
   33 *> matrices and B is bidiagonal.
   34 *>
   35 *> The test ratio to test the reduction is
   36 *>    RESID = norm( A - Q * B * PT ) / ( n * norm(A) * EPS )
   37 *> where PT = P' and EPS is the machine precision.
   38 *> \endverbatim
   39 *
   40 *  Arguments:
   41 *  ==========
   42 *
   43 *> \param[in] M
   44 *> \verbatim
   45 *>          M is INTEGER
   46 *>          The number of rows of the matrices A and Q.
   47 *> \endverbatim
   48 *>
   49 *> \param[in] N
   50 *> \verbatim
   51 *>          N is INTEGER
   52 *>          The number of columns of the matrices A and P'.
   53 *> \endverbatim
   54 *>
   55 *> \param[in] KD
   56 *> \verbatim
   57 *>          KD is INTEGER
   58 *>          If KD = 0, B is diagonal and the array E is not referenced.
   59 *>          If KD = 1, the reduction was performed by xGEBRD; B is upper
   60 *>          bidiagonal if M >= N, and lower bidiagonal if M < N.
   61 *>          If KD = -1, the reduction was performed by xGBBRD; B is
   62 *>          always upper bidiagonal.
   63 *> \endverbatim
   64 *>
   65 *> \param[in] A
   66 *> \verbatim
   67 *>          A is COMPLEX array, dimension (LDA,N)
   68 *>          The m by n matrix A.
   69 *> \endverbatim
   70 *>
   71 *> \param[in] LDA
   72 *> \verbatim
   73 *>          LDA is INTEGER
   74 *>          The leading dimension of the array A.  LDA >= max(1,M).
   75 *> \endverbatim
   76 *>
   77 *> \param[in] Q
   78 *> \verbatim
   79 *>          Q is COMPLEX array, dimension (LDQ,N)
   80 *>          The m by min(m,n) unitary matrix Q in the reduction
   81 *>          A = Q * B * P'.
   82 *> \endverbatim
   83 *>
   84 *> \param[in] LDQ
   85 *> \verbatim
   86 *>          LDQ is INTEGER
   87 *>          The leading dimension of the array Q.  LDQ >= max(1,M).
   88 *> \endverbatim
   89 *>
   90 *> \param[in] D
   91 *> \verbatim
   92 *>          D is REAL array, dimension (min(M,N))
   93 *>          The diagonal elements of the bidiagonal matrix B.
   94 *> \endverbatim
   95 *>
   96 *> \param[in] E
   97 *> \verbatim
   98 *>          E is REAL array, dimension (min(M,N)-1)
   99 *>          The superdiagonal elements of the bidiagonal matrix B if
  100 *>          m >= n, or the subdiagonal elements of B if m < n.
  101 *> \endverbatim
  102 *>
  103 *> \param[in] PT
  104 *> \verbatim
  105 *>          PT is COMPLEX array, dimension (LDPT,N)
  106 *>          The min(m,n) by n unitary matrix P' in the reduction
  107 *>          A = Q * B * P'.
  108 *> \endverbatim
  109 *>
  110 *> \param[in] LDPT
  111 *> \verbatim
  112 *>          LDPT is INTEGER
  113 *>          The leading dimension of the array PT.
  114 *>          LDPT >= max(1,min(M,N)).
  115 *> \endverbatim
  116 *>
  117 *> \param[out] WORK
  118 *> \verbatim
  119 *>          WORK is COMPLEX array, dimension (M+N)
  120 *> \endverbatim
  121 *>
  122 *> \param[out] RWORK
  123 *> \verbatim
  124 *>          RWORK is REAL array, dimension (M)
  125 *> \endverbatim
  126 *>
  127 *> \param[out] RESID
  128 *> \verbatim
  129 *>          RESID is REAL
  130 *>          The test ratio:  norm(A - Q * B * P') / ( n * norm(A) * EPS )
  131 *> \endverbatim
  132 *
  133 *  Authors:
  134 *  ========
  135 *
  136 *> \author Univ. of Tennessee
  137 *> \author Univ. of California Berkeley
  138 *> \author Univ. of Colorado Denver
  139 *> \author NAG Ltd.
  140 *
  141 *> \ingroup complex_eig
  142 *
  143 *  =====================================================================
  144       SUBROUTINE CBDT01( M, N, KD, A, LDA, Q, LDQ, D, E, PT, LDPT, WORK,
  145      $                   RWORK, RESID )
  146 *
  147 *  -- LAPACK test routine --
  148 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  149 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  150 *
  151 *     .. Scalar Arguments ..
  152       INTEGER            KD, LDA, LDPT, LDQ, M, N
  153       REAL               RESID
  154 *     ..
  155 *     .. Array Arguments ..
  156       REAL               D( * ), E( * ), RWORK( * )
  157       COMPLEX            A( LDA, * ), PT( LDPT, * ), Q( LDQ, * ),
  158      $                   WORK( * )
  159 *     ..
  160 *
  161 *  =====================================================================
  162 *
  163 *     .. Parameters ..
  164       REAL               ZERO, ONE
  165       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
  166 *     ..
  167 *     .. Local Scalars ..
  168       INTEGER            I, J
  169       REAL               ANORM, EPS
  170 *     ..
  171 *     .. External Functions ..
  172       REAL               CLANGE, SCASUM, SLAMCH
  173       EXTERNAL           CLANGE, SCASUM, SLAMCH
  174 *     ..
  175 *     .. External Subroutines ..
  176       EXTERNAL           CCOPY, CGEMV
  177 *     ..
  178 *     .. Intrinsic Functions ..
  179       INTRINSIC          CMPLX, MAX, MIN, REAL
  180 *     ..
  181 *     .. Executable Statements ..
  182 *
  183 *     Quick return if possible
  184 *
  185       IF( M.LE.0 .OR. N.LE.0 ) THEN
  186          RESID = ZERO
  187          RETURN
  188       END IF
  189 *
  190 *     Compute A - Q * B * P' one column at a time.
  191 *
  192       RESID = ZERO
  193       IF( KD.NE.0 ) THEN
  194 *
  195 *        B is bidiagonal.
  196 *
  197          IF( KD.NE.0 .AND. M.GE.N ) THEN
  198 *
  199 *           B is upper bidiagonal and M >= N.
  200 *
  201             DO 20 J = 1, N
  202                CALL CCOPY( M, A( 1, J ), 1, WORK, 1 )
  203                DO 10 I = 1, N - 1
  204                   WORK( M+I ) = D( I )*PT( I, J ) + E( I )*PT( I+1, J )
  205    10          CONTINUE
  206                WORK( M+N ) = D( N )*PT( N, J )
  207                CALL CGEMV( 'No transpose', M, N, -CMPLX( ONE ), Q, LDQ,
  208      $                     WORK( M+1 ), 1, CMPLX( ONE ), WORK, 1 )
  209                RESID = MAX( RESID, SCASUM( M, WORK, 1 ) )
  210    20       CONTINUE
  211          ELSE IF( KD.LT.0 ) THEN
  212 *
  213 *           B is upper bidiagonal and M < N.
  214 *
  215             DO 40 J = 1, N
  216                CALL CCOPY( M, A( 1, J ), 1, WORK, 1 )
  217                DO 30 I = 1, M - 1
  218                   WORK( M+I ) = D( I )*PT( I, J ) + E( I )*PT( I+1, J )
  219    30          CONTINUE
  220                WORK( M+M ) = D( M )*PT( M, J )
  221                CALL CGEMV( 'No transpose', M, M, -CMPLX( ONE ), Q, LDQ,
  222      $                     WORK( M+1 ), 1, CMPLX( ONE ), WORK, 1 )
  223                RESID = MAX( RESID, SCASUM( M, WORK, 1 ) )
  224    40       CONTINUE
  225          ELSE
  226 *
  227 *           B is lower bidiagonal.
  228 *
  229             DO 60 J = 1, N
  230                CALL CCOPY( M, A( 1, J ), 1, WORK, 1 )
  231                WORK( M+1 ) = D( 1 )*PT( 1, J )
  232                DO 50 I = 2, M
  233                   WORK( M+I ) = E( I-1 )*PT( I-1, J ) +
  234      $                          D( I )*PT( I, J )
  235    50          CONTINUE
  236                CALL CGEMV( 'No transpose', M, M, -CMPLX( ONE ), Q, LDQ,
  237      $                     WORK( M+1 ), 1, CMPLX( ONE ), WORK, 1 )
  238                RESID = MAX( RESID, SCASUM( M, WORK, 1 ) )
  239    60       CONTINUE
  240          END IF
  241       ELSE
  242 *
  243 *        B is diagonal.
  244 *
  245          IF( M.GE.N ) THEN
  246             DO 80 J = 1, N
  247                CALL CCOPY( M, A( 1, J ), 1, WORK, 1 )
  248                DO 70 I = 1, N
  249                   WORK( M+I ) = D( I )*PT( I, J )
  250    70          CONTINUE
  251                CALL CGEMV( 'No transpose', M, N, -CMPLX( ONE ), Q, LDQ,
  252      $                     WORK( M+1 ), 1, CMPLX( ONE ), WORK, 1 )
  253                RESID = MAX( RESID, SCASUM( M, WORK, 1 ) )
  254    80       CONTINUE
  255          ELSE
  256             DO 100 J = 1, N
  257                CALL CCOPY( M, A( 1, J ), 1, WORK, 1 )
  258                DO 90 I = 1, M
  259                   WORK( M+I ) = D( I )*PT( I, J )
  260    90          CONTINUE
  261                CALL CGEMV( 'No transpose', M, M, -CMPLX( ONE ), Q, LDQ,
  262      $                     WORK( M+1 ), 1, CMPLX( ONE ), WORK, 1 )
  263                RESID = MAX( RESID, SCASUM( M, WORK, 1 ) )
  264   100       CONTINUE
  265          END IF
  266       END IF
  267 *
  268 *     Compute norm(A - Q * B * P') / ( n * norm(A) * EPS )
  269 *
  270       ANORM = CLANGE( '1', M, N, A, LDA, RWORK )
  271       EPS = SLAMCH( 'Precision' )
  272 *
  273       IF( ANORM.LE.ZERO ) THEN
  274          IF( RESID.NE.ZERO )
  275      $      RESID = ONE / EPS
  276       ELSE
  277          IF( ANORM.GE.RESID ) THEN
  278             RESID = ( RESID / ANORM ) / ( REAL( N )*EPS )
  279          ELSE
  280             IF( ANORM.LT.ONE ) THEN
  281                RESID = ( MIN( RESID, REAL( N )*ANORM ) / ANORM ) /
  282      $                 ( REAL( N )*EPS )
  283             ELSE
  284                RESID = MIN( RESID / ANORM, REAL( N ) ) /
  285      $                 ( REAL( N )*EPS )
  286             END IF
  287          END IF
  288       END IF
  289 *
  290       RETURN
  291 *
  292 *     End of CBDT01
  293 *
  294       END