C----------------------------------------------------------------------
C CG:  Standard Conjugate Gradient Method without Preconditioning
C----------------------------------------------------------------------
      subroutine mpi_cg(X, RHS, m,n,s,e,w,
     >              D, G, H, wrk, buf1,buf2,
     >              me,north,south,east,west,
     >              eps, Ceps, iter, arfac, lx, ly, bufsize)
      implicit none
      integer*4 iter,lx,ly,me,bufsize,north,south,east,west
      double precision X,RHS,d,g,h,m,n,s,e,w,
     >                 wrk,buf1,buf2,
     >                 arfac,eps,Ceps,TAU,BETA,delta,delta0,delta1,
     >                 mpi_dot2
      DIMENSION X(lx,ly), RHS(lx,ly),
     >          D(lx,ly), G(lx,ly), H(lx,ly),
     >          m(lx,ly), 
     >          n(lx,ly), s(lx,ly),
     >          e(lx,ly), w(lx,ly),
     >          wrk(0:lx+1,0:ly+1),buf1(bufsize),buf2(bufsize)


      iter = 0
      Ceps = eps
C ---------------------------------------------- initial defect
      call nullv(X,lx,ly)
      call matmul(H, m,n,s,e,w, X, wrk,buf1,buf2,
     >            me,north,south,east,west,lx,ly,bufsize)
      call vvsub(G, H, RHS, lx,ly)
C ---------------------------------------------- initial quasidefect
      call vminv(D, G, lx,ly)
      delta0 = mpi_dot2(G, G, lx,ly)
      Ceps = delta0 * eps
      delta = delta0
      if (me.eq.0) write(*,*) '---------------- delta0:', delta0
      IF (delta0 .LT. eps) GO TO 10

 20   continue

      call matmul(H, m,n,s,e,w, D, wrk,buf1,buf2,
     >            me,north,south,east,west,lx,ly,bufsize)
      iter = iter + 1
c      IF (iter .GT. 100) GOTO 10
      TAU = delta0 / mpi_dot2(D, H, lx, ly)
c ---------------------------------------------- next defect X
      call vpabsv(X, TAU, D, lx, ly)
C ---------------------------------------------- next defect G
      call vpabsv(G, TAU, H, lx, ly)
C ---------------------------------------------- next search direction H
      delta1 = mpi_dot2(G, G, lx, ly)
       if (mod(iter,200) .eq. 0) then
      if (me.eq.0) write(*,*) '------------------ delta1:', iter,delta1
       endif
      IF (delta1 .LT. Ceps) GO TO 10
      beta = delta1 / delta0
      delta0 = delta1
      call vminsv(D, G, beta, D, lx, ly)
      GO TO 20

 10   continue

      if (iter .gt. 0) then
          arfac = (delta1/delta) ** (1.0d0/(2.0d0 * iter))
      else
          arfac = DSQRT(delta)
      endif
      if (me.eq.0) write(*,*) '------------------ cg_iter:', iter
      return
      end
c --------------------------------------------------------------------
      subroutine print(wrk,me,lx,ly,message)
      implicit none

      integer i,j,me,lx,ly
      double precision wrk(lx,ly)
      character*20 message

       write(*,*)
       write(*,*) message,'(',me,')'
       do i=1,lx
          write(*,*) (wrk(i,j), j=1,ly)
       enddo

      return
      end