MODULE cg_alg
  USE cg_precision, ONLY : prec       
CONTAINS
  SUBROUTINE cg(n, b, x, itmax, eps, matr_mult)
    ! Tämä aliohjelma ratkaisee lineaarisen yhtälöryhmän Ax=b
    ! liittogradienttialgoritmilla (CG)
    IMPLICIT NONE
    INTEGER, INTENT(in) :: n               ! Yhtälöryhmän dimensio
    REAL(KIND=prec), INTENT(in) :: b(n)    ! Oikean puolen vektori
    REAL(KIND=prec), INTENT(inout) :: x(n) ! Alkuarvaus ja ratkaisu
    INTEGER, INTENT(in) :: itmax           ! Iteraatioiden lukumäärä
    REAL(KIND=prec), INTENT(in) :: eps     ! Suppenemiskriteeri
    EXTERNAL matr_mult      ! Rutiini, joka laskee matriisi-vektoritulon
    REAL(KIND=prec) :: p(n)        ! Etsintäsuunta
    REAL(KIND=prec) :: ap(n)       ! A*p
    REAL(KIND=prec) :: res(n)      ! Residuaali
    REAL(KIND=prec) :: alpha, beta ! Askelpituudet
    REAL(KIND=prec) :: rnorm       ! Residuaalin normin neliö
    REAL(KIND=prec) :: rnnorm      ! Uuden residuaalin normin neliö
    REAL(KIND=prec) :: rnorm0      ! Alkuresiduaalin normin neliö
    INTEGER :: k                   ! Iteraatiokierros

    ! Laske alkuresiduaali b-Ax
    CALL matr_mult(n, x,res)
    res = b - res
    ! Etsintäsuunta on alkuresiduaali
    p = res

    k = 1
    rnorm = DOT_PRODUCT(res,res)
    rnorm0 = rnorm ! Alkuresiduaalin normin neliö
    DO WHILE( k < itmax .AND. SQRT(rnorm)/SQRT(rnorm0) > eps)
       CALL matr_mult(n, p, ap)
       alpha = rnorm / DOT_PRODUCT(p,ap)
       x = x + alpha * p      ! Päivitä iteraattia
       res = res - alpha * ap ! Päivitä residuaalia
       rnnorm = DOT_PRODUCT(res,res)
       beta = rnnorm / rnorm
       p = res + beta * p     ! Päivitä etsintäsuuntaa
       rnorm = rnnorm
       WRITE(*,*) k, SQRT(rnorm)
       k = k + 1
    END DO
    WRITE (*, '(/,A,I3,A,/)') 'CG suppeni ', k-1, ' iteraatiossa'
  END SUBROUTINE cg
END MODULE cg_alg