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