90 INTEGER LWORK, M, N, L, NB, LDT
98 COMPLEX,
ALLOCATABLE :: AF(:,:), Q(:,:),
99 $ R(:,:), WORK( : ), T(:,:),
100 $ CF(:,:), DF(:,:), A(:,:), C(:,:), D(:,:)
101 REAL,
ALLOCATABLE :: RWORK(:)
106 parameter( zero = 0.0, one = (1.0,0.0), czero=(0.0,0.0) )
109 INTEGER INFO, J, K, M2, NP1
110 REAL ANORM, EPS, RESID, CNORM, DNORM
122 DATA iseed / 1988, 1989, 1990, 1991 /
136 ALLOCATE(a(m2,n),af(m2,n),q(m2,m2),r(m2,m2),rwork(m2),
137 $ work(lwork),t(nb,n),c(m2,n),cf(m2,n),
143 CALL claset(
'Full', m2, n, czero, czero, a, m2 )
144 CALL claset(
'Full', nb, n, czero, czero, t, nb )
146 CALL clarnv( 2, iseed, j, a( 1, j ) )
150 CALL clarnv( 2, iseed, m-l, a( min(n+m,n+1), j ) )
155 CALL clarnv( 2, iseed, min(j,l), a( min(n+m,n+m-l+1), j ) )
161 CALL clacpy(
'Full', m2, n, a, m2, af, m2 )
165 CALL ctpqrt( m,n,l,nb,af,m2,af(np1,1),m2,t,ldt,work,info)
169 CALL claset(
'Full', m2, m2, czero, one, q, m2 )
170 CALL cgemqrt(
'R',
'N', m2, m2, k, nb, af, m2, t, ldt, q, m2,
175 CALL claset(
'Full', m2, n, czero, czero, r, m2 )
176 CALL clacpy(
'Upper', m2, n, af, m2, r, m2 )
180 CALL cgemm(
'C',
'N', m2, n, m2, -one, q, m2, a, m2, one, r, m2 )
181 anorm =
clange(
'1', m2, n, a, m2, rwork )
182 resid =
clange(
'1', m2, n, r, m2, rwork )
183 IF( anorm.GT.zero )
THEN
184 result( 1 ) = resid / (eps*anorm*max(1,m2))
191 CALL claset(
'Full', m2, m2, czero, one, r, m2 )
192 CALL cherk(
'U',
'C', m2, m2, real(-one), q, m2, real(one),
194 resid =
clansy(
'1',
'Upper', m2, r, m2, rwork )
195 result( 2 ) = resid / (eps*max(1,m2))
200 CALL clarnv( 2, iseed, m2, c( 1, j ) )
202 cnorm =
clange(
'1', m2, n, c, m2, rwork)
203 CALL clacpy(
'Full', m2, n, c, m2, cf, m2 )
207 CALL ctpmqrt(
'L',
'N', m,n,k,l,nb,af(np1,1),m2,t,ldt,cf,m2,
208 $ cf(np1,1),m2,work,info)
212 CALL cgemm(
'N',
'N', m2, n, m2, -one, q, m2, c, m2, one, cf, m2 )
213 resid =
clange(
'1', m2, n, cf, m2, rwork )
214 IF( cnorm.GT.zero )
THEN
215 result( 3 ) = resid / (eps*max(1,m2)*cnorm)
222 CALL clacpy(
'Full', m2, n, c, m2, cf, m2 )
226 CALL ctpmqrt(
'L',
'C',m,n,k,l,nb,af(np1,1),m2,t,ldt,cf,m2,
227 $ cf(np1,1),m2,work,info)
231 CALL cgemm(
'C',
'N',m2,n,m2,-one,q,m2,c,m2,one,cf,m2)
232 resid =
clange(
'1', m2, n, cf, m2, rwork )
233 IF( cnorm.GT.zero )
THEN
234 result( 4 ) = resid / (eps*max(1,m2)*cnorm)
242 CALL clarnv( 2, iseed, n, d( 1, j ) )
244 dnorm =
clange(
'1', n, m2, d, n, rwork)
245 CALL clacpy(
'Full', n, m2, d, n, df, n )
249 CALL ctpmqrt(
'R',
'N',n,m,n,l,nb,af(np1,1),m2,t,ldt,df,n,
250 $ df(1,np1),n,work,info)
254 CALL cgemm(
'N',
'N',n,m2,m2,-one,d,n,q,m2,one,df,n)
255 resid =
clange(
'1',n, m2,df,n,rwork )
256 IF( cnorm.GT.zero )
THEN
257 result( 5 ) = resid / (eps*max(1,m2)*dnorm)
264 CALL clacpy(
'Full',n,m2,d,n,df,n )
268 CALL ctpmqrt(
'R',
'C',n,m,n,l,nb,af(np1,1),m2,t,ldt,df,n,
269 $ df(1,np1),n,work,info)
274 CALL cgemm(
'N',
'C', n, m2, m2, -one, d, n, q, m2, one, df, n )
275 resid =
clange(
'1', n, m2, df, n, rwork )
276 IF( cnorm.GT.zero )
THEN
277 result( 6 ) = resid / (eps*max(1,m2)*dnorm)
284 DEALLOCATE ( a, af, q, r, rwork, work, t, c, d, cf, df)