I'm trying to use mkl to compute a equation. But it seem that the array a[] leaks all the time lik this:
44.62 -0.09 -6277438562204192487878988888393020692503707483087375482269988814848.00 -6277438562204192487878988888393020692503707483087375482269988814848.00 -6277438562204192487878988888393020692503707483087375482269988814848.00
-0.09 11.29 -0.09 -6277438562204192487878988888393020692503707483087375482269988814848.00 -6277438562204192487878988888393020692503707483087375482269988814848.00
-6277438562204192487878988888393020692503707483087375482269988814848.00 -0.09 0.18 -0.09 -6277438562204192487878988888393020692503707483087375482269988814848.00
-6277438562204192487878988888393020692503707483087375482269988814848.00 -6277438562204192487878988888393020692503707483087375482269988814848.00 -0.09 11.29 -0.09
-6277438562204192487878988888393020692503707483087375482269988814848.00 -6277438562204192487878988888393020692503707483087375482269988814848.00 -6277438562204192487878988888393020692503707483087375482269988814848.00 -0.09 44.62
And my code is:
#include <stdlib.h>
#include <stdio.h>
#include "mkl_lapacke.h"
/* Auxiliary routines prototypes */
extern void my_print_matrix(char* desc, MKL_INT m, MKL_INT n, double* a, MKL_INT lda, FILE* fpWrite);
extern void print_matrix(char* desc, MKL_INT m, MKL_INT n, double* a, MKL_INT lda);
/* Parameters */
#define N 5//nstep
#define LDA N
#define RMIN -10.0
#define RMAX 10.0
/* Main program */
int main() {
/* Locals */
MKL_INT n = N, lda = LDA, info;
/* Local arrays */
double h = (RMAX - RMIN) / (double(N) + 1.0);;
double xi;
double *w;
double *a;
w= (double*)malloc(sizeof(double) * N);
a = (double*)malloc(sizeof(double) * N*LDA);
for (int i = 0; i < N; i++) {
xi = RMIN + double(1.0+i) * h;
a[i*(N+1)] = 2.0 / h / h+xi * xi;
if (i==0) {
a[1] = -1.0 / h / h;
}
else if (i == N - 1) {
a[LDA * N-2] =- 1.0 / h / h;
}
else {
a[i *(N + 1)+1] = -1.0/h/h;
a[i * (N + 1) - 1] = -1.0/h/h;
}
}
print_matrix("Matrix", n, n, a, lda);
/* Executable statements */
printf("LAPACKE_dsyev (row-major, high-level) Example Program Results\n");
/* Solve eigenproblem */
info = LAPACKE_dsyev(LAPACK_ROW_MAJOR, 'V', 'U', n, a, lda, w);
/* Check for convergence */
if (info > 0) {
printf("The algorithm failed to compute eigenvalues.\n");
exit(1);
}
exit(0);
} /* End of LAPACKE_dsyev Example */
/* Auxiliary routine: printing a matrix */
void print_matrix(char* desc, MKL_INT m, MKL_INT n, double* a, MKL_INT lda) {
MKL_INT i, j;
printf("\n %s\n", desc);
for (i = 0; i < m; i++) {
for (j = 0; j < n; j++) printf(" %6.2f", a[i * lda + j]);
printf("\n");
}
}
The first two elements is seem right but wrong numbers... And the last numbers are a same but very large number. I think it's array leaking,but I don't know how to deal with it. So I ask for help.
And the a[] show is just initialzed,not the result. My problem where is initializing wrong?
Lets take a closer look at how you initialize a:
a[i*(N+1)] = 2.0 / h / h+xi * xi;
if (i==0) {
a[1] = -1.0 / h / h;
}
else if (i == N - 1) {
a[LDA * N-2] =- 1.0 / h / h;
}
else {
a[i *(N + 1)+1] = -1.0/h/h;
a[i * (N + 1) - 1] = -1.0/h/h;
}
Lets take the two cases when i == 0 and i == 1:
i == 0
Here you first do the unconditional initialization
a[i*(N+1)] = 2.0 / h / h+xi * xi;
If we calculate the index i*(N+1) it's 0*(N+1) which is 0. Therefore you will initialize a[0].
Then you have if (i==0) where you initialize a[1].
i == 1
First the unconditional initialization of index i*(N+1), which then is 1*(50+1) which equals 51. So here you initialize a[51].
Then the conditions i==1 and i == N - 1 are both false, so we end up in the final else clause:
a[i *(N + 1)+1] = -1.0/h/h;
a[i * (N + 1) - 1] = -1.0/h/h;
The first index i *(N + 1)+1 will be 1 *(50 + 1)+1 which is 52. So you initialize a[52].
The next index i * (N + 1) - 1 will be 1 * (50 + 1) - 1 which is 50. So you initialize a[50].
This pattern repeats throughout the loop, with ever higher indexes, but never lower.
That means you will never initialize index 2 to 49. These elements will have indeterminate values, and if you're unlucky one of those values could be trap-values which would lead to undefined behavior when using them.
You need to rework your algorithm to initialize all elements of the array a.