1 SUBROUTINE DORT01( ROWCOL, M, N, U, LDU, WORK, LWORK, 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 CHARACTER ROWCOL
9 INTEGER LDU, LWORK, M, N
10 DOUBLE PRECISION RESID
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION U( LDU, * ), WORK( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DORT01 checks that the matrix U is orthogonal by computing the ratio
20 *
21 * RESID = norm( I - U*U' ) / ( n * EPS ), if ROWCOL = 'R',
22 * or
23 * RESID = norm( I - U'*U ) / ( m * EPS ), if ROWCOL = 'C'.
24 *
25 * Alternatively, if there isn't sufficient workspace to form
26 * I - U*U' or I - U'*U, the ratio is computed as
27 *
28 * RESID = abs( I - U*U' ) / ( n * EPS ), if ROWCOL = 'R',
29 * or
30 * RESID = abs( I - U'*U ) / ( m * EPS ), if ROWCOL = 'C'.
31 *
32 * where EPS is the machine precision. ROWCOL is used only if m = n;
33 * if m > n, ROWCOL is assumed to be 'C', and if m < n, ROWCOL is
34 * assumed to be 'R'.
35 *
36 * Arguments
37 * =========
38 *
39 * ROWCOL (input) CHARACTER
40 * Specifies whether the rows or columns of U should be checked
41 * for orthogonality. Used only if M = N.
42 * = 'R': Check for orthogonal rows of U
43 * = 'C': Check for orthogonal columns of U
44 *
45 * M (input) INTEGER
46 * The number of rows of the matrix U.
47 *
48 * N (input) INTEGER
49 * The number of columns of the matrix U.
50 *
51 * U (input) DOUBLE PRECISION array, dimension (LDU,N)
52 * The orthogonal matrix U. U is checked for orthogonal columns
53 * if m > n or if m = n and ROWCOL = 'C'. U is checked for
54 * orthogonal rows if m < n or if m = n and ROWCOL = 'R'.
55 *
56 * LDU (input) INTEGER
57 * The leading dimension of the array U. LDU >= max(1,M).
58 *
59 * WORK (workspace) DOUBLE PRECISION array, dimension (LWORK)
60 *
61 * LWORK (input) INTEGER
62 * The length of the array WORK. For best performance, LWORK
63 * should be at least N*(N+1) if ROWCOL = 'C' or M*(M+1) if
64 * ROWCOL = 'R', but the test will be done even if LWORK is 0.
65 *
66 * RESID (output) DOUBLE PRECISION
67 * RESID = norm( I - U * U' ) / ( n * EPS ), if ROWCOL = 'R', or
68 * RESID = norm( I - U' * U ) / ( m * EPS ), if ROWCOL = 'C'.
69 *
70 * =====================================================================
71 *
72 * .. Parameters ..
73 DOUBLE PRECISION ZERO, ONE
74 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
75 * ..
76 * .. Local Scalars ..
77 CHARACTER TRANSU
78 INTEGER I, J, K, LDWORK, MNMIN
79 DOUBLE PRECISION EPS, TMP
80 * ..
81 * .. External Functions ..
82 LOGICAL LSAME
83 DOUBLE PRECISION DDOT, DLAMCH, DLANSY
84 EXTERNAL LSAME, DDOT, DLAMCH, DLANSY
85 * ..
86 * .. External Subroutines ..
87 EXTERNAL DLASET, DSYRK
88 * ..
89 * .. Intrinsic Functions ..
90 INTRINSIC ABS, DBLE, MAX, MIN
91 * ..
92 * .. Executable Statements ..
93 *
94 RESID = ZERO
95 *
96 * Quick return if possible
97 *
98 IF( M.LE.0 .OR. N.LE.0 )
99 $ RETURN
100 *
101 EPS = DLAMCH( 'Precision' )
102 IF( M.LT.N .OR. ( M.EQ.N .AND. LSAME( ROWCOL, 'R' ) ) ) THEN
103 TRANSU = 'N'
104 K = N
105 ELSE
106 TRANSU = 'T'
107 K = M
108 END IF
109 MNMIN = MIN( M, N )
110 *
111 IF( ( MNMIN+1 )*MNMIN.LE.LWORK ) THEN
112 LDWORK = MNMIN
113 ELSE
114 LDWORK = 0
115 END IF
116 IF( LDWORK.GT.0 ) THEN
117 *
118 * Compute I - U*U' or I - U'*U.
119 *
120 CALL DLASET( 'Upper', MNMIN, MNMIN, ZERO, ONE, WORK, LDWORK )
121 CALL DSYRK( 'Upper', TRANSU, MNMIN, K, -ONE, U, LDU, ONE, WORK,
122 $ LDWORK )
123 *
124 * Compute norm( I - U*U' ) / ( K * EPS ) .
125 *
126 RESID = DLANSY( '1', 'Upper', MNMIN, WORK, LDWORK,
127 $ WORK( LDWORK*MNMIN+1 ) )
128 RESID = ( RESID / DBLE( K ) ) / EPS
129 ELSE IF( TRANSU.EQ.'T' ) THEN
130 *
131 * Find the maximum element in abs( I - U'*U ) / ( m * EPS )
132 *
133 DO 20 J = 1, N
134 DO 10 I = 1, J
135 IF( I.NE.J ) THEN
136 TMP = ZERO
137 ELSE
138 TMP = ONE
139 END IF
140 TMP = TMP - DDOT( M, U( 1, I ), 1, U( 1, J ), 1 )
141 RESID = MAX( RESID, ABS( TMP ) )
142 10 CONTINUE
143 20 CONTINUE
144 RESID = ( RESID / DBLE( M ) ) / EPS
145 ELSE
146 *
147 * Find the maximum element in abs( I - U*U' ) / ( n * EPS )
148 *
149 DO 40 J = 1, M
150 DO 30 I = 1, J
151 IF( I.NE.J ) THEN
152 TMP = ZERO
153 ELSE
154 TMP = ONE
155 END IF
156 TMP = TMP - DDOT( N, U( J, 1 ), LDU, U( I, 1 ), LDU )
157 RESID = MAX( RESID, ABS( TMP ) )
158 30 CONTINUE
159 40 CONTINUE
160 RESID = ( RESID / DBLE( N ) ) / EPS
161 END IF
162 RETURN
163 *
164 * End of DORT01
165 *
166 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 CHARACTER ROWCOL
9 INTEGER LDU, LWORK, M, N
10 DOUBLE PRECISION RESID
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION U( LDU, * ), WORK( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DORT01 checks that the matrix U is orthogonal by computing the ratio
20 *
21 * RESID = norm( I - U*U' ) / ( n * EPS ), if ROWCOL = 'R',
22 * or
23 * RESID = norm( I - U'*U ) / ( m * EPS ), if ROWCOL = 'C'.
24 *
25 * Alternatively, if there isn't sufficient workspace to form
26 * I - U*U' or I - U'*U, the ratio is computed as
27 *
28 * RESID = abs( I - U*U' ) / ( n * EPS ), if ROWCOL = 'R',
29 * or
30 * RESID = abs( I - U'*U ) / ( m * EPS ), if ROWCOL = 'C'.
31 *
32 * where EPS is the machine precision. ROWCOL is used only if m = n;
33 * if m > n, ROWCOL is assumed to be 'C', and if m < n, ROWCOL is
34 * assumed to be 'R'.
35 *
36 * Arguments
37 * =========
38 *
39 * ROWCOL (input) CHARACTER
40 * Specifies whether the rows or columns of U should be checked
41 * for orthogonality. Used only if M = N.
42 * = 'R': Check for orthogonal rows of U
43 * = 'C': Check for orthogonal columns of U
44 *
45 * M (input) INTEGER
46 * The number of rows of the matrix U.
47 *
48 * N (input) INTEGER
49 * The number of columns of the matrix U.
50 *
51 * U (input) DOUBLE PRECISION array, dimension (LDU,N)
52 * The orthogonal matrix U. U is checked for orthogonal columns
53 * if m > n or if m = n and ROWCOL = 'C'. U is checked for
54 * orthogonal rows if m < n or if m = n and ROWCOL = 'R'.
55 *
56 * LDU (input) INTEGER
57 * The leading dimension of the array U. LDU >= max(1,M).
58 *
59 * WORK (workspace) DOUBLE PRECISION array, dimension (LWORK)
60 *
61 * LWORK (input) INTEGER
62 * The length of the array WORK. For best performance, LWORK
63 * should be at least N*(N+1) if ROWCOL = 'C' or M*(M+1) if
64 * ROWCOL = 'R', but the test will be done even if LWORK is 0.
65 *
66 * RESID (output) DOUBLE PRECISION
67 * RESID = norm( I - U * U' ) / ( n * EPS ), if ROWCOL = 'R', or
68 * RESID = norm( I - U' * U ) / ( m * EPS ), if ROWCOL = 'C'.
69 *
70 * =====================================================================
71 *
72 * .. Parameters ..
73 DOUBLE PRECISION ZERO, ONE
74 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
75 * ..
76 * .. Local Scalars ..
77 CHARACTER TRANSU
78 INTEGER I, J, K, LDWORK, MNMIN
79 DOUBLE PRECISION EPS, TMP
80 * ..
81 * .. External Functions ..
82 LOGICAL LSAME
83 DOUBLE PRECISION DDOT, DLAMCH, DLANSY
84 EXTERNAL LSAME, DDOT, DLAMCH, DLANSY
85 * ..
86 * .. External Subroutines ..
87 EXTERNAL DLASET, DSYRK
88 * ..
89 * .. Intrinsic Functions ..
90 INTRINSIC ABS, DBLE, MAX, MIN
91 * ..
92 * .. Executable Statements ..
93 *
94 RESID = ZERO
95 *
96 * Quick return if possible
97 *
98 IF( M.LE.0 .OR. N.LE.0 )
99 $ RETURN
100 *
101 EPS = DLAMCH( 'Precision' )
102 IF( M.LT.N .OR. ( M.EQ.N .AND. LSAME( ROWCOL, 'R' ) ) ) THEN
103 TRANSU = 'N'
104 K = N
105 ELSE
106 TRANSU = 'T'
107 K = M
108 END IF
109 MNMIN = MIN( M, N )
110 *
111 IF( ( MNMIN+1 )*MNMIN.LE.LWORK ) THEN
112 LDWORK = MNMIN
113 ELSE
114 LDWORK = 0
115 END IF
116 IF( LDWORK.GT.0 ) THEN
117 *
118 * Compute I - U*U' or I - U'*U.
119 *
120 CALL DLASET( 'Upper', MNMIN, MNMIN, ZERO, ONE, WORK, LDWORK )
121 CALL DSYRK( 'Upper', TRANSU, MNMIN, K, -ONE, U, LDU, ONE, WORK,
122 $ LDWORK )
123 *
124 * Compute norm( I - U*U' ) / ( K * EPS ) .
125 *
126 RESID = DLANSY( '1', 'Upper', MNMIN, WORK, LDWORK,
127 $ WORK( LDWORK*MNMIN+1 ) )
128 RESID = ( RESID / DBLE( K ) ) / EPS
129 ELSE IF( TRANSU.EQ.'T' ) THEN
130 *
131 * Find the maximum element in abs( I - U'*U ) / ( m * EPS )
132 *
133 DO 20 J = 1, N
134 DO 10 I = 1, J
135 IF( I.NE.J ) THEN
136 TMP = ZERO
137 ELSE
138 TMP = ONE
139 END IF
140 TMP = TMP - DDOT( M, U( 1, I ), 1, U( 1, J ), 1 )
141 RESID = MAX( RESID, ABS( TMP ) )
142 10 CONTINUE
143 20 CONTINUE
144 RESID = ( RESID / DBLE( M ) ) / EPS
145 ELSE
146 *
147 * Find the maximum element in abs( I - U*U' ) / ( n * EPS )
148 *
149 DO 40 J = 1, M
150 DO 30 I = 1, J
151 IF( I.NE.J ) THEN
152 TMP = ZERO
153 ELSE
154 TMP = ONE
155 END IF
156 TMP = TMP - DDOT( N, U( J, 1 ), LDU, U( I, 1 ), LDU )
157 RESID = MAX( RESID, ABS( TMP ) )
158 30 CONTINUE
159 40 CONTINUE
160 RESID = ( RESID / DBLE( N ) ) / EPS
161 END IF
162 RETURN
163 *
164 * End of DORT01
165 *
166 END