PETSC调试

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));