pim.f

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

*  4. rho=alpha=omega=1
      RHO = CONE
      ALPHA = CONE
      OMEGA = CONE

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

*  5. rho=dot(rtilde,r)
          RHO0 = RHO
          DOTS(1) = CDOTC(LOCLEN,WRK(IRTILDE),1,WRK(IR),1)
          CALL PCSUM(1,DOTS)
          RHO = DOTS(1)

*  6. beta=rho*alpha/(rho0*omega)
          KAPPA = RHO0*OMEGA
          IF (KAPPA.EQ.CZERO) THEN
              STATUS = -3
              STEPERR = 6
              GO TO 9999

          END IF

          BETA = RHO*ALPHA/KAPPA

*  7. p=r+beta*(p-omega*v)
          CALL CAXPY(LOCLEN,-OMEGA,WRK(IV),1,WRK(IP),1)
          CALL CCOPY(LOCLEN,WRK(IP),1,WRK(IW),1)
          CALL CCOPY(LOCLEN,WRK(IR),1,WRK(IP),1)
          CALL CAXPY(LOCLEN,BETA,WRK(IW),1,WRK(IP),1)

*  8. v=Q1AQ2p
          IF (PRECONTYPE.EQ.0) THEN
              CALL MATVEC(WRK(IP),WRK(IV),IPAR)

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

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

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

*  9. xi=dot(rtilde,v)
          DOTS(1) = CDOTC(LOCLEN,WRK(IRTILDE),1,WRK(IV),1)
          CALL PCSUM(1,DOTS)
          XI = DOTS(1)

* 10. alpha=rho/xi
          IF (XI.EQ.CZERO) THEN
              STATUS = -3
              STEPERR = 10
              GO TO 9999

          END IF

          ALPHA = RHO/XI

* 11. s=r-alpha*v
          CALL CCOPY(LOCLEN,WRK(IR),1,WRK(IS),1)
          CALL CAXPY(LOCLEN,-ALPHA,WRK(IV),1,WRK(IS),1)

* 12. if ||s||<breaktol then soft-breakdown has occurred
          KAPPA = PSCNRM(LOCLEN,WRK(IS))
          IF (ABS(KAPPA).LT.MACHEPS) THEN
              STATUS = -2
              STEPERR = 12
              CALL CAXPY(LOCLEN,ALPHA,WRK(IP),1,X,1)
              GO TO 9999

          END IF

* 13. t=Q1AQ2s
          IF (PRECONTYPE.EQ.0) THEN
              CALL MATVEC(WRK(IS),WRK(IT),IPAR)

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

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

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

* 14. omega=dot(t,s)/dot(t,t)
          DOTS(1) = CDOTC(LOCLEN,WRK(IT),1,WRK(IT),1)
          DOTS(2) = CDOTC(LOCLEN,WRK(IT),1,WRK(IS),1)

*  Accumulate simultaneously partial values
          CALL PCSUM(2,DOTS)

          IF (DOTS(1).EQ.CZERO) THEN
              STATUS = -3
              STEPERR = 14
              GO TO 9999

          END IF

          OMEGA = DOTS(2)/DOTS(1)

* 15. x=x+alpha*p+omega*s
          CALL CCOPY(LOCLEN,X,1,WRK(IXOLD),1)
          CALL CAXPY(LOCLEN,ALPHA,WRK(IP),1,X,1)
          CALL CAXPY(LOCLEN,OMEGA,WRK(IS),1,X,1)

* 16. r=s-omega*t
          CALL CCOPY(LOCLEN,WRK(IS),1,WRK(IR),1)
          CALL CAXPY(LOCLEN,-OMEGA,WRK(IT),1,WRK(IR),1)

* 17. check stopping criterion
          CALL STOPCRIT(B,WRK(IR),WRK(IZ),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(IZ))

          IF (STATUS.EQ.0) THEN
              GO TO 9999
          END IF
   20 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 PIMCCG(X,B,WRK,IPAR,SPAR,MATVEC,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
*     ..
*     .. Local Scalars ..
      COMPLEX ALPHA,BETA,DELTA,RDOTR,RDOTR0,XI
      REAL EPSILON,EXITNORM,RHSSTOP
      INTEGER BASISDIM,BLKSZ,CNVRTX,IP,IR,IS,ITNO,IW,IXOLD,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,PIMSGETPAR,STOPCRIT
*     ..

      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 = 0

*  Set indices for mapping local vectors into wrk
      IR = 1
      IP = IR + LOCLEN
      IW = IP + LOCLEN
      IZ = IW + LOCLEN
      IS = IZ + LOCLEN
      IXOLD = IS + 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 Qb 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. p=r
      CALL CCOPY(LOCLEN,WRK(IR),1,WRK(IP),1)

*  3. rdotr=dot(r,r)
      DOTS(1) = CDOTC(LOCLEN,WRK(IR),1,WRK(IR),1)
      CALL PCSUM(1,DOTS)
      RDOTR = DOTS(1)

*  4. w=Q1AQ2p
      IF (PRECONTYPE.EQ.0) THEN
          CALL MATVEC(WRK(IP),WRK(IW),IPAR)

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

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

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

*  5. xi=dot(p,w)
      DOTS(1) = CDOTC(LOCLEN,WRK(IP),1,WRK(IW),1)
      CALL PCSUM(1,DOTS)
      XI = DOTS(1)

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

*  6. alpha=rdotr/xi
          IF (XI.EQ.CZERO) THEN
              STATUS = -3
              STEPERR = 6
              GO TO 9999

          END IF

          ALPHA = RDOTR/XI

*  7. x=x+alpha*p
          CALL CCOPY(LOCLEN,X,1,WRK(IXOLD),1)
          CALL CAXPY(LOCLEN,ALPHA,WRK(IP),1,X,1)

*  8. r=r-alpha*w
          CALL CAXPY(LOCLEN,-ALPHA,WRK(IW),1,WRK(IR),1)

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

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

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

* 10. s=Q1AQ2r
          IF (PRECONTYPE.EQ.0) THEN
              CALL MATVEC(WRK(IR),WRK(IS),IPAR)

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

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

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

* 11. rdotr=dot(r,r)
          RDOTR0 = RDOTR
          DOTS(1) = CDOTC(LOCLEN,WRK(IR),1,WRK(IR),1)

* 12. delta=dot(r,s)
          DOTS(2) = CDOTC(LOCLEN,WRK(IR),1,WRK(IS),1)

*  Accumulate simultaneously partial values
          CALL PCSUM(2,DOTS)
          RDOTR = DOTS(1)
          DELTA = DOTS(2)

* 13. beta=rdotr/rdotr0
          IF (RDOTR0.EQ.CZERO) THEN
              STATUS = -3
              STEPERR = 13
              GO TO 9999

          END IF

          BETA = RDOTR/RDOTR0

* 14. p=r+beta*p
          CALL CCOPY(LOCLEN,WRK(IP),1,WRK(IZ),1)
          CALL CCOPY(LOCLEN,WRK(IR),1,WRK(IP),1)
          CALL CAXPY(LOCLEN,BETA,WRK(IZ),1,WRK(IP),1)

* 15. w=s+beta*w
          CALL CCOPY(LOCLEN,WRK(IW),1,WRK(IZ),1)
          CALL CCOPY(LOCLEN,WRK(IS),1,WRK(IW),1)
          CALL CAXPY(LOCLEN,BETA,WRK(IZ),1,WRK(IW),1)

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

   20 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 PIMCCGNE(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