90 INTEGER lwork, m, n, l, nb, ldt
92 DOUBLE PRECISION result(6)
98 COMPLEX*16,
ALLOCATABLE :: af(:,:), q(:,:),
99 $ r(:,:), work( : ), t(:,:),
100 $ cf(:,:), df(:,:), a(:,:), c(:,:), d(:,:)
101 DOUBLE PRECISION,
ALLOCATABLE :: rwork(:)
104 DOUBLE PRECISION zero
105 COMPLEX*16 one, czero
106 parameter( zero = 0.0, one = (1.0,0.0), czero=(0.0,0.0) )
109 INTEGER info, j, k, m2, np1
110 DOUBLE PRECISION 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 zlaset(
'Full', m2, n, czero, czero, a, m2 )
144 CALL zlaset(
'Full', nb, n, czero, czero, t, nb )
146 CALL zlarnv( 2, iseed, j, a( 1, j ) )
150 CALL zlarnv( 2, iseed, m-l, a( min(n+m,n+1), j ) )
155 CALL zlarnv( 2, iseed, min(j,l), a( min(n+m,n+m-l+1), j ) )
161 CALL zlacpy(
'Full', m2, n, a, m2, af, m2 )
165 CALL ztpqrt( m,n,l,nb,af,m2,af(np1,1),m2,t,ldt,work,info)
169 CALL zlaset(
'Full', m2, m2, czero, one, q, m2 )
170 CALL zgemqrt(
'R',
'N', m2, m2, k, nb, af, m2, t, ldt, q, m2,
175 CALL zlaset(
'Full', m2, n, czero, czero, r, m2 )
176 CALL zlacpy(
'Upper', m2, n, af, m2, r, m2 )
180 CALL zgemm(
'C',
'N', m2, n, m2, -one, q, m2, a, m2, one, r, m2 )
181 anorm =
zlange(
'1', m2, n, a, m2, rwork )
182 resid =
zlange(
'1', m2, n, r, m2, rwork )
183 IF( anorm.GT.zero )
THEN
184 result( 1 ) = resid / (eps*anorm*max(1,m2))
191 CALL zlaset(
'Full', m2, m2, czero, one, r, m2 )
192 CALL zherk(
'U',
'C', m2, m2, dreal(-one), q, m2, dreal(one),
194 resid =
zlansy(
'1',
'Upper', m2, r, m2, rwork )
195 result( 2 ) = resid / (eps*max(1,m2))
200 CALL zlarnv( 2, iseed, m2, c( 1, j ) )
202 cnorm =
zlange(
'1', m2, n, c, m2, rwork)
203 CALL zlacpy(
'Full', m2, n, c, m2, cf, m2 )
207 CALL ztpmqrt(
'L',
'N', m,n,k,l,nb,af(np1,1),m2,t,ldt,cf,m2,
208 $ cf(np1,1),m2,work,info)
212 CALL zgemm(
'N',
'N', m2, n, m2, -one, q, m2, c, m2, one, cf, m2 )
213 resid =
zlange(
'1', m2, n, cf, m2, rwork )
214 IF( cnorm.GT.zero )
THEN
215 result( 3 ) = resid / (eps*max(1,m2)*cnorm)
222 CALL zlacpy(
'Full', m2, n, c, m2, cf, m2 )
226 CALL ztpmqrt(
'L',
'C',m,n,k,l,nb,af(np1,1),m2,t,ldt,cf,m2,
227 $ cf(np1,1),m2,work,info)
231 CALL zgemm(
'C',
'N',m2,n,m2,-one,q,m2,c,m2,one,cf,m2)
232 resid =
zlange(
'1', m2, n, cf, m2, rwork )
233 IF( cnorm.GT.zero )
THEN
234 result( 4 ) = resid / (eps*max(1,m2)*cnorm)
242 CALL zlarnv( 2, iseed, n, d( 1, j ) )
244 dnorm =
zlange(
'1', n, m2, d, n, rwork)
245 CALL zlacpy(
'Full', n, m2, d, n, df, n )
249 CALL ztpmqrt(
'R',
'N',n,m,n,l,nb,af(np1,1),m2,t,ldt,df,n,
250 $ df(1,np1),n,work,info)
254 CALL zgemm(
'N',
'N',n,m2,m2,-one,d,n,q,m2,one,df,n)
255 resid =
zlange(
'1',n, m2,df,n,rwork )
256 IF( cnorm.GT.zero )
THEN
257 result( 5 ) = resid / (eps*max(1,m2)*dnorm)
264 CALL zlacpy(
'Full',n,m2,d,n,df,n )
268 CALL ztpmqrt(
'R',
'C',n,m,n,l,nb,af(np1,1),m2,t,ldt,df,n,
269 $ df(1,np1),n,work,info)
274 CALL zgemm(
'N',
'C', n, m2, m2, -one, d, n, q, m2, one, df, n )
275 resid =
zlange(
'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)
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine zlarnv(IDIST, ISEED, N, X)
ZLARNV returns a vector of random numbers from a uniform or normal distribution.
subroutine ztpqrt(M, N, L, NB, A, LDA, B, LDB, T, LDT, WORK, INFO)
ZTPQRT
double precision function dlamch(CMACH)
DLAMCH
subroutine ztpmqrt(SIDE, TRANS, M, N, K, L, NB, V, LDV, T, LDT, A, LDA, B, LDB, WORK, INFO)
ZTPMQRT
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
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...
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, or the element of largest absolute value of a complex symmetric matrix.
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
logical function lsame(CA, CB)
LSAME