求fortran高手帮我看一下这段用高斯乔丹消元法解方程组的代码为什么总是得不到想要的结果

求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)=011    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'              endif12          continue          endif13      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)=dum14        continue          do 15 l=1,m            dum=b(irow,l)            b(irow,l)=b(icol,l)            b(icol,l)=dum15        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)*pivinv16      continue        do 17 l=1,m          b(icol,l)=b(icol,l)*pivinv17      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)*dum18          continue            do 19 l=1,m              b(ll,l)=b(ll,l)-b(icol,l)*dum19          continue          endif21      continue22    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))=dum23        continue        endif24    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)enddoenddo            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
再问： 问题出在哪？
再答： 调用gaussj时，实参与虚参类型不一致。您在主程序中定义数组a和b是双精度8个字节的，而在子程序中却是默认的4个字节。 我把你的子程序的第三行改了，结果就正确了。