参照宋叶志的《Fortran95 2003 科学计算与工程》
1 module tri_eq 2 !---------------------------------------------module comment 3 ! Author : Huang zc 4 ! Date : 2017-6-27 5 !----------------------------------------------------------- 6 ! Description: 7 ! Solve upper triangular system with back-substitution 8 ! Solve lower triangular system with forward-substitution 9 !----------------------------------------------------------- 10 ! 11 contains 12 13 subroutine uptri(A,b) 14 !---------------------------------------- subroutine comment 15 ! Purpose : Solve Ax=b where A is upper triangular 16 ! solution is stored in b 17 !----------------------------------------------------------- 18 ! Input: 19 ! 1.A(N,N) : coefficient matrix 20 ! 2.b(N) : right-hand side vector 21 ! Output: 22 ! 1.b(N) : solution of matrix equation 23 !---------------------------------------------define variable 24 implicit none 25 integer,parameter::iwp = selected_real_kind(15) 26 real(iwp),intent(in),allocatable::A(:,:) 27 real(iwp),intent(inout),allocatable::b(:) 28 integer::i,j,N 29 real(iwp)::tmp 30 !------------------------------------------------------------- 31 N = size(b) 32 !--------------------------------------------print information 33 print*," Subroutine uptri is being called......" 34 print*," Input upper triangular matrix and rhs vector:..." 35 do i = 1,N 36 do j = 1,N 37 write(*,"(f10.4)",advance="no")A(i,j) 38 enddo 39 write(*,"(a3)",advance="no")"|" 40 write(*,"(f10.4)")b(i) 41 enddo 42 !---------------------------------------------------main body 43 b(N) = b(N)/A(N,N) 44 do i = N-1,1,-1 45 tmp = dot_product(A(i,i+1:N),b(i+1:N)) 46 b(i) = b(i)-tmp 47 b(i) = b(i)/A(i,i) 48 enddo 49 !-------------------------------------------------print result 50 print*," Subroutine uptri end......" 51 print*," Solution of the matrix equation:" 52 do i = 1,N 53 write(*,"(f15.8)")b(i) 54 enddo 55 !------------------------------------------------------------- 56 end subroutine uptri 57 58 subroutine lowertri(A,b) 59 !---------------------------------------- subroutine comment 60 ! Purpose : Solve Ax=b where A is lower triangular 61 ! solution is stored in b 62 !----------------------------------------------------------- 63 ! Input: 64 ! 1.A(N,N) : coefficient matrix 65 ! 2.b(N) : right-hand side vector 66 ! Output: 67 ! 1.b(N) : solution of matrix equation 68 !---------------------------------------------define variable 69 implicit none 70 integer,parameter::iwp = selected_real_kind(15) 71 real(iwp),intent(in),allocatable::A(:,:) 72 real(iwp),intent(inout),allocatable::b(:) 73 integer::i,j,N 74 real(iwp)::tmp 75 !------------------------------------------------------------- 76 N = size(b) 77 !--------------------------------------------print information 78 print*," subroutine lowertri is being called......" 79 print*," Input lower triangular matrix and rhs vector:..." 80 do i = 1,N 81 do j = 1,N 82 write(*,"(f10.4)",advance="no")A(i,j) 83 enddo 84 write(*,"(a3)",advance="no")"|" 85 write(*,"(f10.4)")b(i) 86 enddo 87 !---------------------------------------------------main body 88 b(1) = b(1)/A(1,1) 89 do i = 2,N 90 tmp = dot_product(A(i,1:i-1),b(1:i-1)) 91 b(i) = b(i)-tmp 92 b(i) = b(i)/A(i,i) 93 enddo 94 !-------------------------------------------------print result 95 print*," Subroutine lowertri end......" 96 print*," Solution of the matrix equation:" 97 do i = 1,N 98 write(*,"(f15.8)")b(i) 99 enddo 100 !------------------------------------------------------------- 101 end subroutine lowertri 102 !------------------------------------------------------------- 103 end module tri_eq