/*@
DMGlobalToLocalSolve - Solve for the global vector that is mapped to a given local vector by DMGlobalToLocalBegin()/DMGlobalToLocalEnd() with mode
= INSERT_VALUES. It is assumed that the sum of all the local vector sizes is greater than or equal to the global vector size, so the solution is
a least-squares solution. It is also assumed that DMLocalToGlobalBegin()/DMLocalToGlobalEnd() with mode = ADD_VALUES is the adjoint of the
global-to-local map, so that the least-squares solution may be found by the normal equations.
collective
Input Parameters:
+ dm - The DM object
. x - The local vector
- y - The global vector: the input value of globalVec is used as an initial guess
Output Parameters:
. y - The least-squares solution
Level: advanced
Note: If the DM is of type DMPLEX, then y is the solution of L' * D * L * y = L' * D * x, where D is a diagonal mask that is 1 for every point in
the union of the closures of the local cells and 0 otherwise. This difference is only relevant if there are anchor points that are not in the
closure of any local cell (see DMPlexGetAnchors()/DMPlexSetAnchors()).
.seealso: DMGlobalToLocalBegin(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin(), DMLocalToGlobalEnd(), DMPlexGetAnchors(), DMPlexSetAnchors()
@*/
PetscErrorCode DMGlobalToLocalSolve(DM dm, Vec x, Vec y)
{
Mat CtC;
PetscInt n, N, cStart, cEnd, c;
PetscBool isPlex;
KSP ksp;
PC pc;
Vec global, mask=NULL;
projectConstraintsCtx ctx;
PetscErrorCode ierr;
PetscFunctionBegin;
ierr = PetscObjectTypeCompare((PetscObject)dm,DMPLEX,&isPlex);CHKERRQ(ierr);
if (isPlex) {
/* mark points in the closure */
ierr = DMGetLocalVector(dm,&mask);CHKERRQ(ierr);
ierr = VecSet(mask,0.0);CHKERRQ(ierr);
ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr);
if (cEnd > cStart) {
PetscScalar *ones;
PetscInt numValues, i;
ierr = DMPlexVecGetClosure(dm,NULL,mask,cStart,&numValues,NULL);CHKERRQ(ierr);
ierr = PetscMalloc1(numValues,&ones);CHKERRQ(ierr);
for (i = 0; i < numValues; i++) {
ones[i] = 1.;
}
for (c = cStart; c < cEnd; c++) {
ierr = DMPlexVecSetClosure(dm,NULL,mask,c,ones,INSERT_VALUES);CHKERRQ(ierr);
}
ierr = PetscFree(ones);CHKERRQ(ierr);
}
}
ctx.dm = dm;
ctx.mask = mask;
ierr = VecGetSize(y,&N);CHKERRQ(ierr);
ierr = VecGetLocalSize(y,&n);CHKERRQ(ierr);
ierr = MatCreate(PetscObjectComm((PetscObject)dm),&CtC);CHKERRQ(ierr);
ierr = MatSetSizes(CtC,n,n,N,N);CHKERRQ(ierr);
ierr = MatSetType(CtC,MATSHELL);CHKERRQ(ierr);
ierr = MatSetUp(CtC);CHKERRQ(ierr);
ierr = MatShellSetContext(CtC,&ctx);CHKERRQ(ierr);
ierr = MatShellSetOperation(CtC,MATOP_MULT,(void(*)(void))MatMult_GlobalToLocalNormal);CHKERRQ(ierr);
ierr = KSPCreate(PetscObjectComm((PetscObject)dm),&ksp);CHKERRQ(ierr);
ierr = KSPSetOperators(ksp,CtC,CtC);CHKERRQ(ierr);
ierr = KSPSetType(ksp,KSPCG);CHKERRQ(ierr);
ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr);
ierr = KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);CHKERRQ(ierr);
ierr = KSPSetUp(ksp);CHKERRQ(ierr);
ierr = DMGetGlobalVector(dm,&global);CHKERRQ(ierr);
ierr = VecSet(global,0.);CHKERRQ(ierr);
if (mask) {ierr = VecPointwiseMult(x,mask,x);CHKERRQ(ierr);}
ierr = DMLocalToGlobalBegin(dm,x,ADD_VALUES,global);CHKERRQ(ierr);
ierr = DMLocalToGlobalEnd(dm,x,ADD_VALUES,global);CHKERRQ(ierr);
ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
ierr = KSPSolve(ksp,global,y);CHKERRQ(ierr);
ierr = DMRestoreGlobalVector(dm,&global);CHKERRQ(ierr);
/* clean up */
if (isPlex) {
ierr = DMRestoreLocalVector(dm,&mask);CHKERRQ(ierr);
}
ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
ierr = MatDestroy(&CtC);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
/*@C
SNESUpdateCheckJacobian - Checks each Jacobian computed by the nonlinear solver comparing the users function with a finite difference computation.
Options Database:
+ -snes_check_jacobian - use this every time SNESSolve() is called
- -snes_check_jacobian_view - Display difference between Jacobian approximated by finite-differencing and the hand-coded Jacobian
Output:
+ difference - ||J - Jd||, the norm of the difference of the hand-coded Jacobian J and the approximate Jacobian Jd obtained by finite-differencing
the residual,
- ratio - ||J - Jd||/||J||, the ratio of the norms of the above difference and the hand-coded Jacobian.
Notes:
Frobenius norm is used in the above throughout. This check is carried out every SNES iteration.
Level: intermediate
.seealso: SNESTEST, SNESCreate(), SNES, SNESSetType(), SNESNEWTONLS, SNESNEWTONTR, SNESSolve()
@*/
PetscErrorCode SNESUpdateCheckJacobian(SNES snes,PetscInt it)
{
Mat A = snes->jacobian,B;
Vec x = snes->vec_sol,f = snes->vec_func,f1 = snes->vec_sol_update;
PetscErrorCode ierr;
PetscReal nrm,gnorm;
PetscErrorCode (*objective)(SNES,Vec,PetscReal*,void*);
void *ctx;
PetscReal fnorm,f1norm,dnorm;
PetscInt m,n,M,N;
PetscBool complete_print = PETSC_FALSE;
void *functx;
PetscViewer viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes));
PetscFunctionBegin;
ierr = PetscOptionsHasName(((PetscObject)snes)->prefix,"-snes_check_jacobian_view",&complete_print);CHKERRQ(ierr);
if (A != snes->jacobian_pre) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Cannot check Jacobian with alternative preconditioner");
ierr = PetscViewerASCIIAddTab(viewer,((PetscObject)snes)->tablevel);CHKERRQ(ierr);
ierr = PetscViewerASCIIPrintf(viewer," Testing hand-coded Jacobian, if the ratio is O(1.e-8), the hand-coded Jacobian is probably correct.\n");CHKERRQ(ierr);
if (!complete_print) {
ierr = PetscViewerASCIIPrintf(viewer," Run with -snes_check_jacobian_view [viewer][:filename][:format] to show difference of hand-coded and finite difference Jacobian.\n");CHKERRQ(ierr);
}
/* compute both versions of Jacobian */
ierr = SNESComputeJacobian(snes,x,A,A);CHKERRQ(ierr);
ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr);
ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr);
ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr);
ierr = MatSetSizes(B,m,n,M,N);CHKERRQ(ierr);
ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr);
ierr = MatSetUp(B);CHKERRQ(ierr);
ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr);
ierr = SNESGetFunction(snes,NULL,NULL,&functx);CHKERRQ(ierr);
ierr = SNESComputeJacobianDefault(snes,x,B,B,functx);CHKERRQ(ierr);
if (complete_print) {
ierr = PetscViewerASCIIPrintf(viewer," Finite difference Jacobian\n");CHKERRQ(ierr);
ierr = MatViewFromOptions(B,(PetscObject)snes,"-snes_check_jacobian_view");CHKERRQ(ierr);
}
/* compare */
ierr = MatAYPX(B,-1.0,A,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
ierr = MatNorm(B,NORM_FROBENIUS,&nrm);CHKERRQ(ierr);
ierr = MatNorm(A,NORM_FROBENIUS,&gnorm);CHKERRQ(ierr);
if (complete_print) {
ierr = PetscViewerASCIIPrintf(viewer," Hand-coded Jacobian\n");CHKERRQ(ierr);
ierr = MatViewFromOptions(A,(PetscObject)snes,"-snes_check_jacobian_view");CHKERRQ(ierr);
ierr = PetscViewerASCIIPrintf(viewer," Hand-coded minus finite difference Jacobian\n");CHKERRQ(ierr);
ierr = MatViewFromOptions(B,(PetscObject)snes,"-snes_check_jacobian_view");CHKERRQ(ierr);
}
if (!gnorm) gnorm = 1; /* just in case */
ierr = PetscViewerASCIIPrintf(viewer," %g = ||J - Jfd||/||J|| %g = ||J - Jfd||\n",(double)(nrm/gnorm),(double)nrm);CHKERRQ(ierr);
ierr = SNESGetObjective(snes,&objective,&ctx);CHKERRQ(ierr);
if (objective) {
ierr = SNESComputeFunction(snes,x,f);CHKERRQ(ierr);
ierr = VecNorm(f,NORM_2,&fnorm);CHKERRQ(ierr);
if (complete_print) {
ierr = PetscViewerASCIIPrintf(viewer," Hand-coded Objective Function \n");CHKERRQ(ierr);
ierr = VecView(f,viewer);CHKERRQ(ierr);
}
ierr = SNESObjectiveComputeFunctionDefaultFD(snes,x,f1,NULL);CHKERRQ(ierr);
ierr = VecNorm(f1,NORM_2,&f1norm);CHKERRQ(ierr);
if (complete_print) {
ierr = PetscViewerASCIIPrintf(viewer," Finite-Difference Objective Function\n");CHKERRQ(ierr);
ierr = VecView(f1,viewer);CHKERRQ(ierr);
}
/* compare the two */
ierr = VecAXPY(f,-1.0,f1);CHKERRQ(ierr);
ierr = VecNorm(f,NORM_2,&dnorm);CHKERRQ(ierr);
if (!fnorm) fnorm = 1.;
ierr = PetscViewerASCIIPrintf(viewer," %g = Norm of objective function ratio %g = difference\n",dnorm/fnorm,dnorm);CHKERRQ(ierr);
if (complete_print) {
ierr = PetscViewerASCIIPrintf(viewer," Difference\n");CHKERRQ(ierr);
ierr = VecView(f,viewer);CHKERRQ(ierr);
}
}
ierr = PetscViewerASCIISubtractTab(viewer,((PetscObject)snes)->tablevel);CHKERRQ(ierr);
//.........这里部分代码省略.........
//.........这里部分代码省略.........
for (j=0; j<B->ilen[i]; j++) {
PetscInt data,gid1 = aj[B->i[i]+j] + 1;
ierr = PetscTableFind(gid1_lid1,gid1,&data);CHKERRQ(ierr);
if (!data) {
/* one based table */
ierr = PetscTableAdd(gid1_lid1,gid1,++ec,INSERT_VALUES);CHKERRQ(ierr);
}
}
}
/* form array of columns we need */
ierr = PetscMalloc1(ec,&garray);CHKERRQ(ierr);
ierr = PetscTableGetHeadPosition(gid1_lid1,&tpos);CHKERRQ(ierr);
while (tpos) {
ierr = PetscTableGetNext(gid1_lid1,&tpos,&gid,&lid);CHKERRQ(ierr);
gid--; lid--;
garray[lid] = gid;
}
ierr = PetscSortInt(ec,garray);CHKERRQ(ierr);
ierr = PetscTableRemoveAll(gid1_lid1);CHKERRQ(ierr);
for (i=0; i<ec; i++) {
ierr = PetscTableAdd(gid1_lid1,garray[i]+1,i+1,INSERT_VALUES);CHKERRQ(ierr);
}
/* compact out the extra columns in B */
for (i=0; i<B->mbs; i++) {
for (j=0; j<B->ilen[i]; j++) {
PetscInt gid1 = aj[B->i[i] + j] + 1;
ierr = PetscTableFind(gid1_lid1,gid1,&lid);CHKERRQ(ierr);
lid--;
aj[B->i[i]+j] = lid;
}
}
B->nbs = ec;
baij->B->cmap->n = baij->B->cmap->N = ec*mat->rmap->bs;
ierr = PetscLayoutSetUp((baij->B->cmap));CHKERRQ(ierr);
ierr = PetscTableDestroy(&gid1_lid1);CHKERRQ(ierr);
#else
/* For the first stab we make an array as long as the number of columns */
/* mark those columns that are in baij->B */
ierr = PetscCalloc1(Nbs,&indices);CHKERRQ(ierr);
for (i=0; i<B->mbs; i++) {
for (j=0; j<B->ilen[i]; j++) {
if (!indices[aj[B->i[i] + j]]) ec++;
indices[aj[B->i[i] + j]] = 1;
}
}
/* form array of columns we need */
ierr = PetscMalloc1(ec,&garray);CHKERRQ(ierr);
ec = 0;
for (i=0; i<Nbs; i++) {
if (indices[i]) {
garray[ec++] = i;
}
}
/* make indices now point into garray */
for (i=0; i<ec; i++) indices[garray[i]] = i;
/* compact out the extra columns in B */
for (i=0; i<B->mbs; i++) {
for (j=0; j<B->ilen[i]; j++) {
aj[B->i[i] + j] = indices[aj[B->i[i] + j]];
}
}
B->nbs = ec;
baij->B->cmap->n = ec*mat->rmap->bs;
ierr = PetscFree(indices);CHKERRQ(ierr);
#endif
/* create local vector that is used to scatter into */
ierr = VecCreateSeq(PETSC_COMM_SELF,ec*bs,&baij->lvec);CHKERRQ(ierr);
/* create two temporary index sets for building scatter-gather */
ierr = ISCreateBlock(PETSC_COMM_SELF,bs,ec,garray,PETSC_COPY_VALUES,&from);CHKERRQ(ierr);
ierr = PetscMalloc1(ec,&stmp);CHKERRQ(ierr);
for (i=0; i<ec; i++) stmp[i] = i;
ierr = ISCreateBlock(PETSC_COMM_SELF,bs,ec,stmp,PETSC_OWN_POINTER,&to);CHKERRQ(ierr);
/* create temporary global vector to generate scatter context */
ierr = VecCreateMPIWithArray(PetscObjectComm((PetscObject)mat),1,mat->cmap->n,mat->cmap->N,NULL,&gvec);CHKERRQ(ierr);
ierr = VecScatterCreateWithData(gvec,from,baij->lvec,to,&baij->Mvctx);CHKERRQ(ierr);
ierr = VecDestroy(&gvec);CHKERRQ(ierr);
ierr = PetscLogObjectParent((PetscObject)mat,(PetscObject)baij->Mvctx);CHKERRQ(ierr);
ierr = PetscLogObjectParent((PetscObject)mat,(PetscObject)baij->lvec);CHKERRQ(ierr);
ierr = PetscLogObjectParent((PetscObject)mat,(PetscObject)from);CHKERRQ(ierr);
ierr = PetscLogObjectParent((PetscObject)mat,(PetscObject)to);CHKERRQ(ierr);
baij->garray = garray;
ierr = PetscLogObjectMemory((PetscObject)mat,(ec+1)*sizeof(PetscInt));CHKERRQ(ierr);
ierr = ISDestroy(&from);CHKERRQ(ierr);
ierr = ISDestroy(&to);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
PetscErrorCode MatSetUpMultiply_MPISBAIJ(Mat mat)
{
Mat_MPISBAIJ *sbaij = (Mat_MPISBAIJ*)mat->data;
Mat_SeqBAIJ *B = (Mat_SeqBAIJ*)(sbaij->B->data);
PetscErrorCode ierr;
PetscInt Nbs = sbaij->Nbs,i,j,*indices,*aj = B->j,ec = 0,*garray,*sgarray;
PetscInt bs = mat->rmap->bs,*stmp,mbs=sbaij->mbs, vec_size,nt;
IS from,to;
Vec gvec;
PetscMPIInt rank =sbaij->rank,lsize,size=sbaij->size;
PetscInt *owners=sbaij->rangebs,*ec_owner,k;
const PetscInt *sowners;
PetscScalar *ptr;
PetscFunctionBegin;
ierr = VecScatterDestroy(&sbaij->sMvctx);CHKERRQ(ierr);
/* For the first stab we make an array as long as the number of columns */
/* mark those columns that are in sbaij->B */
ierr = PetscCalloc1(Nbs,&indices);CHKERRQ(ierr);
for (i=0; i<mbs; i++) {
for (j=0; j<B->ilen[i]; j++) {
if (!indices[aj[B->i[i] + j]]) ec++;
indices[aj[B->i[i] + j]] = 1;
}
}
/* form arrays of columns we need */
ierr = PetscMalloc1(ec,&garray);CHKERRQ(ierr);
ierr = PetscMalloc2(2*ec,&sgarray,ec,&ec_owner);CHKERRQ(ierr);
ec = 0;
for (j=0; j<size; j++) {
for (i=owners[j]; i<owners[j+1]; i++) {
if (indices[i]) {
garray[ec] = i;
ec_owner[ec] = j;
ec++;
}
}
}
/* make indices now point into garray */
for (i=0; i<ec; i++) indices[garray[i]] = i;
/* compact out the extra columns in B */
for (i=0; i<mbs; i++) {
for (j=0; j<B->ilen[i]; j++) aj[B->i[i] + j] = indices[aj[B->i[i] + j]];
}
B->nbs = ec;
sbaij->B->cmap->n = sbaij->B->cmap->N = ec*mat->rmap->bs;
ierr = PetscLayoutSetUp((sbaij->B->cmap));CHKERRQ(ierr);
ierr = PetscFree(indices);CHKERRQ(ierr);
/* create local vector that is used to scatter into */
ierr = VecCreateSeq(PETSC_COMM_SELF,ec*bs,&sbaij->lvec);CHKERRQ(ierr);
/* create two temporary index sets for building scatter-gather */
ierr = PetscMalloc1(2*ec,&stmp);CHKERRQ(ierr);
ierr = ISCreateBlock(PETSC_COMM_SELF,bs,ec,garray,PETSC_COPY_VALUES,&from);CHKERRQ(ierr);
for (i=0; i<ec; i++) stmp[i] = i;
ierr = ISCreateBlock(PETSC_COMM_SELF,bs,ec,stmp,PETSC_COPY_VALUES,&to);CHKERRQ(ierr);
/* generate the scatter context
-- Mvctx and lvec are not used by MatMult_MPISBAIJ(), but usefule for some applications */
ierr = VecCreateMPIWithArray(PetscObjectComm((PetscObject)mat),1,mat->cmap->n,mat->cmap->N,NULL,&gvec);CHKERRQ(ierr);
ierr = VecScatterCreateWithData(gvec,from,sbaij->lvec,to,&sbaij->Mvctx);CHKERRQ(ierr);
ierr = VecDestroy(&gvec);CHKERRQ(ierr);
sbaij->garray = garray;
ierr = PetscLogObjectParent((PetscObject)mat,(PetscObject)sbaij->Mvctx);CHKERRQ(ierr);
ierr = PetscLogObjectParent((PetscObject)mat,(PetscObject)sbaij->lvec);CHKERRQ(ierr);
ierr = PetscLogObjectParent((PetscObject)mat,(PetscObject)from);CHKERRQ(ierr);
ierr = PetscLogObjectParent((PetscObject)mat,(PetscObject)to);CHKERRQ(ierr);
ierr = ISDestroy(&from);CHKERRQ(ierr);
ierr = ISDestroy(&to);CHKERRQ(ierr);
/* create parallel vector that is used by SBAIJ matrix to scatter from/into */
lsize = (mbs + ec)*bs;
ierr = VecCreateMPI(PetscObjectComm((PetscObject)mat),lsize,PETSC_DETERMINE,&sbaij->slvec0);CHKERRQ(ierr);
ierr = VecDuplicate(sbaij->slvec0,&sbaij->slvec1);CHKERRQ(ierr);
ierr = VecGetSize(sbaij->slvec0,&vec_size);CHKERRQ(ierr);
ierr = VecGetOwnershipRanges(sbaij->slvec0,&sowners);CHKERRQ(ierr);
/* x index in the IS sfrom */
for (i=0; i<ec; i++) {
j = ec_owner[i];
sgarray[i] = garray[i] + (sowners[j]/bs - owners[j]);
}
/* b index in the IS sfrom */
k = sowners[rank]/bs + mbs;
for (i=ec,j=0; i< 2*ec; i++,j++) sgarray[i] = k + j;
ierr = ISCreateBlock(PETSC_COMM_SELF,bs,2*ec,sgarray,PETSC_COPY_VALUES,&from);CHKERRQ(ierr);
/* x index in the IS sto */
//.........这里部分代码省略.........
/*@C
KSPView - Prints the KSP data structure.
Collective on KSP
Input Parameters:
+ ksp - the Krylov space context
- viewer - visualization context
Options Database Keys:
. -ksp_view - print the ksp data structure at the end of a KSPSolve call
Note:
The available visualization contexts include
+ PETSC_VIEWER_STDOUT_SELF - standard output (default)
- PETSC_VIEWER_STDOUT_WORLD - synchronized standard
output where only the first processor opens
the file. All other processors send their
data to the first processor to print.
The user can open an alternative visualization context with
PetscViewerASCIIOpen() - output to a specified file.
Level: beginner
.keywords: KSP, view
.seealso: PCView(), PetscViewerASCIIOpen()
@*/
PetscErrorCode KSPView(KSP ksp,PetscViewer viewer)
{
PetscErrorCode ierr;
PetscBool iascii,isbinary,isdraw;
#if defined(PETSC_HAVE_AMS)
PetscBool isams;
#endif
PetscFunctionBegin;
PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
if (!viewer) viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ksp));
PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2);
PetscCheckSameComm(ksp,1,viewer,2);
ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr);
ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr);
#if defined(PETSC_HAVE_AMS)
ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERAMS,&isams);CHKERRQ(ierr);
#endif
if (iascii) {
ierr = PetscObjectPrintClassNamePrefixType((PetscObject)ksp,viewer);CHKERRQ(ierr);
if (ksp->ops->view) {
ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
ierr = (*ksp->ops->view)(ksp,viewer);CHKERRQ(ierr);
ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
}
if (ksp->guess_zero) {
ierr = PetscViewerASCIIPrintf(viewer," maximum iterations=%D, initial guess is zero\n",ksp->max_it);CHKERRQ(ierr);
} else {
ierr = PetscViewerASCIIPrintf(viewer," maximum iterations=%D\n", ksp->max_it);CHKERRQ(ierr);
}
if (ksp->guess_knoll) {ierr = PetscViewerASCIIPrintf(viewer," using preconditioner applied to right hand side for initial guess\n");CHKERRQ(ierr);}
ierr = PetscViewerASCIIPrintf(viewer," tolerances: relative=%G, absolute=%G, divergence=%G\n",ksp->rtol,ksp->abstol,ksp->divtol);CHKERRQ(ierr);
if (ksp->pc_side == PC_RIGHT) {
ierr = PetscViewerASCIIPrintf(viewer," right preconditioning\n");CHKERRQ(ierr);
} else if (ksp->pc_side == PC_SYMMETRIC) {
ierr = PetscViewerASCIIPrintf(viewer," symmetric preconditioning\n");CHKERRQ(ierr);
} else {
ierr = PetscViewerASCIIPrintf(viewer," left preconditioning\n");CHKERRQ(ierr);
}
if (ksp->guess) {ierr = PetscViewerASCIIPrintf(viewer," using Fischers initial guess method %D with size %D\n",ksp->guess->method,ksp->guess->maxl);CHKERRQ(ierr);}
if (ksp->dscale) {ierr = PetscViewerASCIIPrintf(viewer," diagonally scaled system\n");CHKERRQ(ierr);}
if (ksp->nullsp) {ierr = PetscViewerASCIIPrintf(viewer," has attached null space\n");CHKERRQ(ierr);}
if (!ksp->guess_zero) {ierr = PetscViewerASCIIPrintf(viewer," using nonzero initial guess\n");CHKERRQ(ierr);}
ierr = PetscViewerASCIIPrintf(viewer," using %s norm type for convergence test\n",KSPNormTypes[ksp->normtype]);CHKERRQ(ierr);
} else if (isbinary) {
PetscInt classid = KSP_FILE_CLASSID;
MPI_Comm comm;
PetscMPIInt rank;
char type[256];
ierr = PetscObjectGetComm((PetscObject)ksp,&comm);CHKERRQ(ierr);
ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
if (!rank) {
ierr = PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr);
ierr = PetscStrncpy(type,((PetscObject)ksp)->type_name,256);CHKERRQ(ierr);
ierr = PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR,PETSC_FALSE);CHKERRQ(ierr);
}
if (ksp->ops->view) {
ierr = (*ksp->ops->view)(ksp,viewer);CHKERRQ(ierr);
}
} else if (isdraw) {
PetscDraw draw;
char str[36];
PetscReal x,y,bottom,h;
PetscBool flg;
ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr);
ierr = PetscDrawGetCurrentPoint(draw,&x,&y);CHKERRQ(ierr);
ierr = PetscObjectTypeCompare((PetscObject)ksp,KSPPREONLY,&flg);CHKERRQ(ierr);
//.........这里部分代码省略.........
客服电话
APP下载
官方微信
请发表评论