"Fossies" - the Fresh Open Source Software Archive

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

    1 *> \brief \b CBDT03
    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 CBDT03( UPLO, N, KD, D, E, U, LDU, S, VT, LDVT, WORK,
   12 *                          RESID )
   13 *
   14 *       .. Scalar Arguments ..
   15 *       CHARACTER          UPLO
   16 *       INTEGER            KD, LDU, LDVT, N
   17 *       REAL               RESID
   18 *       ..
   19 *       .. Array Arguments ..
   20 *       REAL               D( * ), E( * ), S( * )
   21 *       COMPLEX            U( LDU, * ), VT( LDVT, * ), WORK( * )
   22 *       ..
   23 *
   24 *
   25 *> \par Purpose:
   26 *  =============
   27 *>
   28 *> \verbatim
   29 *>
   30 *> CBDT03 reconstructs a bidiagonal matrix B from its SVD:
   31 *>    S = U' * B * V
   32 *> where U and V are orthogonal matrices and S is diagonal.
   33 *>
   34 *> The test ratio to test the singular value decomposition is
   35 *>    RESID = norm( B - U * S * VT ) / ( n * norm(B) * EPS )
   36 *> where VT = V' and EPS is the machine precision.
   37 *> \endverbatim
   38 *
   39 *  Arguments:
   40 *  ==========
   41 *
   42 *> \param[in] UPLO
   43 *> \verbatim
   44 *>          UPLO is CHARACTER*1
   45 *>          Specifies whether the matrix B is upper or lower bidiagonal.
   46 *>          = 'U':  Upper bidiagonal
   47 *>          = 'L':  Lower bidiagonal
   48 *> \endverbatim
   49 *>
   50 *> \param[in] N
   51 *> \verbatim
   52 *>          N is INTEGER
   53 *>          The order of the matrix B.
   54 *> \endverbatim
   55 *>
   56 *> \param[in] KD
   57 *> \verbatim
   58 *>          KD is INTEGER
   59 *>          The bandwidth of the bidiagonal matrix B.  If KD = 1, the
   60 *>          matrix B is bidiagonal, and if KD = 0, B is diagonal and E is
   61 *>          not referenced.  If KD is greater than 1, it is assumed to be
   62 *>          1, and if KD is less than 0, it is assumed to be 0.
   63 *> \endverbatim
   64 *>
   65 *> \param[in] D
   66 *> \verbatim
   67 *>          D is REAL array, dimension (N)
   68 *>          The n diagonal elements of the bidiagonal matrix B.
   69 *> \endverbatim
   70 *>
   71 *> \param[in] E
   72 *> \verbatim
   73 *>          E is REAL array, dimension (N-1)
   74 *>          The (n-1) superdiagonal elements of the bidiagonal matrix B
   75 *>          if UPLO = 'U', or the (n-1) subdiagonal elements of B if
   76 *>          UPLO = 'L'.
   77 *> \endverbatim
   78 *>
   79 *> \param[in] U
   80 *> \verbatim
   81 *>          U is COMPLEX array, dimension (LDU,N)
   82 *>          The n by n orthogonal matrix U in the reduction B = U'*A*P.
   83 *> \endverbatim
   84 *>
   85 *> \param[in] LDU
   86 *> \verbatim
   87 *>          LDU is INTEGER
   88 *>          The leading dimension of the array U.  LDU >= max(1,N)
   89 *> \endverbatim
   90 *>
   91 *> \param[in] S
   92 *> \verbatim
   93 *>          S is REAL array, dimension (N)
   94 *>          The singular values from the SVD of B, sorted in decreasing
   95 *>          order.
   96 *> \endverbatim
   97 *>
   98 *> \param[in] VT
   99 *> \verbatim
  100 *>          VT is COMPLEX array, dimension (LDVT,N)
  101 *>          The n by n orthogonal matrix V' in the reduction
  102 *>          B = U * S * V'.
  103 *> \endverbatim
  104 *>
  105 *> \param[in] LDVT
  106 *> \verbatim
  107 *>          LDVT is INTEGER
  108 *>          The leading dimension of the array VT.
  109 *> \endverbatim
  110 *>
  111 *> \param[out] WORK
  112 *> \verbatim
  113 *>          WORK is COMPLEX array, dimension (2*N)
  114 *> \endverbatim
  115 *>
  116 *> \param[out] RESID
  117 *> \verbatim
  118 *>          RESID is REAL
  119 *>          The test ratio:  norm(B - U * S * V') / ( n * norm(A) * EPS )
  120 *> \endverbatim
  121 *
  122 *  Authors:
  123 *  ========
  124 *
  125 *> \author Univ. of Tennessee
  126 *> \author Univ. of California Berkeley
  127 *> \author Univ. of Colorado Denver
  128 *> \author NAG Ltd.
  129 *
  130 *> \ingroup complex_eig
  131 *
  132 *  =====================================================================
  133       SUBROUTINE CBDT03( UPLO, N, KD, D, E, U, LDU, S, VT, LDVT, WORK,
  134      $                   RESID )
  135 *
  136 *  -- LAPACK test routine --
  137 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  138 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  139 *
  140 *     .. Scalar Arguments ..
  141       CHARACTER          UPLO
  142       INTEGER            KD, LDU, LDVT, N
  143       REAL               RESID
  144 *     ..
  145 *     .. Array Arguments ..
  146       REAL               D( * ), E( * ), S( * )
  147       COMPLEX            U( LDU, * ), VT( LDVT, * ), WORK( * )
  148 *     ..
  149 *
  150 * ======================================================================
  151 *
  152 *     .. Parameters ..
  153       REAL               ZERO, ONE
  154       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
  155 *     ..
  156 *     .. Local Scalars ..
  157       INTEGER            I, J
  158       REAL               BNORM, EPS
  159 *     ..
  160 *     .. External Functions ..
  161       LOGICAL            LSAME
  162       INTEGER            ISAMAX
  163       REAL               SCASUM, SLAMCH
  164       EXTERNAL           LSAME, ISAMAX, SCASUM, SLAMCH
  165 *     ..
  166 *     .. External Subroutines ..
  167       EXTERNAL           CGEMV
  168 *     ..
  169 *     .. Intrinsic Functions ..
  170       INTRINSIC          ABS, CMPLX, MAX, MIN, REAL
  171 *     ..
  172 *     .. Executable Statements ..
  173 *
  174 *     Quick return if possible
  175 *
  176       RESID = ZERO
  177       IF( N.LE.0 )
  178      $   RETURN
  179 *
  180 *     Compute B - U * S * V' one column at a time.
  181 *
  182       BNORM = ZERO
  183       IF( KD.GE.1 ) THEN
  184 *
  185 *        B is bidiagonal.
  186 *
  187          IF( LSAME( UPLO, 'U' ) ) THEN
  188 *
  189 *           B is upper bidiagonal.
  190 *
  191             DO 20 J = 1, N
  192                DO 10 I = 1, N
  193                   WORK( N+I ) = S( I )*VT( I, J )
  194    10          CONTINUE
  195                CALL CGEMV( 'No transpose', N, N, -CMPLX( ONE ), U, LDU,
  196      $                     WORK( N+1 ), 1, CMPLX( ZERO ), WORK, 1 )
  197                WORK( J ) = WORK( J ) + D( J )
  198                IF( J.GT.1 ) THEN
  199                   WORK( J-1 ) = WORK( J-1 ) + E( J-1 )
  200                   BNORM = MAX( BNORM, ABS( D( J ) )+ABS( E( J-1 ) ) )
  201                ELSE
  202                   BNORM = MAX( BNORM, ABS( D( J ) ) )
  203                END IF
  204                RESID = MAX( RESID, SCASUM( N, WORK, 1 ) )
  205    20       CONTINUE
  206          ELSE
  207 *
  208 *           B is lower bidiagonal.
  209 *
  210             DO 40 J = 1, N
  211                DO 30 I = 1, N
  212                   WORK( N+I ) = S( I )*VT( I, J )
  213    30          CONTINUE
  214                CALL CGEMV( 'No transpose', N, N, -CMPLX( ONE ), U, LDU,
  215      $                     WORK( N+1 ), 1, CMPLX( ZERO ), WORK, 1 )
  216                WORK( J ) = WORK( J ) + D( J )
  217                IF( J.LT.N ) THEN
  218                   WORK( J+1 ) = WORK( J+1 ) + E( J )
  219                   BNORM = MAX( BNORM, ABS( D( J ) )+ABS( E( J ) ) )
  220                ELSE
  221                   BNORM = MAX( BNORM, ABS( D( J ) ) )
  222                END IF
  223                RESID = MAX( RESID, SCASUM( N, WORK, 1 ) )
  224    40       CONTINUE
  225          END IF
  226       ELSE
  227 *
  228 *        B is diagonal.
  229 *
  230          DO 60 J = 1, N
  231             DO 50 I = 1, N
  232                WORK( N+I ) = S( I )*VT( I, J )
  233    50       CONTINUE
  234             CALL CGEMV( 'No transpose', N, N, -CMPLX( ONE ), U, LDU,
  235      $                  WORK( N+1 ), 1, CMPLX( ZERO ), WORK, 1 )
  236             WORK( J ) = WORK( J ) + D( J )
  237             RESID = MAX( RESID, SCASUM( N, WORK, 1 ) )
  238    60    CONTINUE
  239          J = ISAMAX( N, D, 1 )
  240          BNORM = ABS( D( J ) )
  241       END IF
  242 *
  243 *     Compute norm(B - U * S * V') / ( n * norm(B) * EPS )
  244 *
  245       EPS = SLAMCH( 'Precision' )
  246 *
  247       IF( BNORM.LE.ZERO ) THEN
  248          IF( RESID.NE.ZERO )
  249      $      RESID = ONE / EPS
  250       ELSE
  251          IF( BNORM.GE.RESID ) THEN
  252             RESID = ( RESID / BNORM ) / ( REAL( N )*EPS )
  253          ELSE
  254             IF( BNORM.LT.ONE ) THEN
  255                RESID = ( MIN( RESID, REAL( N )*BNORM ) / BNORM ) /
  256      $                 ( REAL( N )*EPS )
  257             ELSE
  258                RESID = MIN( RESID / BNORM, REAL( N ) ) /
  259      $                 ( REAL( N )*EPS )
  260             END IF
  261          END IF
  262       END IF
  263 *
  264       RETURN
  265 *
  266 *     End of CBDT03
  267 *
  268       END