test Subroutine

private subroutine test(me, m, n, nduplc, npower, damp)

This is an example driver routine for running LSQR. It generates a test problem, solves it, and examines the results. Note that subroutine aprod must be declared external if it is used only in the call to LSQR (and acheck).

History

  • 1982---1991: Various versions implemented.
  • 04 Sep 1991: "wantse" added to argument list of LSQR, making standard errors optional.
  • 10 Feb 1992: Revised for use with lsqrchk fortran.
  • 31 Mar 2005: changed atol = eps*0.666667 to eps0.99 to increase accuracy of the solution. LSQR appears to be successful on all 18 test problems except 5 and 6 (which are over-determined and too ill-conditioned to permit any correct digits). The output from an Intel Xeon system with g77 is in LSQR.LIS. The two "appears to have failed" messages are no cause for alarm. Michael Saunders, Dept of Operations Research, Stanford University.

Type Bound

test_solver

Arguments

Type IntentOptional Attributes Name
class(test_solver), intent(inout) :: me
integer :: m
integer :: n
integer :: nduplc
integer :: npower
real(kind=wp) :: damp

Calls

proc~~test~~CallsGraph proc~test lsqrtest_module::test_solver%test 4 4 proc~test->4 proc~acheck lsqr_module::lsqr_solver%acheck proc~test->proc~acheck proc~dcopy lsqpblas_module::dcopy proc~test->proc~dcopy proc~dnrm2 lsqpblas_module::dnrm2 proc~test->proc~dnrm2 proc~lsqr lsqr_module::lsqr_solver%LSQR proc~test->proc~lsqr proc~lstp lsqrtest_module::test_solver%lstp proc~test->proc~lstp proc~xcheck lsqr_module::lsqr_solver%xcheck proc~test->proc~xcheck proc~acheck->proc~dcopy proc~acheck->proc~dnrm2 aprod aprod proc~acheck->aprod proc~ddot lsqpblas_module::ddot proc~acheck->proc~ddot proc~dscal lsqpblas_module::dscal proc~acheck->proc~dscal proc~lsqr->proc~dcopy proc~lsqr->proc~dnrm2 proc~lsqr->aprod proc~d2norm lsqr_module::d2norm proc~lsqr->proc~d2norm proc~lsqr->proc~dscal proc~lstp->proc~dcopy proc~lstp->proc~dnrm2 proc~aprod1 lsqrtest_module::test_solver%aprod1 proc~lstp->proc~aprod1 proc~lstp->proc~dscal proc~hprod lsqrtest_module::hprod proc~lstp->proc~hprod proc~xcheck->proc~dcopy proc~xcheck->proc~dnrm2 proc~xcheck->aprod proc~xcheck->proc~dscal proc~aprod1->proc~hprod

Called by

proc~~test~~CalledByGraph proc~test lsqrtest_module::test_solver%test proc~lsqr_test lsqrtest_module::lsqr_test proc~lsqr_test->proc~test program~main main program~main->proc~lsqr_test

Source Code

    subroutine test  ( me, m, n, nduplc, npower, damp )

    class(test_solver),intent(inout) :: me
    integer :: m, n, nduplc, npower
    real(wp) :: damp

    integer,parameter :: maxm = 2000
    integer,parameter :: maxn = 2000
    integer,parameter :: mxmn = 2000
    real(wp),parameter :: eps = epsilon(1.0_wp) !! machine precision
    character(len=34),parameter :: line = '----------------------------------'

    integer :: inform, istop, itnlim, j, itn, maxmn, minmn, nprint
    integer :: locd, lochy, lochz, locw, ltotal
    logical :: wantse
    real(wp) :: b(maxm),  u(maxm), &
                v(maxn),  w(mxmn), x(maxn), &
                se(maxn), xtrue(maxn), y(mxmn)
    real(wp) :: atol, btol, conlim, &
                anorm, acond, rnorm, arnorm, &
                enorm, etol, xnorm, test1, test2, test3, wnorm

    if (m > maxm  .or.  n > maxn) then
        ! m or n too large.
        write(me%nout, 8000)
        return
    end if

    ! Set the desired solution xtrue.
    ! For least-squares problems, this is it.
    ! For underdetermined systems, lstp may alter it.

    do j = 1, n
        ! xtrue(j) = one
        xtrue(j) = j * 0.1_wp
    end do

    ! Generate the specified test problem.
    ! The workspace array  rw  is used for the following vectors:
    !    d(minmn), hy(m), hz(n), w(maxmn).
    ! The vectors  d, hy, hz  will define the test matrix A.
    ! w is needed for workspace in aprod1 and aprod2.

    maxmn  = max( m, n )
    minmn  = min( m, n )
    locd   = 1
    lochy  = locd  + minmn
    lochz  = lochy + m
    locw   = lochz + n
    ltotal = locw  + maxmn - 1
    if (ltotal > lenrw) then
        ! Not enough workspace.
        write(me%nout, 9000) ltotal
    end if

    call me%lstp  ( m, n, maxmn, minmn, nduplc, npower, damp, xtrue, &
                b, me%rw(locd), me%rw(lochy), me%rw(lochz), me%rw(locw), &
                acond, rnorm )

    write(me%nout, 1000) line, line, &
                         m, n, nduplc, npower, damp, acond, rnorm, &
                         line, line

    ! Check that aprod generates y + Ax and x + A'y consistently.
    call me%acheck( m, n, me%nout, eps, v, w, x, y, inform )

    if (inform > 0) then
        write(me%nout, '(a)') ' Check eps and power in subroutine acheck'
        stop
    end if

    ! Solve the problem defined by aprod, damp and b.
    ! Copy the rhs vector b into u  (LSQR will overwrite u)
    ! and set the other input parameters for LSQR.
    ! We ask for standard errors only if they are well-defined.

    call dcopy ( m, b, 1, u, 1 )
    !---  wantse = m > n  .or.  damp > zero
    wantse = .false.
    atol   = eps**0.99_wp
    btol   = atol
    conlim = 1000.0_wp * acond
    itnlim = 4*(m + n + 50)

    call me%lsqr(  m, n, damp, wantse, &
                   u, v, w, x, se, &
                   atol, btol, conlim, itnlim, me%nout, &
                   istop, itn, anorm, acond, rnorm, arnorm, xnorm )

    ! Examine the results.

    !-  if (damp == zero) then
    !-     write(nout, 2000)       xnorm, rnorm, arnorm
    !-  else
    !-     write(nout, 2100) damp, xnorm, rnorm, arnorm
    !-  end if

    call me%xcheck( m, n, me%nout, anorm, damp, eps, &
                    b, u, v, w, x, &
                    inform, test1, test2, test3 )

    ! Print the solution and standard error estimates from  LSQR.

    nprint = min( m, n, 8 )
    write(me%nout, 2500)    (j, x(j) , j = 1, nprint)
    if ( wantse ) then
        write(me%nout, 2600) (j, se(j), j = 1, nprint)
    end if

    ! Print a clue about whether the solution looks OK.

    do j = 1, n
        w(j)  = x(j) - xtrue(j)
    end do
    wnorm    = dnrm2 ( n, w    , 1 )
    xnorm    = dnrm2 ( n, xtrue, 1 )
    enorm    = wnorm / (one + xnorm)
    etol     = 0.001_wp
    if (enorm <= etol) then
        write(me%nout, 3000) enorm
    else
        write(me%nout, 3100) enorm
    end if
    return

    1000 format(1p &
    // 1x, 2a &
    /  ' Least-Squares Test Problem      P(', 4i5, e12.2, ' )' &
    // ' Condition no. =', e12.4,  '     Residual function =', e17.9 &
    /  1x, 2a)
    2000 format(1p &
    // ' We are solving    min norm(Ax - b)    with no damping.' &
    // ' Estimates from LSQR:' &
    /  '    norm(x)         =', e10.3, ' = xnorm' &
    /  '    norm(r)         =', e10.3, ' = rnorm' &
    /  '    norm(A''r)       =', e10.3, ' = arnorm')
    2100 format(1p &
    // ' We are solving    min norm(Ax - b)    with damp =', e10.3 &
    /  '                           (damp*x)' &
    // ' Estimates from LSQR:' &
    /  '    norm(x)         =', e10.3, ' = xnorm' &
    /  '    norm(rbar)      =', e10.3, ' = rnorm' &
    /  '    norm(Abar''rbar) =', e10.3, ' = arnorm')
    2500 format(//' Solution  x:' / 4(i6, g14.6))
    2600 format(/ ' Standard errors  se:' / 4(i6, g14.6))
    3000 format(1p / ' LSQR  appears to be successful.', &
            '     Relative error in  x  =', e10.2)
    3100 format(1p / ' LSQR  appears to have failed.  ', &
            '     Relative error in  x  =', e10.2)
    8000 format(/ ' XXX  m or n is too large.')
    9000 format(/ ' XXX  Insufficient workspace.', &
            '  The length of  rw  should be at least', i6)

    end subroutine test