85 DOUBLE PRECISION result(6)
91 DOUBLE PRECISION,
ALLOCATABLE :: af(:,:), q(:,:),
92 $ r(:,:), rwork(:), work( : ), t(:,:),
93 $ cf(:,:), df(:,:), a(:,:), c(:,:), d(:,:)
96 DOUBLE PRECISION one, zero
97 parameter( zero = 0.0, one = 1.0 )
100 INTEGER info,
j, k, l, lwork
101 DOUBLE PRECISION anorm, eps, resid, cnorm, dnorm
115 DATA iseed / 1988, 1989, 1990, 1991 /
120 lwork = max(2,l)*max(2,l)*nb
124 ALLOCATE ( a(m,n), af(m,n), q(m,m), r(m,l), rwork(l),
125 $ work(lwork), t(nb,n), c(m,n), cf(m,n),
132 CALL
dlarnv( 2, iseed, m, a( 1,
j ) )
134 CALL
dlacpy(
'Full', m, n, a, m, af, m )
138 CALL
dgeqrt( m, n, nb, af, m, t, ldt, work, info )
142 CALL
dlaset(
'Full', m, m, zero, one, q, m )
143 CALL
dgemqrt(
'R',
'N', m, m, k, nb, af, m, t, ldt, q, m,
148 CALL
dlaset(
'Full', m, n, zero, zero, r, m )
149 CALL
dlacpy(
'Upper', m, n, af, m, r, m )
153 CALL
dgemm(
'T',
'N', m, n, m, -one, q, m, a, m, one, r, m )
154 anorm =
dlange(
'1', m, n, a, m, rwork )
155 resid =
dlange(
'1', m, n, r, m, rwork )
156 IF( anorm.GT.zero )
THEN
157 result( 1 ) = resid / (eps*max(1,m)*anorm)
164 CALL
dlaset(
'Full', m, m, zero, one, r, m )
165 CALL
dsyrk(
'U',
'C', m, m, -one, q, m, one, r, m )
166 resid =
dlansy(
'1',
'Upper', m, r, m, rwork )
167 result( 2 ) = resid / (eps*max(1,m))
172 CALL
dlarnv( 2, iseed, m, c( 1,
j ) )
174 cnorm =
dlange(
'1', m, n, c, m, rwork)
175 CALL
dlacpy(
'Full', m, n, c, m, cf, m )
179 CALL
dgemqrt(
'L',
'N', m, n, k, nb, af, m, t, nb, cf, m,
184 CALL
dgemm(
'N',
'N', m, n, m, -one, q, m, c, m, one, cf, m )
185 resid =
dlange(
'1', m, n, cf, m, rwork )
186 IF( cnorm.GT.zero )
THEN
187 result( 3 ) = resid / (eps*max(1,m)*cnorm)
194 CALL
dlacpy(
'Full', m, n, c, m, cf, m )
198 CALL
dgemqrt(
'L',
'T', m, n, k, nb, af, m, t, nb, cf, m,
203 CALL
dgemm(
'T',
'N', m, n, m, -one, q, m, c, m, one, cf, m )
204 resid =
dlange(
'1', m, n, cf, m, rwork )
205 IF( cnorm.GT.zero )
THEN
206 result( 4 ) = resid / (eps*max(1,m)*cnorm)
214 CALL
dlarnv( 2, iseed, n, d( 1,
j ) )
216 dnorm =
dlange(
'1', n, m, d, n, rwork)
217 CALL
dlacpy(
'Full', n, m, d, n, df, n )
221 CALL
dgemqrt(
'R',
'N', n, m, k, nb, af, m, t, nb, df, n,
226 CALL
dgemm(
'N',
'N', n, m, m, -one, d, n, q, m, one, df, n )
227 resid =
dlange(
'1', n, m, df, n, rwork )
228 IF( cnorm.GT.zero )
THEN
229 result( 5 ) = resid / (eps*max(1,m)*dnorm)
236 CALL
dlacpy(
'Full', n, m, d, n, df, n )
240 CALL
dgemqrt(
'R',
'T', n, m, k, nb, af, m, t, nb, df, n,
245 CALL
dgemm(
'N',
'T', n, m, m, -one, d, n, q, m, one, df, n )
246 resid =
dlange(
'1', n, m, df, n, rwork )
247 IF( cnorm.GT.zero )
THEN
248 result( 6 ) = resid / (eps*max(1,m)*dnorm)
255 DEALLOCATE ( a, af, q, r, rwork, work, t, c, d, cf, df)
LOGICAL function lsame(CA, CB)
LSAME
subroutine dqrt04(M, N, NB, RESULT)
DQRT04
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
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 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 dgeqrt(M, N, NB, A, LDA, T, LDT, WORK, INFO)
DGEQRT