C++ MatCreateVecs函数代码示例

PetscErrorCode SolveInit(FEMInf fem, int L, PetscScalar *e0, Vec *x) {

  PetscErrorCode ierr;
  Mat H, S;
  ierr = CalcMat(fem, L, &H, &S); CHKERRQ(ierr);

  EPS eps;
  ierr = PrintTimeStamp(fem->comm, "EPS", NULL); CHKERRQ(ierr);
  ierr = EPSCreate(fem->comm, &eps); CHKERRQ(ierr);
  ierr = EPSSetTarget(eps, -0.6); CHKERRQ(ierr);
  ierr = EPSSetWhichEigenpairs(eps, EPS_TARGET_MAGNITUDE);  CHKERRQ(ierr);
  ierr = EPSSetOperators(eps, H, S); CHKERRQ(ierr);
  if(S == NULL) {
    ierr = EPSSetProblemType(eps, EPS_NHEP);  CHKERRQ(ierr);
  } else {
    ierr = EPSSetProblemType(eps, EPS_GNHEP);  CHKERRQ(ierr);
  }
  Vec x0[1]; MatCreateVecs(H, &x0[0], NULL); 
  int num; FEMInfGetSize(fem, &num);
  for(int i = 0; i < num; i++) {
    VecSetValue(x0[0], i, 0.5, INSERT_VALUES);
  }
  VecAssemblyBegin(x0[0]); VecAssemblyEnd(x0[0]);
  EPSSetInitialSpace(eps, 1, x0);
  ierr = EPSSetFromOptions(eps); CHKERRQ(ierr);

  ierr = EPSSolve(eps); CHKERRQ(ierr);
  PetscInt nconv;
  EPSGetConverged(eps, &nconv);

  if(nconv == 0) 
    SETERRQ(fem->comm, 1, "Failed to digonalize in init state\n");

  Vec x_ans;
  MatCreateVecs(H, &x_ans, NULL);
  EPSGetEigenpair(eps, 0, e0, NULL, x_ans, NULL);

  EPSDestroy(&eps);

  PetscScalar v[1]; PetscInt idx[1] = {1};
  VecGetValues(x_ans, 1, idx, v);
  PetscScalar scale_factor = v[0] / cabs(v[0]);
  VecScale( x_ans, 1.0/scale_factor);

  PetscScalar norm0;
  Vec Sx;  MatCreateVecs(S, &Sx, NULL); 
  MatMult(S, x_ans, Sx); VecDot(x_ans, Sx, &norm0);

  VecScale(x_ans, 1.0/sqrt(norm0));

  *x = x_ans;
  return 0;
}

static void ApplyPCPrecPETSCGen(void *x, PRIMME_INT *ldx, void *y,
      PRIMME_INT *ldy, int *blockSize, int trans, PC *pc, MPI_Comm comm) {
   int i;
   Vec xvec, yvec;
   Mat matrix;
   PetscErrorCode ierr;
   PetscInt mLocal, nLocal;
   
   ierr = PCGetOperators(pc[0],&matrix,NULL); CHKERRABORT(comm, ierr);

   assert(sizeof(PetscScalar) == sizeof(SCALAR));   
   ierr = MatGetLocalSize(matrix, &mLocal, &nLocal); CHKERRABORT(comm, ierr);
   assert(mLocal == nLocal && nLocal <= *ldx && mLocal <= *ldy);

   #if PETSC_VERSION_LT(3,6,0)
      ierr = MatGetVecs(matrix, &xvec, &yvec); CHKERRABORT(comm, ierr);
   #else
      ierr = MatCreateVecs(matrix, &xvec, &yvec); CHKERRABORT(comm, ierr);
   #endif
   for (i=0; i<*blockSize; i++) {
      ierr = VecPlaceArray(xvec, ((SCALAR*)x) + (*ldx)*i); CHKERRABORT(comm, ierr);
      ierr = VecPlaceArray(yvec, ((SCALAR*)y) + (*ldy)*i); CHKERRABORT(comm, ierr);
      if (trans == 0) {
         ierr = PCApply(*pc, xvec, yvec); CHKERRABORT(comm, ierr);
      } else if (pc[1]) {
         ierr = PCApply(pc[1], xvec, yvec); CHKERRABORT(comm, ierr);
      } else {
         ierr = PCApplyTranspose(pc[0], xvec, yvec);
      }
      ierr = VecResetArray(xvec); CHKERRABORT(comm, ierr);
      ierr = VecResetArray(yvec); CHKERRABORT(comm, ierr);
   }
   ierr = VecDestroy(&xvec); CHKERRABORT(comm, ierr);
   ierr = VecDestroy(&yvec); CHKERRABORT(comm, ierr);
}

PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping rmapping,ISLocalToGlobalMapping cmapping)
{
  PetscErrorCode ierr;
  PetscInt       n,bs;
  Mat_IS         *is = (Mat_IS*)A->data;
  IS             from,to;
  Vec            global;

  PetscFunctionBegin;
  PetscCheckSameComm(A,1,rmapping,2);
  if (rmapping != cmapping) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"MATIS requires the row and column mappings to be identical");
  if (is->mapping) { /* Currenly destroys the objects that will be created by this routine. Is there anything else that should be checked? */
    ierr = ISLocalToGlobalMappingDestroy(&is->mapping);CHKERRQ(ierr);
    ierr = VecDestroy(&is->x);CHKERRQ(ierr);
    ierr = VecDestroy(&is->y);CHKERRQ(ierr);
    ierr = VecScatterDestroy(&is->ctx);CHKERRQ(ierr);
    ierr = MatDestroy(&is->A);CHKERRQ(ierr);
  }
  ierr = PetscObjectReference((PetscObject)rmapping);CHKERRQ(ierr);
  ierr = ISLocalToGlobalMappingDestroy(&is->mapping);CHKERRQ(ierr);
  is->mapping = rmapping;
/*
  ierr = PetscLayoutSetISLocalToGlobalMapping(A->rmap,rmapping);CHKERRQ(ierr);
  ierr = PetscLayoutSetISLocalToGlobalMapping(A->cmap,cmapping);CHKERRQ(ierr);
*/

  /* Create the local matrix A */
  ierr = ISLocalToGlobalMappingGetSize(rmapping,&n);CHKERRQ(ierr);
  ierr = ISLocalToGlobalMappingGetBlockSize(rmapping,&bs);CHKERRQ(ierr);
  ierr = MatCreate(PETSC_COMM_SELF,&is->A);CHKERRQ(ierr);
  if (bs > 1) {
    ierr = MatSetType(is->A,MATSEQBAIJ);CHKERRQ(ierr);
  } else {
    ierr = MatSetType(is->A,MATSEQAIJ);CHKERRQ(ierr);
  }
  ierr = MatSetSizes(is->A,n,n,n,n);CHKERRQ(ierr);
  ierr = MatSetBlockSize(is->A,bs);CHKERRQ(ierr);
  ierr = MatSetOptionsPrefix(is->A,((PetscObject)A)->prefix);CHKERRQ(ierr);
  ierr = MatAppendOptionsPrefix(is->A,"is_");CHKERRQ(ierr);
  ierr = MatSetFromOptions(is->A);CHKERRQ(ierr);

  /* Create the local work vectors */
  ierr = VecCreate(PETSC_COMM_SELF,&is->x);CHKERRQ(ierr);
  ierr = VecSetBlockSize(is->x,bs);CHKERRQ(ierr);
  ierr = VecSetSizes(is->x,n,n);CHKERRQ(ierr);
  ierr = VecSetOptionsPrefix(is->x,((PetscObject)A)->prefix);CHKERRQ(ierr);
  ierr = VecAppendOptionsPrefix(is->x,"is_");CHKERRQ(ierr);
  ierr = VecSetFromOptions(is->x);CHKERRQ(ierr);
  ierr = VecDuplicate(is->x,&is->y);CHKERRQ(ierr);

  /* setup the global to local scatter */
  ierr = ISCreateStride(PETSC_COMM_SELF,n,0,1,&to);CHKERRQ(ierr);
  ierr = ISLocalToGlobalMappingApplyIS(rmapping,to,&from);CHKERRQ(ierr);
  ierr = MatCreateVecs(A,&global,NULL);CHKERRQ(ierr);
  ierr = VecScatterCreate(global,from,is->x,to,&is->ctx);CHKERRQ(ierr);
  ierr = VecDestroy(&global);CHKERRQ(ierr);
  ierr = ISDestroy(&to);CHKERRQ(ierr);
  ierr = ISDestroy(&from);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

static PetscErrorCode KSPSetUp_TSIRM(KSP ksp)
{
    PetscErrorCode ierr;
    KSP_TSIRM      *tsirm = (KSP_TSIRM*)ksp->data;

    PetscFunctionBegin;
    /* Initialization */
    tsirm->tol_ls     = 1e-40;
    tsirm->size_ls    = 12;
    tsirm->maxiter_ls = 15;
    tsirm->cgls       = 0;

    /* Matrix of the system */
    ierr = KSPGetOperators(ksp,&tsirm->A,NULL);
    CHKERRQ(ierr);    /* Matrix of the system   */
    ierr = MatGetSize(tsirm->A,&tsirm->size,NULL);
    CHKERRQ(ierr); /* Size of the system     */
    ierr = MatGetOwnershipRange(tsirm->A,&tsirm->Istart,&tsirm->Iend);
    CHKERRQ(ierr);

    /* Matrix S of residuals */
    ierr = MatCreate(PETSC_COMM_WORLD,&tsirm->S);
    CHKERRQ(ierr);
    ierr = MatSetSizes(tsirm->S,tsirm->Iend-tsirm->Istart,PETSC_DECIDE,tsirm->size,tsirm->size_ls);
    CHKERRQ(ierr);
    ierr = MatSetType(tsirm->S,MATDENSE);
    CHKERRQ(ierr);
    ierr = MatSetUp(tsirm->S);
    CHKERRQ(ierr);

    /* Residual and vector Alpha computed in the minimization step */
    ierr = MatCreateVecs(tsirm->S,&tsirm->Alpha,&tsirm->r);
    CHKERRQ(ierr);
    PetscFunctionReturn(0);
}

int main(int argc, char **argv)
{
  MPI_Comm       comm;
  Mat            A;    /* A graph */
  Vec            c;    /* A vector giving the component of each vertex */
  Vec            cold; /* The vector c from the last iteration */
  PetscScalar   *carray;
  PetscInt       testnum = 0;
  PetscInt       V, vStart, vEnd, v, n;
  PetscMPIInt    size;
  PetscErrorCode ierr;

  ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
  comm = PETSC_COMM_WORLD;
  ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr);
  /* Use matrix to encode a graph */
  ierr = PetscOptionsGetInt(NULL,NULL, "-testnum", &testnum, NULL);CHKERRQ(ierr);
  ierr = CreateGraph(comm, testnum, &A);CHKERRQ(ierr);
  ierr = MatGetSize(A, &V, NULL);CHKERRQ(ierr);
  /* Replace matrix-vector multiplication with one that calculates the minimum rather than the sum */
  if (size == 1) {
    ierr = MatShellSetOperation(A, MATOP_MULT, (void (*)) MatMultMax_SeqAIJ);CHKERRQ(ierr);
  } else {
    Mat_MPIAIJ *a = (Mat_MPIAIJ *) A->data;

    ierr = MatShellSetOperation(a->A, MATOP_MULT, (void (*)) MatMultMax_SeqAIJ);CHKERRQ(ierr);
    ierr = MatShellSetOperation(a->B, MATOP_MULT, (void (*)) MatMultMax_SeqAIJ);CHKERRQ(ierr);
    ierr = MatShellSetOperation(a->B, MATOP_MULT_ADD, (void (*)) MatMultAddMax_SeqAIJ);CHKERRQ(ierr);
  }
  /* Initialize each vertex as a separate component */
  ierr = MatCreateVecs(A, &c, NULL);CHKERRQ(ierr);
  ierr = MatGetOwnershipRange(A, &vStart, &vEnd);CHKERRQ(ierr);
  ierr = VecGetArray(c, &carray);CHKERRQ(ierr);
  for (v = vStart; v < vEnd; ++v) {
    carray[v-vStart] = v;
  }
  ierr = VecRestoreArray(c, &carray);CHKERRQ(ierr);
  /* Preprocess in parallel to find local components */
  /* Multiply until c does not change */
  ierr = VecDuplicate(c, &cold);CHKERRQ(ierr);
  for (v = 0; v < V; ++v) {
    Vec       cnew = cold;
    PetscBool stop;

    ierr = MatMult(A, c, cnew);CHKERRQ(ierr);
    ierr = VecEqual(c, cnew, &stop);CHKERRQ(ierr);
    if (stop) break;
    cold = c;
    c    = cnew;
  }
  /* Report */
  ierr = VecUniqueEntries(c, &n, NULL);CHKERRQ(ierr);
  ierr = PetscPrintf(comm, "Components: %d Iterations: %d\n", n, v);CHKERRQ(ierr);
  ierr = VecView(c, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  /* Cleanup */
  ierr = VecDestroy(&c);CHKERRQ(ierr);
  ierr = VecDestroy(&cold);CHKERRQ(ierr);
  ierr = PetscFinalize();
  return ierr;
}

int main(int argc,char **args)
{
  Mat            A;
  Vec            x;
  PetscErrorCode ierr;
  PetscViewer    fd;              /* viewer */
  char           file[PETSC_MAX_PATH_LEN]; /* input file name */
  PetscReal      norm;
  PetscBool      flg;

  PetscInitialize(&argc,&args,(char*)0,help);

  /* Determine file from which we read the matrix A */
  ierr = PetscOptionsGetString(NULL,"-f",file,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr);
  if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate binary file with the -f option");

  /* Load matrix A */
  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);CHKERRQ(ierr);
  ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
  ierr = MatLoad(A,fd);CHKERRQ(ierr);
  ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr);
  ierr = MatCreateVecs(A,&x,NULL);CHKERRQ(ierr);
  ierr = MatGetDiagonal(A,x);CHKERRQ(ierr);
  ierr = VecScale(x,-1.0);CHKERRQ(ierr);
  ierr = MatDiagonalSet(A,x,ADD_VALUES);CHKERRQ(ierr);
  ierr = MatGetDiagonal(A,x);CHKERRQ(ierr);
  ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm %g\n",(double)norm);CHKERRQ(ierr);

  /* Free data structures */
  ierr = MatDestroy(&A);CHKERRQ(ierr);
  ierr = VecDestroy(&x);CHKERRQ(ierr);
  ierr = PetscFinalize();
  return 0;
}

/*@
   PCMGGetRScale - Gets the pointwise scaling for the restriction operator from level l to l-1.

   Collective on PC

   Input Parameters:
+  pc - the multigrid context
.  rscale - the scaling
-  l - the level (0 is coarsest) to supply [Do not supply 0]

   Level: advanced

   Notes:
       When evaluating a function on a coarse level one does not want to do F(R * x) one does F(rscale * R * x) where rscale is 1 over the row sums of R.

.keywords: MG, set, multigrid, restriction, level

.seealso: PCMGSetInterpolation(), PCMGGetRestriction()
@*/
PetscErrorCode PCMGGetRScale(PC pc,PetscInt l,Vec *rscale)
{
  PetscErrorCode ierr;
  PC_MG          *mg        = (PC_MG*)pc->data;
  PC_MG_Levels   **mglevels = mg->levels;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc,PC_CLASSID,1);
  if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
  if (l <= 0 || mg->nlevels <= l) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Level %D must be in range {1,...,%D}",l,mg->nlevels-1);
  if (!mglevels[l]->rscale) {
    Mat      R;
    Vec      X,Y,coarse,fine;
    PetscInt M,N;
    ierr = PCMGGetRestriction(pc,l,&R);CHKERRQ(ierr);
    ierr = MatCreateVecs(R,&X,&Y);CHKERRQ(ierr);
    ierr = MatGetSize(R,&M,&N);CHKERRQ(ierr);
    if (M < N) {
      fine = X;
      coarse = Y;
    } else if (N < M) {
      fine = Y; coarse = X;
    } else SETERRQ(PetscObjectComm((PetscObject)R),PETSC_ERR_SUP,"Restriction matrix is square, cannot determine which Vec is coarser");
    ierr = VecSet(fine,1.);CHKERRQ(ierr);
    ierr = MatRestrict(R,fine,coarse);CHKERRQ(ierr);
    ierr = VecDestroy(&fine);CHKERRQ(ierr);
    ierr = VecReciprocal(coarse);CHKERRQ(ierr);
    mglevels[l]->rscale = coarse;
  }
  *rscale = mglevels[l]->rscale;
  PetscFunctionReturn(0);
}

PETSC_EXTERN PetscErrorCode MatMFFDSetBase_MFFD(Mat J,Vec U,Vec F)
{
  PetscErrorCode ierr;
  MatMFFD        ctx = (MatMFFD)J->data;

  PetscFunctionBegin;
  ierr = MatMFFDResetHHistory(J);CHKERRQ(ierr);
  if (!ctx->current_u) {
    ierr = VecDuplicate(U,&ctx->current_u);CHKERRQ(ierr);
    ierr = VecLockReadPush(ctx->current_u);CHKERRQ(ierr);
  }
  ierr = VecLockReadPop(ctx->current_u);CHKERRQ(ierr);
  ierr = VecCopy(U,ctx->current_u);CHKERRQ(ierr);
  ierr = VecLockReadPush(ctx->current_u);CHKERRQ(ierr);
  if (F) {
    if (ctx->current_f_allocated) {ierr = VecDestroy(&ctx->current_f);CHKERRQ(ierr);}
    ctx->current_f           = F;
    ctx->current_f_allocated = PETSC_FALSE;
  } else if (!ctx->current_f_allocated) {
    ierr = MatCreateVecs(J,NULL,&ctx->current_f);CHKERRQ(ierr);

    ctx->current_f_allocated = PETSC_TRUE;
  }
  if (!ctx->w) {
    ierr = VecDuplicate(ctx->current_u,&ctx->w);CHKERRQ(ierr);
  }
  J->assembled = PETSC_TRUE;
  PetscFunctionReturn(0);
}

int main(int argc,char **args)
{
  const PetscInt    inds[]  = {0,1};
  PetscScalar       avals[] = {2,3,5,7};
  Mat               S;
  User              user;
  PetscErrorCode    ierr;
  Vec               base;

  ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
  ierr = PetscNew(&user);CHKERRQ(ierr);
  ierr = MatCreateSeqAIJ(PETSC_COMM_WORLD,2,2,2,NULL,&user->B);CHKERRQ(ierr);
  ierr = MatSetUp(user->B);CHKERRQ(ierr);
  ierr = MatSetValues(user->B,2,inds,2,inds,avals,INSERT_VALUES);CHKERRQ(ierr);
  ierr = MatAssemblyBegin(user->B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(user->B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatCreateVecs(user->B,&base,NULL);CHKERRQ(ierr);
  ierr = MatCreateShell(PETSC_COMM_WORLD,2,2,2,2,user,&S);CHKERRQ(ierr);
  ierr = MatSetUp(S);CHKERRQ(ierr);
  ierr = MatShellSetOperation(S,MATOP_MULT,(void (*)(void))MatMult_User);CHKERRQ(ierr);
  ierr = MatShellSetOperation(S,MATOP_MULT_TRANSPOSE,(void (*)(void))MatMultTranspose_User);CHKERRQ(ierr);

  ierr = MatShellTestMult(S,MyFunction,base,user,NULL);CHKERRQ(ierr);
  ierr = MatShellTestMultTranspose(S,MyFunction,base,user,NULL);CHKERRQ(ierr);

  ierr = VecDestroy(&base);CHKERRQ(ierr);
  ierr = MatDestroy(&user->B);CHKERRQ(ierr);
  ierr = MatDestroy(&S);CHKERRQ(ierr);
  ierr = PetscFree(user);CHKERRQ(ierr);
  ierr = PetscFinalize();
  return ierr;
}

static void PETScMatvecGenNoBlock(void *x, PRIMME_INT ldx, void *y, PRIMME_INT ldy,
      int blockSize, int trans, Mat matrix, MPI_Comm comm) {
   int i;
   Vec xvec, yvec;
   PetscInt m, n, mLocal, nLocal;
   PetscErrorCode ierr;
   
   assert(sizeof(PetscScalar) == sizeof(SCALAR));   
   ierr = MatGetSize(matrix, &m, &n); CHKERRABORT(comm, ierr);
   ierr = MatGetLocalSize(matrix, &mLocal, &nLocal); CHKERRABORT(comm, ierr);

   #if PETSC_VERSION_LT(3,6,0)
      ierr = MatGetVecs(matrix, &xvec, &yvec); CHKERRABORT(comm, ierr);
   #else
      ierr = MatCreateVecs(matrix, &xvec, &yvec); CHKERRABORT(comm, ierr);
   #endif
   if (trans == 1) {
      Vec aux = xvec; xvec = yvec; yvec = aux;
   }
   for (i=0; i<blockSize; i++) {
      ierr = VecPlaceArray(xvec, ((SCALAR*)x) + ldx*i); CHKERRABORT(comm, ierr);
      ierr = VecPlaceArray(yvec, ((SCALAR*)y) + ldy*i); CHKERRABORT(comm, ierr);
      if (trans == 0) {
         ierr = MatMult(matrix, xvec, yvec); CHKERRABORT(comm, ierr);
      } else {
         ierr = MatMultHermitianTranspose(matrix, xvec, yvec); CHKERRABORT(comm, ierr);
      }
      ierr = VecResetArray(xvec); CHKERRABORT(comm, ierr);
      ierr = VecResetArray(yvec); CHKERRABORT(comm, ierr);
   }
   ierr = VecDestroy(&xvec); CHKERRABORT(comm, ierr);
   ierr = VecDestroy(&yvec); CHKERRABORT(comm, ierr);
}

static PetscErrorCode PCBDDCScalingSetUp_Deluxe_Private(PC pc)
{
  PC_BDDC                *pcbddc=(PC_BDDC*)pc->data;
  PCBDDCDeluxeScaling    deluxe_ctx=pcbddc->deluxe_ctx;
  PCBDDCSubSchurs        sub_schurs = pcbddc->sub_schurs;
  PetscErrorCode         ierr;

  PetscFunctionBegin;
  if (!sub_schurs->n_subs) {
    PetscFunctionReturn(0);
  }

  /* Create work vectors for sequential part of deluxe */
  ierr = MatCreateVecs(sub_schurs->S_Ej_all,&deluxe_ctx->seq_work1,&deluxe_ctx->seq_work2);CHKERRQ(ierr);

  /* Compute deluxe sequential scatter */
  if (sub_schurs->reuse_mumps && !sub_schurs->is_dir) {
    PCBDDCReuseMumps reuse_mumps = sub_schurs->reuse_mumps;
    ierr = PetscObjectReference((PetscObject)reuse_mumps->correction_scatter_B);CHKERRQ(ierr);
    deluxe_ctx->seq_scctx = reuse_mumps->correction_scatter_B;
  } else {
    ierr = VecScatterCreate(pcbddc->work_scaling,sub_schurs->is_Ej_all,deluxe_ctx->seq_work1,NULL,&deluxe_ctx->seq_scctx);CHKERRQ(ierr);
  }

  /* Create Mat object for deluxe scaling */
  ierr = PetscObjectReference((PetscObject)sub_schurs->S_Ej_all);CHKERRQ(ierr);
  deluxe_ctx->seq_mat = sub_schurs->S_Ej_all;
  if (sub_schurs->sum_S_Ej_all) { /* if this matrix is present, then we need to create the KSP object to invert it */
    PC               pc_temp;
    MatSolverPackage solver=NULL;
    char             ksp_prefix[256];
    size_t           len;

    ierr = KSPCreate(PETSC_COMM_SELF,&deluxe_ctx->seq_ksp);CHKERRQ(ierr);
    ierr = KSPSetOperators(deluxe_ctx->seq_ksp,sub_schurs->sum_S_Ej_all,sub_schurs->sum_S_Ej_all);CHKERRQ(ierr);
    ierr = KSPSetType(deluxe_ctx->seq_ksp,KSPPREONLY);CHKERRQ(ierr);
    ierr = KSPGetPC(deluxe_ctx->seq_ksp,&pc_temp);CHKERRQ(ierr);
    ierr = PCSetType(pc_temp,PCLU);CHKERRQ(ierr);
    ierr = KSPGetPC(pcbddc->ksp_D,&pc_temp);CHKERRQ(ierr);
    ierr = PCFactorGetMatSolverPackage(pc_temp,(const MatSolverPackage*)&solver);CHKERRQ(ierr);
    if (solver) {
      PC     new_pc;
      PCType type;

      ierr = PCGetType(pc_temp,&type);CHKERRQ(ierr);
      ierr = KSPGetPC(deluxe_ctx->seq_ksp,&new_pc);CHKERRQ(ierr);
      ierr = PCSetType(new_pc,type);CHKERRQ(ierr);
      ierr = PCFactorSetMatSolverPackage(new_pc,solver);CHKERRQ(ierr);
    }
    ierr = PetscStrlen(((PetscObject)(pcbddc->ksp_D))->prefix,&len);CHKERRQ(ierr);
    len -= 10; /* remove "dirichlet_" */
    ierr = PetscStrncpy(ksp_prefix,((PetscObject)(pcbddc->ksp_D))->prefix,len+1);CHKERRQ(ierr);
    ierr = PetscStrcat(ksp_prefix,"deluxe_");CHKERRQ(ierr);
    ierr = KSPSetOptionsPrefix(deluxe_ctx->seq_ksp,ksp_prefix);CHKERRQ(ierr);
    ierr = KSPSetFromOptions(deluxe_ctx->seq_ksp);CHKERRQ(ierr);
    ierr = KSPSetUp(deluxe_ctx->seq_ksp);CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

PetscErrorCode MatMult_SchurComplement(Mat N,Vec x,Vec y)
{
  Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
  PetscErrorCode      ierr;

  PetscFunctionBegin;
  if (!Na->work1) {ierr = MatCreateVecs(Na->A,&Na->work1,NULL);CHKERRQ(ierr);}
  if (!Na->work2) {ierr = MatCreateVecs(Na->A,&Na->work2,NULL);CHKERRQ(ierr);}
  ierr = MatMult(Na->B,x,Na->work1);CHKERRQ(ierr);
  ierr = KSPSolve(Na->ksp,Na->work1,Na->work2);CHKERRQ(ierr);
  ierr = MatMult(Na->C,Na->work2,y);CHKERRQ(ierr);
  ierr = VecScale(y,-1.0);CHKERRQ(ierr);
  if (Na->D) {
    ierr = MatMultAdd(Na->D,x,y,y);CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

PetscErrorCode MatCreateVecs_SchurComplement(Mat N,Vec *right,Vec *left)
{
  Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
  PetscErrorCode      ierr;

  PetscFunctionBegin;
  if (Na->D) {
    ierr = MatCreateVecs(Na->D,right,left);CHKERRQ(ierr);
    PetscFunctionReturn(0);
  }
  if (right) {
    ierr = MatCreateVecs(Na->B,right,NULL);CHKERRQ(ierr);
  }
  if (left) {
    ierr = MatCreateVecs(Na->C,NULL,left);CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

PetscErrorCode MatCreateVecs_Transpose(Mat A,Vec *r, Vec *l)
{
  Mat_Transpose  *Aa = (Mat_Transpose*)A->data;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  ierr = MatCreateVecs(Aa->A,l,r);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PetscErrorCode MatCreateVecs_HT(Mat N,Vec *r, Vec *l)
{
  Mat_HT         *Na = (Mat_HT*)N->data;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  ierr = MatCreateVecs(Na->A,l,r);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}