Relatively Robust Representations algorithmn (dstevr)

y_sigolaev — Thu, 19 Jul 2007 19:59:06 GMT

The big efforts have been spent to repair Relatively Robust Representations algorithmn (dstevr),
but also now it is easy to find examples in which bad results turn out
(see also section 2 on the page devoted to P4 processor: http://www.thesa-store.com/products):

//icl /O2 mkl_dstevr_test.c mkl_c.lib

#include
#include
#include

#include

int main() {
int n = 1001;
double *d__, *e;
double vl, vu;
int il, iu;
double abstol;
int m;
double *w;
double *z__;
int ldz;
int *isuppz;
double *work;
int lwork;
int *iwork;
int liwork;
int info;

int i, j;
double sum, max_modulus;

d__= (double *) malloc (n * sizeof (double));
e = (double *) malloc (n * sizeof (double));
w = (double *) malloc (n * sizeof (double));
z__ = (double *) malloc (n * n * sizeof (double));
isuppz = (int *) malloc (2 * n * sizeof (int));
lwork = 20 * n;
work = (double *) malloc (lwork * sizeof (double));
liwork = 10 * n;
iwork = (int *) malloc (liwork * sizeof (double));
if (! iwork || ! work || ! isuppz || ! z__ || ! e || ! d__ || ! w) {
printf("Not enough memory to allocate buffer ");
exit(1);
}
for (i = 0; i < n; i++) {
if (i <= n / 2) {
d__ = sin((double)(n / 2 - i));
} else {
d__ = sin((double)(i - n / 2));
}
if (i != n - 1) {
e = 1.0;
}
}
abstol = dlamch("S");
ldz = n;
dstevr("V",
"A",
&n,
d__,
e,
&vl,
&vu,
&il,
&iu,
&abstol,
&m,
w,
z__,
&ldz,
isuppz,
work,
&lwork,
iwork,
&liwork,
&info);
if (! info) {
max_modulus = 0.0;
for (i = 0; i < m - 1; i++) {
sum = 0.0;
for (j = 0; j < n; j++) {
sum += z__[i * n + j] * z__[(i + 1) * n + j];
}
if (fabs (sum) > max_modulus) {
max_modulus = fabs (sum);
}
}
printf("Max modulus of scalar dot of contiguous vectors = %e ", max_modulus);
} else {
printf("info_dstevr=%d ", info);
}
free (d__);
free (e);
free (w);
free (z__);
free (isuppz);
free (work);
free (iwork);
return 0;
}

//Max modulus of scalar dot of contiguous vectors = 6.109680e-001

topic Relatively Robust Representations algorithmn (dstevr) in Intel® oneAPI Math Kernel Library

Relatively Robust Representations algorithmn (dstevr)