本文整理汇总了C++中PetscFinalize函数的典型用法代码示例。如果您正苦于以下问题:C++ PetscFinalize函数的具体用法?C++ PetscFinalize怎么用?C++ PetscFinalize使用的例子?那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了PetscFinalize函数的20个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于我们的系统推荐出更棒的C++代码示例。
示例1: main
int main(int argc, char *argv[])
{
PetscViewer viewer;
Vec u;
PetscScalar v;
int VECTOR_GENERATE, VECTOR_READ;
int i, m = 10, rank, size, low, high, ldim, iglobal;
int ierr;
ierr = PetscInitialize(&argc, &argv, NULL, help);CHKERRQ(ierr);
ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);CHKERRQ(ierr);
ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(NULL,NULL, "-m", &m, NULL);CHKERRQ(ierr);
/* PART 1: Generate vector, then write it to Mathematica */
ierr = PetscLogEventRegister("Generate Vector", VEC_CLASSID,&VECTOR_GENERATE);CHKERRQ(ierr);
ierr = PetscLogEventBegin(VECTOR_GENERATE, 0, 0, 0, 0);CHKERRQ(ierr);
/* Generate vector */
ierr = VecCreate(PETSC_COMM_WORLD, &u);CHKERRQ(ierr);
ierr = VecSetSizes(u, PETSC_DECIDE, m);CHKERRQ(ierr);
ierr = VecSetFromOptions(u);CHKERRQ(ierr);
ierr = VecGetOwnershipRange(u, &low, &high);CHKERRQ(ierr);
ierr = VecGetLocalSize(u, &ldim);CHKERRQ(ierr);
for (i = 0; i < ldim; i++) {
iglobal = i + low;
v = (PetscScalar) (i + 100*rank);
ierr = VecSetValues(u, 1, &iglobal, &v, INSERT_VALUES);CHKERRQ(ierr);
}
ierr = VecAssemblyBegin(u);CHKERRQ(ierr);
ierr = VecAssemblyEnd(u);CHKERRQ(ierr);
ierr = VecView(u, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD, "writing vector to Mathematica...\n");CHKERRQ(ierr);
#if 0
ierr = PetscViewerMathematicaOpen(PETSC_COMM_WORLD, 8000, "192.168.119.1", "Connect", &viewer);CHKERRQ(ierr);
ierr = VecView(u, viewer);CHKERRQ(ierr);
#else
ierr = VecView(u, PETSC_VIEWER_MATHEMATICA_WORLD);CHKERRQ(ierr);
#endif
v = 0.0;
ierr = VecSet(u,v);CHKERRQ(ierr);
ierr = PetscLogEventEnd(VECTOR_GENERATE, 0, 0, 0, 0);CHKERRQ(ierr);
/* All processors wait until test vector has been dumped */
ierr = MPI_Barrier(PETSC_COMM_WORLD);CHKERRQ(ierr);
ierr = PetscSleep(10);CHKERRQ(ierr);
/* PART 2: Read in vector in from Mathematica */
ierr = PetscLogEventRegister("Read Vector", VEC_CLASSID,&VECTOR_READ);CHKERRQ(ierr);
ierr = PetscLogEventBegin(VECTOR_READ, 0, 0, 0, 0);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD, "reading vector from Mathematica...\n");CHKERRQ(ierr);
/* Read new vector in binary format */
#if 0
ierr = PetscViewerMathematicaGetVector(viewer, u);CHKERRQ(ierr);
ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
#else
ierr = PetscViewerMathematicaGetVector(PETSC_VIEWER_MATHEMATICA_WORLD, u);CHKERRQ(ierr);
#endif
ierr = PetscLogEventEnd(VECTOR_READ, 0, 0, 0, 0);CHKERRQ(ierr);
ierr = VecView(u, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
/* Free data structures */
ierr = VecDestroy(&u);CHKERRQ(ierr);
ierr = PetscFinalize();
return ierr;
}
开发者ID:tom-klotz,项目名称:petsc,代码行数:69,代码来源:ex15.c
示例2: main
int main(int argc, char* argv[]){
PetscErrorCode ierr;
ierr = PetscInitialize(&argc,&argv,NULL,NULL); CHKERRQ(ierr);
std::map<std::string,bool> commandline;
observeCommand(commandline,"-readLocally");
observeCommand(commandline,"-mshBinary");
observeCommand(commandline,"-outBinary");
observeCommand(commandline,"-cgns");
if(!parseCommand(commandline)){
ierr = PetscFinalize(); CHKERRQ(ierr);
return 0;
}
/******************************************
* START UP
******************************************/
NavierStokesSolver* nsSolver = new NavierStokesSolver;
nsSolver->initSolverParam(); //root only : read param.in and check
nsSolver->readCommand(commandline);
nsSolver->broadcastSolverParam(); //collective fetch param from root
int* elementBuffer = NULL;
double* vertexBuffer = NULL;
int* interfaceBuffer = NULL;
if(!commandline["-readLocally"]){ //transfer geometry through MPI
nsSolver->readAndPartition(getCommand(commandline,"-mshBinary"), getCommand(commandline,"-cgns")); //root only : read msh and partition
nsSolver->broadcastPartitionInfo();
nsSolver->scatterGridFile(&elementBuffer,&vertexBuffer,&interfaceBuffer);//collective
}else{
//read geometry locally
}
//parse the gridfile as original, buffer freeed, boundInfo got;
nsSolver->ReadGridFile(elementBuffer,vertexBuffer,interfaceBuffer);
nsSolver->dataPartition->initPetsc();
//build faces. same sequence as Original;
nsSolver->CreateFaces();
nsSolver->CellFaceInfo();
nsSolver->CheckAndAllocate();
nsSolver->InitFlowField();
/******************************************
* NS_solve
* MAIN CFD
******************************************/
nsSolver->NSSolve();
/******************************************
* Post Process
* VTK, Tecplot, etc.
******************************************/
MPI_Barrier(MPI_COMM_WORLD);
nsSolver->dataPartition->deinit();
delete nsSolver;
//PetscPrintf(MPI_COMM_WORLD,"done\n");
PetscPrintf(MPI_COMM_WORLD,"done\n");
MPI_Barrier(MPI_COMM_WORLD);
ierr = PetscFinalize(); CHKERRQ(ierr);
return 0;
}
开发者ID:Meiyim,项目名称:SYSU_CFD_Parallel,代码行数:68,代码来源:main.cpp
示例3: main
int main(int argc, char *argv[]) {
PetscInitialize(&argc, &argv, NULL, help);
// Get command line options.
Options options;
options.setOptions();
// Get mesh.
Mesh *mesh = Mesh::factory(options);
mesh->read(options);
// Get model.
ExodusModel model(options);
model.initializeParallel();
// Get source.
std::vector<Source*> sources = Source::factory(options);
// Setup reference element.
Square *reference_element = new Acoustic(options);
mesh->setupGlobalDof(reference_element->NumberDofVertex(), reference_element->NumberDofEdge(),
reference_element->NumberDofFace(), reference_element->NumberDofVolume(),
reference_element->NumberDimensions());
mesh->registerFields();
// Clone a list of all local elements.
std::vector<Square *> elements;
for (auto i = 0; i < mesh->NumberElementsLocal(); i++) { elements.push_back(reference_element->clone()); }
// Now things that only local elements are allowed to do.
int element_number = 0;
for (auto &element: elements) {
element->SetLocalElementNumber(element_number++);
element->attachVertexCoordinates(mesh->DistributedMesh());
element->attachSource(sources);
element->interpolateMaterialProperties(model);
element->readOperators();
element->assembleMassMatrix();
element->scatterMassMatrix(mesh);
}
mesh->checkInMassMatrix();
mesh->setUpMovie();
double time = 0;
double timestep = 1e-3;
while (time < 2.0) {
// Pull down fields from global dof.
mesh->checkOutFields();
mesh->zeroFields();
// Compute element-wise terms.
for (auto &element: elements) {
element->SetTime(time);
element->checkOutFields(mesh);
element->computeStiffnessTerm();
element->computeSourceTerm();
element->computeSurfaceTerm();
element->checkInField(mesh);
}
// Sum fields back into global dof.
mesh->checkInFieldsBegin();
mesh->checkInFieldsEnd();
// Invert mass matrix and take time step.
mesh->applyInverseMassMatrix();
mesh->advanceField();
mesh->saveFrame();
time += timestep;
if (!MPI::COMM_WORLD.Get_rank()) std::cout << time << std::endl;
// break;
}
mesh->finalizeMovie();
PetscFinalize();
}
开发者ID:gitter-badger,项目名称:salvus,代码行数:84,代码来源:Main.cpp
示例4: main
//.........这里部分代码省略.........
ierr = MatCreate(PETSC_COMM_WORLD,&user.A);CHKERRQ(ierr);
ierr = MatSetSizes(user.A,PETSC_DECIDE,PETSC_DECIDE,2,2);CHKERRQ(ierr);
ierr = MatSetFromOptions(user.A);CHKERRQ(ierr);
ierr = MatSetUp(user.A);CHKERRQ(ierr);
ierr = MatCreateVecs(user.A,&user.x,NULL);CHKERRQ(ierr);
ierr = MatCreate(PETSC_COMM_WORLD,&user.Jacp);CHKERRQ(ierr);
ierr = MatSetSizes(user.Jacp,PETSC_DECIDE,PETSC_DECIDE,2,1);CHKERRQ(ierr);
ierr = MatSetFromOptions(user.Jacp);CHKERRQ(ierr);
ierr = MatSetUp(user.Jacp);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Create timestepping solver context
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr);
ierr = TSSetType(ts,TSCN);CHKERRQ(ierr);
ierr = TSSetIFunction(ts,NULL,IFunction,&user);CHKERRQ(ierr);
ierr = TSSetIJacobian(ts,user.A,user.A,IJacobian,&user);CHKERRQ(ierr);
ierr = TSSetMaxTime(ts,user.ftime);CHKERRQ(ierr);
ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);CHKERRQ(ierr);
if (monitor) {
ierr = TSMonitorSet(ts,Monitor,&user,NULL);CHKERRQ(ierr);
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Set initial conditions
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = VecGetArray(user.x,&x_ptr);CHKERRQ(ierr);
x_ptr[0] = 2.0; x_ptr[1] = -0.66666654321;
ierr = VecRestoreArray(user.x,&x_ptr);CHKERRQ(ierr);
ierr = TSSetTimeStep(ts,.0001);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Save trajectory of solution so that TSAdjointSolve() may be used
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = TSSetSaveTrajectory(ts);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Set runtime options
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
ierr = TSSolve(ts,user.x);CHKERRQ(ierr);
ierr = TSGetSolveTime(ts,&user.ftime);CHKERRQ(ierr);
ierr = TSGetStepNumber(ts,&user.steps);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Adjoint model starts here
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = MatCreateVecs(user.A,&user.lambda[0],NULL);CHKERRQ(ierr);
/* Set initial conditions for the adjoint integration */
ierr = VecGetArray(user.lambda[0],&y_ptr);CHKERRQ(ierr);
y_ptr[0] = 1.0; y_ptr[1] = 0.0;
ierr = VecRestoreArray(user.lambda[0],&y_ptr);CHKERRQ(ierr);
ierr = MatCreateVecs(user.A,&user.lambda[1],NULL);CHKERRQ(ierr);
ierr = VecGetArray(user.lambda[1],&y_ptr);CHKERRQ(ierr);
y_ptr[0] = 0.0; y_ptr[1] = 1.0;
ierr = VecRestoreArray(user.lambda[1],&y_ptr);CHKERRQ(ierr);
ierr = MatCreateVecs(user.Jacp,&user.mup[0],NULL);CHKERRQ(ierr);
ierr = VecGetArray(user.mup[0],&x_ptr);CHKERRQ(ierr);
x_ptr[0] = 0.0;
ierr = VecRestoreArray(user.mup[0],&x_ptr);CHKERRQ(ierr);
ierr = MatCreateVecs(user.Jacp,&user.mup[1],NULL);CHKERRQ(ierr);
ierr = VecGetArray(user.mup[1],&x_ptr);CHKERRQ(ierr);
x_ptr[0] = 0.0;
ierr = VecRestoreArray(user.mup[1],&x_ptr);CHKERRQ(ierr);
ierr = TSSetCostGradients(ts,2,user.lambda,user.mup);CHKERRQ(ierr);
/* Set RHS JacobianP */
ierr = TSSetRHSJacobianP(ts,user.Jacp,RHSJacobianP,&user);CHKERRQ(ierr);
ierr = TSAdjointSolve(ts);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"\n sensitivity wrt initial conditions: d[y(tf)]/d[y0] d[y(tf)]/d[z0]\n");CHKERRQ(ierr);
ierr = VecView(user.lambda[0],PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"\n sensitivity wrt initial conditions: d[z(tf)]/d[y0] d[z(tf)]/d[z0]\n");CHKERRQ(ierr);
ierr = VecView(user.lambda[1],PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"\n sensitivity wrt parameters: d[y(tf)]/d[mu]\n");CHKERRQ(ierr);
ierr = VecView(user.mup[0],PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"\n sensivitity wrt parameters: d[z(tf)]/d[mu]\n");CHKERRQ(ierr);
ierr = VecView(user.mup[1],PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Free work space. All PETSc objects should be destroyed when they
are no longer needed.
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = MatDestroy(&user.A);CHKERRQ(ierr);
ierr = MatDestroy(&user.Jacp);CHKERRQ(ierr);
ierr = VecDestroy(&user.x);CHKERRQ(ierr);
ierr = VecDestroy(&user.lambda[0]);CHKERRQ(ierr);
ierr = VecDestroy(&user.lambda[1]);CHKERRQ(ierr);
ierr = VecDestroy(&user.mup[0]);CHKERRQ(ierr);
ierr = VecDestroy(&user.mup[1]);CHKERRQ(ierr);
ierr = TSDestroy(&ts);CHKERRQ(ierr);
ierr = PetscFinalize();
return(ierr);
}
开发者ID:firedrakeproject,项目名称:petsc,代码行数:101,代码来源:ex20adj.c
示例5: main
//.........这里部分代码省略.........
SNES snes;
DM da;
Vec x,X,b;
PetscBool youngflg,poissonflg,muflg,lambdaflg,view=PETSC_FALSE,viewline=PETSC_FALSE;
PetscReal poisson=0.2,young=4e4;
char filename[PETSC_MAX_PATH_LEN] = "ex16.vts";
char filename_def[PETSC_MAX_PATH_LEN] = "ex16_def.vts";
ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
ierr = FormElements();CHKERRQ(ierr);
comm = PETSC_COMM_WORLD;
ierr = SNESCreate(comm,&snes);CHKERRQ(ierr);
ierr = DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,11,2,2,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,3,1,NULL,NULL,NULL,&da);CHKERRQ(ierr);
ierr = DMSetFromOptions(da);CHKERRQ(ierr);
ierr = DMSetUp(da);CHKERRQ(ierr);
ierr = SNESSetDM(snes,(DM)da);CHKERRQ(ierr);
ierr = SNESSetNGS(snes,NonlinearGS,&user);CHKERRQ(ierr);
ierr = DMDAGetInfo(da,0,&mx,&my,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);CHKERRQ(ierr);
user.loading = 0.0;
user.arc = PETSC_PI/3.;
user.mu = 4.0;
user.lambda = 1.0;
user.rad = 100.0;
user.height = 3.;
user.width = 1.;
user.ploading = -5e3;
ierr = PetscOptionsGetReal(NULL,NULL,"-arc",&user.arc,NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetReal(NULL,NULL,"-mu",&user.mu,&muflg);CHKERRQ(ierr);
ierr = PetscOptionsGetReal(NULL,NULL,"-lambda",&user.lambda,&lambdaflg);CHKERRQ(ierr);
ierr = PetscOptionsGetReal(NULL,NULL,"-rad",&user.rad,NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetReal(NULL,NULL,"-height",&user.height,NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetReal(NULL,NULL,"-width",&user.width,NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetReal(NULL,NULL,"-loading",&user.loading,NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetReal(NULL,NULL,"-ploading",&user.ploading,NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetReal(NULL,NULL,"-poisson",&poisson,&poissonflg);CHKERRQ(ierr);
ierr = PetscOptionsGetReal(NULL,NULL,"-young",&young,&youngflg);CHKERRQ(ierr);
if ((youngflg || poissonflg) || !(muflg || lambdaflg)) {
/* set the lame' parameters based upon the poisson ratio and young's modulus */
user.lambda = poisson*young / ((1. + poisson)*(1. - 2.*poisson));
user.mu = young/(2.*(1. + poisson));
}
ierr = PetscOptionsGetBool(NULL,NULL,"-view",&view,NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetBool(NULL,NULL,"-view_line",&viewline,NULL);CHKERRQ(ierr);
ierr = DMDASetFieldName(da,0,"x_disp");CHKERRQ(ierr);
ierr = DMDASetFieldName(da,1,"y_disp");CHKERRQ(ierr);
ierr = DMDASetFieldName(da,2,"z_disp");CHKERRQ(ierr);
ierr = DMSetApplicationContext(da,&user);CHKERRQ(ierr);
ierr = DMDASNESSetFunctionLocal(da,INSERT_VALUES,(PetscErrorCode (*)(DMDALocalInfo*,void*,void*,void*))FormFunctionLocal,&user);CHKERRQ(ierr);
ierr = DMDASNESSetJacobianLocal(da,(DMDASNESJacobian)FormJacobianLocal,&user);CHKERRQ(ierr);
ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);
ierr = FormCoordinates(da,&user);CHKERRQ(ierr);
ierr = DMCreateGlobalVector(da,&x);CHKERRQ(ierr);
ierr = DMCreateGlobalVector(da,&b);CHKERRQ(ierr);
ierr = InitialGuess(da,&user,x);CHKERRQ(ierr);
ierr = FormRHS(da,&user,b);CHKERRQ(ierr);
ierr = PetscPrintf(comm,"lambda: %f mu: %f\n",(double)user.lambda,(double)user.mu);CHKERRQ(ierr);
/* show a cross-section of the initial state */
if (viewline) {
ierr = DisplayLine(snes,x);CHKERRQ(ierr);
}
/* get the loaded configuration */
ierr = SNESSolve(snes,b,x);CHKERRQ(ierr);
ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr);
ierr = PetscPrintf(comm,"Number of SNES iterations = %D\n", its);CHKERRQ(ierr);
ierr = SNESGetSolution(snes,&X);CHKERRQ(ierr);
/* show a cross-section of the final state */
if (viewline) {
ierr = DisplayLine(snes,X);CHKERRQ(ierr);
}
if (view) {
PetscViewer viewer;
Vec coords;
ierr = PetscViewerVTKOpen(PETSC_COMM_WORLD,filename,FILE_MODE_WRITE,&viewer);CHKERRQ(ierr);
ierr = VecView(x,viewer);CHKERRQ(ierr);
ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
ierr = DMGetCoordinates(da,&coords);CHKERRQ(ierr);
ierr = VecAXPY(coords,1.0,x);CHKERRQ(ierr);
ierr = PetscViewerVTKOpen(PETSC_COMM_WORLD,filename_def,FILE_MODE_WRITE,&viewer);CHKERRQ(ierr);
ierr = VecView(x,viewer);CHKERRQ(ierr);
ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
}
ierr = VecDestroy(&x);CHKERRQ(ierr);
ierr = VecDestroy(&b);CHKERRQ(ierr);
ierr = DMDestroy(&da);CHKERRQ(ierr);
ierr = SNESDestroy(&snes);CHKERRQ(ierr);
ierr = PetscFinalize();
return ierr;
}
开发者ID:firedrakeproject,项目名称:petsc,代码行数:101,代码来源:ex16.c
示例6: main
int main(int argc,char **argv)
{
PetscErrorCode ierr;
KSP ksp;
PC pc;
Vec x,b;
DM da;
Mat A,Atrans;
PetscInt dof=1,M=-8;
PetscBool flg,trans=PETSC_FALSE;
PetscInitialize(&argc,&argv,(char *)0,help);
ierr = PetscOptionsGetInt(PETSC_NULL,"-dof",&dof,PETSC_NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(PETSC_NULL,"-M",&M,PETSC_NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetBool(PETSC_NULL,"-trans",&trans,PETSC_NULL);CHKERRQ(ierr);
ierr = DMDACreate(PETSC_COMM_WORLD,&da);CHKERRQ(ierr);
ierr = DMDASetDim(da,3);CHKERRQ(ierr);
ierr = DMDASetBoundaryType(da,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE);CHKERRQ(ierr);
ierr = DMDASetStencilType(da,DMDA_STENCIL_STAR);CHKERRQ(ierr);
ierr = DMDASetSizes(da,M,M,M);CHKERRQ(ierr);
ierr = DMDASetNumProcs(da,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
ierr = DMDASetDof(da,dof);CHKERRQ(ierr);
ierr = DMDASetStencilWidth(da,1);CHKERRQ(ierr);
ierr = DMDASetOwnershipRanges(da,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
ierr = DMSetFromOptions(da);CHKERRQ(ierr);
ierr = DMSetUp(da);CHKERRQ(ierr);
ierr = DMCreateGlobalVector(da,&x);CHKERRQ(ierr);
ierr = DMCreateGlobalVector(da,&b);CHKERRQ(ierr);
ierr = ComputeRHS(da,b);CHKERRQ(ierr);
ierr = DMCreateMatrix(da,MATBAIJ,&A);CHKERRQ(ierr);
ierr = ComputeMatrix(da,A);CHKERRQ(ierr);
/* A is non-symmetric. Make A = 0.5*(A + Atrans) symmetric for testing icc and cholesky */
ierr = MatTranspose(A,MAT_INITIAL_MATRIX,&Atrans);CHKERRQ(ierr);
ierr = MatAXPY(A,1.0,Atrans,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
ierr = MatScale(A,0.5);CHKERRQ(ierr);
ierr = MatDestroy(&Atrans);CHKERRQ(ierr);
/* Test sbaij matrix */
flg = PETSC_FALSE;
ierr = PetscOptionsGetBool(PETSC_NULL, "-test_sbaij1", &flg,PETSC_NULL);CHKERRQ(ierr);
if (flg){
Mat sA;
PetscBool issymm;
ierr = MatIsTranspose(A,A,0.0,&issymm);CHKERRQ(ierr);
if (issymm) {
ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr);
} else {
printf("Warning: A is non-symmetric\n");
}
ierr = MatConvert(A,MATSBAIJ,MAT_INITIAL_MATRIX,&sA);CHKERRQ(ierr);
ierr = MatDestroy(&A);CHKERRQ(ierr);
A = sA;
}
ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
ierr = KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
ierr = PCSetDM(pc,(DM)da);CHKERRQ(ierr);
if (trans) {
ierr = KSPSolveTranspose(ksp,b,x);CHKERRQ(ierr);
} else {
ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
}
/* check final residual */
flg = PETSC_FALSE;
ierr = PetscOptionsGetBool(PETSC_NULL, "-check_final_residual", &flg,PETSC_NULL);CHKERRQ(ierr);
if (flg){
Vec b1;
PetscReal norm;
ierr = KSPGetSolution(ksp,&x);CHKERRQ(ierr);
ierr = VecDuplicate(b,&b1);CHKERRQ(ierr);
ierr = MatMult(A,x,b1);CHKERRQ(ierr);
ierr = VecAXPY(b1,-1.0,b);CHKERRQ(ierr);
ierr = VecNorm(b1,NORM_2,&norm);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"Final residual %g\n",norm);CHKERRQ(ierr);
ierr = VecDestroy(&b1);CHKERRQ(ierr);
}
ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
ierr = VecDestroy(&x);CHKERRQ(ierr);
ierr = VecDestroy(&b);CHKERRQ(ierr);
ierr = MatDestroy(&A);CHKERRQ(ierr);
ierr = DMDestroy(&da);CHKERRQ(ierr);
ierr = PetscFinalize();
return 0;
}
开发者ID:Kun-Qu,项目名称:petsc,代码行数:93,代码来源:ex32.c
示例7: main
int main(int argc,char **args)
{
Mat C;
int i,m = 5,rank,size,N,start,end,M;
int ierr,idx[4];
PetscScalar Ke[16];
PetscReal h;
Vec u,b;
KSP ksp;
MatNullSpace nullsp;
ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
ierr = PetscOptionsGetInt(NULL,NULL,"-m",&m,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);
ierr = MatSetUp(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(b,1.0);CHKERRQ(ierr);
ierr = VecSetValue(b,0,1.2,ADD_VALUES);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);CHKERRQ(ierr);
ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
ierr = KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);CHKERRQ(ierr);
ierr = MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,NULL,&nullsp);CHKERRQ(ierr);
/*
The KSP solver will remove this nullspace from the solution at each iteration
*/
ierr = MatSetNullSpace(C,nullsp);CHKERRQ(ierr);
/*
The KSP solver will remove from the right hand side any portion in this nullspace, thus making the linear system consistent.
*/
ierr = MatSetTransposeNullSpace(C,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 ierr;
}
开发者ID:tom-klotz,项目名称:petsc,代码行数:79,代码来源:ex20.c
示例8: main
int main(int argc,char **args)
{
Mat A;
PetscErrorCode ierr;
PetscMPIInt rank,size;
PetscInt *ia,*ja;
MatPartitioning part;
IS is,isn;
PetscInitialize(&argc,&args,(char*)0,help);
ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
if (size != 4) SETERRQ(PETSC_COMM_WORLD,1,"Must run with 4 processors");
ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
ierr = PetscMalloc(5*sizeof(PetscInt),&ia);CHKERRQ(ierr);
ierr = PetscMalloc(16*sizeof(PetscInt),&ja);CHKERRQ(ierr);
if (!rank) {
ja[0] = 1; ja[1] = 4; ja[2] = 0; ja[3] = 2; ja[4] = 5; ja[5] = 1; ja[6] = 3; ja[7] = 6;
ja[8] = 2; ja[9] = 7;
ia[0] = 0; ia[1] = 2; ia[2] = 5; ia[3] = 8; ia[4] = 10;
} else if (rank == 1) {
ja[0] = 0; ja[1] = 5; ja[2] = 8; ja[3] = 1; ja[4] = 4; ja[5] = 6; ja[6] = 9; ja[7] = 2;
ja[8] = 5; ja[9] = 7; ja[10] = 10; ja[11] = 3; ja[12] = 6; ja[13] = 11;
ia[0] = 0; ia[1] = 3; ia[2] = 7; ia[3] = 11; ia[4] = 14;
} else if (rank == 2) {
ja[0] = 4; ja[1] = 9; ja[2] = 12; ja[3] = 5; ja[4] = 8; ja[5] = 10; ja[6] = 13; ja[7] = 6;
ja[8] = 9; ja[9] = 11; ja[10] = 14; ja[11] = 7; ja[12] = 10; ja[13] = 15;
ia[0] = 0; ia[1] = 3; ia[2] = 7; ia[3] = 11; ia[4] = 14;
} else {
ja[0] = 8; ja[1] = 13; ja[2] = 9; ja[3] = 12; ja[4] = 14; ja[5] = 10; ja[6] = 13; ja[7] = 15;
ja[8] = 11; ja[9] = 14;
ia[0] = 0; ia[1] = 2; ia[2] = 5; ia[3] = 8; ia[4] = 10;
}
ierr = MatCreateMPIAdj(PETSC_COMM_WORLD,4,16,ia,ja,NULL,&A);CHKERRQ(ierr);
ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
/*
Partition the graph of the matrix
*/
ierr = MatPartitioningCreate(PETSC_COMM_WORLD,&part);CHKERRQ(ierr);
ierr = MatPartitioningSetAdjacency(part,A);CHKERRQ(ierr);
ierr = MatPartitioningSetFromOptions(part);CHKERRQ(ierr);
/* get new processor owner number of each vertex */
ierr = MatPartitioningApply(part,&is);CHKERRQ(ierr);
/* get new global number of each old global number */
ierr = ISPartitioningToNumbering(is,&isn);CHKERRQ(ierr);
ierr = ISView(isn,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
ierr = ISDestroy(&is);CHKERRQ(ierr);
ierr = ISDestroy(&isn);CHKERRQ(ierr);
ierr = MatPartitioningDestroy(&part);CHKERRQ(ierr);
/*
Free work space. All PETSc objects should be destroyed when they
are no longer needed.
*/
ierr = MatDestroy(&A);CHKERRQ(ierr);
ierr = PetscFinalize();
return 0;
}
开发者ID:feelpp,项目名称:debian-petsc,代码行数:63,代码来源:ex80.c
示例9: indices
//.........这里部分代码省略.........
if (!flg) fromFirst = 2;
ierr = PetscOptionsGetInt(NULL,NULL,"-fromStep",&fromStep,&flg);CHKERRQ(ierr);
if (!flg) fromStep = 2;
if (n>m) {
ierr = PetscPrintf(PETSC_COMM_WORLD,"The vector sizes are %D. The number of elements being scattered is %D\n",m,n);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"Adjust the parameters such that m>=n\n");CHKERRQ(ierr);
} else if (toFirst+(n-1)*toStep >=m) {
ierr = PetscPrintf(PETSC_COMM_WORLD,"The vector sizes are %D. The number of elements being scattered is %D\n",m,n);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"For the Strided Scatter, toFirst=%D and toStep=%D.\n",toFirst,toStep);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"This produces an index (toFirst+(n-1)*toStep)>=m\n");CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"Adjust the parameterrs accordingly with -m, -n, -toFirst, or -toStep\n");CHKERRQ(ierr);
} else if (fromFirst+(n-1)*fromStep>=m) {
ierr = PetscPrintf(PETSC_COMM_WORLD,"The vector sizes are %D. The number of elements being scattered is %D\n",m,n);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"For the Strided Scatter, fromFirst=%D and fromStep=%D.\n",fromFirst,toStep);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"This produces an index (fromFirst+(n-1)*fromStep)>=m\n");CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"Adjust the parameterrs accordingly with -m, -n, -fromFirst, or -fromStep\n");CHKERRQ(ierr);
} else {
ierr = PetscPrintf(PETSC_COMM_WORLD,"m=%D\tn=%D\tfromFirst=%D\tfromStep=%D\ttoFirst=%D\ttoStep=%D\n",m,n,fromFirst,fromStep,toFirst,toStep);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"fromFirst+(n-1)*fromStep=%D\ttoFirst+(n-1)*toStep=%D\n",fromFirst+(n-1)*fromStep,toFirst+(n-1)*toStep);CHKERRQ(ierr);
/* Build the vectors */
ierr = VecCreate(PETSC_COMM_WORLD,&Y);CHKERRQ(ierr);
ierr = VecSetSizes(Y,m,PETSC_DECIDE);CHKERRQ(ierr);
ierr = VecCreate(PETSC_COMM_WORLD,&X);CHKERRQ(ierr);
ierr = VecSetSizes(X,m,PETSC_DECIDE);CHKERRQ(ierr);
ierr = VecSetFromOptions(Y);CHKERRQ(ierr);
ierr = VecSetFromOptions(X);CHKERRQ(ierr);
ierr = VecSet(X,2.0);CHKERRQ(ierr);
ierr = VecSet(Y,1.0);CHKERRQ(ierr);
/* Build the strided index sets */
ierr = ISCreate(PETSC_COMM_WORLD,&toISStrided);CHKERRQ(ierr);
ierr = ISCreate(PETSC_COMM_WORLD,&fromISStrided);CHKERRQ(ierr);
ierr = ISSetType(toISStrided, ISSTRIDE);CHKERRQ(ierr);
ierr = ISSetType(fromISStrided, ISSTRIDE);CHKERRQ(ierr);
ierr = ISStrideSetStride(fromISStrided,n,fromFirst,fromStep);CHKERRQ(ierr);
ierr = ISStrideSetStride(toISStrided,n,toFirst,toStep);CHKERRQ(ierr);
/* Build the general index sets */
ierr = PetscMalloc1(n,&idx);CHKERRQ(ierr);
ierr = PetscMalloc1(n,&idy);CHKERRQ(ierr);
for (i=0; i<n; i++) {
idx[i] = i % m;
idy[i] = (i+m) % m;
}
n1 = n;
n2 = n;
ierr = PetscSortRemoveDupsInt(&n1,idx);CHKERRQ(ierr);
ierr = PetscSortRemoveDupsInt(&n2,idy);CHKERRQ(ierr);
ierr = ISCreateGeneral(PETSC_COMM_WORLD,n1,idx,PETSC_COPY_VALUES,&toISGeneral);CHKERRQ(ierr);
ierr = ISCreateGeneral(PETSC_COMM_WORLD,n2,idy,PETSC_COPY_VALUES,&fromISGeneral);CHKERRQ(ierr);
/* set the mode and the insert/add parameter */
mode = SCATTER_FORWARD;
addv = ADD_VALUES;
/* VecScatter : Seq Strided to Seq Strided */
ierr = VecScatterCreate(X,fromISStrided,Y,toISStrided,&vscatSStoSS);CHKERRQ(ierr);
ierr = VecScatterBegin(vscatSStoSS,X,Y,addv,mode);CHKERRQ(ierr);
ierr = VecScatterEnd(vscatSStoSS,X,Y,addv,mode);CHKERRQ(ierr);
ierr = VecScatterDestroy(&vscatSStoSS);CHKERRQ(ierr);
/* VecScatter : Seq General to Seq Strided */
ierr = VecScatterCreate(Y,fromISGeneral,X,toISStrided,&vscatSGtoSS);CHKERRQ(ierr);
ierr = VecScatterBegin(vscatSGtoSS,Y,X,addv,mode);CHKERRQ(ierr);
ierr = VecScatterEnd(vscatSGtoSS,Y,X,addv,mode);CHKERRQ(ierr);
ierr = VecScatterDestroy(&vscatSGtoSS);CHKERRQ(ierr);
/* VecScatter : Seq General to Seq General */
ierr = VecScatterCreate(X,fromISGeneral,Y,toISGeneral,&vscatSGtoSG);CHKERRQ(ierr);
ierr = VecScatterBegin(vscatSGtoSG,X,Y,addv,mode);CHKERRQ(ierr);
ierr = VecScatterEnd(vscatSGtoSG,X,Y,addv,mode);CHKERRQ(ierr);
ierr = VecScatterDestroy(&vscatSGtoSG);CHKERRQ(ierr);
/* VecScatter : Seq Strided to Seq General */
ierr = VecScatterCreate(Y,fromISStrided,X,toISGeneral,&vscatSStoSG);CHKERRQ(ierr);
ierr = VecScatterBegin(vscatSStoSG,Y,X,addv,mode);CHKERRQ(ierr);
ierr = VecScatterEnd(vscatSStoSG,Y,X,addv,mode);CHKERRQ(ierr);
ierr = VecScatterDestroy(&vscatSStoSG);CHKERRQ(ierr);
/* view the results */
ierr = VecView(Y,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
/* Cleanup */
ierr = VecDestroy(&X);CHKERRQ(ierr);
ierr = VecDestroy(&Y);CHKERRQ(ierr);
ierr = ISDestroy(&toISStrided);CHKERRQ(ierr);
ierr = ISDestroy(&fromISStrided);CHKERRQ(ierr);
ierr = ISDestroy(&toISGeneral);CHKERRQ(ierr);
ierr = ISDestroy(&fromISGeneral);CHKERRQ(ierr);
ierr = PetscFree(idx);CHKERRQ(ierr);
ierr = PetscFree(idy);CHKERRQ(ierr);
}
ierr = PetscFinalize();
return ierr;
}
开发者ID:petsc,项目名称:petsc,代码行数:101,代码来源:ex44.c
示例10: main
Example: mpiexec -n <np> ./ex130 -f <matrix binary file> -mat_solver_type 1 -mat_superlu_equil \n\n";
#include <petscmat.h>
int main(int argc,char **args)
{
Mat A,F;
Vec u,x,b;
PetscErrorCode ierr;
PetscMPIInt rank,size;
PetscInt m,n,nfact,ipack=0;
PetscReal norm,tol=1.e-12,Anorm;
IS perm,iperm;
MatFactorInfo info;
PetscBool flg,testMatSolve=PETSC_TRUE;
PetscViewer fd; /* viewer */
char file[PETSC_MAX_PATH_LEN]; /* input file name */
ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);CHKERRQ(ierr);
ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size);CHKERRQ(ierr);
/* Determine file from which we read the matrix A */
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");
/* Load matrix A */
ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);CHKERRQ(ierr);
ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
ierr = MatLoad(A,fd);CHKERRQ(ierr);
ierr = VecCreate(PETSC_COMM_WORLD,&b);CHKERRQ(ierr);
ierr = VecLoad(b,fd);CHKERRQ(ierr);
ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr);
ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr);
if (m != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%d, %d)", m, n);
ierr = MatNorm(A,NORM_INFINITY,&Anorm);CHKERRQ(ierr);
/* Create vectors */
ierr = VecDuplicate(b,&x);CHKERRQ(ierr);
ierr = VecDuplicate(x,&u);CHKERRQ(ierr); /* save the true solution */
/* Test LU Factorization */
ierr = MatGetOrdering(A,MATORDERINGNATURAL,&perm,&iperm);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(NULL,NULL,"-mat_solver_type",&ipack,NULL);CHKERRQ(ierr);
switch (ipack) {
case 1:
#if defined(PETSC_HAVE_SUPERLU)
if (!rank) printf(" SUPERLU LU:\n");
ierr = MatGetFactor(A,MATSOLVERSUPERLU,MAT_FACTOR_LU,&F);CHKERRQ(ierr);
break;
#endif
case 2:
#if defined(PETSC_HAVE_MUMPS)
if (!rank) printf(" MUMPS LU:\n");
ierr = MatGetFactor(A,MATSOLVERMUMPS,MAT_FACTOR_LU,&F);CHKERRQ(ierr);
{
/* test mumps options */
PetscInt icntl_7 = 5;
ierr = MatMumpsSetIcntl(F,7,icntl_7);CHKERRQ(ierr);
}
break;
#endif
default:
if (!rank) printf(" PETSC LU:\n");
ierr = MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_LU,&F);CHKERRQ(ierr);
}
info.fill = 5.0;
ierr = MatLUFactorSymbolic(F,A,perm,iperm,&info);CHKERRQ(ierr);
for (nfact = 0; nfact < 1; nfact++) {
if (!rank) printf(" %d-the LU numfactorization \n",nfact);
ierr = MatLUFactorNumeric(F,A,&info);CHKERRQ(ierr);
/* Test MatSolve() */
if (testMatSolve) {
ierr = MatSolve(F,b,x);CHKERRQ(ierr);
/* Check the residual */
ierr = MatMult(A,x,u);CHKERRQ(ierr);
ierr = VecAXPY(u,-1.0,b);CHKERRQ(ierr);
ierr = VecNorm(u,NORM_INFINITY,&norm);CHKERRQ(ierr);
if (norm > tol) {
if (!rank) {
ierr = PetscPrintf(PETSC_COMM_SELF,"MatSolve: rel residual %g/%g = %g, LU numfact %d\n",norm,Anorm,norm/Anorm,nfact);CHKERRQ(ierr);
}
}
}
}
/* Free data structures */
ierr = MatDestroy(&A);CHKERRQ(ierr);
ierr = MatDestroy(&F);CHKERRQ(ierr);
ierr = ISDestroy(&perm);CHKERRQ(ierr);
ierr = ISDestroy(&iperm);CHKERRQ(ierr);
ierr = VecDestroy(&x);CHKERRQ(ierr);
ierr = VecDestroy(&b);CHKERRQ(ierr);
ierr = VecDestroy(&u);CHKERRQ(ierr);
ierr = PetscFinalize();
//.........这里部分代码省略.........
开发者ID:firedrakeproject,项目名称:petsc,代码行数:101,代码来源:ex130.c
示例11: main
int main(int argc, char** args)
{
int ntdim = 31;
double* t = (double*)malloc(ntdim*sizeof(double));
for (int i = 0; i<ntdim; i++)
{
t[i] = 0.1*i;
}
int nsdim = 2244;
double* u0 = (double*)malloc(ntdim*nsdim*sizeof(double));
std::string line;
std::ifstream myfile ("vortex_51.txt");
double a;
int linenum = -1;
if(myfile.is_open())
{
while(! myfile.eof())
{
getline(myfile,line);
linenum ++;
if (linenum >= nsdim)
{ break; }
std::istringstream iss(line);
for (int i = 0; i<ntdim; i++)
{
iss >> a;
u0[i*nsdim+linenum] = a;
}
}
}
double* s;
double alpha = 10;
lssSolver A(t,u0,s,ntdim,nsdim,alpha);
MPI_Init(&argc,&args);
PetscInitialize(&argc,&args,(char*)0,NULL);
double* uf = A.lss();
std::ofstream outfile ("uf_31.txt");
if(outfile.is_open())
{
for(int i = 0 ; i < ntdim; i++)
{
for (int j = 0; j < nsdim; j++)
{
outfile << uf[i*nsdim+j] << ", ";
}
outfile << std::endl;
}
}
PetscFinalize();
MPI_Finalize();
return 0;
}
开发者ID:hjbae,项目名称:MTLSS,代码行数:75,代码来源:main.cpp
示例12: main
PetscInt main(PetscInt argc,char **args)
{
PetscErrorCode ierr;
PetscMPIInt rank,size;
PetscInt N0=2048,N1=2048,N2=3,N3=5,N4=5,N=N0*N1;
PetscRandom rdm;
PetscReal enorm;
Vec x,y,z,input,output;
Mat A;
PetscInt DIM, dim[5],vsize;
PetscReal fac;
PetscScalar one=1,two=2,three=3;
ierr = PetscInitialize(&argc,&args,(char*)0,help);CHKERRQ(ierr);
#if defined(PETSC_USE_COMPLEX)
SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP, "This example requires real numbers");
#endif
ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size);CHKERRQ(ierr);
ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);CHKERRQ(ierr);
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 = VecSet(input,one);CHKERRQ(ierr); */
/* ierr = VecSetValue(input,1,two,INSERT_VALUES);CHKERRQ(ierr); */
/* ierr = VecSetValue(input,2,three,INSERT_VALUES);CHKERRQ(ierr); */
/* ierr = VecSetValue(input,3,three,INSERT_VALUES);CHKERRQ(ierr); */
ierr = VecSetRandom(input,rdm);CHKERRQ(ierr);
/* ierr = VecSetRandom(input,rdm);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 = MatGetVecsFFTW(A,&x,&y,&z);CHKERRQ(ierr);
/* ierr = MatGetVecs(A,&x,&y);CHKERRQ(ierr); */
/* ierr = MatGetVecs(A,&z,NULL);CHKERRQ(ierr); */
ierr = VecGetSize(x,&vsize);CHKERRQ(ierr);
printf("The vector size of input from the main routine is %d\n",vsize);
ierr = VecGetSize(z,&vsize);CHKERRQ(ierr);
printf("The vector size of output from the main routine is %d\n",vsize);
ierr = InputTransformFFT(A,input,x);CHKERRQ(ierr);
ierr = MatMult(A,x,y);CHKERRQ(ierr);
ierr = VecAssemblyBegin(y);CHKERRQ(ierr);
ierr = VecAssemblyEnd(y);CHKERRQ(ierr);
ierr = VecView(y,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
ierr = MatMultTranspose(A,y,z);CHKERRQ(ierr);
ierr = OutputTransformFFT(A,z,output);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(&output);CHKERRQ(ierr);
ierr = VecDestroy(&input);CHKERRQ(ierr);
ierr = VecDestroy(&x);CHKERRQ(ierr);
ierr = VecDestroy(&y);CHKERRQ(ierr);
ierr = VecDestroy(&z);CHKERRQ(ierr);
ierr = MatDestroy(&A);CHKERRQ(ierr);
ierr = PetscRandomDestroy(&rdm);CHKERRQ(ierr);
PetscFinalize();
return 0;
}
开发者ID:00liujj,项目名称:petsc,代码行数:84,代码来源:ex157.c
示例13: main
//.........这里部分代码省略.........
if (n == PETSC_DECIDE) SETERRQ(PETSC_COMM_WORLD,1,"Must set -n option with -distribute option");
ierr = PetscMalloc(n*sizeof(PetscInt),&ly);CHKERRQ(ierr);
for (i=0; i<n-1; i++) { ly[i] = 2;}
ly[n-1] = N - 2*(n-1);
}
/* Create distributed array and get vectors */
ierr = DMDACreate2d(PETSC_COMM_WORLD,bx,by,st,M,N,m,n,w,s,lx,ly,&da);CHKERRQ(ierr);
ierr = PetscFree(lx);CHKERRQ(ierr);
ierr = PetscFree(ly);CHKERRQ(ierr);
ierr = DMView(da,viewer);CHKERRQ(ierr);
ierr = DMCreateGlobalVector(da,&global);CHKERRQ(ierr);
ierr = DMCreateLocalVector(da,&local);CHKERRQ(ierr);
/* Set global vector; send ghost points to local vectors */
value = 1;
ierr = VecSet(global,value);CHKERRQ(ierr);
ierr = DMGlobalToLocalBegin(da,global,INSERT_VALUES,local);CHKERRQ(ierr);
ierr = DMGlobalToLocalEnd(da,global,INSERT_VALUES,local);CHKERRQ(ierr);
/* Scale local vectors according to processor rank; pass to global vector */
ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
value = rank;
ierr = VecScale(local,value);CHKERRQ(ierr);
ierr = DMLocalToGlobalBegin(da,local,INSERT_VALUES,global);CHKERRQ(ierr);
ierr = DMLocalToGlobalEnd(da,local,INSERT_VALUES,global);CHKERRQ(ierr);
if (!testorder) { /* turn off printing when testing ordering mappings */
ierr = PetscPrintf (PETSC_COMM_WORLD,"\nGlobal Vectors:\n");CHKERRQ(ierr);
ierr = VecView(global,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
ierr = PetscPrintf (PETSC_COMM_WORLD,"\n\n");CHKERRQ(ierr);
}
/* Send ghost points to local vectors */
ierr = DMGlobalToLocalBegin(da,global,INSERT_VALUES,local);CHKERRQ(ierr);
ierr = DMGlobalToLocalEnd(da,global,INSERT_VALUES,local);CHKERRQ(ierr);
flg = PETSC_FALSE;
ierr = PetscOptionsGetBool(PETSC_NULL,"-local_print",&flg,PETSC_NULL);CHKERRQ(ierr);
if (flg) {
PetscViewer sviewer;
ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"\nLocal Vector: processor %d\n",rank);CHKERRQ(ierr);
ierr = PetscViewerGetSingleton(PETSC_VIEWER_STDOUT_WORLD,&sviewer);CHKERRQ(ierr);
ierr = VecView(local,sviewer);CHKERRQ(ierr);
ierr = PetscViewerRestoreSingleton(PETSC_VIEWER_STDOUT_WORLD,&sviewer);CHKERRQ(ierr);
}
/* Tests mappings betweeen application/PETSc orderings */
if (testorder) {
ierr = DMDAGetGhostCorners(da,&Xs,&Ys,PETSC_NULL,&Xm,&Ym,PETSC_NULL);CHKERRQ(ierr);
ierr = DMDAGetGlobalIndices(da,&nloc,<og);CHKERRQ(ierr);
ierr = DMDAGetAO(da,&ao);CHKERRQ(ierr);
ierr = PetscMalloc(nloc*sizeof(PetscInt),&iglobal);CHKERRQ(ierr);
/* Set iglobal to be global indices for each processor's local and ghost nodes,
using the DMDA ordering of grid points */
kk = 0;
for (j=Ys; j<Ys+Ym; j++) {
for (i=Xs; i<Xs+Xm; i++) {
iloc = w*((j-Ys)*Xm + i-Xs);
for (l=0; l<w; l++) {
iglobal[kk++] = ltog[iloc+l];
}
}
}
/* Map this to the application ordering (which for DMDAs is just the natural ordering
that would be used for 1 processor, numbering most rapidly by x, then y) */
ierr = AOPetscToApplication(ao,nloc,iglobal);CHKERRQ(ierr);
/* Then map the application ordering back to the PETSc DMDA ordering */
ierr = AOApplicationToPetsc(ao,nloc,iglobal);CHKERRQ(ierr);
/* Verify the mappings */
kk=0;
for (j=Ys; j<Ys+Ym; j++) {
for (i=Xs; i<Xs+Xm; i++) {
iloc = w*((j-Ys)*Xm + i-Xs);
for (l=0; l<w; l++) {
if (iglobal[kk] != ltog[iloc+l]) {
ierr = PetscFPrintf(PETSC_COMM_SELF,stdout,"[%d] Problem with mapping: j=%D, i=%D, l=%D, petsc1=%D, petsc2=%D\n",
rank,j,i,l,ltog[iloc+l],iglobal[kk]);}
kk++;
}
}
}
ierr = PetscFree(iglobal);CHKERRQ(ierr);
}
/* Free memory */
ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
ierr = VecDestroy(&local);CHKERRQ(ierr);
ierr = VecDestroy(&global);CHKERRQ(ierr);
ierr = DMDestroy(&da);CHKERRQ(ierr);
ierr = PetscFinalize();
return 0;
}
开发者ID:erdc-cm,项目名称:petsc-dev,代码行数:101,代码来源:ex4.c
|
请发表评论