LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
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, DGELST,
35*> DGETSLS, DGELSS, DGELSY, 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 = 18 )
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_dgelst, lwork_dgetsls,
229 $ lwork_dgelss, 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, dgelst, dgelsy, dgemm,
247 $ dgetsls, dlacpy, dlarnv, dqrt13, dqrt15,
248 $ dqrt16, dscal, xlaenv
249* ..
250* .. Intrinsic Functions ..
251 INTRINSIC dble, int, 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, ?GELST, ?GETSLS, ?GELSY, ?GELSS and ?GELSD
334* routines.
335*
336 DO im = 1, nm
337 m = mval( im )
338 lda = max( 1, m )
339 DO in = 1, nn
340 n = nval( in )
341 mnmin = max(min( m, n ),1)
342 ldb = max( 1, m, n )
343 DO ins = 1, nns
344 nrhs = nsval( ins )
345 DO irank = 1, 2
346 DO iscale = 1, 3
347 itype = ( irank-1 )*3 + iscale
348 IF( dotype( itype ) ) THEN
349 IF( irank.EQ.1 ) THEN
350 DO itran = 1, 2
351 IF( itran.EQ.1 ) THEN
352 trans = 'N'
353 ELSE
354 trans = 'T'
355 END IF
356*
357* Compute workspace needed for DGELS
358 CALL dgels( trans, m, n, nrhs, a, lda,
359 $ b, ldb, wq, -1, info )
360 lwork_dgels = int( wq( 1 ) )
361* Compute workspace needed for DGELST
362 CALL dgelst( trans, m, n, nrhs, a, lda,
363 $ b, ldb, wq, -1, info )
364 lwork_dgelst = int( wq( 1 ) )
365* Compute workspace needed for DGETSLS
366 CALL dgetsls( trans, m, n, nrhs, a, lda,
367 $ b, ldb, wq, -1, info )
368 lwork_dgetsls = int( wq( 1 ) )
369 ENDDO
370 END IF
371* Compute workspace needed for DGELSY
372 CALL dgelsy( m, n, nrhs, a, lda, b, ldb, iwq,
373 $ rcond, crank, wq, -1, info )
374 lwork_dgelsy = int( wq( 1 ) )
375* Compute workspace needed for DGELSS
376 CALL dgelss( m, n, nrhs, a, lda, b, ldb, s,
377 $ rcond, crank, wq, -1 , info )
378 lwork_dgelss = int( wq( 1 ) )
379* Compute workspace needed for DGELSD
380 CALL dgelsd( m, n, nrhs, a, lda, b, ldb, s,
381 $ rcond, crank, wq, -1, iwq, info )
382 lwork_dgelsd = int( wq( 1 ) )
383* Compute LIWORK workspace needed for DGELSY and DGELSD
384 liwork = max( liwork, n, iwq( 1 ) )
385* Compute LWORK workspace needed for all functions
386 lwork = max( lwork, lwork_dgels, lwork_dgelst,
387 $ lwork_dgetsls, lwork_dgelsy,
388 $ lwork_dgelss, lwork_dgelsd )
389 END IF
390 ENDDO
391 ENDDO
392 ENDDO
393 ENDDO
394 ENDDO
395*
396 lwlsy = lwork
397*
398 ALLOCATE( work( lwork ) )
399 ALLOCATE( iwork( liwork ) )
400*
401 DO 150 im = 1, nm
402 m = mval( im )
403 lda = max( 1, m )
404*
405 DO 140 in = 1, nn
406 n = nval( in )
407 mnmin = max(min( m, n ),1)
408 ldb = max( 1, m, n )
409 mb = (mnmin+1)
410*
411 DO 130 ins = 1, nns
412 nrhs = nsval( ins )
413*
414 DO 120 irank = 1, 2
415 DO 110 iscale = 1, 3
416 itype = ( irank-1 )*3 + iscale
417 IF( .NOT.dotype( itype ) )
418 $ GO TO 110
419* =====================================================
420* Begin test DGELS
421* =====================================================
422 IF( irank.EQ.1 ) THEN
423*
424* Generate a matrix of scaling type ISCALE
425*
426 CALL dqrt13( iscale, m, n, copya, lda, norma,
427 $ iseed )
428*
429* Loop for testing different block sizes.
430*
431 DO inb = 1, nnb
432 nb = nbval( inb )
433 CALL xlaenv( 1, nb )
434 CALL xlaenv( 3, nxval( inb ) )
435*
436* Loop for testing non-transposed and transposed.
437*
438 DO itran = 1, 2
439 IF( itran.EQ.1 ) THEN
440 trans = 'N'
441 nrows = m
442 ncols = n
443 ELSE
444 trans = 'T'
445 nrows = n
446 ncols = m
447 END IF
448 ldwork = max( 1, ncols )
449*
450* Set up a consistent rhs
451*
452 IF( ncols.GT.0 ) THEN
453 CALL dlarnv( 2, iseed, ncols*nrhs,
454 $ work )
455 CALL dscal( ncols*nrhs,
456 $ one / dble( ncols ), work,
457 $ 1 )
458 END IF
459 CALL dgemm( trans, 'No transpose', nrows,
460 $ nrhs, ncols, one, copya, lda,
461 $ work, ldwork, zero, b, ldb )
462 CALL dlacpy( 'Full', nrows, nrhs, b, ldb,
463 $ copyb, ldb )
464*
465* Solve LS or overdetermined system
466*
467 IF( m.GT.0 .AND. n.GT.0 ) THEN
468 CALL dlacpy( 'Full', m, n, copya, lda,
469 $ a, lda )
470 CALL dlacpy( 'Full', nrows, nrhs,
471 $ copyb, ldb, b, ldb )
472 END IF
473 srnamt = 'DGELS '
474 CALL dgels( trans, m, n, nrhs, a, lda, b,
475 $ ldb, work, lwork, info )
476 IF( info.NE.0 )
477 $ CALL alaerh( path, 'DGELS ', info, 0,
478 $ trans, m, n, nrhs, -1, nb,
479 $ itype, nfail, nerrs,
480 $ nout )
481*
482* Test 1: Check correctness of results
483* for DGELS, compute the residual:
484* RESID = norm(B - A*X) /
485* / ( max(m,n) * norm(A) * norm(X) * EPS )
486*
487 IF( nrows.GT.0 .AND. nrhs.GT.0 )
488 $ CALL dlacpy( 'Full', nrows, nrhs,
489 $ copyb, ldb, c, ldb )
490 CALL dqrt16( trans, m, n, nrhs, copya,
491 $ lda, b, ldb, c, ldb, work,
492 $ result( 1 ) )
493*
494* Test 2: Check correctness of results
495* for DGELS.
496*
497 IF( ( itran.EQ.1 .AND. m.GE.n ) .OR.
498 $ ( itran.EQ.2 .AND. m.LT.n ) ) THEN
499*
500* Solving LS system, compute:
501* r = norm((B- A*X)**T * A) /
502* / (norm(A)*norm(B)*max(M,N,NRHS)*EPS)
503*
504 result( 2 ) = dqrt17( trans, 1, m, n,
505 $ nrhs, copya, lda, b, ldb,
506 $ copyb, ldb, c, work,
507 $ lwork )
508 ELSE
509*
510* Solving overdetermined system
511*
512 result( 2 ) = dqrt14( trans, m, n,
513 $ nrhs, copya, lda, b, ldb,
514 $ work, lwork )
515 END IF
516*
517* Print information about the tests that
518* did not pass the threshold.
519*
520 DO k = 1, 2
521 IF( result( k ).GE.thresh ) THEN
522 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
523 $ CALL alahd( nout, path )
524 WRITE( nout, fmt = 9999 ) trans, m,
525 $ n, nrhs, nb, itype, k,
526 $ result( k )
527 nfail = nfail + 1
528 END IF
529 END DO
530 nrun = nrun + 2
531 END DO
532 END DO
533 END IF
534* =====================================================
535* End test DGELS
536* =====================================================
537* =====================================================
538* Begin test DGELST
539* =====================================================
540 IF( irank.EQ.1 ) THEN
541*
542* Generate a matrix of scaling type ISCALE
543*
544 CALL dqrt13( iscale, m, n, copya, lda, norma,
545 $ iseed )
546*
547* Loop for testing different block sizes.
548*
549 DO inb = 1, nnb
550 nb = nbval( inb )
551 CALL xlaenv( 1, nb )
552*
553* Loop for testing non-transposed and transposed.
554*
555 DO itran = 1, 2
556 IF( itran.EQ.1 ) THEN
557 trans = 'N'
558 nrows = m
559 ncols = n
560 ELSE
561 trans = 'T'
562 nrows = n
563 ncols = m
564 END IF
565 ldwork = max( 1, ncols )
566*
567* Set up a consistent rhs
568*
569 IF( ncols.GT.0 ) THEN
570 CALL dlarnv( 2, iseed, ncols*nrhs,
571 $ work )
572 CALL dscal( ncols*nrhs,
573 $ one / dble( ncols ), work,
574 $ 1 )
575 END IF
576 CALL dgemm( trans, 'No transpose', nrows,
577 $ nrhs, ncols, one, copya, lda,
578 $ work, ldwork, zero, b, ldb )
579 CALL dlacpy( 'Full', nrows, nrhs, b, ldb,
580 $ copyb, ldb )
581*
582* Solve LS or overdetermined system
583*
584 IF( m.GT.0 .AND. n.GT.0 ) THEN
585 CALL dlacpy( 'Full', m, n, copya, lda,
586 $ a, lda )
587 CALL dlacpy( 'Full', nrows, nrhs,
588 $ copyb, ldb, b, ldb )
589 END IF
590 srnamt = 'DGELST'
591 CALL dgelst( trans, m, n, nrhs, a, lda, b,
592 $ ldb, work, lwork, info )
593 IF( info.NE.0 )
594 $ CALL alaerh( path, 'DGELST', info, 0,
595 $ trans, m, n, nrhs, -1, nb,
596 $ itype, nfail, nerrs,
597 $ nout )
598*
599* Test 3: Check correctness of results
600* for DGELST, compute the residual:
601* RESID = norm(B - A*X) /
602* / ( max(m,n) * norm(A) * norm(X) * EPS )
603*
604 IF( nrows.GT.0 .AND. nrhs.GT.0 )
605 $ CALL dlacpy( 'Full', nrows, nrhs,
606 $ copyb, ldb, c, ldb )
607 CALL dqrt16( trans, m, n, nrhs, copya,
608 $ lda, b, ldb, c, ldb, work,
609 $ result( 3 ) )
610*
611* Test 4: Check correctness of results
612* for DGELST.
613*
614 IF( ( itran.EQ.1 .AND. m.GE.n ) .OR.
615 $ ( itran.EQ.2 .AND. m.LT.n ) ) THEN
616*
617* Solving LS system, compute:
618* r = norm((B- A*X)**T * A) /
619* / (norm(A)*norm(B)*max(M,N,NRHS)*EPS)
620*
621 result( 4 ) = dqrt17( trans, 1, m, n,
622 $ nrhs, copya, lda, b, ldb,
623 $ copyb, ldb, c, work,
624 $ lwork )
625 ELSE
626*
627* Solving overdetermined system
628*
629 result( 4 ) = dqrt14( trans, m, n,
630 $ nrhs, copya, lda, b, ldb,
631 $ work, lwork )
632 END IF
633*
634* Print information about the tests that
635* did not pass the threshold.
636*
637 DO k = 3, 4
638 IF( result( k ).GE.thresh ) THEN
639 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
640 $ CALL alahd( nout, path )
641 WRITE( nout, fmt = 9999 ) trans, m,
642 $ n, nrhs, nb, itype, k,
643 $ result( k )
644 nfail = nfail + 1
645 END IF
646 END DO
647 nrun = nrun + 2
648 END DO
649 END DO
650 END IF
651* =====================================================
652* End test DGELST
653* =====================================================
654* =====================================================
655* Begin test DGETSLS
656* =====================================================
657 IF( irank.EQ.1 ) THEN
658*
659* Generate a matrix of scaling type ISCALE
660*
661 CALL dqrt13( iscale, m, n, copya, lda, norma,
662 $ iseed )
663*
664* Loop for testing different block sizes MB.
665*
666 DO imb = 1, nnb
667 mb = nbval( imb )
668 CALL xlaenv( 1, mb )
669*
670* Loop for testing different block sizes NB.
671*
672 DO inb = 1, nnb
673 nb = nbval( inb )
674 CALL xlaenv( 2, nb )
675*
676* Loop for testing non-transposed
677* and transposed.
678*
679 DO itran = 1, 2
680 IF( itran.EQ.1 ) THEN
681 trans = 'N'
682 nrows = m
683 ncols = n
684 ELSE
685 trans = 'T'
686 nrows = n
687 ncols = m
688 END IF
689 ldwork = max( 1, ncols )
690*
691* Set up a consistent rhs
692*
693 IF( ncols.GT.0 ) THEN
694 CALL dlarnv( 2, iseed, ncols*nrhs,
695 $ work )
696 CALL dscal( ncols*nrhs,
697 $ one / dble( ncols ),
698 $ work, 1 )
699 END IF
700 CALL dgemm( trans, 'No transpose',
701 $ nrows, nrhs, ncols, one,
702 $ copya, lda, work, ldwork,
703 $ zero, b, ldb )
704 CALL dlacpy( 'Full', nrows, nrhs,
705 $ b, ldb, copyb, ldb )
706*
707* Solve LS or overdetermined system
708*
709 IF( m.GT.0 .AND. n.GT.0 ) THEN
710 CALL dlacpy( 'Full', m, n,
711 $ copya, lda, a, lda )
712 CALL dlacpy( 'Full', nrows, nrhs,
713 $ copyb, ldb, b, ldb )
714 END IF
715 srnamt = 'DGETSLS'
716 CALL dgetsls( trans, m, n, nrhs,
717 $ a, lda, b, ldb, work, lwork,
718 $ info )
719 IF( info.NE.0 )
720 $ CALL alaerh( path, 'DGETSLS', info,
721 $ 0, trans, m, n, nrhs,
722 $ -1, nb, itype, nfail,
723 $ nerrs, nout )
724*
725* Test 5: Check correctness of results
726* for DGETSLS, compute the residual:
727* RESID = norm(B - A*X) /
728* / ( max(m,n) * norm(A) * norm(X) * EPS )
729*
730 IF( nrows.GT.0 .AND. nrhs.GT.0 )
731 $ CALL dlacpy( 'Full', nrows, nrhs,
732 $ copyb, ldb, c, ldb )
733 CALL dqrt16( trans, m, n, nrhs,
734 $ copya, lda, b, ldb,
735 $ c, ldb, work,
736 $ result( 5 ) )
737*
738* Test 6: Check correctness of results
739* for DGETSLS.
740*
741 IF( ( itran.EQ.1 .AND. m.GE.n ) .OR.
742 $ ( itran.EQ.2 .AND. m.LT.n ) ) THEN
743*
744* Solving LS system, compute:
745* r = norm((B- A*X)**T * A) /
746* / (norm(A)*norm(B)*max(M,N,NRHS)*EPS)
747*
748 result( 6 ) = dqrt17( trans, 1, m,
749 $ n, nrhs, copya, lda,
750 $ b, ldb, copyb, ldb,
751 $ c, work, lwork )
752 ELSE
753*
754* Solving overdetermined system
755*
756 result( 6 ) = dqrt14( trans, m, n,
757 $ nrhs, copya, lda,
758 $ b, ldb, work, lwork )
759 END IF
760*
761* Print information about the tests that
762* did not pass the threshold.
763*
764 DO k = 5, 6
765 IF( result( k ).GE.thresh ) THEN
766 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
767 $ CALL alahd( nout, path )
768 WRITE( nout, fmt = 9997 ) trans,
769 $ m, n, nrhs, mb, nb, itype,
770 $ k, result( k )
771 nfail = nfail + 1
772 END IF
773 END DO
774 nrun = nrun + 2
775 END DO
776 END DO
777 END DO
778 END IF
779* =====================================================
780* End test DGETSLS
781* =====================================================
782*
783* Generate a matrix of scaling type ISCALE and rank
784* type IRANK.
785*
786 CALL dqrt15( iscale, irank, m, n, nrhs, copya, lda,
787 $ copyb, ldb, copys, rank, norma, normb,
788 $ iseed, work, lwork )
789*
790* workspace used: MAX(M+MIN(M,N),NRHS*MIN(M,N),2*N+M)
791*
792 ldwork = max( 1, m )
793*
794* Loop for testing different block sizes.
795*
796 DO 100 inb = 1, nnb
797 nb = nbval( inb )
798 CALL xlaenv( 1, nb )
799 CALL xlaenv( 3, nxval( inb ) )
800*
801* Test DGELSY
802*
803* DGELSY: Compute the minimum-norm solution X
804* to min( norm( A * X - B ) )
805* using the rank-revealing orthogonal
806* factorization.
807*
808* Initialize vector IWORK.
809*
810 DO 70 j = 1, n
811 iwork( j ) = 0
812 70 CONTINUE
813*
814 CALL dlacpy( 'Full', m, n, copya, lda, a, lda )
815 CALL dlacpy( 'Full', m, nrhs, copyb, ldb, b,
816 $ ldb )
817*
818 srnamt = 'DGELSY'
819 CALL dgelsy( m, n, nrhs, a, lda, b, ldb, iwork,
820 $ rcond, crank, work, lwlsy, info )
821 IF( info.NE.0 )
822 $ CALL alaerh( path, 'DGELSY', info, 0, ' ', m,
823 $ n, nrhs, -1, nb, itype, nfail,
824 $ nerrs, nout )
825*
826* Test 7: Compute relative error in svd
827* workspace: M*N + 4*MIN(M,N) + MAX(M,N)
828*
829 result( 7 ) = dqrt12( crank, crank, a, lda,
830 $ copys, work, lwork )
831*
832* Test 8: Compute error in solution
833* workspace: M*NRHS + M
834*
835 CALL dlacpy( 'Full', m, nrhs, copyb, ldb, work,
836 $ ldwork )
837 CALL dqrt16( 'No transpose', m, n, nrhs, copya,
838 $ lda, b, ldb, work, ldwork,
839 $ work( m*nrhs+1 ), result( 8 ) )
840*
841* Test 9: Check norm of r'*A
842* workspace: NRHS*(M+N)
843*
844 result( 9 ) = zero
845 IF( m.GT.crank )
846 $ result( 9 ) = dqrt17( 'No transpose', 1, m,
847 $ n, nrhs, copya, lda, b, ldb,
848 $ copyb, ldb, c, work, lwork )
849*
850* Test 10: Check if x is in the rowspace of A
851* workspace: (M+NRHS)*(N+2)
852*
853 result( 10 ) = zero
854*
855 IF( n.GT.crank )
856 $ result( 10 ) = dqrt14( 'No transpose', m, n,
857 $ nrhs, copya, lda, b, ldb,
858 $ work, lwork )
859*
860* Test DGELSS
861*
862* DGELSS: Compute the minimum-norm solution X
863* to min( norm( A * X - B ) )
864* using the SVD.
865*
866 CALL dlacpy( 'Full', m, n, copya, lda, a, lda )
867 CALL dlacpy( 'Full', m, nrhs, copyb, ldb, b,
868 $ ldb )
869 srnamt = 'DGELSS'
870 CALL dgelss( m, n, nrhs, a, lda, b, ldb, s,
871 $ rcond, crank, work, lwork, info )
872 IF( info.NE.0 )
873 $ CALL alaerh( path, 'DGELSS', info, 0, ' ', m,
874 $ n, nrhs, -1, nb, itype, nfail,
875 $ nerrs, nout )
876*
877* workspace used: 3*min(m,n) +
878* max(2*min(m,n),nrhs,max(m,n))
879*
880* Test 11: Compute relative error in svd
881*
882 IF( rank.GT.0 ) THEN
883 CALL daxpy( mnmin, -one, copys, 1, s, 1 )
884 result( 11 ) = dasum( mnmin, s, 1 ) /
885 $ dasum( mnmin, copys, 1 ) /
886 $ ( eps*dble( mnmin ) )
887 ELSE
888 result( 11 ) = zero
889 END IF
890*
891* Test 12: Compute error in solution
892*
893 CALL dlacpy( 'Full', m, nrhs, copyb, ldb, work,
894 $ ldwork )
895 CALL dqrt16( 'No transpose', m, n, nrhs, copya,
896 $ lda, b, ldb, work, ldwork,
897 $ work( m*nrhs+1 ), result( 12 ) )
898*
899* Test 13: Check norm of r'*A
900*
901 result( 13 ) = zero
902 IF( m.GT.crank )
903 $ result( 13 ) = dqrt17( 'No transpose', 1, m,
904 $ n, nrhs, copya, lda, b, ldb,
905 $ copyb, ldb, c, work, lwork )
906*
907* Test 14: Check if x is in the rowspace of A
908*
909 result( 14 ) = zero
910 IF( n.GT.crank )
911 $ result( 14 ) = dqrt14( 'No transpose', m, n,
912 $ nrhs, copya, lda, b, ldb,
913 $ work, lwork )
914*
915* Test DGELSD
916*
917* DGELSD: Compute the minimum-norm solution X
918* to min( norm( A * X - B ) ) using a
919* divide and conquer SVD.
920*
921* Initialize vector IWORK.
922*
923 DO 80 j = 1, n
924 iwork( j ) = 0
925 80 CONTINUE
926*
927 CALL dlacpy( 'Full', m, n, copya, lda, a, lda )
928 CALL dlacpy( 'Full', m, nrhs, copyb, ldb, b,
929 $ ldb )
930*
931 srnamt = 'DGELSD'
932 CALL dgelsd( m, n, nrhs, a, lda, b, ldb, s,
933 $ rcond, crank, work, lwork, iwork,
934 $ info )
935 IF( info.NE.0 )
936 $ CALL alaerh( path, 'DGELSD', info, 0, ' ', m,
937 $ n, nrhs, -1, nb, itype, nfail,
938 $ nerrs, nout )
939*
940* Test 15: Compute relative error in svd
941*
942 IF( rank.GT.0 ) THEN
943 CALL daxpy( mnmin, -one, copys, 1, s, 1 )
944 result( 15 ) = dasum( mnmin, s, 1 ) /
945 $ dasum( mnmin, copys, 1 ) /
946 $ ( eps*dble( mnmin ) )
947 ELSE
948 result( 15 ) = zero
949 END IF
950*
951* Test 16: Compute error in solution
952*
953 CALL dlacpy( 'Full', m, nrhs, copyb, ldb, work,
954 $ ldwork )
955 CALL dqrt16( 'No transpose', m, n, nrhs, copya,
956 $ lda, b, ldb, work, ldwork,
957 $ work( m*nrhs+1 ), result( 16 ) )
958*
959* Test 17: Check norm of r'*A
960*
961 result( 17 ) = zero
962 IF( m.GT.crank )
963 $ result( 17 ) = dqrt17( 'No transpose', 1, m,
964 $ n, nrhs, copya, lda, b, ldb,
965 $ copyb, ldb, c, work, lwork )
966*
967* Test 18: Check if x is in the rowspace of A
968*
969 result( 18 ) = zero
970 IF( n.GT.crank )
971 $ result( 18 ) = dqrt14( 'No transpose', m, n,
972 $ nrhs, copya, lda, b, ldb,
973 $ work, lwork )
974*
975* Print information about the tests that did not
976* pass the threshold.
977*
978 DO 90 k = 7, 18
979 IF( result( k ).GE.thresh ) THEN
980 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
981 $ CALL alahd( nout, path )
982 WRITE( nout, fmt = 9998 )m, n, nrhs, nb,
983 $ itype, k, result( k )
984 nfail = nfail + 1
985 END IF
986 90 CONTINUE
987 nrun = nrun + 12
988*
989 100 CONTINUE
990
991
992
993
994
995
996 110 CONTINUE
997 120 CONTINUE
998 130 CONTINUE
999 140 CONTINUE
1000 150 CONTINUE
1001*
1002* Print a summary of the results.
1003*
1004 CALL alasvm( path, nout, nfail, nrun, nerrs )
1005*
1006 9999 FORMAT( ' TRANS=''', a1, ''', M=', i5, ', N=', i5, ', NRHS=', i4,
1007 $ ', NB=', i4, ', type', i2, ', test(', i2, ')=', g12.5 )
1008 9998 FORMAT( ' M=', i5, ', N=', i5, ', NRHS=', i4, ', NB=', i4,
1009 $ ', type', i2, ', test(', i2, ')=', g12.5 )
1010 9997 FORMAT( ' TRANS=''', a1,' M=', i5, ', N=', i5, ', NRHS=', i4,
1011 $ ', MB=', i4,', NB=', i4,', type', i2,
1012 $ ', test(', i2, ')=', g12.5 )
1013*
1014 DEALLOCATE( work )
1015 DEALLOCATE( iwork )
1016 RETURN
1017*
1018* End of DDRVLS
1019*
1020 END
alasvm
subroutine alasvm(type, nout, nfail, nrun, nerrs)
ALASVM
Definition alasvm.f:73
xlaenv
subroutine xlaenv(ispec, nvalue)
XLAENV
Definition xlaenv.f:81
alaerh
subroutine alaerh(path, subnam, info, infoe, opts, m, n, kl, ku, n5, imat, nfail, nerrs, nout)
ALAERH
Definition alaerh.f:147
alahd
subroutine alahd(iounit, path)
ALAHD
Definition alahd.f:107
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
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
dqrt15
subroutine dqrt15(scale, rksel, m, n, nrhs, a, lda, b, ldb, s, rank, norma, normb, iseed, work, lwork)
DQRT15
Definition dqrt15.f:148
dqrt16
subroutine dqrt16(trans, m, n, nrhs, a, lda, x, ldx, b, ldb, rwork, resid)
DQRT16
Definition dqrt16.f:133
daxpy
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
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
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:203
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
dgelst
subroutine dgelst(trans, m, n, nrhs, a, lda, b, ldb, work, lwork, info)
DGELST solves overdetermined or underdetermined systems for GE matrices using QR or LQ factorization ...
Definition dgelst.f:194
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:206
dgemm
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
dgetsls
subroutine dgetsls(trans, m, n, nrhs, a, lda, b, ldb, work, lwork, info)
DGETSLS
Definition dgetsls.f:162
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
dlarnv
subroutine dlarnv(idist, iseed, n, x)
DLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition dlarnv.f:97
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79