- petsc-3.18.1/src/ksp/ksp/tutorials$ ./ex1
- KSP Object: 1 MPI process
- type: gmres
- restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement
- happy breakdown tolerance 1e-30
- maximum iterations=10000, initial guess is zero
- tolerances: relative=1e-05, absolute=1e-50, divergence=10000.
- left preconditioning
- using PRECONDITIONED norm type for convergence test
- PC Object: 1 MPI process
- type: jacobi
- type DIAGONAL
- linear system matrix = precond matrix:
- Mat Object: 1 MPI process
- type: seqaij
- rows=10, cols=10
- total: nonzeros=28, allocated nonzeros=50
- total number of mallocs used during MatSetValues calls=0
- not using I-node routines
- Norm of error 2.41202e-15, Iterations 5
使用的是GMRES。
1) gdb ex1
130 PetscCall(KSPSolve(ksp, b, x));
2)
at /thfs1/home/***/petsc-3.18.1/src/ksp/ksp/interface/itfunc.c:1066
不太1066 PetscFunctionBegin;
- (gdb) p KSP_CLASSID
- $1 = 1211223
- (gdb) p VEC_CLASSID
- $2 = 1211214
不太明白这个ID是啥意思。
实际调用的是KSPSolve_Private()
- (gdb) s
- KSPSolve_Private (ksp=0x5c62a0, b=0x4d3060, x=0x4cdc80)
- at /thfs1/home/monkeycode/tianya/software/petsc-3.18.1/src/ksp/ksp/interface/itfunc.c:791
- 791 PetscBool flg = PETSC_FALSE, inXisinB = PETSC_FALSE, guess_zero;
3)CG可以看ex72.c
-mat_type mpisbaij -ksp_type cg -pc_type eisenstat -ksp_monitor_short -ksp_converged_reason
test:
-mat_type mpisbaij -ksp_type cg -pc_type eisenstat -ksp_monitor_short -ksp_converged_reason
-ksp_cg_single_reduction
-ksp_ksp_type cg -ksp_pc_type bjacobi -ksp_pc_bjacobi_blocks 1
4)
- 125 if (!ksp->guess_zero) {
- 126 PetscCall(KSP_MatMult(ksp, Amat, X, R)); /* r <- b - Ax * /
- 127 PetscCall(VecAYPX(R, -1.0, B));
- 128 } else {
- 129 PetscCall(VecCopy(B, R)); /* r <- b (x is 0) */
- 130 }
-
-
- switch (ksp->normtype) {
- 135 case KSP_NORM_PRECONDITIONED:
- 136 PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
- 137 PetscCall(VecNorm(Z, NORM_2, &dp)); /* dp <- z'*z = e'*A'*B'*B*A*e */
- 138 KSPCheckNorm(ksp, dp);
- 139 break;
- 140 case KSP_NORM_UNPRECONDITIONED:
- 141 PetscCall(VecNorm(R, NORM_2, &dp)); /* dp <- r'*r = e'*A'*A*e */
- 142 KSPCheckNorm(ksp, dp);
- 143 break;
- 144 case KSP_NORM_NATURAL:
- 145 PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
- 146 PetscCall(VecXDot(Z, R, &beta)); /* beta <- z'*r */
- 147 KSPCheckDot(ksp, beta);
- 148 dp = PetscSqrtReal(PetscAbsScalar(beta)); /* dp <- r'*z = r'*B*r = e'*A'*B*A*e */
-
-
- 78 /*
- 79 A macro used in the following KSPSolve_CG and KSPSolve_CG_SingleReduction routine s
- 80 */
- 81 #define VecXDot(x, y, a) (((cg->type) == (KSP_CG_HERMITIAN)) ? VecDot(x, y, a) : VecTD ot(x, y, a))
-
-
- /petsc-3.18.1/src/vec$ grep VecAYPX -rn *
- f90-mod/ftn-auto-interfaces/petscvec.h90:188: subroutine VecAYPX(a,b,c,z)
- f90-mod/ftn-auto-interfaces/petscvec.h90:194: end subroutine VecAYPX
- vec/interface/rvector.c:571:$ VecAYPX(y,beta,x) y = x + beta y
- vec/interface/rvector.c:577:.seealso: `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
- vec/interface/rvector.c:603: VecAYPX - Computes y = x + beta y.
- vec/interface/rvector.c:622:PetscErrorCode VecAYPX(Vec y, PetscScalar beta, Vec x)
- vec/interface/rvector.c:665:.seealso: `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`
-
5)继续gdb ex1
- (gdb)
- 899 PetscUseTypeMethod(ksp, solve);
- (gdb) s
- KSPSolve_GMRES (ksp=0x5c69b0)
- at /thfs1/home/monkeycode/tianya/software/petsc-3.18.1/src/ksp/ksp/impls/gmres/gmres.c:212
- 212 KSP_GMRES *gmres = (KSP_GMRES *)ksp->data;
函数作为参数
有些代码和判断比较长,但是执行次数少,就并不
也可以这样运行:
./ex1 -m 400 -n 400 -pc_type eisenstat -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always
应该核心是:
- 899 PetscUseTypeMethod(ksp, solve);
- (gdb)
- 0 KSP Residual norm 1.52753
- 1 KSP Residual norm 0.437136
- 2 KSP Residual norm 0.22205
- 3 KSP Residual norm 0.137425
- 4 KSP Residual norm 0.0952384
- 5 KSP Residual norm 0.0709072
- 6 KSP Residual norm 0.0554212
- 7 KSP Residual norm 0.0448558
- 8 KSP Residual norm 0.0372695
- 9 KSP Residual norm 0.0316055
- 10 KSP Residual norm 0.0272444
6)
- (gdb)
- 228 PetscCall(KSPGMRESCycle(&its, ksp));
- (gdb)
- 0 KSP Residual norm 1.52753
- 1 KSP Residual norm 0.437136
- 2 KSP Residual norm 0.22205
- 3 KSP Residual norm 0.137425
- 4 KSP Residual norm 0.0952384
- 5 KSP Residual norm 0.0709072
- 6 KSP Residual norm 0.0554212
- 7 KSP Residual norm 0.0448558
- 8 KSP Residual norm 0.0372695
-
- KSPGMRESCycle (itcount=0xffffffff0560, ksp=0x5d0bd0)
- at /thfs1/home/monkeycode/tianya/software/petsc-3.18.1/src/ksp/ksp/impls/gmres/gmres.c:103
这个地方是核心求解?
/* update hessenberg matrix and do Gram-Schmidt */
PetscCall((*gmres->orthog)(ksp, it));
176 PetscCall((*ksp->converged)(ksp, ksp->its, res, &ksp->reason, ksp->cnvP));