OpenRadioss 2025.1.11
OpenRadioss project
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, SGETSLS, SGELSS, SGELSY,
35*> 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 = 16 )
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_sgetsls, lwork_sgelss,
229 $ 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, sgelsy, sgemm, slacpy,
247 $ slarnv, sqrt13, sqrt15, sqrt16, sscal,
248 $ xlaenv, sgetsls
249* ..
250* .. Intrinsic Functions ..
251 INTRINSIC int, log, 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, ?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 SGELS
357 CALL sgels( trans, m, n, nrhs, a, lda,
358 $ b, ldb, wq( 1 ), -1, info )
359 lwork_sgels = int( wq( 1 ) )
360* Compute workspace needed for SGETSLS
361 CALL sgetsls( trans, m, n, nrhs, a, lda,
362 $ b, ldb, wq( 1 ), -1, info )
363 lwork_sgetsls = int( wq( 1 ) )
364 ENDDO
365 END IF
366* Compute workspace needed for SGELSY
367 CALL sgelsy( m, n, nrhs, a, lda, b, ldb, iwq,
368 $ rcond, crank, wq, -1, info )
369 lwork_sgelsy = int( wq( 1 ) )
370* Compute workspace needed for SGELSS
371 CALL sgelss( m, n, nrhs, a, lda, b, ldb, s,
372 $ rcond, crank, wq, -1 , info )
373 lwork_sgelss = int( wq( 1 ) )
374* Compute workspace needed for SGELSD
375 CALL sgelsd( m, n, nrhs, a, lda, b, ldb, s,
376 $ rcond, crank, wq, -1, iwq, info )
377 lwork_sgelsd = int( wq( 1 ) )
378* Compute LIWORK workspace needed for SGELSY and SGELSD
379 liwork = max( liwork, n, iwq( 1 ) )
380* Compute LWORK workspace needed for all functions
381 lwork = max( lwork, lwork_sgels, lwork_sgetsls,
382 $ lwork_sgelsy, lwork_sgelss,
383 $ lwork_sgelsd )
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 SGELS
418*
419* Generate a matrix of scaling type ISCALE
420*
421 CALL sqrt13( 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 slarnv( 2, iseed, ncols*nrhs,
444 $ work )
445 CALL sscal( ncols*nrhs,
446 $ one / real( ncols ), work,
447 $ 1 )
448 END IF
449 CALL sgemm( trans, 'No transpose', nrows,
450 $ nrhs, ncols, one, copya, lda,
451 $ work, ldwork, zero, b, ldb )
452 CALL slacpy( '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 slacpy( 'Full', m, n, copya, lda,
459 $ a, lda )
460 CALL slacpy( 'Full', nrows, nrhs,
461 $ copyb, ldb, b, ldb )
462 END IF
463 srnamt = 'SGELS '
464 CALL sgels( trans, m, n, nrhs, a, lda, b,
465 $ ldb, work, lwork, info )
466 IF( info.NE.0 )
467 $ CALL alaerh( path, 'SGELS ', 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 slacpy( 'Full', nrows, nrhs,
477 $ copyb, ldb, c, ldb )
478 CALL sqrt16( 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 ) = sqrt17( 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 ) = sqrt14( 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 SGETSLS
519*
520* Generate a matrix of scaling type ISCALE
521*
522 CALL sqrt13( 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 slarnv( 2, iseed, ncols*nrhs,
547 $ work )
548 CALL sscal( ncols*nrhs,
549 $ one / real( ncols ), work,
550 $ 1 )
551 END IF
552 CALL sgemm( trans, 'No transpose', nrows,
553 $ nrhs, ncols, one, copya, lda,
554 $ work, ldwork, zero, b, ldb )
555 CALL slacpy( '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 slacpy( 'Full', m, n, copya, lda,
562 $ a, lda )
563 CALL slacpy( 'Full', nrows, nrhs,
564 $ copyb, ldb, b, ldb )
565 END IF
566 srnamt = 'SGETSLS '
567 CALL sgetsls( trans, m, n, nrhs, a,
568 $ lda, b, ldb, work, lwork, info )
569 IF( info.NE.0 )
570 $ CALL alaerh( path, 'SGETSLS ', 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 slacpy( 'Full', nrows, nrhs,
580 $ copyb, ldb, c, ldb )
581 CALL sqrt16( trans, m, n, nrhs, copya,
582 $ lda, b, ldb, c, ldb, work,
583 $ result( 15 ) )
584*
585 IF( ( itran.EQ.1 .AND. m.GE.n ) .OR.
586 $ ( itran.EQ.2 .AND. m.LT.n ) ) THEN
587*
588* Solving LS system
589*
590 result( 16 ) = sqrt17( 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 ) = sqrt14( 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 IF( result( k ).GE.thresh ) THEN
608 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
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 sqrt15( 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 SGELSY
641*
642* SGELSY: 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 slacpy( 'Full', m, n, copya, lda, a, lda )
654 CALL slacpy( 'Full', m, nrhs, copyb, ldb, b,
655 $ ldb )
656*
657 srnamt = 'SGELSY'
658 CALL sgelsy( m, n, nrhs, a, lda, b, ldb, iwork,
659 $ rcond, crank, work, lwlsy, info )
660 IF( info.NE.0 )
661 $ CALL alaerh( path, 'SGELSY', 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 ) = sqrt12( crank, crank, a, lda,
669 $ copys, work, lwork )
670*
671* Test 4: Compute error in solution
672* workspace: M*NRHS + M
673*
674 CALL slacpy( 'Full', m, nrhs, copyb, ldb, work,
675 $ ldwork )
676 CALL sqrt16( '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 IF( m.GT.crank )
685 $ result( 5 ) = sqrt17( '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 IF( n.GT.crank )
695 $ result( 6 ) = sqrt14( 'No transpose', m, n,
696 $ nrhs, copya, lda, b, ldb,
697 $ work, lwork )
698*
699* Test SGELSS
700*
701* SGELSS: Compute the minimum-norm solution X
702* to min( norm( A * X - B ) )
703* using the SVD.
704*
705 CALL slacpy( 'Full', m, n, copya, lda, a, lda )
706 CALL slacpy( 'Full', m, nrhs, copyb, ldb, b,
707 $ ldb )
708 srnamt = 'SGELSS'
709 CALL sgelss( m, n, nrhs, a, lda, b, ldb, s,
710 $ rcond, crank, work, lwork, info )
711 IF( info.NE.0 )
712 $ CALL alaerh( path, 'SGELSS', 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 IF( rank.GT.0 ) THEN
722 CALL saxpy( mnmin, -one, copys, 1, s, 1 )
723 result( 7 ) = sasum( mnmin, s, 1 ) /
724 $ sasum( mnmin, copys, 1 ) /
725 $ ( eps*real( mnmin ) )
726 ELSE
727 result( 7 ) = zero
728 END IF
729*
730* Test 8: Compute error in solution
731*
732 CALL slacpy( 'Full', m, nrhs, copyb, ldb, work,
733 $ ldwork )
734 CALL sqrt16( '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 IF( m.GT.crank )
742 $ result( 9 ) = sqrt17( '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 IF( n.GT.crank )
750 $ result( 10 ) = sqrt14( 'No transpose', m, n,
751 $ nrhs, copya, lda, b, ldb,
752 $ work, lwork )
753*
754* Test SGELSD
755*
756* SGELSD: 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 slacpy( 'Full', m, n, copya, lda, a, lda )
767 CALL slacpy( 'Full', m, nrhs, copyb, ldb, b,
768 $ ldb )
769*
770 srnamt = 'sgelsd'
771 CALL SGELSD( M, N, NRHS, A, LDA, B, LDB, S,
772 $ RCOND, CRANK, WORK, LWORK, IWORK,
773 $ INFO )
774.NE. IF( INFO0 )
775 $ CALL ALAERH( PATH, 'sgelsd', 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 SAXPY( MNMIN, -ONE, COPYS, 1, S, 1 )
783 RESULT( 11 ) = SASUM( MNMIN, S, 1 ) /
784 $ SASUM( MNMIN, COPYS, 1 ) /
785 $ ( EPS*REAL( MNMIN ) )
786 ELSE
787 RESULT( 11 ) = ZERO
788 END IF
789*
790* Test 12: Compute error in solution
791*
792 CALL SLACPY( 'full', M, NRHS, COPYB, LDB, WORK,
793 $ LDWORK )
794 CALL SQRT16( '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 ) = SQRT17( '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 ) = SQRT14( '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 SDRVLS
852*
853 END
slarnv
subroutine slarnv(idist, iseed, n, x)
SLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition slarnv.f:97
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
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
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:210
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
sgetsls
subroutine sgetsls(trans, m, n, nrhs, a, lda, b, ldb, work, lwork, info)
SGETSLS
Definition sgetsls.f:162
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:204
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
sscal
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
saxpy
subroutine saxpy(n, sa, sx, incx, sy, incy)
SAXPY
Definition saxpy.f:89
sgemm
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:187
sqrt16
subroutine sqrt16(trans, m, n, nrhs, a, lda, x, ldx, b, ldb, rwork, resid)
SQRT16
Definition sqrt16.f:133
sqrt15
subroutine sqrt15(scale, rksel, m, n, nrhs, a, lda, b, ldb, s, rank, norma, normb, iseed, work, lwork)
SQRT15
Definition sqrt15.f:148
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
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21