本文整理汇总了C++中VecDestroy函数的典型用法代码示例。如果您正苦于以下问题：C++ VecDestroy函数的具体用法？C++ VecDestroy怎么用？C++ VecDestroy使用的例子？那么, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了VecDestroy函数的15个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: 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);
}示例2: test1_DAInjection3d
PetscErrorCode test1_DAInjection3d(PetscInt mx, PetscInt my, PetscInt mz)
{
PetscErrorCode ierr;
DM dac,daf;
PetscViewer vv;
Vec ac,af;
PetscInt periodicity;
DMBoundaryType bx,by,bz;
PetscFunctionBeginUser;
bx = DM_BOUNDARY_NONE;
by = DM_BOUNDARY_NONE;
bz = DM_BOUNDARY_NONE;
periodicity = 0;
ierr = PetscOptionsGetInt(NULL,"-periodic", &periodicity, NULL);CHKERRQ(ierr);
if (periodicity==1) {
bx = DM_BOUNDARY_PERIODIC;
} else if (periodicity==2) {
by = DM_BOUNDARY_PERIODIC;
} else if (periodicity==3) {
bz = DM_BOUNDARY_PERIODIC;
}
ierr = DMDACreate3d(PETSC_COMM_WORLD, bx,by,bz, DMDA_STENCIL_BOX,
mx+1, my+1,mz+1,
PETSC_DECIDE, PETSC_DECIDE,PETSC_DECIDE,
1, /* 1 dof */
1, /* stencil = 1 */
NULL,NULL,NULL,
&daf);CHKERRQ(ierr);
ierr = DMSetFromOptions(daf);CHKERRQ(ierr);
ierr = DMCoarsen(daf,MPI_COMM_NULL,&dac);CHKERRQ(ierr);
ierr = DMDASetUniformCoordinates(dac, -1.0,1.0, -1.0,1.0, -1.0,1.0);CHKERRQ(ierr);
ierr = DMDASetUniformCoordinates(daf, -1.0,1.0, -1.0,1.0, -1.0,1.0);CHKERRQ(ierr);
{
DM cdaf,cdac;
Vec coordsc,coordsf,coordsf2;
VecScatter inject;
Mat interp;
PetscReal norm;
ierr = DMGetCoordinateDM(dac,&cdac);CHKERRQ(ierr);
ierr = DMGetCoordinateDM(daf,&cdaf);CHKERRQ(ierr);
ierr = DMGetCoordinates(dac,&coordsc);CHKERRQ(ierr);
ierr = DMGetCoordinates(daf,&coordsf);CHKERRQ(ierr);
ierr = DMCreateInjection(cdac,cdaf,&inject);CHKERRQ(ierr);
ierr = VecScatterBegin(inject,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
ierr = VecScatterEnd(inject ,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
ierr = VecScatterDestroy(&inject);CHKERRQ(ierr);
ierr = DMCreateInterpolation(cdac,cdaf,&interp,NULL);CHKERRQ(ierr);
ierr = VecDuplicate(coordsf,&coordsf2);CHKERRQ(ierr);
ierr = MatInterpolate(interp,coordsc,coordsf2);CHKERRQ(ierr);
ierr = VecAXPY(coordsf2,-1.0,coordsf);CHKERRQ(ierr);
ierr = VecNorm(coordsf2,NORM_MAX,&norm);CHKERRQ(ierr);
/* The fine coordinates are only reproduced in certain cases */
if (!bx && !by && !bz && norm > 1.e-10) {ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm %g\n",(double)norm);CHKERRQ(ierr);}
ierr = VecDestroy(&coordsf2);CHKERRQ(ierr);
ierr = MatDestroy(&interp);CHKERRQ(ierr);
}
if (0) {
ierr = DMCreateGlobalVector(dac,&ac);CHKERRQ(ierr);
ierr = VecZeroEntries(ac);CHKERRQ(ierr);
ierr = DMCreateGlobalVector(daf,&af);CHKERRQ(ierr);
ierr = VecZeroEntries(af);CHKERRQ(ierr);
ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_7.vtk", &vv);CHKERRQ(ierr);
ierr = PetscViewerSetFormat(vv, PETSC_VIEWER_ASCII_VTK);CHKERRQ(ierr);
ierr = DMView(dac, vv);CHKERRQ(ierr);
ierr = VecView(ac, vv);CHKERRQ(ierr);
ierr = PetscViewerDestroy(&vv);CHKERRQ(ierr);
ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_7.vtk", &vv);CHKERRQ(ierr);
ierr = PetscViewerSetFormat(vv, PETSC_VIEWER_ASCII_VTK);CHKERRQ(ierr);
ierr = DMView(daf, vv);CHKERRQ(ierr);
ierr = VecView(af, vv);CHKERRQ(ierr);
ierr = PetscViewerDestroy(&vv);CHKERRQ(ierr);
ierr = VecDestroy(&ac);CHKERRQ(ierr);
ierr = VecDestroy(&af);CHKERRQ(ierr);
}
ierr = DMDestroy(&dac);CHKERRQ(ierr);
ierr = DMDestroy(&daf);CHKERRQ(ierr);
PetscFunctionReturn(0);
}示例3: MatGetOrdering_AWBM
//.........这里部分代码省略.........
if (debug) {
ierr = PetscPrintf(PETSC_COMM_SELF, "Replaced match %d -- %d\n", c1, ja[r]);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_SELF, " Added match %d -- %d\n", c, ja[r]);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_SELF, " Added match %d -- %d\n", c1, ja[r1]);CHKERRQ(ierr);
}
match[c] = ja[r];
matchR[ja[r]] = c;
match[c1] = ja[r1];
matchR[ja[r1]] = c1;
break;
}
}
if (match[c] >= 0) break;
}
}
/* Allow matching with non-optimal rows */
ierr = CheckUnmatched(n, match, matchR);CHKERRQ(ierr);
for (c = 0; c < n; ++c) {
if (match[c] >= 0) continue;
for (r = ia[c]; r < ia[c+1]; ++r) {
if (matchR[ja[r]] < 0) {
if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "Matched non-opt %d -- %d\n", c, ja[r]);CHKERRQ(ierr);}
match[c] = ja[r];
matchR[ja[r]] = c;
break;
}
}
}
/* Deal with non-optimal unmatched columns */
ierr = CheckUnmatched(n, match, matchR);CHKERRQ(ierr);
for (c = 0; c < n; ++c) {
if (match[c] >= 0) continue;
for (r = ia[c]; r < ia[c+1]; ++r) {
/* \bar c_ij = 0 and (r, j1) \in M */
c1 = matchR[ja[r]];
for (r1 = ia[c1]; r1 < ia[c1+1]; ++r1) {
if (matchR[ja[r1]] < 0) {
/* (r, c1) in M is replaced by (r, c) and (r1, c1) */
if (debug) {
ierr = PetscPrintf(PETSC_COMM_SELF, "Replaced match %d -- %d\n", c1, ja[r]);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_SELF, " Added match %d -- %d\n", c, ja[r]);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_SELF, " Added match %d -- %d\n", c1, ja[r1]);CHKERRQ(ierr);
}
match[c] = ja[r];
matchR[ja[r]] = c;
match[c1] = ja[r1];
matchR[ja[r1]] = c1;
break;
}
}
if (match[c] >= 0) break;
}
}
/* Complete matching */
ierr = CheckUnmatched(n, match, matchR);CHKERRQ(ierr);
for (c = 0, r = 0; c < n; ++c) {
if (match[c] >= n) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Column %d matched to invalid row %d", c, match[c]);
if (match[c] < 0) {
for (; r < n; ++r) {
if (matchR[r] < 0) {
if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "Matched default %d -- %d\n", c, r);CHKERRQ(ierr);}
match[c] = r;
matchR[r] = c;
break;
}
}
}
}
/* Check matching */
ierr = CheckUnmatched(n, match, matchR);CHKERRQ(ierr);
for (c = 0; c < n; ++c) {
if (match[c] < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Column %d unmatched", c);
if (match[c] >= n) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Column %d matched to invalid row %d", c, match[c]);
}
/* Make permutation */
for (c = 0; c < n; ++c) {p[match[c]] = c;}
ierr = ISCreateGeneral(PETSC_COMM_SELF, n, p, PETSC_OWN_POINTER, permR);CHKERRQ(ierr);
ierr = ISSetPermutation(*permR);CHKERRQ(ierr);
ierr = ISCreateStride(PETSC_COMM_SELF, n, 0, 1, permC);CHKERRQ(ierr);
ierr = ISSetPermutation(*permC);CHKERRQ(ierr);
ierr = PetscFree2(match, matchR);CHKERRQ(ierr);
/* Make scaling */
ierr = VecCreateSeq(PETSC_COMM_SELF, n, scalR);CHKERRQ(ierr);
ierr = VecCreateSeq(PETSC_COMM_SELF, n, scalC);CHKERRQ(ierr);
ierr = VecGetArray(*scalR, &sr);CHKERRQ(ierr);
ierr = VecGetArray(*scalC, &sc);CHKERRQ(ierr);
for (c = 0; c < n; ++c) {
sr[c] = PetscExpReal(v[c])/a_j[c];
sc[c] = PetscExpReal(u[c]);
}
ierr = VecRestoreArray(*scalR, &sr);CHKERRQ(ierr);
ierr = VecRestoreArray(*scalC, &sc);CHKERRQ(ierr);
ierr = VecRestoreArray(colMax, &a_j);CHKERRQ(ierr);
ierr = VecDestroy(&colMax);CHKERRQ(ierr);
ierr = PetscFree3(u,v,weights);CHKERRQ(ierr);
ierr = VecDestroy(scalR);CHKERRQ(ierr);
ierr = VecDestroy(scalC);CHKERRQ(ierr);
PetscFunctionReturn(0);
}示例4: main
//.........这里部分代码省略.........
N1, N2,
PETSC_DECIDE, PETSC_DECIDE,
DOF, NG, PETSC_NULL, PETSC_NULL, &dmda);
// When running in parallel, each process computes from
// [X1Start, X1Start+X1Size] x [X2Start, X2Start+X2Size]
DMDAGetCorners(dmda,
&X1Start, &X2Start, NULL,
&X1Size, &X2Size, NULL);
// Create the solution vector.
DMCreateGlobalVector(dmda, &soln);
// Create the time stepper and link it to the computational grid and the
// residual evaluation function.
TSCreate(PETSC_COMM_WORLD, &ts);
TSSetDM(ts, dmda);
TSSetIFunction(ts, PETSC_NULL, ComputeResidual, NULL);
// OpenCL boilerplate code.
clErr = cl::Platform::get(&platforms);
CheckCLErrors(clErr, "cl::Platform::get");
// Select computation device here.
clErr = platforms.at(1).getDevices(CL_DEVICE_TYPE_CPU, &devices);
CheckCLErrors(clErr, "cl::Platform::getDevices");
context = cl::Context(devices, NULL, NULL, NULL, &clErr);
CheckCLErrors(clErr, "cl::Context::Context");
queue = cl::CommandQueue(context, devices.at(0), 0, &clErr);
CheckCLErrors(clErr, "cl::CommandQueue::CommandQueue");
std::ifstream sourceFile("computeresidual.cl");
std::string sourceCode((std::istreambuf_iterator<char>(sourceFile)),
std::istreambuf_iterator<char>());
cl::Program::Sources source(1, std::make_pair(sourceCode.c_str(),
sourceCode.length()+1));
program = cl::Program(context, source, &clErr);
CheckCLErrors(clErr, "cl::Program::Program");
// Pass in constants to the OpenCL kernel as compiler switches. This is an
// efficient way to handle constants such as domain sizes in OpenCL.
std::string BuildOptions("\
-D X1_SIZE=" +
std::to_string(X1Size) +
" -D X2_SIZE=" +
std::to_string(X2Size) +
" -D TOTAL_X1_SIZE=" +
std::to_string(X1Size+2*NG) +
" -D TOTAL_X2_SIZE=" +
std::to_string(X2Size+2*NG));
// Compile the OpenCL program and extract the kernel.
PetscScalar start = std::clock();
clErr = program.build(devices, BuildOptions.c_str(), NULL, NULL);
const char *buildlog = program.getBuildInfo<CL_PROGRAM_BUILD_LOG>(
devices.at(0),
&clErr).c_str();
PetscPrintf(PETSC_COMM_WORLD, "%s\n", buildlog);
CheckCLErrors(clErr, "cl::Program::build");
PetscScalar end = std::clock();
PetscScalar time = (end - start)/(PetscScalar)CLOCKS_PER_SEC;
PetscPrintf(PETSC_COMM_WORLD,
"Time taken for kernel compilation = %f\n", time);
kernel = cl::Kernel(program, "ComputeResidual", &clErr);
CheckCLErrors(clErr, "cl::Kernel::Kernel");
// How much memory is the kernel using?
cl_ulong localMemSize = kernel.getWorkGroupInfo<CL_KERNEL_LOCAL_MEM_SIZE>(
devices.at(0), &clErr);
cl_ulong privateMemSize = kernel.getWorkGroupInfo<CL_KERNEL_PRIVATE_MEM_SIZE>(
devices.at(0), &clErr);
printf("Local memory used = %llu\n", (unsigned long long)localMemSize);
printf("Private memory used = %llu\n", (unsigned long long)privateMemSize);
// Set initial conditions.
InitialCondition(ts, soln);
TSSetSolution(ts, soln);
TSSetType(ts, TSTHETA);
TSSetFromOptions(ts);
// Finally solve! All time stepping options can be controlled from the
// command line.
TSSolve(ts, soln);
// Delete the data structures in the following order.
DMDestroy(&dmda);
VecDestroy(&soln);
TSDestroy(&ts);
PetscFinalize();
return(0);
}示例5: MatGetSubMatrix
/*@
MatCreateSubMatrix - Creates a composite matrix that acts as a submatrix
Collective on Mat
Input Parameters:
+ A - matrix that we will extract a submatrix of
. isrow - rows to be present in the submatrix
- iscol - columns to be present in the submatrix
Output Parameters:
. newmat - new matrix
Level: developer
Notes:
Most will use MatGetSubMatrix which provides a more efficient representation if it is available.
.seealso: MatGetSubMatrix(), MatSubMatrixUpdate()
@*/
PetscErrorCode MatCreateSubMatrix(Mat A,IS isrow,IS iscol,Mat *newmat)
{
Vec left,right;
PetscInt m,n;
Mat N;
Mat_SubMatrix *Na;
PetscErrorCode ierr;
PetscFunctionBegin;
PetscValidHeaderSpecific(A,MAT_CLASSID,1);
PetscValidHeaderSpecific(isrow,IS_CLASSID,2);
PetscValidHeaderSpecific(iscol,IS_CLASSID,3);
PetscValidPointer(newmat,4);
*newmat = 0;
ierr = MatCreate(PetscObjectComm((PetscObject)A),&N);CHKERRQ(ierr);
ierr = ISGetLocalSize(isrow,&m);CHKERRQ(ierr);
ierr = ISGetLocalSize(iscol,&n);CHKERRQ(ierr);
ierr = MatSetSizes(N,m,n,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr);
ierr = PetscObjectChangeTypeName((PetscObject)N,MATSUBMATRIX);CHKERRQ(ierr);
ierr = PetscNewLog(N,&Na);CHKERRQ(ierr);
N->data = (void*)Na;
ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr);
ierr = PetscObjectReference((PetscObject)isrow);CHKERRQ(ierr);
ierr = PetscObjectReference((PetscObject)iscol);CHKERRQ(ierr);
Na->A = A;
Na->isrow = isrow;
Na->iscol = iscol;
Na->scale = 1.0;
N->ops->destroy = MatDestroy_SubMatrix;
N->ops->mult = MatMult_SubMatrix;
N->ops->multadd = MatMultAdd_SubMatrix;
N->ops->multtranspose = MatMultTranspose_SubMatrix;
N->ops->multtransposeadd = MatMultTransposeAdd_SubMatrix;
N->ops->scale = MatScale_SubMatrix;
N->ops->diagonalscale = MatDiagonalScale_SubMatrix;
ierr = PetscLayoutSetBlockSize(N->rmap,A->rmap->bs);CHKERRQ(ierr);
ierr = PetscLayoutSetBlockSize(N->cmap,A->cmap->bs);CHKERRQ(ierr);
ierr = PetscLayoutSetUp(N->rmap);CHKERRQ(ierr);
ierr = PetscLayoutSetUp(N->cmap);CHKERRQ(ierr);
ierr = MatGetVecs(A,&Na->rwork,&Na->lwork);CHKERRQ(ierr);
ierr = VecCreate(PetscObjectComm((PetscObject)isrow),&left);CHKERRQ(ierr);
ierr = VecCreate(PetscObjectComm((PetscObject)iscol),&right);CHKERRQ(ierr);
ierr = VecSetSizes(left,m,PETSC_DETERMINE);CHKERRQ(ierr);
ierr = VecSetSizes(right,n,PETSC_DETERMINE);CHKERRQ(ierr);
ierr = VecSetUp(left);CHKERRQ(ierr);
ierr = VecSetUp(right);CHKERRQ(ierr);
ierr = VecScatterCreate(Na->lwork,isrow,left,NULL,&Na->lrestrict);CHKERRQ(ierr);
ierr = VecScatterCreate(right,NULL,Na->rwork,iscol,&Na->rprolong);CHKERRQ(ierr);
ierr = VecDestroy(&left);CHKERRQ(ierr);
ierr = VecDestroy(&right);CHKERRQ(ierr);
N->assembled = PETSC_TRUE;
ierr = MatSetUp(N);CHKERRQ(ierr);
*newmat = N;
PetscFunctionReturn(0);
}示例6: LoadTestMatrices
PetscErrorCode LoadTestMatrices(Mat *_A,Vec *_x,Vec *_b,IS *_isu,IS *_isp)
{
Vec f,h,x,b,bX[2];
Mat A,Auu,Aup,Apu,App,bA[2][2];
IS is_u,is_p,bis[2];
PetscInt lnu,lnp,nu,np,i,start_u,end_u,start_p,end_p;
VecScatter *vscat;
PetscMPIInt rank;
PetscErrorCode ierr;
PetscFunctionBeginUser;
/* fetch test matrices and vectors */
ierr = LSCLoadTestOperators(&Auu,&Aup,&Apu,&App,&f,&h);CHKERRQ(ierr);
/* build the mat-nest */
ierr = VecGetSize(f,&nu);CHKERRQ(ierr);
ierr = VecGetSize(h,&np);CHKERRQ(ierr);
ierr = VecGetLocalSize(f,&lnu);CHKERRQ(ierr);
ierr = VecGetLocalSize(h,&lnp);CHKERRQ(ierr);
ierr = VecGetOwnershipRange(f,&start_u,&end_u);CHKERRQ(ierr);
ierr = VecGetOwnershipRange(h,&start_p,&end_p);CHKERRQ(ierr);
ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d] lnu = %D | lnp = %D \n", rank, lnu, lnp);CHKERRQ(ierr);
ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d] s_u = %D | e_u = %D \n", rank, start_u, end_u);CHKERRQ(ierr);
ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d] s_p = %D | e_p = %D \n", rank, start_p, end_p);CHKERRQ(ierr);
ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d] is_u (offset) = %D \n", rank, start_u+start_p);CHKERRQ(ierr);
ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d] is_p (offset) = %D \n", rank, start_u+start_p+lnu);CHKERRQ(ierr);
ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr);
ierr = ISCreateStride(PETSC_COMM_WORLD,lnu,start_u+start_p,1,&is_u);CHKERRQ(ierr);
ierr = ISCreateStride(PETSC_COMM_WORLD,lnp,start_u+start_p+lnu,1,&is_p);CHKERRQ(ierr);
bis[0] = is_u; bis[1] = is_p;
bA[0][0] = Auu; bA[0][1] = Aup;
bA[1][0] = Apu; bA[1][1] = App;
ierr = MatCreateNest(PETSC_COMM_WORLD,2,bis,2,bis,&bA[0][0],&A);CHKERRQ(ierr);
ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
/* Pull f,h into b */
ierr = MatCreateVecs(A,&b,&x);CHKERRQ(ierr);
bX[0] = f; bX[1] = h;
ierr = PetscMalloc1(2,&vscat);CHKERRQ(ierr);
for (i=0; i<2; i++) {
ierr = VecScatterCreateWithData(b,bis[i],bX[i],NULL,&vscat[i]);CHKERRQ(ierr);
ierr = VecScatterBegin(vscat[i],bX[i],b,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
}
for (i=0; i<2; i++) {
ierr = VecScatterEnd(vscat[i],bX[i],b,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
}
/* tidy up */
for (i=0; i<2; i++) {
ierr = VecScatterDestroy(&vscat[i]);CHKERRQ(ierr);
}
ierr = PetscFree(vscat);CHKERRQ(ierr);
ierr = MatDestroy(&Auu);CHKERRQ(ierr);
ierr = MatDestroy(&Aup);CHKERRQ(ierr);
ierr = MatDestroy(&Apu);CHKERRQ(ierr);
ierr = MatDestroy(&App);CHKERRQ(ierr);
ierr = VecDestroy(&f);CHKERRQ(ierr);
ierr = VecDestroy(&h);CHKERRQ(ierr);
*_isu = is_u;
*_isp = is_p;
*_A = A;
*_x = x;
*_b = b;
PetscFunctionReturn(0);
}示例7: main
int main(int argc,char **args)
{
Mat C;
int i,m = 5,rank,size,N,start,end,M;
int ierr,idx[4];
PetscBool flg;
PetscScalar Ke[16];
PetscReal h;
Vec u,b;
KSP ksp;
MatNullSpace nullsp;
PetscInitialize(&argc,&args,(char *)0,help);
ierr = PetscOptionsGetInt(PETSC_NULL,"-m",&m,PETSC_NULL);CHKERRQ(ierr);
N = (m+1)*(m+1); /* dimension of matrix */
M = m*m; /* number of elements */
h = 1.0/m; /* mesh width */
ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
/* Create stiffness matrix */
ierr = MatCreate(PETSC_COMM_WORLD,&C);CHKERRQ(ierr);
ierr = MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,N,N);CHKERRQ(ierr);
ierr = MatSetFromOptions(C);CHKERRQ(ierr);
start = rank*(M/size) + ((M%size) < rank ? (M%size) : rank);
end = start + M/size + ((M%size) > rank);
/* Assemble matrix */
ierr = FormElementStiffness(h*h,Ke); /* element stiffness for Laplacian */
for (i=start; i<end; i++) {
/* location of lower left corner of element */
/* node numbers for the four corners of element */
idx[0] = (m+1)*(i/m) + (i % m);
idx[1] = idx[0]+1; idx[2] = idx[1] + m + 1; idx[3] = idx[2] - 1;
ierr = MatSetValues(C,4,idx,4,idx,Ke,ADD_VALUES);CHKERRQ(ierr);
}
ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
/* Create right-hand-side and solution vectors */
ierr = VecCreate(PETSC_COMM_WORLD,&u);CHKERRQ(ierr);
ierr = VecSetSizes(u,PETSC_DECIDE,N);CHKERRQ(ierr);
ierr = VecSetFromOptions(u);CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject)u,"Approx. Solution");CHKERRQ(ierr);
ierr = VecDuplicate(u,&b);CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject)b,"Right hand side");CHKERRQ(ierr);
ierr = VecSet(u,1.0);CHKERRQ(ierr);
ierr = MatMult(C,u,b);CHKERRQ(ierr);
ierr = VecSet(u,0.0);CHKERRQ(ierr);
/* Solve linear system */
ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
ierr = KSPSetOperators(ksp,C,C,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
ierr = KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);CHKERRQ(ierr);
flg = PETSC_FALSE;
ierr = PetscOptionsGetBool(PETSC_NULL,"-fixnullspace",&flg,PETSC_NULL);CHKERRQ(ierr);
if (flg) {
ierr = MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,PETSC_NULL,&nullsp);CHKERRQ(ierr);
ierr = KSPSetNullSpace(ksp,nullsp);CHKERRQ(ierr);
ierr = MatNullSpaceDestroy(&&nullsp);CHKERRQ(ierr);
}
ierr = KSPSolve(ksp,b,u);CHKERRQ(ierr);
/* Free work space */
ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
ierr = VecDestroy(&u);CHKERRQ(ierr);
ierr = VecDestroy(&b);CHKERRQ(ierr);
ierr = MatDestroy(&C);CHKERRQ(ierr);
ierr = PetscFinalize();
return 0;
}示例8: 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;
}示例9: PetscPrintf
~thing(){
PetscPrintf(PETSC_COMM_WORLD, "x before destroy = %i\n", x);
VecDestroy(&x);
PetscPrintf(PETSC_COMM_WORLD, "x after destroy = %i\n", x);
}示例10: testBSplineHAtom
int testBSplineHAtom() {
PrintTimeStamp(PETSC_COMM_SELF, "H atom", NULL);
MPI_Comm comm = PETSC_COMM_SELF;
BPS bps; BPSCreate(comm, &bps); BPSSetExp(bps, 30.0, 61, 3.0);
int order = 5;
BSS bss; BSSCreate(comm, &bss); BSSSetKnots(bss, order, bps); BSSSetUp(bss);
Mat H; BSSCreateR1Mat(bss, &H);
Mat S; BSSCreateR1Mat(bss, &S);
Mat V; BSSCreateR1Mat(bss, &V);
BSSD2R1Mat(bss, H);
MatScale(H, -0.5);
BSSENR1Mat(bss, 0, 0.0, V);
MatAXPY(H, -1.0, V, DIFFERENT_NONZERO_PATTERN);
BSSSR1Mat(bss, S);
// -- initial space --
Pot psi0; PotCreate(comm, &psi0); PotSetSlater(psi0, 2.0, 1, 1.1);
int n_init_space = 1;
Vec *xs; PetscMalloc1(n_init_space, &xs);
MatCreateVecs(H, &xs[0], NULL);
BSSPotR1Vec(bss, psi0, xs[0]);
EEPS eps; EEPSCreate(comm, &eps);
EEPSSetOperators(eps, H, S);
// EPSSetType(eps->eps, EPSJD);
EPSSetInitialSpace(eps->eps, 1, xs);
EEPSSetTarget(eps, -0.6);
// EPSSetInitialSpace(eps->eps, 1, xs);
EEPSSolve(eps);
int nconv;
PetscScalar kr;
EPSGetConverged(eps->eps, &nconv);
ASSERT_TRUE(nconv > 0);
EPSGetEigenpair(eps->eps, 0, &kr, NULL, NULL, NULL);
ASSERT_DOUBLE_NEAR(-0.5, kr, pow(10.0, -6.0));
Vec cs;
MatCreateVecs(H, &cs, NULL);
EEPSGetEigenvector(eps, 0, cs);
PetscReal x=1.1;
PetscScalar y=0.0;
PetscScalar dy=0.0;
BSSPsiOne(bss, cs, x, &y);
BSSDerivPsiOne(bss, cs, x, &dy);
ASSERT_DOUBLE_NEAR(creal(y), 2.0*x*exp(-x), pow(10.0, -6));
ASSERT_DOUBLE_NEAR(creal(dy), 2.0*exp(-x)-2.0*x*exp(-x), pow(10.0, -6));
VecDestroy(&xs[0]);
PetscFree(xs);
PFDestroy(&psi0);
BSSDestroy(&bss);
MatDestroy(&H);
MatDestroy(&V);
MatDestroy(&S);
EEPSDestroy(&eps);
VecDestroy(&cs);
return 0;
}示例11: testBSplinePsi
int testBSplinePsi() {
PetscErrorCode ierr;
PrintTimeStamp(PETSC_COMM_SELF, "Psi", NULL);
MPI_Comm comm = PETSC_COMM_SELF;
BPS bps; BPSCreate(comm, &bps); BPSSetLine(bps, 5.0, 6);
int order = 3;
BSS bss; BSSCreate(comm, &bss); BSSSetKnots(bss, order, bps);
BSSSetUp(bss);
int n; BSSGetSize(bss, &n);
if(n <= 0){
SETERRQ(comm, 1, "n is 0 or negative");
}
Vec c;
ierr = VecCreate(comm, &c); CHKERRQ(ierr);
ierr = VecSetSizes(c, PETSC_DECIDE, n);CHKERRQ(ierr);
ierr = VecSetUp(c); CHKERRQ(ierr);
PetscScalar *ptr_c;
ierr = VecGetArray(c, &ptr_c); CHKERRQ(ierr);
for(int i = 0; i < n; i++) {
ptr_c[i] = i + 0.1;
}
ierr = VecRestoreArray(c, &ptr_c); CHKERRQ(ierr);
PetscReal x = 1.1;
Vec xs; VecCreate(comm, &xs); VecSetSizes(xs, PETSC_DECIDE, 1);
ierr = VecSetUp(xs); CHKERRQ(ierr);
ierr = VecSetValue(xs, 0, x, INSERT_VALUES); CHKERRQ(ierr);
// Vec xs; VecCreate(comm, &xs); VecSetSizes(xs, PETSC_DEFAULT, 1);
// VecSetUp(xs);
// ierr = VecSetValue(xs, 0, x, INSERT_VALUES); CHKERRQ(ierr);
Vec ys;
ierr = VecDuplicate(xs, &ys); CHKERRQ(ierr);
ierr = BSSPsi(bss, c, xs, ys); CHKERRQ(ierr);
Vec dys;
ierr = VecDuplicate(xs, &dys); CHKERRQ(ierr);
ierr = BSSDerivPsi(bss, c, xs, dys); CHKERRQ(ierr);
PetscScalar y;
ierr = BSSPsiOne(bss, c, x, &y); CHKERRQ(ierr);
PetscScalar *ptr_y;
ierr = VecGetArray(ys, &ptr_y); CHKERRQ(ierr);
ASSERT_SCALAR_EQ(y, ptr_y[0]);
ierr = VecRestoreArray(ys, &ptr_y); CHKERRQ(ierr);
PetscScalar dy;
ierr = BSSDerivPsiOne(bss, c, x, &dy); CHKERRQ(ierr);
PetscScalar *ptr_dy;
ierr = VecGetArray(dys, &ptr_dy); CHKERRQ(ierr);
ASSERT_SCALAR_EQ(dy, ptr_dy[0]);
ierr = VecRestoreArray(dys, &ptr_dy); CHKERRQ(ierr);
BSSDestroy(&bss);
VecRestoreArray(c, &ptr_c);
VecDestroy(&c);
VecDestroy(&xs);
VecRestoreArray(ys, &ptr_y);
VecDestroy(&ys);
VecRestoreArray(dys, &ptr_dy);
VecDestroy(&dys);
return 0;
}示例12: output_mat_and_vec
/**
* output_mat_and_vec
* ------
* Output the matrices and vectors that can be used in MATLAB to solve various problems.
*/
PetscErrorCode output_mat_and_vec(const Mat A, const Vec b, const Vec right_precond, const Mat CF, const GridInfo gi)
{
PetscFunctionBegin;
PetscErrorCode ierr;
char output_name_prefixed[PETSC_MAX_PATH_LEN];
//const char *prefix = "/out/";
const char *ind_extension = "_ind";
const char *A_extension = "_A";
const char *b_extension = "_b";
const char *precond_extension = "_pR";
const char *CF_extension = "_CF";
//ierr = PetscStrcpy(output_name_prefixed, getenv("FD3D_ROOT")); CHKERRQ(ierr);
//ierr = PetscStrcat(output_name_prefixed, prefix); CHKERRQ(ierr);
//ierr = PetscStrcat(output_name_prefixed, gi.output_name); CHKERRQ(ierr);
ierr = PetscStrcpy(output_name_prefixed, gi.output_name); CHKERRQ(ierr);
char ind_file[PETSC_MAX_PATH_LEN];
char A_file[PETSC_MAX_PATH_LEN];
char b_file[PETSC_MAX_PATH_LEN];
char precond_file[PETSC_MAX_PATH_LEN];
char CF_file[PETSC_MAX_PATH_LEN];
ierr = PetscStrcpy(ind_file, output_name_prefixed); CHKERRQ(ierr);
ierr = PetscStrcat(ind_file, ind_extension); CHKERRQ(ierr);
ierr = PetscStrcpy(A_file, output_name_prefixed); CHKERRQ(ierr);
ierr = PetscStrcat(A_file, A_extension); CHKERRQ(ierr);
ierr = PetscStrcpy(b_file, output_name_prefixed); CHKERRQ(ierr);
ierr = PetscStrcat(b_file, b_extension); CHKERRQ(ierr);
ierr = PetscStrcpy(precond_file, output_name_prefixed); CHKERRQ(ierr);
ierr = PetscStrcat(precond_file, precond_extension); CHKERRQ(ierr);
ierr = PetscStrcpy(CF_file, output_name_prefixed); CHKERRQ(ierr);
ierr = PetscStrcat(CF_file, CF_extension); CHKERRQ(ierr);
PetscViewer viewer;
/** It turns out that VecView() shows the DA vector in natural order. Therefore, even though
indApp is constructed in application order, it is shown in natural order by VecView().
On the other hand, indNat reorder indApp in natural order and then distribute the vector to
processors. However, because VecView() reorder a vector before it prints out the content of
the vector, indNat is shown messy by VecView().
Inconsistently, MatView() does not reorder the matrix elements into natural order before it
shows the matrix. Therefore, when the binaries of matrices and vectors are imported in MATLAB,
I need to reorder the matrices but not the vectors. */
Vec indApp;
//ierr = create_index(&indApp, gi); CHKERRQ(ierr);
ierr = createFieldArray(&indApp, set_index_at, gi); CHKERRQ(ierr);
//ierr = VecView(indApp, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
Vec indNat;
ierr = DMDACreateNaturalVector(gi.da, &indNat); CHKERRQ(ierr);
ierr = VecCopy(indApp, indNat); CHKERRQ(ierr);
ierr = VecDestroy(&indApp); CHKERRQ(ierr);
//ierr = VecView(indNat, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
/** Write the index vector ind_app. */
ierr = PetscViewerCreate(PETSC_COMM_WORLD, &viewer); CHKERRQ(ierr);
ierr = PetscViewerSetType(viewer, PETSCVIEWERBINARY); CHKERRQ(ierr);
ierr = PetscViewerFileSetMode(viewer, FILE_MODE_WRITE); CHKERRQ(ierr);
ierr = PetscViewerFileSetName(viewer, ind_file); CHKERRQ(ierr);
ierr = VecView(indNat, viewer); CHKERRQ(ierr);
ierr = VecDestroy(&indNat); CHKERRQ(ierr);
ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
/** Write the coefficient matrix A. */
ierr = PetscViewerCreate(PETSC_COMM_WORLD, &viewer); CHKERRQ(ierr);
ierr = PetscViewerSetType(viewer, PETSCVIEWERBINARY); CHKERRQ(ierr);
ierr = PetscViewerFileSetMode(viewer, FILE_MODE_WRITE); CHKERRQ(ierr);
ierr = PetscViewerFileSetName(viewer, A_file); CHKERRQ(ierr);
ierr = MatView(A, viewer); CHKERRQ(ierr);
ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
/** Write the RHS vector b. */
ierr = PetscViewerCreate(PETSC_COMM_WORLD, &viewer); CHKERRQ(ierr);
ierr = PetscViewerSetType(viewer, PETSCVIEWERBINARY); CHKERRQ(ierr);
ierr = PetscViewerFileSetMode(viewer, FILE_MODE_WRITE); CHKERRQ(ierr);
ierr = PetscViewerFileSetName(viewer, b_file); CHKERRQ(ierr);
ierr = VecView(b, viewer); CHKERRQ(ierr);
ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
/** Write the right preconditioner vector pR. */
ierr = PetscViewerCreate(PETSC_COMM_WORLD, &viewer); CHKERRQ(ierr);
ierr = PetscViewerSetType(viewer, PETSCVIEWERBINARY); CHKERRQ(ierr);
ierr = PetscViewerFileSetMode(viewer, FILE_MODE_WRITE); CHKERRQ(ierr);
ierr = PetscViewerFileSetName(viewer, precond_file); CHKERRQ(ierr);
ierr = VecView(right_precond, viewer); CHKERRQ(ierr);
ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
/** Write the E-to-H converter matrix CF. */
ierr = PetscViewerCreate(PETSC_COMM_WORLD, &viewer); CHKERRQ(ierr);
ierr = PetscViewerSetType(viewer, PETSCVIEWERBINARY); CHKERRQ(ierr);
ierr = PetscViewerFileSetMode(viewer, FILE_MODE_WRITE); CHKERRQ(ierr);
ierr = PetscViewerFileSetName(viewer, CF_file); CHKERRQ(ierr);
ierr = MatView(CF, viewer); CHKERRQ(ierr);
//.........这里部分代码省略.........
示例13: main
int main(int argc,char **args)
{
Mat A[3],B; /* matrix */
PetscViewer fd; /* viewer */
char file[PETSC_MAX_PATH_LEN]; /* input file name */
PetscErrorCode ierr;
PetscBool flg;
Vec x,y,z,work;
PetscReal rnorm;
ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
/*
Determine files from which we read the two linear systems
(matrix and right-hand-side vector).
*/
ierr = PetscOptionsGetString(NULL,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");
/*
Open binary file. Note that we use FILE_MODE_READ to indicate
reading from this file.
*/
ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);CHKERRQ(ierr);
/*
Load the matrix; then destroy the viewer.
*/
ierr = MatCreate(PETSC_COMM_WORLD,&A[0]);CHKERRQ(ierr);
ierr = MatLoad(A[0],fd);CHKERRQ(ierr);
ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr);
ierr = MatDuplicate(A[0],MAT_COPY_VALUES,&A[1]);CHKERRQ(ierr);
ierr = MatDuplicate(A[0],MAT_COPY_VALUES,&A[2]);CHKERRQ(ierr);
ierr = MatShift(A[1],1.0);CHKERRQ(ierr);
ierr = MatShift(A[1],2.0);CHKERRQ(ierr);
ierr = MatCreateVecs(A[0],&x,&y);CHKERRQ(ierr);
ierr = VecDuplicate(y,&work);CHKERRQ(ierr);
ierr = VecDuplicate(y,&z);CHKERRQ(ierr);
ierr = VecSet(x,1.0);CHKERRQ(ierr);
ierr = MatMult(A[0],x,z);CHKERRQ(ierr);
ierr = MatMultAdd(A[1],x,z,z);CHKERRQ(ierr);
ierr = MatMultAdd(A[2],x,z,z);CHKERRQ(ierr);
ierr = MatCreateComposite(PETSC_COMM_WORLD,3,A,&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 > 1.e-10) {
ierr = PetscPrintf(PETSC_COMM_WORLD,"Error with composite add %g\n",(double)rnorm);CHKERRQ(ierr);
}
ierr = MatCreateComposite(PETSC_COMM_WORLD,3,A,&B);CHKERRQ(ierr);
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 > 1.e-10) {
ierr = PetscPrintf(PETSC_COMM_WORLD,"Error with composite add after merge %g\n",(double)rnorm);CHKERRQ(ierr);
}
ierr = VecSet(x,1.0);CHKERRQ(ierr);
ierr = MatMult(A[0],x,z);CHKERRQ(ierr);
ierr = MatMult(A[1],z,work);CHKERRQ(ierr);
ierr = MatMult(A[2],work,z);CHKERRQ(ierr);
ierr = MatCreateComposite(PETSC_COMM_WORLD,3,A,&B);CHKERRQ(ierr);
ierr = MatCompositeSetType(B,MAT_COMPOSITE_MULTIPLICATIVE);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 > 1.e-10) {
ierr = PetscPrintf(PETSC_COMM_WORLD,"Error with composite multiplicative %g\n",(double)rnorm);CHKERRQ(ierr);
}
ierr = MatCreateComposite(PETSC_COMM_WORLD,3,A,&B);CHKERRQ(ierr);
ierr = MatCompositeSetType(B,MAT_COMPOSITE_MULTIPLICATIVE);CHKERRQ(ierr);
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 > 1.e-10) {
ierr = PetscPrintf(PETSC_COMM_WORLD,"Error with composite multiplicative after merge %g\n",(double)rnorm);CHKERRQ(ierr);
}
/*
Free work space. All PETSc objects should be destroyed when they
are no longer needed.
*/
ierr = VecDestroy(&x);CHKERRQ(ierr);
ierr = VecDestroy(&y);CHKERRQ(ierr);
ierr = VecDestroy(&work);CHKERRQ(ierr);
ierr = VecDestroy(&z);CHKERRQ(ierr);
ierr = MatDestroy(&A[0]);CHKERRQ(ierr);
ierr = MatDestroy(&A[1]);CHKERRQ(ierr);
//.........这里部分代码省略.........
示例14: main
int main(int argc,char **argv)
{
PetscMPIInt rank,size;
PetscInt M = 14,time_steps = 20,w=1,s=1,localsize,j,i,mybase,myend,globalsize;
PetscErrorCode ierr;
DM da;
Vec global,local;
PetscScalar *globalptr,*localptr;
PetscReal h,k;
ierr = PetscInitialize(&argc,&argv,(char*)0,help);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(NULL,"-M",&M,NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(NULL,"-time",&time_steps,NULL);CHKERRQ(ierr);
/* Set up the array */
ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,M,w,s,NULL,&da);CHKERRQ(ierr);
ierr = DMCreateGlobalVector(da,&global);CHKERRQ(ierr);
ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
/* Make copy of local array for doing updates */
ierr = DMCreateLocalVector(da,&local);CHKERRQ(ierr);
/* determine starting point of each processor */
ierr = VecGetOwnershipRange(global,&mybase,&myend);CHKERRQ(ierr);
/* Initialize the Array */
ierr = VecGetLocalSize (global,&globalsize);CHKERRQ(ierr);
ierr = VecGetArray (global,&globalptr);CHKERRQ(ierr);
for (i=0; i<globalsize; i++) {
j = i + mybase;
globalptr[i] = PetscSinReal((PETSC_PI*j*6)/((PetscReal)M) + 1.2 * PetscSinReal((PETSC_PI*j*2)/((PetscReal)M))) * 4+4;
}
ierr = VecRestoreArray(global,&localptr);CHKERRQ(ierr);
/* Assign Parameters */
h= 1.0/M;
k= h*h/2.2;
ierr = VecGetLocalSize(local,&localsize);CHKERRQ(ierr);
for (j=0; j<time_steps; j++) {
/* Global to Local */
ierr = DMGlobalToLocalBegin(da,global,INSERT_VALUES,local);CHKERRQ(ierr);
ierr = DMGlobalToLocalEnd(da,global,INSERT_VALUES,local);CHKERRQ(ierr);
/*Extract local array */
ierr = VecGetArray(local,&localptr);CHKERRQ(ierr);
ierr = VecGetArray (global,&globalptr);CHKERRQ(ierr);
/* Update Locally - Make array of new values */
/* Note: I don't do anything for the first and last entry */
for (i=1; i< localsize-1; i++) {
globalptr[i-1] = localptr[i] + (k/(h*h)) * (localptr[i+1]-2.0*localptr[i]+localptr[i-1]);
}
ierr = VecRestoreArray (global,&globalptr);CHKERRQ(ierr);
ierr = VecRestoreArray(local,&localptr);CHKERRQ(ierr);
/* View Wave */
/* Set Up Display to Show Heat Graph */
#if defined(PETSC_USE_SOCKET_VIEWER)
ierr = VecView(global,PETSC_VIEWER_SOCKET_WORLD);CHKERRQ(ierr);
#endif
}
ierr = VecDestroy(&local);CHKERRQ(ierr);
ierr = VecDestroy(&global);CHKERRQ(ierr);
ierr = DMDestroy(&da);CHKERRQ(ierr);
ierr = PetscFinalize();
return 0;
}示例15: main
int main(int argc,char **argv)
{
SNES snes; /* nonlinear solver context */
Vec x,r; /* solution, residual vectors */
Mat J; /* Jacobian matrix */
PetscErrorCode ierr;
PetscInt its;
PetscScalar *xx;
SNESConvergedReason reason;
PetscInitialize(&argc,&argv,(char*)0,help);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Create nonlinear solver context
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Create matrix and vector data structures; set corresponding routines
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/*
Create vectors for solution and nonlinear function
*/
ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr);
ierr = VecSetSizes(x,PETSC_DECIDE,2);CHKERRQ(ierr);
ierr = VecSetFromOptions(x);CHKERRQ(ierr);
ierr = VecDuplicate(x,&r);CHKERRQ(ierr);
/*
Create Jacobian matrix data structure
*/
ierr = MatCreate(PETSC_COMM_WORLD,&J);CHKERRQ(ierr);
ierr = MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2);CHKERRQ(ierr);
ierr = MatSetFromOptions(J);CHKERRQ(ierr);
ierr = MatSetUp(J);CHKERRQ(ierr);
/*
Set function evaluation routine and vector.
*/
ierr = SNESSetFunction(snes,r,FormFunction1,NULL);CHKERRQ(ierr);
/*
Set Jacobian matrix data structure and Jacobian evaluation routine
*/
ierr = SNESSetJacobian(snes,J,J,FormJacobian1,NULL);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Customize nonlinear solver; set runtime options
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Evaluate initial guess; then solve nonlinear system
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = VecGetArray(x,&xx);CHKERRQ(ierr);
xx[0] = -1.2; xx[1] = 1.0;
ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr);
/*
Note: The user should initialize the vector, x, with the initial guess
for the nonlinear solver prior to calling SNESSolve(). In particular,
to employ an initial guess of zero, the user should explicitly set
this vector to zero by calling VecSet().
*/
ierr = SNESSolve(snes,NULL,x);CHKERRQ(ierr);
ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr);
ierr = SNESGetConvergedReason(snes,&reason);CHKERRQ(ierr);
/*
Some systems computes a residual that is identically zero, thus converging
due to FNORM_ABS, others converge due to FNORM_RELATIVE. Here, we only
report the reason if the iteration did not converge so that the tests are
reproducible.
*/
ierr = PetscPrintf(PETSC_COMM_WORLD,"%s number of SNES iterations = %D\n\n",reason>0 ? "CONVERGED" : SNESConvergedReasons[reason],its);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Free work space. All PETSc objects should be destroyed when they
are no longer needed.
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = VecDestroy(&x);CHKERRQ(ierr); ierr = VecDestroy(&r);CHKERRQ(ierr);
ierr = MatDestroy(&J);CHKERRQ(ierr); ierr = SNESDestroy(&snes);CHKERRQ(ierr);
ierr = PetscFinalize();
return 0;
}