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