cgm.for

本fortran程序用共轭梯度法求解实系数的的矩阵方程

      
	
	
	!共轭梯度法度法解稀疏方程子程序
      SUBROUTINE  CGM(A,B,N2,ASUB,ATSUB,X,RSQ)
      PARAMETER (NMAX=1000,EPS=1.E-6)
      DOUBLE PRECISION  A,B,X,G,H,XI,XJ
      DIMENSION A(N2,N2),B(N2),X(N2),G(NMAX),H(NMAX),XI(NMAX),XJ(NMAX)
      EPS2=EPS2*N2
      IRST=0
1     IRST=IRST+1
      CALL  ASUB(A,X,XI,N2)
      RP=0.0
      BSQ=0.0
      DO  11 J=1,N2
        BSQ=BSQ+B(J)**2
	    XI(J)=XI(J)-B(J)
	    RP=RP+XI(J)**2
11    CONTINUE
      CALL  ATSUB(A,XI,G,N2)
       DO  12  J=1,N2
	   G(J)=-G(J)
	   H(J)=G(J)
12     CONTINUE
       DO  19  ITER=1,100*N2
	    CALL ASUB(A,H,XI,N2)
		ANUM=0.0
		ADEN=0.0
		DO  13 J=1,N2
		  ANUM=ANUM+G(J)*H(J)
		  ADEN=ADEN+XI(J)**2
13      CONTINUE
        IF(ADEN.EQ.0.0) PAUSE 'VERY SINGULAR MATRIX'
		ANUM=ANUM/ADEN
		DO  14 J=1,N2
		   XI(J)=X(J)
		   X(J)=X(J)+ANUM*H(J)
14      CONTINUE
        CALL ASUB(A,X,XJ,N2)
		 RSQ=0.0
		DO  15  J=1,N2
		  XJ(J)=XJ(J)-B(J)
		  RSQ=RSQ+XJ(J)**2
15      CONTINUE
        IF((RSQ.EQ.RP).OR.(RSQ.LE.BSQ*EPS2))  RETURN
		IF(RSQ.GT.RP) THEN
		  DO  16 J=1,N2
		    X(J)=XI(J)
16        CONTINUE
       	  IF(IRST.GE.3)  RETURN
		  GOTO  1
		ENDIF
		RP=RSQ
		CALL ATSUB(A,XJ,XI,N2)
		GG=0.0
		DGG=0.0
		DO 17 J=1,N2
		   GG=GG+G(J)**2
		   DGG=DGG+(XI(J)+G(J))*XI(J)
17		CONTINUE
		IF(GG.EQ.0.) RETURN
		GAM=DGG/GG
		DO  18 J=1,N2
		  G(J)=-XI(J)
		  H(J)=G(J)+GAM*H(J)
18      CONTINUE
19    CONTINUE
		PAUSE 'TOO MANY ITERATIONS'
		RETURN
       !计算一个矩阵与一个向量的乘积子程序
       SUBROUTINE    ASUB(A,XIN,XOUT,N2)
       DOUBLE  PRECISION  A,XIN,XOUT
	 DIMENSION  A(N2,N2),XIN(N2),XOUT(N2)
	   DO  10  I=1,N2
	     SUM=0.D0
		 DO  20  J=1,N2
		   SUM=SUM+A(I,J)*XIN(J)
20       CONTINUE
		 XOUT(I)=SUM
10     CONTINUE
       END
    !计算一个矩阵的转置矩阵与一个向量的乘积
       SUBROUTINE    ATSUB(A,XIN,XOUT,N2)
	  DOUBLE  PRECISION  A,XIN,XOUT
	  DIMENSION  A(N2,N2),XIN(N2),XOUT(N2)
	     DO  10  I=1,N2
	       SUM=0.D0
		   DO  20  J=1,N2
		     SUM=SUM+A(J,I)*XIN(J)
20         CONTINUE
		   XOUT(I)=SUM
10       CONTINUE
          END		END