Actual source code: ex83.c

  1: static char help[] = "Test the Fischer-1 initial guess routine with VECNEST.\n\n";

  3: #include <petscksp.h>

  5: int main(int argc, char **args)
  6: {
  7:   /* Test that exceeding the number of stored vectors works correctly - this used to not work with VecNest */
  8:   PetscInt    triangle_size = 10;
  9:   Mat         A, A_nest;
 10:   KSP         ksp;
 11:   KSPGuess    guess;
 12:   Vec         sol, rhs, sol_nest, rhs_nest;
 13:   PetscInt    i, j, indices[] = {0, 1, 2, 3, 4};
 14:   PetscScalar values[] = {1.0, 2.0, 3.0, 4.0, 5.0};

 16:   PetscCall(PetscInitialize(&argc, &args, NULL, help));

 18:   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, triangle_size, triangle_size, NULL, &A));
 19:   for (i = 0; i < triangle_size; ++i) PetscCall(MatSetValue(A, i, i, 1.0, INSERT_VALUES));
 20:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 21:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
 22:   PetscCall(MatCreateNest(PETSC_COMM_SELF, 1, NULL, 1, NULL, &A, &A_nest));
 23:   PetscCall(MatNestSetVecType(A_nest, VECNEST));

 25:   PetscCall(KSPCreate(PETSC_COMM_SELF, &ksp));
 26:   PetscCall(KSPSetOperators(ksp, A_nest, A_nest));
 27:   PetscCall(KSPSetFromOptions(ksp));
 28:   PetscCall(KSPGetGuess(ksp, &guess));
 29:   PetscCall(KSPGuessSetUp(guess));

 31:   for (i = 0; i < 5; ++i) {
 32:     PetscCall(VecCreateSeq(PETSC_COMM_SELF, triangle_size, &sol));
 33:     PetscCall(VecCreateSeq(PETSC_COMM_SELF, triangle_size, &rhs));
 34:     for (j = 0; j < i; ++j) {
 35:       PetscCall(VecSetValue(sol, j, (PetscScalar)j, INSERT_VALUES));
 36:       PetscCall(VecSetValue(rhs, j, (PetscScalar)j, INSERT_VALUES));
 37:     }
 38:     PetscCall(VecAssemblyBegin(sol));
 39:     PetscCall(VecAssemblyBegin(rhs));
 40:     PetscCall(VecAssemblyEnd(sol));
 41:     PetscCall(VecAssemblyEnd(rhs));

 43:     PetscCall(VecCreateNest(PETSC_COMM_SELF, 1, NULL, &sol, &sol_nest));
 44:     PetscCall(VecCreateNest(PETSC_COMM_SELF, 1, NULL, &rhs, &rhs_nest));

 46:     PetscCall(KSPGuessUpdate(guess, rhs_nest, sol_nest));

 48:     PetscCall(VecDestroy(&rhs_nest));
 49:     PetscCall(VecDestroy(&sol_nest));
 50:     PetscCall(VecDestroy(&rhs));
 51:     PetscCall(VecDestroy(&sol));
 52:   }

 54:   PetscCall(VecCreateSeq(PETSC_COMM_SELF, triangle_size, &sol));
 55:   PetscCall(VecCreateSeq(PETSC_COMM_SELF, triangle_size, &rhs));
 56:   PetscCall(VecSetValues(rhs, 5, indices, values, INSERT_VALUES));
 57:   PetscCall(VecAssemblyBegin(rhs));
 58:   PetscCall(VecAssemblyEnd(rhs));

 60:   PetscCall(VecCreateNest(PETSC_COMM_SELF, 1, NULL, &sol, &sol_nest));
 61:   PetscCall(VecCreateNest(PETSC_COMM_SELF, 1, NULL, &rhs, &rhs_nest));

 63:   PetscCall(KSPGuessFormGuess(guess, rhs_nest, sol_nest));
 64:   PetscCall(VecView(sol_nest, PETSC_VIEWER_STDOUT_SELF));

 66:   PetscCall(VecDestroy(&rhs_nest));
 67:   PetscCall(VecDestroy(&sol_nest));
 68:   PetscCall(VecDestroy(&rhs));
 69:   PetscCall(VecDestroy(&sol));

 71:   PetscCall(KSPDestroy(&ksp));

 73:   PetscCall(MatDestroy(&A_nest));
 74:   PetscCall(MatDestroy(&A));

 76:   PetscCall(PetscFinalize());
 77:   return 0;
 78: }

 80: /*TEST

 82:    test:
 83:       args: -ksp_guess_type fischer -ksp_guess_fischer_model 1,3 -ksp_guess_fischer_monitor

 85: TEST*/