1 SUBROUTINE DBDT02( M, N, B, LDB, C, LDC, U, LDU, WORK, RESID )
2 *
3 * -- LAPACK test routine (version 3.1) --
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
5 * November 2006
6 *
7 * .. Scalar Arguments ..
8 INTEGER LDB, LDC, LDU, M, N
9 DOUBLE PRECISION RESID
10 * ..
11 * .. Array Arguments ..
12 DOUBLE PRECISION B( LDB, * ), C( LDC, * ), U( LDU, * ),
13 $ WORK( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DBDT02 tests the change of basis C = U' * B by computing the residual
20 *
21 * RESID = norm( B - U * C ) / ( max(m,n) * norm(B) * EPS ),
22 *
23 * where B and C are M by N matrices, U is an M by M orthogonal matrix,
24 * and EPS is the machine precision.
25 *
26 * Arguments
27 * =========
28 *
29 * M (input) INTEGER
30 * The number of rows of the matrices B and C and the order of
31 * the matrix Q.
32 *
33 * N (input) INTEGER
34 * The number of columns of the matrices B and C.
35 *
36 * B (input) DOUBLE PRECISION array, dimension (LDB,N)
37 * The m by n matrix B.
38 *
39 * LDB (input) INTEGER
40 * The leading dimension of the array B. LDB >= max(1,M).
41 *
42 * C (input) DOUBLE PRECISION array, dimension (LDC,N)
43 * The m by n matrix C, assumed to contain U' * B.
44 *
45 * LDC (input) INTEGER
46 * The leading dimension of the array C. LDC >= max(1,M).
47 *
48 * U (input) DOUBLE PRECISION array, dimension (LDU,M)
49 * The m by m orthogonal matrix U.
50 *
51 * LDU (input) INTEGER
52 * The leading dimension of the array U. LDU >= max(1,M).
53 *
54 * WORK (workspace) DOUBLE PRECISION array, dimension (M)
55 *
56 * RESID (output) DOUBLE PRECISION
57 * RESID = norm( B - U * C ) / ( max(m,n) * norm(B) * EPS ),
58 *
59 * ======================================================================
60 *
61 * .. Parameters ..
62 DOUBLE PRECISION ZERO, ONE
63 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
64 * ..
65 * .. Local Scalars ..
66 INTEGER J
67 DOUBLE PRECISION BNORM, EPS, REALMN
68 * ..
69 * .. External Functions ..
70 DOUBLE PRECISION DASUM, DLAMCH, DLANGE
71 EXTERNAL DASUM, DLAMCH, DLANGE
72 * ..
73 * .. External Subroutines ..
74 EXTERNAL DCOPY, DGEMV
75 * ..
76 * .. Intrinsic Functions ..
77 INTRINSIC DBLE, MAX, MIN
78 * ..
79 * .. Executable Statements ..
80 *
81 * Quick return if possible
82 *
83 RESID = ZERO
84 IF( M.LE.0 .OR. N.LE.0 )
85 $ RETURN
86 REALMN = DBLE( MAX( M, N ) )
87 EPS = DLAMCH( 'Precision' )
88 *
89 * Compute norm( B - U * C )
90 *
91 DO 10 J = 1, N
92 CALL DCOPY( M, B( 1, J ), 1, WORK, 1 )
93 CALL DGEMV( 'No transpose', M, M, -ONE, U, LDU, C( 1, J ), 1,
94 $ ONE, WORK, 1 )
95 RESID = MAX( RESID, DASUM( M, WORK, 1 ) )
96 10 CONTINUE
97 *
98 * Compute norm of B.
99 *
100 BNORM = DLANGE( '1', M, N, B, LDB, WORK )
101 *
102 IF( BNORM.LE.ZERO ) THEN
103 IF( RESID.NE.ZERO )
104 $ RESID = ONE / EPS
105 ELSE
106 IF( BNORM.GE.RESID ) THEN
107 RESID = ( RESID / BNORM ) / ( REALMN*EPS )
108 ELSE
109 IF( BNORM.LT.ONE ) THEN
110 RESID = ( MIN( RESID, REALMN*BNORM ) / BNORM ) /
111 $ ( REALMN*EPS )
112 ELSE
113 RESID = MIN( RESID / BNORM, REALMN ) / ( REALMN*EPS )
114 END IF
115 END IF
116 END IF
117 RETURN
118 *
119 * End of DBDT02
120 *
121 END
2 *
3 * -- LAPACK test routine (version 3.1) --
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
5 * November 2006
6 *
7 * .. Scalar Arguments ..
8 INTEGER LDB, LDC, LDU, M, N
9 DOUBLE PRECISION RESID
10 * ..
11 * .. Array Arguments ..
12 DOUBLE PRECISION B( LDB, * ), C( LDC, * ), U( LDU, * ),
13 $ WORK( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DBDT02 tests the change of basis C = U' * B by computing the residual
20 *
21 * RESID = norm( B - U * C ) / ( max(m,n) * norm(B) * EPS ),
22 *
23 * where B and C are M by N matrices, U is an M by M orthogonal matrix,
24 * and EPS is the machine precision.
25 *
26 * Arguments
27 * =========
28 *
29 * M (input) INTEGER
30 * The number of rows of the matrices B and C and the order of
31 * the matrix Q.
32 *
33 * N (input) INTEGER
34 * The number of columns of the matrices B and C.
35 *
36 * B (input) DOUBLE PRECISION array, dimension (LDB,N)
37 * The m by n matrix B.
38 *
39 * LDB (input) INTEGER
40 * The leading dimension of the array B. LDB >= max(1,M).
41 *
42 * C (input) DOUBLE PRECISION array, dimension (LDC,N)
43 * The m by n matrix C, assumed to contain U' * B.
44 *
45 * LDC (input) INTEGER
46 * The leading dimension of the array C. LDC >= max(1,M).
47 *
48 * U (input) DOUBLE PRECISION array, dimension (LDU,M)
49 * The m by m orthogonal matrix U.
50 *
51 * LDU (input) INTEGER
52 * The leading dimension of the array U. LDU >= max(1,M).
53 *
54 * WORK (workspace) DOUBLE PRECISION array, dimension (M)
55 *
56 * RESID (output) DOUBLE PRECISION
57 * RESID = norm( B - U * C ) / ( max(m,n) * norm(B) * EPS ),
58 *
59 * ======================================================================
60 *
61 * .. Parameters ..
62 DOUBLE PRECISION ZERO, ONE
63 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
64 * ..
65 * .. Local Scalars ..
66 INTEGER J
67 DOUBLE PRECISION BNORM, EPS, REALMN
68 * ..
69 * .. External Functions ..
70 DOUBLE PRECISION DASUM, DLAMCH, DLANGE
71 EXTERNAL DASUM, DLAMCH, DLANGE
72 * ..
73 * .. External Subroutines ..
74 EXTERNAL DCOPY, DGEMV
75 * ..
76 * .. Intrinsic Functions ..
77 INTRINSIC DBLE, MAX, MIN
78 * ..
79 * .. Executable Statements ..
80 *
81 * Quick return if possible
82 *
83 RESID = ZERO
84 IF( M.LE.0 .OR. N.LE.0 )
85 $ RETURN
86 REALMN = DBLE( MAX( M, N ) )
87 EPS = DLAMCH( 'Precision' )
88 *
89 * Compute norm( B - U * C )
90 *
91 DO 10 J = 1, N
92 CALL DCOPY( M, B( 1, J ), 1, WORK, 1 )
93 CALL DGEMV( 'No transpose', M, M, -ONE, U, LDU, C( 1, J ), 1,
94 $ ONE, WORK, 1 )
95 RESID = MAX( RESID, DASUM( M, WORK, 1 ) )
96 10 CONTINUE
97 *
98 * Compute norm of B.
99 *
100 BNORM = DLANGE( '1', M, N, B, LDB, WORK )
101 *
102 IF( BNORM.LE.ZERO ) THEN
103 IF( RESID.NE.ZERO )
104 $ RESID = ONE / EPS
105 ELSE
106 IF( BNORM.GE.RESID ) THEN
107 RESID = ( RESID / BNORM ) / ( REALMN*EPS )
108 ELSE
109 IF( BNORM.LT.ONE ) THEN
110 RESID = ( MIN( RESID, REALMN*BNORM ) / BNORM ) /
111 $ ( REALMN*EPS )
112 ELSE
113 RESID = MIN( RESID / BNORM, REALMN ) / ( REALMN*EPS )
114 END IF
115 END IF
116 END IF
117 RETURN
118 *
119 * End of DBDT02
120 *
121 END