85 DOUBLE PRECISION result(6)
91 DOUBLE PRECISION,
ALLOCATABLE :: af(:,:), q(:,:),
92 $ l(:,:), 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, ll, lwork
101 DOUBLE PRECISION anorm, eps, resid, cnorm, dnorm
115 DATA iseed / 1988, 1989, 1990, 1991 /
120 lwork = max(2,ll)*max(2,ll)*nb
124 ALLOCATE ( a(m,n), af(m,n), q(n,n), l(ll,n), rwork(ll),
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 dgelqt( m, n, nb, af, m, t, ldt, work, info )
142 CALL dlaset(
'Full', n, n, zero, one, q, n )
143 CALL dgemlqt(
'R',
'N', n, n, k, nb, af, m, t, ldt, q, n,
148 CALL dlaset(
'Full', m, n, zero, zero, l, ll )
149 CALL dlacpy(
'Lower', m, n, af, m, l, ll )
153 CALL dgemm(
'N',
'T', m, n, n, -one, a, m, q, n, one, l, ll )
154 anorm =
dlange(
'1', m, n, a, m, rwork )
155 resid =
dlange(
'1', m, n, l, ll, rwork )
156 IF( anorm.GT.zero )
THEN
157 result( 1 ) = resid / (eps*max(1,m)*anorm)
164 CALL dlaset(
'Full', n, n, zero, one, l, ll )
165 CALL dsyrk(
'U',
'C', n, n, -one, q, n, one, l, ll )
166 resid =
dlansy(
'1',
'Upper', n, l, ll, rwork )
167 result( 2 ) = resid / (eps*max(1,n))
172 CALL dlarnv( 2, iseed, n, d( 1, j ) )
174 dnorm =
dlange(
'1', n, m, d, n, rwork)
175 CALL dlacpy(
'Full', n, m, d, n, df, n )
179 CALL dgemlqt(
'L',
'N', n, m, k, nb, af, m, t, nb, df, n,
184 CALL dgemm(
'N',
'N', n, m, n, -one, q, n, d, n, one, df, n )
185 resid =
dlange(
'1', n, m, df, n, rwork )
186 IF( dnorm.GT.zero )
THEN
187 result( 3 ) = resid / (eps*max(1,m)*dnorm)
194 CALL dlacpy(
'Full', n, m, d, n, df, n )
198 CALL dgemlqt(
'L',
'T', n, m, k, nb, af, m, t, nb, df, n,
203 CALL dgemm(
'T',
'N', n, m, n, -one, q, n, d, n, one, df, n )
204 resid =
dlange(
'1', n, m, df, n, rwork )
205 IF( dnorm.GT.zero )
THEN
206 result( 4 ) = resid / (eps*max(1,m)*dnorm)
214 CALL dlarnv( 2, iseed, m, c( 1, j ) )
216 cnorm =
dlange(
'1', m, n, c, m, rwork)
217 CALL dlacpy(
'Full', m, n, c, m, cf, m )
221 CALL dgemlqt(
'R',
'N', m, n, k, nb, af, m, t, nb, cf, m,
226 CALL dgemm(
'N',
'N', m, n, n, -one, c, m, q, n, one, cf, m )
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', m, n, c, m, cf, m )
240 CALL dgemlqt(
'R',
'T', m, n, k, nb, af, m, t, nb, cf, m,
245 CALL dgemm(
'N',
'T', m, n, n, -one, c, m, q, n, one, cf, m )
246 resid =
dlange(
'1', m, n, cf, m, rwork )
247 IF( cnorm.GT.zero )
THEN
248 result( 6 ) = resid / (eps*max(1,m)*dnorm)
255 DEALLOCATE ( a, af, q, l, rwork, work, t, c, d, cf, df)
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
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 dlarnv(IDIST, ISEED, N, X)
DLARNV returns a vector of random numbers from a uniform or normal distribution.
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 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
subroutine dgemlqt(SIDE, TRANS, M, N, K, MB, V, LDV, T, LDT, C, LDC, WORK, INFO)
DGEMLQT
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 ...
logical function lsame(CA, CB)
LSAME
subroutine dgelqt(M, N, MB, A, LDA, T, LDT, WORK, INFO)
DGELQT