使用LAPACK复现Schur分解结果与MATLAB不匹配的技术求助

I've gone through your problem and code, and let's break down the issues you're facing with both the dgehrd/dhseqr pipeline and the dgees sorting option.

1. Why dgehrd + dhseqr Isn't Matching MATLAB #

Your pipeline has a few critical mistakes that are throwing off the results:

a. Incorrect ilo/ihi Values #

LAPACK uses 1-based indexing, so for an n×n matrix, you need to set ilo=1 and ihi=nrows (not nrows-1) to process the entire matrix. Your original setting skipped the last row/column, which broke the decomposition.

b. Missing Proper Orthogonal Transformation Accumulation #

You mentioned calling dorghr but seeing no effect—this is likely because you didn't initialize the Z matrix to the identity before calling dorghr. The dorghr function applies the Householder transformations from dgehrd to an existing matrix; starting with identity ensures you build the full orthogonal matrix U.

c. Potential Storage Order Confusion #

LAPACK uses column-major (Fortran-style) storage, while your C code uses row-major. When copying results from Z (U) and val (T) to your output arrays, you have the transpose right, but make sure your input matrix val is passed in column-major order. If your input A is row-major, you need to transpose it first before feeding it to LAPACK.

Fixed dgehrd + dhseqr Code Snippet #

Here's the corrected pipeline:

int ilo = 1, ihi = nrows; // Correct 1-based bounds for full matrix
double *tau = (double*)calloc(nrows-1, sizeof(double));
double wkopt;
int lwork = -1;
int info;

// Step 1: Generate Householder reductions
LAPACK_dgehrd(&nrows, &ilo, &ihi, val, &nrows, tau, &wkopt, &lwork, &info);
if (info != 0) {
    fprintf(stderr, "dgehrd failed with info=%d\n", info);
    exit(1);
}
lwork = (int)wkopt;
double *work = (double*)malloc(lwork * sizeof(double));
LAPACK_dgehrd(&nrows, &ilo, &ihi, val, &nrows, tau, work, &lwork, &info);

// Step 2: Initialize Z to identity, then accumulate orthogonal transforms
for (int i=0; i<nrows*nrows; i++) Z[i] = 0.0;
for (int i=0; i<nrows; i++) Z[i*nrows + i] = 1.0;
LAPACK_dorghr(&nrows, &ilo, &ihi, Z, &nrows, tau, work, &lwork, &info);
if (info != 0) {
    fprintf(stderr, "dorghr failed with info=%d\n", info);
    exit(1);
}

// Step 3: Compute Schur form T
LAPACK_dhseqr("S", "V", &nrows, &ilo, &ihi, val, &nrows, wr, wi, Z, &nrows, work, &lwork, &info);
if (info > 0) {
    fprintf(stderr, "dhseqr failed to converge at eigenvalue %d\n", info);
    exit(1);
}

// Copy column-major results to row-major U/T (your existing code here is correct)
for (int i=0; i<nrows; i++) {
    for (int j=0; j<nrows; j++) {
        U[i*nrows+j] = Z[j*nrows+i];
        T[i*nrows+j] = val[j*nrows+i];
    }
}

2. Fixing dgees Sorting & Crash Issues #

Your dgees setup has two main problems:

a. Misconfigured sdim Parameter #

The sdim variable is an output parameter—it tells you how many eigenvalues matched your SELECT criteria. You shouldn't set it to nrows-1 upfront; initialize it to 0 instead. Setting it to a non-zero value as input can cause memory corruption and crashes.

b. Incorrect SELECT Function Logic #

MATLAB's default schur(A) sorts Schur blocks by descending real part of eigenvalues (with complex conjugate blocks grouped together). Your current SELECT function only checks if an eigenvalue is exactly (1,0), which doesn't match MATLAB's behavior.

If you want to replicate MATLAB's default order, you need a SELECT function that prioritizes eigenvalues with larger real parts. For example, here's a version that selects eigenvalues with real parts greater than the average real part (adjust based on your specific matrix):

// Global variable to store real parts of all eigenvalues (compute before calling dgees)
double *all_wr;
int nrows_global; // Pass nrows via global to avoid function parameter issues

lapack_logical SELECT(double er, double ei) {
    // Select eigenvalues with real part >= average real part (matches MATLAB's descending order intent)
    double avg = 0.0;
    for (int i=0; i<nrows_global; i++) avg += all_wr[i];
    avg /= nrows_global;
    return (er >= avg) ? 1 : 0;
}

Note: Before calling dgees, compute all eigenvalues with dgeev to populate all_wr, and set nrows_global to your matrix size.

Fixed dgees Call #

Here's the corrected dgees setup:

int sdim = 0; // Initialize to 0 (output parameter)
lapack_logical *bwork = (lapack_logical*)calloc(nrows, sizeof(lapack_logical));
double wkopt;
int lwork = -1;
int info;

// First query optimal workspace size
LAPACK_dgees("V", "S", SELECT, &nrows, val, &nrows, &sdim, wr, wi, Z, &nrows, &wkopt, &lwork, bwork, &info);
lwork = (int)wkopt;
double *work = (double*)malloc(lwork * sizeof(double));

// Actual decomposition with sorting
LAPACK_dgees("V", "S", SELECT, &nrows, val, &nrows, &sdim, wr, wi, Z, &nrows, work, &lwork, bwork, &info);

General Critical Notes #

• Column-Major vs Row-Major: Always ensure your input matrix is in column-major order when passing to LAPACK. If your input is row-major (C-style), transpose it first.
• Workspace Allocation: Never skip the workspace query step (lwork=-1). Using insufficient workspace can lead to crashes or incorrect results.
• Indexing: LAPACK uses 1-based integers for parameters like ilo, ihi, and matrix dimensions—don't mix with C's 0-based indexing.

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