LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dbdt01.f
Go to the documentation of this file.
00001 *> \brief \b DBDT01
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *  Definition:
00009 *  ===========
00010 *
00011 *       SUBROUTINE DBDT01( M, N, KD, A, LDA, Q, LDQ, D, E, PT, LDPT, WORK,
00012 *                          RESID )
00013 * 
00014 *       .. Scalar Arguments ..
00015 *       INTEGER            KD, LDA, LDPT, LDQ, M, N
00016 *       DOUBLE PRECISION   RESID
00017 *       ..
00018 *       .. Array Arguments ..
00019 *       DOUBLE PRECISION   A( LDA, * ), D( * ), E( * ), PT( LDPT, * ),
00020 *      $                   Q( LDQ, * ), WORK( * )
00021 *       ..
00022 *  
00023 *
00024 *> \par Purpose:
00025 *  =============
00026 *>
00027 *> \verbatim
00028 *>
00029 *> DBDT01 reconstructs a general matrix A from its bidiagonal form
00030 *>    A = Q * B * P'
00031 *> where Q (m by min(m,n)) and P' (min(m,n) by n) are orthogonal
00032 *> matrices and B is bidiagonal.
00033 *>
00034 *> The test ratio to test the reduction is
00035 *>    RESID = norm( A - Q * B * PT ) / ( n * norm(A) * EPS )
00036 *> where PT = P' and EPS is the machine precision.
00037 *> \endverbatim
00038 *
00039 *  Arguments:
00040 *  ==========
00041 *
00042 *> \param[in] M
00043 *> \verbatim
00044 *>          M is INTEGER
00045 *>          The number of rows of the matrices A and Q.
00046 *> \endverbatim
00047 *>
00048 *> \param[in] N
00049 *> \verbatim
00050 *>          N is INTEGER
00051 *>          The number of columns of the matrices A and P'.
00052 *> \endverbatim
00053 *>
00054 *> \param[in] KD
00055 *> \verbatim
00056 *>          KD is INTEGER
00057 *>          If KD = 0, B is diagonal and the array E is not referenced.
00058 *>          If KD = 1, the reduction was performed by xGEBRD; B is upper
00059 *>          bidiagonal if M >= N, and lower bidiagonal if M < N.
00060 *>          If KD = -1, the reduction was performed by xGBBRD; B is
00061 *>          always upper bidiagonal.
00062 *> \endverbatim
00063 *>
00064 *> \param[in] A
00065 *> \verbatim
00066 *>          A is DOUBLE PRECISION array, dimension (LDA,N)
00067 *>          The m by n matrix A.
00068 *> \endverbatim
00069 *>
00070 *> \param[in] LDA
00071 *> \verbatim
00072 *>          LDA is INTEGER
00073 *>          The leading dimension of the array A.  LDA >= max(1,M).
00074 *> \endverbatim
00075 *>
00076 *> \param[in] Q
00077 *> \verbatim
00078 *>          Q is DOUBLE PRECISION array, dimension (LDQ,N)
00079 *>          The m by min(m,n) orthogonal matrix Q in the reduction
00080 *>          A = Q * B * P'.
00081 *> \endverbatim
00082 *>
00083 *> \param[in] LDQ
00084 *> \verbatim
00085 *>          LDQ is INTEGER
00086 *>          The leading dimension of the array Q.  LDQ >= max(1,M).
00087 *> \endverbatim
00088 *>
00089 *> \param[in] D
00090 *> \verbatim
00091 *>          D is DOUBLE PRECISION array, dimension (min(M,N))
00092 *>          The diagonal elements of the bidiagonal matrix B.
00093 *> \endverbatim
00094 *>
00095 *> \param[in] E
00096 *> \verbatim
00097 *>          E is DOUBLE PRECISION array, dimension (min(M,N)-1)
00098 *>          The superdiagonal elements of the bidiagonal matrix B if
00099 *>          m >= n, or the subdiagonal elements of B if m < n.
00100 *> \endverbatim
00101 *>
00102 *> \param[in] PT
00103 *> \verbatim
00104 *>          PT is DOUBLE PRECISION array, dimension (LDPT,N)
00105 *>          The min(m,n) by n orthogonal matrix P' in the reduction
00106 *>          A = Q * B * P'.
00107 *> \endverbatim
00108 *>
00109 *> \param[in] LDPT
00110 *> \verbatim
00111 *>          LDPT is INTEGER
00112 *>          The leading dimension of the array PT.
00113 *>          LDPT >= max(1,min(M,N)).
00114 *> \endverbatim
00115 *>
00116 *> \param[out] WORK
00117 *> \verbatim
00118 *>          WORK is DOUBLE PRECISION array, dimension (M+N)
00119 *> \endverbatim
00120 *>
00121 *> \param[out] RESID
00122 *> \verbatim
00123 *>          RESID is DOUBLE PRECISION
00124 *>          The test ratio:  norm(A - Q * B * P') / ( n * norm(A) * EPS )
00125 *> \endverbatim
00126 *
00127 *  Authors:
00128 *  ========
00129 *
00130 *> \author Univ. of Tennessee 
00131 *> \author Univ. of California Berkeley 
00132 *> \author Univ. of Colorado Denver 
00133 *> \author NAG Ltd. 
00134 *
00135 *> \date November 2011
00136 *
00137 *> \ingroup double_eig
00138 *
00139 *  =====================================================================
00140       SUBROUTINE DBDT01( M, N, KD, A, LDA, Q, LDQ, D, E, PT, LDPT, WORK,
00141      $                   RESID )
00142 *
00143 *  -- LAPACK test routine (version 3.4.0) --
00144 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00145 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00146 *     November 2011
00147 *
00148 *     .. Scalar Arguments ..
00149       INTEGER            KD, LDA, LDPT, LDQ, M, N
00150       DOUBLE PRECISION   RESID
00151 *     ..
00152 *     .. Array Arguments ..
00153       DOUBLE PRECISION   A( LDA, * ), D( * ), E( * ), PT( LDPT, * ),
00154      $                   Q( LDQ, * ), WORK( * )
00155 *     ..
00156 *
00157 *  =====================================================================
00158 *
00159 *     .. Parameters ..
00160       DOUBLE PRECISION   ZERO, ONE
00161       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00162 *     ..
00163 *     .. Local Scalars ..
00164       INTEGER            I, J
00165       DOUBLE PRECISION   ANORM, EPS
00166 *     ..
00167 *     .. External Functions ..
00168       DOUBLE PRECISION   DASUM, DLAMCH, DLANGE
00169       EXTERNAL           DASUM, DLAMCH, DLANGE
00170 *     ..
00171 *     .. External Subroutines ..
00172       EXTERNAL           DCOPY, DGEMV
00173 *     ..
00174 *     .. Intrinsic Functions ..
00175       INTRINSIC          DBLE, MAX, MIN
00176 *     ..
00177 *     .. Executable Statements ..
00178 *
00179 *     Quick return if possible
00180 *
00181       IF( M.LE.0 .OR. N.LE.0 ) THEN
00182          RESID = ZERO
00183          RETURN
00184       END IF
00185 *
00186 *     Compute A - Q * B * P' one column at a time.
00187 *
00188       RESID = ZERO
00189       IF( KD.NE.0 ) THEN
00190 *
00191 *        B is bidiagonal.
00192 *
00193          IF( KD.NE.0 .AND. M.GE.N ) THEN
00194 *
00195 *           B is upper bidiagonal and M >= N.
00196 *
00197             DO 20 J = 1, N
00198                CALL DCOPY( M, A( 1, J ), 1, WORK, 1 )
00199                DO 10 I = 1, N - 1
00200                   WORK( M+I ) = D( I )*PT( I, J ) + E( I )*PT( I+1, J )
00201    10          CONTINUE
00202                WORK( M+N ) = D( N )*PT( N, J )
00203                CALL DGEMV( 'No transpose', M, N, -ONE, Q, LDQ,
00204      $                     WORK( M+1 ), 1, ONE, WORK, 1 )
00205                RESID = MAX( RESID, DASUM( M, WORK, 1 ) )
00206    20       CONTINUE
00207          ELSE IF( KD.LT.0 ) THEN
00208 *
00209 *           B is upper bidiagonal and M < N.
00210 *
00211             DO 40 J = 1, N
00212                CALL DCOPY( M, A( 1, J ), 1, WORK, 1 )
00213                DO 30 I = 1, M - 1
00214                   WORK( M+I ) = D( I )*PT( I, J ) + E( I )*PT( I+1, J )
00215    30          CONTINUE
00216                WORK( M+M ) = D( M )*PT( M, J )
00217                CALL DGEMV( 'No transpose', M, M, -ONE, Q, LDQ,
00218      $                     WORK( M+1 ), 1, ONE, WORK, 1 )
00219                RESID = MAX( RESID, DASUM( M, WORK, 1 ) )
00220    40       CONTINUE
00221          ELSE
00222 *
00223 *           B is lower bidiagonal.
00224 *
00225             DO 60 J = 1, N
00226                CALL DCOPY( M, A( 1, J ), 1, WORK, 1 )
00227                WORK( M+1 ) = D( 1 )*PT( 1, J )
00228                DO 50 I = 2, M
00229                   WORK( M+I ) = E( I-1 )*PT( I-1, J ) +
00230      $                          D( I )*PT( I, J )
00231    50          CONTINUE
00232                CALL DGEMV( 'No transpose', M, M, -ONE, Q, LDQ,
00233      $                     WORK( M+1 ), 1, ONE, WORK, 1 )
00234                RESID = MAX( RESID, DASUM( M, WORK, 1 ) )
00235    60       CONTINUE
00236          END IF
00237       ELSE
00238 *
00239 *        B is diagonal.
00240 *
00241          IF( M.GE.N ) THEN
00242             DO 80 J = 1, N
00243                CALL DCOPY( M, A( 1, J ), 1, WORK, 1 )
00244                DO 70 I = 1, N
00245                   WORK( M+I ) = D( I )*PT( I, J )
00246    70          CONTINUE
00247                CALL DGEMV( 'No transpose', M, N, -ONE, Q, LDQ,
00248      $                     WORK( M+1 ), 1, ONE, WORK, 1 )
00249                RESID = MAX( RESID, DASUM( M, WORK, 1 ) )
00250    80       CONTINUE
00251          ELSE
00252             DO 100 J = 1, N
00253                CALL DCOPY( M, A( 1, J ), 1, WORK, 1 )
00254                DO 90 I = 1, M
00255                   WORK( M+I ) = D( I )*PT( I, J )
00256    90          CONTINUE
00257                CALL DGEMV( 'No transpose', M, M, -ONE, Q, LDQ,
00258      $                     WORK( M+1 ), 1, ONE, WORK, 1 )
00259                RESID = MAX( RESID, DASUM( M, WORK, 1 ) )
00260   100       CONTINUE
00261          END IF
00262       END IF
00263 *
00264 *     Compute norm(A - Q * B * P') / ( n * norm(A) * EPS )
00265 *
00266       ANORM = DLANGE( '1', M, N, A, LDA, WORK )
00267       EPS = DLAMCH( 'Precision' )
00268 *
00269       IF( ANORM.LE.ZERO ) THEN
00270          IF( RESID.NE.ZERO )
00271      $      RESID = ONE / EPS
00272       ELSE
00273          IF( ANORM.GE.RESID ) THEN
00274             RESID = ( RESID / ANORM ) / ( DBLE( N )*EPS )
00275          ELSE
00276             IF( ANORM.LT.ONE ) THEN
00277                RESID = ( MIN( RESID, DBLE( N )*ANORM ) / ANORM ) /
00278      $                 ( DBLE( N )*EPS )
00279             ELSE
00280                RESID = MIN( RESID / ANORM, DBLE( N ) ) /
00281      $                 ( DBLE( N )*EPS )
00282             END IF
00283          END IF
00284       END IF
00285 *
00286       RETURN
00287 *
00288 *     End of DBDT01
00289 *
00290       END
 All Files Functions