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