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