I would like to resolve a Ax = b linear system.
My matrix A is tridiagonal. I have been asked to use the UMFPACK package from SuiteSparse. There is a function called umfpack_di_symbolic which uses the CSR matrices:
umfpack_di_symbolic(n, n, ia, ja, a, &Symbolic, Control, Info)
So, I need to be able to convert my tridiagonal matrix into a CSR one.
I tried to write a code in C, here it is:
(NB: the original matrix mtr
is a malloc
array defined with mtr = malloc(n*n*sizeof(double));
, the matrix is filled with the FORTRAN convention (from the top to the bottom, then next column) and not with the C one (from the left to the right, then next line). The function transpo
converts the matrix from one convention to another. I am thus using the C convention all along and then recovering the FORTRAN one.)
int csr(double *mtr, int n, int **ia, int **ja, double **a) {
if (mtr == NULL || n<2) {
printf("ERROR\n");
return 1;
}
int i, j, nnz = 0;
for (j=0;j<n;j++) {
for (i=0;i<n;i++) {
if (mtr[i+j*n] != 0) nnz++;
}
}
*ia = malloc( n * sizeof(int));
*ja = malloc(nnz * sizeof(int));
*a = malloc(nnz * sizeof(double));
if (ia == NULL || ja == NULL || a == NULL) {
printf("ERROR\n");
return 1;
}
if (transpo(mtr, n)) {
printf("failed\n");
return 1;
}
int xia = 0, xja = 0, xa = 0;
int need_ia;
for (j=0;j<n;j++) {
need_ia = 1;
for (i=0;i<n;i++) {
double aij = mtr[i + n*j];
if (aij != 0) {
if (xa >= nnz || xja >= nnz) {
printf("ERROR\n", nnz);
return 1;
}
(*a)[xa++] = aij;
if (need_ia) {
if (xia >= n) {
printf("ERROR\n");
return 1;
}
(*ia)[xia++] = xja;
need_ia = 0;
}
(*ja)[xja++] = i;
}
}
}
FILE *fa = NULL;
fa = fopen("CSRa.txt", "w");
if (fa == NULL) {
printf("ERROR\n");
return 1;
}
for (i=0;i<(nnz);i++) {
fprintf(fa, "%f\n", (*a)[i]);
}
fclose(fa);
FILE *fja = NULL;
fja = fopen("CSRja.txt", "w");
if (fja == NULL) {
printf("ERROR\n");
return 1;
}
for (i=0;i<(nnz);i++) {
fprintf(fja, "%d\n", (*ja)[i]);
}
fclose(fja);
FILE *fia = NULL;
fia = fopen("CSRia.txt", "w");
if (fia == NULL) {
printf("ERROR\n");
return 1;
}
for (i=0;i<(n);i++) {
fprintf(fia, "%d\n", (*ia)[i]);
}
fclose(fia);
if (transpo(mtr, n)) {
printf("failed\n");
return 1;
}
return 0;
}
This code gives me the 3 arrays ia
, ja
and a
(a
contains the value, ja
the column number of each number and ia
shows the first element of each line refered in ja
).
I get:
UMFPACK V5.7.1 (Oct 10, 2014): ERROR: input matrix is invalid
Please help me converting my A matrix to the right CSR model !