|
|
!1. Solving linear equations by QR decomposition/Gram-Schmidt and
! QR decomposition/back-substitution.
!2. Calculating the determinant of a matrix.
program LIN_EQ
implicit none
integer, parameter :: dbl=SELECTED_REAL_KIND(14), n=3, m=2
real(KIND=dbl) :: A(n,m), Q(n,m), R(m,m), QR(n,m), A_old(n,m), &
& b(n), x(m), DET
integer :: i
A(1,1)=1 ; A(1,2)=1
A(2,1)=0 ; A(2,2)=1
A(3,1)=1 ; A(3,2)=1
b(1)=6 ; b(2)=4 ; b(3)=6
A_old = A
CALL QR_GS(A, n, m, Q, R)
!emuch CALL QR_BACK(Q, R, n, m, b, x)
CALL QR_BACK(Q, R, b, x)
!emuch CALL QR_DETERMINANT(R, m, DET)
CALL QR_DETERMINANT(R, DET)
call matout("A is", A, n, m)
call matout("b is", b, n, 1)
call matout("Q is", Q, n, m)
call matout("R is", R, m, m)
QR = MATMUL(Q,R)
call matout("Q*R is", QR, n, m)
call matout("x is", x, m, 1)
if (n/=m) then
print '("The determinant of the matrix R is", f10.6)', DET
else
print '("The determinant of the matrix A is", f10.6)', DET
endif
contains !emuch
!Subroutine for QR decomposition by a modified Gram-Schmidt method
subroutine QR_GS(A, n, m, Q, R)
Implicit none
integer, parameter :: dbl=SELECTED_REAL_KIND(14)
integer :: n, m, i, j, k, l, i1
real(KIND=dbl) :: A(n,m), Q(n,m), R(m,m)
do i = 1, m
R(i,i)=SQRT(DOT_PRODUCT(A(1:n,i),A(1:n,i)))
do j = 1, n
Q(j,i) = A(j,i)/R(i,i)
enddo
do j = i+1, m
R(i,j)=DOT_PRODUCT(Q(1:n,i), A(1:n,j))
do k = 1, n
A(k,j) = A(k,j)-Q(k,i)*R(i,j)
enddo
enddo
enddo
end subroutine QR_GS
!Subroutine for solving linear equations Ax=b => QRx=b by QR back-substitution.
!emuch subroutine QR_BACK(Q, R, n, m, b, x)
subroutine QR_BACK(Q, R, b, x)
implicit none
integer, parameter :: dbl=SELECTED_REAL_KIND(14)
real(KIND=dbl) :: Q(n,m), R(m,m), x(m), b(n), Q_T(m,n), b_ny(m)
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
integer :: i, k
Q_T = TRANSPOSE(Q)
b_ny = MATMUL(Q_T, b)
x = 0.0
do i = m, 1, -1
do k = i+1, m
x(i) = x(i) - R(i,k)*x(k)
enddo
x(i) = (x(i) + b_ny(i))/R(i,i)
enddo
end subroutine QR_BACK
!Subroutine for calculating the determinant of a matrix by QR decomposition.
!emuch subroutine QR_DETERMINANT(R, m, DET)
subroutine QR_DETERMINANT(R, DET)
implicit none
integer, parameter :: dbl=SELECTED_REAL_KIND(14)
real(KIND=dbl) :: R(m,m), DET
integer :: i
DET = 1
do i = 1, m
DET = DET*R(i,i)
enddo
end subroutine QR_DETERMINANT
end program LIN_EQ |
|
回复此楼