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



Reply via email to