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

  3: #include <petscksp.h>

  5: #define SIZE 3

  7: int main(int argc, char **args)
  8: {
  9:   {
 10:     Mat         A;
 11:     PetscInt    indices[SIZE] = {0, 1, 2};
 12:     PetscScalar values[SIZE]  = {1.0, 1.0, 1.0};
 13:     Vec         sol, rhs, newsol, newrhs;

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

 18:     /* common data structures */
 19:     PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, SIZE, SIZE, NULL, &A));
 20:     for (PetscInt i = 0; i < SIZE; ++i) PetscCall(MatSetValue(A, i, i, 1.0, INSERT_VALUES));
 21:     PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 22:     PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

 24:     PetscCall(VecCreateSeq(PETSC_COMM_SELF, SIZE, &sol));
 25:     PetscCall(VecDuplicate(sol, &rhs));
 26:     PetscCall(VecDuplicate(sol, &newrhs));
 27:     PetscCall(VecDuplicate(sol, &newsol));

 29:     PetscCall(VecSetValues(sol, SIZE, indices, values, INSERT_VALUES));
 30:     PetscCall(VecSetValues(rhs, SIZE - 1, indices, values, INSERT_VALUES));
 31:     PetscCall(VecSetValues(newrhs, SIZE - 2, indices, values, INSERT_VALUES));
 32:     PetscCall(VecAssemblyBegin(sol));
 33:     PetscCall(VecAssemblyBegin(rhs));
 34:     PetscCall(VecAssemblyBegin(newrhs));
 35:     PetscCall(VecAssemblyEnd(sol));
 36:     PetscCall(VecAssemblyEnd(rhs));
 37:     PetscCall(VecAssemblyEnd(newrhs));

 39:     /* Test one vector */
 40:     {
 41:       KSP      ksp;
 42:       KSPGuess guess;

 44:       PetscCall(KSPCreate(PETSC_COMM_SELF, &ksp));
 45:       PetscCall(KSPSetOperators(ksp, A, A));
 46:       PetscCall(KSPSetFromOptions(ksp));
 47:       PetscCall(KSPGetGuess(ksp, &guess));
 48:       /* we aren't calling through the KSP so we call this ourselves */
 49:       PetscCall(KSPGuessSetUp(guess));

 51:       PetscCall(KSPGuessUpdate(guess, rhs, sol));
 52:       PetscCall(KSPGuessFormGuess(guess, newrhs, newsol));
 53:       PetscCall(VecView(newsol, PETSC_VIEWER_STDOUT_SELF));

 55:       PetscCall(KSPDestroy(&ksp));
 56:     }

 58:     /* Test a singular projection matrix */
 59:     {
 60:       KSP      ksp;
 61:       KSPGuess guess;

 63:       PetscCall(KSPCreate(PETSC_COMM_SELF, &ksp));
 64:       PetscCall(KSPSetOperators(ksp, A, A));
 65:       PetscCall(KSPSetFromOptions(ksp));
 66:       PetscCall(KSPGetGuess(ksp, &guess));
 67:       PetscCall(KSPGuessSetUp(guess));

 69:       for (PetscInt i = 0; i < 15; ++i) PetscCall(KSPGuessUpdate(guess, rhs, sol));
 70:       PetscCall(KSPGuessFormGuess(guess, newrhs, newsol));
 71:       PetscCall(VecView(newsol, PETSC_VIEWER_STDOUT_SELF));

 73:       PetscCall(KSPDestroy(&ksp));
 74:     }
 75:     PetscCall(VecDestroy(&newsol));
 76:     PetscCall(VecDestroy(&newrhs));
 77:     PetscCall(VecDestroy(&rhs));
 78:     PetscCall(VecDestroy(&sol));

 80:     PetscCall(MatDestroy(&A));
 81:   }

 83:   /* Test something triangular */
 84:   {
 85:     PetscInt triangle_size = 10;
 86:     Mat      A;

 88:     PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, triangle_size, triangle_size, NULL, &A));
 89:     for (PetscInt i = 0; i < triangle_size; ++i) PetscCall(MatSetValue(A, i, i, 1.0, INSERT_VALUES));
 90:     PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 91:     PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

 93:     {
 94:       KSP         ksp;
 95:       KSPGuess    guess;
 96:       Vec         sol, rhs;
 97:       PetscInt    j, indices[] = {0, 1, 2, 3, 4};
 98:       PetscScalar values[] = {1.0, 2.0, 3.0, 4.0, 5.0};

100:       PetscCall(KSPCreate(PETSC_COMM_SELF, &ksp));
101:       PetscCall(KSPSetOperators(ksp, A, A));
102:       PetscCall(KSPSetFromOptions(ksp));
103:       PetscCall(KSPGetGuess(ksp, &guess));
104:       PetscCall(KSPGuessSetUp(guess));

106:       for (PetscInt i = 0; i < 5; ++i) {
107:         PetscCall(VecCreateSeq(PETSC_COMM_SELF, triangle_size, &sol));
108:         PetscCall(VecCreateSeq(PETSC_COMM_SELF, triangle_size, &rhs));
109:         for (j = 0; j < i; ++j) {
110:           PetscCall(VecSetValue(sol, j, (PetscScalar)j, INSERT_VALUES));
111:           PetscCall(VecSetValue(rhs, j, (PetscScalar)j, INSERT_VALUES));
112:         }
113:         PetscCall(VecAssemblyBegin(sol));
114:         PetscCall(VecAssemblyBegin(rhs));
115:         PetscCall(VecAssemblyEnd(sol));
116:         PetscCall(VecAssemblyEnd(rhs));

118:         PetscCall(KSPGuessUpdate(guess, rhs, sol));

120:         PetscCall(VecDestroy(&rhs));
121:         PetscCall(VecDestroy(&sol));
122:       }

124:       PetscCall(VecCreateSeq(PETSC_COMM_SELF, triangle_size, &sol));
125:       PetscCall(VecCreateSeq(PETSC_COMM_SELF, triangle_size, &rhs));
126:       PetscCall(VecSetValues(rhs, 5, indices, values, INSERT_VALUES));
127:       PetscCall(VecAssemblyBegin(sol));
128:       PetscCall(VecAssemblyEnd(sol));

130:       PetscCall(KSPGuessFormGuess(guess, rhs, sol));
131:       PetscCall(VecView(sol, PETSC_VIEWER_STDOUT_SELF));

133:       PetscCall(VecDestroy(&rhs));
134:       PetscCall(VecDestroy(&sol));
135:       PetscCall(KSPDestroy(&ksp));
136:     }
137:     PetscCall(MatDestroy(&A));
138:   }
139:   PetscCall(PetscFinalize());
140:   return 0;
141: }

143: /* The relative tolerance here is strict enough to get rid of all the noise in both single and double precision: values as low as 5e-7 also work */

145: /*TEST

147:    test:
148:       args: -ksp_guess_type fischer -ksp_guess_fischer_model 3,10 -ksp_guess_fischer_monitor -ksp_guess_fischer_tol 1e-6

150: TEST*/