c****************************************************************
c     Iterative solver CG
c****************************************************************
c     Authors: P. C. Chaumet and A. Rahmani
c     Date: 27/10/2010
c     
c     Purpose: iterative solver for linear system Ax=b. There is no
c     condition on the matrix. Notice that the product A x is provided
c     by the user.
c     
c     Reference: if you use this routine in your research, please
c     reference, as appropriate: P. C. Chaumet and A. Rahmani, Efficient
c     discrete dipole approximation for magneto-dielectric scatterers
c     Opt. Lett. 34, 917 (2009). J. Tang, Y. Shen, Y. Zheng, and D. Qiu,
c     Coastal Eng. 51, 143 (2004).
c
c     license: GNU GPL
c     We cannot guarantee correctness of the programs...
c****************************************************************
c****************************************************************
c****************************************************************
c     Main program and example. The main program if provided for testing
c     purposes and should be commented out to use only the routine ZCG      
      implicit none
      integer LDA,NDIM,NLAR,NOU,NSTAT,NLOOP,MAXIT,i,j
      parameter(lda=100,nlar=7)
      double precision NORM,TOL,TOLE
      double complex XI(lda),XR(lda),B(lda),WRK(lda,nlar),ALPHA,RHO ,XXI
     $     ,mat(lda,lda),icomp


c     We want to solve Ax=b
c     mat:  matrix A


c     lda: Input: leading dimension array of the matrix

c     ndim: Input: dimension  of the matrix: ndim.le.lda

c     nlar: Input: size of the work vector: nlar.ge.9

c     MAXIT: Input: Maximum of iteration for the iterative solver

c     NLOOP: Output: number of iteration for the iterative solver Should
c     be initialize to 0.

c     ncompte: number of Ax products.

c     STATUS: Output: if STATUS.eq.-4 the maximum number of iterations
c     has been reached else STATUS.eq.1 the solution has been reached.

c     tol: Input: tolerance requested by the user. At the end we have:
c     r=||Ax-b||/||b||<tol

c     tole: output: relative error obtained

c     b: Input: right-hand side of Ax=b

c     norm: Output: norm of b

c     xi: Input: initial guess; Output: solution of the system

c     xr: Input: xr = A xi

c     NOU: Local integer used by ZCG. Should be initialized to 0.

c     ALPHA,RHO,XXI complex needs for ZCG

c     WRK: local array used by QMR.
      
c     Example 
      ndim=100
      icomp=(0.d0,1.d0)
      do i=1,ndim
         do j=1,ndim
            mat(i,j)=0.d0
         enddo
         mat(i,i)=(1.d0,0.d0)*dble(i)*icomp
         b(i)=(1.d0,0.d0)
         xi(i)=(0.d0,0.d0)
      enddo

      MAXIT=100
      tol=1.d-4
      tole=1.d0
      nloop=0
      nou=0 

 10   call ZCG(XI,XR,B,NORM,WRK,ALPHA,RHO,XXI,LDA,NDIM ,NLAR,NOU
     $     ,NSTAT,NLOOP,MAXIT,TOLE,TOL)

      do i=1,ndim
         xr(i)=0.d0
         do j=1,ndim
            xr(i)=xr(i)+mat(i,j)*xi(j)
         enddo
      enddo

      if (nstat.eq.0) goto 10

      DO I=1,NDIM
         WRITE(*,*) 'SOL',XI(I)
      ENDDO
      write(*,*) 'tolerance',tole,'<',tol,'NLOOP',NLOOP,'MAXIT',MAXIT
     $     ,nstat
      
      end
c****************************************
      subroutine ZCG(XI,XR,B,NORM,WRK,ALPHA,RHO0,XXI,LDA,NDIM,NLAR,NOU
     $     ,NSTAT,NLOOP,MAXIT,TOLE,TOL)
      IMPLICIT NONE

      INTEGER I,NLOOP,MAXIT

C     constant
      DOUBLE PRECISION ONE,ZZERO,ZONE,MZONE
      PARAMETER (ONE=1.0D0,ZZERO= (0.0D0,0.0D0),ZONE= (1.0D0,0.0D0)
     $     ,MZONE= (-1.0D0,0.0D0))
      

C     array
      INTEGER LDA,NDIM,NLAR,NOU,NSTAT
      DOUBLE COMPLEX B(LDA),WRK(LDA,NLAR),XR(LDA),XI(LDA)

c     Local Scalars
      DOUBLE PRECISION TOL,NORM,TOLE
      DOUBLE COMPLEX ALPHA,BETA,DELTA,RHO,RHO0,XXI,TMP1,TMP2

c     enetr with nou=0 and nloop=0

      IF (NOU.EQ.1) GOTO 10
      IF (NOU.EQ.2) GOTO 20
      IF (NOU.EQ.3) GOTO 30
      IF (NOU.EQ.4) GOTO 40
c     NOU=4 si change tolerance en cours de route
      
*  1. r=(b-Ax) CALCUL DU RESIDU
      NOU=1  
*     Loop
      NSTAT=0
c     initial guess
      DO I=1,NDIM
         WRK(I,1)=XI(I)
      ENDDO

     
C     CALCUL Ax AL'EXTERIEUR ET RENVOIE DANS XR
      RETURN
C     CALCUL b-Ax
 10   NORM=0.D0
      DO I=1,NDIM
C     RK-1
         WRK(I,2)=B(I)-XR(I)
C     INITIALISE XK
         NORM=NORM+B(I)*DCONJG(B(I))
      ENDDO
      NORM=dsqrt(NORM)

*     2.  p=r
      DO I=1,NDIM
         WRK(I,3)=WRK(I,2)
      ENDDO

*     3. rdotr=dot(r,r)
      RHO0=0.d0
      DO I=1,NDIM
         RHO0=RHO0+WRK(I,2)*WRK(I,2)
      ENDDO

*  4. w=Ap
      NOU=2
      DO I=1,NDIM
         XI(I)=WRK(I,3)
      ENDDO
      RETURN

 20   DO I=1,NDIM
         WRK(I,5)=XR(I)
      ENDDO

*     5. xi=dot(p,w)     
      XXI=0.d0
      DO I=1,NDIM
         XXI=XXI+WRK(I,3)*WRK(I,5)
      ENDDO
c      write(*,*) 'XXI',XXI
      if (XXI.eq.0.d0) then
         write(*,*) 'probleme avec XXI',XXI
         stop
      endif

*     Loop
     
 50   NLOOP=NLOOP+1

      IF (NLOOP.GE.MAXIT) then
         NSTAT=-4
         NOU=0
         DO I=1,NDIM
            XI(I)=WRK(I,1)
         ENDDO
         RETURN         
      ENDIF
*     6. alpha=rdotr/xi      
      ALPHA=RHO0/XXI

*     7. x=x+alpha*p
      DO I=1,NDIM
         WRK(I,1)= WRK(I,1)+ALPHA*WRK(I,3)
      ENDDO

*     8. r=r-alpha*w      
      DO I=1,NDIM
         WRK(I,2)= WRK(I,2)-ALPHA*WRK(I,5)
      ENDDO

*     9. check stopping criterion
      TMP1=ZZERO
      DO I=1,NDIM
         TMP1=TMP1+WRK(I,2)*DCONJG(WRK(I,2))
      ENDDO
      TOLE=dsqrt(cdabs(TMP1))/NORM
      IF (TOLE.LE.TOL) then
         NSTAT=1
         DO I=1,NDIM
            XI(I)=WRK(I,1)
         ENDDO
         NOU=4
         RETURN         
      ENDIF

*     10. s=Ar
 40   NOU=3
      DO I=1,NDIM
         XI(I)=WRK(I,2)
      ENDDO
      RETURN

 30   DO I=1,NDIM
         WRK(I,4)=XR(I)
      ENDDO

*     11. rdotr=dot(r,r)
      RHO=0.d0
      DO I=1,NDIM
         RHO=RHO+WRK(I,2)*WRK(I,2)
      ENDDO

*     12. delta=dot(r,s)
      DELTA=0.d0
      DO I=1,NDIM
         DELTA=DELTA+WRK(I,2)*WRK(I,4)
      ENDDO

*     13. beta=rdotr/rdotr0
      BETA=RHO/RHO0
      RHO0=RHO

*     14. p=r+beta*p
      DO I=1,NDIM
         WRK(I,3)=WRK(I,2)+BETA*WRK(I,3)
      ENDDO

*     15. w=s+beta*w
      DO I=1,NDIM
         WRK(I,5)=WRK(I,4)+BETA*WRK(I,5)
      ENDDO

*     16. xi=delta-beta^2*xi      
      XXI=DELTA-BETA*BETA*XXI

      GOTO 50


      END
