Relatively Robust Representations algorithmn (dstevr)

[y_sigolaev]

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