! Richardson Residual Minimization
!   Solve r=b-A*x iteratively
!
!   Charles O'Neill Oct 15, 2009
!    charles.oneill@gmail.com
!    www.caselab.okstate.edu
!
!============================================================
! Latex Documentation
!\section{Richardson Residual Minimization}
!The canonical form is\[Ax=b\] which in residual form is\[r=b-Ax\]
!Iterative solutions of the residual equation are used\[x_{i+1}=x_{i}+\alpha_{i}r_{i}\]
!The Richardson Residual Minimization scheme for $\alpha$ is,\[\alpha_{i}=\frac{\sum r_{i}q_{i}}{\sum q_{i}q_{i}}\]
!where $q_{i}$ is calculated as\[q_{i}=Ar_{i}\]
!This method is a linear equation residual minimization along the steepest descent direction.    
!
! Copyright (c) 2009 Charles O'Neill
!
! Permission is hereby granted, free of charge, to any person
! obtaining a copy of this software and associated documentation
! files (the "Software"), to deal in the Software without
! restriction, including without limitation the rights to use,
! copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the
! Software is furnished to do so, subject to the following
! conditions:
!
! The above copyright notice and this permission notice shall be
! included in all copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
! EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
! OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
! NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
! HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
! WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
! FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
! OTHER DEALINGS IN THE SOFTWARE.
!

module RichardsonResMin
  ! Real Single Precision
  integer, parameter,private :: SP = 4 ! 4 gives 32 bits
  ! Real Double Precision
  integer, parameter,private :: DP = 8  ! 4 gives 32 bits, 8 gives 64 bits (double)
  ! Real Working Precision
  integer, parameter,private :: WP = SP   ! User specified precision
  
contains
    
    subroutine RRM(n,imax,RsdTol,a0)
        integer :: n   ! Unknowns Vector Size
        integer :: i, imax
        real(WP) :: alpha
        real(WP) :: r(n), q(n), a0(n)
        real(WP) :: CurRes, RsdTol
        logical :: IterateFlag        
        
        IterateFlag = .True.
        i = 1
        
        ! Compute Residual
        CurRes = sum(r*r)   

        ! Iterate
        do while( IterateFlag )   
    
            ! Compute q=A*r
            q = ExampleAfunction(n,r) ! Replace with your actual "A*r" function
            
            ! Compute Alpha
            alpha = sum(r*q)/sum(q*q) ! Richardson Residual Minimization
    
            ! Update Solution
            a0 = a0 +  alpha * r

            ! Compute Residual
            CurRes = sum(r*r)

            ! Stopping Criteria
            if(i>imax)then
                stop "RRM iteration failure"
            endif
            if(CurRes < RsdTol)then
                exit
            endif

            ! Update iteration
            i=i+1
        enddo
    end subroutine

    ! Bogus Example A*r function.
    function ExampleAFunction(n,r)
        integer :: n
        real(WP) :: ExampleAFunction(n)
        real(WP) :: r(n)
        ExampleAFunction = 1.0_wp * r
    end function

end module