关于Intel MKL 特征值分解 dsyevr 的说明

void dsyevr( const char* jobz, const char* range, const char* uplo, 

             const MKL_INT* n, double* a, const MKL_INT* lda, const double* vl,

             const double* vu, const MKL_INT* il, const MKL_INT* iu,

             const double* abstol, MKL_INT* m, double* w, double* z,

             const MKL_INT* ldz, MKL_INT* isuppz, double* work,

             const MKL_INT* lwork, MKL_INT* iwork, const MKL_INT* liwork,

             MKL_INT* info );

 

首先调用 ?sytrd将矩阵A分解为三对角线矩阵T；

然后尽可能调用 ?stemr 使用相对鲁棒表示计算eigenspectrum。

 ?stemr 通过分而治之算法计算eigenvalues，而通过各种“好的”L*D*LT表示（也称为相对鲁棒表示）计算正交eigen向量。尽可能避免Gram-Schmidt正交化。

 

 

 

?syevr 适合于绝大多数实对称eigenvalue问题，因为其潜在算法快且需要空间更少。

 

jobz = 'V'，则计算eigenvalues 和 eigenvectors；

range = 'A', 则计算所有的eigenvalues；

uplo = 'U', a存储A的上三角部分；

n，A的秩；

lda，a的第一维维数；

vl, vu，如果range = 'A' or 'I', vl和vu没用；

il, iu，如果range = 'A' or 'V', il和iu没用；

abstol，每个所需要的eigenvalues 和 eigenvectors的绝对误差上限（？）；输出的期望精度可以通过输入参数abstol来确定。如果相对高的精度很重要，则设置abstol为 ?lamch('S')。

ldz，输出z的维数；

lwork，work的维数；

iwork，Workspace数组；

liwork，iwork的维数

w存放eigenvalues 

z存放eigenvectors

 

 

#include <stdio.h>

#include <stdlib.h>

#include <string.h>

#include <assert.h>

#include <math.h>

#include "mkl_blas.h"

#include "mkl_service.h"

#include "mkl_lapack.h"

 

#define N 3

static double A[9] = {2, -1, 0, -1, 2, -1, 0, -1, 2};

 

int main (int argc, char **argv)

{

int n = N;

char jobz = 'V', range = 'A', uplo = 'U';

int i, lda, il, iu, m, ldz, *isuppz;

double *a, *w, *z, *work, vl, vu, abstol, dtemp;

int lwork, liwork, *iwork, info, itemp;

 

lda = n;

abstol = dlamch_("S");

ldz = n;

lwork = -1;

liwork = -1;

 

a = (double*) malloc(sizeof(double)*n*n);

memcpy(a, A, sizeof(double)*n*n);

 

w = (double*) malloc(sizeof(double)*n);

z = (double*) malloc(sizeof(double)*ldz*n);

isuppz = (int*) malloc(sizeof(int)*2*n);

 

dsyevr(&jobz, &range, &uplo, &n, a, &lda, &vl, &vu, &il, &iu, &abstol, &m, w, z,

&ldz, isuppz, &dtemp, &lwork, &itemp, &liwork, &info);

assert (info == 0);

lwork = (int)dtemp;

liwork = itemp;

work = (double*) malloc(sizeof(double)*lwork);

iwork = (int*) malloc(sizeof(int)*liwork);

 

dsyevr(&jobz, &range, &uplo, &n, a, &lda, &vl, &vu, &il, &iu, &abstol, &m, w, z,

&ldz, isuppz, work, &lwork, iwork, &liwork, &info);

printf("info = %d/n", info);

printf("# eig found = %d/n/n", m);

for (i = 0; i < m; i++) printf("eigen value [%d] = %f/n", i, w[i]);

free (a);

free (w);

free (z);

free (isuppz);

free (work);

free (iwork);

 

 

{

double x = 1.0, y = 2.0, r;

printf("/n");

r = fmod(x, y);

printf("fmod(%f, %f) = %f/n/n", x, y, r); //

}

return 0;

}

 

////OUTPUT

////info = 0

////# eig found = 3

////

////eigen value [0] = 0.585786

////eigen value [1] = 2.000000

////eigen value [2] = 3.414214

////

////fmod(1.000000, 2.000000) = 1.000000

////