iccg在解有限差分方程中的应用.f90

iccg解有限差分方程 不完全分解共轭梯度法解有限差分线性方程组 行业软件

     SUBROUTINE ICCG(A,N,N1,N2,M1,M2,B,X,D,R,P,Q,EPS,ITR,IER,S)
!***********************************************************************
!  INCOMPLETE CHOLESKY DECOMPOSITION CONJUGATED GRADIENT METHOD        *
!      FOR FINITE DIFFERENCE METHOD.                                   *
!  PARAMETERS                                                          *
!    (1) A: 2-DIM. ARRAY CONTAINING THE MATRIX                         *
!    (2) N: ORDER OF THE MATRIX (A)                                    *
!    (3) N1: UPPER ROW SIZE OF THE ARRAY (A)                           *
!    (4) N2: LOWER ROW SIZE OF THE ARRAY (A)                           *
!    (5) M1: NUMBER OF MESH POINTS ON X-AXIS                           *
!    (6) M2: NUMBER OF MESH POINTS ON Y-AXIS                           *
!    (7) B: 1-DIM. ARRAY CONTAINING RIGHT HAND VECTOROF THE EQUATION   *
!    (8) X: 1-DIM. ARRAY CONTAINING THE SOLUTION VECTOR                *
!    (9) D: 1-DIM. WORKING ARRAY                                       *
!   (10) R: 1-DIM. WORKING ARRAY                                       *
!   (11) P: 1-DIM. WORKING ARRAY                                       *
!   (12) Q: 1-DIM. WORKING ARRAY                                       *
!   (13) EPS: TOLERANCE FOR CONVERGENCE                                *
!   (14) ITR: MAXIMUM NUMBER OF REPITITION                             *
!   (15) IER: ERROR CODE                                               *
!   (16) S: CODE WHICH SPECIFIES A METHOD TO BE USED                   *

!***********************************************************************
       IMPLICIT real*8(A-H,O-Z)
       DIMENSION A(N2:N1,4),B(N),X(N2:N+M2),D(N2:N+M2),R(N2:N+M2),	 
     &   P(N2:N+M2),Q(N2:N+M2)

      IF ((M1 .GE. N) .OR. (M2 .GE. N) .OR. (M1 .GE. N) .OR.		 
     &    (M1 .LE. 1) .OR. (M2 .LE. 1) .OR. (S .LT. 0.0)) THEN
       WRITE(*,*) ' (SUBR. ICCG) INVALID ARGUMENT.'
       IER = 2
       RETURN
      ENDIF

      TH = 1.0D0
      IF (S .GT. 0.0 .AND. S .LT. 1.0) THEN
       TH = S
       S = 1.0D0
       ENDIF
      DO 700 I=1-M2,0
       X(I) = 0.0D0
       P(I) = 0.0D0
       Q(I) = 0.0D0
       D(I) = 0.0D0
       X(I+N+M2) = 0.0D0
       P(I+N+M2) = 0.0D0
       Q(I+N+M2) = 0.0D0
  700 CONTINUE
      DO 710 I=1,N
  710  D(I) = 0.0D0
!  INCOMPLETE CHOLESKY DECOMPOSITION
      IF (S .NE. 0.0) THEN
       DO 720 I=1,N
        W=S*A(I,4)-A(I,3)*(A(I,3)+(A(I-1+M2,1)+A(I-1+M1,2))*TH)*D(I-1)	 
     &          -A(I,2)*(A(I,2)+(A(I-M1+1,3)+A(I-M1+M2,1))*TH)*D(I-M1)	  
     &          -A(I,1)*(A(I,1)+(A(I-M2+M1,2)+A(I-M2+1,3))*TH)*D(I-M2)
  720   D(I) = 1.0D0 / W
      ELSE
       DO 730 I=1,N
        W=A(I,4)-D(I-1)*A(I,3)**2-D(I-M1)*A(I,2)**2-D(I-M2)*A(I,1)**2
  730   D(I) = 1.0D0 / W
      ENDIF

      DO 740 I=1,N
       Q(I) = A(I,1) * X(I-M2) + A(I,2) * X(I-M1) + A(I,3) * X(I-1)		   
     &        + A(I,4) * X(I) + A(I+1,3) * X(I+1) + A(I+M1,2) * X(I+M1)	   
     &        + A(I+M2,1) * X(I+M2)
  740 CONTINUE
      DO 750 I=1,N
  750  R(I) = B(I) - Q(I)

      DO 760 I=1,N
  760  P(I) = D(I) * (R(I) - A(I,3) * P(I-1) - A(I,2) * P(I-M1)			
     &                - A(I,1) * P(I-M2))
      DO 770 I=N,1,-1
  770  P(I) =  D(I) * (P(I) -A(I+1,3) * P(I+1) - A(I+M1,2) * P(I+M1)	 
     &                - A(I+M2,1) * P(I+M2))
      C1 = 0.0D0
      DO 780 I=1,N
  780  C1 = C1 + R(I) * P(I)
!  ITERATION
      DO 880 K=1,ITR
!	WRITE(6,*) K
       DO 810 I=1,N
        Q(I) = A(I,1) * P(I-M2) + A(I,2) * P(I-M1) + A(I,3) * P(I-1)	
     &         + A(I,4) * P(I) + A(I+1,3) * P(I+1)					
     &         + A(I+M1,2) * P(I+M1) + A(I+M2,1) * P(I+M2)
  810  CONTINUE
       C2 = 0.0D0
       DO 820 I=1,N
  820   C2 = C2 + P(I) * Q(I)
       ALPHA = C1 / C2
       X1 = 0.0D0
       X2 = 0.0D0
       DO 830 I=1,N
        Y = X(I)
        X(I) = X(I) + ALPHA * P(I)
        R(I) = R(I) - ALPHA * Q(I)
        X1 = X1 + Y * Y
        X2 = X2 + (X(I) - Y)**2
  830  CONTINUE
!	PAUSE 400
!	WRITE(6,*) X2
!	WRITE(6,*) X1
!	PAUSE 500
       IF (X1 .NE. 0.0) THEN
	
        RES = DSQRT(X2 / X1)
!	WRITE(6,*) RES
        IF (RES .LE. EPS) THEN
         ITRR = K
         EPS = RES
         IER = 0
         IF (TH .NE. 1.0D0) S = TH
         GOTO 900
        ENDIF
       ENDIF
	
       DO 840 I=1,N
  840   Q(I) = D(I) * (R(I) - A(I,3) * Q(I-1) - A(I,2) * Q(I-M1)		
     &                  - A(I,1) * Q(I-M2))
       DO 850 I=N,1,-1
  850   Q(I) = D(I) * (Q(I) - A(I+1,3) * Q(I+1) - A(I+M1,2) * Q(I+M1)	 
     &                       - A(I+M2,1) * Q(I+M2))
       C3 = 0.0D0
       DO 860 I=1,N
  860   C3 = C3 + R(I) * Q(I)
       BETA = C3 / C1
       C1 = C3
       DO 870 I=1,N
  870   P(I) = Q(I) + BETA * P(I)
  880 CONTINUE
      IER = 1
      WRITE(*,*) ' (SUBR. ICCG) NO CONVERGENCE. '
!      EPS = RES
      IF (TH .NE. 1.0D0) S = TH
      RETURN
  900 ITR=ITRR
	WRITE(6,*) 'ITERATION NUMBER'
	WRITE(6,*) K
      RETURN
      END