lixia2.f90

解稀疏对称方程组的ICCG法fortran源程序 有算例

	    MODULE PARA
		PARAMETER(NND=6,NNC=10)
		END MODULE

	  program iccgsolver
      use para
      integer nd (nnd),nc(nnc)
      real*8 AA(NNC),SS(NNC),T(NND),BB(NND)

      OPEN (1,FILE='INPUT1.DAT',STATUS='OLD')
      OPEN (2,FILE='OUTPUT.DAT',STATUS='UNKNOWN')

      READ (1,*)(AA(I),I=1,NNC)
      WRITE (2,10)(AA(I),I=1,NNC)

      READ (1,*)(T(I),I=1,NND)
      WRITE (2,20)(T(I),I=1,NND)

      READ (1,*)(ND(I),I=1,NND)
      WRITE (2,30)(ND(I),I=1,NND)

      READ (1,*)(NC(I),I=1,NNC)
      WRITE (2,40)(NC(I),I=1,NNC)

      DO 5 I=1,NNC
      SS (I)=AA(I)						  
5     CONTINUE
      CALL SILLT (NND,AA,ND,NC)
	  CALL ICCG (AA,SS,BB,T,ND,NC)

	  WRITE (2,50)(BB(I),I=1,NND)
10    FORMAT (5X,'******A******'/(1X,5E13.5))
20    FORMAT (/5X,'******T******'/(1X,5E13.5))
30    FORMAT (/5X,'******ND******'/(1X,10I6))
40    FORMAT (/5X,'******NC******'/(1X,10I6)/)
50    FORMAT (//5X,'******Equations solution ******'//(1X,5E13.5))
      STOP
	  END

	  SUBROUTINE SILLT (N,AA,ND,NC)
	  USE PARA
	  REAL*8 AA(NNC),CA
	  INTEGER ND(NND),NC(NNC)
	  DO 148 I=1 ,N
	    IF (I.EQ.1)MI=1
		IF (I.GT.1)MI=ND(I-1)+1
		NI=ND(I)
		M1=NI-1
		DO 147 K=MI,M1
		IF (MI.GT.M1) GO TO 147
		NS=NC(K)
		KKA=ND(NS)
		IF (ABS(AA(KKA)).LT.1.0E-20) THEN
		WRITE(*,*)' Denominator is Zero (分母为零，请检查原系数矩阵对角线元素)'
		WRITE(*,*)'Line number I=',I
		STOP 1111
		END IF
		CA=AA(K)/AA(KKA)
		K1=K+1
		DO 146 J=K1,NI
		NB=NC(J)
		N2=ND(NB-1)+1
		IF (NC(N2).GT.NS) GO TO 146
		M2=ND(NB)
		DO 145 L=N2,M2
		  IF (NC(L).NE.NS) GO TO 145
		  AA(J)=AA(J)-CA*AA(L)
		  GO TO 146
145       CONTINUE
146     CONTINUE
147     CONTINUE
        CB=ABS(AA(NI))
		IF (CB.GT.1.E-20) GO TO 148
		DO 152 J1=MI,M1
		  CB=CB+ABS(AA(J1))
152     CONTINUE
        AA(NI)=CB
148     CONTINUE
        RETURN
		END
		SUBROUTINE ICCG (AA,SS,BB,T,ND,NC)
		USE PARA
		REAL*8 AA(NNC),SS(NNC)
		REAL*8 P(NND),R(NND),V(NND)
		REAL*8 TF,TR,TN,ARF,BAT,T(NND),BB(NND)
		INTEGER ND(NND),NC(NNC)
		LLL=1
		LL=1
		DO 1 I=1,NND
		  BB(I)=T(I)
1       CONTINUE
        CALL SINA (NND,AA,BB,ND,NC)
		DO 5 I=1,NND
		  C=ABS(T(I))
		  EP0=EP0+C
5       CONTINUE
        WRITE(*,*)'EP0=',EP0
		EPS=EP0*1.E-5
10      CALL SMI (SS,V,BB,ND,NC)
        DO 100 I=1,NND
		  R(I)=T(I)-V(I)
		  P(I)=R(I)
100     CONTINUE
        CALL SINA(NND,AA,P,ND,NC)
		TF=0.
		DO 105 I=1,NND
		TF=TF+R(I)*P(I)
105     CONTINUE
        DO 50 J=1,NND
		  TR=0.
		  EP=0.
		  TN=0.
		CALL SMI (SS,V,P,ND,NC)
		DO 110 I=1,NND
		  TR=TR+P(I)*V(I)
110     CONTINUE
        ARF=TF/TR
		DO 120 I=1,NND
		  BB(I)=BB(I)+ARF*P(I)
		  R(I)=R(I)-ARF*V(I)
		  V(I)=R(I)
		  C=ABS(V(I))
		  EP=EP+C
120     CONTINUE
        IF (LLL.EQ.1.OR.MOD(LLL,1).EQ.0) THEN
		  WRITE(*,*)'LLL=',LLL,'EP=',EP
		  WRITE(2,*)'LLL=',LLL,'EP=',EP
		END IF
256     IF (EP.LT.EPS) GOTO 15
        CALL SINA(NND,AA,V,ND,NC)
		DO 130 I=1,NND
		  TN=TN+R(I)*V(I)
130     CONTINUE
        BAT=TN/TF
		TF=TN
		DO 140 I=1,NND
		  P(I)=BAT*P(I)+V(I)
140     CONTINUE
        LLL=LLL+1
50      CONTINUE
        LL=LL+1
		GO TO 10
15      WRITE(*,*)'LLL=',LLL,'EP=',EP
        END
		
		SUBROUTINE SINA (N,AA,F,ND,NC)
		USE PARA
		REAL*8 AA(NNC)
		REAL*8 F(NND),CA
		INTEGER ND(NND),NC(NNC)
		DO 130 I=1,N
		  IF (I.EQ.1) MI=1
		  IF (I.GT.1) MI=ND(I-1)+1
		  M1=ND(I)-1
		  DO 130 K=MI,M1
		    IF (MI.GT.M1) GO TO 130
			NS=NC(K)
			KKA=ND(NS)
			CA=AA(K)/AA(KKA)
			F(I)=F(I)-CA*F(NS)
130      CONTINUE
         DO 140 I2=1,N
		   I=N-I2+1
		   IF (I.EQ.1) MI=1
		   IF (I.GT.1) MI=ND(I-1)+1
		   NI=ND(I)
		   M1=NI-1
		   F(I)=F(I)/AA(NI)
		   DO 140 K=MI,M1
		     IF (MI.GT.M1) GO TO 140
			 NS=NC(K)
			 F(NS)=F(NS)-AA(K)*F(I)
140     CONTINUE
        RETURN
		END

		SUBROUTINE SMI(SS,Y,X,ND,NC)
		USE PARA
		REAL*8 SS(NNC)
		REAL*8 X(NND),Y(NND)
		INTEGER ND(NND),NC(NNC)
		DO 150 I=1,NND
		  IF (I.EQ.1) MI=1
		  IF (I.GT.1) MI=ND(I-1)+1
		  NI=ND(I)
		  M1=NI-1
		  Y(I)=X(I)*SS(NI)
		  DO 155 K=MI,M1
		    NS=NC(K)
			Y(I)=Y(I)+X(NS)*SS(K)
			Y(NS)=Y(NS)+X(I)*SS(K)
155       CONTINUE
150     CONTINUE
        RETURN
		END