CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   Main CFD Forum (http://www.cfd-online.com/Forums/main/)
-   -   CG, BICGSTAB(2) : problem with matrix operation and boundary conditions (http://www.cfd-online.com/Forums/main/72622-cg-bicgstab-2-problem-matrix-operation-boundary-conditions.html)

moomba February 12, 2010 05:34

CG, BICGSTAB(2) : problem with matrix operation and boundary conditions
 
Hi

I am trying to develop a Parallel Conjugate Gradient (CG) and a Parallel BICGSTAB(2) to solve Dynamic Pressure in a Low Mach DNS solver. I have been using the Van der Vorst paper. :)

The system is working well with perdiodic boundary conditions, but fails to work with Neumann (dP/dxn = 0, converge but bad solution) and Dirichlet (P=constant, do not converge).

Because the problem is exactly the same for both system, I think it comes from the matrix multiplication.

Could someone explain me where I made a mistake ? Here is my fortran 90 code for the BICGSTAB(2).
The subroutine "Matrix_A_Multiply" is made here for periodicity on x1 and x2. The code is of course not optimised yet.

Code:

      module BICGSTAB_2

      use typegrid

      contains


      subroutine BICGSTAB_2_CORE(Grid,b,x)

      implicit none

      type (grid_type) :: Grid

      real(8), dimension(1:Grid%Nx(1),1:Grid%Nx(2),1:Grid%Nx(3)), intent(in) :: b
      real(8), dimension(1:Grid%Nx(1),1:Grid%Nx(2),1:Grid%Nx(3)), intent(out) :: x

      real(8), dimension(1:Grid%Nx(1),1:Grid%Nx(2),1:Grid%Nx(3)) :: xx0
      real(8), dimension(1:Grid%Nx(1),1:Grid%Nx(2),1:Grid%Nx(3)) :: u,v,r,s,t,w

      real(8) :: omega_1,omega_2,rho_1,rho_0,mu,nu,beta,gam,tau,alpha

      integer :: i

      xx0 = 0.0d0

      CALL Matrix_A_Multiply(Grid,xx0,r) ! (x_in,x_out) --> x_out = A * x_in

      r(:,:,:)  = b(:,:,:) - r(:,:,:)

      rho_0    = 1.0d0
      u(:,:,:) = 0.0d0
      alpha    = 0.0d0
      omega_2    = 1.0d0

      do i=0,1000,2

        rho_0 = - omega_2 * rho_0


        ! Even BiCG step :

        rho_1 = Scalar_Product(Grid,r,r)

        beta = alpha * rho_1 / rho_0

        rho_0 = rho_1

        u(:,:,:) = r(:,:,:) - beta * u(:,:,:)

        CALL Matrix_A_Multiply(Grid,u,v)

        gam = Scalar_Product(Grid,v,r) ! return v(1)*r(1) + v(2)*r(2) + v(3)*r(3) + ... + v(n)*r(n)

        alpha = rho_0 / gam

        r(:,:,:) = r(:,:,:) - alpha * v(:,:,:)

        CALL Matrix_A_Multiply(Grid,r,s)

        x(:,:,:) = x(:,:,:) + alpha * u(:,:,:)


        ! Odd BiCG step :

        rho_1 = Scalar_Product(Grid,r,s)

        beta = alpha * rho_1 / rho_0

        rho_0 = rho_1

        v(:,:,:) = s(:,:,:) - beta * v(:,:,:)

        CALL Matrix_A_Multiply(Grid,v,w)

        gam = Scalar_Product(Grid,w,r)

        alpha = rho_0 / gam

        u(:,:,:) = r(:,:,:) - beta * u(:,:,:)

        r(:,:,:) = r(:,:,:) - alpha * v(:,:,:)

        s(:,:,:) = s(:,:,:) - alpha * w(:,:,:)

        CALL Matrix_A_Multiply(Grid,s,t)


        ! GCR(2)-part :

        omega_1 = Scalar_Product(Grid,r,s)

        mu = Scalar_Product(Grid,s,s)

        nu = Scalar_Product(Grid,s,t)

        tau = Scalar_Product(Grid,t,t)

        omega_2 = Scalar_Product(Grid,r,t)

        tau = tau - ( nu * nu / mu )

        omega_2 = ( omega_2 - nu * omega_1 / mu ) / tau
! OPTI ?

        omega_1 = ( omega_1 - nu * omega_2 ) / mu

        x(:,:,:) = x(:,:,:) + omega_1 * r(:,:,:) + omega_2 * s(:,:,:) + alpha * u(:,:,:)

        r(:,:,:) = r(:,:,:) - omega_1 * s(:,:,:) - omega_2 * t(:,:,:)

        ! test of convergence
        if( max( ABS(maxval(r)) , ABS(minval(r)) ) < 1e-10 ) then
            print *,"Converged in",i,"iterations"
            return
        end if

        u(:,:,:) = u(:,:,:) - omega_1 * v(:,:,:) - omega_2 * w(:,:,:)

      end do

      print *,"ERROR, BICGSTAB(2) not converged"

      end subroutine BICGSTAB_2_CORE



      subroutine Matrix_A_Multiply(Grid,p,q)

      implicit none

      type (grid_type) :: Grid
      real(8), dimension(1:Grid%Nx(1),1:Grid%Nx(2),1:Grid%Nx(3)), intent(in) :: p
      real(8), dimension(1:Grid%Nx(1),1:Grid%Nx(2),1:Grid%Nx(3)), intent(out) :: q

      real(8), dimension(0:Grid%Nx(1)+1,0:Grid%Nx(2)+1) :: pl

      real(8), dimension(1:3) :: SDx

      integer :: i,j,n

      pl=77.0
      pl(1:Grid%Nx(1),1:Grid%Nx(2)) = p(1:Grid%Nx(1),1:Grid%Nx(2),1)
! Opti to be made

      pl(          0,1:Grid%Nx(2)) = pl(Grid%Nx(1)-1,1:Grid%Nx(2))
      pl(Grid%Nx(1)+1,1:Grid%Nx(2)) = pl(          2,1:Grid%Nx(2))

      pl(1:Grid%Nx(1),          0) = pl(1:Grid%Nx(1),Grid%Nx(2)-1)
      pl(1:Grid%Nx(1),Grid%Nx(2)+1) = pl(1:Grid%Nx(1),          2)


      do n=1,2
        SDx(n)=(-1.0/(Grid%h(n)*Grid%h(n)))
      end do

      do j=1,Grid%Nx(2)
        do i=1,Grid%Nx(1)
            q(i,j,1)=SDx(1)*(-pl(i-1,j)-pl(i+1,j)+2.0*pl(i,j)) +    &
            &        SDx(2)*(-pl(i,j-1)-pl(i,j+1)+2.0*pl(i,j))
        end do
      end do

      return

      end subroutine Matrix_A_Multiply



      function Scalar_Product(Grid,x1,x2)

      implicit none

      type (grid_type) :: Grid
      real(8), dimension(1:Grid%Nx(1),1:Grid%Nx(2),1:Grid%Nx(3)), intent(in) :: x1,x2
      real(8) :: Scalar_Product
      integer :: i,j,k
     
      Scalar_Product=0.0d0
      do k=1,Grid%Nx(3)
        do j=1,Grid%Nx(2)
            do i=1,Grid%Nx(1)
              Scalar_Product = Scalar_Product + x1(i,j,k) * x2(i,j,k)
            end do
        end do
      end do

      return

      end function Scalar_Product

      end module BICGSTAB_2

In fact, for periodicity on x1, I copy point pl(Nx(1)-1) in ghost pl(0), and the point pl(2) in ghost pl(Nx(1)+1).

For Neumann dP/dxn=0, I copy point pl(1) in ghost pl(0), and pl(Nx(1)) in ghost pl(Nx(1)+1).

And for dirichlet P=0 on Nx(1) for example, I put pl(Nx(1))=0 and the output q(Nx(1))=0.


I of course made a mistake, but where :confused:

moomba February 12, 2010 19:47

I made the problem far more simple : a GC to solve f=x² on [0,1] with Dirichlet on x (0 at x=0, and 1 at x=1 :D), and with dp/dy=0 on y border :

Main :
Code:

      program Main

      use typegrid
      use CG

      implicit none

      type (grid_type)        :: Grid

      real(8), allocatable :: F(:,:),P(:,:)
      real(8) :: PI,xx,yy
      integer :: i,j

      PI=3.141592653589793
      Grid%Ndim=2
      Grid%Nx(1)=100
      Grid%Nx(2)=100
      Grid%Nx(3)=1
      Grid%h(1)=1.0/(Grid%Nx(1)-1)
      Grid%h(2)=1.0/(Grid%Nx(2)-1)
      Grid%h(3)=0.0

      ALLOCATE(F(1:Grid%Nx(1),1:Grid%Nx(2)),P(1:Grid%Nx(1),1:Grid%Nx(2)))

      F(:,:)=2.0

      open(10,file="GC0.dat")
        do i=1,Grid%Nx(1)
            do j=1,Grid%Nx(2)
              Write(10,*)(i-1)*Grid%h(1),(j-1)*Grid%h(2),F(i,j)
            end do
        end do
      close(10)


      CALL CG_CORE(Grid,F,P)


      open(10,file="GC1.dat")
        do i=1,Grid%Nx(1)
            do j=1,Grid%Nx(2)
              xx=(i-1)*Grid%h(1)
              yy=(j-1)*Grid%h(2)
              Write(10,*)(i-1)*Grid%h(1),(j-1)*Grid%h(2),P(i,j)
            end do
        end do
      close(10)

      DEALLOCATE(F,P)

      end program Main

GC :
Code:

      module CG

      use typegrid

      integer :: n01,n11,n02,n12

      contains


      subroutine CG_CORE(Grid,b,x)

      implicit none

      type (grid_type) :: Grid

      real(8), dimension(1:Grid%Nx(1),1:Grid%Nx(2)), intent(in) :: b
      real(8), dimension(1:Grid%Nx(1),1:Grid%Nx(2)), intent(out) :: x

      real(8), dimension(1:Grid%Nx(1),1:Grid%Nx(2)) :: r,p,q,xx0

      real(8) :: rho,rhop1,beta,alpha,res

      integer :: i

      n01 = 2
      n11 = Grid%Nx(1)-1
      n02 = 1
      n12 = Grid%Nx(2)

      xx0(n01:n11,n02:n12) = 0.0d0

      CALL Matrix_A_Multiply(Grid,xx0,r) ! (x_in,x_out) --> x_out = A * x_in

      r(n01:n11,n02:n12) = b(n01:n11,n02:n12) - r(n01:n11,n02:n12)

      p(n01:n11,n02:n12) = 0.0d0

      beta = 0.0d0

      rho = Scalar_Product(Grid,r,r) ! return v(1)*r(1) + v(2)*r(2) + v(3)*r(3) + ... + v(n)*r(n)


      do i=1,1000

        p(n01:n11,n02:n12) = r(n01:n11,n02:n12) + beta * p(n01:n11,n02:n12)

        CALL Matrix_A_Multiply(Grid,p,q)

        alpha = rho / Scalar_Product(Grid,p,q)

        x(n01:n11,n02:n12) = x(n01:n11,n02:n12) + alpha * p(n01:n11,n02:n12)

        r(n01:n11,n02:n12) = r(n01:n11,n02:n12) - alpha * q(n01:n11,n02:n12)

        res = Residut(Grid,r)

        if( res < 1e-6 ) then
            print *, "GC converged in",i,"IT with a residut of",res
            return
        end if

!        print *,res

        rhop1 = Scalar_Product(Grid,r,r)

        beta = rhop1 / rho

        rho = rhop1

      end do

      print *,"ERROR, CG not converged"

      end subroutine CG_CORE



      subroutine Matrix_A_Multiply(Grid,p,q)

      implicit none

      type (grid_type) :: Grid
      real(8), dimension(1:Grid%Nx(1),1:Grid%Nx(2)), intent(in) :: p
      real(8), dimension(1:Grid%Nx(1),1:Grid%Nx(2)), intent(out) :: q

      real(8), dimension(0:Grid%Nx(1)+1,0:Grid%Nx(2)+1) :: pl ! copy with ghost, not optimized

      real(8), dimension(1:3) :: SDx

      integer :: i,j,n

      pl=77.0 ! dummy value
      pl(1:Grid%Nx(1),1:Grid%Nx(2)) = p(1:Grid%Nx(1),1:Grid%Nx(2))

      ! Dirichlet : 0 en x=0, et 1 en x=1
      pl(        1,1:Grid%Nx(2)) = 0.0!((        1-1)*Grid%h(1))**2
      pl(Grid%Nx(1),1:Grid%Nx(2)) = 1.0!((Grid%Nx(1)-1)*Grid%h(1))**2

      ! Neumann dp/dn = 0
      pl(1:Grid%Nx(1),          0) = pl(1:Grid%Nx(1),        1)
      pl(1:Grid%Nx(1),Grid%Nx(2)+1) = pl(1:Grid%Nx(1),Grid%Nx(2))

      do n=1,2
        SDx(n)=(1.0d0/(Grid%h(n)*Grid%h(n)))
      end do

      do j=n02,n12
        do i=n01,n11
            q(i,j) = SDx(1)*(pl(i-1,j)+pl(i+1,j)-2.0d0*pl(i,j)) +    &
            &        SDx(2)*(pl(i,j-1)+pl(i,j+1)-2.0d0*pl(i,j))
        end do
      end do

      return

      end subroutine Matrix_A_Multiply


      function Residut(Grid,r)

      implicit none

      type (grid_type) :: Grid
      real(8), dimension(1:Grid%Nx(1),1:Grid%Nx(2)), intent(in) :: r
      real(8) :: Residut
      integer :: i,j
     
      Residut=0.0d0
      do j=n02,n12
        do i=n01,n11
            Residut = max ( Residut, ABS(r(i,j)) )
        end do
      end do

      return

      end function Residut


      function Scalar_Product(Grid,x1,x2)

      implicit none

      type (grid_type) :: Grid
      real(8), dimension(1:Grid%Nx(1),1:Grid%Nx(2)), intent(in) :: x1,x2
      real(8) :: Scalar_Product
      integer :: i,j

     
      Scalar_Product=0.0d0
      do j=n02,n12
        do i=n01,n11
            Scalar_Product = Scalar_Product + x1(i,j) * x2(i,j)
        end do
      end do

      return

      end function Scalar_Product



      end module CG

Type Grid :
Code:

      module typegrid

      implicit none
        type grid_type
        integer                    :: Ndim
        integer                    :: NPL
        integer, dimension(1:3)    :: Nx
        real, dimension(1:3)        :: X0
        real, dimension(1:3)        :: Lx
        real, dimension(1:3)        :: h
        real, dimension(1:100,1:3)  :: XP

        end type grid_type
      end module typegrid

The CG do not converge, however the result is really near the solution :

http://img534.imageshack.us/img534/4...twindowid0.png


As you can see, the value is a little different near x=1...

Does anyone have an idée ? It would really help me :)

moomba February 17, 2010 04:37

Problem solve, here are the sources if it can help someone.

It is a simple conjugate gradient with neumann, periodic and dirichlet boundary conditions.

Code:


! Neumann at face 1.1
BoundCG(1,1) = 1
! Dirichlet at face 1.2
 BoundCG(1,2) = -1
 PRES(Grid%Nx(1),:,:) = 0.0d0
! periodic on axe y
 BoundCG(2,1) = 12
 BoundCG(2,2) = 12

 CALL CG_CORE(Grid,PRES,BoundCG)

Code:

      module CG 
      use typegrid
 
      contains
 
 
      subroutine CG_CORE(Grid,b,BoundCG)

      use commonGC
 
      implicit none
 
      type (grid_type) :: Grid
 
      real(8), dimension(1:Grid%Nx(1),1:Grid%Nx(2),1:Grid%Nx(3)), intent(inout) :: b
      integer, dimension(1:3,1:2), intent(in) :: BoundCG
 
      real(8), dimension(1:Grid%Nx(1),1:Grid%Nx(2),1:Grid%Nx(3)) :: r,p,q,w

      real(8), dimension(1:3) :: SDx
 
      real(8) :: rho,rhop1,beta,alpha,res,a
 
      integer :: i,n
 
!    BoundCG
!    10 = simple com
!    12 = periodic com
!    1 = inflow (dP/dn = 0)
!    -1 = outflow (P = fixed)

      ! Dirichlet Boundary Conditions, set fixed value in initial guess, this values will not be overwritten during the CG process
      if(BoundCG(1,1) == -1) x(1,:,:) = b(1,:,:)
      if(BoundCG(1,2) == -1) x(Grid%Nx(1),:,:) = b(Grid%Nx(1),:,:)
      if(Grid%ndim > 1) then
        if(BoundCG(2,1) == -1) x(:,1,:) = b(:,1,:)
        if(BoundCG(2,2) == -1) x(:,Grid%Nx(2),:) = b(:,Grid%Nx(2),:)
      end if
      if(Grid%ndim > 2) then
        if(BoundCG(3,1) == -1) x(:,:,1) = b(:,:,1)
        if(BoundCG(3,2) == -1) x(:,:,Grid%Nx(3)) = b(:,:,Grid%Nx(3))
      end if
 
      CALL Matrix_A_Multiply(Grid,x,r,BoundCG) ! (x_in,x_out) --> x_out = A * x_in
 
      r(:,:,:) = b(:,:,:) - r(:,:,:)

      ! Dirichlet Boundary Conditions, set residut at 0.0 (value is already correct)
      if(BoundCG(1,1) == -1) r(1,:,:) = 0.0d0
      if(BoundCG(1,2) == -1) r(Grid%Nx(1),:,:) = 0.0d0
      if(Grid%ndim > 1) then
        if(BoundCG(2,1) == -1) r(:,1,:) = 0.0d0
        if(BoundCG(2,2) == -1) r(:,Grid%Nx(2),:) = 0.0d0
      end if
      if(Grid%ndim > 2) then
        if(BoundCG(3,1) == -1) r(:,:,1) = 0.0d0
        if(BoundCG(3,2) == -1) r(:,:,Grid%Nx(3)) = 0.0d0
      end if
 
      p(:,:,:) = 0.0d0
 
      beta = 0.0d0

      ! PRECONDITIONNEMENT
      do n=1,Grid%ndim
        SDx(n)=(1.0d0/(Grid%h(n)*Grid%h(n)))
      end do
      a=1.0d0/(-2.0d0*SDx(1)-2.0d0*SDx(2))

      w(:,:,:) = a * r(:,:,:)
 
      rho = Scalar_Product(Grid,r,w) ! return v(1)*r(1) + v(2)*r(2) + v(3)*r(3) + ... + v(n)*r(n)
 
 
      do i=1,2000
 
        p(:,:,:) = w(:,:,:) + beta * p(:,:,:)
 
        CALL Matrix_A_Multiply(Grid,p,q,BoundCG)
 
        alpha = rho / Scalar_Product(Grid,p,q)
 
        x(:,:,:) = x(:,:,:) + alpha * p(:,:,:)
 
        r(:,:,:) = r(:,:,:) - alpha * q(:,:,:)
 
        res = Residut(Grid,r)
 
        if( res < 1e-6 ) then
            print *, "GC converged in",i,"IT with a residut of",res
            b(:,:,:) = x(:,:,:)
            return
        end if
 
!        print *,res
        w(:,:,:) = a * r(:,:,:)
 
        rhop1 = Scalar_Product(Grid,r,w)
 
        beta = rhop1 / rho
 
        rho = rhop1
 
      end do
 
      print *,"ERROR, CG not converged"
 
      end subroutine CG_CORE
 
 
 
      subroutine Matrix_A_Multiply(Grid,p,q,BoundCG)
 
      implicit none
 
      type (grid_type) :: Grid
      real(8), dimension(1:Grid%Nx(1),1:Grid%Nx(2),1:Grid%Nx(3)), intent(in) :: p
      real(8), dimension(1:Grid%Nx(1),1:Grid%Nx(2),1:Grid%Nx(3)), intent(out) :: q
      integer, dimension(1:3,1:2), intent(in) :: BoundCG
 
      real(8), dimension(0:Grid%Nx(1)+1,0:Grid%Nx(2)+1,0:Grid%Nx(3)+1) :: pl ! copy with ghost, not optimized
 
      real(8), dimension(1:3) :: SDx
 
      integer :: i,j,n
 
      pl=77.0
      pl(1:Grid%Nx(1),1:Grid%Nx(2),1) = p(1:Grid%Nx(1),1:Grid%Nx(2),1)

      ! Neumann Boundaries, dP/dn = 0 , not optimized (should not copy, but simply not compute)
      if(BoundCG(1,1) == 1) pl(0,:,:) = pl(2,:,:)
      if(BoundCG(1,2) == 1) pl(Grid%Nx(1)+1,:,:) = pl(Grid%Nx(1)-1,:,:)
      if(Grid%ndim > 1) then
        if(BoundCG(2,1) == 1) pl(:,0,:) = pl(:,2,:)
        if(BoundCG(2,2) == 1) pl(:,Grid%Nx(2)+1,:) = pl(:,Grid%Nx(2)-1,:)
      end if
      if(Grid%ndim > 2) then
        if(BoundCG(3,1) == 1) pl(:,:,0) = pl(:,:,2)
        if(BoundCG(3,2) == 1) pl(:,:,Grid%Nx(3)+1) = pl(:,:,Grid%Nx(3)-1)
      end if     
 
      ! Periodic Boundaries
      if(BoundCG(1,1) == 12) pl(0,:,:) = pl(Grid%Nx(1)-1,:,:)
      if(BoundCG(1,2) == 12) pl(Grid%Nx(1)+1,:,:) = pl(2,:,:)
      if(Grid%ndim > 1) then
        if(BoundCG(2,1) == 12) pl(:,0,:) = pl(:,Grid%Nx(2)-1,:)
        if(BoundCG(2,2) == 12) pl(:,Grid%Nx(2)+1,:) = pl(:,2,:)
      end if
      if(Grid%ndim > 2) then
        if(BoundCG(3,1) == 12) pl(:,:,0) = pl(:,:,Grid%Nx(3)-1)
        if(BoundCG(3,2) == 12) pl(:,:,Grid%Nx(3)+1) = pl(:,:,2)
      end if
 
      do n=1,2
        SDx(n)=(1.0/(Grid%h(n)*Grid%h(n)))
      end do
 
      do j=1,Grid%Nx(2)
        do i=1,Grid%Nx(1)
            q(i,j,1)=SDx(1)*(pl(i-1,j,1)+pl(i+1,j,1)-2.0*pl(i,j,1)) +    &
            &        SDx(2)*(pl(i,j-1,1)+pl(i,j+1,1)-2.0*pl(i,j,1))
        end do
      end do

      ! Dirichelt boundaries, set q at 0 (to not alter the x fixed value)
      if(BoundCG(1,1) == -1) q(1,:,:) = 0.0d0
      if(BoundCG(1,2) == -1) q(Grid%Nx(1),:,:) = 0.0d0
      if(Grid%ndim > 1) then
        if(BoundCG(2,1) == -1) q(:,1,:) = 0.0d0
        if(BoundCG(2,2) == -1) q(:,Grid%Nx(2),:) = 0.0d0
      end if
      if(Grid%ndim > 2) then
        if(BoundCG(3,1) == -1) q(:,:,1) = 0.0d0
        if(BoundCG(3,2) == -1) q(:,:,Grid%Nx(3)) = 0.0d0
      end if
 
      return
 
      end subroutine Matrix_A_Multiply
 
 
      function Residut(Grid,r)
 
      implicit none
 
      type (grid_type) :: Grid
      real(8), dimension(1:Grid%Nx(1),1:Grid%Nx(2),1:Grid%Nx(3)), intent(in) :: r
      real(8) :: Residut
      integer :: i,j,k
     
      Residut=0.0d0
      do k=1,Grid%Nx(3)
        do j=1,Grid%Nx(2)
            do i=1,Grid%Nx(1)
              Residut = max ( Residut, ABS(r(i,j,k)) )
            end do
        end do
      end do
 
      return
 
      end function Residut
 
 
      function Scalar_Product(Grid,x1,x2)
 
      implicit none
 
      type (grid_type) :: Grid
      real(8), dimension(1:Grid%Nx(1),1:Grid%Nx(2),1:Grid%Nx(3)), intent(in) :: x1,x2
      real(8) :: Scalar_Product
      integer :: i,j,k
     
      Scalar_Product=0.0d0
      do k=1,Grid%Nx(3)
        do j=1,Grid%Nx(2)
            do i=1,Grid%Nx(1)
              Scalar_Product = Scalar_Product + x1(i,j,k) * x2(i,j,k)
            end do
        end do
      end do
 
      return
 
      end function Scalar_Product
 
 
 
      end module CG

Code:

module commonGC

  implicit none

  real(8), allocatable, dimension(:,:,:), save :: x

end module commonGC

If used, you need to allocate x, and set it to 0.0 before the first time you use the routine, because the CG use the last solution calculated at the previous time as initial guess.

This program is not very efficient, but it is simple and a may help to start with CG.

With my best regards


All times are GMT -4. The time now is 16:56.