97 DOUBLE PRECISION result(6)
103 DOUBLE PRECISION,
ALLOCATABLE :: af(:,:), q(:,:),
104 $ r(:,:), rwork(:), work( : ), t(:),
105 $ cf(:,:), df(:,:), a(:,:), c(:,:), d(:,:), lq(:,:)
108 DOUBLE PRECISION one, zero
109 parameter( zero = 0.0, one = 1.0 )
112 LOGICAL testzeros, ts
113 INTEGER info, j, k, l, lwork, tsize, mnb
114 DOUBLE PRECISION anorm, eps, resid, cnorm, dnorm
118 DOUBLE PRECISION tquery( 5 ), workquery
132 COMMON / srnamc / srnamt
135 DATA iseed / 1988, 1989, 1990, 1991 /
139 ts =
lsame(tssw,
'TS')
153 ALLOCATE ( a(m,n), af(m,n), q(l,l), r(m,l), rwork(l),
155 $ d(n,m), df(n,m), lq(l,n) )
160 CALL dlarnv( 2, iseed, m, a( 1, j ) )
165 CALL dlarnv( 2, iseed, m/2, a( m/4, j ) )
169 CALL dlacpy(
'Full', m, n, a, m, af, m )
175 CALL dgeqr( m, n, af, m, tquery, -1, workquery, -1, info )
176 tsize = int( tquery( 1 ) )
177 lwork = int( workquery )
178 CALL dgemqr(
'L',
'N', m, m, k, af, m, tquery, tsize, cf, m,
179 $ workquery, -1, info)
180 lwork = max( lwork, int( workquery ) )
181 CALL dgemqr(
'L',
'N', m, n, k, af, m, tquery, tsize, cf, m,
182 $ workquery, -1, info)
183 lwork = max( lwork, int( workquery ) )
184 CALL dgemqr(
'L',
'T', m, n, k, af, m, tquery, tsize, cf, m,
185 $ workquery, -1, info)
186 lwork = max( lwork, int( workquery ) )
187 CALL dgemqr(
'R',
'N', n, m, k, af, m, tquery, tsize, df, n,
188 $ workquery, -1, info)
189 lwork = max( lwork, int( workquery ) )
190 CALL dgemqr(
'R',
'T', n, m, k, af, m, tquery, tsize, df, n,
191 $ workquery, -1, info)
192 lwork = max( lwork, int( workquery ) )
193 ALLOCATE ( t( tsize ) )
194 ALLOCATE ( work( lwork ) )
196 CALL dgeqr( m, n, af, m, t, tsize, work, lwork, info )
200 CALL dlaset(
'Full', m, m, zero, one, q, m )
202 CALL dgemqr(
'L',
'N', m, m, k, af, m, t, tsize, q, m,
203 $ work, lwork, info )
207 CALL dlaset(
'Full', m, n, zero, zero, r, m )
208 CALL dlacpy(
'Upper', m, n, af, m, r, m )
212 CALL dgemm(
'T',
'N', m, n, m, -one, q, m, a, m, one, r, m )
213 anorm =
dlange(
'1', m, n, a, m, rwork )
214 resid =
dlange(
'1', m, n, r, m, rwork )
215 IF( anorm.GT.zero )
THEN
216 result( 1 ) = resid / (eps*max(1,m)*anorm)
223 CALL dlaset(
'Full', m, m, zero, one, r, m )
224 CALL dsyrk(
'U',
'C', m, m, -one, q, m, one, r, m )
225 resid =
dlansy(
'1',
'Upper', m, r, m, rwork )
226 result( 2 ) = resid / (eps*max(1,m))
231 CALL dlarnv( 2, iseed, m, c( 1, j ) )
233 cnorm =
dlange(
'1', m, n, c, m, rwork)
234 CALL dlacpy(
'Full', m, n, c, m, cf, m )
239 CALL dgemqr(
'L',
'N', m, n, k, af, m, t, tsize, cf, m,
244 CALL dgemm(
'N',
'N', m, n, m, -one, q, m, c, m, one, cf, m )
245 resid =
dlange(
'1', m, n, cf, m, rwork )
246 IF( cnorm.GT.zero )
THEN
247 result( 3 ) = resid / (eps*max(1,m)*cnorm)
254 CALL dlacpy(
'Full', m, n, c, m, cf, m )
259 CALL dgemqr(
'L',
'T', m, n, k, af, m, t, tsize, cf, m,
264 CALL dgemm(
'T',
'N', m, n, m, -one, q, m, c, m, one, cf, m )
265 resid =
dlange(
'1', m, n, cf, m, rwork )
266 IF( cnorm.GT.zero )
THEN
267 result( 4 ) = resid / (eps*max(1,m)*cnorm)
275 CALL dlarnv( 2, iseed, n, d( 1, j ) )
277 dnorm =
dlange(
'1', n, m, d, n, rwork)
278 CALL dlacpy(
'Full', n, m, d, n, df, n )
283 CALL dgemqr(
'R',
'N', n, m, k, af, m, t, tsize, df, n,
288 CALL dgemm(
'N',
'N', n, m, m, -one, d, n, q, m, one, df, n )
289 resid =
dlange(
'1', n, m, df, n, rwork )
290 IF( dnorm.GT.zero )
THEN
291 result( 5 ) = resid / (eps*max(1,m)*dnorm)
298 CALL dlacpy(
'Full', n, m, d, n, df, n )
302 CALL dgemqr(
'R',
'T', n, m, k, af, m, t, tsize, df, n,
307 CALL dgemm(
'N',
'T', n, m, m, -one, d, n, q, m, one, df, n )
308 resid =
dlange(
'1', n, m, df, n, rwork )
309 IF( cnorm.GT.zero )
THEN
310 result( 6 ) = resid / (eps*max(1,m)*dnorm)
318 CALL dgelq( m, n, af, m, tquery, -1, workquery, -1, info )
319 tsize = int( tquery( 1 ) )
320 lwork = int( workquery )
321 CALL dgemlq(
'R',
'N', n, n, k, af, m, tquery, tsize, q, n,
322 $ workquery, -1, info )
323 lwork = max( lwork, int( workquery ) )
324 CALL dgemlq(
'L',
'N', n, m, k, af, m, tquery, tsize, df, n,
325 $ workquery, -1, info)
326 lwork = max( lwork, int( workquery ) )
327 CALL dgemlq(
'L',
'T', n, m, k, af, m, tquery, tsize, df, n,
328 $ workquery, -1, info)
329 lwork = max( lwork, int( workquery ) )
330 CALL dgemlq(
'R',
'N', m, n, k, af, m, tquery, tsize, cf, m,
331 $ workquery, -1, info)
332 lwork = max( lwork, int( workquery ) )
333 CALL dgemlq(
'R',
'T', m, n, k, af, m, tquery, tsize, cf, m,
334 $ workquery, -1, info)
335 lwork = max( lwork, int( workquery ) )
336 ALLOCATE ( t( tsize ) )
337 ALLOCATE ( work( lwork ) )
339 CALL dgelq( m, n, af, m, t, tsize, work, lwork, info )
344 CALL dlaset(
'Full', n, n, zero, one, q, n )
346 CALL dgemlq(
'R',
'N', n, n, k, af, m, t, tsize, q, n,
347 $ work, lwork, info )
351 CALL dlaset(
'Full', m, n, zero, zero, lq, l )
352 CALL dlacpy(
'Lower', m, n, af, m, lq, l )
356 CALL dgemm(
'N',
'T', m, n, n, -one, a, m, q, n, one, lq, l )
357 anorm =
dlange(
'1', m, n, a, m, rwork )
358 resid =
dlange(
'1', m, n, lq, l, rwork )
359 IF( anorm.GT.zero )
THEN
360 result( 1 ) = resid / (eps*max(1,n)*anorm)
367 CALL dlaset(
'Full', n, n, zero, one, lq, l )
368 CALL dsyrk(
'U',
'C', n, n, -one, q, n, one, lq, l )
369 resid =
dlansy(
'1',
'Upper', n, lq, l, rwork )
370 result( 2 ) = resid / (eps*max(1,n))
375 CALL dlarnv( 2, iseed, n, d( 1, j ) )
377 dnorm =
dlange(
'1', n, m, d, n, rwork)
378 CALL dlacpy(
'Full', n, m, d, n, df, n )
382 CALL dgemlq(
'L',
'N', n, m, k, af, m, t, tsize, df, n,
387 CALL dgemm(
'N',
'N', n, m, n, -one, q, n, d, n, one, df, n )
388 resid =
dlange(
'1', n, m, df, n, rwork )
389 IF( dnorm.GT.zero )
THEN
390 result( 3 ) = resid / (eps*max(1,n)*dnorm)
397 CALL dlacpy(
'Full', n, m, d, n, df, n )
401 CALL dgemlq(
'L',
'T', n, m, k, af, m, t, tsize, df, n,
406 CALL dgemm(
'T',
'N', n, m, n, -one, q, n, d, n, one, df, n )
407 resid =
dlange(
'1', n, m, df, n, rwork )
408 IF( dnorm.GT.zero )
THEN
409 result( 4 ) = resid / (eps*max(1,n)*dnorm)
417 CALL dlarnv( 2, iseed, m, c( 1, j ) )
419 cnorm =
dlange(
'1', m, n, c, m, rwork)
420 CALL dlacpy(
'Full', m, n, c, m, cf, m )
424 CALL dgemlq(
'R',
'N', m, n, k, af, m, t, tsize, cf, m,
429 CALL dgemm(
'N',
'N', m, n, n, -one, c, m, q, n, one, cf, m )
430 resid =
dlange(
'1', n, m, df, n, rwork )
431 IF( cnorm.GT.zero )
THEN
432 result( 5 ) = resid / (eps*max(1,n)*cnorm)
439 CALL dlacpy(
'Full', m, n, c, m, cf, m )
443 CALL dgemlq(
'R',
'T', m, n, k, af, m, t, tsize, cf, m,
448 CALL dgemm(
'N',
'T', m, n, n, -one, c, m, q, n, one, cf, m )
449 resid =
dlange(
'1', m, n, cf, m, rwork )
450 IF( cnorm.GT.zero )
THEN
451 result( 6 ) = resid / (eps*max(1,n)*cnorm)
460 DEALLOCATE ( a, af, q, r, 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.
subroutine dgemlq(SIDE, TRANS, M, N, K, A, LDA, T, TSIZE, C, LDC, WORK, LWORK, INFO)
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 dgeqr(M, N, A, LDA, T, TSIZE, WORK, LWORK, INFO)
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
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
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 dgemqr(SIDE, TRANS, M, N, K, A, LDA, T, TSIZE, C, LDC, WORK, LWORK, INFO)
subroutine dgelq(M, N, A, LDA, T, TSIZE, WORK, LWORK, INFO)