本文整理汇总了C++中KSPCreate函数的典型用法代码示例。如果您正苦于以下问题：C++ KSPCreate函数的具体用法？C++ KSPCreate怎么用？C++ KSPCreate使用的例子？那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了KSPCreate函数的20个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于我们的系统推荐出更棒的C++代码示例。
示例1: petsc_solve
int petsc_solve(int n_, double complex *A_, double complex *b_, double complex *x_)
{
Vec x, b;
Mat A;
KSP ksp;
PC pc;
PetscReal norm, tol=1.e-14;
PetscErrorCode ierr;
PetscInt i, j, n = n_, col[n_], its;
PetscScalar neg_one = -1.0, one = 1.0, value[n_], *x_array;
ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject) x, "Solution");CHKERRQ(ierr);
ierr = VecSetSizes(x,PETSC_DECIDE,n);CHKERRQ(ierr);
ierr = VecSetFromOptions(x);CHKERRQ(ierr);
ierr = VecDuplicate(x,&b);CHKERRQ(ierr);
ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr);
ierr = MatSetFromOptions(A);CHKERRQ(ierr);
ierr = MatSetOption(A,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);CHKERRQ(ierr);
ierr = MatSetUp(A);CHKERRQ(ierr);
for (i=0; i<n; i++) col[i] = i;
printf(" Converting matrix to PETSc\n");
ierr = MatSetValues(A,n,col,n,col,A_,INSERT_VALUES);CHKERRQ(ierr);
printf(" Done\n");
ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
for (i=0; i<n; i++) value[i] = b_[i];
ierr = VecSetValues(b, n, col, value, INSERT_VALUES);CHKERRQ(ierr);
ierr = VecAssemblyBegin(b);CHKERRQ(ierr);
ierr = VecAssemblyEnd(b);CHKERRQ(ierr);
ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
ierr = KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
ierr = PCSetType(pc,PCJACOBI);CHKERRQ(ierr);
ierr = KSPSetTolerances(ksp,1.e-5,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr);
ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
// For a full list, see:
// http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/KSP/KSPType.html
ierr = KSPSetType(ksp, KSPCGS);CHKERRQ(ierr);
printf(" Solving...\n");
ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
printf(" Done\n");
// ierr = KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
ierr = VecGetArray(x, &x_array);CHKERRQ(ierr);
for (i=0; i<n; i++) x_[i] = x_array[i];
ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD, "Iterations %D\n", its);CHKERRQ(ierr);
ierr = VecDestroy(&x);CHKERRQ(ierr);
ierr = VecDestroy(&b);CHKERRQ(ierr); ierr = MatDestroy(&A);CHKERRQ(ierr);
ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
return 0;
}
开发者ID:certik，项目名称:hfsolver，代码行数:66，代码来源:c_petsc.c
示例2: TaylorGalerkinStepIIMassEnergy
/* Notice that this requires the previous momentum solution.
The element stiffness matrix for the identity in linear elements is
1 /2 1 1\
- |1 2 1|
12 \1 1 2/
no matter what the shape of the triangle. */
PetscErrorCode TaylorGalerkinStepIIMassEnergy(DM da, UserContext *user)
{
MPI_Comm comm;
Mat mat;
Vec rhs_m, rhs_e;
PetscScalar identity[9] = {0.16666666667, 0.08333333333, 0.08333333333,
0.08333333333, 0.16666666667, 0.08333333333,
0.08333333333, 0.08333333333, 0.16666666667};
PetscScalar *u_n, *v_n, *p_n, *t_n, *mu_n, *kappa_n;
PetscScalar *rho_n, *rho_u_n, *rho_v_n, *rho_e_n;
PetscScalar *u_phi, *v_phi;
PetscScalar *rho_u_np1, *rho_v_np1;
PetscInt idx[3];
PetscScalar psi_x[3], psi_y[3];
PetscScalar values_m[3];
PetscScalar values_e[3];
PetscScalar phi = user->phi;
PetscScalar mu, kappa, tau_xx, tau_xy, tau_yy, q_x, q_y;
PetscReal hx, hy, area;
KSP ksp;
const PetscInt *necon;
PetscInt j, k, e, ne, nc, mx, my;
PetscErrorCode ierr;
PetscFunctionBeginUser;
ierr = PetscObjectGetComm((PetscObject) da, &comm);CHKERRQ(ierr);
ierr = DMSetMatType(da,MATAIJ);CHKERRQ(ierr);
ierr = DMCreateMatrix(da, &mat);CHKERRQ(ierr);
ierr = MatSetOption(mat,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr);
ierr = DMGetGlobalVector(da, &rhs_m);CHKERRQ(ierr);
ierr = DMGetGlobalVector(da, &rhs_e);CHKERRQ(ierr);
ierr = KSPCreate(comm, &ksp);CHKERRQ(ierr);
ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
ierr = DMDAGetInfo(da, 0, &mx, &my, 0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
hx = 1.0 / (PetscReal)(mx-1);
hy = 1.0 / (PetscReal)(my-1);
area = 0.5*hx*hy;
ierr = VecGetArray(user->sol_n.u, &u_n);CHKERRQ(ierr);
ierr = VecGetArray(user->sol_n.v, &v_n);CHKERRQ(ierr);
ierr = VecGetArray(user->sol_n.p, &p_n);CHKERRQ(ierr);
ierr = VecGetArray(user->sol_n.t, &t_n);CHKERRQ(ierr);
ierr = VecGetArray(user->mu, &mu_n);CHKERRQ(ierr);
ierr = VecGetArray(user->kappa, &kappa_n);CHKERRQ(ierr);
ierr = VecGetArray(user->sol_n.rho, &rho_n);CHKERRQ(ierr);
ierr = VecGetArray(user->sol_n.rho_u, &rho_u_n);CHKERRQ(ierr);
ierr = VecGetArray(user->sol_n.rho_v, &rho_v_n);CHKERRQ(ierr);
ierr = VecGetArray(user->sol_n.rho_e, &rho_e_n);CHKERRQ(ierr);
ierr = VecGetArray(user->sol_phi.u, &u_phi);CHKERRQ(ierr);
ierr = VecGetArray(user->sol_phi.v, &v_phi);CHKERRQ(ierr);
ierr = VecGetArray(user->sol_np1.rho_u, &rho_u_np1);CHKERRQ(ierr);
ierr = VecGetArray(user->sol_np1.rho_v, &rho_v_np1);CHKERRQ(ierr);
ierr = DMDAGetElements(da, &ne, &nc, &necon);CHKERRQ(ierr);
for (e = 0; e < ne; e++) {
for (j = 0; j < 3; j++) {
idx[j] = necon[3*e+j];
values_m[j] = 0.0;
values_e[j] = 0.0;
}
/* Get basis function deriatives (we need the orientation of the element here) */
if (idx[1] > idx[0]) {
psi_x[0] = -hy; psi_x[1] = hy; psi_x[2] = 0.0;
psi_y[0] = -hx; psi_y[1] = 0.0; psi_y[2] = hx;
} else {
psi_x[0] = hy; psi_x[1] = -hy; psi_x[2] = 0.0;
psi_y[0] = hx; psi_y[1] = 0.0; psi_y[2] = -hx;
}
/* <\nabla\psi, F^*>: Divergence of the predicted convective fluxes */
for (j = 0; j < 3; j++) {
values_m[j] += (psi_x[j]*(phi*rho_u_np1[idx[j]] + rho_u_n[idx[j]]) + psi_y[j]*(rho_v_np1[idx[j]] + rho_v_n[idx[j]]))/3.0;
values_e[j] += values_m[j]*((rho_e_n[idx[j]] + p_n[idx[j]]) / rho_n[idx[j]]);
}
/* -<\nabla\psi, F^n_v>: Divergence of the viscous fluxes */
for (j = 0; j < 3; j++) {
/* \tau_{xx} = 2/3 \mu(T) (2 {\partial u\over\partial x} - {\partial v\over\partial y}) */
/* \tau_{xy} = \mu(T) ( {\partial u\over\partial y} + {\partial v\over\partial x}) */
/* \tau_{yy} = 2/3 \mu(T) (2 {\partial v\over\partial y} - {\partial u\over\partial x}) */
/* q_x = -\kappa(T) {\partial T\over\partial x} */
/* q_y = -\kappa(T) {\partial T\over\partial y} */
/* above code commeted out - causing ininitialized variables. */
q_x =0; q_y =0;
mu = 0.0;
kappa = 0.0;
tau_xx = 0.0;
tau_xy = 0.0;
tau_yy = 0.0;
for (k = 0; k < 3; k++) {
mu += mu_n[idx[k]];
kappa += kappa_n[idx[k]];
//.........这里部分代码省略.........
开发者ID:hsahasra，项目名称:petsc-magma-dense-mat，代码行数:101，代码来源:ex31.c
示例3: PCISSetUp
//.........这里部分代码省略.........
ierr = VecScatterCreate(pcis->vec1_N,pcis->is_B_local,pcis->vec1_B,(IS)0,&pcis->N_to_B);CHKERRQ(ierr);
ierr = VecScatterCreate(pcis->vec1_global,pcis->is_B_global,pcis->vec1_B,(IS)0,&pcis->global_to_B);CHKERRQ(ierr);
/* Creating scaling "matrix" D */
ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_is_use_stiffness_scaling",&pcis->use_stiffness_scaling,NULL);CHKERRQ(ierr);
if (!pcis->D) {
ierr = VecDuplicate(pcis->vec1_B,&pcis->D);CHKERRQ(ierr);
if (!pcis->use_stiffness_scaling) {
ierr = VecSet(pcis->D,pcis->scaling_factor);CHKERRQ(ierr);
} else {
ierr = MatGetDiagonal(matis->A,pcis->vec1_N);CHKERRQ(ierr);
ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_N,pcis->D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
ierr = VecScatterEnd (pcis->N_to_B,pcis->vec1_N,pcis->D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
}
}
ierr = VecCopy(pcis->D,pcis->vec1_B);CHKERRQ(ierr);
ierr = MatCreateVecs(pc->pmat,&counter,0);CHKERRQ(ierr); /* temporary auxiliar vector */
ierr = VecSet(counter,0.0);CHKERRQ(ierr);
ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,counter,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
ierr = VecScatterEnd (pcis->global_to_B,pcis->vec1_B,counter,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
ierr = VecScatterBegin(pcis->global_to_B,counter,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
ierr = VecScatterEnd (pcis->global_to_B,counter,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
ierr = VecPointwiseDivide(pcis->D,pcis->D,pcis->vec1_B);CHKERRQ(ierr);
ierr = VecDestroy(&counter);CHKERRQ(ierr);
/* See historical note 01, at the bottom of this file. */
/*
Creating the KSP contexts for the local Dirichlet and Neumann problems.
*/
if (pcis->computesolvers) {
PC pc_ctx;
/* Dirichlet */
ierr = KSPCreate(PETSC_COMM_SELF,&pcis->ksp_D);CHKERRQ(ierr);
ierr = PetscObjectIncrementTabLevel((PetscObject)pcis->ksp_D,(PetscObject)pc,1);CHKERRQ(ierr);
ierr = KSPSetOperators(pcis->ksp_D,pcis->A_II,pcis->A_II);CHKERRQ(ierr);
ierr = KSPSetOptionsPrefix(pcis->ksp_D,"is_localD_");CHKERRQ(ierr);
ierr = KSPGetPC(pcis->ksp_D,&pc_ctx);CHKERRQ(ierr);
ierr = PCSetType(pc_ctx,PCLU);CHKERRQ(ierr);
ierr = KSPSetType(pcis->ksp_D,KSPPREONLY);CHKERRQ(ierr);
ierr = KSPSetFromOptions(pcis->ksp_D);CHKERRQ(ierr);
/* the vectors in the following line are dummy arguments, just telling the KSP the vector size. Values are not used */
ierr = KSPSetUp(pcis->ksp_D);CHKERRQ(ierr);
/* Neumann */
ierr = KSPCreate(PETSC_COMM_SELF,&pcis->ksp_N);CHKERRQ(ierr);
ierr = PetscObjectIncrementTabLevel((PetscObject)pcis->ksp_N,(PetscObject)pc,1);CHKERRQ(ierr);
ierr = KSPSetOperators(pcis->ksp_N,matis->A,matis->A);CHKERRQ(ierr);
ierr = KSPSetOptionsPrefix(pcis->ksp_N,"is_localN_");CHKERRQ(ierr);
ierr = KSPGetPC(pcis->ksp_N,&pc_ctx);CHKERRQ(ierr);
ierr = PCSetType(pc_ctx,PCLU);CHKERRQ(ierr);
ierr = KSPSetType(pcis->ksp_N,KSPPREONLY);CHKERRQ(ierr);
ierr = KSPSetFromOptions(pcis->ksp_N);CHKERRQ(ierr);
{
PetscBool damp_fixed = PETSC_FALSE,
remove_nullspace_fixed = PETSC_FALSE,
set_damping_factor_floating = PETSC_FALSE,
not_damp_floating = PETSC_FALSE,
not_remove_nullspace_floating = PETSC_FALSE;
PetscReal fixed_factor,
floating_factor;
ierr = PetscOptionsGetReal(((PetscObject)pc_ctx)->prefix,"-pc_is_damp_fixed",&fixed_factor,&damp_fixed);CHKERRQ(ierr);
if (!damp_fixed) fixed_factor = 0.0;
ierr = PetscOptionsGetBool(((PetscObject)pc_ctx)->prefix,"-pc_is_damp_fixed",&damp_fixed,NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetBool(((PetscObject)pc_ctx)->prefix,"-pc_is_remove_nullspace_fixed",&remove_nullspace_fixed,NULL);CHKERRQ(ierr);
开发者ID:fengyuqi，项目名称:petsc，代码行数:67，代码来源:pcis.c
示例4: main
int main(int argc, char **argv)
{
PetscErrorCode ierr;
PetscInt n=10000,its,dfid=1;
Vec x,b,u;
Mat A;
KSP ksp;
PC pc,pcnoise;
PCNoise_Ctx ctx={0,NULL};
PetscReal eta=0.1,norm;
PetscScalar(*diagfunc)(PetscInt,PetscInt);
ierr = PetscInitialize(&argc,&argv,(char*)0,help);CHKERRQ(ierr);
/* Process command line options */
ierr = PetscOptionsGetInt(NULL,"-n",&n,NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetReal(NULL,"-eta",&eta,NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(NULL,"-diagfunc",&dfid,NULL);CHKERRQ(ierr);
switch(dfid){
case 1:
diagfunc = diagFunc1;
break;
case 2:
diagfunc = diagFunc2;
break;
case 3:
diagfunc = diagFunc3;
break;
default:
SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Unrecognized diagfunc option");
}
/* Create a diagonal matrix with a given distribution of diagonal elements */
ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr);
ierr = MatSetFromOptions(A);CHKERRQ(ierr);
ierr = MatSetUp(A);CHKERRQ(ierr);
ierr = AssembleDiagonalMatrix(A,diagfunc);CHKERRQ(ierr);
/* Allocate vectors and manufacture an exact solution and rhs */
ierr = MatCreateVecs(A,&x,NULL);CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject)x,"Computed Solution");CHKERRQ(ierr);
ierr = MatCreateVecs(A,&b,NULL);CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject)b,"RHS");CHKERRQ(ierr);
ierr = MatCreateVecs(A,&u,NULL);CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject)u,"Reference Solution");CHKERRQ(ierr);
ierr = VecSet(u,1.0);CHKERRQ(ierr);
ierr = MatMult(A,u,b);CHKERRQ(ierr);
/* Create a KSP object */
ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
/* Set up a composite preconditioner */
ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
ierr = PCSetType(pc,PCCOMPOSITE);CHKERRQ(ierr); /* default composite with single Identity PC */
ierr = PCCompositeSetType(pc,PC_COMPOSITE_ADDITIVE);CHKERRQ(ierr);
ierr = PCCompositeAddPC(pc,PCNONE);CHKERRQ(ierr);
if(eta > 0){
ierr = PCCompositeAddPC(pc,PCSHELL);CHKERRQ(ierr);
ierr = PCCompositeGetPC(pc,1,&pcnoise);CHKERRQ(ierr);
ctx.eta = eta;
ierr = PCShellSetContext(pcnoise,&ctx);CHKERRQ(ierr);
ierr = PCShellSetApply(pcnoise,PCApply_Noise);CHKERRQ(ierr);
ierr = PCShellSetSetUp(pcnoise,PCSetup_Noise);CHKERRQ(ierr);
ierr = PCShellSetDestroy(pcnoise,PCDestroy_Noise);CHKERRQ(ierr);
ierr = PCShellSetName(pcnoise,"Noise PC");CHKERRQ(ierr);
}
/* Set KSP from options (this can override the PC just defined) */
ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
/* Solve */
ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
/* Compute error */
ierr = VecAXPY(x,-1.0,u);CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject)x,"Error");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);
/* Destroy objects and finalize */
ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
ierr = MatDestroy(&A);CHKERRQ(ierr);
ierr = VecDestroy(&x);CHKERRQ(ierr);
ierr = VecDestroy(&b);CHKERRQ(ierr);
ierr = VecDestroy(&u);CHKERRQ(ierr);
PetscFinalize();
return 0;
}
开发者ID:pombredanne，项目名称:petsc，代码行数:92，代码来源:ex60.c
示例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 */
KSP *subksp; /* array of local KSP contexts on this processor */
PC pc; /* PC context */
PC subpc; /* PC context for subdomain */
PetscReal norm; /* norm of solution error */
PetscErrorCode ierr;
PetscInt i,j,Ii,J,*blks,m = 8,n;
PetscMPIInt rank,size;
PetscInt its,nlocal,first,Istart,Iend;
PetscScalar v,one = 1.0,none = -1.0;
PetscBool isbjacobi,flg = PETSC_FALSE;
PetscInitialize(&argc,&args,(char*)0,help);
ierr = PetscOptionsGetInt(NULL,"-m",&m,NULL);CHKERRQ(ierr);
ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
n = m+2;
/* -------------------------------------------------------------------
Compute the matrix and right-hand-side vector that define
the linear system, Ax = b.
------------------------------------------------------------------- */
/*
Create and assemble parallel matrix
*/
ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);CHKERRQ(ierr);
ierr = MatSetFromOptions(A);CHKERRQ(ierr);
ierr = MatMPIAIJSetPreallocation(A,5,NULL,5,NULL);CHKERRQ(ierr);
ierr = MatSeqAIJSetPreallocation(A,5,NULL);CHKERRQ(ierr);
ierr = MatGetOwnershipRange(A,&Istart,&Iend);CHKERRQ(ierr);
for (Ii=Istart; Ii<Iend; Ii++) {
v = -1.0; i = Ii/n; j = Ii - i*n;
if (i>0) {J = Ii - n; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);}
if (i<m-1) {J = Ii + n; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);}
if (j>0) {J = Ii - 1; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);}
if (j<n-1) {J = Ii + 1; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);}
v = 4.0; ierr = MatSetValues(A,1,&Ii,1,&Ii,&v,ADD_VALUES);CHKERRQ(ierr);
}
ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
/*
Create parallel vectors
*/
ierr = VecCreate(PETSC_COMM_WORLD,&u);CHKERRQ(ierr);
ierr = VecSetSizes(u,PETSC_DECIDE,m*n);CHKERRQ(ierr);
ierr = VecSetFromOptions(u);CHKERRQ(ierr);
ierr = VecDuplicate(u,&b);CHKERRQ(ierr);
ierr = VecDuplicate(b,&x);CHKERRQ(ierr);
/*
Set exact solution; then compute right-hand-side vector.
*/
ierr = VecSet(u,one);CHKERRQ(ierr);
ierr = MatMult(A,u,b);CHKERRQ(ierr);
/*
Create linear solver context
*/
ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
/*
Set operators. Here the matrix that defines the linear system
also serves as the preconditioning matrix.
*/
ierr = KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
/*
Set default preconditioner for this program to be block Jacobi.
This choice can be overridden at runtime with the option
-pc_type <type>
*/
ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
ierr = PCSetType(pc,PCBJACOBI);CHKERRQ(ierr);
/* -------------------------------------------------------------------
Define the problem decomposition
------------------------------------------------------------------- */
/*
Call PCBJacobiSetTotalBlocks() to set individually the size of
each block in the preconditioner. This could also be done with
the runtime option
-pc_bjacobi_blocks <blocks>
Also, see the command PCBJacobiSetLocalBlocks() to set the
local blocks.
Note: The default decomposition is 1 block per processor.
*/
ierr = PetscMalloc(m*sizeof(PetscInt),&blks);CHKERRQ(ierr);
for (i=0; i<m; i++) blks[i] = n;
ierr = PCBJacobiSetTotalBlocks(pc,m,blks);CHKERRQ(ierr);
ierr = PetscFree(blks);CHKERRQ(ierr);
//.........这里部分代码省略.........
开发者ID:feelpp，项目名称:debian-petsc，代码行数:101，代码来源:ex7.c
示例6: PCBDDCNullSpaceAssembleCorrection
//.........这里部分代码省略.........
ierr = VecAssemblyBegin(shell_ctx->work_small_1);CHKERRQ(ierr);
ierr = VecAssemblyEnd(shell_ctx->work_small_1);CHKERRQ(ierr);
ierr = MatSolve(small_mat,shell_ctx->work_small_1,shell_ctx->work_small_2);CHKERRQ(ierr);
ierr = VecGetArrayRead(shell_ctx->work_small_2,(const PetscScalar**)&array);CHKERRQ(ierr);
for (i=0;i<basis_size;i++) {
array_mat[i*basis_size+k]=array[i];
}
ierr = VecRestoreArrayRead(shell_ctx->work_small_2,(const PetscScalar**)&array);CHKERRQ(ierr);
}
ierr = MatCreateSeqDense(PETSC_COMM_SELF,basis_size,basis_size,array_mat,&inv_small_mat);CHKERRQ(ierr);
ierr = MatMatMult(shell_ctx->basis_mat,inv_small_mat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&shell_ctx->Lbasis_mat);CHKERRQ(ierr);
ierr = PetscFree(array_mat);CHKERRQ(ierr);
ierr = MatDestroy(&inv_small_mat);CHKERRQ(ierr);
ierr = MatDestroy(&small_mat);CHKERRQ(ierr);
ierr = MatScale(shell_ctx->Kbasis_mat,m_one);CHKERRQ(ierr);
/* Rebuild local PC */
ierr = KSPGetPC(local_ksp,&shell_ctx->local_pc);CHKERRQ(ierr);
ierr = PetscObjectReference((PetscObject)shell_ctx->local_pc);CHKERRQ(ierr);
ierr = PCCreate(PETSC_COMM_SELF,&newpc);CHKERRQ(ierr);
ierr = PCSetOperators(newpc,local_mat,local_mat);CHKERRQ(ierr);
ierr = PCSetType(newpc,PCSHELL);CHKERRQ(ierr);
ierr = PCShellSetContext(newpc,shell_ctx);CHKERRQ(ierr);
ierr = PCShellSetApply(newpc,PCBDDCApplyNullSpaceCorrectionPC);CHKERRQ(ierr);
ierr = PCShellSetDestroy(newpc,PCBDDCDestroyNullSpaceCorrectionPC);CHKERRQ(ierr);
ierr = PCSetUp(newpc);CHKERRQ(ierr);
ierr = KSPSetPC(local_ksp,newpc);CHKERRQ(ierr);
ierr = PCDestroy(&newpc);CHKERRQ(ierr);
ierr = KSPSetUp(local_ksp);CHKERRQ(ierr);
}
/* test */
if (pcbddc->dbg_flag && basis_dofs) {
KSP check_ksp;
PC check_pc;
Mat test_mat;
Vec work3;
PetscReal test_err,lambda_min,lambda_max;
PetscBool setsym,issym=PETSC_FALSE;
PetscInt tabs;
ierr = PetscViewerASCIIGetTab(pcbddc->dbg_viewer,&tabs);CHKERRQ(ierr);
ierr = KSPGetPC(local_ksp,&check_pc);CHKERRQ(ierr);
ierr = VecDuplicate(shell_ctx->work_full_1,&work1);CHKERRQ(ierr);
ierr = VecDuplicate(shell_ctx->work_full_1,&work2);CHKERRQ(ierr);
ierr = VecDuplicate(shell_ctx->work_full_1,&work3);CHKERRQ(ierr);
ierr = VecSetRandom(shell_ctx->work_small_1,NULL);CHKERRQ(ierr);
ierr = MatMult(shell_ctx->basis_mat,shell_ctx->work_small_1,work1);CHKERRQ(ierr);
ierr = VecCopy(work1,work2);CHKERRQ(ierr);
ierr = MatMult(local_mat,work1,work3);CHKERRQ(ierr);
ierr = PCApply(check_pc,work3,work1);CHKERRQ(ierr);
ierr = VecAXPY(work1,m_one,work2);CHKERRQ(ierr);
ierr = VecNorm(work1,NORM_INFINITY,&test_err);CHKERRQ(ierr);
ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d error for nullspace correction for ",PetscGlobalRank);CHKERRQ(ierr);
ierr = PetscViewerASCIIUseTabs(pcbddc->dbg_viewer,PETSC_FALSE);CHKERRQ(ierr);
if (isdir) {
ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Dirichlet ");CHKERRQ(ierr);
} else {
ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Neumann ");CHKERRQ(ierr);
}
ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"solver is :%1.14e\n",test_err);CHKERRQ(ierr);
ierr = PetscViewerASCIISetTab(pcbddc->dbg_viewer,tabs);CHKERRQ(ierr);
ierr = PetscViewerASCIIUseTabs(pcbddc->dbg_viewer,PETSC_TRUE);CHKERRQ(ierr);
ierr = MatTransposeMatMult(shell_ctx->Lbasis_mat,shell_ctx->Kbasis_mat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&test_mat);CHKERRQ(ierr);
ierr = MatShift(test_mat,one);CHKERRQ(ierr);
ierr = MatNorm(test_mat,NORM_INFINITY,&test_err);CHKERRQ(ierr);
ierr = MatDestroy(&test_mat);CHKERRQ(ierr);
ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d error for nullspace matrices is :%1.14e\n",PetscGlobalRank,test_err);CHKERRQ(ierr);
/* Create ksp object suitable for extreme eigenvalues' estimation */
ierr = KSPCreate(PETSC_COMM_SELF,&check_ksp);CHKERRQ(ierr);
ierr = KSPSetErrorIfNotConverged(check_ksp,pc->erroriffailure);CHKERRQ(ierr);
ierr = KSPSetOperators(check_ksp,local_mat,local_mat);CHKERRQ(ierr);
ierr = KSPSetTolerances(check_ksp,1.e-8,1.e-8,PETSC_DEFAULT,basis_dofs);CHKERRQ(ierr);
ierr = KSPSetComputeSingularValues(check_ksp,PETSC_TRUE);CHKERRQ(ierr);
ierr = MatIsSymmetricKnown(pc->pmat,&setsym,&issym);CHKERRQ(ierr);
if (issym) {
ierr = KSPSetType(check_ksp,KSPCG);CHKERRQ(ierr);
}
ierr = KSPSetPC(check_ksp,check_pc);CHKERRQ(ierr);
ierr = KSPSetUp(check_ksp);CHKERRQ(ierr);
ierr = VecSetRandom(work1,NULL);CHKERRQ(ierr);
ierr = MatMult(local_mat,work1,work2);CHKERRQ(ierr);
ierr = KSPSolve(check_ksp,work2,work2);CHKERRQ(ierr);
ierr = VecAXPY(work2,m_one,work1);CHKERRQ(ierr);
ierr = VecNorm(work2,NORM_INFINITY,&test_err);CHKERRQ(ierr);
ierr = KSPComputeExtremeSingularValues(check_ksp,&lambda_max,&lambda_min);CHKERRQ(ierr);
ierr = KSPGetIterationNumber(check_ksp,&k);CHKERRQ(ierr);
ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d error for adapted KSP %1.14e (it %d, eigs %1.6e %1.6e)\n",PetscGlobalRank,test_err,k,lambda_min,lambda_max);CHKERRQ(ierr);
ierr = KSPDestroy(&check_ksp);CHKERRQ(ierr);
ierr = VecDestroy(&work1);CHKERRQ(ierr);
ierr = VecDestroy(&work2);CHKERRQ(ierr);
ierr = VecDestroy(&work3);CHKERRQ(ierr);
}
/* all processes shoud call this, even the void ones */
if (pcbddc->dbg_flag) {
ierr = PetscViewerFlush(pcbddc->dbg_viewer);CHKERRQ(ierr);
}
PetscFunctionReturn(0);
}
开发者ID:plguhur，项目名称:petsc，代码行数:101，代码来源:bddcnullspace.c
示例7: main
//.........这里部分代码省略.........
ierr = MatSetValuesStencil(H, 1, &sxy, 1, &sxy_m, &valconj, ADD_VALUES);CHKERRQ(ierr);
sxy.c = 1; sxy_m.c = 0; /* spin 1, 0 */
val = -PETSC_i*uxy2; valconj = PetscConj(val);
ierr = MatSetValuesStencil(H, 1, &sxy_m, 1, &sxy, &val, ADD_VALUES);CHKERRQ(ierr);
ierr = MatSetValuesStencil(H, 1, &sxy, 1, &sxy_m, &valconj, ADD_VALUES);CHKERRQ(ierr);
sxy.c = 1; sxy_m.c = 1; /* spin 1, 1 */
val = PetscConj(uxy2); valconj = PetscConj(val);
ierr = MatSetValuesStencil(H, 1, &sxy_m, 1, &sxy, &val, ADD_VALUES);CHKERRQ(ierr);
ierr = MatSetValuesStencil(H, 1, &sxy, 1, &sxy_m, &valconj, ADD_VALUES);CHKERRQ(ierr);
}
}
ierr = MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
/* scale H */
ierr = MatScale(H, 1./(2.*h));CHKERRQ(ierr);
/* it looks like H is Hermetian */
/* construct normal equations */
ierr = MatMatMult(H, H, MAT_INITIAL_MATRIX, 1., &HtH);CHKERRQ(ierr);
/* permutation matrix to check whether H and HtH are identical to the ones in the paper */
/* Mat perm; */
/* ierr = DMCreateMatrix(da, &perm);CHKERRQ(ierr); */
/* PetscInt row, col; */
/* PetscScalar one = 1.0; */
/* for (PetscInt i=0; i<n; i++) { */
/* for (PetscInt j=0; j<n; j++) { */
/* row = (i*n+j)*2; col = i*n+j; */
/* ierr = MatSetValues(perm, 1, &row, 1, &col, &one, INSERT_VALUES);CHKERRQ(ierr); */
/* row = (i*n+j)*2+1; col = i*n+j + n*n; */
/* ierr = MatSetValues(perm, 1, &row, 1, &col, &one, INSERT_VALUES);CHKERRQ(ierr); */
/* } */
/* } */
/* ierr = MatAssemblyBegin(perm, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); */
/* ierr = MatAssemblyEnd(perm, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); */
/* Mat Hperm; */
/* ierr = MatPtAP(H, perm, MAT_INITIAL_MATRIX, 1.0, &Hperm);CHKERRQ(ierr); */
/* ierr = PetscPrintf(PETSC_COMM_WORLD, "Matrix H after construction\n");CHKERRQ(ierr); */
/* ierr = MatView(Hperm, PETSC_VIEWER_STDOUT_(PETSC_COMM_WORLD));CHKERRQ(ierr); */
/* Mat HtHperm; */
/* ierr = MatPtAP(HtH, perm, MAT_INITIAL_MATRIX, 1.0, &HtHperm);CHKERRQ(ierr); */
/* ierr = PetscPrintf(PETSC_COMM_WORLD, "Matrix HtH:\n");CHKERRQ(ierr); */
/* ierr = MatView(HtHperm, PETSC_VIEWER_STDOUT_(PETSC_COMM_WORLD));CHKERRQ(ierr); */
/* right hand side */
ierr = DMCreateGlobalVector(da, &b);CHKERRQ(ierr);
ierr = VecSet(b,0.0);CHKERRQ(ierr);
ierr = VecSetValues(b, 1, ix, vals, INSERT_VALUES);CHKERRQ(ierr);
ierr = VecAssemblyBegin(b);CHKERRQ(ierr);
ierr = VecAssemblyEnd(b);CHKERRQ(ierr);
/* ierr = VecSetRandom(b, rctx);CHKERRQ(ierr); */
ierr = VecDuplicate(b, &Htb);CHKERRQ(ierr);
ierr = MatMultTranspose(H, b, Htb);CHKERRQ(ierr);
/* construct solver */
ierr = KSPCreate(PETSC_COMM_WORLD,&kspmg);CHKERRQ(ierr);
ierr = KSPSetType(kspmg, KSPCG);CHKERRQ(ierr);
ierr = KSPGetPC(kspmg,&pcmg);CHKERRQ(ierr);
ierr = PCSetType(pcmg,PCASA);CHKERRQ(ierr);
/* maybe user wants to override some of the choices */
ierr = KSPSetFromOptions(kspmg);CHKERRQ(ierr);
ierr = KSPSetOperators(kspmg, HtH, HtH, DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
ierr = DMDASetRefinementFactor(da, 3, 3, 3);CHKERRQ(ierr);
ierr = PCSetDM(pcmg,da);CHKERRQ(ierr);
ierr = PCASASetTolerances(pcmg, 1.e-6, 1.e-10,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr);
ierr = VecDuplicate(b, &xvec);CHKERRQ(ierr);
ierr = VecSet(xvec, 0.0);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Solve the linear system
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = KSPSolve(kspmg, Htb, xvec);CHKERRQ(ierr);
/* ierr = VecView(xvec, PETSC_VIEWER_STDOUT_(PETSC_COMM_WORLD));CHKERRQ(ierr); */
ierr = KSPDestroy(&kspmg);CHKERRQ(ierr);
ierr = VecDestroy(&xvec);CHKERRQ(ierr);
/* seems to be destroyed by KSPDestroy */
ierr = VecDestroy(&b);CHKERRQ(ierr);
ierr = VecDestroy(&Htb);CHKERRQ(ierr);
ierr = MatDestroy(&HtH);CHKERRQ(ierr);
ierr = MatDestroy(&H);CHKERRQ(ierr);
ierr = DMDestroy(&da);CHKERRQ(ierr);
ierr = PetscRandomDestroy(&rctx);CHKERRQ(ierr);
ierr = PetscFinalize();
return 0;
}
开发者ID:hsahasra，项目名称:petsc-magma-dense-mat，代码行数:101，代码来源:ex39.c
示例8: main
int main(int argc,char **argv)
{
PetscErrorCode ierr;
PetscInt its,n,Nx=PETSC_DECIDE,Ny=PETSC_DECIDE,nlocal;
PetscMPIInt size;
PetscScalar one = 1.0;
PetscInt mx,my;
Mat A;
GridCtx fine_ctx;
KSP ksp;
PetscBool flg;
ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
/* set up discretization matrix for fine grid */
fine_ctx.mx = 9; fine_ctx.my = 9;
ierr = PetscOptionsGetInt(NULL,NULL,"-mx",&mx,&flg);CHKERRQ(ierr);
if (flg) fine_ctx.mx = mx;
ierr = PetscOptionsGetInt(NULL,NULL,"-my",&my,&flg);CHKERRQ(ierr);
if (flg) fine_ctx.my = my;
ierr = PetscPrintf(PETSC_COMM_WORLD,"Fine grid size %D by %D\n",fine_ctx.mx,fine_ctx.my);CHKERRQ(ierr);
n = fine_ctx.mx*fine_ctx.my;
MPI_Comm_size(PETSC_COMM_WORLD,&size);
ierr = PetscOptionsGetInt(NULL,NULL,"-Nx",&Nx,NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(NULL,NULL,"-Ny",&Ny,NULL);CHKERRQ(ierr);
/* Set up distributed array for fine grid */
ierr = DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,fine_ctx.mx,
fine_ctx.my,Nx,Ny,1,1,NULL,NULL,&fine_ctx.da);CHKERRQ(ierr);
ierr = DMCreateGlobalVector(fine_ctx.da,&fine_ctx.x);CHKERRQ(ierr);
ierr = VecDuplicate(fine_ctx.x,&fine_ctx.b);CHKERRQ(ierr);
ierr = VecGetLocalSize(fine_ctx.x,&nlocal);CHKERRQ(ierr);
ierr = DMCreateLocalVector(fine_ctx.da,&fine_ctx.localX);CHKERRQ(ierr);
ierr = VecDuplicate(fine_ctx.localX,&fine_ctx.localF);CHKERRQ(ierr);
ierr = MatCreateAIJ(PETSC_COMM_WORLD,nlocal,nlocal,n,n,5,NULL,3,NULL,&A);CHKERRQ(ierr);
ierr = FormJacobian_Grid(&fine_ctx,&A);CHKERRQ(ierr);
/* create linear solver */
ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
/* set values for rhs vector */
ierr = VecSet(fine_ctx.b,one);CHKERRQ(ierr);
/* set options, then solve system */
ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr); /* calls PCSetFromOptions_ML if 'pc_type=ml' */
ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
ierr = KSPSolve(ksp,fine_ctx.b,fine_ctx.x);CHKERRQ(ierr);
ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %D\n",its);CHKERRQ(ierr);
/* free data structures */
ierr = VecDestroy(&fine_ctx.x);CHKERRQ(ierr);
ierr = VecDestroy(&fine_ctx.b);CHKERRQ(ierr);
ierr = DMDestroy(&fine_ctx.da);CHKERRQ(ierr);
ierr = VecDestroy(&fine_ctx.localX);CHKERRQ(ierr);
ierr = VecDestroy(&fine_ctx.localF);CHKERRQ(ierr);
ierr = MatDestroy(&A);CHKERRQ(ierr);
ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
ierr = PetscFinalize();
return ierr;
}
开发者ID:tom-klotz，项目名称:petsc，代码行数:62，代码来源:ex26.c
示例9: TaoCreate_TRON
PETSC_EXTERN PetscErrorCode TaoCreate_TRON(Tao tao)
{
TAO_TRON *tron;
PetscErrorCode ierr;
const char *morethuente_type = TAOLINESEARCHMT;
PetscFunctionBegin;
tao->ops->setup = TaoSetup_TRON;
tao->ops->solve = TaoSolve_TRON;
tao->ops->view = TaoView_TRON;
tao->ops->setfromoptions = TaoSetFromOptions_TRON;
tao->ops->destroy = TaoDestroy_TRON;
tao->ops->computedual = TaoComputeDual_TRON;
ierr = PetscNewLog(tao,&tron);CHKERRQ(ierr);
tao->data = (void*)tron;
/* Override default settings (unless already changed) */
if (!tao->max_it_changed) tao->max_it = 50;
#if defined(PETSC_USE_REAL_SINGLE)
if (!tao->steptol_changed) tao->steptol = 1.0e-6;
#else
if (!tao->steptol_changed) tao->steptol = 1.0e-12;
#endif
if (!tao->trust0_changed) tao->trust0 = 1.0;
/* Initialize pointers and variables */
tron->n = 0;
tron->maxgpits = 3;
tron->pg_ftol = 0.001;
tron->eta1 = 1.0e-4;
tron->eta2 = 0.25;
tron->eta3 = 0.50;
tron->eta4 = 0.90;
tron->sigma1 = 0.5;
tron->sigma2 = 2.0;
tron->sigma3 = 4.0;
tron->gp_iterates = 0; /* Cumulative number */
tron->total_gp_its = 0;
tron->n_free = 0;
tron->DXFree=NULL;
tron->R=NULL;
tron->X_New=NULL;
tron->G_New=NULL;
tron->Work=NULL;
tron->Free_Local=NULL;
tron->H_sub=NULL;
tron->Hpre_sub=NULL;
tao->subset_type = TAO_SUBSET_SUBVEC;
ierr = TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);CHKERRQ(ierr);
ierr = TaoLineSearchSetType(tao->linesearch,morethuente_type);CHKERRQ(ierr);
ierr = TaoLineSearchUseTaoRoutines(tao->linesearch,tao);CHKERRQ(ierr);
ierr = TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);CHKERRQ(ierr);
ierr = KSPCreate(((PetscObject)tao)->comm, &tao->ksp);CHKERRQ(ierr);
ierr = KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);CHKERRQ(ierr);
ierr = KSPSetType(tao->ksp,KSPSTCG);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
开发者ID:plguhur，项目名称:petsc，代码行数:66，代码来源:tron.c
示例10: main
int main(int argc,char **args)
{
KSP subksp;
Mat A,subA;
Vec x,b,u,subb,subx,subu;
PetscViewer fd;
char file[PETSC_MAX_PATH_LEN];
PetscBool flg;
PetscErrorCode ierr;
PetscInt i,m,n,its;
PetscReal norm;
PetscMPIInt rank,size;
MPI_Comm comm,subcomm;
PetscSubcomm psubcomm;
PetscInt nsubcomm=1,id;
PetscScalar *barray,*xarray,*uarray,*array,one=1.0;
PetscInt type=1;
PetscInitialize(&argc,&args,(char*)0,help);
/* Load the matrix */
ierr = PetscOptionsGetString(NULL,"-f",file,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr);
if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate binary file with the -f option");
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);CHKERRQ(ierr);
ierr = MatLoad(A,fd);CHKERRQ(ierr);
ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr);
ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr);
ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
/* Create rhs vector b */
ierr = MatGetLocalSize(A,&m,NULL);CHKERRQ(ierr);
ierr = VecCreate(PETSC_COMM_WORLD,&b);CHKERRQ(ierr);
ierr = VecSetSizes(b,m,PETSC_DECIDE);CHKERRQ(ierr);
ierr = VecSetFromOptions(b);CHKERRQ(ierr);
ierr = VecSet(b,one);CHKERRQ(ierr);
ierr = VecDuplicate(b,&x);CHKERRQ(ierr);
ierr = VecDuplicate(b,&u);CHKERRQ(ierr);
ierr = VecSet(x,0.0);CHKERRQ(ierr);
/* Test MatGetMultiProcBlock() */
ierr = PetscOptionsGetInt(NULL,"-nsubcomm",&nsubcomm,NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(NULL,"-subcomm_type",&type,NULL);CHKERRQ(ierr);
ierr = PetscSubcommCreate(comm,&psubcomm);CHKERRQ(ierr);
ierr = PetscSubcommSetNumber(psubcomm,nsubcomm);CHKERRQ(ierr);
if (type == PETSC_SUBCOMM_GENERAL) { /* user provides color, subrank and duprank */
PetscMPIInt color,subrank,duprank,subsize;
duprank = size-1 - rank;
subsize = size/nsubcomm;
if (subsize*nsubcomm != size) SETERRQ2(comm,PETSC_ERR_SUP,"This example requires nsubcomm %D divides nproc %D",nsubcomm,size);
color = duprank/subsize;
subrank = duprank - color*subsize;
ierr = PetscSubcommSetTypeGeneral(psubcomm,color,subrank,duprank);CHKERRQ(ierr);
} else if (type == PETSC_SUBCOMM_CONTIGUOUS) {
ierr = PetscSubcommSetType(psubcomm,PETSC_SUBCOMM_CONTIGUOUS);CHKERRQ(ierr);
} else if (type == PETSC_SUBCOMM_INTERLACED) {
ierr = PetscSubcommSetType(psubcomm,PETSC_SUBCOMM_INTERLACED);CHKERRQ(ierr);
} else SETERRQ1(psubcomm->parent,PETSC_ERR_SUP,"PetscSubcommType %D is not supported yet",type);
subcomm = psubcomm->comm;
ierr = PetscOptionsHasName(NULL, "-subcomm_view", &flg);CHKERRQ(ierr);
if (flg) {
PetscMPIInt subsize,subrank,duprank;
ierr = MPI_Comm_size((MPI_Comm)subcomm,&subsize);CHKERRQ(ierr);
ierr = MPI_Comm_rank((MPI_Comm)subcomm,&subrank);CHKERRQ(ierr);
ierr = MPI_Comm_rank((MPI_Comm)psubcomm->dupparent,&duprank);CHKERRQ(ierr);
ierr = PetscSynchronizedPrintf(comm,"[%D], color %D, sub-size %D, sub-rank %D, duprank %D\n",rank,psubcomm->color,subsize,subrank,duprank);
ierr = PetscSynchronizedFlush(comm);CHKERRQ(ierr);
}
/* Create subA */
ierr = MatGetMultiProcBlock(A,subcomm,MAT_INITIAL_MATRIX,&subA);CHKERRQ(ierr);
/* Create sub vectors without arrays. Place b's and x's local arrays into subb and subx */
ierr = MatGetLocalSize(subA,&m,&n);CHKERRQ(ierr);
ierr = VecCreateMPIWithArray(subcomm,1,m,PETSC_DECIDE,NULL,&subb);CHKERRQ(ierr);
ierr = VecCreateMPIWithArray(subcomm,1,n,PETSC_DECIDE,NULL,&subx);CHKERRQ(ierr);
ierr = VecCreateMPIWithArray(subcomm,1,n,PETSC_DECIDE,NULL,&subu);CHKERRQ(ierr);
ierr = VecGetArray(b,&barray);CHKERRQ(ierr);
ierr = VecGetArray(x,&xarray);CHKERRQ(ierr);
ierr = VecGetArray(u,&uarray);CHKERRQ(ierr);
ierr = VecPlaceArray(subb,barray);CHKERRQ(ierr);
ierr = VecPlaceArray(subx,xarray);CHKERRQ(ierr);
ierr = VecPlaceArray(subu,uarray);CHKERRQ(ierr);
/* Create linear solvers associated with subA */
ierr = KSPCreate(subcomm,&subksp);CHKERRQ(ierr);
ierr = KSPSetOperators(subksp,subA,subA,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
ierr = KSPSetFromOptions(subksp);CHKERRQ(ierr);
/* Solve sub systems */
ierr = KSPSolve(subksp,subb,subx);CHKERRQ(ierr);
ierr = KSPGetIterationNumber(subksp,&its);CHKERRQ(ierr);
//.........这里部分代码省略.........
开发者ID:feelpp，项目名称:debian-petsc，代码行数:101，代码来源:ex37.c
示例11: BSSCR_DRIVER_flex
//.........这里部分代码省略.........
useNormInfStoppingConditions = PETSC_FALSE;
PetscOptionsGetTruth( PETSC_NULL ,"-A11_use_norm_inf_stopping_condition", &useNormInfStoppingConditions, &found );
if(useNormInfStoppingConditions)
BSSCR_KSPSetNormInfConvergenceTest( ksp_inner );
useNormInfMonitor = PETSC_FALSE;
PetscOptionsGetTruth( PETSC_NULL, "-A11_ksp_norm_inf_monitor", &useNormInfMonitor, &found );
if(useNormInfMonitor)
KSPMonitorSet( ksp_inner, BSSCR_KSPNormInfMonitor, PETSC_NULL, PETSC_NULL );
useNormInfMonitor = PETSC_FALSE;
PetscOptionsGetTruth( PETSC_NULL, "-A11_ksp_norm_inf_to_norm_2_monitor", &useNormInfMonitor, &found );
if(useNormInfMonitor)
KSPMonitorSet( ksp_inner, BSSCR_KSPNormInfToNorm2Monitor, PETSC_NULL, PETSC_NULL );
/* create right hand side */
/* h_hat = G'*inv(K)*f - h */
MatGetVecs(K,PETSC_NULL,&t);
MatGetVecs( S, PETSC_NULL, &h_hat );
KSPSetOptionsPrefix( ksp_inner, "A11_" );
KSPSetFromOptions( ksp_inner );
KSPSolve(ksp_inner,f,t);/* t=f/K */
//bsscr_writeVec( t, "ts", "Writing t vector");
MatMult(D,t,h_hat);/* G'*t */
VecAXPY(h_hat, -1.0, h);/* h_hat = h_hat - h */
Stg_VecDestroy(&t);
//bsscr_writeVec( h_hat, "h_hat", "Writing h_hat Vector in Solver");
//MatSchurApplyReductionToVecFromBlock( S, stokes_b, h_hat );
/* create solver for S p = h_hat */
KSPCreate( PETSC_COMM_WORLD, &ksp_S );
KSPSetOptionsPrefix( ksp_S, "scr_");
Stg_KSPSetOperators( ksp_S, S,S, SAME_NONZERO_PATTERN );
KSPSetType( ksp_S, "cg" );
/* Build preconditioner for S */
KSPGetPC( ksp_S, &pc_S );
BSSCR_BSSCR_StokesCreatePCSchur2( K,G,D,C,approxS,pc_S,sym, bsscrp_self );
KSPSetFromOptions(ksp_S);
/* Set specific monitor test */
KSPGetTolerances( ksp_S, PETSC_NULL, PETSC_NULL, PETSC_NULL, &max_it );
//BSSCR_KSPLogSetMonitor( ksp_S, max_it, &monitor_index );
/* Pressure / Velocity Solve */
scrSolveTime = MPI_Wtime();
PetscPrintf( PETSC_COMM_WORLD, "\t* Pressure / Velocity Solve \n");
usePreviousGuess = PETSC_FALSE;
if(been_here)
PetscOptionsGetTruth( PETSC_NULL, "-scr_use_previous_guess", &usePreviousGuess, &found );
if(usePreviousGuess) { /* Note this should actually look at checkpoint information */
KSPSetInitialGuessNonzero( ksp_S, PETSC_TRUE );
}
else {
KSPSetInitialGuessNonzero( ksp_S, PETSC_FALSE );
}
//KSPSetRelativeRhsConvergenceTest( ksp_S );
useNormInfStoppingConditions = PETSC_FALSE;
开发者ID:OlympusMonds，项目名称:EarthByte_Underworld，代码行数:67，代码来源:bsscr-driver.c
示例12: main
int main(int argc, char** argv)
{
DM da;
PetscErrorCode ierr;
Vec x, rhs;
Mat A, jac;
ierr = PetscInitialize(&argc, &argv, NULL, NULL);
CHKERRQ(ierr);
ierr = PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Laplacian in 2D", "");
CHKERRQ(ierr);
ierr = PetscOptionsEnd();
CHKERRQ(ierr);
ierr = HpddmRegisterKSP();
CHKERRQ(ierr);
MPI_Barrier(PETSC_COMM_WORLD);
double time = MPI_Wtime();
ierr = DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 10, 10, PETSC_DECIDE, PETSC_DECIDE, 1, 1,
0, 0, &da);
CHKERRQ(ierr);
ierr = DMSetFromOptions(da);
CHKERRQ(ierr);
ierr = DMSetUp(da);
CHKERRQ(ierr);
ierr = DMCreateGlobalVector(da, &rhs);
CHKERRQ(ierr);
ierr = DMCreateGlobalVector(da, &x);
CHKERRQ(ierr);
ierr = DMCreateMatrix(da, &A);
CHKERRQ(ierr);
ierr = DMCreateMatrix(da, &jac);
CHKERRQ(ierr);
ierr = ComputeMatrix(da, jac, A);
CHKERRQ(ierr);
MPI_Barrier(PETSC_COMM_WORLD);
time = MPI_Wtime() - time;
ierr = PetscPrintf(PETSC_COMM_WORLD, "--- Mat assembly = %f\n", time);
CHKERRQ(ierr);
MPI_Barrier(PETSC_COMM_WORLD);
time = MPI_Wtime();
KSP ksp;
ierr = KSPCreate(PETSC_COMM_WORLD, &ksp);
CHKERRQ(ierr);
ierr = KSPSetDM(ksp, da);
CHKERRQ(ierr);
ierr = KSPSetFromOptions(ksp);
CHKERRQ(ierr);
ierr = KSPSetOperators(ksp, A, A);
CHKERRQ(ierr);
ierr = KSPSetDMActive(ksp, PETSC_FALSE);
CHKERRQ(ierr);
ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
CHKERRQ(ierr);
ierr = KSPSetUp(ksp);
CHKERRQ(ierr);
MPI_Barrier(PETSC_COMM_WORLD);
time = MPI_Wtime() - time;
ierr = PetscPrintf(PETSC_COMM_WORLD, "--- PC setup = %f\n", time);
CHKERRQ(ierr);
PetscScalar nus[SIZE_ARRAY_NU] = {0.1, 10.0, 0.001, 100.0};
float t_time[SIZE_ARRAY_NU];
int t_its[SIZE_ARRAY_NU];
int i, j;
for (j = 0; j < 2; ++j) {
{
if (j == 1) {
ierr = KSPSetType(ksp, "hpddm");
CHKERRQ(ierr);
ierr = KSPSetFromOptions(ksp);
CHKERRQ(ierr);
ierr = VecZeroEntries(x);
CHKERRQ(ierr);
}
ierr = KSPSolve(ksp, rhs, x);
CHKERRQ(ierr);
if (j == 1) {
const HpddmOption* const opt = HpddmOptionGet();
int previous = HpddmOptionVal(opt, "krylov_method");
if (previous == HPDDM_KRYLOV_METHOD_GCRODR || previous == HPDDM_KRYLOV_METHOD_BGCRODR) HpddmDestroyRecycling();
}
}
for (i = 0; i < SIZE_ARRAY_NU; ++i) {
ierr = VecZeroEntries(x);
CHKERRQ(ierr);
ierr = ComputeRHS(da, rhs, nus[i]);
CHKERRQ(ierr);
MPI_Barrier(PETSC_COMM_WORLD);
time = MPI_Wtime();
ierr = KSPSolve(ksp, rhs, x);
CHKERRQ(ierr);
MPI_Barrier(PETSC_COMM_WORLD);
t_time[i] = MPI_Wtime() - time;
PetscInt its;
ierr = KSPGetIterationNumber(ksp, &its);
CHKERRQ(ierr);
t_its[i] = its;
ierr = ComputeError(A, rhs, x);
CHKERRQ(ierr);
}
for (i = 0; i < SIZE_ARRAY_NU; ++i) {
ierr = PetscPrintf(PETSC_COMM_WORLD, "%d\t%d\t%f\n", i + 1, t_its[i], t_time[i]);
//.........这里部分代码省略.........
开发者ID:hpddm，项目名称:hpddm，代码行数:101，代码来源:ex32.c
示例13: main
int main(int argc,char **args)
{
Mat C;
PetscMPIInt rank,size;
PetscInt i,m = 5,N,start,end,M,its;
PetscScalar val,Ke[16],r[4];
PetscReal x,y,h,norm,tol=1.e-14;
PetscErrorCode ierr;
PetscInt idx[4],count,*rows;
Vec
|
客服电话
APP下载
官方微信
请发表评论