0

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 !

  • If the matrix is the transpose of the format that you want, you can do the symbolic operation as normal, and then do the numeric factorization using the transpose matrix option. – Juan May 10 '23 at 14:01

0 Answers0