module inv_mat
contains
subroutine inv(A,invA,N)
implicit real*8(a-z)
integer::n
integer::i
real*8::A(n,n),invA(n,n),E(n,n)
E=0
do i=1,n
   E(i,i)=1
end do
call mateq(A,B,X,N,N)
end subroutine inv

subroutine mateq(A,B,X,N,M)
implicit real*8(a-z)
integer::N,M,i
real*8::A(M,N),B(N,M),X(N,M)
real*8::btemp(N),xtemp(N)
do i=1,M
   btemp=B(:,i)
   call elgauss(A,btemp,xtemp,N)
   X(:,i)=xtemp
end do
end subroutine mateq

subroutine elgauss(A,b,x,N)
implicit real*8(a-z)
integer::i,k,N
integer::id_max
real*8::A(N,N),b(N),x(N)
real*8::Aup(N,N),bup(N)
real*8::Ab(N,N+1)
real*8::vtemp1(N+1),vtemp2(N+1)
Ab(1:N,1:N)=A
Ab(:,N+1)=b
do k=1,N-1
   elmax=dabs(Ab(k,k))
   id_max=k
   do i=k+1,n
    if(dabs(Ab(i,k))>elmax)then
        elmax=Ab(i,k)
        id_max=i
        end if
        end do
        vtemp1=Ab(k,:)
        vtemp2=Ab(id_max,:)
        Ab(k,:)=vtemp2
        Ab(id_max,:)=vtemp1
do i=k+1,N
    temp=Ab(i,k)/Ab(k,k)
        Ab(i,:)=Ab(i,:)-temp*Ab(k,:)
end do
end do
Aup(:,:)=Ab(1:N,1:N)
bup(:)=Ab(:,N+1)
call uptri(Aup,bup,x,n)
end subroutine elgauss

subroutine uptri(A,b,x,N)
implicit real*8(a-z)
integer::i,j,N
real*8::A(N,N),b(N),x(N)
x(N)=b(N)/A(N,N)
do i=n-1,1,-1
    x(i)=b(i)
        do j=i+1,N
        x(i)=x(i)-a(i,j)*x(j)
        end do
        x(i)=x(i)/A(i,i)
end do
end subroutine uptri
end module inv_mat

module m_bfs
use inv_mat
contains

subroutine solve(x0,N,tol)
implicit real*8(a-z)
integer::i,n,itmax=50
real*8::x1(n),x0(n),y(n),f0(n),f1(n),dx(n)
real*8::v1(n),v2(n),v3(n)
real*8::H0(n,n),H1(n,n),df(n,n)
real*8::m1(n,n),m2(n,n),m3(n,n)
call jac(df,x0)
call inv(df,H0,N)
write(11,101)
write(12,102)
do i=1,itmax
   call func(f0,x0)
   x1=x0-matmul(H0,f0)
   dx=x1-x0
   call func(f1,x1)
   y=f1-f0
   v1=matmul(H0,y)
   t1=vdot(y,v1,n)
   t2=vdot(dx,y,N)
   u=1d0+t1/t2
   m1=vvmat(dx,dx,n)*u
   m2=vvmat(dx,y,n)
   m2=matmul(m2,H0)
   v3=matmul(H0,y)
   m3=vvmat(v3,dx,n)
   t3=vdot(dx,y,n)
   H1=(m1-m2-m3)/t3
   H1=H0+H1
   write(12,103)i,H0
   x0x=1
   H0=H1
   write(11,104)i,x0
   dx2=dsqrt(dx(1)**2+dx(2)**2)
   if(dx2 end do
101  format(/,T5,'BFS方法计算序列',/)
102  format(/,T5,'BFS方法H矩阵序列为:',/)
103  format(T5,'iter=',I4,/,(F16.10/),/)
104  format(I8,3F16.10)
end subroutine solve

function vdot(a,b,N)
integer::i,N
real*8::a(n),b(n),vdot
vdot=0
do i=1,N
   vdot=vdot+a(i)*b(i)
end do
end function vdot

function vvmat(a,b,n)
integer::n,i,j
real*8::a(n),b(n),vvmat(n,n)
do i=1,n
   do j=1,n
     vvmat(i,j)=a(i)*b(j)
   end do
end do
end function  vvmat

subroutine func(f,x)
implicit real*8(a-z)
real*8::x(3),f(3),pi=3.14159263589793
f(1)=3*x(1)-dcos(x(2)*x(3))-1d0/2
f(2)=x(1)**2-81*(x(2)+0.1d0)**2+dsin(x(3))+1.06d0
f(3)=dexp(-x(1)*x(2))+20*x(3)+(10*pi-3d0)/3
end subroutine func

subroutine jac(df,x)
implicit real*8(a-z)
real*8::x(3),df(3,3)
df(1,1)=3d0
df(1,2)=x(3)*dsin(x(2)*x(3))
df(1,3)=x(2)*dsin(x(2)*x(3))
df(2,1)=2*x(1)
df(2,2)=-162*(x(1)+0.1)
df(2,3)=dcos(x(3))
df(3,1)=-x(2)*dexp(-x(1)*x(2))
df(3,2)=-x(1)*dexp(-x(1)*x(2))
df(3,3)=20d0
end subroutine jac

end module m_bfs

program main
use inv_mat
use m_bfs
implicit real*8(a-z)
integer::N=3
real*8::x0(3)
open(unit=11,file='result.txt')
open(unit=12,file='Hmatrix.txt')
x0=(/0.1d0,0.1d0,-0.1d0/)
call solve(x0,n,1d-8)
end program main

5楼: Originally posted by exabyss916 at 2011-11-19 10:45:59:
能不能简单说一下为什么尽可能不要用 real*4、real*8 这些编译器扩展的用法？

1楼: Originally posted by zyj8119 at 2011-11-11 23:00:16:
[code]module inv_mat
contains
subroutine inv(A,invA,N)
implicit real*8(a-z)
integer::n
integer::i
real*8::A(n,n),invA(n,n),E(n,n)
E=0
do i=1,n
   E(i,i)=1
end do
call mateq(A,B,X,N,N)
e ...

2楼: Originally posted by maomao1210 at 2011-11-11 23:15:23:
大哥；
子程序 inv, X是数组，你忘记定义了。X（n,n）

3楼: Originally posted by zyj8119 at 2011-11-11 23:29:18:
谢谢提醒，修改过来了，编译通过。。