90 INTEGER LWORK, M, N, L, NB, LDT
92 DOUBLE PRECISION RESULT(6)
98 DOUBLE PRECISION,
ALLOCATABLE :: AF(:,:), Q(:,:),
99 $ R(:,:), RWORK(:), WORK( : ), T(:,:),
100 $ CF(:,:), DF(:,:), A(:,:), C(:,:), D(:,:)
103 DOUBLE PRECISION ONE, ZERO
104 parameter( zero = 0.0, one = 1.0 )
107 INTEGER INFO, J, K, N2, NP1,i
108 DOUBLE PRECISION ANORM, EPS, RESID, CNORM, DNORM
114 DOUBLE PRECISION DLAMCH, DLANGE, DLANSY
119 DATA iseed / 1988, 1989, 1990, 1991 /
133 ALLOCATE(a(m,n2),af(m,n2),q(n2,n2),r(n2,n2),rwork(n2),
134 $ work(lwork),t(nb,m),c(n2,m),cf(n2,m),
140 CALL dlaset(
'Full', m, n2, zero, zero, a, m )
141 CALL dlaset(
'Full', nb, m, zero, zero, t, nb )
143 CALL dlarnv( 2, iseed, m-j+1, a( j, j ) )
147 CALL dlarnv( 2, iseed, m, a( 1, min(n+m,m+1) + j - 1 ) )
152 CALL dlarnv( 2, iseed, m-j+1, a( j, min(n+m,n+m-l+1)
159 CALL dlacpy(
'Full', m, n2, a, m, af, m )
163 CALL dtplqt( m,n,l,nb,af,m,af(1,np1),m,t,ldt,work,info)
167 CALL dlaset(
'Full', n2, n2, zero, one, q, n2 )
168 CALL dgemlqt(
'L',
'N', n2, n2, k, nb, af, m, t, ldt, q, n2,
173 CALL dlaset(
'Full', n2, n2, zero, zero, r, n2 )
174 CALL dlacpy(
'Lower', m, n2, af, m, r, n2 )
178 CALL dgemm(
'N',
'T', m, n2, n2, -one, a, m, q, n2, one, r, n2)
179 anorm =
dlange(
'1', m, n2, a, m, rwork )
180 resid =
dlange(
'1', m, n2, r, n2, rwork )
181 IF( anorm.GT.zero )
THEN
182 result( 1 ) = resid / (eps*anorm*max(1,n2))
189 CALL dlaset(
'Full', n2, n2, zero, one, r, n2 )
190 CALL dsyrk(
'U',
'N', n2, n2, -one, q, n2, one, r, n2 )
191 resid =
dlansy(
'1',
'Upper', n2, r, n2, rwork )
192 result( 2 ) = resid / (eps*max(1,n2))
196 CALL dlaset(
'Full', n2, m, zero, one, c, n2 )
198 CALL dlarnv( 2, iseed, n2, c( 1, j ) )
200 cnorm =
dlange(
'1', n2, m, c, n2, rwork)
201 CALL dlacpy(
'Full', n2, m, c, n2, cf, n2 )
205 CALL dtpmlqt(
'L',
'N', n,m,k,l,nb,af(1, np1),m,t,ldt,cf,n2,
206 $ cf(np1,1),n2,work,info)
210 CALL dgemm(
'N',
'N', n2, m, n2, -one, q, n2, c, n2, one, cf, n2 )
211 resid =
dlange(
'1', n2, m, cf, n2, rwork )
212 IF( cnorm.GT.zero )
THEN
213 result( 3 ) = resid / (eps*max(1,n2)*cnorm)
221 CALL dlacpy(
'Full', n2, m, c, n2, cf, n2 )
225 CALL dtpmlqt(
'L',
'T',n,m,k,l,nb,af(1,np1),m,t,ldt,cf,n2,
226 $ cf(np1,1),n2,work,info)
230 CALL dgemm(
'T',
'N',n2,m,n2,-one,q,n2,c,n2,one,cf,n2)
231 resid =
dlange(
'1', n2, m, cf, n2, rwork )
233 IF( cnorm.GT.zero )
THEN
234 result( 4 ) = resid / (eps*max(1,n2)*cnorm)
242 CALL dlarnv( 2, iseed, m, d( 1, j ) )
244 dnorm =
dlange(
'1', m, n2, d, m, rwork)
245 CALL dlacpy(
'Full', m, n2, d, m, df, m )
249 CALL dtpmlqt(
'R',
'N',m,n,k,l,nb,af(1,np1),m,t,ldt,df,m,
250 $ df(1,np1),m,work,info)
254 CALL dgemm(
'N',
'N',m,n2,n2,-one,d,m,q,n2,one,df,m)
255 resid =
dlange(
'1',m, n2,df,m,rwork )
256 IF( cnorm.GT.zero )
THEN
257 result( 5 ) = resid / (eps*max(1,n2)*dnorm)
264 CALL dlacpy(
'Full',m,n2,d,m,df,m )
268 CALL dtpmlqt(
'R',
'T',m,n,k,l,nb,af(1,np1),m,t,ldt,df,m,
269 $ df(1,np1),m,work,info)
274 CALL dgemm(
'N',
'T', m, n2, n2, -one, d, m, q, n2, one, df, m )
275 resid =
dlange(
'1', m, n2, df, m, rwork )
276 IF( cnorm.GT.zero )
THEN
277 result( 6 ) = resid / (eps*max(1,n2)*dnorm)
284 DEALLOCATE ( a, af, q, r, rwork, work, t, c, d, cf, df)