CFD Online Discussion Forums - CG, BICGSTAB(2) : problem with matrix operation and boundary conditions

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)

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:

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 :)

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