Run src/ts/examples/tutorials/advection-diffusion-reaction/ex5.c which provides
only a rhs function and rhsjacobian and requests fully implicit
falls into the chunk of code in TSComputeIJacobian() below repeatedly.
} else if (Arhs) { /* Both IJacobian and RHSJacobian */
MatStructure axpy = DIFFERENT_NONZERO_PATTERN;
if (!ijacobian) { /* No IJacobian provided, but we have a
separate RHS matrix */
ierr = MatZeroEntries(*A);CHKERRQ(ierr);
ierr = MatShift(*A,shift);CHKERRQ(ierr);
if (*A != *B) {
ierr = MatZeroEntries(*B);CHKERRQ(ierr);
ierr = MatShift(*B,shift);CHKERRQ(ierr);
}
}
ierr = MatAXPY(*A,-1,Arhs,axpy);CHKERRQ(ierr);
if (*A != *B) {
ierr = MatAXPY(*B,-1,Brhs,axpy);CHKERRQ(ierr);
}
*flg = PetscMin(*flg,flg2);c
}
Not that axpy is DIFFERENT_NONZERO_PATTERN which means a very slow MatAXPY().
If I understand things correctly after the first time through here *A will have
a nonzero pattern at least as big as Arhs? If so how can the code be organized
so that slow MatAXPY() is not used for an eternity and switches to using subset
for much faster performance?
Thanks
Barry