在 C 中使用 CLapack 计算广义特征值问题

问：

我将向你展示一个示例，如何使用 CLapack 在 C 中计算广义特征值问题。

这里我有两个矩阵。

#define row 5
#define column 5

float A[row * column] = { 9.5605,   6.8163, - 1.9425, - 5.7558,   3.5550,
     0.8163,   9.7506,   0.2625, - 3.7964,   4.7512,
    - 0.9425,   0.2625,   3.4206,   2.9018,   0.4519,
    - 0.7558, - 3.7964,   2.9018,   7.0731, - 2.5701,
     0.5550,   4.7512,   0.4519, - 2.5701,   2.6984 };


float B[row * column] = { 1.4511,   0.2932, - 1.1534, - 0.3546,   0.4427,
     0.2932,   5.8465, - 1.0767,   1.7268,   1.1316,
    - 1.1534, - 1.0767,   7.5729, - 4.1207,   2.6296,
    - 0.3546,   1.7268, - 4.1207,   5.7797, - 3.5050,
     0.4427,   1.1316,   2.6296, - 3.5050,   4.5628 };

/* Eigenvalues real */
float dr[row];

/* Eigenvectors real */
float wr[row * row];

/* Eigenvalues imaginary */
float di[row];

/* Eigenvectors imaginary */
float wi[row * row];

我想使用例程来计算广义特征值问题sggev

    /* Settings */
    integer n = row, lda = row, ldb = row, info, lwork;
    real* beta = (real*)malloc(n * sizeof(float));
    real wkopt;

    /* Eigenvectors */
    real* vl = (real*)malloc(n * lda * sizeof(real));
    real* vr = (real*)malloc(n * lda * sizeof(real));

    /* Copy */
    float* Bcopy = (float*)malloc(row * row * sizeof(float));
    memcpy(Bcopy, B, row * row * sizeof(float));
    float* Acopy = (float*)malloc(row * row * sizeof(float));
    memcpy(Acopy, A, row * row * sizeof(float));

    /* Allocate memory */
    lwork = -1;
    sggev_("V", "N", &n, Acopy, &lda, Bcopy, &ldb, dr, di, beta, vl, &n, vr, &n, &wkopt, &lwork, &info);

    /* Compute */
    lwork = (integer)wkopt;
    real* work = (real*)malloc(lwork * sizeof(real));
    sggev_("V", "N", &n, Acopy, &lda, Bcopy, &ldb, dr, di, beta, vl, &n, vr, &n, work, &lwork, &info);

    /* Compute the eigenvalues */
    size_t i;
    for (i = 0; i < row; i++) {
        /* Check if it's a real eigenvalue */
        if (fabsf(di[i]) < MIN_VALUE) {
            dr[i] = dr[i] / beta[i];
            di[i] = 0.0f;
        }
        else {
            dr[i] = dr[i] / beta[i];
            di[i] = di[i] / beta[i];
        }
    }

    /* Fill the eigenvectors in row major */
    size_t j, s = 0, t = 0;
    memset(wi, 0, row * row * sizeof(float));
    for (i = 0; i < n; i++) {
        if (fabsf(di[i]) < MIN_VALUE) {
            s = t;
            for (j = 0; j < row; j++) {
                wr[j * row + i] = vl[s++];
            }
            t = s;
        }
        else {
            t = s;
            for (j = 0; j < row; j++) {
                wr[j * row + i] = vl[t++];
            }
            for (j = 0; j < row; j++) {
                wi[j * row + i] = vl[t++];
            }
        }
    }

    /* Free */
    free(beta);
    free(vl);
    free(vr);
    free(Bcopy);
    free(Acopy);

    /* Return status */
    bool status = info == 0;

输出为

Eigenvalues real dr:
7.1193342
7.1193337
1.4855472
0.0041028
0.1091108

Eigenvalues imaginary di:
2.2437241
-2.2437241
0.0000000
0.0000000
0.0000000

Eigenvectors real wr:
-0.7547432      -0.7547432      0.6629440       -0.0018632      -0.0072472
0.2152242       0.2152242       -1.0000000      -0.3431580      0.6721919
-0.3558549      -0.3558549      -0.5982113      -0.3981720      -1.0000000
-0.5564440      -0.5564440      -0.6243958      0.3420532       0.7263851
-0.1728931      -0.1728931      -0.6662152      1.0000000       -0.6148371

Eigenvectors imaginary wi:
0.2452568       0.2452568       0.0000000       0.0000000       0.0000000
0.1706368       0.1706368       0.0000000       0.0000000       0.0000000
-0.0627999      -0.0627999      0.0000000       0.0000000       0.0000000
-0.3060447      -0.3060447      0.0000000       0.0000000       0.0000000
-0.2216028      -0.2216028      0.0000000       0.0000000       0.0000000

与MATLAB代码相比：

A = [9.5605   6.8163  -1.9425  -5.7558   3.5550
         0.8163   9.7506   0.2625  -3.7964   4.7512
        -0.9425   0.2625   3.4206   2.9018   0.4519
        -0.7558  -3.7964   2.9018   7.0731  -2.5701
         0.5550   4.7512   0.4519  -2.5701   2.6984];

    B = [1.4511   0.2932  -1.1534  -0.3546   0.4427
         0.2932   5.8465  -1.0767   1.7268   1.1316
        -1.1534  -1.0767   7.5729  -4.1207   2.6296
        -0.3546   1.7268  -4.1207   5.7797  -3.5050
         0.4427   1.1316   2.6296  -3.5050   4.5628];

    [V, D] = eig(A, B)

输出为

V =

 Columns 1 through 3:

     -0.74041 -     0.25959i     -0.74041 +     0.25959i      0.66294 +           0i
      0.21661 -     0.16376i      0.21661 +     0.16376i           -1 +           0i
     -0.35316 +    0.053976i     -0.35316 -    0.053976i     -0.59821 +           0i
     -0.55696 +     0.28984i     -0.55696 -     0.28984i      -0.6244 +           0i
     -0.17593 +      0.2151i     -0.17593 -      0.2151i     -0.66621 +           0i

 Columns 4 and 5:

    0.0018631 +           0i    0.0072472 +           0i
      0.34316 +           0i     -0.67219 +           0i
      0.39817 +           0i            1 +           0i
     -0.34205 +           0i     -0.72638 +           0i
           -1 +           0i      0.61484 +           0i

D =

Diagonal Matrix

 Columns 1 through 3:

       7.1193 +      2.2437i                           0                           0
                           0       7.1193 -      2.2437i                           0
                           0                           0       1.4855 +           0i
                           0                           0                           0
                           0                           0                           0

 Columns 4 and 5:

                           0                           0
                           0                           0
                           0                           0
    0.0041026 +           0i                           0
                           0      0.10911 +           0i

问题：

和的特征值与MATLAB的矩阵相同。 但是特征向量与 和 相比几乎相同，非常接近，但误差约为 0.01。drdiDVwrwi

为什么？我是否未正确使用 Clapack 例程？sggev