黄佳玮 2014301020139 物基二班

8.1 ex0801.f90

要求

编写一个子程序，对任意阶方阵进行上三角化。并编写主程序演示该子程序的功能。

代码

program ex0801

implicit none

real,allocatable::a(:,:)
integer::n,i,j,status
write(*,'(a,$)')'input'
read(*,*) n
allocate(a(n,n),stat=status)
if (status/=0) write(*,*) 'error'
a=0
write(*,*)'input the number'
do i=1,n
  do j=1,n
        read(*,*) a(i,j)
  end do
end do

call upperTriangularMatrix(n,a)
call showResult(n,a)

end

subroutine check(k,n,a)
implicit none
integer::n,i,j,k
real::a(n,n)
if(a(k,k)/=0) then
  return
else
  do i=k,n
    if(a(i,k)/=0) then
      call swap(i,k,n,a)
      return
    end if
  end do
end if
end subroutine check

subroutine swap(i,k,n,a)
implicit none
integer::n,i,j,k
real::a(n,n),temp
do j=1,n
  temp=a(i,j)
  a(i,j)=a(k,j)
  a(k,j)=temp
end do
end subroutine swap

subroutine upperTriangularMatrix(n,a)
implicit none
integer::n,i,j,k
real::a(n,n),temp
do k=1,n-1
  call check(k,n,a)
  do i=k+1,n
    if(a(k,k)/=0) then
      temp=a(i,k)/a(k,k)
      do j=k,n
        a(i,j)=-a(k,j)*temp+a(i,j)
      end do
    end if
  end do
end do
return
end subroutine upperTriangularMatrix

subroutine showResult(n,a)
implicit none
integer::n,i,j
real::a(n,n)
write(*,'(a,i2,a)') 'result'
do i=1,n
    do j=1,n
        write(*,'(f8.1,$)') a(i,j)
    end do
    write(*,*) ' '
end do
return
end subroutine showResult

结果

分析

• 描述算法
  主程序：
  首先建立一个n*n的矩阵，也就是多维数组，动态分配其大小，输入阶数以及矩阵中各元素的值，调用子程序upperTriangularMatrix(n,a)（其中又包括子程序check，check中又包括子程序swap），再调用showResult，结束程序。
  子程序：
  upperTriangularMatrix(n,a)：
  首先检查矩阵对角线上有没有值为零的项， 确认没有之后选取要进行操作的一行，乘以它与下一行的比值后再减去下一行，重复操作，即可得到上三角化的矩阵。
  check(k,n,a)：
  检查矩阵对角线上有没有值为零的项，若没有则回到upperTriangularMatrix，若有，则找到该列上不为零的另一行，进入swap程序。
  swap：
  将check中选定的两行中所对应的每一项交换位置，使对角线上的项不为零。
  showResult(n,a)：
  按照矩阵的形式输出结果。
• 分析缺陷
  输入的形式只能一个数一行，容易出错，如果能在输入数据时就以矩阵的形式读数据会更好一些，也方便校对。另外如果两行数据基本相同，只有一个数字有较小的区别，则计算结果可能会出错。