73 IMPLICIT NONE
74
75
76
77
78
79
80 INTEGER M, N, NB, LDT
81
82 DOUBLE PRECISION RESULT(6)
83
84
85
86
87
88 DOUBLE PRECISION, ALLOCATABLE :: AF(:,:), Q(:,:),
89 $ R(:,:), RWORK(:), WORK( : ), T(:,:),
90 $ CF(:,:), DF(:,:), A(:,:), C(:,:), D(:,:)
91
92
93 DOUBLE PRECISION ONE, ZERO
94 parameter( zero = 0.0, one = 1.0 )
95
96
97 INTEGER INFO, J, K, L, LWORK
98 DOUBLE PRECISION ANORM, EPS, RESID, CNORM, DNORM
99
100
101 INTEGER ISEED( 4 )
102
103
104 DOUBLE PRECISION DLAMCH, DLANGE, DLANSY
105 LOGICAL LSAME
107
108
109 INTRINSIC max, min
110
111
112 DATA iseed / 1988, 1989, 1990, 1991 /
113
115 k = min(m,n)
116 l = max(m,n)
117 lwork = max(2,l)*max(2,l)*nb
118
119
120
121 ALLOCATE ( a(m,n), af(m,n), q(m,m), r(m,l), rwork(l),
122 $ work(lwork), t(nb,n), c(m,n), cf(m,n),
123 $ d(n,m), df(n,m) )
124
125
126
127 ldt=nb
128 DO j=1,n
129 CALL dlarnv( 2, iseed, m, a( 1, j ) )
130 END DO
131 CALL dlacpy(
'Full', m, n, a, m, af, m )
132
133
134
135 CALL dgeqrt( m, n, nb, af, m, t, ldt, work, info )
136
137
138
139 CALL dlaset(
'Full', m, m, zero, one, q, m )
140 CALL dgemqrt(
'R',
'N', m, m, k, nb, af, m, t, ldt, q, m,
141 $ work, info )
142
143
144
145 CALL dlaset(
'Full', m, n, zero, zero, r, m )
146 CALL dlacpy(
'Upper', m, n, af, m, r, m )
147
148
149
150 CALL dgemm(
'T',
'N', m, n, m, -one, q, m, a, m, one, r, m )
151 anorm =
dlange(
'1', m, n, a, m, rwork )
152 resid =
dlange(
'1', m, n, r, m, rwork )
153 IF( anorm.GT.zero ) THEN
154 result( 1 ) = resid / (eps*max(1,m)*anorm)
155 ELSE
156 result( 1 ) = zero
157 END IF
158
159
160
161 CALL dlaset(
'Full', m, m, zero, one, r, m )
162 CALL dsyrk(
'U',
'C', m, m, -one, q, m, one, r, m )
163 resid =
dlansy(
'1',
'Upper', m, r, m, rwork )
164 result( 2 ) = resid / (eps*max(1,m))
165
166
167
168 DO j=1,n
169 CALL dlarnv( 2, iseed, m, c( 1, j ) )
170 END DO
171 cnorm =
dlange(
'1', m, n, c, m, rwork)
172 CALL dlacpy(
'Full', m, n, c, m, cf, m )
173
174
175
176 CALL dgemqrt(
'L',
'N', m, n, k, nb, af, m, t, nb, cf, m,
177 $ work, info)
178
179
180
181 CALL dgemm(
'N',
'N', m, n, m, -one, q, m, c, m, one, cf, m )
182 resid =
dlange(
'1', m, n, cf, m, rwork )
183 IF( cnorm.GT.zero ) THEN
184 result( 3 ) = resid / (eps*max(1,m)*cnorm)
185 ELSE
186 result( 3 ) = zero
187 END IF
188
189
190
191 CALL dlacpy(
'Full', m, n, c, m, cf, m )
192
193
194
195 CALL dgemqrt(
'L',
'T', m, n, k, nb, af, m, t, nb, cf, m,
196 $ work, info)
197
198
199
200 CALL dgemm(
'T',
'N', m, n, m, -one, q, m, c, m, one, cf, m )
201 resid =
dlange(
'1', m, n, cf, m, rwork )
202 IF( cnorm.GT.zero ) THEN
203 result( 4 ) = resid / (eps*max(1,m)*cnorm)
204 ELSE
205 result( 4 ) = zero
206 END IF
207
208
209
210 DO j=1,m
211 CALL dlarnv( 2, iseed, n, d( 1, j ) )
212 END DO
213 dnorm =
dlange(
'1', n, m, d, n, rwork)
214 CALL dlacpy(
'Full', n, m, d, n, df, n )
215
216
217
218 CALL dgemqrt(
'R',
'N', n, m, k, nb, af, m, t, nb, df, n,
219 $ work, info)
220
221
222
223 CALL dgemm(
'N',
'N', n, m, m, -one, d, n, q, m, one, df, n )
224 resid =
dlange(
'1', n, m, df, n, rwork )
225 IF( cnorm.GT.zero ) THEN
226 result( 5 ) = resid / (eps*max(1,m)*dnorm)
227 ELSE
228 result( 5 ) = zero
229 END IF
230
231
232
233 CALL dlacpy(
'Full', n, m, d, n, df, n )
234
235
236
237 CALL dgemqrt(
'R',
'T', n, m, k, nb, af, m, t, nb, df, n,
238 $ work, info)
239
240
241
242 CALL dgemm(
'N',
'T', n, m, m, -one, d, n, q, m, one, df, n )
243 resid =
dlange(
'1', n, m, df, n, rwork )
244 IF( cnorm.GT.zero ) THEN
245 result( 6 ) = resid / (eps*max(1,m)*dnorm)
246 ELSE
247 result( 6 ) = zero
248 END IF
249
250
251
252 DEALLOCATE ( a, af, q, r, rwork, work, t, c, d, cf, df)
253
254 RETURN
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 dgeqrt(m, n, nb, a, lda, t, ldt, work, info)
DGEQRT
subroutine dsyrk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
DSYRK
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 dlange(norm, m, n, a, lda, work)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
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,...
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.
logical function lsame(ca, cb)
LSAME