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:     PetscInitialize(&argc,&args,(char*)0,help);

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

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

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

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

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

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

 57:       KSPDestroy(&ksp);
 58:     }

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

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

 71:       for (i = 0; i < 15; ++i) {
 72:         KSPGuessUpdate(guess,rhs,sol);
 73:       }
 74:       KSPGuessFormGuess(guess,newrhs,newsol);
 75:       VecView(newsol,PETSC_VIEWER_STDOUT_SELF);

 77:       KSPDestroy(&ksp);
 78:     }
 79:     VecDestroy(&newsol);
 80:     VecDestroy(&newrhs);
 81:     VecDestroy(&rhs);
 82:     VecDestroy(&sol);

 84:     MatDestroy(&A);
 85:   }

 87:   /* Test something triangular */
 88:   {
 89:     PetscInt triangle_size = 10;
 90:     Mat      A;

 92:     MatCreateSeqDense(PETSC_COMM_SELF,triangle_size,triangle_size,NULL,&A);
 93:     for (i = 0; i < triangle_size; ++i) {
 94:       MatSetValue(A,i,i,1.0,INSERT_VALUES);
 95:     }
 96:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 97:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

 99:     {
100:       KSP         ksp;
101:       KSPGuess    guess;
102:       Vec         sol,rhs;
103:       PetscInt    j,indices[] = {0,1,2,3,4};
104:       PetscScalar values[] = {1.0,2.0,3.0,4.0,5.0};

106:       KSPCreate(PETSC_COMM_SELF,&ksp);
107:       KSPSetOperators(ksp,A,A);
108:       KSPSetFromOptions(ksp);
109:       KSPGetGuess(ksp,&guess);
110:       KSPGuessSetUp(guess);

112:       for (i = 0; i < 5; ++i) {
113:         VecCreateSeq(PETSC_COMM_SELF,triangle_size,&sol);
114:         VecCreateSeq(PETSC_COMM_SELF,triangle_size,&rhs);
115:         for (j = 0; j < i; ++j) {
116:           VecSetValue(sol,j,(PetscScalar)j,INSERT_VALUES);
117:           VecSetValue(rhs,j,(PetscScalar)j,INSERT_VALUES);
118:         }
119:         VecAssemblyBegin(sol);
120:         VecAssemblyBegin(rhs);
121:         VecAssemblyEnd(sol);
122:         VecAssemblyEnd(rhs);

124:         KSPGuessUpdate(guess,rhs,sol);

126:         VecDestroy(&rhs);
127:         VecDestroy(&sol);
128:       }

130:       VecCreateSeq(PETSC_COMM_SELF,triangle_size,&sol);
131:       VecCreateSeq(PETSC_COMM_SELF,triangle_size,&rhs);
132:       VecSetValues(rhs,5,indices,values,INSERT_VALUES);
133:       VecAssemblyBegin(sol);
134:       VecAssemblyEnd(sol);

136:       KSPGuessFormGuess(guess,rhs,sol);
137:       VecView(sol,PETSC_VIEWER_STDOUT_SELF);

139:       VecDestroy(&rhs);
140:       VecDestroy(&sol);
141:       KSPDestroy(&ksp);
142:     }
143:     MatDestroy(&A);
144:   }
145:   PetscFinalize();
146:   return 0;
147: }

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

151: /*TEST

153:    test:
154:       args: -ksp_guess_type fischer -ksp_guess_fischer_model 3,10 -ksp_guess_fischer_monitor -ksp_guess_fischer_tol 1e-6

156: TEST*/