主题:矩阵求逆出错,方法是参考了根据徐士良编著的Fortran常用算法程序集-第二版上的高斯约当消去法求逆
矩阵为 2 1 这是一个可逆矩阵手算可以得出其逆矩阵为 3/2 -1/2
4 3 -2 1
根据徐士良编著的Fortran常用算法程序集-第二版上的高斯约当消去法求逆写的程序,结果显示不可逆(erro not inv)。求高手相助!不甚感谢啊!程序如下 :
program main
DIMENSION a(2,2),b(2,2),c(2,2),is(2),js(2)
DATA a/2.0,4.0,1.0,3.0/
do i=1,2
do j=1,2
b(i,j)=a(i,j)
enddo
enddo
call brinv(a,2,l,is,js)
if(l.ne.0)then
write(*,10)((a(i,j),j=1,2),i=1,2)
endif
10 format(1x,2f3.1)
end program
subroutine brinv(a,n,l,is,js)
dimension a(n,n),is(n),js(n)
l=1
do k=1,n
d=0.0
do i=k,n
do j=k,n
if(abs(a(i,j)).gt.d) then
d=abs(a(i,j))
is(k)=i
js(k)=j
endif
enddo
enddo
if(d+1.0.eq.1.0)then
l=0
write(*,20)
return
endif
20 format(1x,'erro**not inv')
do j=1,n
t=a(k,j)
a(a,j)=a(is(k),j)
a(is(k),j)=t
enddo
do i=1,n
t=a(i,k)
a(i,k)=a(i,js(k))
a(i,js(k))=t
enddo
a(k,k)=1/a(k,k)
do j=1,n
if(j.ne.k)then
a(k,j)=a(k,j)*a(k,k)
endif
enddo
do i=1,n
if(i.ne.k)then
do j=1,n
if(j.ne.k)then
a(i,j)=a(i,j)-a(i,k)*a(k,j)
endif
enddo
endif
enddo
do i=1,n
if(i.ne.k)then
a(i,k)=-a(i,k)*a(k,k)
end if
enddo
enddo
do k=n,1,-1
do j=1,n
t=a(k,j)
a(k,j)=a(js(k),j)
a(js(k),j)=t
enddo
do i=1,n
t=a(i,k)
a(i,k)=a(i,is(k))
a(i,is(k))=t
enddo
enddo
return
end
4 3 -2 1
根据徐士良编著的Fortran常用算法程序集-第二版上的高斯约当消去法求逆写的程序,结果显示不可逆(erro not inv)。求高手相助!不甚感谢啊!程序如下 :
program main
DIMENSION a(2,2),b(2,2),c(2,2),is(2),js(2)
DATA a/2.0,4.0,1.0,3.0/
do i=1,2
do j=1,2
b(i,j)=a(i,j)
enddo
enddo
call brinv(a,2,l,is,js)
if(l.ne.0)then
write(*,10)((a(i,j),j=1,2),i=1,2)
endif
10 format(1x,2f3.1)
end program
subroutine brinv(a,n,l,is,js)
dimension a(n,n),is(n),js(n)
l=1
do k=1,n
d=0.0
do i=k,n
do j=k,n
if(abs(a(i,j)).gt.d) then
d=abs(a(i,j))
is(k)=i
js(k)=j
endif
enddo
enddo
if(d+1.0.eq.1.0)then
l=0
write(*,20)
return
endif
20 format(1x,'erro**not inv')
do j=1,n
t=a(k,j)
a(a,j)=a(is(k),j)
a(is(k),j)=t
enddo
do i=1,n
t=a(i,k)
a(i,k)=a(i,js(k))
a(i,js(k))=t
enddo
a(k,k)=1/a(k,k)
do j=1,n
if(j.ne.k)then
a(k,j)=a(k,j)*a(k,k)
endif
enddo
do i=1,n
if(i.ne.k)then
do j=1,n
if(j.ne.k)then
a(i,j)=a(i,j)-a(i,k)*a(k,j)
endif
enddo
endif
enddo
do i=1,n
if(i.ne.k)then
a(i,k)=-a(i,k)*a(k,k)
end if
enddo
enddo
do k=n,1,-1
do j=1,n
t=a(k,j)
a(k,j)=a(js(k),j)
a(js(k),j)=t
enddo
do i=1,n
t=a(i,k)
a(i,k)=a(i,is(k))
a(i,is(k))=t
enddo
enddo
return
end