73 IMPLICIT NONE
74
75
76
77
78
79
80 INTEGER M, N, NB
81
82 DOUBLE PRECISION RESULT(6)
83
84
85
86
87
88 COMPLEX*16, ALLOCATABLE :: AF(:,:), Q(:,:),
89 $ L(:,:), WORK( : ), T(:,:),
90 $ CF(:,:), DF(:,:), A(:,:), C(:,:), D(:,:)
91 DOUBLE PRECISION, ALLOCATABLE :: RWORK(:)
92
93
94 DOUBLE PRECISION ZERO
95 COMPLEX*16 ONE, CZERO
96 parameter( zero = 0.0)
97 parameter( one = (1.0,0.0), czero=(0.0,0.0) )
98
99
100 INTEGER INFO, J, K, LL, LWORK, LDT
101 DOUBLE PRECISION ANORM, EPS, RESID, CNORM, DNORM
102
103
104 INTEGER ISEED( 4 )
105
106
107 DOUBLE PRECISION DLAMCH
108 DOUBLE PRECISION ZLANGE, ZLANSY
109 LOGICAL LSAME
111
112
113 INTRINSIC max, min
114
115
116 DATA iseed / 1988, 1989, 1990, 1991 /
117
119 k = min(m,n)
120 ll = max(m,n)
121 lwork = max(2,ll)*max(2,ll)*nb
122
123
124
125 ALLOCATE ( a(m,n), af(m,n), q(n,n), l(ll,n), rwork(ll),
126 $ work(lwork), t(nb,n), c(m,n), cf(m,n),
127 $ d(n,m), df(n,m) )
128
129
130
131 ldt=nb
132 DO j=1,n
133 CALL zlarnv( 2, iseed, m, a( 1, j ) )
134 END DO
135 CALL zlacpy(
'Full', m, n, a, m, af, m )
136
137
138
139 CALL zgelqt( m, n, nb, af, m, t, ldt, work, info )
140
141
142
143 CALL zlaset(
'Full', n, n, czero, one, q, n )
144 CALL zgemlqt(
'R',
'N', n, n, k, nb, af, m, t, ldt, q, n,
145 $ work, info )
146
147
148
149 CALL zlaset(
'Full', ll, n, czero, czero, l, ll )
150 CALL zlacpy(
'Lower', m, n, af, m, l, ll )
151
152
153
154 CALL zgemm(
'N',
'C', m, n, n, -one, a, m, q, n, one, l, ll )
155 anorm =
zlange(
'1', m, n, a, m, rwork )
156 resid =
zlange(
'1', m, n, l, ll, rwork )
157 IF( anorm.GT.zero ) THEN
158 result( 1 ) = resid / (eps*max(1,m)*anorm)
159 ELSE
160 result( 1 ) = zero
161 END IF
162
163
164
165 CALL zlaset(
'Full', n, n, czero, one, l, ll )
166 CALL zherk(
'U',
'C', n, n, dreal(-one), q, n, dreal(one), l, ll)
167 resid =
zlansy(
'1',
'Upper', n, l, ll, rwork )
168 result( 2 ) = resid / (eps*max(1,n))
169
170
171
172 DO j=1,m
173 CALL zlarnv( 2, iseed, n, d( 1, j ) )
174 END DO
175 dnorm =
zlange(
'1', n, m, d, n, rwork)
176 CALL zlacpy(
'Full', n, m, d, n, df, n )
177
178
179
180 CALL zgemlqt(
'L',
'N', n, m, k, nb, af, m, t, nb, df, n,
181 $ work, info)
182
183
184
185 CALL zgemm(
'N',
'N', n, m, n, -one, q, n, d, n, one, df, n )
186 resid =
zlange(
'1', n, m, df, n, rwork )
187 IF( dnorm.GT.zero ) THEN
188 result( 3 ) = resid / (eps*max(1,m)*dnorm)
189 ELSE
190 result( 3 ) = zero
191 END IF
192
193
194
195 CALL zlacpy(
'Full', n, m, d, n, df, n )
196
197
198
199 CALL zgemlqt(
'L',
'C', n, m, k, nb, af, m, t, nb, df, n,
200 $ work, info)
201
202
203
204 CALL zgemm(
'C',
'N', n, m, n, -one, q, n, d, n, one, df, n )
205 resid =
zlange(
'1', n, m, df, n, rwork )
206 IF( dnorm.GT.zero ) THEN
207 result( 4 ) = resid / (eps*max(1,m)*dnorm)
208 ELSE
209 result( 4 ) = zero
210 END IF
211
212
213
214 DO j=1,n
215 CALL zlarnv( 2, iseed, m, c( 1, j ) )
216 END DO
217 cnorm =
zlange(
'1', m, n, c, m, rwork)
218 CALL zlacpy(
'Full', m, n, c, m, cf, m )
219
220
221
222 CALL zgemlqt(
'R',
'N', m, n, k, nb, af, m, t, nb, cf, m,
223 $ work, info)
224
225
226
227 CALL zgemm(
'N',
'N', m, n, n, -one, c, m, q, n, one, cf, m )
228 resid =
zlange(
'1', n, m, df, n, rwork )
229 IF( cnorm.GT.zero ) THEN
230 result( 5 ) = resid / (eps*max(1,m)*dnorm)
231 ELSE
232 result( 5 ) = zero
233 END IF
234
235
236
237 CALL zlacpy(
'Full', m, n, c, m, cf, m )
238
239
240
241 CALL zgemlqt(
'R',
'C', m, n, k, nb, af, m, t, nb, cf, m,
242 $ work, info)
243
244
245
246 CALL zgemm(
'N',
'C', m, n, n, -one, c, m, q, n, one, cf, m )
247 resid =
zlange(
'1', m, n, cf, m, rwork )
248 IF( cnorm.GT.zero ) THEN
249 result( 6 ) = resid / (eps*max(1,m)*dnorm)
250 ELSE
251 result( 6 ) = zero
252 END IF
253
254
255
256 DEALLOCATE ( a, af, q, l, rwork, work, t, c, d, cf, df)
257
258 RETURN
subroutine zgelqt(m, n, mb, a, lda, t, ldt, work, info)
ZGELQT
subroutine zgemlqt(side, trans, m, n, k, mb, v, ldv, t, ldt, c, ldc, work, info)
ZGEMLQT
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
subroutine zherk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
ZHERK
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
double precision function dlamch(cmach)
DLAMCH
double precision function zlange(norm, m, n, a, lda, work)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
double precision function zlansy(norm, uplo, n, a, lda, work)
ZLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
subroutine zlarnv(idist, iseed, n, x)
ZLARNV returns a vector of random numbers from a uniform or normal distribution.
subroutine zlaset(uplo, m, n, alpha, beta, a, lda)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
logical function lsame(ca, cb)
LSAME