How to allocate memory for successive iterative solutions with potentially different non-zero structure?

Question

Background

I am solving the unsteady heat equation in 3D using an alternating direction implicit (ADI) method. This means that I am solving three different tridiagonal systems within a single timestep (once for each cartesian axis). The size of these matrices is on the order of 10^5 x 10^5.

Although each direction's A matrix has all non-zeros along the three diagonals, there are situations where one direction will have a zero value along one of the diagonals where another direction does not.

My original plan was to set up an AIJ matrix and preallocate for 3 non-zero values and reuse the same memory allocation for each of the three directions. This works fine for solving the first direction, but when I hit the second direction PETSc complains because I am attempting to add a new non-zero value where there wasn't one before (related to differences in boundary conditions).

Questions

What is the most efficient way to proceed:

• Is there a way to allocate differently that will allow me to have values change between zero and non-zero?
• OR should I create a new matrix and new allocation for each direction? (This I know works, but I'm worried about the overhead of destroying and creating a new matrix for every solution.)
• Other options?

In pseudo code, my program looks like this:

// Create and preallocate A matrix:
MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,M,M,3,NULL,3,NULL,&Amat);
    
// Create x vector
VecCreate(PETSC_COMM_WORLD,&x);

// Create b vector
VecDuplicate(x,&b);

// Solve for each timestep:
for (each timestep)
{
    for (X,Y,Z)
    {
        MatSetValues(Amat,...,INSERT_VALUES);
        VecSetValues(b,...,INSERT_VALUES);

        MatAssemblyBegin(Amat,MAT_FINAL_ASSEMBLY);
        MatAssemblyEnd(Amat,MAT_FINAL_ASSEMBLY);
        
        VecAssemblyBegin(b);
        VecAssemblyEnd(b);   
     
        KSPSetOperators(ksp, Amat, Amat, DIFFERENT_NONZERO_PATTERN);
        KSPSetUp(ksp);
        KSPSolve(ksp,b,x);

        MatZeroEntries(Amat);
    }
}

Answer 1

1. ADI is not a very good parallel algorithm. You should seriously consider formulating the problem in 3D and solving with multigrid. You could get a start with src/ksp/ksp/examples/tutorials/ex45.c.

2. If you insist on using ADI, you should seriously consider allocating the three matrices separately. It's more memory, but then you won't have to reassemble on each iteration, so you amortize setup costs. This uses more memory.

3. If you insist on ADI and on reusing the same matrix for each direction, you should allocate for three nonzeros per row and fill in explicit zeros even when the direction is missing an entry. Otherwise the first assembly compacts out the unused entries and they can't be cheaply inserted later.