I'm going to show you an example how to compute generalized eigenvalue problem in C with CLapack.
Here I have two matrices.
#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];
And I want to use the routine sggev to compute the generalized eigenvalue problem
/* 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;
The output is
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
Compared to MATLAB code:
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)
The outputs are
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
Question:
The eigenvalues from dr and di are the same as D matrix of MATLAB.
But the eigenvectors V compared to wr and wi is almost the same, very close, but it's about 0.01 in error.
Why? Am I using the Clapack routine sggev not correctly?
