LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
sdrvls.f
Go to the documentation of this file.
1*> \brief \b SDRVLS
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 SDRVLS( 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* REAL THRESH
19* ..
20* .. Array Arguments ..
21* LOGICAL DOTYPE( * )
22* INTEGER MVAL( * ), NBVAL( * ), NSVAL( * ),
23* $ NVAL( * ), NXVAL( * )
24* REAL A( * ), B( * ), C( * ), COPYA( * ), COPYB( * ),
25* $ COPYS( * ), S( * )
26* ..
27*
28*
29*> \par Purpose:
30* =============
31*>
32*> \verbatim
33*>
34*> SDRVLS tests the least squares driver routines SGELS, SGELST,
35*> SGETSLS, SGELSS, SGELSY and SGELSD.
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 REAL
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 REAL 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 REAL array, dimension (MMAX*NMAX)
141*> \endverbatim
142*>
143*> \param[out] B
144*> \verbatim
145*> B is REAL 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 REAL array, dimension (MMAX*NSMAX)
153*> \endverbatim
154*>
155*> \param[out] C
156*> \verbatim
157*> C is REAL array, dimension (MMAX*NSMAX)
158*> \endverbatim
159*>
160*> \param[out] S
161*> \verbatim
162*> S is REAL array, dimension
163*> (min(MMAX,NMAX))
164*> \endverbatim
165*>
166*> \param[out] COPYS
167*> \verbatim
168*> COPYS is REAL 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 single_lin
187*
188* =====================================================================
189 SUBROUTINE sdrvls( 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 REAL THRESH
201* ..
202* .. Array Arguments ..
203 LOGICAL DOTYPE( * )
204 INTEGER MVAL( * ), NBVAL( * ), NSVAL( * ),
205 $ nval( * ), nxval( * )
206 REAL 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 REAL ONE, TWO, ZERO
218 parameter( one = 1.0e0, two = 2.0e0, zero = 0.0e0 )
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_sgels, lwork_sgelst, lwork_sgetsls,
229 $ lwork_sgelss, lwork_sgelsy, lwork_sgelsd
230 REAL EPS, NORMA, NORMB, RCOND
231* ..
232* .. Local Arrays ..
233 INTEGER ISEED( 4 ), ISEEDY( 4 ), IWQ( 1 )
234 REAL RESULT( NTESTS ), WQ( 1 )
235* ..
236* .. Allocatable Arrays ..
237 REAL, ALLOCATABLE :: WORK (:)
238 INTEGER, ALLOCATABLE :: IWORK (:)
239* ..
240* .. External Functions ..
241 REAL SASUM, SLAMCH, SQRT12, SQRT14, SQRT17
242 EXTERNAL SASUM, SLAMCH, SQRT12, SQRT14, SQRT17
243* ..
244* .. External Subroutines ..
245 EXTERNAL alaerh, alahd, alasvm, saxpy, serrls, sgels,
246 $ sgelsd, sgelss, sgelst, sgelsy, sgemm,
247 $ sgetsls, slacpy, slarnv, sqrt13, sqrt15,
248 $ sqrt16, sscal, xlaenv
249* ..
250* .. Intrinsic Functions ..
251 INTRINSIC int, max, min, real, 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 ) = 'SINGLE 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 = slamch( '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 serrls( 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* SQRT14, SQRT17 (two side cases), SQRT15 and SQRT12
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 SGELS
358 CALL sgels( trans, m, n, nrhs, a, lda,
359 $ b, ldb, wq( 1 ), -1, info )
360 lwork_sgels = int( wq( 1 ) )
361* Compute workspace needed for SGELST
362 CALL sgelst( trans, m, n, nrhs, a, lda,
363 $ b, ldb, wq, -1, info )
364 lwork_sgelst = int( wq( 1 ) )
365* Compute workspace needed for SGETSLS
366 CALL sgetsls( trans, m, n, nrhs, a, lda,
367 $ b, ldb, wq( 1 ), -1, info )
368 lwork_sgetsls = int( wq( 1 ) )
369 ENDDO
370 END IF
371* Compute workspace needed for SGELSY
372 CALL sgelsy( m, n, nrhs, a, lda, b, ldb, iwq,
373 $ rcond, crank, wq, -1, info )
374 lwork_sgelsy = int( wq( 1 ) )
375* Compute workspace needed for SGELSS
376 CALL sgelss( m, n, nrhs, a, lda, b, ldb, s,
377 $ rcond, crank, wq, -1 , info )
378 lwork_sgelss = int( wq( 1 ) )
379* Compute workspace needed for SGELSD
380 CALL sgelsd( m, n, nrhs, a, lda, b, ldb, s,
381 $ rcond, crank, wq, -1, iwq, info )
382 lwork_sgelsd = int( wq( 1 ) )
383* Compute LIWORK workspace needed for SGELSY and SGELSD
384 liwork = max( liwork, n, iwq( 1 ) )
385* Compute LWORK workspace needed for all functions
386 lwork = max( lwork, lwork_sgels, lwork_sgelst,
387 $ lwork_sgetsls, lwork_sgelsy,
388 $ lwork_sgelss, lwork_sgelsd )
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 SGELS
421* =====================================================
422 IF( irank.EQ.1 ) THEN
423*
424* Generate a matrix of scaling type ISCALE
425*
426 CALL sqrt13( 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 slarnv( 2, iseed, ncols*nrhs,
454 $ work )
455 CALL sscal( ncols*nrhs,
456 $ one / real( ncols ), work,
457 $ 1 )
458 END IF
459 CALL sgemm( trans, 'No transpose', nrows,
460 $ nrhs, ncols, one, copya, lda,
461 $ work, ldwork, zero, b, ldb )
462 CALL slacpy( '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 slacpy( 'Full', m, n, copya, lda,
469 $ a, lda )
470 CALL slacpy( 'Full', nrows, nrhs,
471 $ copyb, ldb, b, ldb )
472 END IF
473 srnamt = 'SGELS '
474 CALL sgels( trans, m, n, nrhs, a, lda, b,
475 $ ldb, work, lwork, info )
476 IF( info.NE.0 )
477 $ CALL alaerh( path, 'SGELS ', 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 SGELS, 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 slacpy( 'Full', nrows, nrhs,
489 $ copyb, ldb, c, ldb )
490 CALL sqrt16( 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 SGELS.
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 ) = sqrt17( 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 ) = sqrt14( 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 SGELS
536* =====================================================
537* =====================================================
538* Begin test SGELST
539* =====================================================
540 IF( irank.EQ.1 ) THEN
541*
542* Generate a matrix of scaling type ISCALE
543*
544 CALL sqrt13( 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 slarnv( 2, iseed, ncols*nrhs,
571 $ work )
572 CALL sscal( ncols*nrhs,
573 $ one / real( ncols ), work,
574 $ 1 )
575 END IF
576 CALL sgemm( trans, 'No transpose', nrows,
577 $ nrhs, ncols, one, copya, lda,
578 $ work, ldwork, zero, b, ldb )
579 CALL slacpy( '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 slacpy( 'Full', m, n, copya, lda,
586 $ a, lda )
587 CALL slacpy( 'Full', nrows, nrhs,
588 $ copyb, ldb, b, ldb )
589 END IF
590 srnamt = 'SGELST'
591 CALL sgelst( trans, m, n, nrhs, a, lda, b,
592 $ ldb, work, lwork, info )
593 IF( info.NE.0 )
594 $ CALL alaerh( path, 'SGELST', 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 SGELST, 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 slacpy( 'Full', nrows, nrhs,
606 $ copyb, ldb, c, ldb )
607 CALL sqrt16( 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 SGELST.
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 ) = sqrt17( 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 ) = sqrt14( 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 SGELST
653* =====================================================
654* =====================================================
655* Begin test SGETSLS
656* =====================================================
657 IF( irank.EQ.1 ) THEN
658*
659* Generate a matrix of scaling type ISCALE
660*
661 CALL sqrt13( 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 slarnv( 2, iseed, ncols*nrhs,
695 $ work )
696 CALL sscal( ncols*nrhs,
697 $ one / real( ncols ),
698 $ work, 1 )
699 END IF
700 CALL sgemm( trans, 'No transpose',
701 $ nrows, nrhs, ncols, one,
702 $ copya, lda, work, ldwork,
703 $ zero, b, ldb )
704 CALL slacpy( '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 slacpy( 'Full', m, n,
711 $ copya, lda, a, lda )
712 CALL slacpy( 'Full', nrows, nrhs,
713 $ copyb, ldb, b, ldb )
714 END IF
715 srnamt = 'SGETSLS'
716 CALL sgetsls( trans, m, n, nrhs,
717 $ a, lda, b, ldb, work, lwork,
718 $ info )
719 IF( info.NE.0 )
720 $ CALL alaerh( path, 'SGETSLS', 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 SGETSLS, 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 slacpy( 'Full', nrows, nrhs,
732 $ copyb, ldb, c, ldb )
733 CALL sqrt16( 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 SGETSLS.
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 ) = sqrt17( 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 ) = sqrt14( 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 SGETSLS
781* =====================================================
782*
783* Generate a matrix of scaling type ISCALE and rank
784* type IRANK.
785*
786 CALL sqrt15( 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 SGELSY
802*
803* SGELSY: 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 slacpy( 'Full', m, n, copya, lda, a, lda )
815 CALL slacpy( 'Full', m, nrhs, copyb, ldb, b,
816 $ ldb )
817*
818 srnamt = 'SGELSY'
819 CALL sgelsy( m, n, nrhs, a, lda, b, ldb, iwork,
820 $ rcond, crank, work, lwlsy, info )
821 IF( info.NE.0 )
822 $ CALL alaerh( path, 'SGELSY', 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 ) = sqrt12( crank, crank, a, lda,
830 $ copys, work, lwork )
831*
832* Test 8: Compute error in solution
833* workspace: M*NRHS + M
834*
835 CALL slacpy( 'Full', m, nrhs, copyb, ldb, work,
836 $ ldwork )
837 CALL sqrt16( '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 ) = sqrt17( '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 ) = sqrt14( 'No transpose', m, n,
857 $ nrhs, copya, lda, b, ldb,
858 $ work, lwork )
859*
860* Test SGELSS
861*
862* SGELSS: Compute the minimum-norm solution X
863* to min( norm( A * X - B ) )
864* using the SVD.
865*
866 CALL slacpy( 'Full', m, n, copya, lda, a, lda )
867 CALL slacpy( 'Full', m, nrhs, copyb, ldb, b,
868 $ ldb )
869 srnamt = 'SGELSS'
870 CALL sgelss( m, n, nrhs, a, lda, b, ldb, s,
871 $ rcond, crank, work, lwork, info )
872 IF( info.NE.0 )
873 $ CALL alaerh( path, 'SGELSS', 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 saxpy( mnmin, -one, copys, 1, s, 1 )
884 result( 11 ) = sasum( mnmin, s, 1 ) /
885 $ sasum( mnmin, copys, 1 ) /
886 $ ( eps*real( mnmin ) )
887 ELSE
888 result( 11 ) = zero
889 END IF
890*
891* Test 12: Compute error in solution
892*
893 CALL slacpy( 'Full', m, nrhs, copyb, ldb, work,
894 $ ldwork )
895 CALL sqrt16( '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 ) = sqrt17( '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 ) = sqrt14( 'No transpose', m, n,
912 $ nrhs, copya, lda, b, ldb,
913 $ work, lwork )
914*
915* Test SGELSD
916*
917* SGELSD: 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 slacpy( 'Full', m, n, copya, lda, a, lda )
928 CALL slacpy( 'Full', m, nrhs, copyb, ldb, b,
929 $ ldb )
930*
931 srnamt = 'SGELSD'
932 CALL sgelsd( 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, 'SGELSD', 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 saxpy( mnmin, -one, copys, 1, s, 1 )
944 result( 15 ) = sasum( mnmin, s, 1 ) /
945 $ sasum( mnmin, copys, 1 ) /
946 $ ( eps*real( mnmin ) )
947 ELSE
948 result( 15 ) = zero
949 END IF
950*
951* Test 16: Compute error in solution
952*
953 CALL slacpy( 'Full', m, nrhs, copyb, ldb, work,
954 $ ldwork )
955 CALL sqrt16( '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 ) = sqrt17( '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 ) = sqrt14( '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 110 CONTINUE
991 120 CONTINUE
992 130 CONTINUE
993 140 CONTINUE
994 150 CONTINUE
995*
996* Print a summary of the results.
997*
998 CALL alasvm( path, nout, nfail, nrun, nerrs )
999*
1000 9999 FORMAT( ' TRANS=''', a1, ''', M=', i5, ', N=', i5, ', NRHS=', i4,
1001 $ ', NB=', i4, ', type', i2, ', test(', i2, ')=', g12.5 )
1002 9998 FORMAT( ' M=', i5, ', N=', i5, ', NRHS=', i4, ', NB=', i4,
1003 $ ', type', i2, ', test(', i2, ')=', g12.5 )
1004 9997 FORMAT( ' TRANS=''', a1,' M=', i5, ', N=', i5, ', NRHS=', i4,
1005 $ ', MB=', i4,', NB=', i4,', type', i2,
1006 $ ', test(', i2, ')=', g12.5 )
1007*
1008 DEALLOCATE( work )
1009 DEALLOCATE( iwork )
1010 RETURN
1011*
1012* End of SDRVLS
1013*
1014 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
saxpy
subroutine saxpy(n, sa, sx, incx, sy, incy)
SAXPY
Definition saxpy.f:89
sgels
subroutine sgels(trans, m, n, nrhs, a, lda, b, ldb, work, lwork, info)
SGELS solves overdetermined or underdetermined systems for GE matrices
Definition sgels.f:183
sgelsd
subroutine sgelsd(m, n, nrhs, a, lda, b, ldb, s, rcond, rank, work, lwork, iwork, info)
SGELSD computes the minimum-norm solution to a linear least squares problem for GE matrices
Definition sgelsd.f:204
sgelss
subroutine sgelss(m, n, nrhs, a, lda, b, ldb, s, rcond, rank, work, lwork, info)
SGELSS solves overdetermined or underdetermined systems for GE matrices
Definition sgelss.f:172
sgelst
subroutine sgelst(trans, m, n, nrhs, a, lda, b, ldb, work, lwork, info)
SGELST solves overdetermined or underdetermined systems for GE matrices using QR or LQ factorization ...
Definition sgelst.f:194
sgelsy
subroutine sgelsy(m, n, nrhs, a, lda, b, ldb, jpvt, rcond, rank, work, lwork, info)
SGELSY solves overdetermined or underdetermined systems for GE matrices
Definition sgelsy.f:206
sgemm
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:188
sgetsls
subroutine sgetsls(trans, m, n, nrhs, a, lda, b, ldb, work, lwork, info)
SGETSLS
Definition sgetsls.f:162
slacpy
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
Definition slacpy.f:103
slarnv
subroutine slarnv(idist, iseed, n, x)
SLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition slarnv.f:97
sscal
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
sdrvls
subroutine sdrvls(dotype, nm, mval, nn, nval, nns, nsval, nnb, nbval, nxval, thresh, tsterr, a, copya, b, copyb, c, s, copys, nout)
SDRVLS
Definition sdrvls.f:192
serrls
subroutine serrls(path, nunit)
SERRLS
Definition serrls.f:55
sqrt13
subroutine sqrt13(scale, m, n, a, lda, norma, iseed)
SQRT13
Definition sqrt13.f:91
sqrt15
subroutine sqrt15(scale, rksel, m, n, nrhs, a, lda, b, ldb, s, rank, norma, normb, iseed, work, lwork)
SQRT15
Definition sqrt15.f:148
sqrt16
subroutine sqrt16(trans, m, n, nrhs, a, lda, x, ldx, b, ldb, rwork, resid)
SQRT16
Definition sqrt16.f:133