1: static char help[] = "Tests solving linear system on 0 by 0 matrix, and KSPLSQR convergence test handling.\n\n";

  3: #include <petscksp.h>

  5: static PetscErrorCode GetConvergenceTestName(KSPConvergenceTestFn *converged, char name[], size_t n)
  6: {
  7:   PetscFunctionBegin;
  8:   if (converged == KSPConvergedDefault) {
  9:     PetscCall(PetscStrncpy(name, "default", n));
 10:   } else if (converged == KSPConvergedSkip) {
 11:     PetscCall(PetscStrncpy(name, "skip", n));
 12:   } else if (converged == KSPLSQRConvergedDefault) {
 13:     PetscCall(PetscStrncpy(name, "lsqr", n));
 14:   } else {
 15:     PetscCall(PetscStrncpy(name, "other", n));
 16:   }
 17:   PetscFunctionReturn(PETSC_SUCCESS);
 18: }

 20: int main(int argc, char **args)
 21: {
 22:   Mat       C;
 23:   PetscInt  N = 0;
 24:   Vec       u, b, x;
 25:   KSP       ksp;
 26:   PetscReal norm;
 27:   PetscBool flg = PETSC_FALSE;

 29:   PetscFunctionBeginUser;
 30:   PetscCall(PetscInitialize(&argc, &args, NULL, help));

 32:   /* create stiffness matrix */
 33:   PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
 34:   PetscCall(MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, N, N));
 35:   PetscCall(MatSetFromOptions(C));
 36:   PetscCall(MatSetUp(C));
 37:   PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
 38:   PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));

 40:   /* create right-hand side and solution */
 41:   PetscCall(VecCreate(PETSC_COMM_WORLD, &u));
 42:   PetscCall(VecSetSizes(u, PETSC_DECIDE, N));
 43:   PetscCall(VecSetFromOptions(u));
 44:   PetscCall(VecDuplicate(u, &b));
 45:   PetscCall(VecDuplicate(u, &x));

 47:   PetscCall(VecAssemblyBegin(b));
 48:   PetscCall(VecAssemblyEnd(b));

 50:   /* solve linear system */
 51:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
 52:   PetscCall(KSPSetOperators(ksp, C, C));
 53:   PetscCall(KSPSetFromOptions(ksp));
 54:   PetscCall(KSPSolve(ksp, b, u));

 56:   /* test proper handling of convergence test by KSPLSQR */
 57:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_lsqr", &flg, NULL));
 58:   if (flg) {
 59:     char                 *type;
 60:     char                  convtestname[16];
 61:     PetscBool             islsqr;
 62:     KSPConvergenceTestFn *converged, *converged1;
 63:     PetscCtxDestroyFn    *destroy, *destroy1;
 64:     void                 *ctx, *ctx1;

 66:     {
 67:       const char *typeP;
 68:       PetscCall(KSPGetType(ksp, &typeP));
 69:       PetscCall(PetscStrallocpy(typeP, &type));
 70:     }
 71:     PetscCall(PetscStrcmp(type, KSPLSQR, &islsqr));
 72:     PetscCall(KSPGetConvergenceTest(ksp, &converged, &ctx, &destroy));
 73:     PetscCall(GetConvergenceTestName(converged, convtestname, 16));
 74:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "convergence test: %s\n", convtestname));
 75:     PetscCall(KSPSetType(ksp, KSPLSQR));
 76:     PetscCall(KSPGetConvergenceTest(ksp, &converged1, &ctx1, &destroy1));
 77:     PetscCheck(converged1 == KSPLSQRConvergedDefault, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "convergence test should be KSPLSQRConvergedDefault");
 78:     PetscCheck(destroy1 == KSPConvergedDefaultDestroy, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "convergence test destroy function should be KSPConvergedDefaultDestroy");
 79:     if (islsqr) {
 80:       PetscCheck(converged1 == converged, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "convergence test should be kept");
 81:       PetscCheck(destroy1 == destroy, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "convergence test destroy function should be kept");
 82:       PetscCheck(ctx1 == ctx, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "convergence test context should be kept");
 83:     }
 84:     PetscCall(GetConvergenceTestName(converged1, convtestname, 16));
 85:     PetscCall(KSPViewFromOptions(ksp, NULL, "-ksp1_view"));
 86:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "convergence test: %s\n", convtestname));
 87:     PetscCall(KSPSetType(ksp, type));
 88:     PetscCall(KSPGetConvergenceTest(ksp, &converged1, &ctx1, &destroy1));
 89:     PetscCheck(converged1 == converged, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "convergence test not reverted properly");
 90:     PetscCheck(destroy1 == destroy, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "convergence test destroy function not reverted properly");
 91:     PetscCheck(ctx1 == ctx, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "convergence test context not reverted properly");
 92:     PetscCall(GetConvergenceTestName(converged1, convtestname, 16));
 93:     PetscCall(KSPViewFromOptions(ksp, NULL, "-ksp2_view"));
 94:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "convergence test: %s\n", convtestname));
 95:     PetscCall(PetscFree(type));
 96:   }

 98:   PetscCall(MatMult(C, u, x));
 99:   PetscCall(VecAXPY(x, -1.0, b));
100:   PetscCall(VecNorm(x, NORM_2, &norm));

102:   PetscCall(KSPDestroy(&ksp));
103:   PetscCall(VecDestroy(&u));
104:   PetscCall(VecDestroy(&x));
105:   PetscCall(VecDestroy(&b));
106:   PetscCall(MatDestroy(&C));
107:   PetscCall(PetscFinalize());
108:   return 0;
109: }

111: /*TEST

113:     test:
114:       args: -pc_type jacobi -ksp_monitor -ksp_orthogonalization_cgs_refinement_type refine_always

116:     test:
117:       suffix: 2
118:       nsize: 2
119:       args: -pc_type jacobi -ksp_monitor -ksp_orthogonalization_cgs_refinement_type refine_always

121:     test:
122:       suffix: 3
123:       args: -pc_type sor -pc_sor_symmetric -ksp_monitor -ksp_orthogonalization_cgs_refinement_type refine_always

125:     test:
126:       suffix: 5
127:       args: -pc_type eisenstat -ksp_monitor -ksp_orthogonalization_cgs_refinement_type refine_always

129:     testset:
130:       args: -test_lsqr -ksp{,1,2}_view -pc_type jacobi
131:       filter: grep -E "(^  type:|preconditioning|norm type|convergence test:)"
132:       test:
133:         suffix: lsqr_0
134:         args: -ksp_convergence_test {{default skip}separate output}
135:       test:
136:         suffix: lsqr_1
137:         args: -ksp_type cg -ksp_convergence_test {{default skip}separate output}
138:       test:
139:         suffix: lsqr_2
140:         args: -ksp_type lsqr

142: TEST*/