# 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);
}
}


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.