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