使用dgeqrf进行QR分解的Fortran程序问题及秩相关疑问

Hi there, let's tackle your two questions about QR decomposition in Fortran using LAPACK's dgeqrf subroutine step by step.

Fixing the Runtime Crash with dgeqrf #

Your code compiles but crashes at runtime due to a few critical mistakes. Let's break them down and fix the code:

Key Issues in Your Original Code #

1. Wrong Parameter Order in dgeqrf Call
   The standard dgeqrf signature requires the leading dimension LDA right after the matrix A. You skipped this argument, which shifts all subsequent parameters into the wrong positions—this is the main cause of your runtime error. The correct call order is:

   call dgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
   

2. Incorrect Array Dimensions

   • TAU should be a 1D array (not 2D) with length min(M, N)—it stores scalar factors for the elementary reflectors used in the decomposition.
   • WORK is also a 1D workspace array (not 2D) used by LAPACK for temporary calculations.

3. Unallocated Arrays
   You declared TAU and WORK as allocatable but never allocated memory for them. For WORK, it's best practice to first query the optimal workspace size (by setting LWORK=-1) to avoid undersizing.

Corrected Working Code #

program decomposition
    implicit none
    real,dimension(2,2) :: A_mat, R, A
    integer :: M, N, LDA, LWORK, INFO
    real,allocatable, dimension(:) :: TAU, WORK

    external dgeqrf

    ! Initialize matrix dimensions
    M = 2
    N = 2
    LDA = 2
    INFO = 0

    ! Populate input matrix
    A_mat(1,1) = 4
    A_mat(1,2) = 1
    A_mat(2,1) = 3
    A_mat(2,2) = 1
    A = A_mat  ! Make a copy to preserve original matrix

    ! Allocate TAU (size = min(M,N) for square matrices)
    allocate(TAU(min(M,N)))

    ! Query optimal workspace size for WORK
    LWORK = -1
    allocate(WORK(1))
    call dgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
    LWORK = int(WORK(1))
    deallocate(WORK)
    allocate(WORK(LWORK))

    ! Perform actual QR decomposition
    call dgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)

    ! Extract upper triangular matrix R from the output A
    ! Note: dgeqrf stores reflector data in the lower triangle of A
    R = 0.0
    do j = 1, N
        do i = 1, min(j, M)
            R(i,j) = A(i,j)
        end do
    end do

    ! Print results
    print *, "Upper Triangular Matrix R:"
    print *, R(1,:)
    print *, R(2,:)
    print *, "INFO Code (0 = Success):", INFO

    ! Clean up allocated arrays
    deallocate(TAU, WORK)
end program decomposition

Quick Notes on the Fixes #

• Added implicit none to catch undeclared variables (always recommended in Fortran).
• Fixed the dgeqrf parameter order to include LDA.
• Adjusted TAU and WORK to 1D arrays and properly allocated them, including querying the optimal workspace size for efficiency.
• Explicitly extracted the upper triangular R matrix (since dgeqrf overwrites A with both reflector data and the upper triangular R).

QR Decomposition Behavior When Matrix A Has Zero Elements #

To answer your second question clearly:

• Isolated zero elements in A do not cause a "zero decomposition" or rank reduction. QR decomposition (A = Q*R) exists for any real matrix (square or rectangular), no matter how many individual zeros it has.
• Rank reduction only happens if A itself is rank-deficient (i.e., its columns are linearly dependent). For example:
  • If A has an entire column of zeros, the corresponding column in R will also be all zeros, making R rank-deficient (matching the rank of A).
  • If A has linearly dependent columns (even without explicit zeros), R will have one or more zero (or numerically near-zero) diagonal elements, indicating reduced rank.
• Q is always an orthogonal matrix (for square A) or has orthonormal columns (for rectangular A), so Q will never be singular or "zeroed out". R remains upper triangular, with a rank equal to that of A.

内容的提问来源于stack exchange，提问作者Nobody