lstp Subroutine

private subroutine lstp(me, m, n, maxmn, minmn, nduplc, npower, damp, x, b, d, hy, hz, w, acond, rnorm)

lstp generate a sparse least-squares test problem of the form ( A )x = ( b ) ( dampI ) ( 0 ) for solution by LSQR, or a sparse underdetermined system Ax + damps = b for solution by CRAIG. The matrix A is m by n and is constructed in the form A = HYD*HZ, where D is an m by n diagonal matrix, and HY and HZ are Householder transformations.

m and n may contain any positive values. The first 8 parameters are input to lstp. The last 8 are output. If m >= n or damp = 0, the true solution is x as given. Otherwise, x is modified to contain the true solution.

Type Bound

test_solver

Arguments

Type IntentOptional Attributes Name
class(test_solver), intent(inout) :: me
integer :: m
integer :: n
integer :: maxmn
integer :: minmn
integer :: nduplc
integer :: npower
real(kind=wp) :: damp
real(kind=wp) :: x(n)
real(kind=wp) :: b(m)
real(kind=wp) :: d(minmn)
real(kind=wp) :: hy(m)
real(kind=wp) :: hz(n)
real(kind=wp) :: w(maxmn)
real(kind=wp) :: acond
real(kind=wp) :: rnorm

Calls

proc~~lstp~~CallsGraph proc~lstp lsqrtest_module::test_solver%lstp proc~aprod1 lsqrtest_module::test_solver%aprod1 proc~lstp->proc~aprod1 proc~dcopy lsqpblas_module::dcopy proc~lstp->proc~dcopy proc~dnrm2 lsqpblas_module::dnrm2 proc~lstp->proc~dnrm2 proc~dscal lsqpblas_module::dscal proc~lstp->proc~dscal proc~hprod lsqrtest_module::hprod proc~lstp->proc~hprod proc~aprod1->proc~hprod

Called by

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

Source Code

    subroutine lstp  ( me, m, n, maxmn, minmn, nduplc, npower, damp, x, &
                       b, d, hy, hz, w, acond, rnorm )

    class(test_solver),intent(inout) :: me
    integer :: m, n, maxmn, minmn, nduplc, npower
    real(wp) :: damp, acond, rnorm
    real(wp) :: b(m), x(n), d(minmn), hy(m), hz(n), w(maxmn)

    integer :: i, j
    real(wp) alfa, beta, dampsq, t

    real(wp),parameter :: fourpi = 4.0_wp * acos(-1.0_wp)   !4.0*3.141592

    ! ------------------------------------------------------------------
    ! Make two vectors of norm 1.0 for the Householder transformations.
    ! fourpi  need not be exact.
    ! ------------------------------------------------------------------
    minmn  = min( m, n )
    dampsq = damp**2
    alfa   = fourpi / m
    beta   = fourpi / n

    do i = 1, m
        hy(i) = sin( i * alfa )
    end do

    do i = 1, n
        hz(i) = cos( i * beta )
    end do

    alfa   = dnrm2 ( m, hy, 1 )
    beta   = dnrm2 ( n, hz, 1 )
    call dscal ( m, (- one / alfa), hy, 1 )
    call dscal ( n, (- one / beta), hz, 1 )

    ! ------------------------------------------------------------------
    ! Set the diagonal matrix  D.  These are the singular values of  A.
    ! ------------------------------------------------------------------
    do i = 1, minmn
        j     = (i - 1 + nduplc) / nduplc
        t     =  j * nduplc
        t     =  t / minmn
        d(i)  =  t**npower
    end do

    acond  = (d(minmn)**2 + dampsq) / (d(1)**2 + dampsq)
    acond  = sqrt( acond )

    ! ------------------------------------------------------------------
    ! Set the true solution   x.
    ! It must be of the form  x = Z ( w )  for some  w.
    !                               ( 0 )
    ! ------------------------------------------------------------------
    call hprod ( n, hz, x, w )

    do i = m + 1, n
        w(i)  = zero
    end do

    call hprod ( n, hz, w, x )

    ! Solve D r1bar = dampsq x1bar
    ! where r1bar and x1bar are both in w.

    do i = 1, minmn
        w(i)  = dampsq * w(i) / d(i)
    end do

    ! Set r2bar to be anything.  (It is empty if m <= n)
    ! Then form r = Y rbar (again in w).

    do i = minmn+1, m
        w(i)  = one
    end do

    call Hprod ( m, hy, w, w )

    ! Compute the rhs    b = r  +  Ax.

    rnorm  = dnrm2 ( m, w, 1 )
    call dcopy ( m, w, 1, b, 1 )
    call me%aprod1( m, n, maxmn, minmn, x, b, d, hy, hz, w )

    end subroutine lstp