C++ MatDuplicate函数代码示例

int main(int argc,char **args)
{
  Mat                 A1,A2,A3,A4,nest;
  Mat                 mata[4];
  Mat                 aij;
  MPI_Comm            comm;
  PetscInt            m,n,istart,iend,ii,i,J,j;
  PetscScalar         v;
  PetscMPIInt         size;
  PetscErrorCode      ierr;

  ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
  comm = PETSC_COMM_WORLD;
  ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
  /*
     Assemble the matrix for the five point stencil, YET AGAIN
  */
  ierr = MatCreate(comm,&A1);CHKERRQ(ierr);
  m=2,n=2;
  ierr = MatSetSizes(A1,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);CHKERRQ(ierr);
  ierr = MatSetFromOptions(A1);CHKERRQ(ierr);
  ierr = MatSetUp(A1);CHKERRQ(ierr);
  ierr = MatGetOwnershipRange(A1,&istart,&iend);CHKERRQ(ierr);
  for (ii=istart; ii<iend; ii++) {
    v = -1.0; i = ii/n; j = ii - i*n;
    if (i>0)   {J = ii - n; ierr = MatSetValues(A1,1,&ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);}
    if (i<m-1) {J = ii + n; ierr = MatSetValues(A1,1,&ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);}
    if (j>0)   {J = ii - 1; ierr = MatSetValues(A1,1,&ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);}
    if (j<n-1) {J = ii + 1; ierr = MatSetValues(A1,1,&ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);}
    v = 4.0; ierr = MatSetValues(A1,1,&ii,1,&ii,&v,INSERT_VALUES);CHKERRQ(ierr);
  }
  ierr = MatAssemblyBegin(A1,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(A1,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatView(A1,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  ierr = MatDuplicate(A1,MAT_COPY_VALUES,&A2);CHKERRQ(ierr);
  ierr = MatDuplicate(A1,MAT_COPY_VALUES,&A3);CHKERRQ(ierr);
  ierr = MatDuplicate(A1,MAT_COPY_VALUES,&A4);CHKERRQ(ierr);
  /*create a nest matrix */
  ierr = MatCreate(comm,&nest);CHKERRQ(ierr);
  ierr = MatSetType(nest,MATNEST);CHKERRQ(ierr);
  mata[0]=A1,mata[1]=A2,mata[2]=A3,mata[3]=A4;
  ierr = MatNestSetSubMats(nest,2,NULL,2,NULL,mata);CHKERRQ(ierr);
  ierr = MatSetUp(nest);CHKERRQ(ierr);
  ierr = MatConvert(nest,MATAIJ,MAT_INITIAL_MATRIX,&aij);CHKERRQ(ierr);
  ierr = MatView(aij,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  ierr = MatDestroy(&nest);CHKERRQ(ierr);
  ierr = MatDestroy(&aij);CHKERRQ(ierr);
  ierr = MatDestroy(&A1);CHKERRQ(ierr);
  ierr = MatDestroy(&A2);CHKERRQ(ierr);
  ierr = MatDestroy(&A3);CHKERRQ(ierr);
  ierr = MatDestroy(&A4);CHKERRQ(ierr);
  ierr = PetscFinalize();
  return ierr;
}

PETSC_EXTERN PetscErrorCode MatConvert_MPISBAIJ_MPISBSTRM(Mat A,MatType type,MatReuse reuse,Mat *newmat)
{
  PetscErrorCode ierr;
  Mat            B = *newmat;
  Mat_SeqSBSTRM  *sbstrm;

  PetscFunctionBegin;
  if (reuse == MAT_INITIAL_MATRIX) {
    ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr);
  }
  /* printf(" --- in MatConvert_MPISBAIJ_MPISBSTRM  -- 1 \n"); */

  ierr     = PetscNewLog(B,   Mat_SeqSBSTRM,&sbstrm);CHKERRQ(ierr);
  B->spptr = (void*)sbstrm;

  /* Set function pointers for methods that we inherit from AIJ but override.
     B->ops->duplicate   = MatDuplicate_SBSTRM;
     B->ops->mult        = MatMult_SBSTRM;
     B->ops->destroy     = MatDestroy_MPISBSTRM;
  */
  B->ops->assemblyend = MatAssemblyEnd_MPISBSTRM;

  /* If A has already been assembled, compute the permutation. */
  if (A->assembled) {
    ierr = MPISBSTRM_create_sbstrm(B);CHKERRQ(ierr);
  }

  ierr = PetscObjectChangeTypeName((PetscObject) B, MATMPISBSTRM);CHKERRQ(ierr);
  ierr = PetscObjectComposeFunction((PetscObject)B,"MatMPISBAIJSetPreallocation_C",MatMPISBAIJSetPreallocation_MPISBSTRM);CHKERRQ(ierr);
  *newmat = B;
  PetscFunctionReturn(0);
}

PETSC_EXTERN PetscErrorCode MatConvert_MPIAIJ_MPIAIJCRL(Mat A,MatType type,MatReuse reuse,Mat *newmat)
{
  PetscErrorCode ierr;
  Mat            B = *newmat;
  Mat_AIJCRL     *aijcrl;

  PetscFunctionBegin;
  if (reuse == MAT_INITIAL_MATRIX) {
    ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr);
  }

  ierr     = PetscNewLog(B,Mat_AIJCRL,&aijcrl);CHKERRQ(ierr);
  B->spptr = (void*) aijcrl;

  /* Set function pointers for methods that we inherit from AIJ but override. */
  B->ops->duplicate   = MatDuplicate_AIJCRL;
  B->ops->assemblyend = MatAssemblyEnd_MPIAIJCRL;
  B->ops->destroy     = MatDestroy_MPIAIJCRL;
  B->ops->mult        = MatMult_AIJCRL;

  /* If A has already been assembled, compute the permutation. */
  if (A->assembled) {
    ierr = MatMPIAIJCRL_create_aijcrl(B);CHKERRQ(ierr);
  }
  ierr    = PetscObjectChangeTypeName((PetscObject)B,MATMPIAIJCRL);CHKERRQ(ierr);
  *newmat = B;
  PetscFunctionReturn(0);
}

PetscErrorCode SNESMonitorJacUpdateSpectrum(SNES snes,PetscInt it,PetscReal fnorm,void *ctx)
{
#if defined(PETSC_MISSING_LAPACK_GEEV)
  SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_SUP,"GEEV - Lapack routine is unavailable\nNot able to provide eigen values.");
#elif defined(PETSC_HAVE_ESSL)
  SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_SUP,"GEEV - No support for ESSL Lapack Routines");
#else
  Vec            X;
  Mat            J,dJ,dJdense;
  PetscErrorCode ierr;
  PetscErrorCode (*func)(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
  PetscInt       n,i;
  PetscBLASInt   nb,lwork;
  PetscReal      *eigr,*eigi;
  MatStructure   flg = DIFFERENT_NONZERO_PATTERN;
  PetscScalar    *work;
  PetscScalar    *a;

  PetscFunctionBegin;
  if (it == 0) PetscFunctionReturn(0);
  /* create the difference between the current update and the current jacobian */
  ierr = SNESGetSolution(snes,&X);CHKERRQ(ierr);
  ierr = SNESGetJacobian(snes,&J,NULL,&func,NULL);CHKERRQ(ierr);
  ierr = MatDuplicate(J,MAT_COPY_VALUES,&dJ);CHKERRQ(ierr);
  ierr = SNESComputeJacobian(snes,X,&dJ,&dJ,&flg);CHKERRQ(ierr);
  ierr = MatAXPY(dJ,-1.0,J,SAME_NONZERO_PATTERN);CHKERRQ(ierr);

  /* compute the spectrum directly */
  ierr  = MatConvert(dJ,MATSEQDENSE,MAT_INITIAL_MATRIX,&dJdense);CHKERRQ(ierr);
  ierr  = MatGetSize(dJ,&n,NULL);CHKERRQ(ierr);
  ierr  = PetscBLASIntCast(n,&nb);CHKERRQ(ierr);
  lwork = 3*nb;
  ierr  = PetscMalloc(n*sizeof(PetscReal),&eigr);CHKERRQ(ierr);
  ierr  = PetscMalloc(n*sizeof(PetscReal),&eigi);CHKERRQ(ierr);
  ierr  = PetscMalloc(lwork*sizeof(PetscScalar),&work);CHKERRQ(ierr);
  ierr  = MatDenseGetArray(dJdense,&a);CHKERRQ(ierr);
#if !defined(PETSC_USE_COMPLEX)
  {
    PetscBLASInt lierr;
    ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
    PetscStackCall("LAPACKgeev",LAPACKgeev_("N","N",&nb,a,&nb,eigr,eigi,NULL,&nb,NULL,&nb,work,&lwork,&lierr));
    if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"geev() error %d",lierr);
    ierr = PetscFPTrapPop();CHKERRQ(ierr);
  }
#else
  SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not coded for complex");
#endif
  PetscPrintf(PetscObjectComm((PetscObject)snes),"Eigenvalues of J_%d - J_%d:\n",it,it-1);CHKERRQ(ierr);
  for (i=0;i<n;i++) {
    PetscPrintf(PetscObjectComm((PetscObject)snes),"%5d: %20.5g + %20.5gi\n",i,eigr[i],eigi[i]);CHKERRQ(ierr);
  }
  ierr = MatDenseRestoreArray(dJdense,&a);CHKERRQ(ierr);
  ierr = MatDestroy(&dJ);CHKERRQ(ierr);
  ierr = MatDestroy(&dJdense);CHKERRQ(ierr);
  ierr = PetscFree(eigr);CHKERRQ(ierr);
  ierr = PetscFree(eigi);CHKERRQ(ierr);
  ierr = PetscFree(work);CHKERRQ(ierr);
  PetscFunctionReturn(0);
#endif
}

PETSC_EXTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqSBSTRM(Mat A,MatType type,MatReuse reuse,Mat *newmat)
{
  PetscErrorCode ierr;

  Mat           B = *newmat;
  Mat_SeqSBSTRM *sbstrm;
  /* PetscInt       bs = A->rmap->bs; */

  PetscFunctionBegin;
  if (reuse == MAT_INITIAL_MATRIX) {
    ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr);
  }


  ierr     = PetscNewLog(B,&sbstrm);CHKERRQ(ierr);
  B->spptr = (void*) sbstrm;

  /* Set function pointers for methods that we inherit from BAIJ but override. */
  B->ops->duplicate   = MatDuplicate_SeqSBSTRM;
  B->ops->assemblyend = MatAssemblyEnd_SeqSBSTRM;
  B->ops->destroy     = MatDestroy_SeqSBSTRM;
  /*B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_sbstrm;
    B->ops->choleskyfactornumeric  = MatCholeskyFactorNumeric_sbstrm; */

  /* If A has already been assembled, compute the permutation. */
  if (A->assembled) {
    ierr = SeqSBSTRM_create_sbstrm(B);CHKERRQ(ierr);
  }
  ierr    = PetscObjectChangeTypeName((PetscObject)B,MATSEQSBSTRM);CHKERRQ(ierr);
  *newmat = B;
  PetscFunctionReturn(0);
}

/*@
  MatSchurComplementComputeExplicitOperator - Compute the Schur complement matrix explicitly

  Collective on Mat

  Input Parameter:
. M - the matrix obtained with MatCreateSchurComplement()

  Output Parameter:
. S - the Schur complement matrix

  Note: This can be expensive, so it is mainly for testing

  Level: advanced

.seealso: MatCreateSchurComplement(), MatSchurComplementUpdate()
@*/
PetscErrorCode MatSchurComplementComputeExplicitOperator(Mat M, Mat *S)
{
  Mat            B, C, D;
  KSP            ksp;
  PC             pc;
  PetscBool      isLU, isILU;
  PetscReal      fill = 2.0;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  ierr = MatSchurComplementGetSubMatrices(M, NULL, NULL, &B, &C, &D);CHKERRQ(ierr);
  ierr = MatSchurComplementGetKSP(M, &ksp);CHKERRQ(ierr);
  ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr);
  ierr = PetscObjectTypeCompare((PetscObject) pc, PCLU, &isLU);CHKERRQ(ierr);
  ierr = PetscObjectTypeCompare((PetscObject) pc, PCILU, &isILU);CHKERRQ(ierr);
  if (isLU || isILU) {
    Mat       fact, Bd, AinvB, AinvBd;
    PetscReal eps = 1.0e-10;

    /* This can be sped up for banded LU */
    ierr = KSPSetUp(ksp);CHKERRQ(ierr);
    ierr = PCFactorGetMatrix(pc, &fact);CHKERRQ(ierr);
    ierr = MatConvert(B, MATDENSE, MAT_INITIAL_MATRIX, &Bd);CHKERRQ(ierr);
    ierr = MatDuplicate(Bd, MAT_DO_NOT_COPY_VALUES, &AinvBd);CHKERRQ(ierr);
    ierr = MatMatSolve(fact, Bd, AinvBd);CHKERRQ(ierr);
    ierr = MatDestroy(&Bd);CHKERRQ(ierr);
    ierr = MatChop(AinvBd, eps);CHKERRQ(ierr);
    ierr = MatConvert(AinvBd, MATAIJ, MAT_INITIAL_MATRIX, &AinvB);CHKERRQ(ierr);
    ierr = MatDestroy(&AinvBd);CHKERRQ(ierr);
    ierr = MatMatMult(C, AinvB, MAT_INITIAL_MATRIX, fill, S);CHKERRQ(ierr);
    ierr = MatDestroy(&AinvB);CHKERRQ(ierr);
  } else {
    Mat Ainvd, Ainv;

    ierr = PCComputeExplicitOperator(pc, &Ainvd);CHKERRQ(ierr);
    ierr = MatConvert(Ainvd, MATAIJ, MAT_INITIAL_MATRIX, &Ainv);CHKERRQ(ierr);
    ierr = MatDestroy(&Ainvd);CHKERRQ(ierr);
#if 0
    /* Symmetric version */
    ierr = MatPtAP(Ainv, B, MAT_INITIAL_MATRIX, fill, S);CHKERRQ(ierr);
#else
    /* Nonsymmetric version */
    ierr = MatMatMatMult(C, Ainv, B, MAT_INITIAL_MATRIX, fill, S);CHKERRQ(ierr);
#endif
    ierr = MatDestroy(&Ainv);CHKERRQ(ierr);
  }

  ierr = PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr);
  ierr = MatView(*S, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  ierr = PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

  if (D) {
    MatInfo info;

    ierr = MatGetInfo(D, MAT_GLOBAL_SUM, &info);CHKERRQ(ierr);
    if (info.nz_used) SETERRQ(PetscObjectComm((PetscObject) M), PETSC_ERR_SUP, "Not yet implemented");
  }
  PetscFunctionReturn(0);
}

 PetscMatrix<Scalar>* PetscMatrix<Scalar>::duplicate() const
 {
   PetscMatrix<Scalar>*ptscmatrix = new PetscMatrix<Scalar>();
   MatDuplicate(matrix, MAT_COPY_VALUES, &(ptscmatrix->matrix));
   ptscmatrix->size = this->size;
   ptscmatrix->nnz = nnz;
   return ptscmatrix;
 };

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);
}

PetscErrorCode NEPSolve_Interpol(NEP nep)
{
  PetscErrorCode ierr;
  NEP_INTERPOL   *ctx = (NEP_INTERPOL*)nep->data;
  Mat            *A;   /*T=nep->function,Tp=nep->jacobian;*/
  PetscScalar    *x,*fx,t;
  PetscReal      *cs,a,b,s;
  PetscInt       i,j,k,deg=ctx->deg;

  PetscFunctionBegin;
  ierr = PetscMalloc4(deg+1,&A,(deg+1)*(deg+1),&cs,deg+1,&x,(deg+1)*nep->nt,&fx);CHKERRQ(ierr);
  ierr = RGIntervalGetEndpoints(nep->rg,&a,&b,NULL,NULL);CHKERRQ(ierr);
  ierr = ChebyshevNodes(deg,a,b,x,cs);CHKERRQ(ierr);
  for (j=0;j<nep->nt;j++) {
    for (i=0;i<=deg;i++) {
      ierr = FNEvaluateFunction(nep->f[j],x[i],&fx[i+j*(deg+1)]);CHKERRQ(ierr);
    }
  }

  /* Polynomial coefficients */
  for (k=0;k<=deg;k++) {
    ierr = MatDuplicate(nep->A[0],MAT_COPY_VALUES,&A[k]);CHKERRQ(ierr);
    t = 0.0;
    for (i=0;i<deg+1;i++) t += fx[i]*cs[i*(deg+1)+k];
    t *= 2.0/(deg+1); 
    if (k==0) t /= 2.0;
    ierr = MatScale(A[k],t);CHKERRQ(ierr);
    for (j=1;j<nep->nt;j++) {
      t = 0.0;
      for (i=0;i<deg+1;i++) t += fx[i+j*(deg+1)]*cs[i*(deg+1)+k];
      t *= 2.0/(deg+1); 
      if (k==0) t /= 2.0;
      ierr = MatAXPY(A[k],t,nep->A[j],SUBSET_NONZERO_PATTERN);CHKERRQ(ierr);
    }
  }

  ierr = PEPSetOperators(ctx->pep,deg+1,A);CHKERRQ(ierr);
  for (k=0;k<=deg;k++) {
    ierr = MatDestroy(&A[k]);CHKERRQ(ierr);
  }
  ierr = PetscFree4(A,cs,x,fx);CHKERRQ(ierr);

  /* Solve polynomial eigenproblem */
  ierr = PEPSolve(ctx->pep);CHKERRQ(ierr);
  ierr = PEPGetConverged(ctx->pep,&nep->nconv);CHKERRQ(ierr);
  ierr = PEPGetIterationNumber(ctx->pep,&nep->its);CHKERRQ(ierr);
  ierr = PEPGetConvergedReason(ctx->pep,(PEPConvergedReason*)&nep->reason);CHKERRQ(ierr);
  s = 2.0/(b-a);
  for (i=0;i<nep->nconv;i++) {
    ierr = PEPGetEigenpair(ctx->pep,i,&nep->eigr[i],&nep->eigi[i],NULL,NULL);CHKERRQ(ierr);
    nep->eigr[i] /= s;
    nep->eigr[i] += (a+b)/2.0;
    nep->eigi[i] /= s;
  }
  nep->state = NEP_STATE_EIGENVECTORS;
  PetscFunctionReturn(0);
}

PetscErrorCode MatDuplicate_Transpose(Mat N, MatDuplicateOption op, Mat* m)
{
  Mat_Transpose  *Na = (Mat_Transpose*)N->data;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  if (op == MAT_COPY_VALUES) {
    ierr = MatTranspose(Na->A,MAT_INITIAL_MATRIX,m);CHKERRQ(ierr);
  } else if (op == MAT_DO_NOT_COPY_VALUES) {
    ierr = MatDuplicate(Na->A,MAT_DO_NOT_COPY_VALUES,m);CHKERRQ(ierr);
    ierr = MatTranspose(*m,MAT_INPLACE_MATRIX,m);CHKERRQ(ierr);
  } else SETERRQ(PetscObjectComm((PetscObject)N),PETSC_ERR_SUP,"MAT_SHARE_NONZERO_PATTERN not supported for this matrix type");
  PetscFunctionReturn(0);
}

static PetscErrorCode MatCopy_User(Mat A,Mat B,MatStructure str)
{
  User           userA,userB;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  ierr = MatShellGetContext(A,&userA);CHKERRQ(ierr);
  if (userA) {
    ierr = PetscNew(&userB);CHKERRQ(ierr);
    ierr = MatDuplicate(userA->A,MAT_COPY_VALUES,&userB->A);CHKERRQ(ierr);
    ierr = MatShellSetContext(B, userB);CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

/*@C
  TaoMatGetSubMat - Gets a submatrix using the IS

  Input Parameters:
+ M - the full matrix (n x n)
. is - the index set for the submatrix (both row and column index sets need to be the same)
. v1 - work vector of dimension n, needed for TAO_SUBSET_MASK option
- subset_type - the method TAO is using for subsetting (TAO_SUBSET_SUBVEC, TAO_SUBSET_MASK,
  TAO_SUBSET_MATRIXFREE)

  Output Parameters:
. Msub - the submatrix
@*/
PetscErrorCode TaoMatGetSubMat(Mat M, IS is, Vec v1, TaoSubsetType subset_type, Mat *Msub)
{
  PetscErrorCode ierr;
  IS             iscomp;
  PetscBool      flg = PETSC_FALSE;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(M,MAT_CLASSID,1);
  PetscValidHeaderSpecific(is,IS_CLASSID,2);
  ierr = MatDestroy(Msub);CHKERRQ(ierr);
  switch (subset_type) {
  case TAO_SUBSET_SUBVEC:
    ierr = MatGetSubMatrix(M, is, is, MAT_INITIAL_MATRIX, Msub);CHKERRQ(ierr);
    break;

  case TAO_SUBSET_MASK:
    /* Get Reduced Hessian
     Msub[i,j] = M[i,j] if i,j in Free_Local or i==j
     Msub[i,j] = 0      if i!=j and i or j not in Free_Local
     */
    ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)M),NULL,NULL,NULL);CHKERRQ(ierr);
    ierr = PetscOptionsBool("-different_submatrix","use separate hessian matrix when computing submatrices","TaoSubsetType",flg,&flg,NULL);CHKERRQ(ierr);
    ierr = PetscOptionsEnd();CHKERRQ(ierr);
    if (flg) {
      ierr = MatDuplicate(M, MAT_COPY_VALUES, Msub);CHKERRQ(ierr);
    } else {
      /* Act on hessian directly (default) */
      ierr = PetscObjectReference((PetscObject)M);CHKERRQ(ierr);
      *Msub = M;
    }
    /* Save the diagonal to temporary vector */
    ierr = MatGetDiagonal(*Msub,v1);CHKERRQ(ierr);

    /* Zero out rows and columns */
    ierr = ISComplementVec(is,v1,&iscomp);CHKERRQ(ierr);

    /* Use v1 instead of 0 here because of PETSc bug */
    ierr = MatZeroRowsColumnsIS(*Msub,iscomp,1.0,v1,v1);CHKERRQ(ierr);

    ierr = ISDestroy(&iscomp);CHKERRQ(ierr);
    break;
  case TAO_SUBSET_MATRIXFREE:
    ierr = ISComplementVec(is,v1,&iscomp);CHKERRQ(ierr);
    ierr = MatCreateSubMatrixFree(M,iscomp,iscomp,Msub);CHKERRQ(ierr);
    ierr = ISDestroy(&iscomp);CHKERRQ(ierr);
    break;
  }
  PetscFunctionReturn(0);
}

int main(int argc,char **args)
{
  PetscErrorCode ierr;
  Mat            A,B,C;
  PetscBool      different=PETSC_FALSE,skip=PETSC_FALSE;
  PetscInt       m0,m1,n=128,i;

  PetscInitialize(&argc,&args,(char*)0,help);
  ierr = PetscOptionsGetBool(NULL,"-different",&different,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsGetBool(NULL,"-skip",&skip,NULL);CHKERRQ(ierr);
  /*
     Create matrices
     A = tridiag(1,-2,1) and B = diag(7);
  */
  ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
  ierr = MatCreate(PETSC_COMM_WORLD,&B);CHKERRQ(ierr);
  ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr);
  ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr);
  ierr = MatSetFromOptions(A);CHKERRQ(ierr);
  ierr = MatSetFromOptions(B);CHKERRQ(ierr);
  ierr = MatSetUp(A);CHKERRQ(ierr);
  ierr = MatSetUp(B);CHKERRQ(ierr);
  ierr = MatGetOwnershipRange(A,&m0,&m1);CHKERRQ(ierr);
  for (i=m0;i<m1;i++) {
    if (i>0) { ierr = MatSetValue(A,i,i-1,-1.0,INSERT_VALUES);CHKERRQ(ierr); }
    if (i<n-1) { ierr = MatSetValue(A,i,i+1,-1.0,INSERT_VALUES);CHKERRQ(ierr); }
    ierr = MatSetValue(A,i,i,2.0,INSERT_VALUES);CHKERRQ(ierr);
    ierr = MatSetValue(B,i,i,7.0,INSERT_VALUES);CHKERRQ(ierr);
  }
  ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

  ierr = MatDuplicate(A,MAT_COPY_VALUES,&C);CHKERRQ(ierr);
  /* Add B */
  ierr = MatAXPY(C,1.0,B,(different)?DIFFERENT_NONZERO_PATTERN:SUBSET_NONZERO_PATTERN);CHKERRQ(ierr);
  /* Add A */
  if (!skip) { ierr = MatAXPY(C,1.0,A,SUBSET_NONZERO_PATTERN);CHKERRQ(ierr); }

  /*
     Free memory
  */
  ierr = MatDestroy(&A);CHKERRQ(ierr);
  ierr = MatDestroy(&B);CHKERRQ(ierr);
  ierr = MatDestroy(&C);CHKERRQ(ierr);
  ierr = PetscFinalize();
  return 0;
}

PetscErrorCode TestMatZeroRows_with_no_allocation(Mat A,IS is,PetscScalar diag)
{
  Mat            B;
  PetscErrorCode ierr;

  /* Now copy A into B, and test it with MatZeroRows() */
  ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr);
  /* Set this flag after assembly. This way, it affects only MatZeroRows() */
  ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);

  ierr = MatZeroRowsIS(B,is,diag,0,0);CHKERRQ(ierr);
  ierr = MatView(B,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  ierr = MatDestroy(&B);CHKERRQ(ierr);
  return 0;
}

PETSC_INTERN PetscErrorCode MatConvert_MPIAIJ_MPIAIJMKL(Mat A,MatType type,MatReuse reuse,Mat *newmat)
{
  PetscErrorCode ierr;
  Mat            B = *newmat;

  PetscFunctionBegin;
  if (reuse == MAT_INITIAL_MATRIX) {
    ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr);
  }

  ierr = PetscObjectChangeTypeName((PetscObject) B, MATMPIAIJMKL);CHKERRQ(ierr);
  ierr = PetscObjectComposeFunction((PetscObject)B,"MatMPIAIJSetPreallocation_C",MatMPIAIJSetPreallocation_MPIAIJMKL);CHKERRQ(ierr);
  *newmat = B;
  PetscFunctionReturn(0);
}