- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
Link Copied
- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page