1 static char help[] = "Test the Fischer-3 initial guess routine.\n\n";
2
3 #include <petscksp.h>
4
5 #define SIZE 3
6
main(int argc,char ** args)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;
15
16 PetscFunctionBeginUser;
17 PetscCall(PetscInitialize(&argc, &args, NULL, help));
18
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));
24
25 PetscCall(VecCreateSeq(PETSC_COMM_SELF, SIZE, &sol));
26 PetscCall(VecDuplicate(sol, &rhs));
27 PetscCall(VecDuplicate(sol, &newrhs));
28 PetscCall(VecDuplicate(sol, &newsol));
29
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));
39
40 /* Test one vector */
41 {
42 KSP ksp;
43 KSPGuess guess;
44
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));
51
52 PetscCall(KSPGuessUpdate(guess, rhs, sol));
53 PetscCall(KSPGuessFormGuess(guess, newrhs, newsol));
54 PetscCall(VecView(newsol, PETSC_VIEWER_STDOUT_SELF));
55
56 PetscCall(KSPDestroy(&ksp));
57 }
58
59 /* Test a singular projection matrix */
60 {
61 KSP ksp;
62 KSPGuess guess;
63
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));
69
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));
73
74 PetscCall(KSPDestroy(&ksp));
75 }
76 PetscCall(VecDestroy(&newsol));
77 PetscCall(VecDestroy(&newrhs));
78 PetscCall(VecDestroy(&rhs));
79 PetscCall(VecDestroy(&sol));
80
81 PetscCall(MatDestroy(&A));
82 }
83
84 /* Test something triangular */
85 {
86 PetscInt triangle_size = 10;
87 Mat A;
88
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));
93
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};
100
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));
106
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));
118
119 PetscCall(KSPGuessUpdate(guess, rhs, sol));
120
121 PetscCall(VecDestroy(&rhs));
122 PetscCall(VecDestroy(&sol));
123 }
124
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));
130
131 PetscCall(KSPGuessFormGuess(guess, rhs, sol));
132 PetscCall(VecView(sol, PETSC_VIEWER_STDOUT_SELF));
133
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 }
143
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 */
145
146 /*TEST
147
148 test:
149 args: -ksp_guess_type fischer -ksp_guess_fischer_model 3,10 -ksp_guess_fischer_monitor -ksp_guess_fischer_tol 1e-6
150
151 TEST*/
152