![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
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