"Fossies" - the Fresh Open Source Software Archive

Member "lapack-3.9.1/TESTING/EIG/cbdt05.f" (25 Mar 2021, 5735 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 "cbdt05.f": 3.9.0_vs_3.9.1.

    1 *> \brief \b CBDT05
    2 *  =========== DOCUMENTATION ===========
    3 *
    4 * Online html documentation available at
    5 *            http://www.netlib.org/lapack/explore-html/
    6 *
    7 *  Definition:
    8 *  ===========
    9 *
   10 *       SUBROUTINE CBDT05( M, N, A, LDA, S, NS, U, LDU,
   11 *                          VT, LDVT, WORK, RESID )
   12 *
   13 *       .. Scalar Arguments ..
   14 *       INTEGER            LDA, LDU, LDVT, N, NS
   15 *       REAL               RESID
   16 *       ..
   17 *       .. Array Arguments ..
   18 *      REAL               S( * )
   19 *      COMPLEX            A( LDA, * ), U( * ), VT( LDVT, * ), WORK( * )
   20 *       ..
   21 *
   22 *> \par Purpose:
   23 *  =============
   24 *>
   25 *> \verbatim
   26 *>
   27 *> CBDT05 reconstructs a bidiagonal matrix B from its (partial) SVD:
   28 *>    S = U' * B * V
   29 *> where U and V are orthogonal matrices and S is diagonal.
   30 *>
   31 *> The test ratio to test the singular value decomposition is
   32 *>    RESID = norm( S - U' * B * V ) / ( n * norm(B) * EPS )
   33 *> where VT = V' and EPS is the machine precision.
   34 *> \endverbatim
   35 *
   36 *  Arguments:
   37 *  ==========
   38 *
   39 *> \param[in] M
   40 *> \verbatim
   41 *>          M is INTEGER
   42 *>          The number of rows of the matrices A and U.
   43 *> \endverbatim
   44 *>
   45 *> \param[in] N
   46 *> \verbatim
   47 *>          N is INTEGER
   48 *>          The number of columns of the matrices A and VT.
   49 *> \endverbatim
   50 *>
   51 *> \param[in] A
   52 *> \verbatim
   53 *>          A is COMPLEX array, dimension (LDA,N)
   54 *>          The m by n matrix A.
   55 *> \endverbatim
   56 *>
   57 *> \param[in] LDA
   58 *> \verbatim
   59 *>          LDA is INTEGER
   60 *>          The leading dimension of the array A.  LDA >= max(1,M).
   61 *> \endverbatim
   62 *>
   63 *> \param[in] S
   64 *> \verbatim
   65 *>          S is REAL array, dimension (NS)
   66 *>          The singular values from the (partial) SVD of B, sorted in
   67 *>          decreasing order.
   68 *> \endverbatim
   69 *>
   70 *> \param[in] NS
   71 *> \verbatim
   72 *>          NS is INTEGER
   73 *>          The number of singular values/vectors from the (partial)
   74 *>          SVD of B.
   75 *> \endverbatim
   76 *>
   77 *> \param[in] U
   78 *> \verbatim
   79 *>          U is COMPLEX array, dimension (LDU,NS)
   80 *>          The n by ns orthogonal matrix U in S = U' * B * V.
   81 *> \endverbatim
   82 *>
   83 *> \param[in] LDU
   84 *> \verbatim
   85 *>          LDU is INTEGER
   86 *>          The leading dimension of the array U.  LDU >= max(1,N)
   87 *> \endverbatim
   88 *>
   89 *> \param[in] VT
   90 *> \verbatim
   91 *>          VT is COMPLEX array, dimension (LDVT,N)
   92 *>          The n by ns orthogonal matrix V in S = U' * B * V.
   93 *> \endverbatim
   94 *>
   95 *> \param[in] LDVT
   96 *> \verbatim
   97 *>          LDVT is INTEGER
   98 *>          The leading dimension of the array VT.
   99 *> \endverbatim
  100 *>
  101 *> \param[out] WORK
  102 *> \verbatim
  103 *>          WORK is COMPLEX array, dimension (M,N)
  104 *> \endverbatim
  105 *>
  106 *> \param[out] RESID
  107 *> \verbatim
  108 *>          RESID is REAL
  109 *>          The test ratio:  norm(S - U' * A * V) / ( n * norm(A) * EPS )
  110 *> \endverbatim
  111 *
  112 *  Authors:
  113 *  ========
  114 *
  115 *> \author Univ. of Tennessee
  116 *> \author Univ. of California Berkeley
  117 *> \author Univ. of Colorado Denver
  118 *> \author NAG Ltd.
  119 *
  120 *> \ingroup double_eig
  121 *
  122 *  =====================================================================
  123       SUBROUTINE CBDT05( M, N, A, LDA, S, NS, U, LDU,
  124      $                    VT, LDVT, WORK, RESID )
  125 *
  126 *  -- LAPACK test routine --
  127 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  128 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  129 *
  130 *     .. Scalar Arguments ..
  131       INTEGER            LDA, LDU, LDVT, M, N, NS
  132       REAL               RESID
  133 *     ..
  134 *     .. Array Arguments ..
  135       REAL               S( * )
  136       COMPLEX            A( LDA, * ), U( * ), VT( LDVT, * ), WORK( * )
  137 *     ..
  138 *
  139 * ======================================================================
  140 *
  141 *     .. Parameters ..
  142       REAL               ZERO, ONE
  143       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
  144       COMPLEX            CZERO, CONE
  145       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
  146      $                   CONE = ( 1.0E+0, 0.0E+0 ) )
  147 *     ..
  148 *     .. Local Scalars ..
  149       INTEGER            I, J
  150       REAL               ANORM, EPS
  151 *     ..
  152 *     .. Local Arrays ..
  153       REAL               DUM( 1 )
  154 *     ..
  155 *     .. External Functions ..
  156       LOGICAL            LSAME
  157       INTEGER            ISAMAX
  158       REAL               SASUM, SLAMCH, CLANGE
  159       EXTERNAL           LSAME, ISAMAX, SASUM, SLAMCH, CLANGE
  160       REAL               SCASUM
  161 *     ..
  162 *     .. External Subroutines ..
  163       EXTERNAL           CGEMM
  164 *     ..
  165 *     .. Intrinsic Functions ..
  166       INTRINSIC          ABS, REAL, MAX, MIN
  167 *     ..
  168 *     .. Executable Statements ..
  169 *
  170 *     Quick return if possible.
  171 *
  172       RESID = ZERO
  173       IF( MIN( M, N ).LE.0 .OR. NS.LE.0 )
  174      $   RETURN
  175 *
  176       EPS = SLAMCH( 'Precision' )
  177       ANORM = CLANGE( 'M', M, N, A, LDA, DUM )
  178 *
  179 *     Compute U' * A * V.
  180 *
  181       CALL CGEMM( 'N', 'C', M, NS, N, CONE, A, LDA, VT,
  182      $            LDVT, CZERO, WORK( 1+NS*NS ), M )
  183       CALL CGEMM( 'C', 'N', NS, NS, M, -CONE, U, LDU, WORK( 1+NS*NS ),
  184      $            M, CZERO, WORK, NS )
  185 *
  186 *     norm(S - U' * B * V)
  187 *
  188       J = 0
  189       DO 10 I = 1, NS
  190          WORK( J+I ) =  WORK( J+I ) + CMPLX( S( I ), ZERO )
  191          RESID = MAX( RESID, SCASUM( NS, WORK( J+1 ), 1 ) )
  192          J = J + NS
  193    10 CONTINUE
  194 *
  195       IF( ANORM.LE.ZERO ) THEN
  196          IF( RESID.NE.ZERO )
  197      $      RESID = ONE / EPS
  198       ELSE
  199          IF( ANORM.GE.RESID ) THEN
  200             RESID = ( RESID / ANORM ) / ( REAL( N )*EPS )
  201          ELSE
  202             IF( ANORM.LT.ONE ) THEN
  203                RESID = ( MIN( RESID, REAL( N )*ANORM ) / ANORM ) /
  204      $                 ( REAL( N )*EPS )
  205             ELSE
  206                RESID = MIN( RESID / ANORM, REAL( N ) ) /
  207      $                 ( REAL( N )*EPS )
  208             END IF
  209          END IF
  210       END IF
  211 *
  212       RETURN
  213 *
  214 *     End of CBDT05
  215 *
  216       END