pim.f

利用离散偶极近似方法计算散射体的电磁场。 DDA 方法

*  Check consistency of stop types
      IF ((STOPTYPE.NE.1) .AND. (STOPTYPE.NE.2) .AND.
     +    (STOPTYPE.NE.7)) THEN
          ITNO = 0
          STATUS = -6
          STEPERR = 0
          GO TO 9999

      END IF

*  Does not need conversion Y=Q2X for residual
      CNVRTX = 0

*  Set indices for mapping local vectors into wrk
      IW = 1
      IK = IW + LOCLEN
      IZ = IK + LOCLEN
      IR = IZ + LOCLEN
      IXOLD = IR + LOCLEN

*  Set rhs of stopping criteria
      RHSSTOP = SCSETRHSSTOP(B,WRK(IR),EPSILON,IPAR,PRECONL,PSCNRM)

*  1. Set parameters for iteration
      IF ((SPAR(3).EQ.ZERO) .AND. (SPAR(4).EQ.ZERO) .AND.
     +    (SPAR(5).EQ.ZERO) .AND. (SPAR(6).EQ.ZERO)) THEN
          STATUS = -7
          STEPERR = 1
          GO TO 9999

      ELSE IF (SPAR(5).EQ.SPAR(6)) THEN
*     Eigenvalues are contained in the interval [SPAR(3),SPAR(4)] on
*     the real axis:
*         sigma=(dpar(4)-dpar(3))/(2-dpar(4)-dpar(3))
*         gamma=2/(2-dpar(4)-dpar(3))
          SIGMA = (SPAR(4)-SPAR(3))/ (TWO-SPAR(4)-SPAR(3))
          SIGMASQ = SIGMA*SIGMA
          GAMMA = TWO/ (TWO-SPAR(4)-SPAR(3))

      ELSE IF (SPAR(3).EQ.SPAR(4)) THEN
*     Eigenvalues are contained in the interval [SPAR(5),SPAR(6)] on
*     the imaginary axis:
*         sigma^2=-max(dpar(5),dpar(6))
*         gamma=1
          SIGMASQ = -MAX(SPAR(5),SPAR(6))
          GAMMA = ONE

      ELSE
*     Eigenvalues are complex and contained in the box
*     SPAR(3)<= Real(e) <= SPAR(4) and SPAR(5)<= Imag(e) <= SPAR(6).
*     Compute the minimum bounding ellipse that circumscribes the box;
*     this is defined by its axes a=sqrt(2)*(dpar(4)-dpar(3))/2 (along
*     the real axis) and b=sqrt(2)*(dpar(6)-dpar(5))/2 (along the
*     imaginary axis). The center of the ellipse is d.
*         sigma^2=(a^2+b^2)/(1-d)^2
*         gamma=1/(1-d)
          LENGTHR = (SPAR(4)-SPAR(3))/TWO
          LENGTHI = (SPAR(6)-SPAR(5))/TWO
          AXISRSQ = LENGTHR*LENGTHR*TWO
          AXISISQ = LENGTHI*LENGTHI*TWO
          D = (SPAR(6)+SPAR(5))/TWO
          SIGMASQ = (AXISRSQ-AXISISQ)/ (ONE-D)**2
          GAMMA = ONE/ (ONE-D)

      END IF

*  2. k=gamma*Q1b
      IF (PRECONTYPE.EQ.0) THEN
          CALL CCOPY(LOCLEN,B,1,WRK(IK),1)
          CALL CSCAL(LOCLEN,CMPLX(GAMMA),WRK(IK),1)

      ELSE IF ((PRECONTYPE.EQ.1) .OR. (PRECONTYPE.EQ.3)) THEN
          CALL PRECONL(B,WRK(IK),IPAR)
          CALL CSCAL(LOCLEN,CMPLX(GAMMA),WRK(IK),1)
      END IF

*    xold=x
      CALL CCOPY(LOCLEN,X,1,WRK(IXOLD),1)

*  Loop
      STATUS = 0
      EXITNORM = -ONE
      STEPERR = -1
      DO 10 ITNO = 1,MAXIT

*  3. rho
          IF (ITNO.EQ.1) THEN
              RHO = ONE

          ELSE IF (ITNO.EQ.2) THEN
              RHO = ONE/ (ONE-SIGMASQ/TWO)

          ELSE
              RHO = ONE/ (ONE-RHO*SIGMASQ/4.0)
          END IF

*  4. w=(I-Q1AQ2)x
          IF (PRECONTYPE.EQ.0) THEN
              CALL MATVEC(X,WRK(IZ),IPAR)

          ELSE IF (PRECONTYPE.EQ.1) THEN
              CALL MATVEC(X,WRK(IW),IPAR)
              CALL PRECONL(WRK(IW),WRK(IZ),IPAR)

          ELSE IF (PRECONTYPE.EQ.2) THEN
              CALL PRECONR(X,WRK(IW),IPAR)
              CALL MATVEC(WRK(IW),WRK(IZ),IPAR)

          ELSE IF (PRECONTYPE.EQ.3) THEN
              CALL PRECONR(X,WRK(IZ),IPAR)
              CALL MATVEC(WRK(IZ),WRK(IW),IPAR)
              CALL PRECONL(WRK(IW),WRK(IZ),IPAR)
          END IF

          CALL CCOPY(LOCLEN,X,1,WRK(IW),1)
          CALL CAXPY(LOCLEN,-ZONE,WRK(IZ),1,WRK(IW),1)

*  5. x=rho*(gamma*((I-Q1A)x+Q1b)+(1-gamma)*x)+(1-rho)*xold
          DELTA = RHO*GAMMA
          CALL CSCAL(LOCLEN,CMPLX(ONE-RHO),WRK(IXOLD),1)
          CALL CAXPY(LOCLEN,CMPLX(RHO),WRK(IK),1,WRK(IXOLD),1)
          CALL CAXPY(LOCLEN,CMPLX(RHO-DELTA),X,1,WRK(IXOLD),1)
          CALL CAXPY(LOCLEN,CMPLX(DELTA),WRK(IW),1,WRK(IXOLD),1)
          CALL CSWAP(LOCLEN,WRK(IXOLD),1,X,1)

*  6. check stopping criterion
          CALL STOPCRIT(B,WRK(IZ),WRK(IR),X,WRK(IXOLD),WRK(IW),RHSSTOP,
     +                  CNVRTX,EXITNORM,STATUS,IPAR,MATVEC,MATVEC,
     +                  PRECONR,PCSUM,PSCNRM)

*  Call monitoring routine
          CALL PROGRESS(LOCLEN,ITNO,EXITNORM,X,WRK(IR),WRK(IR))

          IF (STATUS.EQ.0) THEN
              GO TO 9999
          END IF

   10 CONTINUE

      IF (ITNO.GT.MAXIT) THEN
          STATUS = -1
          ITNO = MAXIT
      END IF

 9999 CONTINUE

      IF ((PRECONTYPE.EQ.2) .OR. (PRECONTYPE.EQ.3)) THEN
          CALL CCOPY(LOCLEN,X,1,WRK(IZ),1)
          CALL PRECONR(WRK(IZ),X,IPAR)
      END IF

*  Set output parameters
      IPAR(11) = ITNO
      IPAR(12) = STATUS
      IPAR(13) = STEPERR
      SPAR(2) = EXITNORM

      RETURN

      END
      SUBROUTINE PIMCQMR(X,B,WRK,IPAR,SPAR,MATVEC,TMATVEC,PRECONL,
     +                   PRECONR,PCSUM,PSCNRM,PROGRESS)
      IMPLICIT NONE

*           PIM -- The Parallel Iterative Methods package
*           ---------------------------------------------
*
*                      Rudnei Dias da Cunha
*     National Supercomputing Centre and Mathematics Institute
*         Universidade Federal do Rio Grande do Sul, Brasil
*
*                          Tim Hopkins
*     Computing Laboratory, University of Kent at Canterbury, U.K.
*
* ----------------------------------------------------------------------
*
*     .. Parameters ..
      REAL ONE
      PARAMETER (ONE=1.0E0)
      COMPLEX CZERO
      PARAMETER (CZERO= (0.0E0,0.0E0))
      COMPLEX CONE
      PARAMETER (CONE= (1.0E0,0.0E0))
      INTEGER IPARSIZ
      PARAMETER (IPARSIZ=13)
      INTEGER SPARSIZ
      PARAMETER (SPARSIZ=2)
*     ..
*     .. Array Arguments ..
      COMPLEX B(*),WRK(*),X(*)
      REAL SPAR(SPARSIZ)
      INTEGER IPAR(IPARSIZ)
*     ..
*     .. Function Arguments ..
      REAL PSCNRM
      EXTERNAL PSCNRM
*     ..
*     .. Subroutine Arguments ..
      EXTERNAL MATVEC,PCSUM,PRECONL,PRECONR,PROGRESS,TMATVEC
*     ..
*     .. Local Scalars ..
      COMPLEX BETA,C,C0,DE,DELTA,DI,EPS,ETA,KSI,OMEGA,RHO,RHO0,THETA,
     +        THETA0
      REAL EPSILON,EXITNORM,OVERFLOW,RHSSTOP
      INTEGER BASISDIM,BLKSZ,CNVRTX,ID,IP,IQ,IR,ITNO,IV,IVTILDE,IW,
     +        IWTILDE,IXOLD,IY,IZ,LDA,LOCLEN,MAXIT,N,NPROCS,PRECONTYPE,
     +        PROCID,STATUS,STEPERR,STOPTYPE
*     ..
*     .. Local Arrays ..
      COMPLEX DOTS(2)
*     ..
*     .. External Functions ..
      COMPLEX CDOTC
      REAL SCSETRHSSTOP
      EXTERNAL CDOTC,SCSETRHSSTOP
*     ..
*     .. External Subroutines ..
      EXTERNAL CAXPY,CCOPY,CINIT,CSCAL,PIMSGETPAR,SMACHCONS,STOPCRIT
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC ABS,SQRT
*     ..
      CALL SMACHCONS('O',OVERFLOW)

      CALL PIMSGETPAR(IPAR,SPAR,LDA,N,BLKSZ,LOCLEN,BASISDIM,NPROCS,
     +                PROCID,PRECONTYPE,STOPTYPE,MAXIT,ITNO,STATUS,
     +                STEPERR,EPSILON,EXITNORM)

*  Check consistency of preconditioning and stop types
      IF (((PRECONTYPE.EQ.0).OR. (PRECONTYPE.EQ.2)) .AND.
     +    (STOPTYPE.EQ.6)) THEN
          ITNO = 0
          STATUS = -4
          STEPERR = 0
          GO TO 9999

      END IF

*  Does not need conversion Y=Q2X for residual
      CNVRTX = 1

*  Set indices for mapping local vectors into wrk
      IR = 1
      IV = IR + LOCLEN
      IW = IV + LOCLEN
      IP = IW + LOCLEN
      IQ = IP + LOCLEN
      ID = IQ + LOCLEN
      IVTILDE = ID + LOCLEN
      IWTILDE = IVTILDE + LOCLEN
      IXOLD = IWTILDE + LOCLEN
      IZ = IXOLD + LOCLEN
      IY = IZ + LOCLEN

*  Set RHS of stopping criteria
      RHSSTOP = SCSETRHSSTOP(B,WRK(IR),EPSILON,IPAR,PRECONL,PSCNRM)

*  1. r=Q1(b-AQ2x)
      IF (STOPTYPE.NE.6) THEN
          IF (PRECONTYPE.EQ.0) THEN
*     r=b-Ax
              CALL CCOPY(LOCLEN,B,1,WRK(IR),1)
              CALL MATVEC(X,WRK(IW),IPAR)
              CALL CAXPY(LOCLEN,-CONE,WRK(IW),1,WRK(IR),1)

          ELSE IF (PRECONTYPE.EQ.1) THEN
*     r=Q1(b-Ax)
              CALL CCOPY(LOCLEN,B,1,WRK(IZ),1)
              CALL MATVEC(X,WRK(IW),IPAR)
              CALL CAXPY(LOCLEN,-CONE,WRK(IW),1,WRK(IZ),1)
              CALL PRECONL(WRK(IZ),WRK(IR),IPAR)

          ELSE IF (PRECONTYPE.EQ.2) THEN
*     r=b-AQ2x
              CALL CCOPY(LOCLEN,B,1,WRK(IR),1)
              CALL PRECONR(X,WRK(IW),IPAR)
              CALL MATVEC(WRK(IW),WRK(IZ),IPAR)
              CALL CAXPY(LOCLEN,-CONE,WRK(IZ),1,WRK(IR),1)

          ELSE IF (PRECONTYPE.EQ.3) THEN
*     r=Q1(b-AQ2x)
              CALL CCOPY(LOCLEN,B,1,WRK(IP),1)
              CALL PRECONR(X,WRK(IW),IPAR)
              CALL MATVEC(WRK(IW),WRK(IZ),IPAR)
              CALL CAXPY(LOCLEN,-CONE,WRK(IZ),1,WRK(IP),1)
              CALL PRECONL(WRK(IP),WRK(IR),IPAR)
          END IF

      ELSE
*     r has been set to Q1b in the call to dsetrhsstop
          IF (PRECONTYPE.EQ.1) THEN
*     r=Q1(b-Ax)
              CALL MATVEC(X,WRK(IW),IPAR)
              CALL PRECONL(WRK(IW),WRK(IZ),IPAR)
              CALL CAXPY(LOCLEN,-CONE,WRK(IZ),1,WRK(IR),1)

          ELSE IF (PRECONTYPE.EQ.3) THEN
*     r=Q1(b-AQ2x)
              CALL PRECONR(X,WRK(IZ),IPAR)
              CALL MATVEC(WRK(IZ),WRK(IW),IPAR)
              CALL PRECONL(WRK(IW),WRK(IZ),IPAR)
              CALL CAXPY(LOCLEN,-CONE,WRK(IZ),1,WRK(IR),1)
          END IF

      END IF

*  2. rho=||r||_{2}
      RHO = PSCNRM(LOCLEN,WRK(IR))

*  3. v=r/rho
      IF (RHO.EQ.CZERO) THEN
          STATUS = -3
          STEPERR = 3
          GO TO 9999

      END IF

      CALL CCOPY(LOCLEN,WRK(IR),1,WRK(IV),1)
      CALL CSCAL(LOCLEN,CONE/RHO,WRK(IV),1)

*  4. w=-r/rho
      CALL CCOPY(LOCLEN,WRK(IR),1,WRK(IW),1)
      CALL CSCAL(LOCLEN,-CONE/RHO,WRK(IW),1)

*  5. p=q=d=0
      CALL CINIT(LOCLEN,CZERO,WRK(IP),1)
      CALL CINIT(LOCLEN,CZERO,WRK(IQ),1)
      CALL CINIT(LOCLEN,CZERO,WRK(ID),1)

*  6. c0=1, eps=1, ksi=1, theta0=0, eta=-1, omega=1
      C0 = CONE
      EPS = CONE
      KSI = CONE
      THETA0 = CZERO
      ETA = -CONE
      OMEGA = CONE

*  Loop
      STATUS = 0
      EXITNORM = -ONE
      STEPERR = -1
      DO 10 ITNO = 1,MAXIT

*  7. delta=w^{T}v
          DOTS(1) = CDOTC(LOCLEN,WRK(IW),1,WRK(IV),1)
          CALL PCSUM(1,DOTS)
          DELTA = DOTS(1)

*  8. if eps=0 then breakdown
          IF (EPS.EQ.CZERO) THEN
              STATUS = -3
              STEPERR = 8
              GO TO 9999

          END IF

*  9. if delta=0 then breakdown
          IF (DELTA.EQ.CZERO) THEN
              STATUS = -3
              STEPERR = 9
              GO TO 9999

          END IF

* 10. p=v-(ksi*delta/eps)*p
          DE = DELTA/EPS
          CALL CCOPY(LOCLEN,WRK(IP),1,WRK(IZ),1)
          CALL CCOPY(LOCLEN,WRK(IV),1,WRK(IP),1)
          CALL CAXPY(LOCLEN,-KSI*DE,WRK(IZ),1,WRK(IP),1)

* 11. q=w-(rho*delta/eps)*q
          CALL CCOPY(LOCLEN,WRK(IQ),1,WRK(IZ),1)
          CALL CCOPY(LOCLEN,WRK(IW),1,WRK(IQ),1)
          CALL CAXPY(LOCLEN,-RHO*DE,WRK(IZ),1,WRK(IQ),1)

* 12. vtilde=Q1AQ2p
          IF (PRECONTYPE.EQ.0) THEN
              CALL MATVEC(WRK(IP),WRK(IVTILDE),IPAR)

          ELSE IF (PRECONTYPE.EQ.1) THEN
              CALL MATVEC(WRK(IP),WRK(IZ),IPAR)
              CALL PRECONL(WRK(IZ),WRK(IVTILDE),IPAR)

          ELSE IF (PRECONTYPE.EQ.2) THEN
              CALL PRECONR(WRK(IP),WRK(IZ),IPAR)
              CALL MATVEC(WRK(IZ),WRK(IVTILDE),IPAR)

          ELSE IF (PRECONTYPE.EQ.3) THEN
              CALL PRECONR(WRK(IP),WRK(IVTILDE),IPAR)
              CALL MATVEC(WRK(IVTILDE),WRK(IZ),IPAR)
              CALL PRECONL(WRK(IZ),WRK(IVTILDE),IPAR)
          END IF

* 13. eps=q^{T}(Q1AQ2p)=q^{T}vtilde
          DOTS(1) = CDOTC(LOCLEN,WRK(IQ),1,WRK(IVTILDE),1)
          CALL PCSUM(1,DOTS)
          EPS = DOTS(1)

* 14. beta=eps/delta
          BETA = EPS/DELTA

* 15. vtilde=vtilde-beta*v
          CALL CAXPY(LOCLEN,-BETA,WRK(IV),1,WRK(IVTILDE),1)

* 16. wtilde=Q1A^{T}Q2q-beta*w
          IF (PRECONTYPE.EQ.0) THEN
              CALL TMATVEC(WRK(IQ),WRK(IWTILDE),IPAR)

          ELSE IF (PRECONTYPE.EQ.1) THEN
              CALL TMATVEC(WRK(IQ),WRK(IZ),IPAR)
              CALL PRECONL(WRK(IZ),WRK(IWTILDE),IPAR)

          ELSE IF (PRECONTYPE.EQ.2) THEN
              CALL PRECONR(WRK(IQ),WRK(IZ),IPAR)
              CALL TMATVEC(WRK(IZ),WRK(IWTILDE),IPAR)

          ELSE IF (PRECONTYPE.EQ.3) THEN
              CALL PRECONR(WRK(IQ),WRK(IWTILDE),IPAR)
              CALL TMATVEC(WRK(IWTILDE),WRK(IZ),IPAR)
              CALL PRECONL(WRK(IZ),WRK(IWTILDE),IPAR)
          END IF

          CALL CAXPY(LOCLEN,-BETA,WRK(IW),1,WRK(IWTILDE),1)

* 17. rho=norm(vtilde)
          RHO0 = RHO
          DOTS(1) = CDOTC(LOCLEN,WRK(IVTILDE),1,WRK(IVTILDE),1)

* 18. ksi=norm(wtilde)
          DOTS(2) = CDOTC(LOCLEN,WRK(IWTILDE),1,WRK(IWTILDE),1)

*  Accumulate simultaneously partial values
          CALL PCSUM(2,DOTS)
          IF (ABS(DOTS(1)).GE.OVERFLOW) THEN
              STATUS = -3
              STEPERR = 17
              GO TO 9999

          END IF
          RHO = SQRT(DOTS(1))

          IF (ABS(DOTS(2)).GE.OVERFLOW) THEN
              STATUS = -3
              STEPERR = 18
              GO TO 9999

          END IF
          KSI = SQRT(DOTS(2))

* 19. theta=(omega*rho)/(omega*c0*abs(beta))
          DI = OMEGA*C0*ABS(BETA)
          IF (DI.EQ.CZERO) THEN
              STATUS = -3
              STEPERR = 19