作业帮求fortran高手帮我看一下这段用高斯乔丹消元法解方程组的代码为什么总是得不到想要的结果这里的线性方程组对应的系数矩阵是希伯特矩阵这是我的f90文件下载地址 希望有人能回答一下这
来源：学生作业帮助网 编辑：作业帮 时间：2024/04/25 21:25:11
求fortran高手帮我看一下这段用高斯乔丹消元法解方程组的代码为什么总是得不到想要的结果这里的线性方程组对应的系数矩阵是希伯特矩阵这是我的f90文件下载地址 希望有人能回答一下这
求fortran高手帮我看一下这段用高斯乔丹消元法解方程组的代码为什么总是得不到想要的结果
这里的线性方程组对应的系数矩阵是希伯特矩阵
这是我的f90文件下载地址 希望有人能回答一下这个问题 感激不尽
求fortran高手帮我看一下这段用高斯乔丹消元法解方程组的代码为什么总是得不到想要的结果这里的线性方程组对应的系数矩阵是希伯特矩阵这是我的f90文件下载地址 希望有人能回答一下这
主要问题：
调用gaussj时,实参与虚参类型不一致.您在主程序中定义数组a和b是双精度8个字节的,而在子程序中却是默认的4个字节.
我输入x(1:3)=(5.5 ,5.5, 5.5)时输出结果如下图：
SUBROUTINE gaussj(a,n,np,b,m,mp) !引入GaussJ消元法
INTEGER m,mp,n,np,NMAX
REAL*8 a(np,np),b(np,mp)
PARAMETER (NMAX=50)
INTEGER i,icol,irow,j,k,l,ll,indxc(NMAX),indxr(NMAX),ipiv(NMAX)
REAL big,dum,pivinv
do 11 j=1,n
ipiv(j)=0
11 continue
do 22 i=1,n
big=0.
do 13 j=1,n
if(ipiv(j).ne.1)then
do 12 k=1,n
if (ipiv(k).eq.0) then
if (abs(a(j,k)).ge.big)then
big=abs(a(j,k))
irow=j
icol=k
endif
else if (ipiv(k).gt.1) then
pause 'singular matrix in gaussj'
endif
12 continue
endif
13 continue
ipiv(icol)=ipiv(icol)+1
if (irow.ne.icol) then
do 14 l=1,n
dum=a(irow,l)
a(irow,l)=a(icol,l)
a(icol,l)=dum
14 continue
do 15 l=1,m
dum=b(irow,l)
b(irow,l)=b(icol,l)
b(icol,l)=dum
15 continue
endif
indxr(i)=irow
indxc(i)=icol
if (a(icol,icol).eq.0.) pause 'singular matrix in gaussj'
pivinv=1./a(icol,icol)
a(icol,icol)=1.
do 16 l=1,n
a(icol,l)=a(icol,l)*pivinv
16 continue
do 17 l=1,m
b(icol,l)=b(icol,l)*pivinv
17 continue
do 21 ll=1,n
if(ll.ne.icol)then
dum=a(ll,icol)
a(ll,icol)=0.
do 18 l=1,n
a(ll,l)=a(ll,l)-a(icol,l)*dum
18 continue
do 19 l=1,m
b(ll,l)=b(ll,l)-b(icol,l)*dum
19 continue
endif
21 continue
22 continue
do 24 l=n,1,-1
if(indxr(l).ne.indxc(l))then
do 23 k=1,n
dum=a(k,indxr(l))
a(k,indxr(l))=a(k,indxc(l))
a(k,indxc(l))=dum
23 continue
endif
24 continue
return
end
program main
implicit none
integer,parameter :: row=3
integer,parameter :: col=3
integer nn !关于x
integer mm !关于b
integer ii
integer jj
integer kk
real*8:: HI(row,col)
real*8 AI(row,col) !储存矩阵HI的值
real*8 x(row)
real*8 b(row)
real*4 iii !整数型与浮点型的转换
real*4 jjj
do ii=1,row,1 !给H矩阵赋值
do jj=1,col,1
iii=ii
jjj=jj
HI(ii,jj)=1/(iii+jjj-1)
AI(ii,jj)=HI(ii,jj) !将初始矩阵H的值赋给替换矩阵A
write(*,"(A3,I1,I1,A3,F9.4)")"H(",ii,jj,")=",HI(ii,jj)
enddo
enddo
write(*,*) "The substitude matrix is: " !显示替换矩阵A
do ii=1,row,1
do jj=1,col,1
write(*,"(A3,I1,I1,A3,F9.4)")"A(",ii,jj,")=",AI(ii,jj)
enddo
enddo
write(*,*) "Enter the x: " !给x赋值
do nn=1,row,1
read(*,*) x(nn)
enddo
do mm=1,row,1 !求b的值
do jj=1,col,1
b(mm)=b(mm)+HI(mm,jj)*x(jj)
enddo
write(*,"(A3,I1,A3,F9.4)")"b(",mm,")=",b(mm)
enddo !以上代码将矩阵H、解X、常数b均已确定
call gaussj(HI,row,row,b,row,row)
write(*,*) ""
write(*,"(A20)") "The result is : " !显示经gaussj消元法得到的矩阵H
do ii=1,row,1
do jj=1,col,1
write(*,"(A3,I1,I1,A3,F9.4)")"H(",ii,jj,")=",HI(ii,jj)
enddo
enddo
do mm=1,row,1
write(*,*)"b(",mm,")=",b(mm)
enddo
do ii=1,row,1
x(ii)=0
enddo
write(*,*) "The solution vectors are: " !显示解
do ii=1,row,1
do kk=1,col,1
x(ii)=x(ii)+AI(ii,kk)*b(kk)
enddo
write(*,*)"x(",ii,")=",x(ii)
enddo
stop
end