![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CBDT01 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 CBDT01( 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 * REAL RESID 00017 * .. 00018 * .. Array Arguments .. 00019 * REAL D( * ), E( * ), RWORK( * ) 00020 * COMPLEX A( LDA, * ), PT( LDPT, * ), Q( LDQ, * ), 00021 * $ WORK( * ) 00022 * .. 00023 * 00024 * 00025 *> \par Purpose: 00026 * ============= 00027 *> 00028 *> \verbatim 00029 *> 00030 *> CBDT01 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 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 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 REAL 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 REAL 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 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 array, dimension (M+N) 00120 *> \endverbatim 00121 *> 00122 *> \param[out] RWORK 00123 *> \verbatim 00124 *> RWORK is REAL array, dimension (M) 00125 *> \endverbatim 00126 *> 00127 *> \param[out] RESID 00128 *> \verbatim 00129 *> RESID is REAL 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 complex_eig 00144 * 00145 * ===================================================================== 00146 SUBROUTINE CBDT01( 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 REAL RESID 00157 * .. 00158 * .. Array Arguments .. 00159 REAL D( * ), E( * ), RWORK( * ) 00160 COMPLEX A( LDA, * ), PT( LDPT, * ), Q( LDQ, * ), 00161 $ WORK( * ) 00162 * .. 00163 * 00164 * ===================================================================== 00165 * 00166 * .. Parameters .. 00167 REAL ZERO, ONE 00168 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00169 * .. 00170 * .. Local Scalars .. 00171 INTEGER I, J 00172 REAL ANORM, EPS 00173 * .. 00174 * .. External Functions .. 00175 REAL CLANGE, SCASUM, SLAMCH 00176 EXTERNAL CLANGE, SCASUM, SLAMCH 00177 * .. 00178 * .. External Subroutines .. 00179 EXTERNAL CCOPY, CGEMV 00180 * .. 00181 * .. Intrinsic Functions .. 00182 INTRINSIC CMPLX, MAX, MIN, REAL 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 CCOPY( 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 CGEMV( 'No transpose', M, N, -CMPLX( ONE ), Q, LDQ, 00211 $ WORK( M+1 ), 1, CMPLX( ONE ), WORK, 1 ) 00212 RESID = MAX( RESID, SCASUM( 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 CCOPY( 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 CGEMV( 'No transpose', M, M, -CMPLX( ONE ), Q, LDQ, 00225 $ WORK( M+1 ), 1, CMPLX( ONE ), WORK, 1 ) 00226 RESID = MAX( RESID, SCASUM( M, WORK, 1 ) ) 00227 40 CONTINUE 00228 ELSE 00229 * 00230 * B is lower bidiagonal. 00231 * 00232 DO 60 J = 1, N 00233 CALL CCOPY( 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 CGEMV( 'No transpose', M, M, -CMPLX( ONE ), Q, LDQ, 00240 $ WORK( M+1 ), 1, CMPLX( ONE ), WORK, 1 ) 00241 RESID = MAX( RESID, SCASUM( 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 CCOPY( 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 CGEMV( 'No transpose', M, N, -CMPLX( ONE ), Q, LDQ, 00255 $ WORK( M+1 ), 1, CMPLX( ONE ), WORK, 1 ) 00256 RESID = MAX( RESID, SCASUM( M, WORK, 1 ) ) 00257 80 CONTINUE 00258 ELSE 00259 DO 100 J = 1, N 00260 CALL CCOPY( 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 CGEMV( 'No transpose', M, M, -CMPLX( ONE ), Q, LDQ, 00265 $ WORK( M+1 ), 1, CMPLX( ONE ), WORK, 1 ) 00266 RESID = MAX( RESID, SCASUM( 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 = CLANGE( '1', M, N, A, LDA, RWORK ) 00274 EPS = SLAMCH( '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 ) / ( REAL( N )*EPS ) 00282 ELSE 00283 IF( ANORM.LT.ONE ) THEN 00284 RESID = ( MIN( RESID, REAL( N )*ANORM ) / ANORM ) / 00285 $ ( REAL( N )*EPS ) 00286 ELSE 00287 RESID = MIN( RESID / ANORM, REAL( N ) ) / 00288 $ ( REAL( N )*EPS ) 00289 END IF 00290 END IF 00291 END IF 00292 * 00293 RETURN 00294 * 00295 * End of CBDT01 00296 * 00297 END