（原）使用mkl計算特征值和特征向量

轉載請注明出處：

http://www.cnblogs.com/darkknightzh/p/5585271.html

參考文檔：mkl官方文檔

 

lapack_int LAPACKE_sgeev(int matrix_layout, char jobvl, char jobvr, lapack_int n, float* a, lapack_int lda, float* wr, float* wi, float* vl, lapack_int ldvl, float* vr, lapack_int ldvr);

說明：

用於計算n*n實/復非對稱矩陣A的特征值和左/右特征向量。

A的右特征值v滿足：A*v = λ*v，λ為特征值。

A的左特征值u滿足：${{u}^{H}}*A=\lambda *{{u}^{H}}$，λ為特征值。 ${{u}^{H}}$為u的共軛轉置。The computed eigenvectors are normalized to have Euclidean norm equal to 1 and largest component real。

輸入：

matrix_layout：數據存儲格式，行優先- LAPACK_ROW_MAJOR，列優先- LAPACK_COL_MAJOR

jobvl：‘N’，不計算A的左特征值；‘V’，計算A的左特征值。

jobvr：‘N’，不計算A的右特征值；‘V’，計算A的右特征值。

n：矩陣A的階數(n≥ 0).

a：(大小最小 max(1, lda*n)) ，包含n*n 的矩陣 A.

lda：數組a的leading dimension. 最小max(1, n).

ldvl, ldvr：輸出數組vl、vr的leading dimensions

    ldvl≥ 1; ldvr≥ 1.
    If jobvl = 'V', ldvl≥ max(1, n);
    If jobvr = 'V', ldvr≥ max(1, n).

輸出：

a：該緩沖區會被覆蓋

wr, wi：最小max (1, n) 的緩沖區。包含計算到的特征值的實部和虛部。

    Complex conjugate pairs of eigenvalues appear consecutively with the eigenvalue having positive imaginary part first.

vl, vr：

    vl (最小max(1, ldvl*n)) .

    jobvl = 'N'時，vl不使用

    計算到的特征值為實數時，對於第j個特征值：列優先時，第j個特征向量的第i個元素存儲在vl[(i - 1) + (j - 1)*ldvl]（每行存儲對應的特征向量）；行優先時，在vl[(i - 1)*ldvl + (j - 1)]（每列存儲對應的特征向量）。

    計算到的特征值為復數及復矩陣的情況見官方文檔。

    jobvr = 'N'時，不使用vr

    計算到的特征值為實數時，對於第j個特征值：列優先時，第j個特征向量的第i個元素存儲在vr[(i - 1) + (j - 1)*ldvr]（每行存儲對應的特征向量）；行優先時，在vr[(i - 1)*ldvr + (j - 1)]（每列存儲對應的特征向量）。

    計算到的特征值為復數及復矩陣的情況見官方文檔。

返回值：

0：計算成功

-i：第i個參數不合法

i：QR算法未能成功計算所有特征值，並且未計算任何特征向量。elements i+1:n of wr and wi (for real flavors) or w (for complex flavors) contain those eigenvalues which have converged.

 

lapack_int LAPACKE_ssyevd(int matrix_layout, char jobz, char uplo, lapack_int n, float* a, lapack_int lda, float* w);

說明：

計算實對稱矩陣A所有特征值，（特征向量可選）。只計算特征值時，使用Pal-Walker-Kahan variant of the QL or QR算法。需要計算特征向量時，使用divide and conquer算法。

注意：

對於實對稱問題，大部分情況下，默認的選擇是syevr函數，該函數更快，並且內存使用更少。syevd函數需要更多內存，但是更快，特別是對於大矩陣。

輸入：

matrix_layout：數據存儲格式，行優先- LAPACK_ROW_MAJOR，列優先- LAPACK_COL_MAJOR

jobz：‘N’，只計算特征值；‘V’，計算特征值和特征向量。

uplo：‘U’，a存儲了A的上三角部分；‘L’，a存儲了A的下三角部分。

n：矩陣A的階數(n≥ 0).

a：(大小max(1, lda*n)) ，包含實對稱矩陣 A的上三角或下三角部分。通過參數uplo指定。

lda：數組a的leading dimension。最小max(1, n).

ldvl, ldvr：輸出數組vl、vr的leading dimensions

    ldvl≥ 1; ldvr≥ 1.
    If jobvl = 'V', ldvl≥ max(1, n);
    If jobvr = 'V', ldvr≥ max(1, n).

輸出：

w：數組，至少max(1, n)。當返回值為0時，包含升序的特征值。

a：當jobz = 'V'時，該參數會被正交的特征向量覆蓋。

返回值：

0：計算成功

-i：第i個參數不合法

i：jobz = 'N'時，算法未能收斂；i indicates the number of off-diagonal elements of an intermediate tridiagonal form which did not converge to zero；

    jobz = 'V'時，the algorithm failed to compute an eigenvalue while working on the submatrix lying in rows and columns info/(n+1) through mod(info,n+1).

 

https://software.intel.com/en-us/articles/intel-mkl-111-release-notes

中寫道：If LAPACKE_ssyevd fails to evaluate a matrix , workaround is to use LAPACKE_ssyev

lapack_int LAPACKE_ssyev(int matrix_layout, char jobz, char uplo, lapack_int n, float* a, lapack_int lda, float* w);

說明：

計算實對稱矩陣A所有特征值，（特征向量可選）。

注意：

對於實對稱問題，大部分情況下，默認的選擇是syevr函數，該函數更快，並且內存使用更少。syevd函數需要更多內存，但是更快，特別是對於大矩陣。

輸入：

matrix_layout：數據存儲格式，行優先- LAPACK_ROW_MAJOR，列優先- LAPACK_COL_MAJOR

jobz：‘N’，只計算特征值；‘V’，計算特征值和特征向量。

uplo：‘U’，a存儲了A的上三角部分；‘L’，a存儲了A的下三角部分。

n：矩陣A的階數(n≥ 0).

a：(大小max(1, lda*n)) ，包含實對稱矩陣 A的上三角或下三角部分。通過參數uplo指定。

lda：數組a的leading dimension。最小max(1, n).

輸出：

w：數組，至少max(1, n)。當返回值為0時，包含升序的特征值。

a：當jobz = 'V'時，如果info = 0，a中包含正交的特征向量。當jobz = 'N'時，包含對角線的下三角(uplo = 'L')或上三角(uplo = 'U')元素會被覆蓋。

返回值：

0：計算成功

-i：第i個參數不合法

i：算法未能收斂。i indicates the number of elements of an intermediate tridiagonal form which did not converge to zero.

 

lapack_int LAPACKE_ssyevr(int matrix_layout, char jobz, char range, char uplo, lapack_int n, float* a, lapack_int lda, float vl, float vu, lapack_int il, lapack_int iu, float abstol, lapack_int* m, float* w, float* z, lapack_int ldz, lapack_int* isuppz);

說明：

計算實對稱矩陣A所有特征值，（特征向量可選）。

可通過設定特征值的范圍或者特征值的位置來選定特征值（Eigenvalues and eigenvectors can be selected by specifying either a range of values or a range of indices for the desired eigenvalues.）

輸入：

matrix_layout：數據存儲格式，行優先- LAPACK_ROW_MAJOR，列優先- LAPACK_COL_MAJOR

jobz：‘N’，只計算特征值；‘V’，計算特征值和特征向量。

range：'A'，計算所有特征值；'V'，計算vl < w[i]≤vu在半開區間的特征值；'I'，計算索引從il到iu之內的特征值。

uplo：‘U’，a存儲了A的上三角部分；‘L’，a存儲了A的下三角部分。

n：矩陣A的階數(n≥ 0).

a：(大小max(1, lda*n)) ，包含實對稱矩陣 A的上三角或下三角部分。通過參數uplo指定。

lda：數組a的leading dimension。最小max(1, n).

vl, vu：range = 'V'時，計算到的特征值的左右邊界（vl< vu）；range = 'A'或'I'，這兩個參數不使用。

il, iu：當range = 'I'時，最小和最大特征值的索引。

    n > 0時，1 ≤il≤iu≤n

    n = 0時，il=1，iu=0

    range = 'A'或'V'時，不使用這兩個參數。

abstol：當jobz = 'V'時，the eigenvalues and eigenvectors output have residual norms bounded by abstol, and the dot products between different eigenvectors are bounded by abstol.

    If abstol < n *eps*||T||, then n *eps*||T|| is used instead, where eps is the machine precision, and ||T|| is the 1-norm of the matrix T. The eigenvalues are computed to an accuracy of eps*||T|| irrespective of abstol.

ldz：輸出緩沖區z的leading dimension。jobz = 'N'時，ldz≥ 1；jobz = 'V'時，列優先，ldz≥ max(1, n)；行優先，ldz≥ max(1, m)

輸出：

a：包含對角線的下三角(uplo = 'L')或上三角(uplo = 'U')元素會被覆蓋。

m：總共找到的特征值數量，0 ≤m≤n。range = 'A'時，m = n；range = 'I'時，m = iu-il+1；range ='V'時，m未知（the exact value of m is not known in advance）

w：數組，至少max(1, n)。，包含升序的特征值，存儲在w[0]w[m - 1]。

z：列優先時，大小max(1, ldz*m)；行優先時，大小max(1, ldz*n)。jobz = 'V'時，如果返回值為0，z的前m列包含對應於特征值的正交的特征向量（with the i-th column of z holding the eigenvector associated with w[i - 1]）；jobz = 'N'時，不使用z。

isuppz：數組，最小 2 *max(1, m)。The support of the eigenvectors in z, i.e., the indices indicating the nonzero elements in z. The i-th eigenvector is nonzero only in elements isuppz[2i - 2] through isuppz[2i - 1]. Referenced only if eigenvectors are needed (jobz = 'V') and all eigenvalues are needed, that is, range = 'A' or range = 'I' and il = 1 and iu = n.

返回值：

0：計算成功

-i：第i個參數不合法

i：發生了內部錯誤。

 

程序：

下面的程序將參數固定了，實際上可自行調整。

// 計算矩陣的特征值和特征向量
int SEigen(float* pEigVal, float* pEigVec, const float* pSrc, int dim)
{
    float* eigValImag = new float[dim];
    float* eigVecVl = new float[dim * dim];
    float* pSrcBak = new float[dim * dim];
    memcpy(pSrcBak, pSrc, sizeof(float) * dim * dim);

    int nRetVal = LAPACKE_sgeev(LAPACK_ROW_MAJOR, 'N', 'V', dim, pSrcBak, dim,
        pEigVal, eigValImag, eigVecVl, dim, pEigVec, dim);   // 計算特征值和特征向量

    delete[] eigValImag;
    eigValImag = nullptr;
    delete[] eigVecVl;
    eigVecVl = nullptr;
    delete[] pSrcBak;
    pSrcBak = nullptr;

    return nRetVal;
}

// 計算實對稱矩陣的特征值和特征向量
int SRealSymEigen1(float* pEigVal, float* pEigVec, const float* pSrc, int dim)
{
    memcpy(pEigVec, pSrc, sizeof(float)* dim * dim);
    return LAPACKE_ssyevd(LAPACK_ROW_MAJOR, 'V', 'U', dim, pEigVec, dim, pEigVal);   // 計算特征值和特征向量
}

// 計算實對稱矩陣的特征值和特征向量
int SRealSymEigen2(float* pEigVal, float* pEigVec, const float* pSrc, int dim)
{
    memcpy(pEigVec, pSrc, sizeof(float)* dim * dim);
    return LAPACKE_ssyev(LAPACK_ROW_MAJOR, 'V', 'U', dim, pEigVec, dim, pEigVal);  // 計算特征值和特征向量
}

 

測試程序：

    const int nDim = 3;
    float pa[nDim * nDim] = { 1, 2, 3, 6, 5, 4, 7, 5, 8 };
    float eigVal[nDim];
    float eigVec[nDim * nDim];

    memset(eigVal, 0, sizeof(float)* nDim);
    memset(eigVec, 0, sizeof(float)* nDim * nDim);
    int ret1 = SEigen(eigVal, eigVec, pa, nDim);  // 計算非對稱矩陣的特征值和特征向量

    memset(eigVal, 0, sizeof(float)* nDim);
    memset(eigVec, 0, sizeof(float)* nDim * nDim);
    int ret2 = SRealSymEigen1(eigVal, eigVec, pa, nDim);  // 計算實對稱矩陣的特征值和特征向量

    int ret3 = SRealSymEigen2(eigVal, eigVec, pa, nDim);  // 計算實對稱矩陣的特征值和特征向量

 

SRealSymEigen1中LAPACKE_sgeev結果：

特征值：

特征向量：

matlab結果：

注意的是，C++中每列存儲對應的特征向量。除了正負號之外，和matlab結果一致。

 

SrealSymEigen2中LAPACKE_ssyevd結果：

特征值：

特征向量：

matlab結果：

可見結果一致。

注意：測試程序中，pa並不是對稱矩陣，但是使用SRealSymEigen1后，仍舊能得到正確的結果。說明，LAPACKE_ssyevd不會校驗輸入矩陣是否對稱，直接按照輸入要求，把對應的上三角或下三角作為A矩陣進行計算。

 

SrealSymEigen3中LAPACKE_ssyev結果：

特征值：

特征向量：

可見和LAPACKE_ssyevd特征值一樣，特征向量在誤差允許的范圍內一致。

 

對於最后一個函數LAPACKE_ssyevr：

// 計算實對稱矩陣的特征值和特征向量
int SRealSymEigen3()
{
    const int nDim = 3;
    float a[nDim * nDim] = { 1, 2, 3, 6, 5, 4, 7, 5, 8 };
    float eigVal[nDim];
    float eigVec[nDim * nDim];
    int isuppz[nDim * 2] = { 0 };
    int m = 0;

    float aBakUp[nDim * nDim];
    memcpy(aBakUp, a, sizeof(float) * nDim * nDim);

    memset(eigVal, 0, sizeof(float)* nDim);
    memset(eigVec, 0, sizeof(float)* nDim * nDim);
    LAPACKE_ssyevr(LAPACK_ROW_MAJOR, 'V', 'A', 'U', nDim, aBakUp,
        nDim, -2, 5, 2, 3, 10, &m, eigVal, eigVec, nDim, isuppz);  // step1

    memcpy(aBakUp, a, sizeof(float)* nDim * nDim);
    memset(eigVal, 0, sizeof(float)* nDim);
    memset(eigVec, 0, sizeof(float)* nDim * nDim);
    LAPACKE_ssyevr(LAPACK_ROW_MAJOR, 'V', 'V', 'U', nDim, aBakUp,
         nDim, -2, 5, 2, 3, 10, &m, eigVal, eigVec, nDim, isuppz); // step2

    memcpy(aBakUp, a, sizeof(float)* nDim * nDim);
    memset(eigVal, 0, sizeof(float)* nDim);
    memset(eigVec, 0, sizeof(float)* nDim * nDim);
    LAPACKE_ssyevr(LAPACK_ROW_MAJOR, 'V', 'I', 'U', nDim, aBakUp,
         nDim, -2, 5, 2, 3, 10, &m, eigVal, eigVec, nDim, isuppz); // step3

     return 0;
}

執行完step1后：

特征值：

特征向量：

此時，所有特征值和特征向量都計算。

 

執行完step2后：

特征值：

特征向量：

此時，只篩選了符合要求的特征值（特征向量全部保存），但是特征值也不完全對應，不懂。

 

執行完step3后：

特征值：

特征向量：

此時，特征值和特征向量都不對應，更不懂了。。。