本文整理汇总了C++中MatGetVecs函数的典型用法代码示例。如果您正苦于以下问题:C++ MatGetVecs函数的具体用法?C++ MatGetVecs怎么用?C++ MatGetVecs使用的例子?那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了MatGetVecs函数的20个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于我们的系统推荐出更棒的C++代码示例。
示例1: bsscr_DGMiGtD
PetscErrorCode bsscr_DGMiGtD( Mat *_K2, Mat K, Mat G, Mat M){
Mat K2;
Vec diag;
Vec Mdiag;
Mat MinvGt;
Mat Gtrans;
PetscErrorCode ierr;
PetscFunctionBegin;
MatGetVecs( K, &diag, PETSC_NULL );
MatGetDiagonal( K, diag );
VecSqrt(diag);
//VecReciprocal(diag);/* trying something different here */
MatGetVecs( M, &Mdiag, PETSC_NULL );
MatGetDiagonal( M, Mdiag );
VecReciprocal(Mdiag);
#if( PETSC_VERSION_MAJOR <= 2 )
ierr=MatTranspose(G, &Gtrans);CHKERRQ(ierr);
#else
ierr=MatTranspose(G, MAT_INITIAL_MATRIX,&Gtrans);CHKERRQ(ierr);
#endif
ierr=MatConvert(Gtrans, MATSAME, MAT_INITIAL_MATRIX, &MinvGt);CHKERRQ(ierr);/* copy Gtrans -> MinvGt */
MatDiagonalScale(MinvGt, Mdiag, PETSC_NULL);/* Minv*Gtrans */
/* MAT_INITIAL_MATRIX -> creates K2 matrix : PETSC_DEFAULT for fill ratio: run with -info to find what it should be*/
ierr=MatMatMult( G, MinvGt, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &K2);CHKERRQ(ierr);/* K2 = G*Minv*Gtrans */
MatDiagonalScale(K2, diag, diag );/* K2 = D*K2*D = D*G*Minv*Gtrans*D */
Stg_MatDestroy(&Gtrans);
Stg_MatDestroy(&MinvGt);
Stg_VecDestroy(&diag);
Stg_VecDestroy(&Mdiag);
*_K2=K2;
PetscFunctionReturn(0);
}
开发者ID:dansand,项目名称:underworld2,代码行数:33,代码来源:createK2.c
示例2: KSPBuildPressure_CB_Nullspace_BSSCR
PetscErrorCode KSPBuildPressure_CB_Nullspace_BSSCR(KSP ksp)
{
KSP_BSSCR *bsscr = (KSP_BSSCR *)ksp->data;
FeEquationNumber *eq_num = bsscr->solver->st_sle->pSolnVec->feVariable->eqNum;
FeMesh *feMesh = bsscr->solver->st_sle->pSolnVec->feVariable->feMesh; /* is the pressure mesh */
unsigned ijk[3];
Vec t, v;
int numLocalNodes, globalNodeNumber, i, j, eq;
MatStokesBlockScaling BA = bsscr->BA;
PetscErrorCode ierr;
Mat Amat,Pmat, G;
MatStructure pflag;
PetscFunctionBegin;
/* get G matrix from Amat matrix operator on ksp */
ierr=Stg_PCGetOperators(ksp->pc,&Amat,&Pmat,&pflag);CHKERRQ(ierr);
MatNestGetSubMat( Amat, 0,1, &G );/* G should always exist */
/* now create Vecs t and v to match size of G: i.e. pressure */ /* NOTE: not using "h" vector from ksp->vec_rhs because this part of the block vector doesn't always exist */
MatGetVecs( G, &t, PETSC_NULL );/* t and v are destroyed in KSPDestroy_BSSCR */
MatGetVecs( G, &v, PETSC_NULL );/* t and v such that can do G*t */
numLocalNodes = Mesh_GetLocalSize( feMesh, MT_VERTEX); /* number of nodes on current proc not counting any shadow nodes */
for(j=0;j<numLocalNodes;j++){
i = globalNodeNumber = Mesh_DomainToGlobal( feMesh, MT_VERTEX, j);
RegularMeshUtils_Element_1DTo3D(feMesh, i, ijk);
eq = eq_num->destinationArray[j][0];/* get global equation number -- 2nd arg is always 0 because pressure has only one dof */
if(eq != -1){
if( (ijk[0]+ijk[1]+ijk[2])%2 ==0 ){
VecSetValue(t,eq,1.0,INSERT_VALUES);
}
else{
VecSetValue(v,eq,1.0,INSERT_VALUES); }}}
VecAssemblyBegin( t );
VecAssemblyEnd( t );
VecAssemblyBegin( v );
VecAssemblyEnd( v );
/* Scaling the null vectors here because it easier at the moment *//* maybe should do this in the original scaling function */
if( BA->scaling_exists == PETSC_TRUE ){
Vec R2;
/* Get the scalings out the block mat data */
VecNestGetSubVec( BA->Rz, 1, &R2 );
VecPointwiseDivide( t, t, R2); /* x <- x * 1/R2 */
VecPointwiseDivide( v, v, R2);
}
bsscr_writeVec( t, "t", "Writing t vector");
bsscr_writeVec( v, "v", "Writing v vector");
bsscr->t=t;
bsscr->v=v;
PetscFunctionReturn(0);
}
开发者ID:AngelValverdePerez,项目名称:underworld2,代码行数:55,代码来源:ksp_pressure_nullspace.c
示例3: MatGetVecs
/*@C
KSPGetVecs - Gets a number of work vectors.
Input Parameters:
+ ksp - iterative context
. rightn - number of right work vectors
- leftn - number of left work vectors to allocate
Output Parameter:
+ right - the array of vectors created
- left - the array of left vectors
Note: The right vector has as many elements as the matrix has columns. The left
vector has as many elements as the matrix has rows.
Level: advanced
.seealso: MatGetVecs()
@*/
PetscErrorCode KSPGetVecs(KSP ksp,PetscInt rightn, Vec **right,PetscInt leftn,Vec **left)
{
PetscErrorCode ierr;
Vec vecr,vecl;
PetscFunctionBegin;
if (rightn) {
if (!right) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"You asked for right vectors but did not pass a pointer to hold them");
if (ksp->vec_sol) vecr = ksp->vec_sol;
else {
if (ksp->dm) {
ierr = DMGetGlobalVector(ksp->dm,&vecr);CHKERRQ(ierr);
} else {
Mat mat;
if (!ksp->pc) {ierr = KSPGetPC(ksp,&ksp->pc);CHKERRQ(ierr);}
ierr = PCGetOperators(ksp->pc,&mat,NULL,NULL);CHKERRQ(ierr);
ierr = MatGetVecs(mat,&vecr,NULL);CHKERRQ(ierr);
}
}
ierr = VecDuplicateVecs(vecr,rightn,right);CHKERRQ(ierr);
if (!ksp->vec_sol) {
if (ksp->dm) {
ierr = DMRestoreGlobalVector(ksp->dm,&vecr);CHKERRQ(ierr);
} else {
ierr = VecDestroy(&vecr);CHKERRQ(ierr);
}
}
}
if (leftn) {
if (!left) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"You asked for left vectors but did not pass a pointer to hold them");
if (ksp->vec_rhs) vecl = ksp->vec_rhs;
else {
if (ksp->dm) {
ierr = DMGetGlobalVector(ksp->dm,&vecl);CHKERRQ(ierr);
} else {
Mat mat;
if (!ksp->pc) {ierr = KSPGetPC(ksp,&ksp->pc);CHKERRQ(ierr);}
ierr = PCGetOperators(ksp->pc,&mat,NULL,NULL);CHKERRQ(ierr);
ierr = MatGetVecs(mat,NULL,&vecl);CHKERRQ(ierr);
}
}
ierr = VecDuplicateVecs(vecl,leftn,left);CHKERRQ(ierr);
if (!ksp->vec_rhs) {
if (ksp->dm) {
ierr = DMRestoreGlobalVector(ksp->dm,&vecl);CHKERRQ(ierr);
} else {
ierr = VecDestroy(&vecl);CHKERRQ(ierr);
}
}
}
PetscFunctionReturn(0);
}
开发者ID:ZJLi2013,项目名称:petsc,代码行数:72,代码来源:iterativ.c
示例4: PCSetUp_Eisenstat
static PetscErrorCode PCSetUp_Eisenstat(PC pc)
{
PetscErrorCode ierr;
PetscInt M,N,m,n;
PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
PetscFunctionBegin;
if (!pc->setupcalled) {
ierr = MatGetSize(pc->mat,&M,&N);CHKERRQ(ierr);
ierr = MatGetLocalSize(pc->mat,&m,&n);CHKERRQ(ierr);
ierr = MatCreate(((PetscObject)pc)->comm,&eis->shell);CHKERRQ(ierr);
ierr = MatSetSizes(eis->shell,m,n,M,N);CHKERRQ(ierr);
ierr = MatSetType(eis->shell,MATSHELL);CHKERRQ(ierr);
ierr = MatSetUp(eis->shell);CHKERRQ(ierr);
ierr = MatShellSetContext(eis->shell,(void*)pc);CHKERRQ(ierr);
ierr = PetscLogObjectParent(pc,eis->shell);CHKERRQ(ierr);
ierr = MatShellSetOperation(eis->shell,MATOP_MULT,(void(*)(void))PCMult_Eisenstat);CHKERRQ(ierr);
}
if (!eis->usediag) PetscFunctionReturn(0);
if (!pc->setupcalled) {
ierr = MatGetVecs(pc->pmat,&eis->diag,0);CHKERRQ(ierr);
ierr = PetscLogObjectParent(pc,eis->diag);CHKERRQ(ierr);
}
ierr = MatGetDiagonal(pc->pmat,eis->diag);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
开发者ID:erdc-cm,项目名称:petsc-dev,代码行数:26,代码来源:eisen.c
示例5: ApplyPCPrecPETSCGen
static void ApplyPCPrecPETSCGen(void *x, PRIMME_INT *ldx, void *y,
PRIMME_INT *ldy, int *blockSize, int trans, PC *pc, MPI_Comm comm) {
int i;
Vec xvec, yvec;
Mat matrix;
PetscErrorCode ierr;
PetscInt mLocal, nLocal;
ierr = PCGetOperators(pc[0],&matrix,NULL); CHKERRABORT(comm, ierr);
assert(sizeof(PetscScalar) == sizeof(SCALAR));
ierr = MatGetLocalSize(matrix, &mLocal, &nLocal); CHKERRABORT(comm, ierr);
assert(mLocal == nLocal && nLocal <= *ldx && mLocal <= *ldy);
#if PETSC_VERSION_LT(3,6,0)
ierr = MatGetVecs(matrix, &xvec, &yvec); CHKERRABORT(comm, ierr);
#else
ierr = MatCreateVecs(matrix, &xvec, &yvec); CHKERRABORT(comm, ierr);
#endif
for (i=0; i<*blockSize; i++) {
ierr = VecPlaceArray(xvec, ((SCALAR*)x) + (*ldx)*i); CHKERRABORT(comm, ierr);
ierr = VecPlaceArray(yvec, ((SCALAR*)y) + (*ldy)*i); CHKERRABORT(comm, ierr);
if (trans == 0) {
ierr = PCApply(*pc, xvec, yvec); CHKERRABORT(comm, ierr);
} else if (pc[1]) {
ierr = PCApply(pc[1], xvec, yvec); CHKERRABORT(comm, ierr);
} else {
ierr = PCApplyTranspose(pc[0], xvec, yvec);
}
ierr = VecResetArray(xvec); CHKERRABORT(comm, ierr);
ierr = VecResetArray(yvec); CHKERRABORT(comm, ierr);
}
ierr = VecDestroy(&xvec); CHKERRABORT(comm, ierr);
ierr = VecDestroy(&yvec); CHKERRABORT(comm, ierr);
}
开发者ID:primme,项目名称:primme,代码行数:35,代码来源:petscw.c
示例6: createOuterContext
void createOuterContext(OuterContext* & ctx) {
int rank;
int npes;
MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
MPI_Comm_size(PETSC_COMM_WORLD, &npes);
ctx = new OuterContext;
ctx->data = NULL;
ctx->root = NULL;
ctx->outerMat = PETSC_NULL;
ctx->outerKsp = PETSC_NULL;
ctx->outerPC = PETSC_NULL;
ctx->outerSol = PETSC_NULL;
ctx->outerRhs = PETSC_NULL;
createLocalData(ctx->data);
createRSDtree(ctx->root, rank, npes);
createOuterMat(ctx);
createOuterPC(ctx);
createOuterKsp(ctx);
MatGetVecs(ctx->outerMat, &(ctx->outerSol), &(ctx->outerRhs));
}
开发者ID:hsundar,项目名称:schur,代码行数:28,代码来源:schurSetup.C
示例7: BSSCR_MatStokesMVBlockDefaultBuildScaling
// updated
PetscErrorCode BSSCR_MatStokesMVBlockDefaultBuildScaling( MatStokesBlockScaling BA, Mat A, Vec b, Vec x, PetscTruth sym )
{
Mat K,G,D,C;
Vec rG;
PetscScalar rg2, rg, ra;
PetscInt N;
Vec rA, rC;
Vec L1,L2, R1,R2;
Mat S;
VecNestGetSubVec( BA->Lz, 0, &L1 );
VecNestGetSubVec( BA->Lz, 1, &L2 );
VecNestGetSubVec( BA->Rz, 0, &R1 );
VecNestGetSubVec( BA->Rz, 1, &R2 );
rA = L1;
rC = L2;
MatNestGetSubMat( A, 0,0, &K );
MatNestGetSubMat( A, 0,1, &G );
MatNestGetSubMat( A, 1,0, &D );
MatNestGetSubMat( A, 1,1, &C );
VecDuplicate( rA, &rG );
/* Get magnitude of K */
//px_MatGetAbsRowSum( K, rA );
//MatGetRowMax( K, rA, PETSC_NULL );
MatGetDiagonal( K, rA);
VecSqrt( rA );
VecReciprocal( rA );
/* VecDot( rA,rA, &ra ); */
/* VecGetSize( rA, &N ); */
/* ra = PetscSqrtScalar( ra/N ); */
/* Get magnitude of G */
//px_MatGetAbsRowSum( G, rG );
//MatGetRowMax( G, rG, PETSC_NULL );
Mat A21_cpy;
Mat Shat;
Vec diag; /* same as rA*rA */
MatGetVecs( K, &diag, PETSC_NULL );
MatGetDiagonal( K, diag );
VecReciprocal( diag );
if( sym )
#if( PETSC_VERSION_MAJOR <= 2 )
MatTranspose( G, &A21_cpy );
#else
MatTranspose( G, MAT_INITIAL_MATRIX, &A21_cpy );
#endif
else {
开发者ID:AngelValverdePerez,项目名称:underworld2,代码行数:61,代码来源:stokes_mvblock_scaling.c
示例8: NonlocalCollection
NonlocalCollection(Mat depends_on, IS interest_set) {
find_influences(depends_on, interest_set, &nodes);
PetscInt local_size;
ISGetLocalSize(nodes, &local_size);
VecCreateMPI(PETSC_COMM_WORLD, local_size, PETSC_DECIDE, &vec);
IS onto_index_set;
describe_partition(vec, &onto_index_set);
PetscInt begin;
PetscInt end;
VecGetOwnershipRange(vec, &begin, &end);
PetscInt *indicies;
ISGetIndices(nodes, &indicies);
assert(local_size == end-begin);
for (int ii=0; ii<local_size; ii++) {
map[indicies[ii]] = ii+begin;
}
ISRestoreIndices(nodes, &indicies);
Vec w;
MatGetVecs(depends_on, PETSC_NULL, &w);
VecScatterCreate(w, nodes, vec, onto_index_set, &scatter);
VecDestroy(w);
ISDestroy(onto_index_set);
}
开发者ID:rblake,项目名称:petsc_amg,代码行数:25,代码来源:mglib.c
示例9: MatMultTranspose_Composite_Multiplicative
PetscErrorCode MatMultTranspose_Composite_Multiplicative(Mat A,Vec x,Vec y)
{
Mat_Composite *shell = (Mat_Composite*)A->data;
Mat_CompositeLink tail = shell->tail;
PetscErrorCode ierr;
Vec in,out;
PetscFunctionBegin;
if (!tail) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()");
in = x;
if (shell->left) {
if (!shell->leftwork) {
ierr = VecDuplicate(shell->left,&shell->leftwork);CHKERRQ(ierr);
}
ierr = VecPointwiseMult(shell->leftwork,shell->left,in);CHKERRQ(ierr);
in = shell->leftwork;
}
while (tail->prev) {
if (!tail->prev->work) { /* should reuse previous work if the same size */
ierr = MatGetVecs(tail->mat,NULL,&tail->prev->work);CHKERRQ(ierr);
}
out = tail->prev->work;
ierr = MatMultTranspose(tail->mat,in,out);CHKERRQ(ierr);
in = out;
tail = tail->prev;
}
ierr = MatMultTranspose(tail->mat,in,y);CHKERRQ(ierr);
if (shell->right) {
ierr = VecPointwiseMult(y,shell->right,y);CHKERRQ(ierr);
}
ierr = VecScale(y,shell->scale);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
开发者ID:feelpp,项目名称:debian-petsc,代码行数:33,代码来源:mcomposite.c
示例10: PETScMatvecGenNoBlock
static void PETScMatvecGenNoBlock(void *x, PRIMME_INT ldx, void *y, PRIMME_INT ldy,
int blockSize, int trans, Mat matrix, MPI_Comm comm) {
int i;
Vec xvec, yvec;
PetscInt m, n, mLocal, nLocal;
PetscErrorCode ierr;
assert(sizeof(PetscScalar) == sizeof(SCALAR));
ierr = MatGetSize(matrix, &m, &n); CHKERRABORT(comm, ierr);
ierr = MatGetLocalSize(matrix, &mLocal, &nLocal); CHKERRABORT(comm, ierr);
#if PETSC_VERSION_LT(3,6,0)
ierr = MatGetVecs(matrix, &xvec, &yvec); CHKERRABORT(comm, ierr);
#else
ierr = MatCreateVecs(matrix, &xvec, &yvec); CHKERRABORT(comm, ierr);
#endif
if (trans == 1) {
Vec aux = xvec; xvec = yvec; yvec = aux;
}
for (i=0; i<blockSize; i++) {
ierr = VecPlaceArray(xvec, ((SCALAR*)x) + ldx*i); CHKERRABORT(comm, ierr);
ierr = VecPlaceArray(yvec, ((SCALAR*)y) + ldy*i); CHKERRABORT(comm, ierr);
if (trans == 0) {
ierr = MatMult(matrix, xvec, yvec); CHKERRABORT(comm, ierr);
} else {
ierr = MatMultHermitianTranspose(matrix, xvec, yvec); CHKERRABORT(comm, ierr);
}
ierr = VecResetArray(xvec); CHKERRABORT(comm, ierr);
ierr = VecResetArray(yvec); CHKERRABORT(comm, ierr);
}
ierr = VecDestroy(&xvec); CHKERRABORT(comm, ierr);
ierr = VecDestroy(&yvec); CHKERRABORT(comm, ierr);
}
开发者ID:primme,项目名称:primme,代码行数:33,代码来源:petscw.c
示例11: computeMaxEigVal
PetscErrorCode computeMaxEigVal(Mat A, PetscInt its, PetscScalar *eig)
{
PetscErrorCode ierr;
PetscRandom rctx; /* random number generator context */
Vec x0, x, x_1, tmp;
PetscScalar lambda_its, lambda_its_1;
PetscInt i;
PetscFunctionBeginUser;
ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rctx);CHKERRQ(ierr);
ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr);
ierr = MatGetVecs(A, &x_1, &x);CHKERRQ(ierr);
ierr = VecSetRandom(x, rctx);CHKERRQ(ierr);
ierr = VecDuplicate(x, &x0);CHKERRQ(ierr);
ierr = VecCopy(x, x0);CHKERRQ(ierr);
ierr = MatMult(A, x, x_1);CHKERRQ(ierr);
for (i=0; i<its; i++) {
tmp = x; x = x_1; x_1 = tmp;
ierr = MatMult(A, x, x_1);CHKERRQ(ierr);
}
ierr = VecDot(x0, x, &lambda_its);CHKERRQ(ierr);
ierr = VecDot(x0, x_1, &lambda_its_1);CHKERRQ(ierr);
*eig = lambda_its_1/lambda_its;
ierr = VecDestroy(&x0);CHKERRQ(ierr);
ierr = VecDestroy(&x);CHKERRQ(ierr);
ierr = VecDestroy(&x_1);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
开发者ID:hsahasra,项目名称:petsc-magma-dense-mat,代码行数:31,代码来源:ex39.c
示例12: level
/*@
PCMGGetRScale - Gets the pointwise scaling for the restriction operator from level l to l-1.
Collective on PC
Input Parameters:
+ pc - the multigrid context
. rscale - the scaling
- l - the level (0 is coarsest) to supply [Do not supply 0]
Level: advanced
Notes:
When evaluating a function on a coarse level one does not want to do F(R * x) one does F(rscale * R * x) where rscale is 1 over the row sums of R.
.keywords: MG, set, multigrid, restriction, level
.seealso: PCMGSetInterpolation(), PCMGGetRestriction()
@*/
PetscErrorCode PCMGGetRScale(PC pc,PetscInt l,Vec *rscale)
{
PetscErrorCode ierr;
PC_MG *mg = (PC_MG*)pc->data;
PC_MG_Levels **mglevels = mg->levels;
PetscFunctionBegin;
PetscValidHeaderSpecific(pc,PC_CLASSID,1);
if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
if (l <= 0 || mg->nlevels <= l) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Level %D must be in range {1,...,%D}",l,mg->nlevels-1);
if (!mglevels[l]->rscale) {
Mat R;
Vec X,Y,coarse,fine;
PetscInt M,N;
ierr = PCMGGetRestriction(pc,l,&R);CHKERRQ(ierr);
ierr = MatGetVecs(R,&X,&Y);CHKERRQ(ierr);
ierr = MatGetSize(R,&M,&N);CHKERRQ(ierr);
if (M < N) {
fine = X;
coarse = Y;
} else if (N < M) {
fine = Y; coarse = X;
} else SETERRQ(PetscObjectComm((PetscObject)R),PETSC_ERR_SUP,"Restriction matrix is square, cannot determine which Vec is coarser");
ierr = VecSet(fine,1.);CHKERRQ(ierr);
ierr = MatRestrict(R,fine,coarse);CHKERRQ(ierr);
ierr = VecDestroy(&fine);CHKERRQ(ierr);
ierr = VecReciprocal(coarse);CHKERRQ(ierr);
mglevels[l]->rscale = coarse;
}
*rscale = mglevels[l]->rscale;
PetscFunctionReturn(0);
}
开发者ID:prbrune,项目名称:petsc,代码行数:51,代码来源:mgfunc.c
示例13: MatMult_SchurComplement
PetscErrorCode MatMult_SchurComplement(Mat N,Vec x,Vec y)
{
Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
PetscErrorCode ierr;
PetscFunctionBegin;
if (!Na->work1) {ierr = MatGetVecs(Na->A,&Na->work1,PETSC_NULL);CHKERRQ(ierr);}
if (!Na->work2) {ierr = MatGetVecs(Na->A,&Na->work2,PETSC_NULL);CHKERRQ(ierr);}
ierr = MatMult(Na->B,x,Na->work1);CHKERRQ(ierr);
ierr = KSPSolve(Na->ksp,Na->work1,Na->work2);CHKERRQ(ierr);
ierr = MatMult(Na->C,Na->work2,y);CHKERRQ(ierr);
ierr = VecScale(y,-1.0);CHKERRQ(ierr);
if (Na->D) {
ierr = MatMultAdd(Na->D,x,y,y);CHKERRQ(ierr);
}
PetscFunctionReturn(0);
}
开发者ID:erdc-cm,项目名称:petsc-dev,代码行数:17,代码来源:schurm.c
示例14: z
/*
PEPLinearExtract_Norm - Auxiliary routine that copies the solution of the
linear eigenproblem to the PEP object. The eigenvector of the generalized
problem is supposed to be
z = [ x ]
[ l*x ]
If |l|<1.0, the eigenvector is taken from z(1:n), otherwise from z(n+1:2*n).
Finally, x is normalized so that ||x||_2 = 1.
*/
static PetscErrorCode PEPLinearExtract_Norm(PEP pep,EPS eps)
{
PetscErrorCode ierr;
PetscInt i,offset;
PetscScalar *px;
Vec xr,xi,w,vi;
#if !defined(PETSC_USE_COMPLEX)
Vec vi1;
#endif
Mat A;
PetscFunctionBegin;
ierr = EPSGetOperators(eps,&A,NULL);CHKERRQ(ierr);
ierr = MatGetVecs(A,&xr,NULL);CHKERRQ(ierr);
ierr = VecDuplicate(xr,&xi);CHKERRQ(ierr);
ierr = VecCreateMPIWithArray(PetscObjectComm((PetscObject)pep),1,pep->nloc,pep->n,NULL,&w);CHKERRQ(ierr);
for (i=0;i<pep->nconv;i++) {
ierr = EPSGetEigenpair(eps,i,&pep->eigr[i],&pep->eigi[i],xr,xi);CHKERRQ(ierr);
pep->eigr[i] *= pep->sfactor;
pep->eigi[i] *= pep->sfactor;
if (SlepcAbsEigenvalue(pep->eigr[i],pep->eigi[i])>1.0) offset = pep->nloc;
else offset = 0;
#if !defined(PETSC_USE_COMPLEX)
if (pep->eigi[i]>0.0) { /* first eigenvalue of a complex conjugate pair */
ierr = VecGetArray(xr,&px);CHKERRQ(ierr);
ierr = VecPlaceArray(w,px+offset);CHKERRQ(ierr);
ierr = BVInsertVec(pep->V,i,w);CHKERRQ(ierr);
ierr = VecResetArray(w);CHKERRQ(ierr);
ierr = VecRestoreArray(xr,&px);CHKERRQ(ierr);
ierr = VecGetArray(xi,&px);CHKERRQ(ierr);
ierr = VecPlaceArray(w,px+offset);CHKERRQ(ierr);
ierr = BVInsertVec(pep->V,i+1,w);CHKERRQ(ierr);
ierr = VecResetArray(w);CHKERRQ(ierr);
ierr = VecRestoreArray(xi,&px);CHKERRQ(ierr);
ierr = BVGetColumn(pep->V,i,&vi);CHKERRQ(ierr);
ierr = BVGetColumn(pep->V,i+1,&vi1);CHKERRQ(ierr);
ierr = SlepcVecNormalize(vi,vi1,PETSC_TRUE,NULL);CHKERRQ(ierr);
ierr = BVRestoreColumn(pep->V,i,&vi);CHKERRQ(ierr);
ierr = BVRestoreColumn(pep->V,i+1,&vi1);CHKERRQ(ierr);
} else if (pep->eigi[i]==0.0) /* real eigenvalue */
#endif
{
ierr = VecGetArray(xr,&px);CHKERRQ(ierr);
ierr = VecPlaceArray(w,px+offset);CHKERRQ(ierr);
ierr = BVInsertVec(pep->V,i,w);CHKERRQ(ierr);
ierr = VecResetArray(w);CHKERRQ(ierr);
ierr = VecRestoreArray(xr,&px);CHKERRQ(ierr);
ierr = BVGetColumn(pep->V,i,&vi);CHKERRQ(ierr);
ierr = SlepcVecNormalize(vi,NULL,PETSC_FALSE,NULL);CHKERRQ(ierr);
ierr = BVRestoreColumn(pep->V,i,&vi);CHKERRQ(ierr);
}
}
ierr = VecDestroy(&w);CHKERRQ(ierr);
ierr = VecDestroy(&xr);CHKERRQ(ierr);
ierr = VecDestroy(&xi);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
开发者ID:OpenCMISS-Dependencies,项目名称:slepc,代码行数:66,代码来源:linear.c
示例15: MatGetVecs_SchurComplement
PetscErrorCode MatGetVecs_SchurComplement(Mat N,Vec *right,Vec *left)
{
Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
PetscErrorCode ierr;
PetscFunctionBegin;
if (Na->D) {
ierr = MatGetVecs(Na->D,right,left);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
if (right) {
ierr = MatGetVecs(Na->B,right,PETSC_NULL);CHKERRQ(ierr);
}
if (left) {
ierr = MatGetVecs(Na->C,PETSC_NULL,left);CHKERRQ(ierr);
}
PetscFunctionReturn(0);
}
开发者ID:erdc-cm,项目名称:petsc-dev,代码行数:18,代码来源:schurm.c
示例16: RunTest
PetscErrorCode RunTest(void)
{
PetscInt N = 100;
PetscBool draw = PETSC_FALSE;
PetscReal rnorm;
Mat A;
Vec b,x,r;
KSP ksp;
PC pc;
PetscErrorCode ierr;
PetscFunctionBegin;
ierr = PetscOptionsGetInt(0,"-N",&N,NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetBool(0,"-draw",&draw,NULL);CHKERRQ(ierr);
ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N);CHKERRQ(ierr);
ierr = MatSetType(A,MATPYTHON);CHKERRQ(ierr);
ierr = MatPythonSetType(A,"example1.py:Laplace1D");CHKERRQ(ierr);
ierr = MatSetUp(A);CHKERRQ(ierr);
ierr = MatGetVecs(A,&x,&b);CHKERRQ(ierr);
ierr = VecSet(b,1);CHKERRQ(ierr);
ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
ierr = KSPSetType(ksp,KSPPYTHON);CHKERRQ(ierr);
ierr = KSPPythonSetType(ksp,"example1.py:ConjGrad");CHKERRQ(ierr);
ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
ierr = PCSetType(pc,PCPYTHON);CHKERRQ(ierr);
ierr = PCPythonSetType(pc,"example1.py:Jacobi");CHKERRQ(ierr);
ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
ierr = VecDuplicate(b,&r);CHKERRQ(ierr);
ierr = MatMult(A,x,r);CHKERRQ(ierr);
ierr = VecAYPX(r,-1,b);CHKERRQ(ierr);
ierr = VecNorm(r,NORM_2,&rnorm);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"error norm = %g\n",rnorm);CHKERRQ(ierr);
if (draw) {
ierr = VecView(x,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);
ierr = PetscSleep(2);CHKERRQ(ierr);
}
ierr = VecDestroy(&x);CHKERRQ(ierr);
ierr = VecDestroy(&b);CHKERRQ(ierr);
ierr = VecDestroy(&r);CHKERRQ(ierr);
ierr = MatDestroy(&A);CHKERRQ(ierr);
ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
开发者ID:00liujj,项目名称:petsc,代码行数:56,代码来源:ex1.c
示例17: MatSetLocalToGlobalMapping_IS
PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping rmapping,ISLocalToGlobalMapping cmapping)
{
PetscErrorCode ierr;
PetscInt n,bs;
Mat_IS *is = (Mat_IS*)A->data;
IS from,to;
Vec global;
PetscFunctionBegin;
PetscCheckSameComm(A,1,rmapping,2);
if (rmapping != cmapping) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"MATIS requires the row and column mappings to be identical");
if (is->mapping) { /* Currenly destroys the objects that will be created by this routine. Is there anything else that should be checked? */
ierr = ISLocalToGlobalMappingDestroy(&is->mapping);CHKERRQ(ierr);
ierr = VecDestroy(&is->x);CHKERRQ(ierr);
ierr = VecDestroy(&is->y);CHKERRQ(ierr);
ierr = VecScatterDestroy(&is->ctx);CHKERRQ(ierr);
ierr = MatDestroy(&is->A);CHKERRQ(ierr);
}
ierr = PetscObjectReference((PetscObject)rmapping);CHKERRQ(ierr);
ierr = ISLocalToGlobalMappingDestroy(&is->mapping);CHKERRQ(ierr);
is->mapping = rmapping;
/*
ierr = PetscLayoutSetISLocalToGlobalMapping(A->rmap,rmapping);CHKERRQ(ierr);
ierr = PetscLayoutSetISLocalToGlobalMapping(A->cmap,cmapping);CHKERRQ(ierr);
*/
/* Create the local matrix A */
ierr = ISLocalToGlobalMappingGetSize(rmapping,&n);CHKERRQ(ierr);
ierr = MatGetBlockSize(A,&bs);CHKERRQ(ierr);
ierr = MatCreate(PETSC_COMM_SELF,&is->A);CHKERRQ(ierr);
ierr = MatSetSizes(is->A,n,n,n,n);CHKERRQ(ierr);
ierr = MatSetBlockSize(is->A,bs);CHKERRQ(ierr);
ierr = MatSetOptionsPrefix(is->A,((PetscObject)A)->prefix);CHKERRQ(ierr);
ierr = MatAppendOptionsPrefix(is->A,"is_");CHKERRQ(ierr);
ierr = MatSetFromOptions(is->A);CHKERRQ(ierr);
/* Create the local work vectors */
ierr = VecCreate(PETSC_COMM_SELF,&is->x);CHKERRQ(ierr);
ierr = VecSetBlockSize(is->x,bs);CHKERRQ(ierr);
ierr = VecSetSizes(is->x,n,n);CHKERRQ(ierr);
ierr = VecSetOptionsPrefix(is->x,((PetscObject)A)->prefix);CHKERRQ(ierr);
ierr = VecAppendOptionsPrefix(is->x,"is_");CHKERRQ(ierr);
ierr = VecSetFromOptions(is->x);CHKERRQ(ierr);
ierr = VecDuplicate(is->x,&is->y);CHKERRQ(ierr);
/* setup the global to local scatter */
ierr = ISCreateStride(PETSC_COMM_SELF,n,0,1,&to);CHKERRQ(ierr);
ierr = ISLocalToGlobalMappingApplyIS(rmapping,to,&from);CHKERRQ(ierr);
ierr = MatGetVecs(A,&global,NULL);CHKERRQ(ierr);
ierr = VecScatterCreate(global,from,is->x,to,&is->ctx);CHKERRQ(ierr);
ierr = VecDestroy(&global);CHKERRQ(ierr);
ierr = ISDestroy(&to);CHKERRQ(ierr);
ierr = ISDestroy(&from);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
开发者ID:OpenCMISS-Dependencies,项目名称:petsc,代码行数:55,代码来源:matis.c
示例18: MatComputeExplicitOperator
/*@
MatComputeExplicitOperator - Computes the explicit matrix
Collective on Mat
Input Parameter:
. inmat - the matrix
Output Parameter:
. mat - the explict preconditioned operator
Notes:
This computation is done by applying the operators to columns of the
identity matrix.
Currently, this routine uses a dense matrix format when 1 processor
is used and a sparse format otherwise. This routine is costly in general,
and is recommended for use only with relatively small systems.
Level: advanced
.keywords: Mat, compute, explicit, operator
@*/
PetscErrorCode MatComputeExplicitOperator(Mat inmat,Mat *mat)
{
Vec in,out;
PetscErrorCode ierr;
PetscInt i,m,n,M,N,*rows,start,end;
MPI_Comm comm;
PetscScalar *array,zero = 0.0,one = 1.0;
PetscMPIInt size;
PetscFunctionBegin;
PetscValidHeaderSpecific(inmat,MAT_CLASSID,1);
PetscValidPointer(mat,2);
ierr = PetscObjectGetComm((PetscObject)inmat,&comm);CHKERRQ(ierr);
ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
ierr = MatGetLocalSize(inmat,&m,&n);CHKERRQ(ierr);
ierr = MatGetSize(inmat,&M,&N);CHKERRQ(ierr);
ierr = MatGetVecs(inmat,&in,&out);CHKERRQ(ierr);
ierr = VecSetOption(in,VEC_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr);
ierr = VecGetOwnershipRange(out,&start,&end);CHKERRQ(ierr);
ierr = PetscMalloc(m*sizeof(PetscInt),&rows);CHKERRQ(ierr);
for (i=0; i<m; i++) rows[i] = start + i;
ierr = MatCreate(comm,mat);CHKERRQ(ierr);
ierr = MatSetSizes(*mat,m,n,M,N);CHKERRQ(ierr);
if (size == 1) {
ierr = MatSetType(*mat,MATSEQDENSE);CHKERRQ(ierr);
ierr = MatSeqDenseSetPreallocation(*mat,NULL);CHKERRQ(ierr);
} else {
ierr = MatSetType(*mat,MATMPIAIJ);CHKERRQ(ierr);
ierr = MatMPIAIJSetPreallocation(*mat,n,NULL,N-n,NULL);CHKERRQ(ierr);
}
for (i=0; i<N; i++) {
ierr = VecSet(in,zero);CHKERRQ(ierr);
ierr = VecSetValues(in,1,&i,&one,INSERT_VALUES);CHKERRQ(ierr);
ierr = VecAssemblyBegin(in);CHKERRQ(ierr);
ierr = VecAssemblyEnd(in);CHKERRQ(ierr);
ierr = MatMult(inmat,in,out);CHKERRQ(ierr);
ierr = VecGetArray(out,&array);CHKERRQ(ierr);
ierr = MatSetValues(*mat,m,rows,1,&i,array,INSERT_VALUES);CHKERRQ(ierr);
ierr = VecRestoreArray(out,&array);CHKERRQ(ierr);
}
ierr = PetscFree(rows);CHKERRQ(ierr);
ierr = VecDestroy(&out);CHKERRQ(ierr);
ierr = VecDestroy(&in);CHKERRQ(ierr);
ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
开发者ID:feelpp,项目名称:debian-petsc,代码行数:78,代码来源:axpy.c
示例19: PCSetUp_Jacobi_NonSymmetric
static PetscErrorCode PCSetUp_Jacobi_NonSymmetric(PC pc)
{
PetscErrorCode ierr;
PC_Jacobi *jac = (PC_Jacobi*)pc->data;
PetscFunctionBegin;
ierr = MatGetVecs(pc->pmat,&jac->diag,0);CHKERRQ(ierr);
ierr = PetscLogObjectParent(pc,jac->diag);CHKERRQ(ierr);
ierr = PCSetUp_Jacobi(pc);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
开发者ID:Kun-Qu,项目名称:petsc,代码行数:11,代码来源:jacobi.c
示例20: PCGAMGOptProl_Classical_Jacobi
PetscErrorCode PCGAMGOptProl_Classical_Jacobi(PC pc,const Mat A,Mat *P)
{
PetscErrorCode ierr;
PetscInt f,s,n,cf,cs,i,idx;
PetscInt *coarserows;
PetscInt ncols;
const PetscInt *pcols;
const PetscScalar *pvals;
Mat Pnew;
Vec diag;
PC_MG *mg = (PC_MG*)pc->data;
PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
PC_GAMG_Classical *cls = (PC_GAMG_Classical*)pc_gamg->subctx;
PetscFunctionBegin;
if (cls->nsmooths == 0) {
ierr = PCGAMGTruncateProlongator_Private(pc,P);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
ierr = MatGetOwnershipRange(*P,&s,&f);CHKERRQ(ierr);
n = f-s;
ierr = MatGetOwnershipRangeColumn(*P,&cs,&cf);CHKERRQ(ierr);
ierr = PetscMalloc(sizeof(PetscInt)*n,&coarserows);CHKERRQ(ierr);
/* identify the rows corresponding to coarse unknowns */
idx = 0;
for (i=s;i<f;i++) {
ierr = MatGetRow(*P,i,&ncols,&pcols,&pvals);CHKERRQ(ierr);
/* assume, for now, that it's a coarse unknown if it has a single unit entry */
if (ncols == 1) {
if (pvals[0] == 1.) {
coarserows[idx] = i;
idx++;
}
}
ierr = MatRestoreRow(*P,i,&ncols,&pcols,&pvals);CHKERRQ(ierr);
}
ierr = MatGetVecs(A,&diag,0);CHKERRQ(ierr);
ierr = MatGetDiagonal(A,diag);CHKERRQ(ierr);
ierr = VecReciprocal(diag);CHKERRQ(ierr);
for (i=0;i<cls->nsmooths;i++) {
ierr = MatMatMult(A,*P,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&Pnew);CHKERRQ(ierr);
ierr = MatZeroRows(Pnew,idx,coarserows,0.,NULL,NULL);CHKERRQ(ierr);
ierr = MatDiagonalScale(Pnew,diag,0);CHKERRQ(ierr);
ierr = MatAYPX(Pnew,-1.0,*P,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
ierr = MatDestroy(P);CHKERRQ(ierr);
*P = Pnew;
Pnew = NULL;
}
ierr = VecDestroy(&diag);CHKERRQ(ierr);
ierr = PetscFree(coarserows);CHKERRQ(ierr);
ierr = PCGAMGTruncateProlongator_Private(pc,P);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
开发者ID:ZJLi2013,项目名称:petsc,代码行数:54,代码来源:classical.c
注:本文中的MatGetVecs函数示例由纯净天空整理自Github/MSDocs等源码及文档管理平台,相关代码片段筛选自各路编程大神贡献的开源项目,源码版权归原作者所有,传播和使用请参考对应项目的License;未经允许,请勿转载。 |
请发表评论