1
$\begingroup$

Was trying to implement a poisson 2d solver using Conjugate Gradient Method, so from 10x10 grid the matrix becomes 100x100 (since we have 100 nodes to find the values at), 100x100 grid goes to 10000x10000 matrix, but matrix is sparse with 5 elements on each row (not pentadiagonal!).

Hence, receiving segmentation fault, perhaps due to bad memory allocation, but still, quite a big matrix to store and use, because zeros are still considered as doubles, which must potentially slow down the solver.

Thought of implementing data structure for sparse matrix to overcome the memory issues. Is there any other ideas on the memory management in solvers?

Watching these lectures And reading the book of the video lecturer "Numerical Methods for Partial Differential Equations: Finite Difference and Finite Volume Methods" by Sandip Mazumder.

Below is the draft C code to fill up the coefficient matrix for all nodes into dense matrix (not sparse).

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

 

int main(){
    double L = 1;
    int N = 4, M = 3;
    int kmax = M*N;
    double sum = 0;
    int fi = 1;
    double dx = L/(N-1);
    double dy = L/(M-1);
    double dx2 = dx*dx; 
    double dy2 = dy*dy;
    double *S; 
    S = (double *) malloc(N*M * sizeof(double));
    double *Q; 
    Q = (double *) malloc(N*M * sizeof(double));

    int** Acoord = (int**)malloc(kmax * sizeof(int*));
    for (int i = 0; i < kmax; i++)
        Acoord[i] = (int*)malloc(5 * sizeof(int));
    double** Aval = (double**)malloc(kmax * sizeof(double*));
    for (int i = 0; i < kmax; i++)
        Aval[i] = (double*)malloc(5 * sizeof(double));

    // fill S(ource) rhs f function
    for (int i = 0; i< N; i++){
        for (int j = 0; j < M ; j++){
            int k = j*N+i;
            S[k] = -exp(-pow(i*dx-L/2,2)-pow(j*dy-L/2,2));
            sum = sum + S[k];
        }
    }
    sum = sum/N/M;
    for (int i = 0; i < N; i++){
        for (int j = 0; j < M; j++){
            int k = j*N+i;
            S[k] = S[k] - sum;
        }
    }

    for (int i = 1; i < N-1; i++){
        for (int j = 1; j < M-1; j++){
            int k = j*N+i;

            Aval[k][2] = -(2/dx2+2/dy2);
            Aval[k][1] = 1/dx2;
            Aval[k][3] = 1/dx2;
            Aval[k][0] = 1/dy2;
            Aval[k][4] = 1/dy2;

            Acoord[k][0] = k - N;
            Acoord[k][1] = k - 1;
            Acoord[k][2] = k;
            Acoord[k][3] = k + 1;
            Acoord[k][4] = k + N;
            Q[k] = S[k];
        }
    }
    


    
    // fill boundary part
    int i = N-1; // RIGHT
    double JR = 0; // FLUX
    for (int j = 1; j < M-1; j++){ 
        int k = j*N+i;

        Aval[k][0] = 1/dy2;
        Aval[k][1] = 1/dx2;
        Aval[k][2] = -1/dx2-2/dy2;
        Aval[k][3] = 0;
        Aval[k][4] = 1/dy2;
        

        Acoord[k][0] = k - N;
        Acoord[k][1] = k - 1;
        Acoord[k][2] = k;
        Acoord[k][3] = -1;
        Acoord[k][4] = k + N;


        Q[k] = S[k];
    }

    i = 0; // LEFT
    double JL = 0; // FLUX
    for (int j = 1; j < M-1; j++){
        int k = j*N+i;

        Aval[k][0] = 1/dy2;
        Aval[k][1] = 0;
        Aval[k][2] = -1/dx2-2/dy2;
        Aval[k][3] = 1/dx2;
        Aval[k][4] = 1/dy2;

        Acoord[k][0] = k - N;
        Acoord[k][1] = - 1;
        Acoord[k][2] = k;
        Acoord[k][3] = k + 1;
        Acoord[k][4] = k + N;
        
        Q[k] = S[k];

    }

    int j = 0; // BOT
    double JB = 0; // FLUX
    for (int i = 1; i < N-1; i++){
        int k = j*N+i;

        Aval[k][0] = 0;
        Aval[k][1] = 1/dx2;
        Aval[k][2] = -2/dx2-1/dy2;
        Aval[k][3] = 1/dx2;
        Aval[k][4] = 1/dy2;

        Acoord[k][0] = -1;
        Acoord[k][1] = k - 1;
        Acoord[k][2] = k;
        Acoord[k][3] = k + 1;
        Acoord[k][4] = k + N;

        Q[k] = S[k];
        // printf("k = %d kmax = %d\n",k,kmax);
    }

    j = M-1; // TOP
    double JT = 0; // FLUX
    for (int i = 1; i < N-1; i++){
        int k = j*N+i;

        Aval[k][0] = 1/dy2;
        Aval[k][1] = 1/dx2;
        Aval[k][2] = -2/dx2-1/dy2;
        Aval[k][3] = 1/dx2;
        Aval[k][4] = 0;

        Acoord[k][0] = k - N;
        Acoord[k][1] = k - 1;
        Acoord[k][2] = k;
        Acoord[k][3] = k + 1;
        Acoord[k][4] = -1;


        Q[k] = S[k];
    }

    i = 0; // LEFT
    j = 0; // BOT 
    int k = j*N+i;

    Aval[k][0] = 0;
    Aval[k][1] = 0;
    Aval[k][2] = -1/dx2 - 1/dy2;
    Aval[k][3] = 1/dx2;
    Aval[k][4] = 1/dy2;

    Acoord[k][0] = -1;
    Acoord[k][1] = -1;
    Acoord[k][2] = k;
    Acoord[k][3] = k + 1;
    Acoord[k][4] = k + N;

    Q[k] = S[k]; 

    i = 0; // LEFT
    j = M-1; // TOP 
    k = j*N+i;

    Aval[k][0] = 1/dy2;
    Aval[k][1] = 0;
    Aval[k][2] = -1/dx2 - 1/dy2;
    Aval[k][3] = 1/dx2;
    Aval[k][4] = 0;

    Acoord[k][0] = k - N;
    Acoord[k][1] = -1;
    Acoord[k][2] = k;
    Acoord[k][3] = k + 1;
    Acoord[k][4] = -1;

    Q[k] = S[k]; 

    i = N-1; // RIGHT
    j = M-1; // TOP 
    k = j*N+i;

    Aval[k][0] = 1/dy2;
    Aval[k][1] = 1/dx2;
    Aval[k][2] = -1/dx2 - 1/dy2;
    Aval[k][3] = 0;
    Aval[k][4] = 0;

    Acoord[k][0] = k - N;
    Acoord[k][1] = k - 1;
    Acoord[k][2] = k;
    Acoord[k][3] = -1;
    Acoord[k][4] = -1;

    Q[k] = S[k]; 

    i = N-1; // RIGHT
    j = 0; // BOT 
    k = (j)*M+i;

    Aval[k][0] = 0;
    Aval[k][1] = 1/dx2;
    Aval[k][2] = -1/dx2 - 1/dy2;
    Aval[k][3] = 0;
    Aval[k][4] = 1/dy2;

    Acoord[k][0] = -1;
    Acoord[k][1] = k - 1;
    Acoord[k][2] = k;
    Acoord[k][3] = -1;
    Acoord[k][4] = k + N;

    Q[k] = S[k];     

    // CJ solver



    // CJ solver ends

    FILE *f1 = fopen("testCJ_coord.txt", "w+");
    if (f1 == NULL)
    {
        printf("Error opening file!\n");
        exit(1);
    }

    for (int i = 0; i < kmax; i++){
        fprintf(f1,"%d\t",i);
        for (int j = 0; j < 5; j++){
                fprintf(f1,"%d\t",Acoord[i][j]);
        }
        fprintf(f1,"\n");
    }
    fclose(f1);

    FILE *f2 = fopen("testCJ_val.txt", "w+");
    if (f2 == NULL)
    {
        printf("Error opening file!\n");
        exit(1);
    }

    for (int i = 0; i < kmax; i++){
        fprintf(f2,"%d\t",i);
        for (int j = 0; j < 5; j++){
                fprintf(f2,"%.4f\t",Aval[i][j]);

        }
        fprintf(f2,"\n");
    }
    fclose(f2);

    FILE *f3 = fopen("testCJ_Q.txt", "w+");
    if (f3 == NULL)
    {
        printf("Error opening file!\n");
        exit(1);
    }

    for (int i = 0; i < kmax; i++){
        fprintf(f3,"%d\t",i);
        fprintf(f3,"%12.5f",Q[i]);
        fprintf(f3,"\n");
    }
    fclose(f3);
    

    free(S);
    free(Q);
    for (int i = 0; i < kmax; i++)
        free(Aval[i]);
    free(Aval);

    for (int i = 0; i < kmax; i++)
        free(Acoord[i]);
    free(Acoord);
    
    return 0;
}
Improve this question
$\endgroup$

1 Answer 1

Reset to default
3
$\begingroup$

As you pointed out, your matrices are sparse which means that the number of non-zeros is small compared to the number of zeros. There are several formats to store such matrices, e.g. COO, CSC, CSR etc. These formats store only the non-zero entries and assume that every other entry is zero. Furthermore there are solvers and routines (sparse matrix vector product) which makes use of the formats to save calculation time and memory. I would suggest to perform research on the formats mentioned to find out which fits best for your case and then build everything else around it.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.