本文整理汇总了C++中VecDuplicate函数的典型用法代码示例。如果您正苦于以下问题：C++ VecDuplicate函数的具体用法？C++ VecDuplicate怎么用？C++ VecDuplicate使用的例子？那么, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了VecDuplicate函数的15个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: TaoSolve_NTR
static PetscErrorCode TaoSolve_NTR(Tao tao)
{
TAO_NTR *tr = (TAO_NTR *)tao->data;
PC pc;
KSPConvergedReason ksp_reason;
TaoConvergedReason reason;
PetscReal fmin, ftrial, prered, actred, kappa, sigma, beta;
PetscReal tau, tau_1, tau_2, tau_max, tau_min, max_radius;
PetscReal f, gnorm;
PetscReal delta;
PetscReal norm_d;
PetscErrorCode ierr;
PetscInt iter = 0;
PetscInt bfgsUpdates = 0;
PetscInt needH;
PetscInt i_max = 5;
PetscInt j_max = 1;
PetscInt i, j, N, n, its;
PetscFunctionBegin;
if (tao->XL || tao->XU || tao->ops->computebounds) {
ierr = PetscPrintf(((PetscObject)tao)->comm,"WARNING: Variable bounds have been set but will be ignored by ntr algorithm\n");CHKERRQ(ierr);
}
tao->trust = tao->trust0;
/* Modify the radius if it is too large or small */
tao->trust = PetscMax(tao->trust, tr->min_radius);
tao->trust = PetscMin(tao->trust, tr->max_radius);
if (NTR_PC_BFGS == tr->pc_type && !tr->M) {
ierr = VecGetLocalSize(tao->solution,&n);CHKERRQ(ierr);
ierr = VecGetSize(tao->solution,&N);CHKERRQ(ierr);
ierr = MatCreateLMVM(((PetscObject)tao)->comm,n,N,&tr->M);CHKERRQ(ierr);
ierr = MatLMVMAllocateVectors(tr->M,tao->solution);CHKERRQ(ierr);
}
/* Check convergence criteria */
ierr = TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);CHKERRQ(ierr);
ierr = VecNorm(tao->gradient,NORM_2,&gnorm);CHKERRQ(ierr);
if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
needH = 1;
ierr = TaoMonitor(tao, iter, f, gnorm, 0.0, 1.0, &reason);CHKERRQ(ierr);
if (reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0);
/* Create vectors for the limited memory preconditioner */
if ((NTR_PC_BFGS == tr->pc_type) &&
(BFGS_SCALE_BFGS != tr->bfgs_scale_type)) {
if (!tr->Diag) {
ierr = VecDuplicate(tao->solution, &tr->Diag);CHKERRQ(ierr);
}
}
switch(tr->ksp_type) {
case NTR_KSP_NASH:
ierr = KSPSetType(tao->ksp, KSPNASH);CHKERRQ(ierr);
if (tao->ksp->ops->setfromoptions) {
(*tao->ksp->ops->setfromoptions)(tao->ksp);
}
break;
case NTR_KSP_STCG:
ierr = KSPSetType(tao->ksp, KSPSTCG);CHKERRQ(ierr);
if (tao->ksp->ops->setfromoptions) {
(*tao->ksp->ops->setfromoptions)(tao->ksp);
}
break;
default:
ierr = KSPSetType(tao->ksp, KSPGLTR);CHKERRQ(ierr);
if (tao->ksp->ops->setfromoptions) {
(*tao->ksp->ops->setfromoptions)(tao->ksp);
}
break;
}
/* Modify the preconditioner to use the bfgs approximation */
ierr = KSPGetPC(tao->ksp, &pc);CHKERRQ(ierr);
switch(tr->pc_type) {
case NTR_PC_NONE:
ierr = PCSetType(pc, PCNONE);CHKERRQ(ierr);
if (pc->ops->setfromoptions) {
(*pc->ops->setfromoptions)(pc);
}
break;
case NTR_PC_AHESS:
ierr = PCSetType(pc, PCJACOBI);CHKERRQ(ierr);
if (pc->ops->setfromoptions) {
(*pc->ops->setfromoptions)(pc);
}
ierr = PCJacobiSetUseAbs(pc);CHKERRQ(ierr);
break;
case NTR_PC_BFGS:
ierr = PCSetType(pc, PCSHELL);CHKERRQ(ierr);
//.........这里部分代码省略.........
示例2: output
/**
* output
* ------
* Output the E and H fields.
*/
PetscErrorCode output(char *output_name, const Vec x, const Mat CF, const Vec conjParam, const Vec conjSrc, const GridInfo gi)
{
PetscFunctionBegin;
PetscErrorCode ierr;
char output_name_prefixed[PETSC_MAX_PATH_LEN];
//const char *prefix = "/out/";
const char *h_extension = ".H.h5";
const char *e_extension = ".E.h5";
//ierr = PetscStrcpy(output_name_prefixed, getenv("FD3D_ROOT")); CHKERRQ(ierr);
//ierr = PetscStrcat(output_name_prefixed, prefix); CHKERRQ(ierr);
//ierr = PetscStrcat(output_name_prefixed, output_name); CHKERRQ(ierr);
ierr = PetscStrcpy(output_name_prefixed, output_name); CHKERRQ(ierr);
char h_file[PETSC_MAX_PATH_LEN];
char e_file[PETSC_MAX_PATH_LEN];
ierr = PetscStrcpy(h_file, output_name_prefixed); CHKERRQ(ierr);
ierr = PetscStrcat(h_file, h_extension); CHKERRQ(ierr);
ierr = PetscStrcpy(e_file, output_name_prefixed); CHKERRQ(ierr);
ierr = PetscStrcat(e_file, e_extension); CHKERRQ(ierr);
Vec y; // H field vector if x_type == Etype
ierr = VecDuplicate(gi.vecTemp, &y); CHKERRQ(ierr);
ierr = VecCopy(conjSrc, y); CHKERRQ(ierr);
if (gi.x_type==Etype) {
ierr = MatMultAdd(CF, x, y, y); CHKERRQ(ierr);
ierr = VecScale(y, -1.0/PETSC_i/gi.omega); CHKERRQ(ierr);
} else {
ierr = VecScale(y, -1.0); CHKERRQ(ierr);
ierr = MatMultAdd(CF, x, y, y); CHKERRQ(ierr);
ierr = VecScale(y, 1.0/PETSC_i/gi.omega); CHKERRQ(ierr);
}
ierr = VecPointwiseDivide(y, y, conjParam);
PetscViewer viewer;
//viewer = PETSC_VIEWER_STDOUT_WORLD;
//ierr = PetscViewerHDF5Open(PETSC_COMM_WORLD, h_file, FILE_MODE_WRITE, &viewer); CHKERRQ(ierr);
/** Write the E-field file. */
ierr = PetscViewerHDF5Open(PETSC_COMM_WORLD, e_file, FILE_MODE_WRITE, &viewer); CHKERRQ(ierr);
ierr = PetscViewerHDF5PushGroup(viewer, "/");
if (gi.x_type==Etype) {
ierr = PetscObjectSetName((PetscObject) x, "E"); CHKERRQ(ierr);
ierr = VecView(x, viewer); CHKERRQ(ierr);
} else {
assert(gi.x_type==Htype);
ierr = PetscObjectSetName((PetscObject) y, "E"); CHKERRQ(ierr);
ierr = VecView(y, viewer); CHKERRQ(ierr);
}
ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
/** Write the H-field file. */
ierr = PetscViewerHDF5Open(PETSC_COMM_WORLD, h_file, FILE_MODE_WRITE, &viewer); CHKERRQ(ierr);
ierr = PetscViewerHDF5PushGroup(viewer, "/");
if (gi.x_type==Etype) {
ierr = PetscObjectSetName((PetscObject) y, "H"); CHKERRQ(ierr);
ierr = VecView(y, viewer); CHKERRQ(ierr);
} else {
assert(gi.x_type==Htype);
ierr = PetscObjectSetName((PetscObject) x, "H"); CHKERRQ(ierr);
ierr = VecView(x, viewer); CHKERRQ(ierr);
}
ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
ierr = VecDestroy(&y); CHKERRQ(ierr);
PetscFunctionReturn(0);
}示例3: main
PetscInt main(PetscInt argc,char **args)
{
PetscErrorCode ierr;
PetscMPIInt rank,size;
PetscInt N0=4096,N1=4096,N2=256,N3=10,N4=10,N=N0*N1;
PetscRandom rdm;
PetscReal enorm;
Vec x,y,z,input,output;
Mat A;
PetscInt DIM, dim[5],vsize,row,col;
PetscReal fac;
ierr = PetscInitialize(&argc,&args,(char*)0,help);CHKERRQ(ierr);
ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size);CHKERRQ(ierr);
ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);CHKERRQ(ierr);
#if defined(PETSC_USE_COMPLEX)
SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP, "Example for Real DFT. Your current data type is complex!");
#endif
ierr = PetscRandomCreate(PETSC_COMM_WORLD, &rdm);CHKERRQ(ierr);
ierr = PetscRandomSetFromOptions(rdm);CHKERRQ(ierr);
ierr = VecCreate(PETSC_COMM_WORLD,&input);CHKERRQ(ierr);
ierr = VecSetSizes(input,PETSC_DECIDE,N);CHKERRQ(ierr);
ierr = VecSetFromOptions(input);CHKERRQ(ierr);
ierr = VecSetRandom(input,rdm);CHKERRQ(ierr);
ierr = VecDuplicate(input,&output);
DIM = 2; dim[0] = N0; dim[1] = N1; dim[2] = N2; dim[3] = N3; dim[4] = N4;
ierr = MatCreateFFT(PETSC_COMM_WORLD,DIM,dim,MATFFTW,&A);CHKERRQ(ierr);
ierr = MatGetLocalSize(A,&row,&col);CHKERRQ(ierr);
printf("The Matrix size is %d and %d from process %d\n",row,col,rank);
ierr = MatCreateVecsFFTW(A,&x,&y,&z);CHKERRQ(ierr);
ierr = VecGetSize(x,&vsize);CHKERRQ(ierr);
ierr = VecGetSize(z,&vsize);CHKERRQ(ierr);
printf("The vector size of output from the main routine is %d\n",vsize);
ierr = VecScatterPetscToFFTW(A,input,x);CHKERRQ(ierr);
/*ierr = VecDestroy(&input);CHKERRQ(ierr);*/
ierr = MatMult(A,x,y);CHKERRQ(ierr);
ierr = MatMultTranspose(A,y,z);CHKERRQ(ierr);
ierr = VecScatterFFTWToPetsc(A,z,output);CHKERRQ(ierr);
/*ierr = VecDestroy(&z);CHKERRQ(ierr);*/
fac = 1.0/(PetscReal)N;
ierr = VecScale(output,fac);CHKERRQ(ierr);
ierr = VecAssemblyBegin(input);CHKERRQ(ierr);
ierr = VecAssemblyEnd(input);CHKERRQ(ierr);
ierr = VecAssemblyBegin(output);CHKERRQ(ierr);
ierr = VecAssemblyEnd(output);CHKERRQ(ierr);
/* ierr = VecView(input,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);*/
/* ierr = VecView(output,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);*/
ierr = VecAXPY(output,-1.0,input);CHKERRQ(ierr);
ierr = VecNorm(output,NORM_1,&enorm);CHKERRQ(ierr);
if (enorm > 1.e-14) {
ierr = PetscPrintf(PETSC_COMM_SELF," Error norm of |x - z| %e\n",enorm);CHKERRQ(ierr);
}
ierr = VecDestroy(&x);CHKERRQ(ierr);
ierr = VecDestroy(&y);CHKERRQ(ierr);
ierr = VecDestroy(&z);CHKERRQ(ierr);
ierr = VecDestroy(&output);CHKERRQ(ierr);
ierr = VecDestroy(&input);CHKERRQ(ierr);
ierr = MatDestroy(&A);CHKERRQ(ierr);
ierr = PetscRandomDestroy(&rdm);CHKERRQ(ierr);
PetscFinalize();
return 0;
}示例4: main
int main(int argc,char **args)
{
Mat *A,B; /* matrix */
PetscErrorCode ierr;
Vec x,y,v,v2,z;
PetscReal rnorm;
PetscInt n = 20; /* size of the matrix */
PetscInt nmat = 3; /* number of matrices */
PetscInt i;
PetscRandom rctx;
MatCompositeType type;
ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
ierr = PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(NULL,NULL,"-nmat",&nmat,NULL);CHKERRQ(ierr);
/*
Create random matrices
*/
ierr = PetscMalloc1(nmat+3,&A);CHKERRQ(ierr);
ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rctx);CHKERRQ(ierr);
ierr = MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n/2,3,NULL,3,NULL,&A[0]);CHKERRQ(ierr);
for (i = 1; i < nmat+1; i++) {
ierr = MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,3,NULL,3,NULL,&A[i]);CHKERRQ(ierr);
}
ierr = MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n/2,n,3,NULL,3,NULL,&A[nmat+1]);CHKERRQ(ierr);
for (i = 0; i < nmat+2; i++) {
ierr = MatSetRandom(A[i],rctx);CHKERRQ(ierr);
}
ierr = MatCreateVecs(A[1],&x,&y);CHKERRQ(ierr);
ierr = VecDuplicate(y,&z);CHKERRQ(ierr);
ierr = MatCreateVecs(A[0],&v,NULL);CHKERRQ(ierr);
ierr = VecDuplicate(v,&v2);CHKERRQ(ierr);
ierr = VecSet(x,1.0);CHKERRQ(ierr);
ierr = VecSet(y,0.0);CHKERRQ(ierr);
ierr = MatMult(A[1],x,z);CHKERRQ(ierr);
for (i = 2; i < nmat+1; i++) {
ierr = MatMultAdd(A[i],x,z,z);CHKERRQ(ierr);
}
ierr = MatCreateComposite(PETSC_COMM_WORLD,nmat,A+1,&B);CHKERRQ(ierr);
ierr = MatMultAdd(B,x,y,y);CHKERRQ(ierr);
ierr = VecAXPY(y,-1.0,z);CHKERRQ(ierr);
ierr = VecNorm(y,NORM_2,&rnorm);CHKERRQ(ierr);
if (rnorm > 10000.0*PETSC_MACHINE_EPSILON) {
ierr = PetscPrintf(PETSC_COMM_WORLD,"Error with composite add %g\n",(double)rnorm);CHKERRQ(ierr);
}
ierr = MatCompositeSetMatStructure(B,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); /* default */
ierr = MatCompositeMerge(B);CHKERRQ(ierr);
ierr = MatMult(B,x,y);CHKERRQ(ierr);
ierr = MatDestroy(&B);CHKERRQ(ierr);
ierr = VecAXPY(y,-1.0,z);CHKERRQ(ierr);
ierr = VecNorm(y,NORM_2,&rnorm);CHKERRQ(ierr);
if (rnorm > 10000.0*PETSC_MACHINE_EPSILON) {
ierr = PetscPrintf(PETSC_COMM_WORLD,"Error with composite add after merge %g\n",(double)rnorm);CHKERRQ(ierr);
}
/*
Test n x n/2 multiplicative composite
*/
ierr = VecSet(v,1.0);CHKERRQ(ierr);
ierr = MatMult(A[0],v,z);CHKERRQ(ierr);
for (i = 1; i < nmat; i++) {
ierr = MatMult(A[i],z,y);CHKERRQ(ierr);
ierr = VecCopy(y,z);CHKERRQ(ierr);
}
ierr = MatCreateComposite(PETSC_COMM_WORLD,nmat,A,&B);CHKERRQ(ierr);
ierr = MatCompositeSetType(B,MAT_COMPOSITE_MULTIPLICATIVE);CHKERRQ(ierr);
ierr = MatCompositeSetMergeType(B,MAT_COMPOSITE_MERGE_LEFT);CHKERRQ(ierr);
ierr = MatSetFromOptions(B);CHKERRQ(ierr);
ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); /* do MatCompositeMerge() if -mat_composite_merge 1 */
ierr = MatMult(B,v,y);CHKERRQ(ierr);
ierr = MatDestroy(&B);CHKERRQ(ierr);
ierr = VecAXPY(y,-1.0,z);CHKERRQ(ierr);
ierr = VecNorm(y,NORM_2,&rnorm);CHKERRQ(ierr);
if (rnorm > 10000.0*PETSC_MACHINE_EPSILON) {
ierr = PetscPrintf(PETSC_COMM_WORLD,"Error with composite multiplicative %g\n",(double)rnorm);CHKERRQ(ierr);
}
/*
Test n/2 x n multiplicative composite
*/
ierr = VecSet(x,1.0);CHKERRQ(ierr);
ierr = MatMult(A[2],x,z);CHKERRQ(ierr);
for (i = 3; i < nmat+1; i++) {
ierr = MatMult(A[i],z,y);CHKERRQ(ierr);
ierr = VecCopy(y,z);CHKERRQ(ierr);
}
ierr = MatMult(A[nmat+1],z,v);CHKERRQ(ierr);
ierr = MatCreateComposite(PETSC_COMM_WORLD,nmat,A+2,&B);CHKERRQ(ierr);
ierr = MatCompositeSetType(B,MAT_COMPOSITE_MULTIPLICATIVE);CHKERRQ(ierr);
ierr = MatSetFromOptions(B);CHKERRQ(ierr);
ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); /* do MatCompositeMerge() if -mat_composite_merge 1 */
//.........这里部分代码省略.........
示例5: main
int main(int argc,char **args)
{
Vec x,b,u; /* approx solution, RHS, exact solution */
Mat A; /* linear system matrix */
KSP ksp; /* KSP context */
PetscErrorCode ierr;
PetscInt i,n = 10,col[3],its,i1,i2;
PetscScalar none = -1.0,value[3],avalue;
PetscReal norm;
PC pc;
ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
ierr = PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);CHKERRQ(ierr);
/* Create vectors */
ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr);
ierr = VecSetSizes(x,PETSC_DECIDE,n);CHKERRQ(ierr);
ierr = VecSetFromOptions(x);CHKERRQ(ierr);
ierr = VecDuplicate(x,&b);CHKERRQ(ierr);
ierr = VecDuplicate(x,&u);CHKERRQ(ierr);
/* create a solution that is orthogonal to the constants */
ierr = VecGetOwnershipRange(u,&i1,&i2);CHKERRQ(ierr);
for (i=i1; i<i2; i++) {
avalue = i;
VecSetValues(u,1,&i,&avalue,INSERT_VALUES);
}
ierr = VecAssemblyBegin(u);CHKERRQ(ierr);
ierr = VecAssemblyEnd(u);CHKERRQ(ierr);
ierr = VecSum(u,&avalue);CHKERRQ(ierr);
avalue = -avalue/(PetscReal)n;
ierr = VecShift(u,avalue);CHKERRQ(ierr);
/* Create and assemble matrix */
ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr);
ierr = MatSetFromOptions(A);CHKERRQ(ierr);
value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
for (i=1; i<n-1; i++) {
col[0] = i-1; col[1] = i; col[2] = i+1;
ierr = MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);CHKERRQ(ierr);
}
i = n - 1; col[0] = n - 2; col[1] = n - 1; value[1] = 1.0;
ierr = MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr);
i = 0; col[0] = 0; col[1] = 1; value[0] = 1.0; value[1] = -1.0;
ierr = MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr);
ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatMult(A,u,b);CHKERRQ(ierr);
/* Create KSP context; set operators and options; solve linear system */
ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
/* Insure that preconditioner has same null space as matrix */
/* currently does not do anything */
ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
/* ierr = KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */
/* Check error */
ierr = VecAXPY(x,none,u);CHKERRQ(ierr);
ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr);
ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g, Iterations %D\n",(double)norm,its);CHKERRQ(ierr);
/* Free work space */
ierr = VecDestroy(&x);CHKERRQ(ierr);ierr = VecDestroy(&u);CHKERRQ(ierr);
ierr = VecDestroy(&b);CHKERRQ(ierr);ierr = MatDestroy(&A);CHKERRQ(ierr);
ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
ierr = PetscFinalize();
return ierr;
}示例6: TaoSetup_GPCG
static PetscErrorCode TaoSetup_GPCG(Tao tao)
{
PetscErrorCode ierr;
TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
PetscFunctionBegin;
/* Allocate some arrays */
if (!tao->gradient) {
ierr = VecDuplicate(tao->solution, &tao->gradient);CHKERRQ(ierr);
}
if (!tao->stepdirection) {
ierr = VecDuplicate(tao->solution, &tao->stepdirection);CHKERRQ(ierr);
}
if (!tao->XL) {
ierr = VecDuplicate(tao->solution,&tao->XL);CHKERRQ(ierr);
ierr = VecSet(tao->XL,PETSC_NINFINITY);CHKERRQ(ierr);
}
if (!tao->XU) {
ierr = VecDuplicate(tao->solution,&tao->XU);CHKERRQ(ierr);
ierr = VecSet(tao->XU,PETSC_INFINITY);CHKERRQ(ierr);
}
ierr=VecDuplicate(tao->solution,&gpcg->B);CHKERRQ(ierr);
ierr=VecDuplicate(tao->solution,&gpcg->Work);CHKERRQ(ierr);
ierr=VecDuplicate(tao->solution,&gpcg->X_New);CHKERRQ(ierr);
ierr=VecDuplicate(tao->solution,&gpcg->G_New);CHKERRQ(ierr);
ierr=VecDuplicate(tao->solution,&gpcg->DXFree);CHKERRQ(ierr);
ierr=VecDuplicate(tao->solution,&gpcg->R);CHKERRQ(ierr);
ierr=VecDuplicate(tao->solution,&gpcg->PG);CHKERRQ(ierr);
/*
if (gpcg->ksp_type == GPCG_KSP_NASH) {
ierr = KSPSetType(tao->ksp,KSPNASH);CHKERRQ(ierr);
} else if (gpcg->ksp_type == GPCG_KSP_STCG) {
ierr = KSPSetType(tao->ksp,KSPSTCG);CHKERRQ(ierr);
} else {
ierr = KSPSetType(tao->ksp,KSPGLTR);CHKERRQ(ierr);
}
if (tao->ksp->ops->setfromoptions) {
(*tao->ksp->ops->setfromoptions)(tao->ksp);
}
}
*/
PetscFunctionReturn(0);
}示例7: PCBDDCNullSpaceAssembleCorrection
PetscErrorCode PCBDDCNullSpaceAssembleCorrection(PC pc,IS local_dofs)
{
PC_BDDC *pcbddc = (PC_BDDC*)pc->data;
PC_IS *pcis = (PC_IS*)pc->data;
Mat_IS* matis = (Mat_IS*)pc->pmat->data;
KSP *local_ksp;
PC newpc;
NullSpaceCorrection_ctx shell_ctx;
Mat local_mat,local_pmat,small_mat,inv_small_mat;
Vec work1,work2;
const Vec *nullvecs;
VecScatter scatter_ctx;
IS is_aux;
MatFactorInfo matinfo;
PetscScalar *basis_mat,*Kbasis_mat,*array,*array_mat;
PetscScalar one = 1.0,zero = 0.0, m_one = -1.0;
PetscInt basis_dofs,basis_size,nnsp_size,i,k,n_I,n_R;
PetscBool nnsp_has_cnst;
PetscErrorCode ierr;
PetscFunctionBegin;
/* Infer the local solver */
ierr = ISGetSize(local_dofs,&basis_dofs);CHKERRQ(ierr);
ierr = VecGetSize(pcis->vec1_D,&n_I);CHKERRQ(ierr);
ierr = VecGetSize(pcbddc->vec1_R,&n_R);CHKERRQ(ierr);
if (basis_dofs == n_I) {
/* Dirichlet solver */
local_ksp = &pcbddc->ksp_D;
} else if (basis_dofs == n_R) {
/* Neumann solver */
local_ksp = &pcbddc->ksp_R;
} else {
SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in %s: unknown local IS size %d. n_I=%d, n_R=%d)\n",__FUNCT__,basis_dofs,n_I,n_R);
}
ierr = KSPGetOperators(*local_ksp,&local_mat,&local_pmat);CHKERRQ(ierr);
/* Get null space vecs */
ierr = MatNullSpaceGetVecs(pcbddc->NullSpace,&nnsp_has_cnst,&nnsp_size,&nullvecs);CHKERRQ(ierr);
basis_size = nnsp_size;
if (nnsp_has_cnst) {
basis_size++;
}
if (basis_dofs) {
/* Create shell ctx */
ierr = PetscMalloc(sizeof(*shell_ctx),&shell_ctx);CHKERRQ(ierr);
/* Create work vectors in shell context */
ierr = VecCreate(PETSC_COMM_SELF,&shell_ctx->work_small_1);CHKERRQ(ierr);
ierr = VecSetSizes(shell_ctx->work_small_1,basis_size,basis_size);CHKERRQ(ierr);
ierr = VecSetType(shell_ctx->work_small_1,VECSEQ);CHKERRQ(ierr);
ierr = VecDuplicate(shell_ctx->work_small_1,&shell_ctx->work_small_2);CHKERRQ(ierr);
ierr = VecCreate(PETSC_COMM_SELF,&shell_ctx->work_full_1);CHKERRQ(ierr);
ierr = VecSetSizes(shell_ctx->work_full_1,basis_dofs,basis_dofs);CHKERRQ(ierr);
ierr = VecSetType(shell_ctx->work_full_1,VECSEQ);CHKERRQ(ierr);
ierr = VecDuplicate(shell_ctx->work_full_1,&shell_ctx->work_full_2);CHKERRQ(ierr);
/* Allocate workspace */
ierr = MatCreateSeqDense(PETSC_COMM_SELF,basis_dofs,basis_size,NULL,&shell_ctx->basis_mat );CHKERRQ(ierr);
ierr = MatCreateSeqDense(PETSC_COMM_SELF,basis_dofs,basis_size,NULL,&shell_ctx->Kbasis_mat);CHKERRQ(ierr);
ierr = MatDenseGetArray(shell_ctx->basis_mat,&basis_mat);CHKERRQ(ierr);
ierr = MatDenseGetArray(shell_ctx->Kbasis_mat,&Kbasis_mat);CHKERRQ(ierr);
/* Restrict local null space on selected dofs (Dirichlet or Neumann)
and compute matrices N and K*N */
ierr = VecDuplicate(shell_ctx->work_full_1,&work1);CHKERRQ(ierr);
ierr = VecDuplicate(shell_ctx->work_full_1,&work2);CHKERRQ(ierr);
ierr = VecScatterCreate(pcis->vec1_N,local_dofs,work1,(IS)0,&scatter_ctx);CHKERRQ(ierr);
}
for (k=0;k<nnsp_size;k++) {
ierr = VecScatterBegin(matis->ctx,nullvecs[k],pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
ierr = VecScatterEnd(matis->ctx,nullvecs[k],pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
if (basis_dofs) {
ierr = VecPlaceArray(work1,(const PetscScalar*)&basis_mat[k*basis_dofs]);CHKERRQ(ierr);
ierr = VecScatterBegin(scatter_ctx,pcis->vec1_N,work1,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
ierr = VecScatterEnd(scatter_ctx,pcis->vec1_N,work1,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
ierr = VecPlaceArray(work2,(const PetscScalar*)&Kbasis_mat[k*basis_dofs]);CHKERRQ(ierr);
ierr = MatMult(local_mat,work1,work2);CHKERRQ(ierr);
ierr = VecResetArray(work1);CHKERRQ(ierr);
ierr = VecResetArray(work2);CHKERRQ(ierr);
}
}
if (basis_dofs) {
if (nnsp_has_cnst) {
ierr = VecPlaceArray(work1,(const PetscScalar*)&basis_mat[k*basis_dofs]);CHKERRQ(ierr);
ierr = VecSet(work1,one);CHKERRQ(ierr);
ierr = VecPlaceArray(work2,(const PetscScalar*)&Kbasis_mat[k*basis_dofs]);CHKERRQ(ierr);
ierr = MatMult(local_mat,work1,work2);CHKERRQ(ierr);
ierr = VecResetArray(work1);CHKERRQ(ierr);
ierr = VecResetArray(work2);CHKERRQ(ierr);
}
ierr = VecDestroy(&work1);CHKERRQ(ierr);
ierr = VecDestroy(&work2);CHKERRQ(ierr);
ierr = VecScatterDestroy(&scatter_ctx);CHKERRQ(ierr);
ierr = MatDenseRestoreArray(shell_ctx->basis_mat,&basis_mat);CHKERRQ(ierr);
ierr = MatDenseRestoreArray(shell_ctx->Kbasis_mat,&Kbasis_mat);CHKERRQ(ierr);
/* Assemble another Mat object in shell context */
//.........这里部分代码省略.........
示例8: main
int main(int argc, char **argv)
{
DM dm; /* Problem specification */
DM dmAux; /* Material specification */
SNES snes; /* nonlinear solver */
Vec u,r; /* solution, residual vectors */
Mat A,J; /* Jacobian matrix */
MatNullSpace nullSpace; /* May be necessary for Neumann conditions */
AppCtx user; /* user-defined work context */
JacActionCtx userJ; /* context for Jacobian MF action */
PetscInt its; /* iterations for convergence */
PetscReal error = 0.0; /* L_2 error in the solution */
PetscInt numComponents;
PetscErrorCode ierr;
ierr = PetscInitialize(&argc, &argv, NULL, help);CHKERRQ(ierr);
ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr);
ierr = SNESCreate(PETSC_COMM_WORLD, &snes);CHKERRQ(ierr);
ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr);
ierr = SNESSetDM(snes, dm);CHKERRQ(ierr);
ierr = DMSetApplicationContext(dm, &user);CHKERRQ(ierr);
user.fem.fe = user.fe;
ierr = PetscFECreateDefault(dm, user.dim, 1, PETSC_TRUE, NULL, -1, &user.fe[0]);CHKERRQ(ierr);
if (user.bcType != NEUMANN) {
user.fem.feBd = NULL; user.feBd[0] = NULL;
} else {
user.fem.feBd = user.feBd;
ierr = PetscFECreateDefault(dm, user.dim-1, 1, PETSC_TRUE, "bd_", -1, &user.feBd[0]);CHKERRQ(ierr);
}
ierr = PetscFEGetNumComponents(user.fe[0], &numComponents);CHKERRQ(ierr);
ierr = PetscMalloc(NUM_FIELDS * sizeof(void (*)(const PetscReal[], PetscScalar *, void *)), &user.exactFuncs);CHKERRQ(ierr);
ierr = SetupProblem(dm, &user);CHKERRQ(ierr);
ierr = SetupSection(dm, &user);CHKERRQ(ierr);
/* Setup Material */
ierr = DMClone(dm, &dmAux);CHKERRQ(ierr);
ierr = DMPlexCopyCoordinates(dm, dmAux);CHKERRQ(ierr);
if (user.variableCoefficient != COEFF_FIELD) {
user.fem.feAux = NULL; user.feAux[0] = NULL;
} else {
PetscQuadrature q;
user.fem.feAux = user.feAux;
ierr = PetscFECreateDefault(dmAux, user.dim, 1, PETSC_TRUE, "mat_", -1, &user.feAux[0]);CHKERRQ(ierr);
ierr = PetscFEGetQuadrature(user.fe[0], &q);CHKERRQ(ierr);
ierr = PetscFESetQuadrature(user.feAux[0], q);CHKERRQ(ierr);
}
ierr = SetupMaterialSection(dmAux, &user);CHKERRQ(ierr);
ierr = SetupMaterial(dm, dmAux, &user);CHKERRQ(ierr);
ierr = DMDestroy(&dmAux);CHKERRQ(ierr);
ierr = DMCreateGlobalVector(dm, &u);CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject) u, "potential");CHKERRQ(ierr);
ierr = VecDuplicate(u, &r);CHKERRQ(ierr);
ierr = DMSetMatType(dm,MATAIJ);CHKERRQ(ierr);
ierr = DMCreateMatrix(dm, &J);CHKERRQ(ierr);
if (user.jacobianMF) {
PetscInt M, m, N, n;
ierr = MatGetSize(J, &M, &N);CHKERRQ(ierr);
ierr = MatGetLocalSize(J, &m, &n);CHKERRQ(ierr);
ierr = MatCreate(PETSC_COMM_WORLD, &A);CHKERRQ(ierr);
ierr = MatSetSizes(A, m, n, M, N);CHKERRQ(ierr);
ierr = MatSetType(A, MATSHELL);CHKERRQ(ierr);
ierr = MatSetUp(A);CHKERRQ(ierr);
#if 0
ierr = MatShellSetOperation(A, MATOP_MULT, (void (*)(void))FormJacobianAction);CHKERRQ(ierr);
#endif
userJ.dm = dm;
userJ.J = J;
userJ.user = &user;
ierr = DMCreateLocalVector(dm, &userJ.u);CHKERRQ(ierr);
ierr = DMPlexProjectFunctionLocal(dm, user.fe, user.exactFuncs, NULL, INSERT_BC_VALUES, userJ.u);CHKERRQ(ierr);
ierr = MatShellSetContext(A, &userJ);CHKERRQ(ierr);
} else {
A = J;
}
if (user.bcType == NEUMANN) {
ierr = MatNullSpaceCreate(PetscObjectComm((PetscObject) dm), PETSC_TRUE, 0, NULL, &nullSpace);CHKERRQ(ierr);
ierr = MatSetNullSpace(J, nullSpace);CHKERRQ(ierr);
if (A != J) {
ierr = MatSetNullSpace(A, nullSpace);CHKERRQ(ierr);
}
}
ierr = DMSNESSetFunctionLocal(dm, (PetscErrorCode (*)(DM,Vec,Vec,void*)) DMPlexComputeResidualFEM, &user);CHKERRQ(ierr);
ierr = DMSNESSetJacobianLocal(dm, (PetscErrorCode (*)(DM,Vec,Mat,Mat,void*)) DMPlexComputeJacobianFEM, &user);CHKERRQ(ierr);
ierr = SNESSetJacobian(snes, A, J, NULL, NULL);CHKERRQ(ierr);
ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);
ierr = DMPlexProjectFunction(dm, user.fe, user.exactFuncs, NULL, INSERT_ALL_VALUES, u);CHKERRQ(ierr);
if (user.checkpoint) {
#if defined(PETSC_HAVE_HDF5)
ierr = PetscViewerHDF5PushGroup(user.checkpoint, "/fields");CHKERRQ(ierr);
ierr = VecLoad(u, user.checkpoint);CHKERRQ(ierr);
ierr = PetscViewerHDF5PopGroup(user.checkpoint);CHKERRQ(ierr);
//.........这里部分代码省略.........
示例9: main
//.........这里部分代码省略.........
pvfmm::Profile::Toc();
pvfmm::Profile::Tic("FMMCreateShell",&comm,true);
Mat A;
{ // Create Matrix. A
FMMCreateShell(&fmm_data, &A);
MatShellSetOperation(A,MATOP_MULT,(void(*)(void))mult);
}
pvfmm::Profile::Toc();
pvfmm::Profile::Tic("Phi",&comm,true);
{ // Compute phi(x) = phi_0(x) -\int G(x,y)eta(y)phi_0(y)dy
CompPhiUsingBorn(phi, fmm_data);
fmm_data.tree->Write2File("results/phi",0);
// tree2vec(fmm_data,b);
}
pvfmm::Profile::Toc();
pvfmm::Profile::Tic("Right_hand_side",&comm,true);
{ // Compute phi_0(x) - phi(x)
// After this block phi becomes the RHS
ierr = VecAXPY(phi,-1,phi_0);
CHKERRQ(ierr);
vec2tree(phi,fmm_data);
//fmm_data.tree->RunFMM();
//fmm_data.tree->Copy_FMMOutput();
fmm_data.tree->Write2File("results/rhs",0);
// tree2vec(fmm_data,b);
}
pvfmm::Profile::Toc();
// Create solution vector
pvfmm::Profile::Tic("Initial_Vector_eta_comp",&comm,true);
ierr = VecDuplicate(pt_sources,&eta_comp);CHKERRQ(ierr);
pvfmm::Profile::Toc();
// Create linear solver context
pvfmm::Profile::Tic("KSPCreate",&comm,true);
KSP ksp ; ierr = KSPCreate(PETSC_COMM_WORLD,&ksp );CHKERRQ(ierr);
pvfmm::Profile::Toc();
// Set operators. Here the matrix that defines the linear system
// also serves as the preconditioning matrix.
pvfmm::Profile::Tic("KSPSetOperators",&comm,true);
ierr = KSPSetOperators(ksp ,A ,A ,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
pvfmm::Profile::Toc();
// Set runtime options
KSPSetType(ksp ,KSPGMRES);
KSPSetNormType(ksp , KSP_NORM_UNPRECONDITIONED);
KSPSetTolerances(ksp ,GMRES_TOL ,PETSC_DEFAULT,PETSC_DEFAULT,MAX_ITER );
//KSPGMRESSetRestart(ksp , MAX_ITER );
KSPGMRESSetRestart(ksp , 100 );
ierr = KSPSetFromOptions(ksp );CHKERRQ(ierr);
// -------------------------------------------------------------------
// Solve the linear system
// -------------------------------------------------------------------
pvfmm::Profile::Tic("KSPSolve",&comm,true);
time_ksp=-omp_get_wtime();
ierr = KSPSolve(ksp,phi,eta_comp);CHKERRQ(ierr);
MPI_Barrier(comm);
time_ksp+=omp_get_wtime();
pvfmm::Profile::Toc();
// View info about the solver
KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
// -------------------------------------------------------------------
// Check solution and clean up
// -------------------------------------------------------------------
// Iterations
PetscInt its;
ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"Iterations %D\n",its);CHKERRQ(ierr);
iter_ksp=its;
{ // Write output
vec2tree(eta_comp, fmm_data);
fmm_data.tree->Write2File("results/eta_comp",VTK_ORDER);
}
// Free work space. All PETSc objects should be destroyed when they
// are no longer needed.
ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
ierr = VecDestroy(&eta_comp);CHKERRQ(ierr);
ierr = VecDestroy(&phi_0);CHKERRQ(ierr);
ierr = MatDestroy(&A);CHKERRQ(ierr);
// Delete fmm data
FMMDestroy(&fmm_data);
pvfmm::Profile::print(&comm);
}
ierr = PetscFinalize();
return 0;
}示例10: mult
//.........这里部分代码省略.........
ierr = VecGetLocalSize(U, &U_size);
ierr = VecGetLocalSize(phi_0_vec, &phi_0_size);
int data_dof=U_size/(n_coeff3*nlist.size());
assert(data_dof*n_coeff3*nlist.size()==U_size);
PetscScalar *U_ptr;
PetscScalar* phi_0_ptr;
ierr = VecGetArray(U, &U_ptr);
ierr = VecGetArray(phi_0_vec, &phi_0_ptr);
#pragma omp parallel for
for(size_t tid=0;tid<omp_p;tid++){
size_t i_start=(nlist.size()* tid )/omp_p;
size_t i_end =(nlist.size()*(tid+1))/omp_p;
pvfmm::Vector<double> coeff_vec(n_coeff3*data_dof);
pvfmm::Vector<double> val_vec(n_nodes3*data_dof);
pvfmm::Vector<double> phi_0_part(n_nodes3*data_dof);
for(size_t i=i_start;i<i_end;i++){
double s=std::pow(2.0,COORD_DIM*nlist[i]->Depth()*0.5*SCAL_EXP);
{ // coeff_vec: Cheb coeff data for this node
size_t U_offset=i*n_coeff3*data_dof;
size_t phi_0_offset = i*n_nodes3*data_dof;
for(size_t j=0;j<n_coeff3*data_dof;j++){
coeff_vec[j]=PetscRealPart(U_ptr[j+U_offset])*s;
}
for(size_t j=0;j<n_nodes3*data_dof;j++){
phi_0_part[j]=PetscRealPart(phi_0_ptr[j+phi_0_offset]);
}
}
// val_vec: Evaluate coeff_vec at Chebyshev node points
cheb_eval(coeff_vec, cheb_deg, cheb_node_coord1, cheb_node_coord1, cheb_node_coord1, val_vec);
{// phi_0_part*val_vec
for(size_t j0=0;j0<data_dof;j0++){
double* vec=&val_vec[j0*n_nodes3];
double* phi_0=&phi_0_part[j0*n_nodes3];
for(size_t j1=0;j1<n_nodes3;j1++) vec[j1]*=phi_0[j1];
}
}
{ // Compute Chebyshev approx
pvfmm::Vector<double>& coeff_vec=nlist[i]->ChebData();
if(coeff_vec.Dim()!=(data_dof*(cheb_deg+1)*(cheb_deg+2)*(cheb_deg+3))/6){
coeff_vec.ReInit((data_dof*(cheb_deg+1)*(cheb_deg+2)*(cheb_deg+3))/6);
}
pvfmm::cheb_approx<double,double>(&val_vec[0], cheb_deg, data_dof, &coeff_vec[0]);
nlist[i]->DataDOF()=data_dof;
}
}
}
}
pvfmm::Profile::Toc();
// Run FMM ( Compute: G[ \eta * u ] )
tree->ClearFMMData();
tree->RunFMM();
// Copy data from tree to Y
pvfmm::Profile::Tic("tree2vec",comm,true);
{
PetscInt Y_size;
ierr = VecGetLocalSize(Y, &Y_size);
int data_dof=Y_size/(n_coeff3*nlist.size());
PetscScalar *Y_ptr;
ierr = VecGetArray(Y, &Y_ptr);
#pragma omp parallel for
for(size_t tid=0;tid<omp_p;tid++){
size_t i_start=(nlist.size()* tid )/omp_p;
size_t i_end =(nlist.size()*(tid+1))/omp_p;
for(size_t i=i_start;i<i_end;i++){
pvfmm::Vector<double>& coeff_vec=((FMM_Mat_t::FMMData*) nlist[i]->FMMData())->cheb_out;
double s=std::pow(0.5,COORD_DIM*nlist[i]->Depth()*0.5*SCAL_EXP);
size_t Y_offset=i*n_coeff3*data_dof;
for(size_t j=0;j<n_coeff3*data_dof;j++) Y_ptr[j+Y_offset]=coeff_vec[j]*s;
}
}
ierr = VecRestoreArray(Y, &Y_ptr);
}
pvfmm::Profile::Toc();
// Regularize
Vec alpha;
PetscScalar sca = (PetscScalar).00001;
VecDuplicate(Y,&alpha);
VecSet(alpha,sca);
ierr = VecPointwiseMult(alpha,alpha,U);
ierr = VecAXPY(Y,1,alpha);
// Output Vector ( Compute: U + G[ \eta * U ] )
pvfmm::Profile::Tic("FMM_Output",comm,true);
ierr = VecAXPY(Y,1,U);CHKERRQ(ierr);
pvfmm::Profile::Toc();
ierr = VecDestroy(&alpha); CHKERRQ(ierr);
pvfmm::Profile::Toc();
return 0;
}示例11: SNESSetUp_FAS
//.........这里部分代码省略.........
}
/* set the injection from the DM */
if (!fas->inject) {
ierr = DMCreateInjection(next->dm, snes->dm, &fas->inject);CHKERRQ(ierr);
}
}
}
/*pass the smoother, function, and jacobian up to the next level if it's not user set already */
if (fas->galerkin) {
if (next) {
ierr = SNESSetFunction(next, NULL, SNESFASGalerkinDefaultFunction, next);CHKERRQ(ierr);
}
if (fas->smoothd && fas->level != fas->levels - 1) {
ierr = SNESSetFunction(fas->smoothd, NULL, SNESFASGalerkinDefaultFunction, snes);CHKERRQ(ierr);
}
if (fas->smoothu && fas->level != fas->levels - 1) {
ierr = SNESSetFunction(fas->smoothu, NULL, SNESFASGalerkinDefaultFunction, snes);CHKERRQ(ierr);
}
}
/* sets the down (pre) smoother's default norm and sets it from options */
if (fas->smoothd) {
if (fas->level == 0 && fas->levels != 1) {
ierr = SNESSetNormSchedule(fas->smoothd, SNES_NORM_NONE);CHKERRQ(ierr);
} else {
ierr = SNESSetNormSchedule(fas->smoothd, SNES_NORM_FINAL_ONLY);CHKERRQ(ierr);
}
ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)fas->smoothd);CHKERRQ(ierr);
ierr = SNESSetFromOptions(fas->smoothd);CHKERRQ(ierr);
ierr = SNESGetLineSearch(snes,&linesearch);CHKERRQ(ierr);
ierr = SNESGetLineSearch(fas->smoothd,&slinesearch);CHKERRQ(ierr);
ierr = SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);CHKERRQ(ierr);
ierr = SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);CHKERRQ(ierr);
ierr = SNESLineSearchSetPreCheck(slinesearch,precheck,lsprectx);CHKERRQ(ierr);
ierr = SNESLineSearchSetPostCheck(slinesearch,postcheck,lspostctx);CHKERRQ(ierr);
ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)slinesearch);CHKERRQ(ierr);
fas->smoothd->vec_sol = snes->vec_sol;
ierr = PetscObjectReference((PetscObject)snes->vec_sol);CHKERRQ(ierr);
fas->smoothd->vec_sol_update = snes->vec_sol_update;
ierr = PetscObjectReference((PetscObject)snes->vec_sol_update);CHKERRQ(ierr);
fas->smoothd->vec_func = snes->vec_func;
ierr = PetscObjectReference((PetscObject)snes->vec_func);CHKERRQ(ierr);
if (fas->eventsmoothsetup) {ierr = PetscLogEventBegin(fas->eventsmoothsetup,0,0,0,0);CHKERRQ(ierr);}
ierr = SNESSetUp(fas->smoothd);CHKERRQ(ierr);
if (fas->eventsmoothsetup) {ierr = PetscLogEventEnd(fas->eventsmoothsetup,0,0,0,0);CHKERRQ(ierr);}
}
/* sets the up (post) smoother's default norm and sets it from options */
if (fas->smoothu) {
if (fas->level != fas->levels - 1) {
ierr = SNESSetNormSchedule(fas->smoothu, SNES_NORM_NONE);CHKERRQ(ierr);
} else {
ierr = SNESSetNormSchedule(fas->smoothu, SNES_NORM_FINAL_ONLY);CHKERRQ(ierr);
}
ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)fas->smoothu);CHKERRQ(ierr);
ierr = SNESSetFromOptions(fas->smoothu);CHKERRQ(ierr);
ierr = SNESGetLineSearch(snes,&linesearch);CHKERRQ(ierr);
ierr = SNESGetLineSearch(fas->smoothu,&slinesearch);CHKERRQ(ierr);
ierr = SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);CHKERRQ(ierr);
ierr = SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);CHKERRQ(ierr);
ierr = SNESLineSearchSetPreCheck(slinesearch,precheck,lsprectx);CHKERRQ(ierr);
ierr = SNESLineSearchSetPostCheck(slinesearch,postcheck,lspostctx);CHKERRQ(ierr);
ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)slinesearch);CHKERRQ(ierr);
fas->smoothu->vec_sol = snes->vec_sol;
ierr = PetscObjectReference((PetscObject)snes->vec_sol);CHKERRQ(ierr);
fas->smoothu->vec_sol_update = snes->vec_sol_update;
ierr = PetscObjectReference((PetscObject)snes->vec_sol_update);CHKERRQ(ierr);
fas->smoothu->vec_func = snes->vec_func;
ierr = PetscObjectReference((PetscObject)snes->vec_func);CHKERRQ(ierr);
if (fas->eventsmoothsetup) {ierr = PetscLogEventBegin(fas->eventsmoothsetup,0,0,0,0);CHKERRQ(ierr);}
ierr = SNESSetUp(fas->smoothu);CHKERRQ(ierr);
if (fas->eventsmoothsetup) {ierr = PetscLogEventEnd(fas->eventsmoothsetup,0,0,0,0);CHKERRQ(ierr);}
}
if (next) {
/* gotta set up the solution vector for this to work */
if (!next->vec_sol) {ierr = SNESFASCreateCoarseVec(snes,&next->vec_sol);CHKERRQ(ierr);}
if (!next->vec_rhs) {ierr = SNESFASCreateCoarseVec(snes,&next->vec_rhs);CHKERRQ(ierr);}
ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)next);CHKERRQ(ierr);
ierr = SNESGetLineSearch(snes,&linesearch);CHKERRQ(ierr);
ierr = SNESGetLineSearch(fas->next,&slinesearch);CHKERRQ(ierr);
ierr = SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);CHKERRQ(ierr);
ierr = SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);CHKERRQ(ierr);
ierr = SNESLineSearchSetPreCheck(slinesearch,precheck,lsprectx);CHKERRQ(ierr);
ierr = SNESLineSearchSetPostCheck(slinesearch,postcheck,lspostctx);CHKERRQ(ierr);
ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)slinesearch);CHKERRQ(ierr);
ierr = SNESSetUp(next);CHKERRQ(ierr);
}
/* setup FAS work vectors */
if (fas->galerkin) {
ierr = VecDuplicate(snes->vec_sol, &fas->Xg);CHKERRQ(ierr);
ierr = VecDuplicate(snes->vec_sol, &fas->Fg);CHKERRQ(ierr);
}
PetscFunctionReturn(0);
}示例12: SetupIterativeSolver
void VentilationProblem::SolveIterativelyFromPressure()
{
if (mTerminalInteractionMatrix == NULL)
{
SetupIterativeSolver();
}
// double start = Timer::GetElapsedTime();
/* Now use the pressure boundary conditions to determine suitable flux boundary conditions
* and iteratively update them until we are done
*/
assert(mPressure[mOutletNodeIndex] == mPressureCondition[mOutletNodeIndex]);
unsigned max_iterations=1000;
unsigned num_terminals = mMesh.GetNumBoundaryNodes()-1u;
double pressure_tolerance = 1e-4;
if (mLengthScaling < 1e-2)
{
// Using SI units
pressure_tolerance = 1e-7;
}
bool converged=false;
double last_norm_pressure_change;
Vec old_terminal_pressure_change;
VecDuplicate(mTerminalPressureChangeVector, &old_terminal_pressure_change);
for (unsigned iteration = 0; iteration < max_iterations && converged==false; iteration++)
{
for (unsigned terminal=0; terminal<num_terminals; terminal++)
{
unsigned node_index = mTerminalToNodeIndex[terminal];
// How far we are away from matching this boundary condition.
double delta_pressure = mPressure[node_index] - mPressureCondition[node_index];
// Offset the first iteration
if (iteration == 0)
{
delta_pressure += mPressureCondition[mOutletNodeIndex];
}
VecSetValue(mTerminalPressureChangeVector, terminal, delta_pressure, INSERT_VALUES);
}
double norm_pressure_change;
VecNorm(mTerminalPressureChangeVector, NORM_2, &last_norm_pressure_change);
if (last_norm_pressure_change < pressure_tolerance)
{
converged = true;
break;
}
VecCopy(mTerminalPressureChangeVector, old_terminal_pressure_change);
KSPSolve(mTerminalKspSolver, mTerminalPressureChangeVector, mTerminalFluxChangeVector);
double* p_terminal_flux_change_vector;
VecGetArray(mTerminalFluxChangeVector, &p_terminal_flux_change_vector);
for (unsigned terminal=0; terminal<num_terminals; terminal++)
{
double estimated_terminal_flux_change=p_terminal_flux_change_vector[terminal];
unsigned edge_index = mTerminalToEdgeIndex[terminal];
mFlux[edge_index] += estimated_terminal_flux_change;
}
SolveDirectFromFlux();
/* Look at the magnitude of the response */
for (unsigned terminal=0; terminal<num_terminals; terminal++)
{
unsigned node_index = mTerminalToNodeIndex[terminal];
// How far we are away from matching this boundary condition.
double delta_pressure = mPressure[node_index] - mPressureCondition[node_index];
// Offset the first iteration
if (iteration == 0)
{
delta_pressure += mPressureCondition[mOutletNodeIndex];
}
VecSetValue(mTerminalPressureChangeVector, terminal, delta_pressure, INSERT_VALUES);
}
VecNorm(mTerminalPressureChangeVector, NORM_2, &norm_pressure_change);
if (norm_pressure_change < pressure_tolerance)
{
converged = true;
break;
}
double pressure_change_dot_product;
VecDot(mTerminalPressureChangeVector, old_terminal_pressure_change, &pressure_change_dot_product);
if (pressure_change_dot_product < 0.0)
{
/* The pressure correction has changed sign
* * so we have overshot the root
* * back up by setting a correction factor
*/
double terminal_flux_correction = last_norm_pressure_change / (last_norm_pressure_change + norm_pressure_change) - 1.0;
///\todo some sqrt response?
// if (mDynamicResistance)
// {
// terminal_flux_correction *= 0.99;
// }
for (unsigned terminal=0; terminal<num_terminals; terminal++)
{
double estimated_terminal_flux_change=p_terminal_flux_change_vector[terminal];
unsigned edge_index = mTerminalToEdgeIndex[terminal];
mFlux[edge_index] += terminal_flux_correction*estimated_terminal_flux_change;
//.........这里部分代码省略.........
示例13: PCBDDCNullSpaceAdaptGlobal
PetscErrorCode PCBDDCNullSpaceAdaptGlobal(PC pc)
{
PC_IS* pcis = (PC_IS*) (pc->data);
PC_BDDC* pcbddc = (PC_BDDC*)(pc->data);
KSP inv_change;
PC pc_change;
const Vec *nsp_vecs;
Vec *new_nsp_vecs;
PetscInt i,nsp_size,new_nsp_size,start_new;
PetscBool nsp_has_cnst;
MatNullSpace new_nsp;
PetscErrorCode ierr;
MPI_Comm comm;
PetscFunctionBegin;
/* create KSP for change of basis */
ierr = KSPCreate(PETSC_COMM_SELF,&inv_change);CHKERRQ(ierr);
ierr = KSPSetOperators(inv_change,pcbddc->ChangeOfBasisMatrix,pcbddc->ChangeOfBasisMatrix);CHKERRQ(ierr);
ierr = KSPSetType(inv_change,KSPPREONLY);CHKERRQ(ierr);
ierr = KSPGetPC(inv_change,&pc_change);CHKERRQ(ierr);
ierr = PCSetType(pc_change,PCLU);CHKERRQ(ierr);
ierr = KSPSetUp(inv_change);CHKERRQ(ierr);
/* get nullspace and transform it */
ierr = MatNullSpaceGetVecs(pcbddc->NullSpace,&nsp_has_cnst,&nsp_size,&nsp_vecs);CHKERRQ(ierr);
new_nsp_size = nsp_size;
if (nsp_has_cnst) {
new_nsp_size++;
}
ierr = PetscMalloc1(new_nsp_size,&new_nsp_vecs);CHKERRQ(ierr);
for (i=0;i<new_nsp_size;i++) {
ierr = VecDuplicate(pcis->vec1_global,&new_nsp_vecs[i]);CHKERRQ(ierr);
}
start_new = 0;
if (nsp_has_cnst) {
start_new = 1;
ierr = VecSet(new_nsp_vecs[0],1.0);CHKERRQ(ierr);
ierr = VecSet(pcis->vec1_B,1.0);CHKERRQ(ierr);
ierr = KSPSolve(inv_change,pcis->vec1_B,pcis->vec1_B);
ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,new_nsp_vecs[0],INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
ierr = VecScatterEnd (pcis->global_to_B,pcis->vec1_B,new_nsp_vecs[0],INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
}
for (i=0;i<nsp_size;i++) {
ierr = VecCopy(nsp_vecs[i],new_nsp_vecs[i+start_new]);CHKERRQ(ierr);
ierr = VecScatterBegin(pcis->global_to_B,nsp_vecs[i],pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
ierr = VecScatterEnd (pcis->global_to_B,nsp_vecs[i],pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
ierr = KSPSolve(inv_change,pcis->vec1_B,pcis->vec1_B);
ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,new_nsp_vecs[i+start_new],INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
ierr = VecScatterEnd (pcis->global_to_B,pcis->vec1_B,new_nsp_vecs[i+start_new],INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
}
ierr = PCBDDCOrthonormalizeVecs(new_nsp_size,new_nsp_vecs);CHKERRQ(ierr);
#if 0
PetscBool nsp_t=PETSC_FALSE;
ierr = MatNullSpaceTest(pcbddc->NullSpace,pc->pmat,&nsp_t);CHKERRQ(ierr);
printf("Original Null Space test: %d\n",nsp_t);
Mat temp_mat;
Mat_IS* matis = (Mat_IS*)pc->pmat->data;
temp_mat = matis->A;
matis->A = pcbddc->local_mat;
pcbddc->local_mat = temp_mat;
ierr = MatNullSpaceTest(pcbddc->NullSpace,pc->pmat,&nsp_t);CHKERRQ(ierr);
printf("Original Null Space, mat changed test: %d\n",nsp_t);
{
PetscReal test_norm;
for (i=0;i<new_nsp_size;i++) {
ierr = MatMult(pc->pmat,new_nsp_vecs[i],pcis->vec1_global);CHKERRQ(ierr);
ierr = VecNorm(pcis->vec1_global,NORM_2,&test_norm);CHKERRQ(ierr);
if (test_norm > 1.e-12) {
printf("------------ERROR VEC %d------------------\n",i);
ierr = VecView(pcis->vec1_global,PETSC_VIEWER_STDOUT_WORLD);
printf("------------------------------------------\n");
}
}
}
#endif
ierr = KSPDestroy(&inv_change);CHKERRQ(ierr);
ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr);
ierr = MatNullSpaceCreate(comm,PETSC_FALSE,new_nsp_size,new_nsp_vecs,&new_nsp);CHKERRQ(ierr);
ierr = PCBDDCSetNullSpace(pc,new_nsp);CHKERRQ(ierr);
ierr = MatNullSpaceDestroy(&new_nsp);CHKERRQ(ierr);
#if 0
ierr = MatNullSpaceTest(pcbddc->NullSpace,pc->pmat,&nsp_t);CHKERRQ(ierr);
printf("New Null Space, mat changed: %d\n",nsp_t);
temp_mat = matis->A;
matis->A = pcbddc->local_mat;
pcbddc->local_mat = temp_mat;
ierr = MatNullSpaceTest(pcbddc->NullSpace,pc->pmat,&nsp_t);CHKERRQ(ierr);
printf("New Null Space, mat original: %d\n",nsp_t);
#endif
for (i=0;i<new_nsp_size;i++) {
ierr = VecDestroy(&new_nsp_vecs[i]);CHKERRQ(ierr);
}
ierr = PetscFree(new_nsp_vecs);CHKERRQ(ierr);
PetscFunctionReturn(0);
}示例14: MatFDColoringApply_AIJ
PetscErrorCode MatFDColoringApply_AIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx)
{
PetscErrorCode (*f)(void*,Vec,Vec,void*) = (PetscErrorCode (*)(void*,Vec,Vec,void *))coloring->f;
PetscErrorCode ierr;
PetscInt k,start,end,l,row,col,srow,**vscaleforrow;
PetscScalar dx,*y,*xx,*w3_array;
PetscScalar *vscale_array;
PetscReal epsilon = coloring->error_rel,umin = coloring->umin,unorm;
Vec w1=coloring->w1,w2=coloring->w2,w3;
void *fctx = coloring->fctx;
PetscBool flg = PETSC_FALSE;
PetscInt ctype=coloring->ctype,N,col_start=0,col_end=0;
Vec x1_tmp;
PetscFunctionBegin;
PetscValidHeaderSpecific(J,MAT_CLASSID,1);
PetscValidHeaderSpecific(coloring,MAT_FDCOLORING_CLASSID,2);
PetscValidHeaderSpecific(x1,VEC_CLASSID,3);
if (!f) SETERRQ(((PetscObject)J)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must call MatFDColoringSetFunction()");
ierr = PetscLogEventBegin(MAT_FDColoringApply,coloring,J,x1,0);CHKERRQ(ierr);
ierr = MatSetUnfactored(J);CHKERRQ(ierr);
ierr = PetscOptionsGetBool(PETSC_NULL,"-mat_fd_coloring_dont_rezero",&flg,PETSC_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);
}
}
x1_tmp = x1;
if (!coloring->vscale){
ierr = VecDuplicate(x1_tmp,&coloring->vscale);CHKERRQ(ierr);
}
if (coloring->htype[0] == 'w') { /* tacky test; need to make systematic if we add other approaches to computing h*/
ierr = VecNorm(x1_tmp,NORM_2,&unorm);CHKERRQ(ierr);
}
ierr = VecGetOwnershipRange(w1,&start,&end);CHKERRQ(ierr); /* OwnershipRange is used by ghosted x! */
/* Set w1 = F(x1) */
if (!coloring->fset) {
ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr);
ierr = (*f)(sctx,x1_tmp,w1,fctx);CHKERRQ(ierr);
ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr);
} else {
coloring->fset = PETSC_FALSE;
}
if (!coloring->w3) {
ierr = VecDuplicate(x1_tmp,&coloring->w3);CHKERRQ(ierr);
ierr = PetscLogObjectParent(coloring,coloring->w3);CHKERRQ(ierr);
}
w3 = coloring->w3;
/* Compute all the local scale factors, including ghost points */
ierr = VecGetLocalSize(x1_tmp,&N);CHKERRQ(ierr);
ierr = VecGetArray(x1_tmp,&xx);CHKERRQ(ierr);
ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr);
if (ctype == IS_COLORING_GHOSTED){
col_start = 0; col_end = N;
} else if (ctype == IS_COLORING_GLOBAL){
xx = xx - start;
vscale_array = vscale_array - start;
col_start = start; col_end = N + start;
}
for (col=col_start; col<col_end; col++){
/* Loop over each local column, vscale[col] = 1./(epsilon*dx[col]) */
if (coloring->htype[0] == 'w') {
dx = 1.0 + unorm;
} else {
dx = xx[col];
}
if (dx == (PetscScalar)0.0) dx = 1.0;
if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin;
else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin;
dx *= epsilon;
vscale_array[col] = (PetscScalar)1.0/dx;
}
if (ctype == IS_COLORING_GLOBAL) vscale_array = vscale_array + start;
ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr);
if (ctype == IS_COLORING_GLOBAL){
ierr = VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
ierr = VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
}
if (coloring->vscaleforrow) {
vscaleforrow = coloring->vscaleforrow;
} else SETERRQ(((PetscObject)J)->comm,PETSC_ERR_ARG_NULL,"Null Object: coloring->vscaleforrow");
/*
Loop over each color
*/
ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr);
for (k=0; k<coloring->ncolors; k++) {
coloring->currentcolor = k;
ierr = VecCopy(x1_tmp,w3);CHKERRQ(ierr);
//.........这里部分代码省略.........
示例15: PCBDDCNullSpaceAssembleCoarse
PetscErrorCode PCBDDCNullSpaceAssembleCoarse(PC pc, Mat coarse_mat, MatNullSpace* CoarseNullSpace)
{
PC_BDDC *pcbddc = (PC_BDDC*)pc->data;
Mat_IS *matis = (Mat_IS*)pc->pmat->data;
MatNullSpace tempCoarseNullSpace;
const Vec *nsp_vecs;
Vec *coarse_nsp_vecs,local_vec,local_primal_vec;
PetscInt nsp_size,coarse_nsp_size,i;
PetscBool nsp_has_cnst;
PetscReal test_null;
PetscErrorCode ierr;
PetscFunctionBegin;
tempCoarseNullSpace = 0;
coarse_nsp_size = 0;
coarse_nsp_vecs = 0;
ierr = MatNullSpaceGetVecs(pcbddc->NullSpace,&nsp_has_cnst,&nsp_size,&nsp_vecs);CHKERRQ(ierr);
if (coarse_mat) {
ierr = PetscMalloc1((nsp_size+1),&coarse_nsp_vecs);CHKERRQ(ierr);
for (i=0;i<nsp_size+1;i++) {
ierr = VecDuplicate(pcbddc->coarse_vec,&coarse_nsp_vecs[i]);CHKERRQ(ierr);
}
}
ierr = MatGetVecs(pcbddc->ConstraintMatrix,&local_vec,&local_primal_vec);CHKERRQ(ierr);
if (nsp_has_cnst) {
ierr = VecSet(local_vec,1.0);CHKERRQ(ierr);
ierr = MatMult(pcbddc->ConstraintMatrix,local_vec,local_primal_vec);CHKERRQ(ierr);
ierr = PCBDDCScatterCoarseDataBegin(pc,local_primal_vec,pcbddc->coarse_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
ierr = PCBDDCScatterCoarseDataEnd(pc,local_primal_vec,pcbddc->coarse_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
if (coarse_mat) {
if (pcbddc->dbg_flag) {
ierr = MatMult(coarse_mat,pcbddc->coarse_vec,pcbddc->coarse_rhs);CHKERRQ(ierr);
ierr = VecNorm(pcbddc->coarse_rhs,NORM_INFINITY,&test_null);CHKERRQ(ierr);
ierr = PetscViewerASCIIPrintf(pcbddc->dbg_viewer,"Constant coarse null space error % 1.14e\n",test_null);CHKERRQ(ierr);
}
ierr = VecCopy(pcbddc->coarse_vec,coarse_nsp_vecs[coarse_nsp_size]);CHKERRQ(ierr);
coarse_nsp_size++;
}
}
for (i=0;i<nsp_size;i++) {
ierr = VecScatterBegin(matis->ctx,nsp_vecs[i],local_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
ierr = VecScatterEnd(matis->ctx,nsp_vecs[i],local_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
ierr = MatMult(pcbddc->ConstraintMatrix,local_vec,local_primal_vec);CHKERRQ(ierr);
ierr = PCBDDCScatterCoarseDataBegin(pc,local_primal_vec,pcbddc->coarse_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
ierr = PCBDDCScatterCoarseDataEnd(pc,local_primal_vec,pcbddc->coarse_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
if (coarse_mat) {
if (pcbddc->dbg_flag) {
ierr = MatMult(coarse_mat,pcbddc->coarse_vec,pcbddc->coarse_rhs);CHKERRQ(ierr);
ierr = VecNorm(pcbddc->coarse_rhs,NORM_2,&test_null);CHKERRQ(ierr);
ierr = PetscViewerASCIIPrintf(pcbddc->dbg_viewer,"Vec %d coarse null space error % 1.14e\n",i,test_null);CHKERRQ(ierr);
}
ierr = VecCopy(pcbddc->coarse_vec,coarse_nsp_vecs[coarse_nsp_size]);CHKERRQ(ierr);
coarse_nsp_size++;
}
}
if (coarse_nsp_size > 0) {
ierr = PCBDDCOrthonormalizeVecs(coarse_nsp_size,coarse_nsp_vecs);CHKERRQ(ierr);
ierr = MatNullSpaceCreate(PetscObjectComm((PetscObject)coarse_mat),PETSC_FALSE,coarse_nsp_size,coarse_nsp_vecs,&tempCoarseNullSpace);CHKERRQ(ierr);
for (i=0;i<nsp_size+1;i++) {
ierr = VecDestroy(&coarse_nsp_vecs[i]);CHKERRQ(ierr);
}
}
ierr = PetscFree(coarse_nsp_vecs);CHKERRQ(ierr);
ierr = VecDestroy(&local_vec);CHKERRQ(ierr);
ierr = VecDestroy(&local_primal_vec);CHKERRQ(ierr);
*CoarseNullSpace = tempCoarseNullSpace;
PetscFunctionReturn(0);
}