我正在尝试构建一个程序来进行矩阵的 QR 分解
我试过这个：
program QR_Factorization
implicit none
integer, parameter :: n = 3
real(8), dimension(1:9) :: AA =(/1, 6, -4, -51, 167, 24, 4, -68, -41/)
real(8), dimension(n,n) :: A, Q, R
real(8), dimension(n) :: v
real(8) :: alpha
integer :: i, j, k
! Initialize matrix A (replace with your own matrix)
A = reshape( AA, (/n,n/))
print *, "Matrix A:"
do i = 1, n
print *, A(i, :)
end do
! Initialize Q as an identity matrix
Q = identity_matrix(n)
print *, "Matrix Q:"
do i = 1, n
print *, Q(i, :)
end do
! Perform QR Factorization
do k = 1, n-1
call householder(A(:, k:n),k , v)
A(k:n, k:n) = A(k:n, k:n) - 2.0 * outer_product(v, v)* A(k:n, k:n)
Q(:, k:n) = Q(k:n, k:n) - 2.0 * outer_product(v, v)* Q(k:n, k:n)
end do
! Display results
print *, "Matrix A:"
do i = 1, n
print *, A(i, :)
end do
print *, "Matrix Q:"
do i = 1, n
print *, Q(i, :)
end do
print *, "Matrix R:"
R = matmul(transpose(Q), A)
do i = 1, n
print *, R(i, :)
end do
contains
subroutine householder(A, k, v)
real(8), dimension(:,:) :: A
real(8), dimension(:) :: v
real(8) :: alpha, beta
integer, intent(in) :: k
integer :: i
alpha = sqrt(sum(A(k:n, k)**2))
if (A(k, k) < 0.0) alpha = -alpha
beta = 1.0 / (alpha * (alpha - A(k, k)))
v = 0.0
v(k) = 1.0
v(k+1:) = A(k+1:, k) / A(k, k)
v = beta * v
end subroutine householder
function identity_matrix(n) result(mat)
integer, intent(in) :: n
real(8), dimension(n, n) :: mat
integer :: i, j
mat = 0.0
do i = 1, n
mat(i, i) = 1.0
end do
end function identity_matrix
function outer_product(u, v) result(mat)
real(8), dimension(:), intent(in) :: u, v
real(8), dimension(size(u), size(v)) :: mat
integer :: i, j
mat = 0.0
do i = 1, size(u)
do j = 1, size(v)
mat(i, j) = u(i) * v(j)
end do
end do
end function outer_product
end program QR_Factorization
得到了这个：
Matrix A:
1.0000000000000000 -51.000000000000000 4.0000000000000000
6.0000000000000000 167.00000000000000 -68.000000000000000
-4.0000000000000000 24.000000000000000 -41.000000000000000
Matrix Q:
1.0000000000000000 0.0000000000000000 0.0000000000000000
0.0000000000000000 1.0000000000000000 0.0000000000000000
0.0000000000000000 0.0000000000000000 1.0000000000000000
Matrix A:
5.9655462561413701 -50.707220229258688 4.0153087461825523
-4.0153052049990920 161.24773862190591 -69.561492110620350
-4.0153087461825523 24.551114862571886 -40.372235704486464
Matrix Q:
0.99904320336359043 0.96555532108925701 0.0000000000000000
0.0000000000000000 0.0000000000000000 0.98468867571918195
0.0000000000000000 0.0000000000000000 4.2439851328916695E-314
Matrix R:
5.9598384415491488 -50.658703731501653 4.0114669122800590
5.7600649308213958 -48.960626310005544 3.8770027256927961
-3.9538255649188945 158.77882220631733 -68.496413547457081
我知道这是完全错误的。对于愚蠢的错误，我提前表示歉意，我对 Fortran 很陌生。
编辑：
好吧，所以我做了更改（添加辅助矩阵 B 来进行计算并使用 real64），但我仍然得到错误的矩阵：
program QR_Factorization
use iso_Fortran_env, only: real64
implicit none
integer, parameter :: n = 3
real(real64), dimension(1:9) :: AA =(/1, 6, -4, -51, 167, 24, 4, -68, -41/)
real(real64), dimension(n,n) :: A, Q, R, B
real(real64), dimension(n) :: v
real(real64) :: alpha
integer :: i, j, k
! Initialize matrix A (replace with your own matrix)
A = reshape( AA, (/n,n/))
B = A
print *, "Matrix A:"
do i = 1, n
print *, A(i, :)
end do
! Initialize Q as an identity matrix
Q = identity_matrix(n)
print *, "Matrix Q:"
do i = 1, n
print *, Q(i, :)
end do
! Perform QR Factorization
do k = 1, n-1
call householder(B(:, k:n),k , v)
B(k:n, k:n) = B(k:n, k:n) - 2._real64 * outer_product(v, v)* B(k:n, k:n)
Q(:, k:n) = Q(k:n, k:n) - 2._real64 * outer_product(v, v)* Q(k:n, k:n)
end do
! Display results
print *, "Matrix A:"
do i = 1, n
print *, A(i, :)
end do
print *, "Matrix Q:"
do i = 1, n
print *, Q(i, :)
end do
print *, "Matrix R:"
R = matmul(transpose(Q), A)
do i = 1, n
print *, R(i, :)
end do
contains
subroutine householder(A, k, v)
real(real64), dimension(:,:) :: A
real(real64), dimension(:) :: v
real(real64) :: alpha, beta
integer, intent(in) :: k
integer :: i
alpha = sqrt(sum(A(k:n, k)**2))
if (A(k, k) < 0._real64) alpha = -alpha
beta = 1.0 / (alpha * (alpha - A(k, k)))
v = 0._real64
v(k) = 1._real64
v(k+1:) = A(k+1:, k) / A(k, k)
v = beta * v
end subroutine householder
function identity_matrix(n) result(mat)
integer, intent(in) :: n
real(real64), dimension(n, n) :: mat
integer :: i, j
mat = 0._real64
do i = 1, n
mat(i, i) = 1._real64
end do
end function identity_matrix
function outer_product(u, v) result(mat)
real(real64), dimension(:), intent(in) :: u, v
real(real64), dimension(size(u), size(v)) :: mat
integer :: i, j
mat = 0._real64
do i = 1, size(u)
do j = 1, size(v)
mat(i, j) = u(i) * v(j)
end do
end do
end function outer_product
end program QR_Factorization
有人有一个不使用 Lapack 的 QR 分解代码我可以看一下吗？
我整理了你的主程序和v向量的评估。还添加了一个辅助例程来打印矩阵。我没有理会你的其他日常事务。
确保您进行了适当的检查：即 Q 是正交的 (QT.Q=I)，R 是上三角的（您可能仍然需要清理它 - 有些元素很小，而不是 0）并且 QR 等于您的原始矩阵.
以下是基本版本。
program QR_Factorization
use iso_Fortran_env, only: real64
implicit none
integer, parameter :: n = 3
real(real64), dimension(1:9) :: AA = (/1, 6, -4, -51, 167, 24, 4, -68, -41 /)
real(real64), dimension(n,n) :: A, Q, R, Qk
real(real64), dimension(n) :: v
integer k
character( len=* ), parameter :: FMT = "( *( f8.3 ) )"
! Initialize matrix A (replace with your own matrix)
A = reshape( AA, (/ n, n /) )
! Initialise Q and R
Q = identity_matrix( n )
R = A
! Perform QR Factorization
do k = 1, n-1
call householder( R, k, v )
Qk = identity_matrix(n)
Qk(k:n,k:n) = Qk(k:n,k:n) - 2.0 * outer_product( v(k:n), v(k:n) )
Q = matmul( Qk, Q )
R = matmul( Qk, R )
end do
Q = transpose( Q )
! Display results
call printMatrix( "Original matrix ", A, FMT )
call printMatrix( "CHECK QR: ", matmul( Q, R ), FMT )
call printMatrix( "CHECK ORTHOGONALITY (QT.Q): ", matmul( transpose( Q ), Q ), FMT )
call printMatrix( "Q: ", Q, FMT )
call printMatrix( "R: ", R, FMT )
contains
subroutine printMatrix( text, M, fmt )
character(len=*), intent( in ) :: text
real(real64), dimension(:,:), intent(in) :: M
character(len=*), intent( in ) :: fmt
integer i
print *
print *, text
do i = 1, size( M, 1 )
print fmt, M(i,:)
end do
end subroutine printMatrix
subroutine householder( R, k, v )
real(real64), dimension(:,:), intent(in) :: R
real(real64), dimension(:), intent(out) :: v
integer, intent(in) :: k
real(real64) vnorm
v = 0
v(k) = R(k,k) - norm2( R(k:n,k) )
v(k+1:) = R(k+1:,k)
vnorm = norm2( v )
if ( vnorm > 1.0e-10 ) v = v / vnorm
end subroutine householder
function identity_matrix(n) result(mat)
integer, intent(in) :: n
real(real64), dimension(n,n) :: mat
integer :: i
mat = 0._real64
do i = 1, n
mat(i, i) = 1._real64
end do
end function identity_matrix
function outer_product(u, v) result(mat)
real(real64), dimension(:), intent(in) :: u, v
real(real64), dimension(size(u), size(v)) :: mat
integer :: i, j
mat = 0._real64
do i = 1, size(u)
do j = 1, size(v)
mat(i, j) = u(i) * v(j)
end do
end do
end function outer_product
end program QR_Factorization
输出：
Original matrix
1.000 -51.000 4.000
6.000 167.000 -68.000
-4.000 24.000 -41.000
CHECK QR:
1.000 -51.000 4.000
6.000 167.000 -68.000
-4.000 24.000 -41.000
CHECK ORTHOGONALITY (QT.Q):
1.000 -0.000 0.000
-0.000 1.000 -0.000
0.000 -0.000 1.000
Q:
0.137 -0.511 0.849
0.824 0.534 0.188
-0.549 0.674 0.494
R:
7.280 117.443 -32.967
0.000 131.427 -65.986
0.000 0.000 -29.666