c****************************************************************
c     Iterative solver QMR
c****************************************************************
c     Authors: P. C. Chaumet and A. Rahmani
c     Date: 04/02/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 QMR

      implicit none
      integer LDA,NDIM,NLAR,NOU,ITNO,MAXIT,i,j,NT,STATUS,STEPERR,ncompte
      parameter(lda=200,nlar=9)
      double precision NORM,TOL,TOLE
      double complex XS(lda),XI(lda),XR(lda),B(lda),WRK(lda,nlar)
     $     ,mat(lda,lda),icomp
      DOUBLE COMPLEX GAMMA,EPSILON ,KAPPA,KSI,LAMBDA,MU,RHO,TAU,THETA
     $     ,DOTS(4)
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     ITNO: 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 GPBICG
c     STATUS=1 the requested tolerance has been reached or the maximum number of
c     iterations has been reached

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

c     NT: output. If NT=1 compute A x in the main program if NT=2
c     compute A^t x

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 QMR. Should be initialized to 0.

c     GAMMA,EPSILON ,KAPPA,KSI,LAMBDA,MU,RHO,TAU,THETA ,DOTS(4) local
c     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)=1.d0
         enddo
         mat(i,i)=(1.d0,0.d0)*dble(i)*icomp
         b(i)=(1.d0,0.d0)*icomp
         xi(i)=0.d0
      enddo

      tol=1.d-6
      tole=1.d0
      nou=0
      MAXIT=100
      ncompte=0

 10   call PIMZQMR(XS,XI,XR,B,WRK,NORM,LDA,NDIM,NLAR,LAMBDA,KAPPA ,THETA
     $     ,GAMMA,KSI,RHO,EPSILON,MU,TAU,DOTS,NOU,NT,ITNO,MAXIT,TOLE
     $     ,TOL ,STATUS ,STEPERR)

      if (status.lt.0) then
         write(*,*) 'stop nstat',status,STEPERR
         stop
      endif
      ncompte=ncompte+1

      if (nt.eq.1) then
         do i=1,ndim
            xr(i)=0.d0
            do j=1,ndim
               xr(i)=xr(i)+mat(i,j)*xi(j)
            enddo
         enddo      
      elseif (nt.eq.2) then
         do i=1,ndim
            xr(i)=0.d0
            do j=1,ndim
               xr(i)=xr(i)+(mat(j,i))*xi(j)
            enddo
         enddo
      endif

      if (status.ne.1) goto 10

      if (STEPERR.eq.0) then
         write(*,*) 'ITNO has reached MAXIT',ITNO,MAXIT
      endif

      DO I=1,NDIM
         WRITE(*,*) 'SOL',XS(I)
      ENDDO

c     relative error
      write(*,*) 'relative error',tole,'<',tol
      write(*,*) 'number of iteration',ITNO,MAXIT
      write(*,*) 'number of product Ax',ncompte


      end
c****************************************
      SUBROUTINE PIMZQMR(XS,XI,XR,B,WRK,NORM,LDA,NDIM,NLAR,LAMBDA,KAPPA
     $     ,THETA,GAMMA,KSI,RHO,EPSILON,MU,TAU,DOTS,NOU,NT,ITNO,MAXIT
     $     ,TOLE ,TOL ,STATUS ,STEPERR)

      IMPLICIT NONE

*     .. Array Arguments ..
      DOUBLE COMPLEX B(LDA),WRK(LDA,NLAR),XI(LDA),XR(LDA),XS(LDA)
*     ..
*     .. Local Scalars ..
      DOUBLE COMPLEX ABSGAMMA2,ABSTAU02,DEN,GAMMA,EPSILON,EPSILON0,
     +     GAMMA0,KAPPA,KAPPA0,KSI,KSI0,LAMBDA,LAMBDA0,MU,MU0,
     +     RHO,RHO0,TAU,TAU0,THETA,TMP1
      DOUBLE PRECISION NORM, TOL, TOLE
      INTEGER NDIM, NLAR, NOU, I, ITNO, LDA ,MAXIT,N, STATUS,STEPERR,NT
*     ..
*     .. Local Arrays ..
      DOUBLE COMPLEX DOTS(4)
*     ..

      if (nou.eq.0) goto 10
      if (nou.eq.1) goto 20
      if (nou.eq.2) goto 30
      if (nou.eq.3) goto 40
      if (nou.eq.4) goto 50
      if (nou.eq.5) goto 100

*     1. lambda=1, kappa=-1, theta=-1
 10   LAMBDA = (1.d0,0.d0)
      KAPPA = -(1.d0,0.d0)
      THETA = -(1.d0,0.d0)
      NORM=0.d0
      DO I=1,NDIM
         NORM=NORM+B(I)*DCONJG(B(I))
      ENDDO
      NORM=dsqrt(NORM)

*     Loop
      STATUS = 0
      STEPERR = -1
      ITNO=0

*     2. wtilde=vtilde=r=b-Ax
*     r=b-Ax

c     A*x=wrk(3)
      nou=1
      nt=1
      do i=1,ndim
         xs(i)=xi(i)
      enddo

c     compute A*xi

      Return

 20   do i=1,ndim
         wrk(i,3)=xr(i)
         wrk(i,1)=B(i)-wrk(i,3)
         wrk(i,7)=wrk(i,1)
         wrk(i,8)=wrk(i,1)
      enddo	

*     3. p=q=d=s=0
      do i=1,ndim
	 wrk(i,2)=0.d0
	 wrk(i,4)=0.d0
	 wrk(i,5)=0.d0
	 wrk(i,6)=0.d0
      enddo

*     4. gamma=||vtilde||_{2}, ksi=||wtilde||_{2},
*     rho=wtilde^{T}vtilde, epsilon=(A^{T}wtilde)^{T}vtilde, mu=0

      DOTS(1) = 0.d0
      DOTS(2) = 0.d0
      DOTS(3) = 0.d0

      do i=1,ndim
         DOTS(1) = DOTS(1)+wrk(i,7)*dconjg(wrk(i,7))
         DOTS(2) = DOTS(2)+wrk(i,8)*dconjg(wrk(i,8))
         DOTS(3) = DOTS(3)+wrk(i,7)*wrk(i,8)
      enddo
*     Compute A^{T}wtilde
c     CALL TMATVEC(WRK(IWTILDE),WRK(IATWTILDE),IPAR)
      DO I=1,NDIM
         XI(I)= WRK(I,8)
      ENDDO
      
      nt=2
      nou=2
      return


 30   DOTS(4) = 0.d0
      do i=1,ndim
	 wrk(i,9)=xr(i)
         DOTS(4) = DOTS(4)+wrk(i,7)*wrk(i,9)
      enddo

*     Accumulate simultaneously partial inner-products
c     CALL PDZSUM(4,DOTS)

      GAMMA = CDSQRT(DOTS(1))
      KSI = CDSQRT(DOTS(2))
      RHO = DOTS(3)
      EPSILON = DOTS(4)
      MU = 0.d0
*     5. tau=epsilon/rho
      IF (RHO.EQ.0.d0) THEN
         ITNO = 0
         STATUS = -3
         STEPERR = 5
         GO TO 200
      END IF

      TAU = EPSILON/RHO
 100  ITNO=ITNO+1
      
*     6. p=1/gamma*vtilde-mu*p

      IF (GAMMA.EQ.0.d0) THEN
         STATUS = -3
         STEPERR = 6
         GO TO 200

      END IF
      do i=1,ndim
         wrk(i,2)=wrk(i,7)/gamma-mu*wrk(i,2)
      enddo

*     7. q=1/ksi*A^{T}wtilde-(gamma*mu)/ksi*q
      IF (KSI.EQ.0.d0) THEN
         STATUS = -3
         STEPERR = 7
         GO TO 200
      END IF

      do i=1,ndim
         wrk(i,4)=(wrk(i,9)-gamma*mu*wrk(i,4))/ksi
      enddo
      
*     8. vtilde=Ap-tau/gamma*vtilde
      DO I=1,NDIM
         XI(I)= WRK(I,2)
      ENDDO
      
      nt=1
      nou=3
      return

 40   do i=1,ndim
         wrk(i,3)=xr(i)
         wrk(i,7)=wrk(i,3)-tau/gamma*wrk(i,7)
      enddo

*     9. wtilde=q-tau/ksi*wtilde

      do i=1,ndim
         wrk(i,8)=wrk(i,4)-tau/ksi*wrk(i,8)
      enddo

*     11. gamma=||vtilde||_{2}, ksi=||wtilde||_{2},
*     rho=wtilde^{T}vtilde, epsilon=(A^{T}wtilde)^{T}vtilde

      DOTS(1) = 0.d0
      DOTS(2) = 0.d0
      DOTS(3) = 0.d0

      do i=1,ndim
         DOTS(1) = DOTS(1)+wrk(i,7)*dconjg(wrk(i,7))
         DOTS(2) = DOTS(2)+wrk(i,8)*dconjg(wrk(i,8))
         DOTS(3) = DOTS(3)+wrk(i,7)*wrk(i,8)
      enddo

*     Compute A^{T}wtilde
      DO I=1,NDIM
         XI(I)= WRK(I,8)
      ENDDO
      
      nt=2
      nou=4
      return
 50   DOTS(4) = 0.d0
      do i=1,ndim
	 wrk(i,9)=xr(i)
         DOTS(4) = DOTS(4)+wrk(i,7)*wrk(i,9)
      enddo

*     Accumulate simultaneously partial inner-products
c     CALL PDZSUM(4,DOTS)

      GAMMA0 = GAMMA
      GAMMA = CDSQRT(DOTS(1))
      KSI0 = KSI
      KSI = CDSQRT(DOTS(2))
      RHO0 = RHO
      RHO = DOTS(3)
      EPSILON0 = EPSILON
      EPSILON = DOTS(4)
*     12. mu=(gamma0*ksi0*rho)/(gamma*tau*rho0)

      DEN = GAMMA*TAU*RHO0
      IF (DEN.EQ.0.d0) THEN
         STATUS = -3
         STEPERR = 12
         GO TO 200

      END IF
      MU0 = MU
      MU = (GAMMA0*KSI0*RHO)/DEN

*     13. tau=epsilon/rho-gamma*mu
      IF (RHO.EQ.0.d0) THEN
         STATUS = -3
         STEPERR = 13
         GO TO 200

      END IF
      TAU0 = TAU
      TAU = EPSILON/RHO - GAMMA*MU
*     14. theta=(|tau0|^2*(1-lambda))/(lambda*|tau|^2+|gamma|^2)

      ABSTAU02 = CDABS(TAU0)**2.D0
      ABSGAMMA2 = CDABS(GAMMA)**2.D0
      DEN = LAMBDA*ABSTAU02 + ABSGAMMA2
      IF (DEN.EQ.0.d0) THEN
         STATUS = -3
         STEPERR = 14
         GO TO 200

      END IF
      THETA = (ABSTAU02* ((1.d0,0.d0)-LAMBDA))/DEN

*     15. kappa=(-gamma0*CONJG(tau0)*kappa0)/(gamma0*|tau|^2+|gamma|^2)
      KAPPA0 = KAPPA
      KAPPA = - (GAMMA0*DCONJG(TAU0)*KAPPA0)/DEN

*     16. lambda=(lambda0*|tau0|^2)/(gamma0*|tau|^2+|gamma|^2)
      LAMBDA0 = LAMBDA
      LAMBDA = LAMBDA0*ABSTAU02/DEN

*     17. d=theta*d+kappa*p

      do i=1,ndim
         wrk(i,5)=theta*wrk(i,5)+kappa*wrk(i,2)
      enddo

*     18. s=theta*s+kappa*A*p

      do i=1,ndim
         wrk(i,6)=theta*wrk(i,6)+kappa*wrk(i,3)
      enddo

*     19. x=x+d

      do i=1,ndim
         xs(i)=xs(i)+wrk(i,5) 
      enddo

*     20. r=r-s

      do i=1,ndim
         wrk(i,1)=wrk(i,1)-wrk(i,6)
      enddo

c     criterion to stop
      TMP1=0.d0
      DO I=1,NDIM
         TMP1=TMP1+WRK(I,1)*DCONJG(WRK(I,1))
      ENDDO
      TOLE=dsqrt(cdabs(TMP1))/NORM

      IF (TOLE.LE.TOL) then
         NOU=5
         STATUS=1
         RETURN         
      ENDIF

      IF (ITNO.GT.MAXIT) THEN
         STATUS = 1
         STEPERR = 0
         RETURN
      END IF
      GOTO 100

 200  RETURN

      END

