acheck Subroutine

private subroutine acheck(me, m, n, nout, eps, v, w, x, y, inform)

Checks the two modes of aprod for LSQR.

acheck may be called to test the user-written subroutine aprod required by LSQR and CRAIG. For some m x n matrix A, aprod with mode = 1 and 2 supplies LSQR and CRAIG with products of the form y := y + Ax and x := x + A'y respectively, where A' means A(transpose). acheck tries to verify that A and A' refer to the same matrix.

Method

We cook up some "unlikely" vectors x and y of unit length and test if y'(y + Ax) = x'(x + A'y).

History

  • 04 Sep 1991 Initial design and code. Michael Saunders, Dept of Operations Research, Stanford University
  • 10 Feb 1992 aprod and eps added as parameters. tol defined via power.

Type Bound

lsqr_solver

Arguments

Type IntentOptional Attributes Name
class(lsqr_solver), intent(inout) :: me
integer, intent(in) :: m

No. of rows of A.
integer, intent(in) :: n

No. of columns of A.
integer, intent(in) :: nout

A file number for printed output.
real(kind=wp), intent(in) :: eps

The machine precision.
real(kind=wp) :: v(n)
real(kind=wp) :: w(m)
real(kind=wp) :: x(n)
real(kind=wp) :: y(m)
integer, intent(out) :: inform

Error indicator. inform = 0 if aprod seems to be consistent. inform = 1 otherwise.

Calls

proc~~acheck~~CallsGraph proc~acheck lsqr_module::lsqr_solver%acheck aprod aprod proc~acheck->aprod proc~dcopy lsqpblas_module::dcopy proc~acheck->proc~dcopy proc~ddot lsqpblas_module::ddot proc~acheck->proc~ddot proc~dnrm2 lsqpblas_module::dnrm2 proc~acheck->proc~dnrm2 proc~dscal lsqpblas_module::dscal proc~acheck->proc~dscal

Source Code

   subroutine acheck( me, m, n, nout, eps, &
                     v, w, x, y, inform )

   implicit none

   class(lsqr_solver),intent(inout) :: me
   integer,intent(in)   :: m       !! No. of rows of A.
   integer,intent(in)   :: n       !! No. of columns of A.
   integer,intent(in)   :: nout    !! A file number for printed output.
   integer,intent(out)  :: inform  !! Error indicator.
                                   !! inform = 0 if aprod seems to be
                                   !! consistent.
                                   !! inform = 1 otherwise.
   real(wp),intent(in)  :: eps     !! The machine precision.
   real(wp)             :: v(n)
   real(wp)             :: w(m)
   real(wp)             :: x(n)
   real(wp)             :: y(m)

   real(wp),parameter :: power = 0.5_wp   !! eps**power is used as the tolerance for judging
                                          !! whether `y'(y + Ax)  =  x'(x + A'y)`
                                          !! to sufficient accuracy.
                                          !! power should be in the range (0.25, 0.9) say.
                                          !! For example, power = 0.75 means that we are happy
                                          !! if three quarters of the available digits agree.
                                          !! power = 0.5 seems a reasonable requirement
                                          !! (asking for half the digits to agree).

   integer :: i, j
   real(wp) :: alfa, beta, t, test1, test2, test3, tol

   tol = eps**power
   if (nout /= 0) write(nout, '(//A)') &
      'Enter acheck. Test of aprod for LSQR and CRAIG'

   ! ==================================================================
   ! Cook up some "unlikely" vectors x and y of unit length.
   ! ==================================================================
   t = one
   do j = 1, n
      t    = t + one
      x(j) = sqrt( t )
   end do

   t = one
   do i = 1, m
      t    = t + one
      y(i) = one / sqrt( t )
   end do

   alfa = dnrm2 ( n, x, 1 )
   beta = dnrm2 ( m, y, 1 )
   call dscal ( n, (one / alfa), x, 1 )
   call dscal ( m, (one / beta), y, 1 )

   ! ==================================================================
   ! Test if y'(y + Ax) = x'(x + A'y).
   ! ==================================================================

   ! First set w = y + Ax, v = x + A'y.

   call dcopy ( m, y, 1, w, 1 )
   call dcopy ( n, x, 1, v, 1 )
   call me%aprod ( 1, m, n, x, w )
   call me%aprod ( 2, m, n, v, y )

   ! Now set alfa = y'w, beta = x'v.

   alfa  = ddot  ( m, y, 1, w, 1 )
   beta  = ddot  ( n, x, 1, v, 1 )
   test1 = abs( alfa - beta )
   test2 = one  +  abs( alfa )  +  abs( beta )
   test3 = test1 / test2

   ! See if alfa and beta are essentially the same.

   if ( test3 <= tol ) then
      inform = 0
      if (nout /= 0) write(nout, '(1P,A,1X,E10.1)') &
         'aprod seems OK. Relative error =', test3
   else
      inform = 1
      if (nout /= 0) write(nout, '(1P,A,1X,E10.1)') &
         'aprod seems incorrect. Relative error =', test3
   end if

   end subroutine acheck