Hello, when I look at IJacobian function in the example:
http://www.mcs.anl.gov/petsc/petsc-current/src/ts/examples/tutorials/ex19.c
static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal
a,Mat *A,Mat *B,MatStructure *flag,void *ctx)
{
PetscErrorCode ierr;
PetscInt rowcol[] = {0,1};
PetscScalar *x,J[2][2];
PetscFunctionBeginUser;
ierr = VecGetArray(X,&x);CHKERRQ(ierr);
J[0][0] = a; J[0][1] = -1.;
J[1][0] = 1.; J[1][1] = -1. + x[1]*x[1];
ierr =
MatSetValues(*B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);CHKERRQ(ierr);
ierr = VecRestoreArray(X,&x);CHKERRQ(ierr);
ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
if (*A != *B) {
ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
}
*flag = SAME_NONZERO_PATTERN;
PetscFunctionReturn(0);
}
There's a PetscReal a. Is it a constant? And if the value of a known outside of
the IJacobian function, saying the main() part.
If it is, then as shown in this example, J[0][0] = a; J[0][1] = -1.; J[1][0] =
1.; are constant through all the iterations. How can I use
MatRetrieveValues(*A) and MatStoreValues(*A) to reuse them?
I'm asking this question because I have a large Jacobian matrix, but half of
the matrix contains constant values. However, I don't know what's "PetscReal a"
when I try to use it to compute my constant elements of Jacobian matrix which
doesn't depend on the value of x at all.
Thanks,
Shuangshuang