c****************************************************************
c     Iterative solver TFQMR
c****************************************************************
c     Authors: P. C. Chaumet and A. Rahmani
c     Date: 26/10/2010
c     
c     Purpose: iterative solver for linear system Ax=b. There is no
c     condition on the matrix. Notice that the products A x and At x are
c     provided 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). R. D. Da Cunha and T. Hopkins,
c     Appl. Numer. Math. 19, 33 (1995).

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 TFQMR
      implicit none
      integer LDA,NDIM,NLAR,NOU,STATUS,STEPERR,NLOOP,MAXIT,i,j,IM
     $     ,ncompte
      parameter(lda=200,nlar=9)
      double precision NORM,TOL,TOLE
      double complex XI(lda),XR(lda),xs(lda),B(lda),WRK(lda,nlar)
     $     ,mat(lda,lda),ALPHA,BETA,RHO,TAU,THETA,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.lt.0 a problem occured in TFQMR
c     STATUS=1 the requested tolerance has been reached or the maximum number of
c     iterations has been reached

c     STEPERR: Output: if STEPERR.gt.0: indicates where the problem
c     occurs in TFQMR. STEPERR=0 the maximum number of iterations has
c     been reached.  STEPERR=-1 routine completed without error

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;

c     xr: Input: xr = A xi

c     xs: Output: solution of the system

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

c     ALPHA,BETA,RHO,TAU,THETA local complex needs for GPBICG

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
      
      tol=1.d-9
      tole=1.d0
      nloop=0
      nou=0
      MAXIT=50
      ncompte=0
      
 10   call TFQMR(Xs,Xi,XR,B,lda,ndim,nlar,nou,WRK,NLOOP,MAXIT,TOL,NORM
     $     ,ALPHA,BETA,RHO,TAU,THETA,STATUS,STEPERR)
      ncompte=ncompte+1
      if (STATUS.lt.0) then
         write(*,*) 'stop nstat',STATUS,STEPERR
         stop
      endif
      
      do i=1,ndim
         xr(i)=0.d0
         do j=1,ndim
            xr(i)=xr(i)+mat(i,j)*xi(j)
         enddo
      enddo
      
      if (STATUS.ne.1) goto 10
  
      
      DO I=1,NDIM
         WRITE(*,*) 'solution',XS(I)
      ENDDO
      
c     calcul erreur relative
      tole=0.d0
      do i=1,ndim
         tole=tole+cdabs(xr(i)-b(i))**2.d0
      enddo
      tole=dsqrt(tole)/norm
      write(*,*) 'convergence',STEPERR
      write(*,*) 'relative error',tole,'<',tol
      write(*,*) 'number of iteration',NLOOP,MAXIT
      write(*,*) 'number of product Ax',ncompte


  
      end

c*************************************************
      SUBROUTINE TFQMR(X,Xi,XR,B,lda,ndim,nlar,nou,WRK,ITNO,MAXIT ,TOL
     $     ,NORM,ALPHA,BETA,RHO,TAU,THETA,STATUS,STEPERR)
      IMPLICIT NONE

*     .. Array Arguments ..
      INTEGER ii,nou,ndim,IM,ITNO,LDA,MAXIT,STATUS,STEPERR,ILOOP,nlar
      DOUBLE COMPLEX B(lda),WRK(lda,nlar),X(lda),xi(lda),xr(lda)
 
*     .. Local Scalars ..
      DOUBLE COMPLEX ALPHA,BETA,C,ETA,RHO,RHO0,SIGMA,TAU,THETA
      DOUBLE PRECISION TOL,NORM,KAPPA

      IF (NOU.EQ.1) GOTO 10
      IF (NOU.EQ.2) GOTO 20
      IF (NOU.EQ.3) GOTO 30
      IF (NOU.EQ.4) GOTO 40

      STATUS = 0
      STEPERR = -1

c     1:u
c     2:w
c     3:r
c     4:v
c     5:d
c     6:rtilde
c     7: u+1
c     8: Au
c
c     1
c     Compute norm and Ax0
      NORM=0.d0
      do ii=1,ndim
         x(ii)=xi(ii)
         NORM=NORM+cdabs(B(ii))**2.d0
      enddo
      NORM=dsqrt(NORM)
      nou=1
      return

c     initialisze w0=u0=r0=b-Ax0
 10   do ii=1,ndim
         WRK(ii,3)=B(ii)-xr(ii)
         WRK(ii,2)=WRK(ii,3)
         WRK(ii,1)=WRK(ii,3)
         xi(ii)=WRK(ii,1)
      enddo
      nou=2
      return
c     initialize v0 et d0
 20   do ii=1,ndim
         WRK(ii,4)=xr(ii)
         WRK(ii,5)=0.d0
         WRK(ii,8)=WRK(ii,4)
      enddo
c     2
c     initialize tau,theta,eta
      TAU=(1.d0,0.d0)*NORM
      THETA=0.d0
      ETA=0.d0

c     3
c     rtilde
      RHO=0.d0
      do ii=1,ndim
         WRK(ii,6)=WRK(ii,3)
         RHO=RHO+dconjg(WRK(ii,6))*WRK(ii,3)
      enddo

c     4
c     begin the loop
      ITNO=-1
 100  ITNO=ITNO+1

c     5
      if (mod(ITNO,2).eq.0) then
c     6
c     compute alpha
         ALPHA=0.d0
         do ii=1,ndim
            ALPHA=ALPHA+dconjg(WRK(ii,6))*WRK(ii,4)
         enddo
         if (ALPHA.eq.0.d0) then
            STATUS=-1
            STEPERR=1
            return 
         endif

         ALPHA=RHO/ALPHA
c     7
c     compute of u+1
         do ii=1,ndim
            WRK(ii,7)=WRK(ii,1)-ALPHA*WRK(ii,4)
         enddo
      endif
c     8
c     9
c     compute w        
      do ii=1,ndim
         WRK(ii,2)=WRK(ii,2)-ALPHA*WRK(ii,8)
      enddo

c     10
      do ii=1,ndim
         WRK(ii,5)=WRK(ii,1)+THETA*THETA/ALPHA*ETA*WRK(ii,5)
      enddo
c     11
c     compute THETA and C
      THETA=0.d0
      do ii=1,ndim
         THETA=THETA+dconjg(WRK(ii,2))*WRK(ii,2)
      enddo
      if (TAU.eq.0.d0) then
         STATUS=-1
         STEPERR=2
         return 
      endif
      THETA=cdsqrt(THETA)/TAU
      if (1.d0+THETA*THETA.eq.0.d0) then
         STATUS=-1
         STEPERR=3
         return 
      endif
      C=(1.d0,0.d0)/cdsqrt(1.d0+THETA*THETA)
c     12
c     compute TAU and ETA
      TAU=TAU*THETA*C
      ETA=C*C*ALPHA
c     13
c     compute X
      do ii=1,ndim
         x(ii)=x(ii)+ETA*WRK(ii,5)
      enddo
c     stop the loop?
      if (dsqrt(dble(ITNO+1))*cdabs(TAU)/NORM.le.TOL) then
         STATUS=1
         do ii=1,ndim
            xi(ii)=x(ii)
         enddo
         nou=4
         return
 40   endif


c     14
      if (mod(ITNO,2).eq.1) then
c     15
         RHO0=RHO
         RHO=0.d0
         do ii=1,ndim
            RHO=RHO+dconjg(WRK(ii,6))*WRK(ii,2)
         enddo
         BETA=RHO/RHO0
c     16
         do ii=1,ndim
            WRK(ii,7)=WRK(ii,2)+BETA*WRK(ii,1)
            WRK(ii,9)=WRK(ii,8)
         enddo
      endif

c     compute Au
      do ii=1,ndim
         xi(ii)=WRK(ii,7)
         WRK(ii,1)=WRK(ii,7)
      enddo
      nou=3
      return
         
 30   do ii=1,ndim
         WRK(ii,8)=xr(ii)
      enddo
      
      if (mod(ITNO,2).eq.1) then
         do ii=1,ndim
            WRK(ii,4)=WRK(ii,8)+BETA*(WRK(ii,9)+BETA*WRK(ii,4))
         enddo
      endif

      if (ITNO.le.MAXIT) goto 100

      STATUS=1
      STEPERR=0
      do ii=1,ndim
         xi(ii)=x(ii)
      enddo

      END
