87 INTEGER LWORK, M, N, L, NB, LDT
95 COMPLEX,
ALLOCATABLE :: AF(:,:), Q(:,:),
96 $ R(:,:), WORK( : ), T(:,:),
97 $ CF(:,:), DF(:,:), A(:,:), C(:,:), D(:,:)
98 REAL,
ALLOCATABLE :: RWORK(:)
103 parameter( zero = 0.0, one = (1.0,0.0), czero=(0.0,0.0) )
106 INTEGER INFO, J, K, N2, NP1,i
107 REAL ANORM, EPS, RESID, CNORM, DNORM
116 EXTERNAL slamch, clange, clansy, lsame
119 DATA iseed / 1988, 1989, 1990, 1991 /
121 eps = slamch(
'Epsilon' )
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 claset(
'Full', m, n2, czero, czero, a, m )
141 CALL claset(
'Full', nb, m, czero, czero, t, nb )
143 CALL clarnv( 2, iseed, m-j+1, a( j, j ) )
147 CALL clarnv( 2, iseed, m, a( 1, min(n+m,m+1) + j - 1 ) )
152 CALL clarnv( 2, iseed, m-j+1, a( j, min(n+m,n+m-l+1)
159 CALL clacpy(
'Full', m, n2, a, m, af, m )
163 CALL ctplqt( m,n,l,nb,af,m,af(1,np1),m,t,ldt,work,info)
167 CALL claset(
'Full', n2, n2, czero, one, q, n2 )
168 CALL cgemlqt(
'L',
'N', n2, n2, k, nb, af, m, t, ldt, q, n2,
173 CALL claset(
'Full', n2, n2, czero, czero, r, n2 )
174 CALL clacpy(
'Lower', m, n2, af, m, r, n2 )
178 CALL cgemm(
'N',
'C', m, n2, n2, -one, a, m, q, n2, one, r, n2)
179 anorm = clange(
'1', m, n2, a, m, rwork )
180 resid = clange(
'1', m, n2, r, n2, rwork )
181 IF( anorm.GT.zero )
THEN
182 result( 1 ) = resid / (eps*anorm*max(1,n2))
189 CALL claset(
'Full', n2, n2, czero, one, r, n2 )
190 CALL cherk(
'U',
'N', n2, n2, real(-one), q, n2, real(one),
192 resid = clansy(
'1',
'Upper', n2, r, n2, rwork )
193 result( 2 ) = resid / (eps*max(1,n2))
197 CALL claset(
'Full', n2, m, czero, one, c, n2 )
199 CALL clarnv( 2, iseed, n2, c( 1, j ) )
201 cnorm = clange(
'1', n2, m, c, n2, rwork)
202 CALL clacpy(
'Full', n2, m, c, n2, cf, n2 )
206 CALL ctpmlqt(
'L',
'N', n,m,k,l,nb,af(1, np1),m,t,ldt,cf,n2,
207 $ cf(np1,1),n2,work,info)
211 CALL cgemm(
'N',
'N', n2, m, n2, -one, q, n2, c, n2, one, cf, n2 )
212 resid = clange(
'1', n2, m, cf, n2, rwork )
213 IF( cnorm.GT.zero )
THEN
214 result( 3 ) = resid / (eps*max(1,n2)*cnorm)
222 CALL clacpy(
'Full', n2, m, c, n2, cf, n2 )
226 CALL ctpmlqt(
'L',
'C',n,m,k,l,nb,af(1,np1),m,t,ldt,cf,n2,
227 $ cf(np1,1),n2,work,info)
231 CALL cgemm(
'C',
'N',n2,m,n2,-one,q,n2,c,n2,one,cf,n2)
232 resid = clange(
'1', n2, m, cf, n2, rwork )
234 IF( cnorm.GT.zero )
THEN
235 result( 4 ) = resid / (eps*max(1,n2)*cnorm)
243 CALL clarnv( 2, iseed, m, d( 1, j ) )
245 dnorm = clange(
'1', m, n2, d, m, rwork)
246 CALL clacpy(
'Full', m, n2, d, m, df, m )
250 CALL ctpmlqt(
'R',
'N',m,n,k,l,nb,af(1,np1),m,t,ldt,df,m,
251 $ df(1,np1),m,work,info)
255 CALL cgemm(
'N',
'N',m,n2,n2,-one,d,m,q,n2,one,df,m)
256 resid = clange(
'1',m, n2,df,m,rwork )
257 IF( cnorm.GT.zero )
THEN
258 result( 5 ) = resid / (eps*max(1,n2)*dnorm)
265 CALL clacpy(
'Full',m,n2,d,m,df,m )
269 CALL ctpmlqt(
'R',
'C',m,n,k,l,nb,af(1,np1),m,t,ldt,df,m,
270 $ df(1,np1),m,work,info)
275 CALL cgemm(
'N',
'C', m, n2, n2, -one, d, m, q, n2, one, df, m )
276 resid = clange(
'1', m, n2, df, m, rwork )
277 IF( cnorm.GT.zero )
THEN
278 result( 6 ) = resid / (eps*max(1,n2)*dnorm)
285 DEALLOCATE ( a, af, q, r, rwork, work, t, c, d, cf, df)