LAPACK 3.3.0

sdrvls.f

Go to the documentation of this file.
00001       SUBROUTINE SDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
00002      $                   NBVAL, NXVAL, THRESH, TSTERR, A, COPYA, B,
00003      $                   COPYB, C, S, COPYS, WORK, IWORK, NOUT )
00004 *
00005 *  -- LAPACK test routine (version 3.1.1) --
00006 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00007 *     January 2007
00008 *
00009 *     .. Scalar Arguments ..
00010       LOGICAL            TSTERR
00011       INTEGER            NM, NN, NNB, NNS, NOUT
00012       REAL               THRESH
00013 *     ..
00014 *     .. Array Arguments ..
00015       LOGICAL            DOTYPE( * )
00016       INTEGER            IWORK( * ), MVAL( * ), NBVAL( * ), NSVAL( * ),
00017      $                   NVAL( * ), NXVAL( * )
00018       REAL               A( * ), B( * ), C( * ), COPYA( * ), COPYB( * ),
00019      $                   COPYS( * ), S( * ), WORK( * )
00020 *     ..
00021 *
00022 *  Purpose
00023 *  =======
00024 *
00025 *  SDRVLS tests the least squares driver routines SGELS, SGELSS, SGELSX,
00026 *  SGELSY and SGELSD.
00027 *
00028 *  Arguments
00029 *  =========
00030 *
00031 *  DOTYPE  (input) LOGICAL array, dimension (NTYPES)
00032 *          The matrix types to be used for testing.  Matrices of type j
00033 *          (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
00034 *          .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
00035 *          The matrix of type j is generated as follows:
00036 *          j=1: A = U*D*V where U and V are random orthogonal matrices
00037 *               and D has random entries (> 0.1) taken from a uniform 
00038 *               distribution (0,1). A is full rank.
00039 *          j=2: The same of 1, but A is scaled up.
00040 *          j=3: The same of 1, but A is scaled down.
00041 *          j=4: A = U*D*V where U and V are random orthogonal matrices
00042 *               and D has 3*min(M,N)/4 random entries (> 0.1) taken
00043 *               from a uniform distribution (0,1) and the remaining
00044 *               entries set to 0. A is rank-deficient. 
00045 *          j=5: The same of 4, but A is scaled up.
00046 *          j=6: The same of 5, but A is scaled down.
00047 *
00048 *  NM      (input) INTEGER
00049 *          The number of values of M contained in the vector MVAL.
00050 *
00051 *  MVAL    (input) INTEGER array, dimension (NM)
00052 *          The values of the matrix row dimension M.
00053 *
00054 *  NN      (input) INTEGER
00055 *          The number of values of N contained in the vector NVAL.
00056 *
00057 *  NVAL    (input) INTEGER array, dimension (NN)
00058 *          The values of the matrix column dimension N.
00059 *
00060 *  NNS     (input) INTEGER
00061 *          The number of values of NRHS contained in the vector NSVAL.
00062 *
00063 *  NSVAL   (input) INTEGER array, dimension (NNS)
00064 *          The values of the number of right hand sides NRHS.
00065 *
00066 *  NNB     (input) INTEGER
00067 *          The number of values of NB and NX contained in the
00068 *          vectors NBVAL and NXVAL.  The blocking parameters are used
00069 *          in pairs (NB,NX).
00070 *
00071 *  NBVAL   (input) INTEGER array, dimension (NNB)
00072 *          The values of the blocksize NB.
00073 *
00074 *  NXVAL   (input) INTEGER array, dimension (NNB)
00075 *          The values of the crossover point NX.
00076 *
00077 *  THRESH  (input) REAL
00078 *          The threshold value for the test ratios.  A result is
00079 *          included in the output file if RESULT >= THRESH.  To have
00080 *          every test ratio printed, use THRESH = 0.
00081 *
00082 *  TSTERR  (input) LOGICAL
00083 *          Flag that indicates whether error exits are to be tested.
00084 *
00085 *  A       (workspace) REAL array, dimension (MMAX*NMAX)
00086 *          where MMAX is the maximum value of M in MVAL and NMAX is the
00087 *          maximum value of N in NVAL.
00088 *
00089 *  COPYA   (workspace) REAL array, dimension (MMAX*NMAX)
00090 *
00091 *  B       (workspace) REAL array, dimension (MMAX*NSMAX)
00092 *          where MMAX is the maximum value of M in MVAL and NSMAX is the
00093 *          maximum value of NRHS in NSVAL.
00094 *
00095 *  COPYB   (workspace) REAL array, dimension (MMAX*NSMAX)
00096 *
00097 *  C       (workspace) REAL array, dimension (MMAX*NSMAX)
00098 *
00099 *  S       (workspace) REAL array, dimension
00100 *                      (min(MMAX,NMAX))
00101 *
00102 *  COPYS   (workspace) REAL array, dimension
00103 *                      (min(MMAX,NMAX))
00104 *
00105 *  WORK    (workspace) REAL array,
00106 *                      dimension (MMAX*NMAX + 4*NMAX + MMAX).
00107 *
00108 *  IWORK   (workspace) INTEGER array, dimension (15*NMAX)
00109 *
00110 *  NOUT    (input) INTEGER
00111 *          The unit number for output.
00112 *
00113 *  =====================================================================
00114 *
00115 *     .. Parameters ..
00116       INTEGER            NTESTS
00117       PARAMETER          ( NTESTS = 18 )
00118       INTEGER            SMLSIZ
00119       PARAMETER          ( SMLSIZ = 25 )
00120       REAL               ONE, TWO, ZERO
00121       PARAMETER          ( ONE = 1.0E0, TWO = 2.0E0, ZERO = 0.0E0 )
00122 *     ..
00123 *     .. Local Scalars ..
00124       CHARACTER          TRANS
00125       CHARACTER*3        PATH
00126       INTEGER            CRANK, I, IM, IN, INB, INFO, INS, IRANK, 
00127      $                   ISCALE, ITRAN, ITYPE, J, K, LDA, LDB, LDWORK, 
00128      $                   LWLSY, LWORK, M, MNMIN, N, NB, NCOLS, NERRS, 
00129      $                   NFAIL, NLVL, NRHS, NROWS, NRUN, RANK
00130       REAL               EPS, NORMA, NORMB, RCOND
00131 *     ..
00132 *     .. Local Arrays ..
00133       INTEGER            ISEED( 4 ), ISEEDY( 4 )
00134       REAL               RESULT( NTESTS )
00135 *     ..
00136 *     .. External Functions ..
00137       REAL               SASUM, SLAMCH, SQRT12, SQRT14, SQRT17
00138       EXTERNAL           SASUM, SLAMCH, SQRT12, SQRT14, SQRT17
00139 *     ..
00140 *     .. External Subroutines ..
00141       EXTERNAL           ALAERH, ALAHD, ALASVM, SAXPY, SERRLS, SGELS,
00142      $                   SGELSD, SGELSS, SGELSX, SGELSY, SGEMM, SLACPY,
00143      $                   SLARNV, SQRT13, SQRT15, SQRT16, SSCAL,
00144      $                   XLAENV
00145 *     ..
00146 *     .. Intrinsic Functions ..
00147       INTRINSIC          INT, LOG, MAX, MIN, REAL, SQRT
00148 *     ..
00149 *     .. Scalars in Common ..
00150       LOGICAL            LERR, OK
00151       CHARACTER*32       SRNAMT
00152       INTEGER            INFOT, IOUNIT
00153 *     ..
00154 *     .. Common blocks ..
00155       COMMON             / INFOC / INFOT, IOUNIT, OK, LERR
00156       COMMON             / SRNAMC / SRNAMT
00157 *     ..
00158 *     .. Data statements ..
00159       DATA               ISEEDY / 1988, 1989, 1990, 1991 /
00160 *     ..
00161 *     .. Executable Statements ..
00162 *
00163 *     Initialize constants and the random number seed.
00164 *
00165       PATH( 1: 1 ) = 'Single precision'
00166       PATH( 2: 3 ) = 'LS'
00167       NRUN = 0
00168       NFAIL = 0
00169       NERRS = 0
00170       DO 10 I = 1, 4
00171          ISEED( I ) = ISEEDY( I )
00172    10 CONTINUE
00173       EPS = SLAMCH( 'Epsilon' )
00174 *
00175 *     Threshold for rank estimation
00176 *
00177       RCOND = SQRT( EPS ) - ( SQRT( EPS )-EPS ) / 2
00178 *
00179 *     Test the error exits
00180 *
00181       CALL XLAENV( 2, 2 )
00182       CALL XLAENV( 9, SMLSIZ )
00183       IF( TSTERR )
00184      $   CALL SERRLS( PATH, NOUT )
00185 *
00186 *     Print the header if NM = 0 or NN = 0 and THRESH = 0.
00187 *
00188       IF( ( NM.EQ.0 .OR. NN.EQ.0 ) .AND. THRESH.EQ.ZERO )
00189      $   CALL ALAHD( NOUT, PATH )
00190       INFOT = 0
00191 *
00192       DO 150 IM = 1, NM
00193          M = MVAL( IM )
00194          LDA = MAX( 1, M )
00195 *
00196          DO 140 IN = 1, NN
00197             N = NVAL( IN )
00198             MNMIN = MIN( M, N )
00199             LDB = MAX( 1, M, N )
00200 *
00201             DO 130 INS = 1, NNS
00202                NRHS = NSVAL( INS )
00203                NLVL = MAX( INT( LOG( MAX( ONE, REAL( MNMIN ) ) /
00204      $                REAL( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1, 0 )
00205                LWORK = MAX( 1, ( M+NRHS )*( N+2 ), ( N+NRHS )*( M+2 ),
00206      $                 M*N+4*MNMIN+MAX( M, N ), 12*MNMIN+2*MNMIN*SMLSIZ+
00207      $                 8*MNMIN*NLVL+MNMIN*NRHS+(SMLSIZ+1)**2 )
00208 *
00209                DO 120 IRANK = 1, 2
00210                   DO 110 ISCALE = 1, 3
00211                      ITYPE = ( IRANK-1 )*3 + ISCALE
00212                      IF( .NOT.DOTYPE( ITYPE ) )
00213      $                  GO TO 110
00214 *
00215                      IF( IRANK.EQ.1 ) THEN
00216 *
00217 *                       Test SGELS
00218 *
00219 *                       Generate a matrix of scaling type ISCALE
00220 *
00221                         CALL SQRT13( ISCALE, M, N, COPYA, LDA, NORMA,
00222      $                               ISEED )
00223                         DO 40 INB = 1, NNB
00224                            NB = NBVAL( INB )
00225                            CALL XLAENV( 1, NB )
00226                            CALL XLAENV( 3, NXVAL( INB ) )
00227 *
00228                            DO 30 ITRAN = 1, 2
00229                               IF( ITRAN.EQ.1 ) THEN
00230                                  TRANS = 'N'
00231                                  NROWS = M
00232                                  NCOLS = N
00233                               ELSE
00234                                  TRANS = 'T'
00235                                  NROWS = N
00236                                  NCOLS = M
00237                               END IF
00238                               LDWORK = MAX( 1, NCOLS )
00239 *
00240 *                             Set up a consistent rhs
00241 *
00242                               IF( NCOLS.GT.0 ) THEN
00243                                  CALL SLARNV( 2, ISEED, NCOLS*NRHS,
00244      $                                        WORK )
00245                                  CALL SSCAL( NCOLS*NRHS,
00246      $                                       ONE / REAL( NCOLS ), WORK,
00247      $                                       1 )
00248                               END IF
00249                               CALL SGEMM( TRANS, 'No transpose', NROWS,
00250      $                                    NRHS, NCOLS, ONE, COPYA, LDA,
00251      $                                    WORK, LDWORK, ZERO, B, LDB )
00252                               CALL SLACPY( 'Full', NROWS, NRHS, B, LDB,
00253      $                                     COPYB, LDB )
00254 *
00255 *                             Solve LS or overdetermined system
00256 *
00257                               IF( M.GT.0 .AND. N.GT.0 ) THEN
00258                                  CALL SLACPY( 'Full', M, N, COPYA, LDA,
00259      $                                        A, LDA )
00260                                  CALL SLACPY( 'Full', NROWS, NRHS,
00261      $                                        COPYB, LDB, B, LDB )
00262                               END IF
00263                               SRNAMT = 'SGELS '
00264                               CALL SGELS( TRANS, M, N, NRHS, A, LDA, B,
00265      $                                    LDB, WORK, LWORK, INFO )
00266                               IF( INFO.NE.0 )
00267      $                           CALL ALAERH( PATH, 'SGELS ', INFO, 0,
00268      $                                        TRANS, M, N, NRHS, -1, NB,
00269      $                                        ITYPE, NFAIL, NERRS,
00270      $                                        NOUT )
00271 *
00272 *                             Check correctness of results
00273 *
00274                               LDWORK = MAX( 1, NROWS )
00275                               IF( NROWS.GT.0 .AND. NRHS.GT.0 )
00276      $                           CALL SLACPY( 'Full', NROWS, NRHS,
00277      $                                        COPYB, LDB, C, LDB )
00278                               CALL SQRT16( TRANS, M, N, NRHS, COPYA,
00279      $                                     LDA, B, LDB, C, LDB, WORK,
00280      $                                     RESULT( 1 ) )
00281 *
00282                               IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR.
00283      $                            ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN
00284 *
00285 *                                Solving LS system
00286 *
00287                                  RESULT( 2 ) = SQRT17( TRANS, 1, M, N,
00288      $                                         NRHS, COPYA, LDA, B, LDB,
00289      $                                         COPYB, LDB, C, WORK,
00290      $                                         LWORK )
00291                               ELSE
00292 *
00293 *                                Solving overdetermined system
00294 *
00295                                  RESULT( 2 ) = SQRT14( TRANS, M, N,
00296      $                                         NRHS, COPYA, LDA, B, LDB,
00297      $                                         WORK, LWORK )
00298                               END IF
00299 *
00300 *                             Print information about the tests that
00301 *                             did not pass the threshold.
00302 *
00303                               DO 20 K = 1, 2
00304                                  IF( RESULT( K ).GE.THRESH ) THEN
00305                                     IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
00306      $                                 CALL ALAHD( NOUT, PATH )
00307                                     WRITE( NOUT, FMT = 9999 )TRANS, M,
00308      $                                 N, NRHS, NB, ITYPE, K,
00309      $                                 RESULT( K )
00310                                     NFAIL = NFAIL + 1
00311                                  END IF
00312    20                         CONTINUE
00313                               NRUN = NRUN + 2
00314    30                      CONTINUE
00315    40                   CONTINUE
00316                      END IF
00317 *
00318 *                    Generate a matrix of scaling type ISCALE and rank
00319 *                    type IRANK.
00320 *
00321                      CALL SQRT15( ISCALE, IRANK, M, N, NRHS, COPYA, LDA,
00322      $                            COPYB, LDB, COPYS, RANK, NORMA, NORMB,
00323      $                            ISEED, WORK, LWORK )
00324 *
00325 *                    workspace used: MAX(M+MIN(M,N),NRHS*MIN(M,N),2*N+M)
00326 *
00327 *                    Initialize vector IWORK.
00328 *
00329                      DO 50 J = 1, N
00330                         IWORK( J ) = 0
00331    50                CONTINUE
00332                      LDWORK = MAX( 1, M )
00333 *
00334 *                    Test SGELSX
00335 *
00336 *                    SGELSX:  Compute the minimum-norm solution X
00337 *                    to min( norm( A * X - B ) ) using a complete
00338 *                    orthogonal factorization.
00339 *
00340                      CALL SLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
00341                      CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, B, LDB )
00342 *
00343                      SRNAMT = 'SGELSX'
00344                      CALL SGELSX( M, N, NRHS, A, LDA, B, LDB, IWORK,
00345      $                            RCOND, CRANK, WORK, INFO )
00346                      IF( INFO.NE.0 )
00347      $                  CALL ALAERH( PATH, 'SGELSX', INFO, 0, ' ', M, N,
00348      $                               NRHS, -1, NB, ITYPE, NFAIL, NERRS,
00349      $                               NOUT )
00350 *
00351 *                    workspace used: MAX( MNMIN+3*N, 2*MNMIN+NRHS )
00352 *
00353 *                    Test 3:  Compute relative error in svd
00354 *                             workspace: M*N + 4*MIN(M,N) + MAX(M,N)
00355 *
00356                      RESULT( 3 ) = SQRT12( CRANK, CRANK, A, LDA, COPYS,
00357      $                             WORK, LWORK )
00358 *
00359 *                    Test 4:  Compute error in solution
00360 *                             workspace:  M*NRHS + M
00361 *
00362                      CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
00363      $                            LDWORK )
00364                      CALL SQRT16( 'No transpose', M, N, NRHS, COPYA,
00365      $                            LDA, B, LDB, WORK, LDWORK,
00366      $                            WORK( M*NRHS+1 ), RESULT( 4 ) )
00367 *
00368 *                    Test 5:  Check norm of r'*A
00369 *                             workspace: NRHS*(M+N)
00370 *
00371                      RESULT( 5 ) = ZERO
00372                      IF( M.GT.CRANK )
00373      $                  RESULT( 5 ) = SQRT17( 'No transpose', 1, M, N,
00374      $                                NRHS, COPYA, LDA, B, LDB, COPYB,
00375      $                                LDB, C, WORK, LWORK )
00376 *
00377 *                    Test 6:  Check if x is in the rowspace of A
00378 *                             workspace: (M+NRHS)*(N+2)
00379 *
00380                      RESULT( 6 ) = ZERO
00381 *
00382                      IF( N.GT.CRANK )
00383      $                  RESULT( 6 ) = SQRT14( 'No transpose', M, N,
00384      $                                NRHS, COPYA, LDA, B, LDB, WORK,
00385      $                                LWORK )
00386 *
00387 *                    Print information about the tests that did not
00388 *                    pass the threshold.
00389 *
00390                      DO 60 K = 3, 6
00391                         IF( RESULT( K ).GE.THRESH ) THEN
00392                            IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
00393      $                        CALL ALAHD( NOUT, PATH )
00394                            WRITE( NOUT, FMT = 9998 )M, N, NRHS, NB,
00395      $                        ITYPE, K, RESULT( K )
00396                            NFAIL = NFAIL + 1
00397                         END IF
00398    60                CONTINUE
00399                      NRUN = NRUN + 4
00400 *
00401 *                    Loop for testing different block sizes.
00402 *
00403                      DO 100 INB = 1, NNB
00404                         NB = NBVAL( INB )
00405                         CALL XLAENV( 1, NB )
00406                         CALL XLAENV( 3, NXVAL( INB ) )
00407 *
00408 *                       Test SGELSY
00409 *
00410 *                       SGELSY:  Compute the minimum-norm solution X
00411 *                       to min( norm( A * X - B ) )
00412 *                       using the rank-revealing orthogonal
00413 *                       factorization.
00414 *
00415 *                       Initialize vector IWORK.
00416 *
00417                         DO 70 J = 1, N
00418                            IWORK( J ) = 0
00419    70                   CONTINUE
00420 *
00421 *                       Set LWLSY to the adequate value.
00422 *
00423                         LWLSY = MAX( 1, MNMIN+2*N+NB*( N+1 ),
00424      $                          2*MNMIN+NB*NRHS )
00425 *
00426                         CALL SLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
00427                         CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, B,
00428      $                               LDB )
00429 *
00430                         SRNAMT = 'SGELSY'
00431                         CALL SGELSY( M, N, NRHS, A, LDA, B, LDB, IWORK,
00432      $                               RCOND, CRANK, WORK, LWLSY, INFO )
00433                         IF( INFO.NE.0 )
00434      $                     CALL ALAERH( PATH, 'SGELSY', INFO, 0, ' ', M,
00435      $                                  N, NRHS, -1, NB, ITYPE, NFAIL,
00436      $                                  NERRS, NOUT )
00437 *
00438 *                       Test 7:  Compute relative error in svd
00439 *                                workspace: M*N + 4*MIN(M,N) + MAX(M,N)
00440 *
00441                         RESULT( 7 ) = SQRT12( CRANK, CRANK, A, LDA,
00442      $                                COPYS, WORK, LWORK )
00443 *
00444 *                       Test 8:  Compute error in solution
00445 *                                workspace:  M*NRHS + M
00446 *
00447                         CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
00448      $                               LDWORK )
00449                         CALL SQRT16( 'No transpose', M, N, NRHS, COPYA,
00450      $                               LDA, B, LDB, WORK, LDWORK,
00451      $                               WORK( M*NRHS+1 ), RESULT( 8 ) )
00452 *
00453 *                       Test 9:  Check norm of r'*A
00454 *                                workspace: NRHS*(M+N)
00455 *
00456                         RESULT( 9 ) = ZERO
00457                         IF( M.GT.CRANK )
00458      $                     RESULT( 9 ) = SQRT17( 'No transpose', 1, M,
00459      $                                   N, NRHS, COPYA, LDA, B, LDB,
00460      $                                   COPYB, LDB, C, WORK, LWORK )
00461 *
00462 *                       Test 10:  Check if x is in the rowspace of A
00463 *                                workspace: (M+NRHS)*(N+2)
00464 *
00465                         RESULT( 10 ) = ZERO
00466 *
00467                         IF( N.GT.CRANK )
00468      $                     RESULT( 10 ) = SQRT14( 'No transpose', M, N,
00469      $                                    NRHS, COPYA, LDA, B, LDB,
00470      $                                    WORK, LWORK )
00471 *
00472 *                       Test SGELSS
00473 *
00474 *                       SGELSS:  Compute the minimum-norm solution X
00475 *                       to min( norm( A * X - B ) )
00476 *                       using the SVD.
00477 *
00478                         CALL SLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
00479                         CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, B,
00480      $                               LDB )
00481                         SRNAMT = 'SGELSS'
00482                         CALL SGELSS( M, N, NRHS, A, LDA, B, LDB, S,
00483      $                               RCOND, CRANK, WORK, LWORK, INFO )
00484                         IF( INFO.NE.0 )
00485      $                     CALL ALAERH( PATH, 'SGELSS', INFO, 0, ' ', M,
00486      $                                  N, NRHS, -1, NB, ITYPE, NFAIL,
00487      $                                  NERRS, NOUT )
00488 *
00489 *                       workspace used: 3*min(m,n) +
00490 *                                       max(2*min(m,n),nrhs,max(m,n))
00491 *
00492 *                       Test 11:  Compute relative error in svd
00493 *
00494                         IF( RANK.GT.0 ) THEN
00495                            CALL SAXPY( MNMIN, -ONE, COPYS, 1, S, 1 )
00496                            RESULT( 11 ) = SASUM( MNMIN, S, 1 ) /
00497      $                                    SASUM( MNMIN, COPYS, 1 ) /
00498      $                                    ( EPS*REAL( MNMIN ) )
00499                         ELSE
00500                            RESULT( 11 ) = ZERO
00501                         END IF
00502 *
00503 *                       Test 12:  Compute error in solution
00504 *
00505                         CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
00506      $                               LDWORK )
00507                         CALL SQRT16( 'No transpose', M, N, NRHS, COPYA,
00508      $                               LDA, B, LDB, WORK, LDWORK,
00509      $                               WORK( M*NRHS+1 ), RESULT( 12 ) )
00510 *
00511 *                       Test 13:  Check norm of r'*A
00512 *
00513                         RESULT( 13 ) = ZERO
00514                         IF( M.GT.CRANK )
00515      $                     RESULT( 13 ) = SQRT17( 'No transpose', 1, M,
00516      $                                    N, NRHS, COPYA, LDA, B, LDB,
00517      $                                    COPYB, LDB, C, WORK, LWORK )
00518 *
00519 *                       Test 14:  Check if x is in the rowspace of A
00520 *
00521                         RESULT( 14 ) = ZERO
00522                         IF( N.GT.CRANK )
00523      $                     RESULT( 14 ) = SQRT14( 'No transpose', M, N,
00524      $                                    NRHS, COPYA, LDA, B, LDB,
00525      $                                    WORK, LWORK )
00526 *
00527 *                       Test SGELSD
00528 *
00529 *                       SGELSD:  Compute the minimum-norm solution X
00530 *                       to min( norm( A * X - B ) ) using a
00531 *                       divide and conquer SVD.
00532 *
00533 *                       Initialize vector IWORK.
00534 *
00535                         DO 80 J = 1, N
00536                            IWORK( J ) = 0
00537    80                   CONTINUE
00538 *
00539                         CALL SLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
00540                         CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, B,
00541      $                               LDB )
00542 *
00543                         SRNAMT = 'SGELSD'
00544                         CALL SGELSD( M, N, NRHS, A, LDA, B, LDB, S,
00545      $                               RCOND, CRANK, WORK, LWORK, IWORK,
00546      $                               INFO )
00547                         IF( INFO.NE.0 )
00548      $                     CALL ALAERH( PATH, 'SGELSD', INFO, 0, ' ', M,
00549      $                                  N, NRHS, -1, NB, ITYPE, NFAIL,
00550      $                                  NERRS, NOUT )
00551 *
00552 *                       Test 15:  Compute relative error in svd
00553 *
00554                         IF( RANK.GT.0 ) THEN
00555                            CALL SAXPY( MNMIN, -ONE, COPYS, 1, S, 1 )
00556                            RESULT( 15 ) = SASUM( MNMIN, S, 1 ) /
00557      $                                    SASUM( MNMIN, COPYS, 1 ) /
00558      $                                    ( EPS*REAL( MNMIN ) )
00559                         ELSE
00560                            RESULT( 15 ) = ZERO
00561                         END IF
00562 *
00563 *                       Test 16:  Compute error in solution
00564 *
00565                         CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
00566      $                               LDWORK )
00567                         CALL SQRT16( 'No transpose', M, N, NRHS, COPYA,
00568      $                               LDA, B, LDB, WORK, LDWORK,
00569      $                               WORK( M*NRHS+1 ), RESULT( 16 ) )
00570 *
00571 *                       Test 17:  Check norm of r'*A
00572 *
00573                         RESULT( 17 ) = ZERO
00574                         IF( M.GT.CRANK )
00575      $                     RESULT( 17 ) = SQRT17( 'No transpose', 1, M,
00576      $                                    N, NRHS, COPYA, LDA, B, LDB,
00577      $                                    COPYB, LDB, C, WORK, LWORK )
00578 *
00579 *                       Test 18:  Check if x is in the rowspace of A
00580 *
00581                         RESULT( 18 ) = ZERO
00582                         IF( N.GT.CRANK )
00583      $                     RESULT( 18 ) = SQRT14( 'No transpose', M, N,
00584      $                                    NRHS, COPYA, LDA, B, LDB,
00585      $                                    WORK, LWORK )
00586 *
00587 *                       Print information about the tests that did not
00588 *                       pass the threshold.
00589 *
00590                         DO 90 K = 7, NTESTS
00591                            IF( RESULT( K ).GE.THRESH ) THEN
00592                               IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
00593      $                           CALL ALAHD( NOUT, PATH )
00594                               WRITE( NOUT, FMT = 9998 )M, N, NRHS, NB,
00595      $                           ITYPE, K, RESULT( K )
00596                               NFAIL = NFAIL + 1
00597                            END IF
00598    90                   CONTINUE
00599                         NRUN = NRUN + 12 
00600 *
00601   100                CONTINUE
00602   110             CONTINUE
00603   120          CONTINUE
00604   130       CONTINUE
00605   140    CONTINUE
00606   150 CONTINUE
00607 *
00608 *     Print a summary of the results.
00609 *
00610       CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS )
00611 *
00612  9999 FORMAT( ' TRANS=''', A1, ''', M=', I5, ', N=', I5, ', NRHS=', I4,
00613      $      ', NB=', I4, ', type', I2, ', test(', I2, ')=', G12.5 )
00614  9998 FORMAT( ' M=', I5, ', N=', I5, ', NRHS=', I4, ', NB=', I4,
00615      $      ', type', I2, ', test(', I2, ')=', G12.5 )
00616       RETURN
00617 *
00618 *     End of SDRVLS
00619 *
00620       END
 All Files Functions