xcheck Subroutine

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

   implicit none

   class(lsqr_solver),intent(inout) :: me
   integer,intent(in)   :: m  !! The number of rows in A.
   integer,intent(in)   :: n  !! The number of columns in A.
   integer,intent(in)   :: nout  !! A file number for printed output.
                                 !! If nout = 0, nothing is printed.
   integer,intent(out)  :: inform   !! inform = 0 if b = 0 and x = 0.
                                    !! inform = 1, 2 or 3 if x seems to
                                    !! solve systems 1 2 or 3 above.
   real(wp),intent(in)  :: anorm    !! An estimate of norm(A) or
                                    !! norm( A, delta*I ) if delta > 0.
                                    !! Normally this will be available
                                    !! from LSQR or CRAIG.
   real(wp),intent(in)  :: damp     !! Possibly defines a damped problem.
   real(wp),intent(in)  :: eps      !! Machine precision.
   real(wp),intent(out) :: test1    !! These are dimensionless quantities
                                    !! that should be "small" if x does
                                    !! seem to solve one of the systems.
                                    !! "small" means less than
                                    !! tol = eps**power, where power is
                                    !! defined as a parameter below.
   real(wp),intent(out) :: test2    !! These are dimensionless quantities
                                    !! that should be "small" if x does
                                    !! seem to solve one of the systems.
                                    !! "small" means less than
                                    !! tol = eps**power, where power is
                                    !! defined as a parameter below.
   real(wp),intent(out) :: test3    !! These are dimensionless quantities
                                    !! that should be "small" if x does
                                    !! seem to solve one of the systems.
                                    !! "small" means less than
                                    !! tol = eps**power, where power is
                                    !! defined as a parameter below.
   real(wp),intent(in)  :: b(m)     !! The right-hand side of Ax = b etc.
   real(wp),intent(out) :: u(m)     !! On exit, u = r (where r = b - Ax).
   real(wp),intent(out) :: v(n)     !! On exit, v = A'r.
   real(wp),intent(out) :: w(n)     !! On exit, w = A'r - damp^2 x.
   real(wp),intent(in)  :: x(n)     !! The given estimate of a solution.

   real(wp),parameter :: power = 0.5_wp

   integer :: j
   real(wp) :: bnorm, dampsq, rho1, rho2, sigma1, sigma2, &
               tol, snorm, xnorm, xsnorm
   real(wp),dimension(n) :: xtmp

   dampsq = damp**2
   tol    = eps**power
   xtmp   = x

   ! Compute u = b - Ax via u = -b + Ax, u = -u.
   ! This is usual residual vector r.

   call dcopy ( m, b, 1, u, 1 )
   call dscal ( m, (-one), u, 1 )
   call me%aprod ( 1, m, n, xtmp, u )
   call dscal ( m, (-one), u, 1 )

   ! Compute v = A'u via v = 0, v = v + A'u.

   do j = 1, n
      v(j)  = zero
   end do
   call me%aprod ( 2, m, n, v, u )

   ! Compute w = A'u - damp**2 * x.
   ! This will be close to zero in all cases
   ! if x is close to a solution.

   call dcopy ( n, v, 1, w, 1 )
   if (damp /= zero) then
      do j = 1, n
         w(j)  = w(j)  -  dampsq * x(j)
      end do
   end if

   ! Compute the norms associated with b, x, u, v, w.

   bnorm  = dnrm2 ( m, b, 1 )
   xnorm  = dnrm2 ( n, x, 1 )
   rho1   = dnrm2 ( m, u, 1 )
   sigma1 = dnrm2 ( n, v, 1 )
   if (nout /= 0) then
      write(nout, '(//A)') 'Enter xcheck. Does x solve Ax = b, etc?'
      write(nout, '(1P,A,E10.3)')      ' damp            =', damp
      write(nout, '(1P,A,E10.3)')      ' norm(x)         =', xnorm
      write(nout, '(1P,A,E15.8,A)')    ' norm(r)         =', rho1,   ' = rho1'
      write(nout, '(1P,A,E10.3,5X,A)') ' norm(A''r)       =', sigma1, ' = sigma1'
   end if

   if (damp == zero) then
      rho2   = rho1
      sigma2 = sigma1
   else
      rho2   = sqrt( rho1**2  +  dampsq * xnorm**2 )
      sigma2 = dnrm2 ( n, w, 1 )
      snorm  = rho1 / damp
      xsnorm = rho2 / damp
      if (nout /= 0) then
         write(nout, '(1P/A,E10.3)')      ' norm(s)         =', snorm
         write(nout, '(1P,A,E10.3)')      ' norm(x,s)       =', xsnorm
         write(nout, '(1P,A,E15.8,A)')    ' norm(rbar)      =', rho2,    ' = rho2'
         write(nout, '(1P,A,E10.3,5X,A)') ' norm(Abar''rbar) =', sigma2, ' = sigma2'
      end if
   end if

   ! See if x seems to solve Ax = b or min norm(Ax - b)
   ! or the damped least-squares system.

   if (bnorm == zero  .and.  xnorm == zero) then
      inform = 0
      test1  = zero
      test2  = zero
      test3  = zero
   else
      inform = 4
      test1  = rho1 / (bnorm  +  anorm * xnorm)
      test2  = zero
      if (rho1 > zero) test2  = sigma1 / (anorm * rho1)
      test3  = test2
      if (rho2 > zero) test3  = sigma2 / (anorm * rho2)
      if (test3 <= tol) inform = 3
      if (test2 <= tol) inform = 2
      if (test1 <= tol) inform = 1
   end if

   if (nout /= 0) then
      write(nout, '(/A,I2)')        ' inform          =', inform
      write(nout, '(1P,A,E10.3)')   ' tol             =', tol
      write(nout, '(1P,A,E10.3,A)') ' test1           =', test1, ' (Ax = b)'
      write(nout, '(1P,A,E10.3,A)') ' test2           =', test2, ' (least-squares)'
      write(nout, '(1P,A,E10.3,A)') ' test3           =', test3, ' (damped least-squares)'
   end if

   end subroutine xcheck