80 IMPLICIT NONE
81
82
83
84
85
86
87 INTEGER LWORK, M, N, L, NB, LDT
88
89 DOUBLE PRECISION RESULT(6)
90
91
92
93
94
95 COMPLEX*16, ALLOCATABLE :: AF(:,:), Q(:,:),
96 $ R(:,:), WORK( : ), T(:,:),
97 $ CF(:,:), DF(:,:), A(:,:), C(:,:), D(:,:)
98 DOUBLE PRECISION, ALLOCATABLE :: RWORK(:)
99
100
101 DOUBLE PRECISION ZERO
102 COMPLEX*16 ONE, CZERO
103 parameter( zero = 0.0, one = (1.0,0.0), czero=(0.0,0.0) )
104
105
106 INTEGER INFO, J, K, M2, NP1
107 DOUBLE PRECISION ANORM, EPS, RESID, CNORM, DNORM
108
109
110 INTEGER ISEED( 4 )
111
112
113 DOUBLE PRECISION DLAMCH
114 DOUBLE PRECISION ZLANGE, ZLANSY
115 LOGICAL LSAME
117
118
119 DATA iseed / 1988, 1989, 1990, 1991 /
120
122 k = n
123 m2 = m+n
124 IF( m.GT.0 ) THEN
125 np1 = n+1
126 ELSE
127 np1 = 1
128 END IF
129 lwork = m2*m2*nb
130
131
132
133 ALLOCATE(a(m2,n),af(m2,n),q(m2,m2),r(m2,m2),rwork(m2),
134 $ work(lwork),t(nb,n),c(m2,n),cf(m2,n),
135 $ d(n,m2),df(n,m2) )
136
137
138
139 ldt=nb
140 CALL zlaset(
'Full', m2, n, czero, czero, a, m2 )
141 CALL zlaset(
'Full', nb, n, czero, czero, t, nb )
142 DO j=1,n
143 CALL zlarnv( 2, iseed, j, a( 1, j ) )
144 END DO
145 IF( m.GT.0 ) THEN
146 DO j=1,n
147 CALL zlarnv( 2, iseed, m-l, a( min(n+m,n+1), j ) )
148 END DO
149 END IF
150 IF( l.GT.0 ) THEN
151 DO j=1,n
152 CALL zlarnv( 2, iseed, min(j,l), a( min(n+m,n+m-l+1), j ) )
153 END DO
154 END IF
155
156
157
158 CALL zlacpy(
'Full', m2, n, a, m2, af, m2 )
159
160
161
162 CALL ztpqrt( m,n,l,nb,af,m2,af(np1,1),m2,t,ldt,work,info)
163
164
165
166 CALL zlaset(
'Full', m2, m2, czero, one, q, m2 )
167 CALL zgemqrt(
'R',
'N', m2, m2, k, nb, af, m2, t, ldt, q, m2,
168 $ work, info )
169
170
171
172 CALL zlaset(
'Full', m2, n, czero, czero, r, m2 )
173 CALL zlacpy(
'Upper', m2, n, af, m2, r, m2 )
174
175
176
177 CALL zgemm(
'C',
'N', m2, n, m2, -one, q, m2, a, m2, one, r, m2 )
178 anorm =
zlange(
'1', m2, n, a, m2, rwork )
179 resid =
zlange(
'1', m2, n, r, m2, rwork )
180 IF( anorm.GT.zero ) THEN
181 result( 1 ) = resid / (eps*anorm*max(1,m2))
182 ELSE
183 result( 1 ) = zero
184 END IF
185
186
187
188 CALL zlaset(
'Full', m2, m2, czero, one, r, m2 )
189 CALL zherk(
'U',
'C', m2, m2, dreal(-one), q, m2, dreal(one),
190 $ r, m2 )
191 resid =
zlansy(
'1',
'Upper', m2, r, m2, rwork )
192 result( 2 ) = resid / (eps*max(1,m2))
193
194
195
196 DO j=1,n
197 CALL zlarnv( 2, iseed, m2, c( 1, j ) )
198 END DO
199 cnorm =
zlange(
'1', m2, n, c, m2, rwork)
200 CALL zlacpy(
'Full', m2, n, c, m2, cf, m2 )
201
202
203
204 CALL ztpmqrt(
'L',
'N', m,n,k,l,nb,af(np1,1),m2,t,ldt,cf,m2,
205 $ cf(np1,1),m2,work,info)
206
207
208
209 CALL zgemm(
'N',
'N', m2, n, m2, -one, q, m2, c, m2, one, cf, m2 )
210 resid =
zlange(
'1', m2, n, cf, m2, rwork )
211 IF( cnorm.GT.zero ) THEN
212 result( 3 ) = resid / (eps*max(1,m2)*cnorm)
213 ELSE
214 result( 3 ) = zero
215 END IF
216
217
218
219 CALL zlacpy(
'Full', m2, n, c, m2, cf, m2 )
220
221
222
223 CALL ztpmqrt(
'L',
'C',m,n,k,l,nb,af(np1,1),m2,t,ldt,cf,m2,
224 $ cf(np1,1),m2,work,info)
225
226
227
228 CALL zgemm(
'C',
'N',m2,n,m2,-one,q,m2,c,m2,one,cf,m2)
229 resid =
zlange(
'1', m2, n, cf, m2, rwork )
230 IF( cnorm.GT.zero ) THEN
231 result( 4 ) = resid / (eps*max(1,m2)*cnorm)
232 ELSE
233 result( 4 ) = zero
234 END IF
235
236
237
238 DO j=1,m2
239 CALL zlarnv( 2, iseed, n, d( 1, j ) )
240 END DO
241 dnorm =
zlange(
'1', n, m2, d, n, rwork)
242 CALL zlacpy(
'Full', n, m2, d, n, df, n )
243
244
245
246 CALL ztpmqrt(
'R',
'N',n,m,n,l,nb,af(np1,1),m2,t,ldt,df,n,
247 $ df(1,np1),n,work,info)
248
249
250
251 CALL zgemm(
'N',
'N',n,m2,m2,-one,d,n,q,m2,one,df,n)
252 resid =
zlange(
'1',n, m2,df,n,rwork )
253 IF( cnorm.GT.zero ) THEN
254 result( 5 ) = resid / (eps*max(1,m2)*dnorm)
255 ELSE
256 result( 5 ) = zero
257 END IF
258
259
260
261 CALL zlacpy(
'Full',n,m2,d,n,df,n )
262
263
264
265 CALL ztpmqrt(
'R',
'C',n,m,n,l,nb,af(np1,1),m2,t,ldt,df,n,
266 $ df(1,np1),n,work,info)
267
268
269
270
271 CALL zgemm(
'N',
'C', n, m2, m2, -one, d, n, q, m2, one, df, n )
272 resid =
zlange(
'1', n, m2, df, n, rwork )
273 IF( cnorm.GT.zero ) THEN
274 result( 6 ) = resid / (eps*max(1,m2)*dnorm)
275 ELSE
276 result( 6 ) = zero
277 END IF
278
279
280
281 DEALLOCATE ( a, af, q, r, rwork, work, t, c, d, cf, df)
282 RETURN
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
subroutine zgemqrt(side, trans, m, n, k, nb, v, ldv, t, ldt, c, ldc, work, info)
ZGEMQRT
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
subroutine ztpmqrt(side, trans, m, n, k, l, nb, v, ldv, t, ldt, a, lda, b, ldb, work, info)
ZTPMQRT
subroutine ztpqrt(m, n, l, nb, a, lda, b, ldb, t, ldt, work, info)
ZTPQRT