pim.f

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

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

*  Needs conversion Y=Q2X for residual
      CNVRTX = 1

*  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-AA^{T}Q2x)
      IF (STOPTYPE.NE.6) THEN
          IF (PRECONTYPE.EQ.0) THEN
*     r=b-Ax
              CALL CCOPY(LOCLEN,B,1,WRK(IR),1)
              CALL TMATVEC(X,WRK(IP),IPAR)
              CALL MATVEC(WRK(IP),WRK(IW),IPAR)
              CALL CAXPY(LOCLEN,-CONE,WRK(IW),1,WRK(IR),1)

          ELSE IF (PRECONTYPE.EQ.1) THEN
*     r=Q1(b-AA^{T})
              CALL CCOPY(LOCLEN,B,1,WRK(IZ),1)
              CALL TMATVEC(X,WRK(IP),IPAR)
              CALL MATVEC(WRK(IP),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-AA^{T}Q2x
              CALL CCOPY(LOCLEN,B,1,WRK(IR),1)
              CALL PRECONR(X,WRK(IP),IPAR)
              CALL TMATVEC(WRK(IP),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-AA^{T}Q2x)
              CALL CCOPY(LOCLEN,B,1,WRK(IP),1)
              CALL PRECONR(X,WRK(IW),IPAR)
              CALL TMATVEC(WRK(IW),WRK(IZ),IPAR)
              CALL MATVEC(WRK(IZ),WRK(IW),IPAR)
              CALL CAXPY(LOCLEN,-CONE,WRK(IW),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-AA^{T}x)
              CALL TMATVEC(X,WRK(IW),IPAR)
              CALL MATVEC(WRK(IW),WRK(IP),IPAR)
              CALL PRECONL(WRK(IP),WRK(IZ),IPAR)
              CALL CAXPY(LOCLEN,-CONE,WRK(IZ),1,WRK(IR),1)

          ELSE IF (PRECONTYPE.EQ.3) THEN
*     r=Q1(b-AA^{T}Q2x)
              CALL PRECONR(X,WRK(IZ),IPAR)
              CALL TMATVEC(WRK(IZ),WRK(IW),IPAR)
              CALL MATVEC(WRK(IW),WRK(IZ),IPAR)
              CALL PRECONL(WRK(IZ),WRK(IW),IPAR)
              CALL CAXPY(LOCLEN,-CONE,WRK(IW),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=Q1AA^{T}Q2p
      IF (PRECONTYPE.EQ.0) THEN
          CALL TMATVEC(WRK(IP),WRK(IZ),IPAR)
          CALL MATVEC(WRK(IZ),WRK(IW),IPAR)

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

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

      ELSE IF (PRECONTYPE.EQ.3) THEN
          CALL PRECONR(WRK(IP),WRK(IZ),IPAR)
          CALL TMATVEC(WRK(IZ),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,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=Q1AA^{T}Q2r
          IF (PRECONTYPE.EQ.0) THEN
              CALL TMATVEC(WRK(IR),WRK(IZ),IPAR)
              CALL MATVEC(WRK(IZ),WRK(IS),IPAR)

          ELSE IF (PRECONTYPE.EQ.1) THEN
              CALL TMATVEC(WRK(IR),WRK(IS),IPAR)
              CALL MATVEC(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 TMATVEC(WRK(IS),WRK(IZ),IPAR)
              CALL MATVEC(WRK(IZ),WRK(IS),IPAR)

          ELSE IF (PRECONTYPE.EQ.3) THEN
              CALL PRECONR(WRK(IR),WRK(IZ),IPAR)
              CALL TMATVEC(WRK(IZ),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 = 10
              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),WRK(IXOLD),IPAR)
          CALL TMATVEC(WRK(IXOLD),X,IPAR)
      ELSE
          CALL CCOPY(LOCLEN,X,1,WRK(IZ),1)
          CALL TMATVEC(WRK(IZ),X,IPAR)
      END IF

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

      RETURN

      END
      SUBROUTINE PIMCCGNR(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 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
*  On (P)CGNR, the rhs vector is actually A^{T}b
      CALL TMATVEC(B,WRK(IZ),IPAR)
      RHSSTOP = SCSETRHSSTOP(WRK(IZ),WRK(IR),EPSILON,IPAR,PRECONL,
     +          PSCNRM)

*  1. r=Q1(A^{T}b-A^{T}AQ2x)
      IF (STOPTYPE.NE.6) THEN
          IF (PRECONTYPE.EQ.0) THEN
*     r=b-Ax
              CALL MATVEC(X,WRK(IP),IPAR)
              CALL TMATVEC(WRK(IP),WRK(IW),IPAR)
              CALL CAXPY(LOCLEN,-CONE,WRK(IW),1,WRK(IZ),1)
              CALL CCOPY(LOCLEN,WRK(IZ),1,WRK(IR),1)

          ELSE IF (PRECONTYPE.EQ.1) THEN
*     r=Q1(A^{T}b-A^{T}Ax)
              CALL MATVEC(X,WRK(IP),IPAR)
              CALL TMATVEC(WRK(IP),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=A^{T}b-A^{T}AQ2x
              CALL PRECONR(X,WRK(IW),IPAR)
              CALL MATVEC(WRK(IW),WRK(IP),IPAR)
              CALL TMATVEC(WRK(IP),WRK(IW),IPAR)
              CALL CAXPY(LOCLEN,-CONE,WRK(IW),1,WRK(IZ),1)
              CALL CCOPY(LOCLEN,WRK(IZ),1,WRK(IR),1)

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

      ELSE
*     r has been set to Q(A^{T}b) in the call to dsetrhsstop
          IF (PRECONTYPE.EQ.1) THEN
*     r=Q1(A^{T}b-A^{T}Ax)
              CALL MATVEC(X,WRK(IW),IPAR)
              CALL TMATVEC(WRK(IW),WRK(IP),IPAR)
              CALL PRECONL(WRK(IP),WRK(IW),IPAR)
              CALL CAXPY(LOCLEN,-CONE,WRK(IW),1,WRK(IR),1)

          ELSE IF (PRECONTYPE.EQ.3) THEN
*     r=Q1(A^{T}b-A^{T}AQ2x)
              CALL PRECONR(X,WRK(IW),IPAR)
              CALL MATVEC(WRK(IW),WRK(IP),IPAR)
              CALL TMATVEC(WRK(IP),WRK(IW),IPAR)
              CALL PRECONL(WRK(IW),WRK(IP),IPAR)
              CALL CAXPY(LOCLEN,-CONE,WRK(IP),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=Q1A^{T}AQ2p
      IF (PRECONTYPE.EQ.0) THEN
          CALL MATVEC(WRK(IP),WRK(IZ),IPAR)
          CALL TMATVEC(WRK(IZ),WRK(IW),IPAR)

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

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

      ELSE IF (PRECONTYPE.EQ.3) THEN
          CALL PRECONR(WRK(IP),WRK(IZ),IPAR)
          CALL MATVEC(WRK(IZ),WRK(IW),IPAR)
          CALL TMATVEC(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