要用到之前发的解上三角矩阵和下三角矩阵方程的模块tri_eq.f90。
博客园代码不支持fortran格式。。。
1 module lina_gauss 2 !---------------------------------------------module comment 3 ! Author : Huang zc 4 ! Date : 2017-6-27 5 !----------------------------------------------------------- 6 ! Description: 7 ! Solve matrix equations by variable gauss methods 8 !----------------------------------------------------------- 9 ! 10 use tri_eq 11 contains 12 13 subroutine gauss(A,b) 14 !---------------------------------------- subroutine comment 15 ! Purpose : Solve Ax=b by gauss method 16 ! column pivoting is used 17 ! solution is stored in b 18 !----------------------------------------------------------- 19 ! Input: 20 ! 1.A(N,N) : coefficient matrix 21 ! 2.b(N) : right-hand side vector 22 ! Output: 23 ! 1.b(N) : solution of matrix equation 24 !---------------------------------------------define variable 25 implicit none 26 integer,parameter::iwp = selected_real_kind(15) 27 real(iwp),parameter::eps = 3.0d-15 28 real(iwp),intent(inout),allocatable::A(:,:) 29 real(iwp),intent(inout),allocatable::b(:) 30 real(iwp),allocatable::vtemp(:) 31 real(iwp)::tmp 32 integer::i,j,k(1),N 33 !------------------------------------------------------------- 34 N = size(b) 35 allocate(vtemp(N)) 36 !--------------------------------------------print information 37 print*," Subroutine gauss is being called......" 38 print*," Input coefficient matrix and rhs vector:..." 39 do i = 1,N 40 do j = 1,N 41 write(*,"(f10.4)",advance="no")A(i,j) 42 enddo 43 write(*,"(a3)",advance="no")"|" 44 write(*,"(f10.4)")b(i) 45 enddo 46 !---------------------------------------------------main body 47 do j = 1,N-1 48 !---------------------------------pivoting within columns 49 k = maxloc(dabs(A(j:N,j))) 50 k(1) = k(1)+j-1 51 if (k(1) /= j) then 52 vtemp = A(j,:);A(j,:) = A(k(1),:);A(k(1),:) = vtemp 53 tmp = b(j);b(j) = b(k(1));b(k(1)) = tmp 54 print*," Pivoting has been used......" 55 endif 56 !----------------------------------------check if singular 57 if (abs(A(j,j)) < eps) then 58 print*,"The input coefficient matrix seems to be singular!" 59 write(*,110)j 60 stop "Program stop at gauss subroutine!" 61 endif 62 !-----------------------------------eliminate to upper tri 63 do i = j+1,N 64 tmp = -A(i,j)/A(j,j) 65 A(i,j:) = A(i,j:) + tmp*A(j,j:) 66 b(i) = b(i) + tmp*b(j) 67 enddo 68 enddo 69 print*," The coefficient matrix has been changed to upper tri..." 70 call uptri(A,b) 71 !-------------------------------------------------print result 72 print*," Subroutine gauss end......" 73 110 format(t2,'The',i3,' th pivot element is close to zero!') 74 !------------------------------------------------------------- 75 end subroutine gauss 76 77 subroutine gauss2(A,b,x) 78 !---------------------------------------- subroutine comment 79 ! Purpose : Solve Ax=b by gauss method with augment matix 80 ! column pivoting is used 81 ! solution is stored in x 82 !----------------------------------------------------------- 83 ! Input: 84 ! 1.A(N,N) : coefficient matrix 85 ! 2.b(N) : right-hand side vector 86 ! Output: 87 ! 1.x(N) : solution of matrix equation 88 !---------------------------------------------define variable 89 implicit none 90 integer,parameter::iwp = selected_real_kind(15) 91 real(iwp),parameter::eps = 3.0d-15 92 real(iwp),intent(in),allocatable::A(:,:) 93 real(iwp),intent(in),allocatable::b(:) 94 real(iwp),intent(inout),allocatable::x(:) 95 real(iwp),allocatable::Ab(:,:),Aup(:,:),bup(:),vtemp(:) 96 integer::i,j,k(1),N 97 real(iwp)::tmp 98 !------------------------------------------------------------- 99 N = size(b) 100 allocate(Ab(N,N+1),Aup(N,N),bup(N),vtemp(N+1)) 101 Ab(:,1:N) = A;Ab(:,N+1) = b 102 !--------------------------------------------print information 103 print*," Subroutine gauss2 is being called......" 104 print*," Input coefficient matrix and rhs vector:..." 105 do i = 1,N 106 do j = 1,N 107 write(*,"(f10.4)",advance="no")A(i,j) 108 enddo 109 write(*,"(a3)",advance="no")"|" 110 write(*,"(f10.4)")b(i) 111 enddo 112 !---------------------------------------------------main body 113 do j = 1,N-1 114 !---------------------------------pivoting within columns 115 k = maxloc(dabs(Ab(j:N,j))) 116 k(1) = k(1)+j-1 117 if (k(1) /= j) then 118 vtemp = Ab(j,:);Ab(j,:) = Ab(k(1),:);Ab(k(1),:) = vtemp 119 print*," Pivoting has been used......" 120 endif 121 !----------------------------------------check if singular 122 if (abs(Ab(j,j)) < eps) then 123 print*,"The input coefficient matrix seems to be singular!" 124 write(*,110)j 125 stop "Program stop at gauss2 subroutine!" 126 endif 127 !-----------------------------------eliminate to upper tri 128 do i = j+1,N 129 tmp = -Ab(i,j)/Ab(j,j) 130 Ab(i,j:) = Ab(i,j:) + tmp*Ab(j,j:) 131 enddo 132 enddo 133 Aup = Ab(:,1:N);bup = Ab(:,N+1) 134 call uptri(Aup,bup) 135 x = bup 136 !-------------------------------------------------print result 137 print*," Subroutine gauss2 end......" 138 110 format(t2,'The',i3,' th pivot element is close to zero!') 139 !------------------------------------------------------------- 140 end subroutine gauss2 141 142 subroutine iter_solu(A,b,x) 143 !---------------------------------------- subroutine comment 144 ! Purpose : Solve Ax=b by gauss2 method and advance the accuracy 145 ! by iteration of residual equation 146 ! column pivoting is used 147 ! solution is stored in b 148 !----------------------------------------------------------- 149 ! Input: 150 ! 1.A(N,N) : coefficient matrix 151 ! 2.b(N) : right-hand side vector 152 ! Output: 153 ! 1.b(N) : solution of matrix equation 154 !---------------------------------------------define variable 155 implicit none 156 integer,parameter::iwp = selected_real_kind(15) 157 integer,parameter::itmax = 6 158 real(iwp),parameter::eps = 3.0d-15 159 real(iwp),intent(in),allocatable::A(:,:) 160 real(iwp),intent(in),allocatable::b(:) 161 real(iwp),intent(inout),allocatable::x(:) 162 real(iwp),allocatable::x1(:),dx(:),x2(:),db(:) 163 integer::i,N 164 !------------------------------------------------------------- 165 N = size(b) 166 allocate(dx(N),x2(N),db(N)) 167 !--------------------------------------------print information 168 print*," Subroutine iter_solu is being called......" 169 !---------------------------------------------------main body 170 call gauss2(A,b,x1) 171 do i = 1,itmax 172 db = matmul(A,x1) - b 173 call gauss2(A,db,dx) 174 x2 = x1-dx 175 x1 = x2 176 enddo 177 x = x1 178 !-------------------------------------------------print result 179 print*," Solution of the final iteration:" 180 do i = 1,N 181 write(*,"(f15.8)")x(i) 182 enddo 183 print*," Subroutine iter_solu end......" 184 !------------------------------------------------------------- 185 end subroutine iter_solu 186 187 end module lina_gauss