C++ MatGetVecs函数代码示例

PetscErrorCode bsscr_DGMiGtD( Mat *_K2, Mat K, Mat G, Mat M){
    Mat K2;
    Vec diag;
    Vec Mdiag;
    Mat MinvGt;
    Mat Gtrans;
    PetscErrorCode ierr;

    PetscFunctionBegin;
    MatGetVecs( K, &diag, PETSC_NULL );
    MatGetDiagonal( K, diag );
    VecSqrt(diag);
    //VecReciprocal(diag);/* trying something different here */
    MatGetVecs( M, &Mdiag, PETSC_NULL );
    MatGetDiagonal( M, Mdiag );
    VecReciprocal(Mdiag);
    #if( PETSC_VERSION_MAJOR <= 2 )
    ierr=MatTranspose(G, &Gtrans);CHKERRQ(ierr);
    #else
    ierr=MatTranspose(G, MAT_INITIAL_MATRIX,&Gtrans);CHKERRQ(ierr);
    #endif
    ierr=MatConvert(Gtrans, MATSAME, MAT_INITIAL_MATRIX, &MinvGt);CHKERRQ(ierr);/* copy Gtrans -> MinvGt */
    MatDiagonalScale(MinvGt, Mdiag, PETSC_NULL);/* Minv*Gtrans */
    /* MAT_INITIAL_MATRIX -> creates K2 matrix : PETSC_DEFAULT for fill ratio: run with -info to find what it should be*/
    ierr=MatMatMult( G, MinvGt, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &K2);CHKERRQ(ierr);/* K2 = G*Minv*Gtrans */
    MatDiagonalScale(K2, diag, diag );/* K2 = D*K2*D  = D*G*Minv*Gtrans*D */
    Stg_MatDestroy(&Gtrans);
    Stg_MatDestroy(&MinvGt);
    Stg_VecDestroy(&diag);
    Stg_VecDestroy(&Mdiag);
    *_K2=K2;
    PetscFunctionReturn(0);
}

PetscErrorCode KSPBuildPressure_CB_Nullspace_BSSCR(KSP ksp)
{
    KSP_BSSCR        *bsscr = (KSP_BSSCR *)ksp->data;
    FeEquationNumber *eq_num = bsscr->solver->st_sle->pSolnVec->feVariable->eqNum;
    FeMesh           *feMesh = bsscr->solver->st_sle->pSolnVec->feVariable->feMesh; /* is the pressure mesh */
    unsigned          ijk[3];
    Vec               t, v;
    int numLocalNodes, globalNodeNumber, i, j, eq;
    MatStokesBlockScaling BA = bsscr->BA;
    PetscErrorCode ierr;
    Mat            Amat,Pmat, G;
    MatStructure   pflag;

    PetscFunctionBegin;
    /* get G matrix from Amat matrix operator on ksp */
    ierr=Stg_PCGetOperators(ksp->pc,&Amat,&Pmat,&pflag);CHKERRQ(ierr);
    MatNestGetSubMat( Amat, 0,1, &G );/* G should always exist */
    /* now create Vecs t and v to match size of G: i.e. pressure */ /* NOTE: not using "h" vector from ksp->vec_rhs because this part of the block vector doesn't always exist */
    MatGetVecs( G, &t, PETSC_NULL );/* t and v are destroyed in KSPDestroy_BSSCR */
    MatGetVecs( G, &v, PETSC_NULL );/* t and v such that can do G*t */

    numLocalNodes = Mesh_GetLocalSize( feMesh,  MT_VERTEX); /* number of nodes on current proc not counting any shadow nodes */
    for(j=0;j<numLocalNodes;j++){
	i = globalNodeNumber = Mesh_DomainToGlobal( feMesh, MT_VERTEX, j);
	RegularMeshUtils_Element_1DTo3D(feMesh, i, ijk);
	eq = eq_num->destinationArray[j][0];/* get global equation number -- 2nd arg is always 0 because pressure has only one dof */
	if(eq != -1){
	    if( (ijk[0]+ijk[1]+ijk[2])%2 ==0 ){
		VecSetValue(t,eq,1.0,INSERT_VALUES);
	    }
	    else{
		VecSetValue(v,eq,1.0,INSERT_VALUES);	    }}}

    VecAssemblyBegin( t );
    VecAssemblyEnd( t );
    VecAssemblyBegin( v );
    VecAssemblyEnd( v );

    /* Scaling the null vectors here because it easier at the moment *//* maybe should do this in the original scaling function */
    if( BA->scaling_exists == PETSC_TRUE ){
	Vec R2;
	/* Get the scalings out the block mat data */
	VecNestGetSubVec( BA->Rz, 1, &R2 );
	VecPointwiseDivide( t, t, R2); /* x <- x * 1/R2 */
	VecPointwiseDivide( v, v, R2);
    }

    bsscr_writeVec( t, "t", "Writing t vector");
    bsscr_writeVec( v, "v", "Writing v vector");

    bsscr->t=t;
    bsscr->v=v;

    PetscFunctionReturn(0);
}

/*@C
  KSPGetVecs - Gets a number of work vectors.

  Input Parameters:
+ ksp  - iterative context
. rightn  - number of right work vectors
- leftn   - number of left work vectors to allocate

  Output Parameter:
+  right - the array of vectors created
-  left - the array of left vectors

   Note: The right vector has as many elements as the matrix has columns. The left
     vector has as many elements as the matrix has rows.

   Level: advanced

.seealso:   MatGetVecs()

@*/
PetscErrorCode KSPGetVecs(KSP ksp,PetscInt rightn, Vec **right,PetscInt leftn,Vec **left)
{
  PetscErrorCode ierr;
  Vec            vecr,vecl;

  PetscFunctionBegin;
  if (rightn) {
    if (!right) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"You asked for right vectors but did not pass a pointer to hold them");
    if (ksp->vec_sol) vecr = ksp->vec_sol;
    else {
      if (ksp->dm) {
        ierr = DMGetGlobalVector(ksp->dm,&vecr);CHKERRQ(ierr);
      } else {
        Mat mat;
        if (!ksp->pc) {ierr = KSPGetPC(ksp,&ksp->pc);CHKERRQ(ierr);}
        ierr = PCGetOperators(ksp->pc,&mat,NULL,NULL);CHKERRQ(ierr);
        ierr = MatGetVecs(mat,&vecr,NULL);CHKERRQ(ierr);
      }
    }
    ierr = VecDuplicateVecs(vecr,rightn,right);CHKERRQ(ierr);
    if (!ksp->vec_sol) {
      if (ksp->dm) {
        ierr = DMRestoreGlobalVector(ksp->dm,&vecr);CHKERRQ(ierr);
      } else {
        ierr = VecDestroy(&vecr);CHKERRQ(ierr);
      }
    }
  }
  if (leftn) {
    if (!left) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"You asked for left vectors but did not pass a pointer to hold them");
    if (ksp->vec_rhs) vecl = ksp->vec_rhs;
    else {
      if (ksp->dm) {
        ierr = DMGetGlobalVector(ksp->dm,&vecl);CHKERRQ(ierr);
      } else {
        Mat mat;
        if (!ksp->pc) {ierr = KSPGetPC(ksp,&ksp->pc);CHKERRQ(ierr);}
        ierr = PCGetOperators(ksp->pc,&mat,NULL,NULL);CHKERRQ(ierr);
        ierr = MatGetVecs(mat,NULL,&vecl);CHKERRQ(ierr);
      }
    }
    ierr = VecDuplicateVecs(vecl,leftn,left);CHKERRQ(ierr);
    if (!ksp->vec_rhs) {
      if (ksp->dm) {
        ierr = DMRestoreGlobalVector(ksp->dm,&vecl);CHKERRQ(ierr);
      } else {
        ierr = VecDestroy(&vecl);CHKERRQ(ierr);
      }
    }
  }
  PetscFunctionReturn(0);
}

static PetscErrorCode PCSetUp_Eisenstat(PC pc)
{
  PetscErrorCode ierr;
  PetscInt       M,N,m,n;
  PC_Eisenstat   *eis = (PC_Eisenstat*)pc->data;

  PetscFunctionBegin;
  if (!pc->setupcalled) {
    ierr = MatGetSize(pc->mat,&M,&N);CHKERRQ(ierr);
    ierr = MatGetLocalSize(pc->mat,&m,&n);CHKERRQ(ierr);
    ierr = MatCreate(((PetscObject)pc)->comm,&eis->shell);CHKERRQ(ierr);
    ierr = MatSetSizes(eis->shell,m,n,M,N);CHKERRQ(ierr);
    ierr = MatSetType(eis->shell,MATSHELL);CHKERRQ(ierr);
    ierr = MatSetUp(eis->shell);CHKERRQ(ierr);
    ierr = MatShellSetContext(eis->shell,(void*)pc);CHKERRQ(ierr);
    ierr = PetscLogObjectParent(pc,eis->shell);CHKERRQ(ierr);
    ierr = MatShellSetOperation(eis->shell,MATOP_MULT,(void(*)(void))PCMult_Eisenstat);CHKERRQ(ierr);
  }
  if (!eis->usediag) PetscFunctionReturn(0);
  if (!pc->setupcalled) {
    ierr = MatGetVecs(pc->pmat,&eis->diag,0);CHKERRQ(ierr);
    ierr = PetscLogObjectParent(pc,eis->diag);CHKERRQ(ierr);
  }
  ierr = MatGetDiagonal(pc->pmat,eis->diag);CHKERRQ(ierr);
  PetscFunctionReturn(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);
}

void createOuterContext(OuterContext* & ctx) {
  int rank;
  int npes;
  MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
  MPI_Comm_size(PETSC_COMM_WORLD, &npes);

  ctx = new OuterContext;

  ctx->data = NULL;
  ctx->root = NULL;
  ctx->outerMat = PETSC_NULL;
  ctx->outerKsp = PETSC_NULL;
  ctx->outerPC = PETSC_NULL;
  ctx->outerSol = PETSC_NULL;
  ctx->outerRhs = PETSC_NULL;

  createLocalData(ctx->data);

  createRSDtree(ctx->root, rank, npes);

  createOuterMat(ctx);

  createOuterPC(ctx);

  createOuterKsp(ctx);

  MatGetVecs(ctx->outerMat, &(ctx->outerSol), &(ctx->outerRhs));
}

// updated
PetscErrorCode BSSCR_MatStokesMVBlockDefaultBuildScaling( MatStokesBlockScaling BA, Mat A, Vec b, Vec x, PetscTruth sym )
{
	Mat K,G,D,C;
	Vec rG;
	PetscScalar rg2, rg, ra;  
	PetscInt N;
	Vec rA, rC;
	Vec L1,L2, R1,R2;
	Mat S;
	
	VecNestGetSubVec( BA->Lz, 0, &L1 );
	VecNestGetSubVec( BA->Lz, 1, &L2 );
	
	VecNestGetSubVec( BA->Rz, 0, &R1 );
	VecNestGetSubVec( BA->Rz, 1, &R2 );
	
	rA = L1;
	rC = L2;
	
	MatNestGetSubMat( A, 0,0, &K );
	MatNestGetSubMat( A, 0,1, &G );
	MatNestGetSubMat( A, 1,0, &D );
	MatNestGetSubMat( A, 1,1, &C );
	
	VecDuplicate( rA, &rG );
	
	
	/* Get magnitude of K */  
	//px_MatGetAbsRowSum( K, rA );
        //MatGetRowMax( K, rA, PETSC_NULL );
	MatGetDiagonal( K, rA);

	VecSqrt( rA );  
	VecReciprocal( rA );
	
	/* VecDot( rA,rA, &ra ); */
	/* VecGetSize( rA, &N ); */
	/* ra = PetscSqrtScalar( ra/N ); */
	
	
	/* Get magnitude of G */
	//px_MatGetAbsRowSum( G, rG );
	//MatGetRowMax( G, rG, PETSC_NULL );

	Mat A21_cpy;
	Mat Shat;
	Vec diag; /* same as rA*rA */
	
	
	MatGetVecs( K, &diag, PETSC_NULL );
	MatGetDiagonal( K, diag );
	VecReciprocal( diag );

	if( sym )
#if( PETSC_VERSION_MAJOR <= 2 )
	    MatTranspose( G, &A21_cpy );
#else
	    MatTranspose( G, MAT_INITIAL_MATRIX, &A21_cpy );
#endif
	else {

    NonlocalCollection(Mat depends_on, IS interest_set) {
	find_influences(depends_on, interest_set, &nodes);

	PetscInt local_size;
	ISGetLocalSize(nodes, &local_size);
	VecCreateMPI(PETSC_COMM_WORLD, local_size, PETSC_DECIDE, &vec);
	IS onto_index_set;
	describe_partition(vec, &onto_index_set);
	PetscInt begin;
	PetscInt end;
	VecGetOwnershipRange(vec, &begin, &end);
	PetscInt *indicies;
	ISGetIndices(nodes, &indicies);
	assert(local_size == end-begin);
	for (int ii=0; ii<local_size; ii++) {
	    map[indicies[ii]] = ii+begin;
	}
	ISRestoreIndices(nodes, &indicies);
	Vec w;
	MatGetVecs(depends_on, PETSC_NULL, &w);
	VecScatterCreate(w, nodes, vec, onto_index_set, &scatter);
	VecDestroy(w);

	ISDestroy(onto_index_set);
    }

PetscErrorCode MatMultTranspose_Composite_Multiplicative(Mat A,Vec x,Vec y)
{
  Mat_Composite     *shell = (Mat_Composite*)A->data;
  Mat_CompositeLink tail   = shell->tail;
  PetscErrorCode    ierr;
  Vec               in,out;

  PetscFunctionBegin;
  if (!tail) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()");
  in = x;
  if (shell->left) {
    if (!shell->leftwork) {
      ierr = VecDuplicate(shell->left,&shell->leftwork);CHKERRQ(ierr);
    }
    ierr = VecPointwiseMult(shell->leftwork,shell->left,in);CHKERRQ(ierr);
    in   = shell->leftwork;
  }
  while (tail->prev) {
    if (!tail->prev->work) { /* should reuse previous work if the same size */
      ierr = MatGetVecs(tail->mat,NULL,&tail->prev->work);CHKERRQ(ierr);
    }
    out  = tail->prev->work;
    ierr = MatMultTranspose(tail->mat,in,out);CHKERRQ(ierr);
    in   = out;
    tail = tail->prev;
  }
  ierr = MatMultTranspose(tail->mat,in,y);CHKERRQ(ierr);
  if (shell->right) {
    ierr = VecPointwiseMult(y,shell->right,y);CHKERRQ(ierr);
  }
  ierr = VecScale(y,shell->scale);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

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

PetscErrorCode computeMaxEigVal(Mat A, PetscInt its, PetscScalar *eig)
{
  PetscErrorCode ierr;
  PetscRandom    rctx;      /* random number generator context */
  Vec            x0, x, x_1, tmp;
  PetscScalar    lambda_its, lambda_its_1;
  PetscInt       i;

  PetscFunctionBeginUser;
  ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rctx);CHKERRQ(ierr);
  ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr);
  ierr = MatGetVecs(A, &x_1, &x);CHKERRQ(ierr);
  ierr = VecSetRandom(x, rctx);CHKERRQ(ierr);
  ierr = VecDuplicate(x, &x0);CHKERRQ(ierr);
  ierr = VecCopy(x, x0);CHKERRQ(ierr);

  ierr = MatMult(A, x, x_1);CHKERRQ(ierr);
  for (i=0; i<its; i++) {
    tmp  = x; x = x_1; x_1 = tmp;
    ierr = MatMult(A, x, x_1);CHKERRQ(ierr);
  }
  ierr = VecDot(x0, x, &lambda_its);CHKERRQ(ierr);
  ierr = VecDot(x0, x_1, &lambda_its_1);CHKERRQ(ierr);

  *eig = lambda_its_1/lambda_its;

  ierr = VecDestroy(&x0);CHKERRQ(ierr);
  ierr = VecDestroy(&x);CHKERRQ(ierr);
  ierr = VecDestroy(&x_1);CHKERRQ(ierr);
  PetscFunctionReturn(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 = MatGetVecs(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);
}

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

  PetscFunctionBegin;
  if (!Na->work1) {ierr = MatGetVecs(Na->A,&Na->work1,PETSC_NULL);CHKERRQ(ierr);}
  if (!Na->work2) {ierr = MatGetVecs(Na->A,&Na->work2,PETSC_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);
}

/*
   PEPLinearExtract_Norm - Auxiliary routine that copies the solution of the
   linear eigenproblem to the PEP object. The eigenvector of the generalized
   problem is supposed to be
                               z = [  x  ]
                                   [ l*x ]
   If |l|<1.0, the eigenvector is taken from z(1:n), otherwise from z(n+1:2*n).
   Finally, x is normalized so that ||x||_2 = 1.
*/
static PetscErrorCode PEPLinearExtract_Norm(PEP pep,EPS eps)
{
  PetscErrorCode ierr;
  PetscInt       i,offset;
  PetscScalar    *px;
  Vec            xr,xi,w,vi;
#if !defined(PETSC_USE_COMPLEX)
  Vec            vi1;
#endif
  Mat            A;

  PetscFunctionBegin;
  ierr = EPSGetOperators(eps,&A,NULL);CHKERRQ(ierr);
  ierr = MatGetVecs(A,&xr,NULL);CHKERRQ(ierr);
  ierr = VecDuplicate(xr,&xi);CHKERRQ(ierr);
  ierr = VecCreateMPIWithArray(PetscObjectComm((PetscObject)pep),1,pep->nloc,pep->n,NULL,&w);CHKERRQ(ierr);
  for (i=0;i<pep->nconv;i++) {
    ierr = EPSGetEigenpair(eps,i,&pep->eigr[i],&pep->eigi[i],xr,xi);CHKERRQ(ierr);
    pep->eigr[i] *= pep->sfactor;
    pep->eigi[i] *= pep->sfactor;
    if (SlepcAbsEigenvalue(pep->eigr[i],pep->eigi[i])>1.0) offset = pep->nloc;
    else offset = 0;
#if !defined(PETSC_USE_COMPLEX)
    if (pep->eigi[i]>0.0) {   /* first eigenvalue of a complex conjugate pair */
      ierr = VecGetArray(xr,&px);CHKERRQ(ierr);
      ierr = VecPlaceArray(w,px+offset);CHKERRQ(ierr);
      ierr = BVInsertVec(pep->V,i,w);CHKERRQ(ierr);
      ierr = VecResetArray(w);CHKERRQ(ierr);
      ierr = VecRestoreArray(xr,&px);CHKERRQ(ierr);
      ierr = VecGetArray(xi,&px);CHKERRQ(ierr);
      ierr = VecPlaceArray(w,px+offset);CHKERRQ(ierr);
      ierr = BVInsertVec(pep->V,i+1,w);CHKERRQ(ierr);
      ierr = VecResetArray(w);CHKERRQ(ierr);
      ierr = VecRestoreArray(xi,&px);CHKERRQ(ierr);
      ierr = BVGetColumn(pep->V,i,&vi);CHKERRQ(ierr);
      ierr = BVGetColumn(pep->V,i+1,&vi1);CHKERRQ(ierr);
      ierr = SlepcVecNormalize(vi,vi1,PETSC_TRUE,NULL);CHKERRQ(ierr);
      ierr = BVRestoreColumn(pep->V,i,&vi);CHKERRQ(ierr);
      ierr = BVRestoreColumn(pep->V,i+1,&vi1);CHKERRQ(ierr);
    } else if (pep->eigi[i]==0.0)   /* real eigenvalue */
#endif
    {
      ierr = VecGetArray(xr,&px);CHKERRQ(ierr);
      ierr = VecPlaceArray(w,px+offset);CHKERRQ(ierr);
      ierr = BVInsertVec(pep->V,i,w);CHKERRQ(ierr);
      ierr = VecResetArray(w);CHKERRQ(ierr);
      ierr = VecRestoreArray(xr,&px);CHKERRQ(ierr);
      ierr = BVGetColumn(pep->V,i,&vi);CHKERRQ(ierr);
      ierr = SlepcVecNormalize(vi,NULL,PETSC_FALSE,NULL);CHKERRQ(ierr);
      ierr = BVRestoreColumn(pep->V,i,&vi);CHKERRQ(ierr);
    }
  }
  ierr = VecDestroy(&w);CHKERRQ(ierr);
  ierr = VecDestroy(&xr);CHKERRQ(ierr);
  ierr = VecDestroy(&xi);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

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

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