LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ sdrvls()

 subroutine sdrvls ( logical, dimension( * ) dotype, integer nm, integer, dimension( * ) mval, integer nn, integer, dimension( * ) nval, integer nns, integer, dimension( * ) nsval, integer nnb, integer, dimension( * ) nbval, integer, dimension( * ) nxval, real thresh, logical tsterr, real, dimension( * ) a, real, dimension( * ) copya, real, dimension( * ) b, real, dimension( * ) copyb, real, dimension( * ) c, real, dimension( * ) s, real, dimension( * ) copys, integer nout )

SDRVLS

Purpose:
``` SDRVLS tests the least squares driver routines SGELS, SGELST,
SGETSLS, SGELSS, SGELSY and SGELSD.```
Parameters
 [in] DOTYPE ``` DOTYPE is LOGICAL array, dimension (NTYPES) The matrix types to be used for testing. Matrices of type j (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) = .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used. The matrix of type j is generated as follows: j=1: A = U*D*V where U and V are random orthogonal matrices and D has random entries (> 0.1) taken from a uniform distribution (0,1). A is full rank. j=2: The same of 1, but A is scaled up. j=3: The same of 1, but A is scaled down. j=4: A = U*D*V where U and V are random orthogonal matrices and D has 3*min(M,N)/4 random entries (> 0.1) taken from a uniform distribution (0,1) and the remaining entries set to 0. A is rank-deficient. j=5: The same of 4, but A is scaled up. j=6: The same of 5, but A is scaled down.``` [in] NM ``` NM is INTEGER The number of values of M contained in the vector MVAL.``` [in] MVAL ``` MVAL is INTEGER array, dimension (NM) The values of the matrix row dimension M.``` [in] NN ``` NN is INTEGER The number of values of N contained in the vector NVAL.``` [in] NVAL ``` NVAL is INTEGER array, dimension (NN) The values of the matrix column dimension N.``` [in] NNS ``` NNS is INTEGER The number of values of NRHS contained in the vector NSVAL.``` [in] NSVAL ``` NSVAL is INTEGER array, dimension (NNS) The values of the number of right hand sides NRHS.``` [in] NNB ``` NNB is INTEGER The number of values of NB and NX contained in the vectors NBVAL and NXVAL. The blocking parameters are used in pairs (NB,NX).``` [in] NBVAL ``` NBVAL is INTEGER array, dimension (NNB) The values of the blocksize NB.``` [in] NXVAL ``` NXVAL is INTEGER array, dimension (NNB) The values of the crossover point NX.``` [in] THRESH ``` THRESH is REAL The threshold value for the test ratios. A result is included in the output file if RESULT >= THRESH. To have every test ratio printed, use THRESH = 0.``` [in] TSTERR ``` TSTERR is LOGICAL Flag that indicates whether error exits are to be tested.``` [out] A ``` A is REAL array, dimension (MMAX*NMAX) where MMAX is the maximum value of M in MVAL and NMAX is the maximum value of N in NVAL.``` [out] COPYA ` COPYA is REAL array, dimension (MMAX*NMAX)` [out] B ``` B is REAL array, dimension (MMAX*NSMAX) where MMAX is the maximum value of M in MVAL and NSMAX is the maximum value of NRHS in NSVAL.``` [out] COPYB ` COPYB is REAL array, dimension (MMAX*NSMAX)` [out] C ` C is REAL array, dimension (MMAX*NSMAX)` [out] S ``` S is REAL array, dimension (min(MMAX,NMAX))``` [out] COPYS ``` COPYS is REAL array, dimension (min(MMAX,NMAX))``` [in] NOUT ``` NOUT is INTEGER The unit number for output.```

Definition at line 189 of file sdrvls.f.

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,
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*
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
real function sasum(n, sx, incx)
SASUM
Definition sasum.f:72
subroutine saxpy(n, sa, sx, incx, sy, incy)
SAXPY
Definition saxpy.f:89
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
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
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
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
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
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:188
subroutine sgetsls(trans, m, n, nrhs, a, lda, b, ldb, work, lwork, info)
SGETSLS
Definition sgetsls.f:162
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
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
subroutine slarnv(idist, iseed, n, x)
SLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition slarnv.f:97
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
subroutine serrls(path, nunit)
SERRLS
Definition serrls.f:55
real function sqrt12(m, n, a, lda, s, work, lwork)
SQRT12
Definition sqrt12.f:89
subroutine sqrt13(scale, m, n, a, lda, norma, iseed)
SQRT13
Definition sqrt13.f:91
real function sqrt14(trans, m, n, nrhs, a, lda, x, ldx, work, lwork)
SQRT14
Definition sqrt14.f:116
subroutine sqrt15(scale, rksel, m, n, nrhs, a, lda, b, ldb, s, rank, norma, normb, iseed, work, lwork)
SQRT15
Definition sqrt15.f:148
subroutine sqrt16(trans, m, n, nrhs, a, lda, x, ldx, b, ldb, rwork, resid)
SQRT16
Definition sqrt16.f:133
real function sqrt17(trans, iresid, m, n, nrhs, a, lda, x, ldx, b, ldb, c, work, lwork)
SQRT17
Definition sqrt17.f:153
Here is the call graph for this function:
Here is the caller graph for this function: