C++ MatZeroEntries函数代码示例

static PetscErrorCode DMCreateMatrix_Shell(DM dm,Mat *J)
{
    PetscErrorCode ierr;
    DM_Shell       *shell = (DM_Shell*)dm->data;
    Mat            A;

    PetscFunctionBegin;
    PetscValidHeaderSpecific(dm,DM_CLASSID,1);
    PetscValidPointer(J,3);
    if (!shell->A) {
        if (shell->Xglobal) {
            PetscInt m,M;
            ierr = PetscInfo(dm,"Naively creating matrix using global vector distribution without preallocation\n");
            CHKERRQ(ierr);
            ierr = VecGetSize(shell->Xglobal,&M);
            CHKERRQ(ierr);
            ierr = VecGetLocalSize(shell->Xglobal,&m);
            CHKERRQ(ierr);
            ierr = MatCreate(PetscObjectComm((PetscObject)dm),&shell->A);
            CHKERRQ(ierr);
            ierr = MatSetSizes(shell->A,m,m,M,M);
            CHKERRQ(ierr);
            ierr = MatSetType(shell->A,dm->mattype);
            CHKERRQ(ierr);
            ierr = MatSetUp(shell->A);
            CHKERRQ(ierr);
        } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Must call DMShellSetMatrix(), DMShellSetCreateMatrix(), or provide a vector");
    }
    A = shell->A;
    /* the check below is tacky and incomplete */
    if (dm->mattype) {
        PetscBool flg,aij,seqaij,mpiaij;
        ierr = PetscObjectTypeCompare((PetscObject)A,dm->mattype,&flg);
        CHKERRQ(ierr);
        ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&seqaij);
        CHKERRQ(ierr);
        ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&mpiaij);
        CHKERRQ(ierr);
        ierr = PetscStrcmp(dm->mattype,MATAIJ,&aij);
        CHKERRQ(ierr);
        if (!flg) {
            if (!(aij && (seqaij || mpiaij))) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_NOTSAMETYPE,"Requested matrix of type %s, but only %s available",dm->mattype,((PetscObject)A)->type_name);
        }
    }
    if (((PetscObject)A)->refct < 2) { /* We have an exclusive reference so we can give it out */
        ierr = PetscObjectReference((PetscObject)A);
        CHKERRQ(ierr);
        ierr = MatZeroEntries(A);
        CHKERRQ(ierr);
        *J   = A;
    } else {                      /* Need to create a copy, could use MAT_SHARE_NONZERO_PATTERN in most cases */
        ierr = MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,J);
        CHKERRQ(ierr);
        ierr = MatZeroEntries(*J);
        CHKERRQ(ierr);
    }
    PetscFunctionReturn(0);
}

void PETScLinearSolver::Initialize( )
{

   VecSet(b, 0.0);
   VecSet(x, 0.0);
   MatZeroEntries(A);
} 

void PetscSparseStorage::atPutZeros( int row, int col,
				     int rowExtent, int colExtent )
{
  int m, n;
  this->getSize( m, n );
  if( 0 == row && 0 == col &&
      m == rowExtent && n == colExtent ) {
    int ierr;
    ierr = MatZeroEntries( M ); assert( ierr  == 0);
  } else {
    double * zeros = new double[colExtent];
    double * A     = new double[colExtent];
    int    * jcol  = new int   [colExtent];
	
    int i, nnz, info;
    for ( i = 0; i < colExtent; i++ ) {
      zeros[i] = 0.0;
    }

    for( i = row; i < row + rowExtent; i++ ) {
      // Both calls will always succeed.
      this->fromGetSpRow( i, col, A, colExtent, jcol, nnz,
			  colExtent, info );
      this->atPutSpRow( i, zeros, nnz, jcol, info );
    }

    delete [] jcol;
    delete [] A;
    delete [] zeros;
  }
}

static PetscErrorCode RHSJacobianP(TS ts,PetscReal t,Vec X,Mat A,void *ctx0)
{
  PetscErrorCode ierr;
  PetscScalar    a;
  PetscInt       row,col;
  Userctx        *ctx=(Userctx*)ctx0;
  PetscInt       idx=0;
  Vec            Xgen,Xnet;
  PetscScalar    *xgen,*xnet;

  PetscScalar    Eqp,Edp;
  PetscScalar    Id,Iq;

  PetscFunctionBeginUser;

  /*if (ctx->jacp_flg) { delete me */
  ierr = MatZeroEntries(A);CHKERRQ(ierr);

    //recompute since this changes
  M[0] = 2*H[0]/w_s; M[1] = 2*H[1]/w_s; M[2] = 2*H[2]/w_s;
  D[0] = 0.1*M[0]; D[1] = 0.1*M[1]; D[2] = 0.1*M[2];

  ierr = DMCompositeGetLocalVectors(ctx->dmpgrid,&Xgen,&Xnet);CHKERRQ(ierr);
  ierr = DMCompositeScatter(ctx->dmpgrid,X,Xgen,Xnet);CHKERRQ(ierr);

  /* Generator subsystem initialization */
  ierr = VecGetArray(Xgen,&xgen);CHKERRQ(ierr);
  ierr = VecGetArray(Xnet,&xnet);CHKERRQ(ierr);

  for (col=0;col<ngen;col++) {

    Eqp   = xgen[idx];
    Edp   = xgen[idx+1];
    Id    = xgen[idx+4];
    Iq    = xgen[idx+5];
    TM[col] = PG[col];

    a = -1.0 *(TM[col] - Edp*Id - Eqp*Iq - (Xqp[col]-Xdp[col])*Id*Iq)/M[col]/H[col];

    row  = 9*col+3; // E is in the 4th equation hence +3
    idx = idx + 9;  //9 equations per generator, hence move to next gen

    ierr = MatSetValues(A,1,&row,1,&col,&a,INSERT_VALUES);CHKERRQ(ierr);
 }

  /* ctx->jacp_flg = PETSC_FALSE; delete me */
  ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

  ierr = VecRestoreArray(Xgen,&xgen);CHKERRQ(ierr);
  ierr = VecRestoreArray(Xnet,&xnet);CHKERRQ(ierr);

  ierr = DMCompositeGather(ctx->dmpgrid,X,INSERT_VALUES,Xgen,Xnet);CHKERRQ(ierr);
  ierr = DMCompositeRestoreLocalVectors(ctx->dmpgrid,&Xgen,&Xnet);CHKERRQ(ierr);

  //MatView(A,PETSC_VIEWER_STDOUT_SELF);

  /* }delete me */
  PetscFunctionReturn(0);
}

/*
   FormHessian - Evaluates Hessian matrix.

   Input Parameters:
.  tao   - the Tao context
.  x     - input vector
.  ptr   - optional user-defined context, as set by TaoSetHessian()

   Output Parameters:
.  H     - Hessian matrix

   Note:  Providing the Hessian may not be necessary.  Only some solvers
   require this matrix.
*/
PetscErrorCode FormHessian(Tao tao,Vec X,Mat H, Mat Hpre, void *ptr)
{
  AppCtx         *user = (AppCtx*)ptr;
  PetscErrorCode ierr;
  PetscInt       i, ind[2];
  PetscReal      alpha=user->alpha;
  PetscReal      v[2][2],*x;
  PetscBool      assembled;

  /* Zero existing matrix entries */
  ierr = MatAssembled(H,&assembled);CHKERRQ(ierr);
  if (assembled){ierr = MatZeroEntries(H); CHKERRQ(ierr);}

  /* Get a pointer to vector data */
  ierr = VecGetArray(X,&x);CHKERRQ(ierr);

  /* Compute H(X) entries */
  for (i=0; i<user->n/2; i++){
    v[1][1] = 2*alpha;
    v[0][0] = -4*alpha*(x[2*i+1]-3*x[2*i]*x[2*i]) + 2;
    v[1][0] = v[0][1] = -4.0*alpha*x[2*i];
    ind[0]=2*i; ind[1]=2*i+1;
    ierr = MatSetValues(H,2,ind,2,ind,v[0],INSERT_VALUES);CHKERRQ(ierr);
  }
  ierr = VecRestoreArray(X,&x);CHKERRQ(ierr);

  /* Assemble matrix */
  ierr = MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = PetscLogFlops(9.0*user->n/2.0);CHKERRQ(ierr);
  return 0;
}

double SFieldSolveFor(SField sfv, double *Y, unsigned int yCount) {
    mySField sf = static_cast<mySField>(sfv);
    assert(yCount <= sf->maxN);
    assert(Y);
    assert(sf->running);
    sf->Y = Y;
    sf->curN = yCount;

    // -------------- SOLVE
    PetscErrorCode ierr;
    PetscLogDouble tic,toc;
    PetscTime(&tic);
    int pt[sf->d];

    ierr = MatZeroEntries(sf->J); CHKERRQ(ierr);
    JacobianOnD(sf->J, sf->F, 0, pt, sf);

    ierr = MatAssemblyBegin(sf->J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    ierr = MatAssemblyEnd(sf->J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);

    PetscTime(&toc);
    sf->timeAssembly += toc-tic;
    PetscTime(&tic);
    ierr = VecZeroEntries(sf->U);  CHKERRQ(ierr);
    ierr = KSPSetOperators(sf->ksp, sf->J, sf->J); CHKERRQ(ierr);
    ierr = KSPSetUp(sf->ksp);  CHKERRQ(ierr);
    ierr = KSPSolve(sf->ksp,sf->F,sf->U); CHKERRQ(ierr);
    PetscTime(&toc);
    sf->timeSolver  += toc-tic;
    return Integrate(sf->U,pt,0,sf);
}

/*@
    MatFDColoringApply - Given a matrix for which a MatFDColoring context
    has been created, computes the Jacobian for a function via finite differences.

    Collective on MatFDColoring

    Input Parameters:
+   mat - location to store Jacobian
.   coloring - coloring context created with MatFDColoringCreate()
.   x1 - location at which Jacobian is to be computed
-   sctx - context required by function, if this is being used with the SNES solver then it is SNES object, otherwise it is null

    Options Database Keys:
+    -mat_fd_type - "wp" or "ds"  (see MATMFFD_WP or MATMFFD_DS)
.    -mat_fd_coloring_view - Activates basic viewing or coloring
.    -mat_fd_coloring_view draw - Activates drawing of coloring
-    -mat_fd_coloring_view ::ascii_info - Activates viewing of coloring info

    Level: intermediate

.seealso: MatFDColoringCreate(), MatFDColoringDestroy(), MatFDColoringView(), MatFDColoringSetFunction()

.keywords: coloring, Jacobian, finite differences
@*/
PetscErrorCode  MatFDColoringApply(Mat J,MatFDColoring coloring,Vec x1,void *sctx)
{
  PetscErrorCode ierr;
  PetscBool      flg = PETSC_FALSE;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(J,MAT_CLASSID,1);
  PetscValidHeaderSpecific(coloring,MAT_FDCOLORING_CLASSID,2);
  PetscValidHeaderSpecific(x1,VEC_CLASSID,3);
  if (!coloring->f) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_WRONGSTATE,"Must call MatFDColoringSetFunction()");
  if (!J->ops->fdcoloringapply) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not supported for this matrix type %s",((PetscObject)J)->type_name);
  if (!coloring->setupcalled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call MatFDColoringSetUp()");

  ierr = MatSetUnfactored(J);CHKERRQ(ierr);
  ierr = PetscOptionsGetBool(NULL,"-mat_fd_coloring_dont_rezero",&flg,NULL);CHKERRQ(ierr);
  if (flg) {
    ierr = PetscInfo(coloring,"Not calling MatZeroEntries()\n");CHKERRQ(ierr);
  } else {
    PetscBool assembled;
    ierr = MatAssembled(J,&assembled);CHKERRQ(ierr);
    if (assembled) {
      ierr = MatZeroEntries(J);CHKERRQ(ierr);
    }
  }

  ierr = PetscLogEventBegin(MAT_FDColoringApply,coloring,J,x1,0);CHKERRQ(ierr);
  ierr = (*J->ops->fdcoloringapply)(J,coloring,x1,sctx);CHKERRQ(ierr);
  ierr = PetscLogEventEnd(MAT_FDColoringApply,coloring,J,x1,0);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PetscErrorCode _RHS_time_dep_ham_p(TS ts,PetscReal t,Vec X,Mat AA,Mat BB,void *ctx){
  double time_dep_val;
  PetscScalar time_dep_scalar;
  int i,j;
  operator op;

  MatZeroEntries(AA);

  MatCopy(full_A,AA,SAME_NONZERO_PATTERN);

  for (i=0;i<_num_time_dep;i++){
    time_dep_val = _time_dep_list[i].time_dep_func(t);
    _add_ops_to_mat_ham(time_dep_val,AA,_time_dep_list[i].num_ops,_time_dep_list[i].ops);
  }

  for (i=0;i<_num_time_dep_lin;i++){
    time_dep_val = _time_dep_list_lin[i].time_dep_func(t);
    _add_ops_to_mat_lin(time_dep_val,AA,_time_dep_list_lin[i].num_ops,_time_dep_list_lin[i].ops);
  }

  MatAssemblyBegin(AA,MAT_FINAL_ASSEMBLY);
  MatAssemblyEnd(AA,MAT_FINAL_ASSEMBLY);

  if(AA!=BB) {
    MatAssemblyBegin(AA,MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(AA,MAT_FINAL_ASSEMBLY);
  }

  PetscFunctionReturn(0);
}

/*
   FormHessian - Evaluates Hessian matrix.

   Input Parameters:
.  tao   - the Tao context
.  x     - input vector
.  ptr   - optional user-defined context, as set by TaoSetHessian()

   Output Parameters:
.  H     - Hessian matrix

   Note:  Providing the Hessian may not be necessary.  Only some solvers
   require this matrix.
*/
PetscErrorCode FormHessian(Tao tao,Vec X,Mat H, Mat Hpre, void *ptr)
{
  AppCtx            *user = (AppCtx*)ptr;
  PetscErrorCode    ierr;
  PetscInt          i, ind[2];
  PetscReal         alpha=user->alpha;
  PetscReal         v[2][2];
  const PetscScalar *x;
  PetscBool         assembled;

  PetscFunctionBeginUser;
  /* Zero existing matrix entries */
  ierr = MatAssembled(H,&assembled);CHKERRQ(ierr);
  if (assembled){ierr = MatZeroEntries(H); CHKERRQ(ierr);}

  /* Get a pointer to vector data */
  ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr);

  /* Compute H(X) entries */
  if (user->chained) {
    ierr = MatZeroEntries(H);CHKERRQ(ierr);
    for (i=0; i<user->n-1; i++) {
      PetscScalar t1 = x[i+1] - x[i]*x[i];
      v[0][0] = 2 + 2*alpha*(t1*(-2) - 2*x[i]);
      v[0][1] = 2*alpha*(-2*x[i]);
      v[1][0] = 2*alpha*(-2*x[i]);
      v[1][1] = 2*alpha*t1;
      ind[0] = i; ind[1] = i+1;
      ierr = MatSetValues(H,2,ind,2,ind,v[0],ADD_VALUES);CHKERRQ(ierr);
    }
  } else {
    for (i=0; i<user->n/2; i++){
      v[1][1] = 2*alpha;
      v[0][0] = -4*alpha*(x[2*i+1]-3*x[2*i]*x[2*i]) + 2;
      v[1][0] = v[0][1] = -4.0*alpha*x[2*i];
      ind[0]=2*i; ind[1]=2*i+1;
      ierr = MatSetValues(H,2,ind,2,ind,v[0],INSERT_VALUES);CHKERRQ(ierr);
    }
  }
  ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);

  /* Assemble matrix */
  ierr = MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = PetscLogFlops(9.0*user->n/2.0);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PetscErrorCode MatZeroEntries_IS(Mat A)
{
  Mat_IS         *a = (Mat_IS*)A->data;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  ierr = MatZeroEntries(a->A);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

static PetscErrorCode DRDPJacobianTranspose(TS ts,PetscReal t,Vec U,Mat DRDP,AppCtx *ctx)
{
  PetscErrorCode    ierr;

  PetscFunctionBegin;
  ierr = MatZeroEntries(DRDP);CHKERRQ(ierr);
  ierr = MatAssemblyBegin(DRDP,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(DRDP,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

static PetscErrorCode FormRHSJacobian(TS ts,PetscReal t,Vec X,Mat Amat,Mat Pmat,void *ptr)
{
  User              user = (User)ptr;
  PetscErrorCode    ierr;
  const PetscScalar **x;
  PetscInt          M = user->Nspec+1,i,j,xs,xm;;
  DM                dm;

  PetscFunctionBeginUser;
  if (user->reactions) {
    ierr = TSGetDM(ts,&dm);CHKERRQ(ierr);
    ierr = MatZeroEntries(Pmat);CHKERRQ(ierr);
    ierr = MatSetOption(Pmat,MAT_ROW_ORIENTED,PETSC_FALSE);CHKERRQ(ierr);
    ierr = MatSetOption(Pmat,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);CHKERRQ(ierr);
    ierr = DMDAVecGetArrayDOFRead(dm,X,&x);CHKERRQ(ierr);
    ierr = DMDAGetCorners(dm,&xs,NULL,NULL,&xm,NULL,NULL);CHKERRQ(ierr);

    for (i=xs; i<xs+xm; i++) {
      ierr = PetscMemcpy(user->tchemwork,x[i],(user->Nspec+1)*sizeof(x[xs][0]));CHKERRQ(ierr);
      user->tchemwork[0] *= user->Tini;  /* Dimensionalize temperature (first row) because that is what Tchem wants */
      ierr = TC_getJacTYN(user->tchemwork,user->Nspec,user->Jdense,1);CHKERRQ(ierr);

      for (j=0; j<M; j++) user->Jdense[j + 0*M] /= user->Tini; /* Non-dimensionalize first column */
      for (j=0; j<M; j++) user->Jdense[0 + j*M] /= user->Tini; /* Non-dimensionalize first row */
      for (j=0; j<M; j++) user->rows[j] = i*M+j;
      ierr = MatSetValues(Pmat,M,user->rows,M,user->rows,user->Jdense,INSERT_VALUES);CHKERRQ(ierr);
    }
    ierr = DMDAVecRestoreArrayDOFRead(dm,X,&x);CHKERRQ(ierr);
    ierr = MatAssemblyBegin(Pmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
    ierr = MatAssemblyEnd(Pmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  } else {
    ierr = MatZeroEntries(Pmat);CHKERRQ(ierr);
  }
  if (user->diffusion) {
    ierr = FormDiffusionJacobian(ts,t,X,Amat,Pmat,ptr);CHKERRQ(ierr);
  }
  if (Amat != Pmat) {
    ierr = MatAssemblyBegin(Amat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
    ierr = MatAssemblyEnd(Amat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

/*
   FormJacobian1 - Evaluates Jacobian matrix.

   Input Parameters:
.  snes - the SNES context
.  x - input vector
.  dummy - optional user-defined context (not used here)

   Output Parameters:
.  jac - Jacobian matrix
.  B - optionally different preconditioning matrix
.  flag - flag indicating matrix structure
*/
PetscErrorCode FormJacobian1(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *ictx)
{
  PetscScalar    *xx;
  PetscErrorCode ierr;
  PetscInt       i;
  Ctx            *ctx = (Ctx*)ictx;

  ierr = MatZeroEntries(*B);CHKERRQ(ierr);
  /*
     Get pointer to vector data
  */
  ierr = VecGetArray(x,&xx);CHKERRQ(ierr);

  /*
     Compute Jacobian entries and insert into matrix.
      - Since this is such a small problem, we set all entries for
        the matrix at once.
  */
  ierr = MatSetValue(*B,0,0, 2.0 + 1200.0*xx[0]*xx[0] - 400.0*xx[1],ADD_VALUES);CHKERRQ(ierr);
  ierr = MatSetValue(*B,0,1,-400.0*xx[0],ADD_VALUES);CHKERRQ(ierr);

  for (i=1; i<ctx->p+1; i++) {
    ierr = MatSetValue(*B,i,i-1, -400.0*xx[i-1],ADD_VALUES);CHKERRQ(ierr);
    ierr = MatSetValue(*B,i,i, 2.0 + 1200.0*xx[i]*xx[i] - 400.0*xx[i+1] + 200.0,ADD_VALUES);CHKERRQ(ierr);
    ierr = MatSetValue(*B,i,i+1,-400.0*xx[i],ADD_VALUES);CHKERRQ(ierr);
  }

  ierr = MatSetValue(*B,ctx->p+1,ctx->p, -400.0*xx[ctx->p],ADD_VALUES);CHKERRQ(ierr);
  ierr = MatSetValue(*B,ctx->p+1,ctx->p+1,200,ADD_VALUES);CHKERRQ(ierr);

  *flag = SAME_NONZERO_PATTERN;

  for (i=ctx->p+2; i<2+ctx->p+ctx->n; i++) {
    ierr = MatSetValue(*B,i,i,1.0,ADD_VALUES);CHKERRQ(ierr);
    ierr = MatSetValue(*B,i,0,-1.0,ADD_VALUES);CHKERRQ(ierr);
    ierr = MatSetValue(*B,i,1,.7,ADD_VALUES);CHKERRQ(ierr);
    ierr = MatSetValue(*B,i,i-1,-.4*xx[i-1],ADD_VALUES);CHKERRQ(ierr);
  }
  /*
     Restore vector
  */
  ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr);

  /*
     Assemble matrix
  */
  ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  if (*jac != *B) {
    ierr = MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
    ierr = MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  }
  return 0;
}

void
PetscSparseMtrx :: zero()
{
    // test if receiver is already assembled
    PetscBool assembled;
    MatAssembled(this->mtrx, & assembled);
    if ( assembled ) {
        MatZeroEntries(this->mtrx);
    }
    this->newValues = true;
}

PetscErrorCode SNESComputeJacobian_MyShell(SNES snes,Vec X,Mat A,Mat B,void *ctx)
{
  static PetscInt fail = 0;
  PetscErrorCode  ierr;

  PetscFunctionBegin;
  ierr = SNESComputeJacobian_DMDA(snes,X,A,B,ctx);CHKERRQ(ierr);
  if (fail++ > 0) {
    ierr = MatZeroEntries(A);CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}