要用到之前發的解上三角矩陣和下三角矩陣方程的模塊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

列選主元的高斯消元法的Fortran程序