pim.f

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

      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,TMATVEC,
     +                  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=Q1A^{T}AQ2r
          IF (PRECONTYPE.EQ.0) THEN
              CALL MATVEC(WRK(IR),WRK(IZ),IPAR)
              CALL TMATVEC(WRK(IZ),WRK(IS),IPAR)

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

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

          ELSE IF (PRECONTYPE.EQ.3) THEN
              CALL PRECONR(WRK(IR),WRK(IZ),IPAR)
              CALL MATVEC(WRK(IZ),WRK(IS),IPAR)
              CALL TMATVEC(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 PIMCCGS(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,RHO,RHO0,XI
      REAL EPSILON,EXITNORM,RHSSTOP
      INTEGER BASISDIM,BLKSZ,CNVRTX,IP,IQ,IR,IRTILDE,IS,IT,ITNO,IU,IW,
     +        IXOLD,IZ,LDA,LOCLEN,MAXIT,N,NPROCS,PRECONTYPE,PROCID,
     +        STATUS,STEPERR,STOPTYPE
*     ..
*     .. Local Arrays ..
      COMPLEX DOTS(1)
*     ..
*     .. 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
      IRTILDE = IR + LOCLEN
      IP = IRTILDE + LOCLEN
      IQ = IP + LOCLEN
      IS = IQ + LOCLEN
      IT = IS + LOCLEN
      IU = IT + LOCLEN
      IW = IU + LOCLEN
      IZ = IW + LOCLEN
      IXOLD = 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 Qb in the call to ssetrhsstop
          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=s=rtilde=r
      CALL CCOPY(LOCLEN,WRK(IR),1,WRK(IRTILDE),1)
      CALL CCOPY(LOCLEN,WRK(IR),1,WRK(IP),1)
      CALL CCOPY(LOCLEN,WRK(IR),1,WRK(IS),1)

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

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

*  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(rtilde,w)
          DOTS(1) = CDOTC(LOCLEN,WRK(IRTILDE),1,WRK(IW),1)
          CALL PCSUM(1,DOTS)
          XI = DOTS(1)

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

          END IF

          ALPHA = RHO/XI

*  7. t=s-alpha*w
          CALL CCOPY(LOCLEN,WRK(IS),1,WRK(IT),1)
          CALL CAXPY(LOCLEN,-ALPHA,WRK(IW),1,WRK(IT),1)

*  8. w=s+t
          CALL CCOPY(LOCLEN,WRK(IS),1,WRK(IW),1)
          CALL CAXPY(LOCLEN,CONE,WRK(IT),1,WRK(IW),1)

*  9. x=x+alpha*w
          CALL CCOPY(LOCLEN,X,1,WRK(IXOLD),1)
          CALL CAXPY(LOCLEN,ALPHA,WRK(IW),1,X,1)

* 10. r=r-alpha*Q1AQ2w
          IF (PRECONTYPE.EQ.0) THEN
              CALL MATVEC(WRK(IW),WRK(IU),IPAR)

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

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

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

          CALL CAXPY(LOCLEN,-ALPHA,WRK(IU),1,WRK(IR),1)

* 11. check stopping criterion
          CALL STOPCRIT(B,WRK(IR),WRK(IZ),X,WRK(IXOLD),WRK(IU),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
* 12. rho=dot(rtilde0,r)
          RHO0 = RHO
          DOTS(1) = CDOTC(LOCLEN,WRK(IRTILDE),1,WRK(IR),1)
          CALL PCSUM(1,DOTS)
          RHO = DOTS(1)

* 13. beta=rho/rho0
          IF (RHO0.EQ.CZERO) THEN
              STATUS = -3
              STEPERR = 13
              GO TO 9999

          END IF

          BETA = RHO/RHO0

* 14. s=r+beta*t
          CALL CCOPY(LOCLEN,WRK(IR),1,WRK(IS),1)
          CALL CAXPY(LOCLEN,BETA,WRK(IT),1,WRK(IS),1)

* 15. w=t+beta*p
          CALL CCOPY(LOCLEN,WRK(IT),1,WRK(IW),1)
          CALL CAXPY(LOCLEN,BETA,WRK(IP),1,WRK(IW),1)

* 16. p=s+beta*w
          CALL CCOPY(LOCLEN,WRK(IS),1,WRK(IP),1)
          CALL CAXPY(LOCLEN,BETA,WRK(IW),1,WRK(IP),1)

   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 PIMCCHEBYSHEV(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 ZERO
      PARAMETER (ZERO=0.0E0)
      COMPLEX ZONE
      PARAMETER (ZONE= (1.0E0,0.0E0))
      REAL ONE
      PARAMETER (ONE=1.0E0)
      REAL TWO
      PARAMETER (TWO=2.0)
      INTEGER IPARSIZ
      PARAMETER (IPARSIZ=13)
      INTEGER SPARSIZ
      PARAMETER (SPARSIZ=6)
*     ..
*     .. 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 ..
      REAL AXISISQ,AXISRSQ,D,DELTA,EPSILON,EXITNORM,GAMMA,LENGTHI,
     +     LENGTHR,RHO,RHSSTOP,SIGMA,SIGMASQ
      INTEGER BASISDIM,BLKSZ,CNVRTX,IK,IR,ITNO,IW,IXOLD,IZ,LDA,LOCLEN,
     +        MAXIT,N,NPROCS,PRECONTYPE,PROCID,STATUS,STEPERR,STOPTYPE
*     ..
*     .. External Functions ..
      REAL SCSETRHSSTOP
      EXTERNAL SCSETRHSSTOP
*     ..
*     .. External Subroutines ..
      EXTERNAL CAXPY,CCOPY,CSCAL,CSWAP,PIMSGETPAR,STOPCRIT
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC CMPLX,MAX
*     ..
      CALL PIMSGETPAR(IPAR,SPAR,LDA,N,BLKSZ,LOCLEN,BASISDIM,NPROCS,
     +                PROCID,PRECONTYPE,STOPTYPE,MAXIT,ITNO,STATUS,
     +                STEPERR,EPSILON,EXITNORM)