OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
ddrvls.f
Go to the documentation of this file.
1*> \brief \b DDRVLS
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8* Definition:
9* ===========
10*
11* SUBROUTINE DDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
12* NBVAL, NXVAL, THRESH, TSTERR, A, COPYA, B,
13* COPYB, C, S, COPYS, NOUT )
14*
15* .. Scalar Arguments ..
16* LOGICAL TSTERR
17* INTEGER NM, NN, NNB, NNS, NOUT
18* DOUBLE PRECISION THRESH
19* ..
20* .. Array Arguments ..
21* LOGICAL DOTYPE( * )
22* INTEGER MVAL( * ), NBVAL( * ), NSVAL( * ),
23* $ NVAL( * ), NXVAL( * )
24* DOUBLE PRECISION A( * ), B( * ), C( * ), COPYA( * ), COPYB( * ),
25* $ COPYS( * ), S( * )
26* ..
27*
28*
29*> \par Purpose:
30* =============
31*>
32*> \verbatim
33*>
34*> DDRVLS tests the least squares driver routines DGELS, DGETSLS, DGELSS, DGELSY,
35*> and DGELSD.
36*> \endverbatim
37*
38* Arguments:
39* ==========
40*
41*> \param[in] DOTYPE
42*> \verbatim
43*> DOTYPE is LOGICAL array, dimension (NTYPES)
44*> The matrix types to be used for testing. Matrices of type j
45*> (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
46*> .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
47*> The matrix of type j is generated as follows:
48*> j=1: A = U*D*V where U and V are random orthogonal matrices
49*> and D has random entries (> 0.1) taken from a uniform
50*> distribution (0,1). A is full rank.
51*> j=2: The same of 1, but A is scaled up.
52*> j=3: The same of 1, but A is scaled down.
53*> j=4: A = U*D*V where U and V are random orthogonal matrices
54*> and D has 3*min(M,N)/4 random entries (> 0.1) taken
55*> from a uniform distribution (0,1) and the remaining
56*> entries set to 0. A is rank-deficient.
57*> j=5: The same of 4, but A is scaled up.
58*> j=6: The same of 5, but A is scaled down.
59*> \endverbatim
60*>
61*> \param[in] NM
62*> \verbatim
63*> NM is INTEGER
64*> The number of values of M contained in the vector MVAL.
65*> \endverbatim
66*>
67*> \param[in] MVAL
68*> \verbatim
69*> MVAL is INTEGER array, dimension (NM)
70*> The values of the matrix row dimension M.
71*> \endverbatim
72*>
73*> \param[in] NN
74*> \verbatim
75*> NN is INTEGER
76*> The number of values of N contained in the vector NVAL.
77*> \endverbatim
78*>
79*> \param[in] NVAL
80*> \verbatim
81*> NVAL is INTEGER array, dimension (NN)
82*> The values of the matrix column dimension N.
83*> \endverbatim
84*>
85*> \param[in] NNS
86*> \verbatim
87*> NNS is INTEGER
88*> The number of values of NRHS contained in the vector NSVAL.
89*> \endverbatim
90*>
91*> \param[in] NSVAL
92*> \verbatim
93*> NSVAL is INTEGER array, dimension (NNS)
94*> The values of the number of right hand sides NRHS.
95*> \endverbatim
96*>
97*> \param[in] NNB
98*> \verbatim
99*> NNB is INTEGER
100*> The number of values of NB and NX contained in the
101*> vectors NBVAL and NXVAL. The blocking parameters are used
102*> in pairs (NB,NX).
103*> \endverbatim
104*>
105*> \param[in] NBVAL
106*> \verbatim
107*> NBVAL is INTEGER array, dimension (NNB)
108*> The values of the blocksize NB.
109*> \endverbatim
110*>
111*> \param[in] NXVAL
112*> \verbatim
113*> NXVAL is INTEGER array, dimension (NNB)
114*> The values of the crossover point NX.
115*> \endverbatim
116*>
117*> \param[in] THRESH
118*> \verbatim
119*> THRESH is DOUBLE PRECISION
120*> The threshold value for the test ratios. A result is
121*> included in the output file if RESULT >= THRESH. To have
122*> every test ratio printed, use THRESH = 0.
123*> \endverbatim
124*>
125*> \param[in] TSTERR
126*> \verbatim
127*> TSTERR is LOGICAL
128*> Flag that indicates whether error exits are to be tested.
129*> \endverbatim
130*>
131*> \param[out] A
132*> \verbatim
133*> A is DOUBLE PRECISION array, dimension (MMAX*NMAX)
134*> where MMAX is the maximum value of M in MVAL and NMAX is the
135*> maximum value of N in NVAL.
136*> \endverbatim
137*>
138*> \param[out] COPYA
139*> \verbatim
140*> COPYA is DOUBLE PRECISION array, dimension (MMAX*NMAX)
141*> \endverbatim
142*>
143*> \param[out] B
144*> \verbatim
145*> B is DOUBLE PRECISION array, dimension (MMAX*NSMAX)
146*> where MMAX is the maximum value of M in MVAL and NSMAX is the
147*> maximum value of NRHS in NSVAL.
148*> \endverbatim
149*>
150*> \param[out] COPYB
151*> \verbatim
152*> COPYB is DOUBLE PRECISION array, dimension (MMAX*NSMAX)
153*> \endverbatim
154*>
155*> \param[out] C
156*> \verbatim
157*> C is DOUBLE PRECISION array, dimension (MMAX*NSMAX)
158*> \endverbatim
159*>
160*> \param[out] S
161*> \verbatim
162*> S is DOUBLE PRECISION array, dimension
163*> (min(MMAX,NMAX))
164*> \endverbatim
165*>
166*> \param[out] COPYS
167*> \verbatim
168*> COPYS is DOUBLE PRECISION array, dimension
169*> (min(MMAX,NMAX))
170*> \endverbatim
171*>
172*> \param[in] NOUT
173*> \verbatim
174*> NOUT is INTEGER
175*> The unit number for output.
176*> \endverbatim
177*
178* Authors:
179* ========
180*
181*> \author Univ. of Tennessee
182*> \author Univ. of California Berkeley
183*> \author Univ. of Colorado Denver
184*> \author NAG Ltd.
185*
186*> \ingroup double_lin
187*
188* =====================================================================
189 SUBROUTINE ddrvls( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
190 $ NBVAL, NXVAL, THRESH, TSTERR, A, COPYA, B,
191 $ COPYB, C, S, COPYS, NOUT )
192*
193* -- LAPACK test routine --
194* -- LAPACK is a software package provided by Univ. of Tennessee, --
195* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
196*
197* .. Scalar Arguments ..
198 LOGICAL TSTERR
199 INTEGER NM, NN, NNB, NNS, NOUT
200 DOUBLE PRECISION THRESH
201* ..
202* .. Array Arguments ..
203 LOGICAL DOTYPE( * )
204 INTEGER MVAL( * ), NBVAL( * ), NSVAL( * ),
205 $ nval( * ), nxval( * )
206 DOUBLE PRECISION A( * ), B( * ), C( * ), COPYA( * ), COPYB( * ),
207 $ COPYS( * ), S( * )
208* ..
209*
210* =====================================================================
211*
212* .. Parameters ..
213 INTEGER NTESTS
214 PARAMETER ( NTESTS = 16 )
215 INTEGER SMLSIZ
216 parameter( smlsiz = 25 )
217 DOUBLE PRECISION ONE, TWO, ZERO
218 parameter( one = 1.0d0, two = 2.0d0, zero = 0.0d0 )
219* ..
220* .. Local Scalars ..
221 CHARACTER TRANS
222 CHARACTER*3 PATH
223 INTEGER CRANK, I, IM, IMB, IN, INB, INFO, INS, IRANK,
224 $ iscale, itran, itype, j, k, lda, ldb, ldwork,
225 $ lwlsy, lwork, m, mnmin, n, nb, ncols, nerrs,
226 $ nfail, nrhs, nrows, nrun, rank, mb,
227 $ mmax, nmax, nsmax, liwork,
228 $ lwork_dgels, lwork_dgetsls, lwork_dgelss,
229 $ lwork_dgelsy, lwork_dgelsd
230 DOUBLE PRECISION EPS, NORMA, NORMB, RCOND
231* ..
232* .. Local Arrays ..
233 INTEGER ISEED( 4 ), ISEEDY( 4 ), IWQ( 1 )
234 DOUBLE PRECISION RESULT( NTESTS ), WQ( 1 )
235* ..
236* .. Allocatable Arrays ..
237 DOUBLE PRECISION, ALLOCATABLE :: WORK (:)
238 INTEGER, ALLOCATABLE :: IWORK (:)
239* ..
240* .. External Functions ..
241 DOUBLE PRECISION DASUM, DLAMCH, DQRT12, DQRT14, DQRT17
242 EXTERNAL DASUM, DLAMCH, DQRT12, DQRT14, DQRT17
243* ..
244* .. External Subroutines ..
245 EXTERNAL alaerh, alahd, alasvm, daxpy, derrls, dgels,
246 $ dgelsd, dgelss, dgelsy, dgemm, dlacpy,
247 $ dlarnv, dlasrt, dqrt13, dqrt15, dqrt16, dscal,
248 $ xlaenv
249* ..
250* .. Intrinsic Functions ..
251 INTRINSIC dble, int, log, max, min, sqrt
252* ..
253* .. Scalars in Common ..
254 LOGICAL LERR, OK
255 CHARACTER*32 SRNAMT
256 INTEGER INFOT, IOUNIT
257* ..
258* .. Common blocks ..
259 COMMON / infoc / infot, iounit, ok, lerr
260 COMMON / srnamc / srnamt
261* ..
262* .. Data statements ..
263 DATA iseedy / 1988, 1989, 1990, 1991 /
264* ..
265* .. Executable Statements ..
266*
267* Initialize constants and the random number seed.
268*
269 path( 1: 1 ) = 'Double precision'
270 path( 2: 3 ) = 'LS'
271 nrun = 0
272 nfail = 0
273 nerrs = 0
274 DO 10 i = 1, 4
275 iseed( i ) = iseedy( i )
276 10 CONTINUE
277 eps = dlamch( 'Epsilon' )
278*
279* Threshold for rank estimation
280*
281 rcond = sqrt( eps ) - ( sqrt( eps )-eps ) / 2
282*
283* Test the error exits
284*
285 CALL xlaenv( 2, 2 )
286 CALL xlaenv( 9, smlsiz )
287 IF( tsterr )
288 $ CALL derrls( path, nout )
289*
290* Print the header if NM = 0 or NN = 0 and THRESH = 0.
291*
292 IF( ( nm.EQ.0 .OR. nn.EQ.0 ) .AND. thresh.EQ.zero )
293 $ CALL alahd( nout, path )
294 infot = 0
295 CALL xlaenv( 2, 2 )
296 CALL xlaenv( 9, smlsiz )
297*
298* Compute maximal workspace needed for all routines
299*
300 nmax = 0
301 mmax = 0
302 nsmax = 0
303 DO i = 1, nm
304 IF ( mval( i ).GT.mmax ) THEN
305 mmax = mval( i )
306 END IF
307 ENDDO
308 DO i = 1, nn
309 IF ( nval( i ).GT.nmax ) THEN
310 nmax = nval( i )
311 END IF
312 ENDDO
313 DO i = 1, nns
314 IF ( nsval( i ).GT.nsmax ) THEN
315 nsmax = nsval( i )
316 END IF
317 ENDDO
318 m = mmax
319 n = nmax
320 nrhs = nsmax
321 mnmin = max( min( m, n ), 1 )
322*
323* Compute workspace needed for routines
324* DQRT14, DQRT17 (two side cases), DQRT15 and DQRT12
325*
326 lwork = max( 1, ( m+n )*nrhs,
327 $ ( n+nrhs )*( m+2 ), ( m+nrhs )*( n+2 ),
328 $ max( m+mnmin, nrhs*mnmin,2*n+m ),
329 $ max( m*n+4*mnmin+max(m,n), m*n+2*mnmin+4*n ) )
330 liwork = 1
331*
332* Iterate through all test cases and compute necessary workspace
333* sizes for ?GELS, ?GETSLS, ?GELSY, ?GELSS and ?GELSD routines.
334*
335 DO im = 1, nm
336 m = mval( im )
337 lda = max( 1, m )
338 DO in = 1, nn
339 n = nval( in )
340 mnmin = max(min( m, n ),1)
341 ldb = max( 1, m, n )
342 DO ins = 1, nns
343 nrhs = nsval( ins )
344 DO irank = 1, 2
345 DO iscale = 1, 3
346 itype = ( irank-1 )*3 + iscale
347 IF( dotype( itype ) ) THEN
348 IF( irank.EQ.1 ) THEN
349 DO itran = 1, 2
350 IF( itran.EQ.1 ) THEN
351 trans = 'N'
352 ELSE
353 trans = 'T'
354 END IF
355*
356* Compute workspace needed for DGELS
357 CALL dgels( trans, m, n, nrhs, a, lda,
358 $ b, ldb, wq, -1, info )
359 lwork_dgels = int( wq( 1 ) )
360* Compute workspace needed for DGETSLS
361 CALL dgetsls( trans, m, n, nrhs, a, lda,
362 $ b, ldb, wq, -1, info )
363 lwork_dgetsls = int( wq( 1 ) )
364 ENDDO
365 END IF
366* Compute workspace needed for DGELSY
367 CALL dgelsy( m, n, nrhs, a, lda, b, ldb, iwq,
368 $ rcond, crank, wq, -1, info )
369 lwork_dgelsy = int( wq( 1 ) )
370* Compute workspace needed for DGELSS
371 CALL dgelss( m, n, nrhs, a, lda, b, ldb, s,
372 $ rcond, crank, wq, -1 , info )
373 lwork_dgelss = int( wq( 1 ) )
374* Compute workspace needed for DGELSD
375 CALL dgelsd( m, n, nrhs, a, lda, b, ldb, s,
376 $ rcond, crank, wq, -1, iwq, info )
377 lwork_dgelsd = int( wq( 1 ) )
378* Compute LIWORK workspace needed for DGELSY and DGELSD
379 liwork = max( liwork, n, iwq( 1 ) )
380* Compute LWORK workspace needed for all functions
381 lwork = max( lwork, lwork_dgels, lwork_dgetsls,
382 $ lwork_dgelsy, lwork_dgelss,
383 $ lwork_dgelsd )
384 END IF
385 ENDDO
386 ENDDO
387 ENDDO
388 ENDDO
389 ENDDO
390*
391 lwlsy = lwork
392*
393 ALLOCATE( work( lwork ) )
394 ALLOCATE( iwork( liwork ) )
395*
396 DO 150 im = 1, nm
397 m = mval( im )
398 lda = max( 1, m )
399*
400 DO 140 in = 1, nn
401 n = nval( in )
402 mnmin = max(min( m, n ),1)
403 ldb = max( 1, m, n )
404 mb = (mnmin+1)
405*
406 DO 130 ins = 1, nns
407 nrhs = nsval( ins )
408*
409 DO 120 irank = 1, 2
410 DO 110 iscale = 1, 3
411 itype = ( irank-1 )*3 + iscale
412 IF( .NOT.dotype( itype ) )
413 $ GO TO 110
414*
415 IF( irank.EQ.1 ) THEN
416*
417* Test DGELS
418*
419* Generate a matrix of scaling type ISCALE
420*
421 CALL dqrt13( iscale, m, n, copya, lda, norma,
422 $ iseed )
423 DO 40 inb = 1, nnb
424 nb = nbval( inb )
425 CALL xlaenv( 1, nb )
426 CALL xlaenv( 3, nxval( inb ) )
427*
428 DO 30 itran = 1, 2
429 IF( itran.EQ.1 ) THEN
430 trans = 'N'
431 nrows = m
432 ncols = n
433 ELSE
434 trans = 'T'
435 nrows = n
436 ncols = m
437 END IF
438 ldwork = max( 1, ncols )
439*
440* Set up a consistent rhs
441*
442 IF( ncols.GT.0 ) THEN
443 CALL dlarnv( 2, iseed, ncols*nrhs,
444 $ work )
445 CALL dscal( ncols*nrhs,
446 $ one / dble( ncols ), work,
447 $ 1 )
448 END IF
449 CALL dgemm( trans, 'No transpose', nrows,
450 $ nrhs, ncols, one, copya, lda,
451 $ work, ldwork, zero, b, ldb )
452 CALL dlacpy( 'Full', nrows, nrhs, b, ldb,
453 $ copyb, ldb )
454*
455* Solve LS or overdetermined system
456*
457 IF( m.GT.0 .AND. n.GT.0 ) THEN
458 CALL dlacpy( 'Full', m, n, copya, lda,
459 $ a, lda )
460 CALL dlacpy( 'Full', nrows, nrhs,
461 $ copyb, ldb, b, ldb )
462 END IF
463 srnamt = 'DGELS '
464 CALL dgels( trans, m, n, nrhs, a, lda, b,
465 $ ldb, work, lwork, info )
466 IF( info.NE.0 )
467 $ CALL alaerh( path, 'DGELS ', info, 0,
468 $ trans, m, n, nrhs, -1, nb,
469 $ itype, nfail, nerrs,
470 $ nout )
471*
472* Check correctness of results
473*
474 ldwork = max( 1, nrows )
475 IF( nrows.GT.0 .AND. nrhs.GT.0 )
476 $ CALL dlacpy( 'Full', nrows, nrhs,
477 $ copyb, ldb, c, ldb )
478 CALL dqrt16( trans, m, n, nrhs, copya,
479 $ lda, b, ldb, c, ldb, work,
480 $ result( 1 ) )
481*
482 IF( ( itran.EQ.1 .AND. m.GE.n ) .OR.
483 $ ( itran.EQ.2 .AND. m.LT.n ) ) THEN
484*
485* Solving LS system
486*
487 result( 2 ) = dqrt17( trans, 1, m, n,
488 $ nrhs, copya, lda, b, ldb,
489 $ copyb, ldb, c, work,
490 $ lwork )
491 ELSE
492*
493* Solving overdetermined system
494*
495 result( 2 ) = dqrt14( trans, m, n,
496 $ nrhs, copya, lda, b, ldb,
497 $ work, lwork )
498 END IF
499*
500* Print information about the tests that
501* did not pass the threshold.
502*
503 DO 20 k = 1, 2
504 IF( result( k ).GE.thresh ) THEN
505 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
506 $ CALL alahd( nout, path )
507 WRITE( nout, fmt = 9999 )trans, m,
508 $ n, nrhs, nb, itype, k,
509 $ result( k )
510 nfail = nfail + 1
511 END IF
512 20 CONTINUE
513 nrun = nrun + 2
514 30 CONTINUE
515 40 CONTINUE
516*
517*
518* Test DGETSLS
519*
520* Generate a matrix of scaling type ISCALE
521*
522 CALL dqrt13( iscale, m, n, copya, lda, norma,
523 $ iseed )
524 DO 65 inb = 1, nnb
525 mb = nbval( inb )
526 CALL xlaenv( 1, mb )
527 DO 62 imb = 1, nnb
528 nb = nbval( imb )
529 CALL xlaenv( 2, nb )
530*
531 DO 60 itran = 1, 2
532 IF( itran.EQ.1 ) THEN
533 trans = 'N'
534 nrows = m
535 ncols = n
536 ELSE
537 trans = 'T'
538 nrows = n
539 ncols = m
540 END IF
541 ldwork = max( 1, ncols )
542*
543* Set up a consistent rhs
544*
545 IF( ncols.GT.0 ) THEN
546 CALL dlarnv( 2, iseed, ncols*nrhs,
547 $ work )
548 CALL dscal( ncols*nrhs,
549 $ one / dble( ncols ), work,
550 $ 1 )
551 END IF
552 CALL dgemm( trans, 'No transpose', nrows,
553 $ nrhs, ncols, one, copya, lda,
554 $ work, ldwork, zero, b, ldb )
555 CALL dlacpy( 'Full', nrows, nrhs, b, ldb,
556 $ copyb, ldb )
557*
558* Solve LS or overdetermined system
559*
560 IF( m.GT.0 .AND. n.GT.0 ) THEN
561 CALL dlacpy( 'Full', m, n, copya, lda,
562 $ a, lda )
563 CALL dlacpy( 'Full', nrows, nrhs,
564 $ copyb, ldb, b, ldb )
565 END IF
566 srnamt = 'DGETSLS '
567 CALL dgetsls( trans, m, n, nrhs, a,
568 $ lda, b, ldb, work, lwork, info )
569 IF( info.NE.0 )
570 $ CALL alaerh( path, 'DGETSLS ', info, 0,
571 $ trans, m, n, nrhs, -1, nb,
572 $ itype, nfail, nerrs,
573 $ nout )
574*
575* Check correctness of results
576*
577 ldwork = max( 1, nrows )
578 IF( nrows.GT.0 .AND. nrhs.GT.0 )
579 $ CALL dlacpy( 'full', NROWS, NRHS,
580 $ COPYB, LDB, C, LDB )
581 CALL DQRT16( TRANS, M, N, NRHS, COPYA,
582 $ LDA, B, LDB, C, LDB, WORK,
583 $ RESULT( 15 ) )
584*
585.EQ..AND..GE..OR. IF( ( ITRAN1 MN )
586.EQ..AND..LT. $ ( ITRAN2 MN ) ) THEN
587*
588* Solving LS system
589*
590 RESULT( 16 ) = DQRT17( TRANS, 1, M, N,
591 $ NRHS, COPYA, LDA, B, LDB,
592 $ COPYB, LDB, C, WORK,
593 $ LWORK )
594 ELSE
595*
596* Solving overdetermined system
597*
598 RESULT( 16 ) = DQRT14( TRANS, M, N,
599 $ NRHS, COPYA, LDA, B, LDB,
600 $ WORK, LWORK )
601 END IF
602*
603* Print information about the tests that
604* did not pass the threshold.
605*
606 DO 50 K = 15, 16
607.GE. IF( RESULT( K )THRESH ) THEN
608.EQ..AND..EQ. IF( NFAIL0 NERRS0 )
609 $ CALL ALAHD( NOUT, PATH )
610 WRITE( NOUT, FMT = 9997 )TRANS, M,
611 $ N, NRHS, MB, NB, ITYPE, K,
612 $ RESULT( K )
613 NFAIL = NFAIL + 1
614 END IF
615 50 CONTINUE
616 NRUN = NRUN + 2
617 60 CONTINUE
618 62 CONTINUE
619 65 CONTINUE
620 END IF
621*
622* Generate a matrix of scaling type ISCALE and rank
623* type IRANK.
624*
625 CALL DQRT15( ISCALE, IRANK, M, N, NRHS, COPYA, LDA,
626 $ COPYB, LDB, COPYS, RANK, NORMA, NORMB,
627 $ ISEED, WORK, LWORK )
628*
629* workspace used: MAX(M+MIN(M,N),NRHS*MIN(M,N),2*N+M)
630*
631 LDWORK = MAX( 1, M )
632*
633* Loop for testing different block sizes.
634*
635 DO 100 INB = 1, NNB
636 NB = NBVAL( INB )
637 CALL XLAENV( 1, NB )
638 CALL XLAENV( 3, NXVAL( INB ) )
639*
640* Test DGELSY
641*
642* DGELSY: Compute the minimum-norm solution X
643* to min( norm( A * X - B ) )
644* using the rank-revealing orthogonal
645* factorization.
646*
647* Initialize vector IWORK.
648*
649 DO 70 J = 1, N
650 IWORK( J ) = 0
651 70 CONTINUE
652*
653 CALL DLACPY( 'full', M, N, COPYA, LDA, A, LDA )
654 CALL DLACPY( 'full', M, NRHS, COPYB, LDB, B,
655 $ LDB )
656*
657 SRNAMT = 'dgelsy'
658 CALL DGELSY( M, N, NRHS, A, LDA, B, LDB, IWORK,
659 $ RCOND, CRANK, WORK, LWLSY, INFO )
660.NE. IF( INFO0 )
661 $ CALL ALAERH( PATH, 'dgelsy', INFO, 0, ' ', M,
662 $ N, NRHS, -1, NB, ITYPE, NFAIL,
663 $ NERRS, NOUT )
664*
665* Test 3: Compute relative error in svd
666* workspace: M*N + 4*MIN(M,N) + MAX(M,N)
667*
668 RESULT( 3 ) = DQRT12( CRANK, CRANK, A, LDA,
669 $ COPYS, WORK, LWORK )
670*
671* Test 4: Compute error in solution
672* workspace: M*NRHS + M
673*
674 CALL DLACPY( 'full', M, NRHS, COPYB, LDB, WORK,
675 $ LDWORK )
676 CALL DQRT16( 'no transpose', M, N, NRHS, COPYA,
677 $ LDA, B, LDB, WORK, LDWORK,
678 $ WORK( M*NRHS+1 ), RESULT( 4 ) )
679*
680* Test 5: Check norm of r'*A
681* workspace: NRHS*(M+N)
682*
683 RESULT( 5 ) = ZERO
684.GT. IF( MCRANK )
685 $ RESULT( 5 ) = DQRT17( 'no transpose', 1, M,
686 $ N, NRHS, COPYA, LDA, B, LDB,
687 $ COPYB, LDB, C, WORK, LWORK )
688*
689* Test 6: Check if x is in the rowspace of A
690* workspace: (M+NRHS)*(N+2)
691*
692 RESULT( 6 ) = ZERO
693*
694.GT. IF( NCRANK )
695 $ RESULT( 6 ) = DQRT14( 'no transpose', M, N,
696 $ NRHS, COPYA, LDA, B, LDB,
697 $ WORK, LWORK )
698*
699* Test DGELSS
700*
701* DGELSS: Compute the minimum-norm solution X
702* to min( norm( A * X - B ) )
703* using the SVD.
704*
705 CALL DLACPY( 'full', M, N, COPYA, LDA, A, LDA )
706 CALL DLACPY( 'full', M, NRHS, COPYB, LDB, B,
707 $ LDB )
708 SRNAMT = 'dgelss'
709 CALL DGELSS( M, N, NRHS, A, LDA, B, LDB, S,
710 $ RCOND, CRANK, WORK, LWORK, INFO )
711.NE. IF( INFO0 )
712 $ CALL ALAERH( PATH, 'dgelss', INFO, 0, ' ', M,
713 $ N, NRHS, -1, NB, ITYPE, NFAIL,
714 $ NERRS, NOUT )
715*
716* workspace used: 3*min(m,n) +
717* max(2*min(m,n),nrhs,max(m,n))
718*
719* Test 7: Compute relative error in svd
720*
721.GT. IF( RANK0 ) THEN
722 CALL DAXPY( MNMIN, -ONE, COPYS, 1, S, 1 )
723 RESULT( 7 ) = DASUM( MNMIN, S, 1 ) /
724 $ DASUM( MNMIN, COPYS, 1 ) /
725 $ ( EPS*DBLE( MNMIN ) )
726 ELSE
727 RESULT( 7 ) = ZERO
728 END IF
729*
730* Test 8: Compute error in solution
731*
732 CALL DLACPY( 'full', M, NRHS, COPYB, LDB, WORK,
733 $ LDWORK )
734 CALL DQRT16( 'no transpose', M, N, NRHS, COPYA,
735 $ LDA, B, LDB, WORK, LDWORK,
736 $ WORK( M*NRHS+1 ), RESULT( 8 ) )
737*
738* Test 9: Check norm of r'*A
739*
740 RESULT( 9 ) = ZERO
741.GT. IF( MCRANK )
742 $ RESULT( 9 ) = DQRT17( 'no transpose', 1, M,
743 $ N, NRHS, COPYA, LDA, B, LDB,
744 $ COPYB, LDB, C, WORK, LWORK )
745*
746* Test 10: Check if x is in the rowspace of A
747*
748 RESULT( 10 ) = ZERO
749.GT. IF( NCRANK )
750 $ RESULT( 10 ) = DQRT14( 'no transpose', M, N,
751 $ NRHS, COPYA, LDA, B, LDB,
752 $ WORK, LWORK )
753*
754* Test DGELSD
755*
756* DGELSD: Compute the minimum-norm solution X
757* to min( norm( A * X - B ) ) using a
758* divide and conquer SVD.
759*
760* Initialize vector IWORK.
761*
762 DO 80 J = 1, N
763 IWORK( J ) = 0
764 80 CONTINUE
765*
766 CALL DLACPY( 'full', M, N, COPYA, LDA, A, LDA )
767 CALL DLACPY( 'full', M, NRHS, COPYB, LDB, B,
768 $ LDB )
769*
770 SRNAMT = 'dgelsd'
771 CALL DGELSD( M, N, NRHS, A, LDA, B, LDB, S,
772 $ RCOND, CRANK, WORK, LWORK, IWORK,
773 $ INFO )
774.NE. IF( INFO0 )
775 $ CALL ALAERH( PATH, 'dgelsd', INFO, 0, ' ', M,
776 $ N, NRHS, -1, NB, ITYPE, NFAIL,
777 $ NERRS, NOUT )
778*
779* Test 11: Compute relative error in svd
780*
781.GT. IF( RANK0 ) THEN
782 CALL DAXPY( MNMIN, -ONE, COPYS, 1, S, 1 )
783 RESULT( 11 ) = DASUM( MNMIN, S, 1 ) /
784 $ DASUM( MNMIN, COPYS, 1 ) /
785 $ ( EPS*DBLE( MNMIN ) )
786 ELSE
787 RESULT( 11 ) = ZERO
788 END IF
789*
790* Test 12: Compute error in solution
791*
792 CALL DLACPY( 'full', M, NRHS, COPYB, LDB, WORK,
793 $ LDWORK )
794 CALL DQRT16( 'no transpose', M, N, NRHS, COPYA,
795 $ LDA, B, LDB, WORK, LDWORK,
796 $ WORK( M*NRHS+1 ), RESULT( 12 ) )
797*
798* Test 13: Check norm of r'*A
799*
800 RESULT( 13 ) = ZERO
801.GT. IF( MCRANK )
802 $ RESULT( 13 ) = DQRT17( 'no transpose', 1, M,
803 $ N, NRHS, COPYA, LDA, B, LDB,
804 $ COPYB, LDB, C, WORK, LWORK )
805*
806* Test 14: Check if x is in the rowspace of A
807*
808 RESULT( 14 ) = ZERO
809.GT. IF( NCRANK )
810 $ RESULT( 14 ) = DQRT14( 'no transpose', M, N,
811 $ NRHS, COPYA, LDA, B, LDB,
812 $ WORK, LWORK )
813*
814* Print information about the tests that did not
815* pass the threshold.
816*
817 DO 90 K = 3, 14
818.GE. IF( RESULT( K )THRESH ) THEN
819.EQ..AND..EQ. IF( NFAIL0 NERRS0 )
820 $ CALL ALAHD( NOUT, PATH )
821 WRITE( NOUT, FMT = 9998 )M, N, NRHS, NB,
822 $ ITYPE, K, RESULT( K )
823 NFAIL = NFAIL + 1
824 END IF
825 90 CONTINUE
826 NRUN = NRUN + 12
827*
828 100 CONTINUE
829 110 CONTINUE
830 120 CONTINUE
831 130 CONTINUE
832 140 CONTINUE
833 150 CONTINUE
834*
835* Print a summary of the results.
836*
837 CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS )
838*
839 9999 FORMAT( ' trans=''', A1, ''', m=', I5, ', n=', I5, ', nrhs=', I4,
840 $ ', nb=', I4, ', type', I2, ', test(', I2, ')=', G12.5 )
841 9998 FORMAT( ' m=', I5, ', n=', I5, ', nrhs=', I4, ', nb=', I4,
842 $ ', type', I2, ', test(', I2, ')=', G12.5 )
843 9997 FORMAT( ' trans=''', A1,' m=', I5, ', n=', I5, ', nrhs=', I4,
844 $ ', mb=', I4,', nb=', I4,', type', I2,
845 $ ', test(', I2, ')=', G12.5 )
846*
847 DEALLOCATE( WORK )
848 DEALLOCATE( IWORK )
849 RETURN
850*
851* End of DDRVLS
852*
853 END
dlarnv
subroutine dlarnv(idist, iseed, n, x)
DLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition dlarnv.f:97
dlacpy
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:103
dlasrt
subroutine dlasrt(id, n, d, info)
DLASRT sorts numbers in increasing or decreasing order.
Definition dlasrt.f:88
xlaenv
subroutine xlaenv(ispec, nvalue)
XLAENV
Definition xlaenv.f:81
alasvm
subroutine alasvm(type, nout, nfail, nrun, nerrs)
ALASVM
Definition alasvm.f:73
alahd
subroutine alahd(iounit, path)
ALAHD
Definition alahd.f:107
alaerh
subroutine alaerh(path, subnam, info, infoe, opts, m, n, kl, ku, n5, imat, nfail, nerrs, nout)
ALAERH
Definition alaerh.f:147
dgelsd
subroutine dgelsd(m, n, nrhs, a, lda, b, ldb, s, rcond, rank, work, lwork, iwork, info)
DGELSD computes the minimum-norm solution to a linear least squares problem for GE matrices
Definition dgelsd.f:209
dgetsls
subroutine dgetsls(trans, m, n, nrhs, a, lda, b, ldb, work, lwork, info)
DGETSLS
Definition dgetsls.f:162
dgelsy
subroutine dgelsy(m, n, nrhs, a, lda, b, ldb, jpvt, rcond, rank, work, lwork, info)
DGELSY solves overdetermined or underdetermined systems for GE matrices
Definition dgelsy.f:204
dgels
subroutine dgels(trans, m, n, nrhs, a, lda, b, ldb, work, lwork, info)
DGELS solves overdetermined or underdetermined systems for GE matrices
Definition dgels.f:183
dgelss
subroutine dgelss(m, n, nrhs, a, lda, b, ldb, s, rcond, rank, work, lwork, info)
DGELSS solves overdetermined or underdetermined systems for GE matrices
Definition dgelss.f:172
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
daxpy
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
dgemm
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:187
dqrt15
subroutine dqrt15(scale, rksel, m, n, nrhs, a, lda, b, ldb, s, rank, norma, normb, iseed, work, lwork)
DQRT15
Definition dqrt15.f:148
ddrvls
subroutine ddrvls(dotype, nm, mval, nn, nval, nns, nsval, nnb, nbval, nxval, thresh, tsterr, a, copya, b, copyb, c, s, copys, nout)
DDRVLS
Definition ddrvls.f:192
dqrt16
subroutine dqrt16(trans, m, n, nrhs, a, lda, x, ldx, b, ldb, rwork, resid)
DQRT16
Definition dqrt16.f:133
derrls
subroutine derrls(path, nunit)
DERRLS
Definition derrls.f:55
dqrt13
subroutine dqrt13(scale, m, n, a, lda, norma, iseed)
DQRT13
Definition dqrt13.f:91
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21