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;

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

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

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

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

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

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

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

 56:       KSPDestroy(&ksp);
 57:     }

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

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

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

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

 81:     MatDestroy(&A);
 82:   }

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

 89:     MatCreateSeqDense(PETSC_COMM_SELF, triangle_size, triangle_size, NULL, &A);
 90:     for (i = 0; i < triangle_size; ++i) MatSetValue(A, i, i, 1.0, INSERT_VALUES);
 91:     MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 92:     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:       KSPCreate(PETSC_COMM_SELF, &ksp);
102:       KSPSetOperators(ksp, A, A);
103:       KSPSetFromOptions(ksp);
104:       KSPGetGuess(ksp, &guess);
105:       KSPGuessSetUp(guess);

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

119:         KSPGuessUpdate(guess, rhs, sol);

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

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

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

134:       VecDestroy(&rhs);
135:       VecDestroy(&sol);
136:       KSPDestroy(&ksp);
137:     }
138:     MatDestroy(&A);
139:   }
140:   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*/