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, m2, np1
108 DOUBLE PRECISION anorm, eps, resid, cnorm, dnorm
119 DATA iseed / 1988, 1989, 1990, 1991 /
133 ALLOCATE(a(m2,n),af(m2,n),q(m2,m2),r(m2,m2),rwork(m2),
134 $ work(lwork),t(nb,n),c(m2,n),cf(m2,n),
140 CALL dlaset(
'Full', m2, n, zero, zero, a, m2 )
141 CALL dlaset(
'Full', nb, n, zero, zero, t, nb )
143 CALL dlarnv( 2, iseed, j, a( 1, j ) )
147 CALL dlarnv( 2, iseed, m-l, a( min(n+m,n+1), j ) )
152 CALL dlarnv( 2, iseed, min(j,l), a( min(n+m,n+m-l+1), j ) )
158 CALL dlacpy(
'Full', m2, n, a, m2, af, m2 )
162 CALL dtpqrt( m,n,l,nb,af,m2,af(np1,1),m2,t,ldt,work,info)
166 CALL dlaset(
'Full', m2, m2, zero, one, q, m2 )
167 CALL dgemqrt(
'R',
'N', m2, m2, k, nb, af, m2, t, ldt, q, m2,
172 CALL dlaset(
'Full', m2, n, zero, zero, r, m2 )
173 CALL dlacpy(
'Upper', m2, n, af, m2, r, m2 )
177 CALL dgemm(
'T',
'N', m2, n, m2, -one, q, m2, a, m2, one, r, m2 )
178 anorm =
dlange(
'1', m2, n, a, m2, rwork )
179 resid =
dlange(
'1', m2, n, r, m2, rwork )
180 IF( anorm.GT.zero )
THEN
181 result( 1 ) = resid / (eps*anorm*max(1,m2))
188 CALL dlaset(
'Full', m2, m2, zero, one, r, m2 )
189 CALL dsyrk(
'U',
'C', m2, m2, -one, q, m2, one, r, m2 )
190 resid =
dlansy(
'1',
'Upper', m2, r, m2, rwork )
191 result( 2 ) = resid / (eps*max(1,m2))
196 CALL dlarnv( 2, iseed, m2, c( 1, j ) )
198 cnorm =
dlange(
'1', m2, n, c, m2, rwork)
199 CALL dlacpy(
'Full', m2, n, c, m2, cf, m2 )
203 CALL dtpmqrt(
'L',
'N', m,n,k,l,nb,af(np1,1),m2,t,ldt,cf,m2,
204 $ cf(np1,1),m2,work,info)
208 CALL dgemm(
'N',
'N', m2, n, m2, -one, q, m2, c, m2, one, cf, m2 )
209 resid =
dlange(
'1', m2, n, cf, m2, rwork )
210 IF( cnorm.GT.zero )
THEN
211 result( 3 ) = resid / (eps*max(1,m2)*cnorm)
218 CALL dlacpy(
'Full', m2, n, c, m2, cf, m2 )
222 CALL dtpmqrt(
'L',
'T',m,n,k,l,nb,af(np1,1),m2,t,ldt,cf,m2,
223 $ cf(np1,1),m2,work,info)
227 CALL dgemm(
'T',
'N',m2,n,m2,-one,q,m2,c,m2,one,cf,m2)
228 resid =
dlange(
'1', m2, n, cf, m2, rwork )
229 IF( cnorm.GT.zero )
THEN
230 result( 4 ) = resid / (eps*max(1,m2)*cnorm)
238 CALL dlarnv( 2, iseed, n, d( 1, j ) )
240 dnorm =
dlange(
'1', n, m2, d, n, rwork)
241 CALL dlacpy(
'Full', n, m2, d, n, df, n )
245 CALL dtpmqrt(
'R',
'N',n,m,n,l,nb,af(np1,1),m2,t,ldt,df,n,
246 $ df(1,np1),n,work,info)
250 CALL dgemm(
'N',
'N',n,m2,m2,-one,d,n,q,m2,one,df,n)
251 resid =
dlange(
'1',n, m2,df,n,rwork )
252 IF( cnorm.GT.zero )
THEN
253 result( 5 ) = resid / (eps*max(1,m2)*dnorm)
260 CALL dlacpy(
'Full',n,m2,d,n,df,n )
264 CALL dtpmqrt(
'R',
'T',n,m,n,l,nb,af(np1,1),m2,t,ldt,df,n,
265 $ df(1,np1),n,work,info)
270 CALL dgemm(
'N',
'T', n, m2, m2, -one, d, n, q, m2, one, df, n )
271 resid =
dlange(
'1', n, m2, df, n, rwork )
272 IF( cnorm.GT.zero )
THEN
273 result( 6 ) = resid / (eps*max(1,m2)*dnorm)
280 DEALLOCATE ( a, af, q, r, rwork, work, t, c, d, cf, df)
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
double precision function dlansy(NORM, UPLO, N, A, LDA, WORK)
DLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a real symmetric matrix.
double precision function dlamch(CMACH)
DLAMCH
subroutine dtpmqrt(SIDE, TRANS, M, N, K, L, NB, V, LDV, T, LDT, A, LDA, B, LDB, WORK, INFO)
DTPMQRT
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
subroutine dtpqrt(M, N, L, NB, A, LDA, B, LDB, T, LDT, WORK, INFO)
DTPQRT
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
subroutine dsyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
DSYRK
double precision function dlange(NORM, M, N, A, LDA, WORK)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine dgemqrt(SIDE, TRANS, M, N, K, NB, V, LDV, T, LDT, C, LDC, WORK, INFO)
DGEMQRT
subroutine dlarnv(IDIST, ISEED, N, X)
DLARNV returns a vector of random numbers from a uniform or normal distribution.
logical function lsame(CA, CB)
LSAME