90 INTEGER LWORK, M, N, L, NB, LDT
92 DOUBLE PRECISION RESULT(6)
98 COMPLEX*16,
ALLOCATABLE :: AF(:,:), Q(:,:),
99 $ R(:,:), RWORK(:), WORK( : ), T(:,:),
100 $ CF(:,:), DF(:,:), A(:,:), C(:,:), D(:,:)
103 DOUBLE PRECISION ZERO
104 COMPLEX*16 ONE, CZERO
105 parameter( zero = 0.0, one = (1.0,0.0), czero=(0.0,0.0) )
108 INTEGER INFO, J, K, N2, NP1,i
109 DOUBLE PRECISION ANORM, EPS, RESID, CNORM, DNORM
115 DOUBLE PRECISION DLAMCH
116 DOUBLE PRECISION ZLANGE, ZLANSY
121 DATA iseed / 1988, 1989, 1990, 1991 /
135 ALLOCATE(a(m,n2),af(m,n2),q(n2,n2),r(n2,n2),rwork(n2),
136 $ work(lwork),t(nb,m),c(n2,m),cf(n2,m),
142 CALL zlaset(
'Full', m, n2, czero, czero, a, m )
143 CALL zlaset(
'Full', nb, m, czero, czero, t, nb )
145 CALL zlarnv( 2, iseed, m-j+1, a( j, j ) )
149 CALL zlarnv( 2, iseed, m, a( 1, min(n+m,m+1) + j - 1 ) )
154 CALL zlarnv( 2, iseed, m-j+1, a( j, min(n+m,n+m-l+1)
161 CALL zlacpy(
'Full', m, n2, a, m, af, m )
165 CALL ztplqt( m,n,l,nb,af,m,af(1,np1),m,t,ldt,work,info)
169 CALL zlaset(
'Full', n2, n2, czero, one, q, n2 )
170 CALL zgemlqt(
'L',
'N', n2, n2, k, nb, af, m, t, ldt, q, n2,
175 CALL zlaset(
'Full', n2, n2, czero, czero, r, n2 )
176 CALL zlacpy(
'Lower', m, n2, af, m, r, n2 )
180 CALL zgemm(
'N',
'C', m, n2, n2, -one, a, m, q, n2, one, r, n2)
181 anorm =
zlange(
'1', m, n2, a, m, rwork )
182 resid =
zlange(
'1', m, n2, r, n2, rwork )
183 IF( anorm.GT.zero )
THEN
184 result( 1 ) = resid / (eps*anorm*max(1,n2))
191 CALL zlaset(
'Full', n2, n2, czero, one, r, n2 )
192 CALL zherk(
'U',
'N', n2, n2, dreal(-one), q, n2, dreal(one),
194 resid =
zlansy(
'1',
'Upper', n2, r, n2, rwork )
195 result( 2 ) = resid / (eps*max(1,n2))
199 CALL zlaset(
'Full', n2, m, czero, one, c, n2 )
201 CALL zlarnv( 2, iseed, n2, c( 1, j ) )
203 cnorm =
zlange(
'1', n2, m, c, n2, rwork)
204 CALL zlacpy(
'Full', n2, m, c, n2, cf, n2 )
208 CALL ztpmlqt(
'L',
'N', n,m,k,l,nb,af(1, np1),m,t,ldt,cf,n2,
209 $ cf(np1,1),n2,work,info)
213 CALL zgemm(
'N',
'N', n2, m, n2, -one, q, n2, c, n2, one, cf, n2 )
214 resid =
zlange(
'1', n2, m, cf, n2, rwork )
215 IF( cnorm.GT.zero )
THEN
216 result( 3 ) = resid / (eps*max(1,n2)*cnorm)
224 CALL zlacpy(
'Full', n2, m, c, n2, cf, n2 )
228 CALL ztpmlqt(
'L',
'C',n,m,k,l,nb,af(1,np1),m,t,ldt,cf,n2,
229 $ cf(np1,1),n2,work,info)
233 CALL zgemm(
'C',
'N',n2,m,n2,-one,q,n2,c,n2,one,cf,n2)
234 resid =
zlange(
'1', n2, m, cf, n2, rwork )
236 IF( cnorm.GT.zero )
THEN
237 result( 4 ) = resid / (eps*max(1,n2)*cnorm)
245 CALL zlarnv( 2, iseed, m, d( 1, j ) )
247 dnorm =
zlange(
'1', m, n2, d, m, rwork)
248 CALL zlacpy(
'Full', m, n2, d, m, df, m )
252 CALL ztpmlqt(
'R',
'N',m,n,k,l,nb,af(1,np1),m,t,ldt,df,m,
253 $ df(1,np1),m,work,info)
257 CALL zgemm(
'N',
'N',m,n2,n2,-one,d,m,q,n2,one,df,m)
258 resid =
zlange(
'1',m, n2,df,m,rwork )
259 IF( cnorm.GT.zero )
THEN
260 result( 5 ) = resid / (eps*max(1,n2)*dnorm)
267 CALL zlacpy(
'Full',m,n2,d,m,df,m )
271 CALL ztpmlqt(
'R',
'C',m,n,k,l,nb,af(1,np1),m,t,ldt,df,m,
272 $ df(1,np1),m,work,info)
277 CALL zgemm(
'N',
'C', m, n2, n2, -one, d, m, q, n2, one, df, m )
278 resid =
zlange(
'1', m, n2, df, m, rwork )
279 IF( cnorm.GT.zero )
THEN
280 result( 6 ) = resid / (eps*max(1,n2)*dnorm)
287 DEALLOCATE ( a, af, q, r, rwork, work, t, c, d, cf, df)