LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ 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, 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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

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 = 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,
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  IF( info.NE.0 )
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  IF( rank.GT.0 ) 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  IF( m.GT.crank )
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  IF( n.GT.crank )
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  IF( result( k ).GE.thresh ) THEN
819  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
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 *
subroutine slarnv(IDIST, ISEED, N, X)
SLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition: slarnv.f:97
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
subroutine alasvm(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASVM
Definition: alasvm.f:73
subroutine xlaenv(ISPEC, NVALUE)
XLAENV
Definition: xlaenv.f:81
subroutine alahd(IOUNIT, PATH)
ALAHD
Definition: alahd.f:107
subroutine alaerh(PATH, SUBNAM, INFO, INFOE, OPTS, M, N, KL, KU, N5, IMAT, NFAIL, NERRS, NOUT)
ALAERH
Definition: alaerh.f:147
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 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 sgetsls(TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK, INFO)
SGETSLS
Definition: sgetsls.f:162
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
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
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
Definition: saxpy.f:89
real function sasum(N, SX, INCX)
SASUM
Definition: sasum.f:72
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:187
subroutine serrls(PATH, NUNIT)
SERRLS
Definition: serrls.f:55
real function sqrt17(TRANS, IRESID, M, N, NRHS, A, LDA, X, LDX, B, LDB, C, WORK, LWORK)
SQRT17
Definition: sqrt17.f:153
real function sqrt14(TRANS, M, N, NRHS, A, LDA, X, LDX, WORK, LWORK)
SQRT14
Definition: sqrt14.f:116
subroutine sqrt13(SCALE, M, N, A, LDA, NORMA, ISEED)
SQRT13
Definition: sqrt13.f:91
subroutine sqrt15(SCALE, RKSEL, M, N, NRHS, A, LDA, B, LDB, S, RANK, NORMA, NORMB, ISEED, WORK, LWORK)
SQRT15
Definition: sqrt15.f:148
real function sqrt12(M, N, A, LDA, S, WORK, LWORK)
SQRT12
Definition: sqrt12.f:89
subroutine sqrt16(TRANS, M, N, NRHS, A, LDA, X, LDX, B, LDB, RWORK, RESID)
SQRT16
Definition: sqrt16.f:133
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: