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)
LOGICAL function lsame(CA, CB)
LSAME
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.
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
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.
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 dsyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
DSYRK
subroutine dtpqrt(M, N, L, NB, A, LDA, B, LDB, T, LDT, WORK, INFO)
DTPQRT
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
DOUBLE PRECISION function dlamch(CMACH)
DLAMCH
subroutine dqrt05(M, N, L, NB, RESULT)
DQRT05
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...
subroutine dtpmqrt(SIDE, TRANS, M, N, K, L, NB, V, LDV, T, LDT, A, LDA, B, LDB, WORK, INFO)
DTPMQRT