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,
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
subroutine alasvm(type, nout, nfail, nrun, nerrs)
ALASVM
Definition alasvm.f:73
subroutine xlaenv(ispec, nvalue)
XLAENV
Definition xlaenv.f:81
subroutine alaerh(path, subnam, info, infoe, opts, m, n, kl, ku, n5, imat, nfail, nerrs, nout)
ALAERH
Definition alaerh.f:147
subroutine alahd(iounit, path)
ALAHD
Definition alahd.f:107
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
subroutine derrls(path, nunit)
DERRLS
Definition derrls.f:55
subroutine dqrt13(scale, m, n, a, lda, norma, iseed)
DQRT13
Definition dqrt13.f:91
subroutine dqrt15(scale, rksel, m, n, nrhs, a, lda, b, ldb, s, rank, norma, normb, iseed, work, lwork)
DQRT15
Definition dqrt15.f:148
subroutine dqrt16(trans, m, n, nrhs, a, lda, x, ldx, b, ldb, rwork, resid)
DQRT16
Definition dqrt16.f:133
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
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
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
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
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
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
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
subroutine dgetsls(trans, m, n, nrhs, a, lda, b, ldb, work, lwork, info)
DGETSLS
Definition dgetsls.f:162
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
subroutine dlarnv(idist, iseed, n, x)
DLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition dlarnv.f:97
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79