[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