LAPACK: TESTING/LIN/ddrvls.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DDRVLS( 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       DOUBLE PRECISION   THRESH
00013 *     ..
00014 *     .. Array Arguments ..
00015       LOGICAL            DOTYPE( * )
00016       INTEGER            IWORK( * ), MVAL( * ), NBVAL( * ), NSVAL( * ),
00017      $                   NVAL( * ), NXVAL( * )
00018       DOUBLE PRECISION   A( * ), B( * ), C( * ), COPYA( * ), COPYB( * ),
00019      $                   COPYS( * ), S( * ), WORK( * )
00020 *     ..
00021 *
00022 *  Purpose
00023 *  =======
00024 *
00025 *  DDRVLS tests the least squares driver routines DGELS, DGELSS, DGELSX,
00026 *  DGELSY and DGELSD.
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) DOUBLE PRECISION
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) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension (MMAX*NMAX)
00090 *
00091 *  B       (workspace) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension (MMAX*NSMAX)
00096 *
00097 *  C       (workspace) DOUBLE PRECISION array, dimension (MMAX*NSMAX)
00098 *
00099 *  S       (workspace) DOUBLE PRECISION array, dimension
00100 *                      (min(MMAX,NMAX))
00101 *
00102 *  COPYS   (workspace) DOUBLE PRECISION array, dimension
00103 *                      (min(MMAX,NMAX))
00104 *
00105 *  WORK    (workspace) DOUBLE PRECISION 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       DOUBLE PRECISION   ONE, TWO, ZERO
00121       PARAMETER          ( ONE = 1.0D0, TWO = 2.0D0, ZERO = 0.0D0 )
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       DOUBLE PRECISION   EPS, NORMA, NORMB, RCOND
00131 *     ..
00132 *     .. Local Arrays ..
00133       INTEGER            ISEED( 4 ), ISEEDY( 4 )
00134       DOUBLE PRECISION   RESULT( NTESTS )
00135 *     ..
00136 *     .. External Functions ..
00137       DOUBLE PRECISION   DASUM, DLAMCH, DQRT12, DQRT14, DQRT17
00138       EXTERNAL           DASUM, DLAMCH, DQRT12, DQRT14, DQRT17
00139 *     ..
00140 *     .. External Subroutines ..
00141       EXTERNAL           ALAERH, ALAHD, ALASVM, DAXPY, DERRLS, DGELS,
00142      $                   DGELSD, DGELSS, DGELSX, DGELSY, DGEMM, DLACPY,
00143      $                   DLARNV, DLASRT, DQRT13, DQRT15, DQRT16, DSCAL,
00144      $                   XLAENV
00145 *     ..
00146 *     .. Intrinsic Functions ..
00147       INTRINSIC          DBLE, INT, LOG, MAX, MIN, 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 ) = 'Double 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 = DLAMCH( '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 DERRLS( 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       CALL XLAENV( 2, 2 )
00192       CALL XLAENV( 9, SMLSIZ )
00193 *
00194       DO 150 IM = 1, NM
00195          M = MVAL( IM )
00196          LDA = MAX( 1, M )
00197 *
00198          DO 140 IN = 1, NN
00199             N = NVAL( IN )
00200             MNMIN = MIN( M, N )
00201             LDB = MAX( 1, M, N )
00202 *
00203             DO 130 INS = 1, NNS
00204                NRHS = NSVAL( INS )
00205                NLVL = MAX( INT( LOG( MAX( ONE, DBLE( MNMIN ) ) /
00206      $                DBLE( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1, 0 )
00207                LWORK = MAX( 1, ( M+NRHS )*( N+2 ), ( N+NRHS )*( M+2 ),
00208      $                 M*N+4*MNMIN+MAX( M, N ), 12*MNMIN+2*MNMIN*SMLSIZ+
00209      $                 8*MNMIN*NLVL+MNMIN*NRHS+(SMLSIZ+1)**2 )
00210 *
00211                DO 120 IRANK = 1, 2
00212                   DO 110 ISCALE = 1, 3
00213                      ITYPE = ( IRANK-1 )*3 + ISCALE
00214                      IF( .NOT.DOTYPE( ITYPE ) )
00215      $                  GO TO 110
00216 *
00217                      IF( IRANK.EQ.1 ) THEN
00218 *
00219 *                       Test DGELS
00220 *
00221 *                       Generate a matrix of scaling type ISCALE
00222 *
00223                         CALL DQRT13( ISCALE, M, N, COPYA, LDA, NORMA,
00224      $                               ISEED )
00225                         DO 40 INB = 1, NNB
00226                            NB = NBVAL( INB )
00227                            CALL XLAENV( 1, NB )
00228                            CALL XLAENV( 3, NXVAL( INB ) )
00229 *
00230                            DO 30 ITRAN = 1, 2
00231                               IF( ITRAN.EQ.1 ) THEN
00232                                  TRANS = 'N'
00233                                  NROWS = M
00234                                  NCOLS = N
00235                               ELSE
00236                                  TRANS = 'T'
00237                                  NROWS = N
00238                                  NCOLS = M
00239                               END IF
00240                               LDWORK = MAX( 1, NCOLS )
00241 *
00242 *                             Set up a consistent rhs
00243 *
00244                               IF( NCOLS.GT.0 ) THEN
00245                                  CALL DLARNV( 2, ISEED, NCOLS*NRHS,
00246      $                                        WORK )
00247                                  CALL DSCAL( NCOLS*NRHS,
00248      $                                       ONE / DBLE( NCOLS ), WORK,
00249      $                                       1 )
00250                               END IF
00251                               CALL DGEMM( TRANS, 'No transpose', NROWS,
00252      $                                    NRHS, NCOLS, ONE, COPYA, LDA,
00253      $                                    WORK, LDWORK, ZERO, B, LDB )
00254                               CALL DLACPY( 'Full', NROWS, NRHS, B, LDB,
00255      $                                     COPYB, LDB )
00256 *
00257 *                             Solve LS or overdetermined system
00258 *
00259                               IF( M.GT.0 .AND. N.GT.0 ) THEN
00260                                  CALL DLACPY( 'Full', M, N, COPYA, LDA,
00261      $                                        A, LDA )
00262                                  CALL DLACPY( 'Full', NROWS, NRHS,
00263      $                                        COPYB, LDB, B, LDB )
00264                               END IF
00265                               SRNAMT = 'DGELS '
00266                               CALL DGELS( TRANS, M, N, NRHS, A, LDA, B,
00267      $                                    LDB, WORK, LWORK, INFO )
00268                               IF( INFO.NE.0 )
00269      $                           CALL ALAERH( PATH, 'DGELS ', INFO, 0,
00270      $                                        TRANS, M, N, NRHS, -1, NB,
00271      $                                        ITYPE, NFAIL, NERRS,
00272      $                                        NOUT )
00273 *
00274 *                             Check correctness of results
00275 *
00276                               LDWORK = MAX( 1, NROWS )
00277                               IF( NROWS.GT.0 .AND. NRHS.GT.0 )
00278      $                           CALL DLACPY( 'Full', NROWS, NRHS,
00279      $                                        COPYB, LDB, C, LDB )
00280                               CALL DQRT16( TRANS, M, N, NRHS, COPYA,
00281      $                                     LDA, B, LDB, C, LDB, WORK,
00282      $                                     RESULT( 1 ) )
00283 *
00284                               IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR.
00285      $                            ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN
00286 *
00287 *                                Solving LS system
00288 *
00289                                  RESULT( 2 ) = DQRT17( TRANS, 1, M, N,
00290      $                                         NRHS, COPYA, LDA, B, LDB,
00291      $                                         COPYB, LDB, C, WORK,
00292      $                                         LWORK )
00293                               ELSE
00294 *
00295 *                                Solving overdetermined system
00296 *
00297                                  RESULT( 2 ) = DQRT14( TRANS, M, N,
00298      $                                         NRHS, COPYA, LDA, B, LDB,
00299      $                                         WORK, LWORK )
00300                               END IF
00301 *
00302 *                             Print information about the tests that
00303 *                             did not pass the threshold.
00304 *
00305                               DO 20 K = 1, 2
00306                                  IF( RESULT( K ).GE.THRESH ) THEN
00307                                     IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
00308      $                                 CALL ALAHD( NOUT, PATH )
00309                                     WRITE( NOUT, FMT = 9999 )TRANS, M,
00310      $                                 N, NRHS, NB, ITYPE, K,
00311      $                                 RESULT( K )
00312                                     NFAIL = NFAIL + 1
00313                                  END IF
00314    20                         CONTINUE
00315                               NRUN = NRUN + 2
00316    30                      CONTINUE
00317    40                   CONTINUE
00318                      END IF
00319 *
00320 *                    Generate a matrix of scaling type ISCALE and rank
00321 *                    type IRANK.
00322 *
00323                      CALL DQRT15( ISCALE, IRANK, M, N, NRHS, COPYA, LDA,
00324      $                            COPYB, LDB, COPYS, RANK, NORMA, NORMB,
00325      $                            ISEED, WORK, LWORK )
00326 *
00327 *                    workspace used: MAX(M+MIN(M,N),NRHS*MIN(M,N),2*N+M)
00328 *
00329 *                    Initialize vector IWORK.
00330 *
00331                      DO 50 J = 1, N
00332                         IWORK( J ) = 0
00333    50                CONTINUE
00334                      LDWORK = MAX( 1, M )
00335 *
00336 *                    Test DGELSX
00337 *
00338 *                    DGELSX:  Compute the minimum-norm solution X
00339 *                    to min( norm( A * X - B ) ) using a complete
00340 *                    orthogonal factorization.
00341 *
00342                      CALL DLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
00343                      CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, B, LDB )
00344 *
00345                      SRNAMT = 'DGELSX'
00346                      CALL DGELSX( M, N, NRHS, A, LDA, B, LDB, IWORK,
00347      $                            RCOND, CRANK, WORK, INFO )
00348                      IF( INFO.NE.0 )
00349      $                  CALL ALAERH( PATH, 'DGELSX', INFO, 0, ' ', M, N,
00350      $                               NRHS, -1, NB, ITYPE, NFAIL, NERRS,
00351      $                               NOUT )
00352 *
00353 *                    workspace used: MAX( MNMIN+3*N, 2*MNMIN+NRHS )
00354 *
00355 *                    Test 3:  Compute relative error in svd
00356 *                             workspace: M*N + 4*MIN(M,N) + MAX(M,N)
00357 *
00358                      RESULT( 3 ) = DQRT12( CRANK, CRANK, A, LDA, COPYS,
00359      $                             WORK, LWORK )
00360 *
00361 *                    Test 4:  Compute error in solution
00362 *                             workspace:  M*NRHS + M
00363 *
00364                      CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
00365      $                            LDWORK )
00366                      CALL DQRT16( 'No transpose', M, N, NRHS, COPYA,
00367      $                            LDA, B, LDB, WORK, LDWORK,
00368      $                            WORK( M*NRHS+1 ), RESULT( 4 ) )
00369 *
00370 *                    Test 5:  Check norm of r'*A
00371 *                             workspace: NRHS*(M+N)
00372 *
00373                      RESULT( 5 ) = ZERO
00374                      IF( M.GT.CRANK )
00375      $                  RESULT( 5 ) = DQRT17( 'No transpose', 1, M, N,
00376      $                                NRHS, COPYA, LDA, B, LDB, COPYB,
00377      $                                LDB, C, WORK, LWORK )
00378 *
00379 *                    Test 6:  Check if x is in the rowspace of A
00380 *                             workspace: (M+NRHS)*(N+2)
00381 *
00382                      RESULT( 6 ) = ZERO
00383 *
00384                      IF( N.GT.CRANK )
00385      $                  RESULT( 6 ) = DQRT14( 'No transpose', M, N,
00386      $                                NRHS, COPYA, LDA, B, LDB, WORK,
00387      $                                LWORK )
00388 *
00389 *                    Print information about the tests that did not
00390 *                    pass the threshold.
00391 *
00392                      DO 60 K = 3, 6
00393                         IF( RESULT( K ).GE.THRESH ) THEN
00394                            IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
00395      $                        CALL ALAHD( NOUT, PATH )
00396                            WRITE( NOUT, FMT = 9998 )M, N, NRHS, NB,
00397      $                        ITYPE, K, RESULT( K )
00398                            NFAIL = NFAIL + 1
00399                         END IF
00400    60                CONTINUE
00401                      NRUN = NRUN + 4
00402 *
00403 *                    Loop for testing different block sizes.
00404 *
00405                      DO 100 INB = 1, NNB
00406                         NB = NBVAL( INB )
00407                         CALL XLAENV( 1, NB )
00408                         CALL XLAENV( 3, NXVAL( INB ) )
00409 *
00410 *                       Test DGELSY
00411 *
00412 *                       DGELSY:  Compute the minimum-norm solution X
00413 *                       to min( norm( A * X - B ) )
00414 *                       using the rank-revealing orthogonal
00415 *                       factorization.
00416 *
00417 *                       Initialize vector IWORK.
00418 *
00419                         DO 70 J = 1, N
00420                            IWORK( J ) = 0
00421    70                   CONTINUE
00422 *
00423 *                       Set LWLSY to the adequate value.
00424 *
00425                         LWLSY = MAX( 1, MNMIN+2*N+NB*( N+1 ),
00426      $                          2*MNMIN+NB*NRHS )
00427 *
00428                         CALL DLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
00429                         CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, B,
00430      $                               LDB )
00431 *
00432                         SRNAMT = 'DGELSY'
00433                         CALL DGELSY( M, N, NRHS, A, LDA, B, LDB, IWORK,
00434      $                               RCOND, CRANK, WORK, LWLSY, INFO )
00435                         IF( INFO.NE.0 )
00436      $                     CALL ALAERH( PATH, 'DGELSY', INFO, 0, ' ', M,
00437      $                                  N, NRHS, -1, NB, ITYPE, NFAIL,
00438      $                                  NERRS, NOUT )
00439 *
00440 *                       Test 7:  Compute relative error in svd
00441 *                                workspace: M*N + 4*MIN(M,N) + MAX(M,N)
00442 *
00443                         RESULT( 7 ) = DQRT12( CRANK, CRANK, A, LDA,
00444      $                                COPYS, WORK, LWORK )
00445 *
00446 *                       Test 8:  Compute error in solution
00447 *                                workspace:  M*NRHS + M
00448 *
00449                         CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
00450      $                               LDWORK )
00451                         CALL DQRT16( 'No transpose', M, N, NRHS, COPYA,
00452      $                               LDA, B, LDB, WORK, LDWORK,
00453      $                               WORK( M*NRHS+1 ), RESULT( 8 ) )
00454 *
00455 *                       Test 9:  Check norm of r'*A
00456 *                                workspace: NRHS*(M+N)
00457 *
00458                         RESULT( 9 ) = ZERO
00459                         IF( M.GT.CRANK )
00460      $                     RESULT( 9 ) = DQRT17( 'No transpose', 1, M,
00461      $                                   N, NRHS, COPYA, LDA, B, LDB,
00462      $                                   COPYB, LDB, C, WORK, LWORK )
00463 *
00464 *                       Test 10:  Check if x is in the rowspace of A
00465 *                                workspace: (M+NRHS)*(N+2)
00466 *
00467                         RESULT( 10 ) = ZERO
00468 *
00469                         IF( N.GT.CRANK )
00470      $                     RESULT( 10 ) = DQRT14( 'No transpose', M, N,
00471      $                                    NRHS, COPYA, LDA, B, LDB,
00472      $                                    WORK, LWORK )
00473 *
00474 *                       Test DGELSS
00475 *
00476 *                       DGELSS:  Compute the minimum-norm solution X
00477 *                       to min( norm( A * X - B ) )
00478 *                       using the SVD.
00479 *
00480                         CALL DLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
00481                         CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, B,
00482      $                               LDB )
00483                         SRNAMT = 'DGELSS'
00484                         CALL DGELSS( M, N, NRHS, A, LDA, B, LDB, S,
00485      $                               RCOND, CRANK, WORK, LWORK, INFO )
00486                         IF( INFO.NE.0 )
00487      $                     CALL ALAERH( PATH, 'DGELSS', INFO, 0, ' ', M,
00488      $                                  N, NRHS, -1, NB, ITYPE, NFAIL,
00489      $                                  NERRS, NOUT )
00490 *
00491 *                       workspace used: 3*min(m,n) +
00492 *                                       max(2*min(m,n),nrhs,max(m,n))
00493 *
00494 *                       Test 11:  Compute relative error in svd
00495 *
00496                         IF( RANK.GT.0 ) THEN
00497                            CALL DAXPY( MNMIN, -ONE, COPYS, 1, S, 1 )
00498                            RESULT( 11 ) = DASUM( MNMIN, S, 1 ) /
00499      $                                    DASUM( MNMIN, COPYS, 1 ) /
00500      $                                    ( EPS*DBLE( MNMIN ) )
00501                         ELSE
00502                            RESULT( 11 ) = ZERO
00503                         END IF
00504 *
00505 *                       Test 12:  Compute error in solution
00506 *
00507                         CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
00508      $                               LDWORK )
00509                         CALL DQRT16( 'No transpose', M, N, NRHS, COPYA,
00510      $                               LDA, B, LDB, WORK, LDWORK,
00511      $                               WORK( M*NRHS+1 ), RESULT( 12 ) )
00512 *
00513 *                       Test 13:  Check norm of r'*A
00514 *
00515                         RESULT( 13 ) = ZERO
00516                         IF( M.GT.CRANK )
00517      $                     RESULT( 13 ) = DQRT17( 'No transpose', 1, M,
00518      $                                    N, NRHS, COPYA, LDA, B, LDB,
00519      $                                    COPYB, LDB, C, WORK, LWORK )
00520 *
00521 *                       Test 14:  Check if x is in the rowspace of A
00522 *
00523                         RESULT( 14 ) = ZERO
00524                         IF( N.GT.CRANK )
00525      $                     RESULT( 14 ) = DQRT14( 'No transpose', M, N,
00526      $                                    NRHS, COPYA, LDA, B, LDB,
00527      $                                    WORK, LWORK )
00528 *
00529 *                       Test DGELSD
00530 *
00531 *                       DGELSD:  Compute the minimum-norm solution X
00532 *                       to min( norm( A * X - B ) ) using a
00533 *                       divide and conquer SVD.
00534 *
00535 *                       Initialize vector IWORK.
00536 *
00537                         DO 80 J = 1, N
00538                            IWORK( J ) = 0
00539    80                   CONTINUE
00540 *
00541                         CALL DLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
00542                         CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, B,
00543      $                               LDB )
00544 *
00545                         SRNAMT = 'DGELSD'
00546                         CALL DGELSD( M, N, NRHS, A, LDA, B, LDB, S,
00547      $                               RCOND, CRANK, WORK, LWORK, IWORK,
00548      $                               INFO )
00549                         IF( INFO.NE.0 )
00550      $                     CALL ALAERH( PATH, 'DGELSD', INFO, 0, ' ', M,
00551      $                                  N, NRHS, -1, NB, ITYPE, NFAIL,
00552      $                                  NERRS, NOUT )
00553 *
00554 *                       Test 15:  Compute relative error in svd
00555 *
00556                         IF( RANK.GT.0 ) THEN
00557                            CALL DAXPY( MNMIN, -ONE, COPYS, 1, S, 1 )
00558                            RESULT( 15 ) = DASUM( MNMIN, S, 1 ) /
00559      $                                    DASUM( MNMIN, COPYS, 1 ) /
00560      $                                    ( EPS*DBLE( MNMIN ) )
00561                         ELSE
00562                            RESULT( 15 ) = ZERO
00563                         END IF
00564 *
00565 *                       Test 16:  Compute error in solution
00566 *
00567                         CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
00568      $                               LDWORK )
00569                         CALL DQRT16( 'No transpose', M, N, NRHS, COPYA,
00570      $                               LDA, B, LDB, WORK, LDWORK,
00571      $                               WORK( M*NRHS+1 ), RESULT( 16 ) )
00572 *
00573 *                       Test 17:  Check norm of r'*A
00574 *
00575                         RESULT( 17 ) = ZERO
00576                         IF( M.GT.CRANK )
00577      $                     RESULT( 17 ) = DQRT17( 'No transpose', 1, M,
00578      $                                    N, NRHS, COPYA, LDA, B, LDB,
00579      $                                    COPYB, LDB, C, WORK, LWORK )
00580 *
00581 *                       Test 18:  Check if x is in the rowspace of A
00582 *
00583                         RESULT( 18 ) = ZERO
00584                         IF( N.GT.CRANK )
00585      $                     RESULT( 18 ) = DQRT14( 'No transpose', M, N,
00586      $                                    NRHS, COPYA, LDA, B, LDB,
00587      $                                    WORK, LWORK )
00588 *
00589 *                       Print information about the tests that did not
00590 *                       pass the threshold.
00591 *
00592                         DO 90 K = 7, NTESTS
00593                            IF( RESULT( K ).GE.THRESH ) THEN
00594                               IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
00595      $                           CALL ALAHD( NOUT, PATH )
00596                               WRITE( NOUT, FMT = 9998 )M, N, NRHS, NB,
00597      $                           ITYPE, K, RESULT( K )
00598                               NFAIL = NFAIL + 1
00599                            END IF
00600    90                   CONTINUE
00601                         NRUN = NRUN + 12 
00602 *
00603   100                CONTINUE
00604   110             CONTINUE
00605   120          CONTINUE
00606   130       CONTINUE
00607   140    CONTINUE
00608   150 CONTINUE
00609 *
00610 *     Print a summary of the results.
00611 *
00612       CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS )
00613 *
00614  9999 FORMAT( ' TRANS=''', A1, ''', M=', I5, ', N=', I5, ', NRHS=', I4,
00615      $      ', NB=', I4, ', type', I2, ', test(', I2, ')=', G12.5 )
00616  9998 FORMAT( ' M=', I5, ', N=', I5, ', NRHS=', I4, ', NB=', I4,
00617      $      ', type', I2, ', test(', I2, ')=', G12.5 )
00618       RETURN
00619 *
00620 *     End of DDRVLS
00621 *
00622       END