|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE ZBDT02( M, N, B, LDB, C, LDC, U, LDU, WORK, RWORK, 00002 $ RESID ) 00003 * 00004 * -- LAPACK test routine (version 3.1) -- 00005 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00006 * November 2006 00007 * 00008 * .. Scalar Arguments .. 00009 INTEGER LDB, LDC, LDU, M, N 00010 DOUBLE PRECISION RESID 00011 * .. 00012 * .. Array Arguments .. 00013 DOUBLE PRECISION RWORK( * ) 00014 COMPLEX*16 B( LDB, * ), C( LDC, * ), U( LDU, * ), 00015 $ WORK( * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * ZBDT02 tests the change of basis C = U' * B by computing the residual 00022 * 00023 * RESID = norm( B - U * C ) / ( max(m,n) * norm(B) * EPS ), 00024 * 00025 * where B and C are M by N matrices, U is an M by M orthogonal matrix, 00026 * and EPS is the machine precision. 00027 * 00028 * Arguments 00029 * ========= 00030 * 00031 * M (input) INTEGER 00032 * The number of rows of the matrices B and C and the order of 00033 * the matrix Q. 00034 * 00035 * N (input) INTEGER 00036 * The number of columns of the matrices B and C. 00037 * 00038 * B (input) COMPLEX*16 array, dimension (LDB,N) 00039 * The m by n matrix B. 00040 * 00041 * LDB (input) INTEGER 00042 * The leading dimension of the array B. LDB >= max(1,M). 00043 * 00044 * C (input) COMPLEX*16 array, dimension (LDC,N) 00045 * The m by n matrix C, assumed to contain U' * B. 00046 * 00047 * LDC (input) INTEGER 00048 * The leading dimension of the array C. LDC >= max(1,M). 00049 * 00050 * U (input) COMPLEX*16 array, dimension (LDU,M) 00051 * The m by m orthogonal matrix U. 00052 * 00053 * LDU (input) INTEGER 00054 * The leading dimension of the array U. LDU >= max(1,M). 00055 * 00056 * WORK (workspace) COMPLEX*16 array, dimension (M) 00057 * 00058 * RWORK (workspace) DOUBLE PRECISION array, dimension (M) 00059 * 00060 * RESID (output) DOUBLE PRECISION 00061 * RESID = norm( B - U * C ) / ( max(m,n) * norm(B) * EPS ), 00062 * 00063 * ====================================================================== 00064 * 00065 * .. Parameters .. 00066 DOUBLE PRECISION ZERO, ONE 00067 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00068 * .. 00069 * .. Local Scalars .. 00070 INTEGER J 00071 DOUBLE PRECISION BNORM, EPS, REALMN 00072 * .. 00073 * .. External Functions .. 00074 DOUBLE PRECISION DLAMCH, DZASUM, ZLANGE 00075 EXTERNAL DLAMCH, DZASUM, ZLANGE 00076 * .. 00077 * .. External Subroutines .. 00078 EXTERNAL ZCOPY, ZGEMV 00079 * .. 00080 * .. Intrinsic Functions .. 00081 INTRINSIC DBLE, DCMPLX, MAX, MIN 00082 * .. 00083 * .. Executable Statements .. 00084 * 00085 * Quick return if possible 00086 * 00087 RESID = ZERO 00088 IF( M.LE.0 .OR. N.LE.0 ) 00089 $ RETURN 00090 REALMN = DBLE( MAX( M, N ) ) 00091 EPS = DLAMCH( 'Precision' ) 00092 * 00093 * Compute norm( B - U * C ) 00094 * 00095 DO 10 J = 1, N 00096 CALL ZCOPY( M, B( 1, J ), 1, WORK, 1 ) 00097 CALL ZGEMV( 'No transpose', M, M, -DCMPLX( ONE ), U, LDU, 00098 $ C( 1, J ), 1, DCMPLX( ONE ), WORK, 1 ) 00099 RESID = MAX( RESID, DZASUM( M, WORK, 1 ) ) 00100 10 CONTINUE 00101 * 00102 * Compute norm of B. 00103 * 00104 BNORM = ZLANGE( '1', M, N, B, LDB, RWORK ) 00105 * 00106 IF( BNORM.LE.ZERO ) THEN 00107 IF( RESID.NE.ZERO ) 00108 $ RESID = ONE / EPS 00109 ELSE 00110 IF( BNORM.GE.RESID ) THEN 00111 RESID = ( RESID / BNORM ) / ( REALMN*EPS ) 00112 ELSE 00113 IF( BNORM.LT.ONE ) THEN 00114 RESID = ( MIN( RESID, REALMN*BNORM ) / BNORM ) / 00115 $ ( REALMN*EPS ) 00116 ELSE 00117 RESID = MIN( RESID / BNORM, REALMN ) / ( REALMN*EPS ) 00118 END IF 00119 END IF 00120 END IF 00121 RETURN 00122 * 00123 * End of ZBDT02 00124 * 00125 END
1.7.3 
