juliang.f

MOM 巨量法 初学巨量法编程很好的参考。

      IMPLICIT NONE 
	INTEGER N, JS, IS 
      PARAMETER (n=6)
      DIMENSION L(n,n), g(n), IS(N), JS(N), A(N)
	REAL A, F, H 	
      DOUBLE PRECISION l, G
	CALL getmat1(l,n)
	CALL brinv(l,n,h,is,js)	
	CALL getmat2(g,n)  	
 	CALL mul(l,g,n,a)	
	CALL getnum(a,n,f)
	END
  
      subroutine getmat1(l,m)
	dimension l(m,m)
	double precision l 	
 	real i,j
	  do 10 j=1,m
	   do 20 i=1,m
	    t=i*(i+1)*(j/(m+1))**(i-1) 
	    l(i,j)=t
20       continue  
10      continue					      
      end	 
		
	subroutine brinv(a,m,h,is,js)
	dimension a(m,m),is(m),js(m)	
      double precision a,t,d
	h=1.0
	do 100 k=1,m
		d=0.0
		do 10 i=k,m
		do 10 j=k,m
	         if (abs(a(i,j)).gt.d) then
			  d=abs(a(i,j))
			  is(k)=i
			  js(k)=j
			 endif 
10           continue
        if(d+1.0.eq.1.0) then
	 h=0.0
	  write(*,20)
	   return
	  endif
20      format(1x,'err**not inv')
           do 30 j=1,m
	       t=a(k,j)
	       a(k,j)=a(is(k),j)
	       a(is(k),j)=t
30         continue
          do 40 i=1,m
	       t=a(i,k)
	       a(i,k)=a(i,js(k))
	       a(i,js(k))=t
40        continue
        a(k,k)=1/a(k,K)
	   do 50 j=1,m
	     if(j.ne.k) then
	     a(k,j)=a(k,j)*a(k,k)
	     endif
50       continue
          do 70 i=1,m
	      if	(i.ne.k) then
	     do 60 j=1,m
	       if(j.ne.k) then
	        a(i,j)=a(i,j)-a(i,k)*a(k,j)
	       endif
60         continue
            endif
70        continue
        do 80 i=1,m
	     if(i.ne.k) then
	      a(i,k)=-a(i,k)*a(k,k)
	     endif
80      continue
100   continue
        do 130 k=m,1,-1
	    do 110 j=1,m
	      t=a(k,j)
	      a(k,j)=a(js(k),j)
	      a(js(k),j)=t
110       continue
         do 120 i=1,m
	     t=a(i,k)
	     a(i,k)=a(i,is(k))
	     a(i,is(k))=t
120      continue
130     continue
        return
	end			
			
	subroutine getmat2(g,m)
	dimension g(m)
	real i
      double precision t,g
        do 20 i=1,m
	     t=1+4*(i/(m+1))**2 
	     g(i)=t	
20      continue 
      end
									      
      subroutine mul(l,g,m,a)
	dimension g(m),l(m,m),a(m)
      double precision g,l
	real a		 	
	 do 20 i=1,m
	    t=0.0
	   do 30 j=1,m
		 t=t+l(i,j)*g(j)
30       continue
         a(i)=t 	   			
20     continue	
    	end	     
     	
     	subroutine getnum(a,m,f)
	dimension a(m),t(m)	   
	real x,f,a 
	open (6,file='tab6.dat')  
       do 20 x=0.0,1.0,0.01
            f=0.0
	    do 10 i=1,m 
	      t(i)=x-x**(i+1)
	      f=f+a(i)*t(i)	
10        continue 
         write(6,100) x,f 
20     continue
        close(6)
100    format(f12.5,f12.5)
       end