Actual source code: ex80.c

  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:   PetscInt i;
 10:   {
 11:     Mat         A;
 12:     PetscInt    indices[SIZE] = {0, 1, 2};
 13:     PetscScalar values[SIZE]  = {1.0, 1.0, 1.0};
 14:     Vec         sol, rhs, newsol, newrhs;

 16:     PetscFunctionBeginUser;
 17:     PetscCall(PetscInitialize(&argc, &args, (char *)0, help));

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

144: /* 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 */

146: /*TEST

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

151: TEST*/