[petsc-dev] TS confusion
Barry Smith
bsmith at mcs.anl.gov
Sat Sep 14 21:55:55 CDT 2013
If no ijacobian is provided and not imex TS still seems to produce two matrices
Mat Arhs = NULL,Brhs = NULL;
MatStructure flg2;
if (rhsjacobian) {
ierr = TSGetRHSMats_Private(ts,&Arhs,&Brhs);CHKERRQ(ierr);
ierr = TSComputeRHSJacobian(ts,t,U,&Arhs,&Brhs,&flg2);CHKERRQ(ierr);
}
while TSComputeRHSJacobian
ierr = TSGetIJacobian(ts,&A,&B,NULL,NULL);CHKERRQ(ierr);
if (Arhs) {
if (!ts->Arhs) {
ierr = MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&ts->Arhs);CHKERRQ(ierr);
}
*Arhs = ts->Arhs;
}
} 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);
The end result is two sets of matrices. A and Arhs and then MatAXPY costs me
MatAXPY 47 1.0 7.1141e+00 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 22 0 0 0 0 22 0 0 0 0 0
when it doesn't even need to be done. The code needs some refactoring but I don't understand enough about its convoluted logic.
What is going on here. Why does it feel it needs to create the non RHS Jacobian matrix?
Barry
More information about the petsc-dev
mailing list