1 const char help[] = "Test validity of different approaches to specifying a MatLMVM";
2
3 #include <petscksp.h>
4
5 enum {
6 TEST_INIT_NONE,
7 TEST_INIT_VECS,
8 TEST_INIT_DIAG,
9 TEST_INIT_MAT,
10 TEST_INIT_PC,
11 TEST_INIT_KSP,
12 TEST_INIT_COUNT
13 };
14
15 typedef PetscInt TestInitType;
16
17 enum {
18 TEST_SIZE_NONE,
19 TEST_SIZE_LOCAL,
20 TEST_SIZE_GLOBAL,
21 TEST_SIZE_BOTH,
22 TEST_SIZE_COUNT
23 };
24
25 typedef PetscInt TestSizeType;
26
CreateMatWithTestSizes(MPI_Comm comm,MatType mat_type,PetscInt n,PetscInt N,TestSizeType size_type,PetscBool call_setup,Mat * B)27 static PetscErrorCode CreateMatWithTestSizes(MPI_Comm comm, MatType mat_type, PetscInt n, PetscInt N, TestSizeType size_type, PetscBool call_setup, Mat *B)
28 {
29 PetscFunctionBegin;
30 PetscCall(MatCreate(comm, B));
31 switch (size_type) {
32 case TEST_SIZE_LOCAL:
33 PetscCall(MatSetSizes(*B, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
34 break;
35 case TEST_SIZE_GLOBAL:
36 PetscCall(MatSetSizes(*B, PETSC_DECIDE, PETSC_DECIDE, N, N));
37 break;
38 case TEST_SIZE_BOTH:
39 PetscCall(MatSetSizes(*B, n, n, N, N));
40 break;
41 default:
42 break;
43 }
44 PetscCall(MatSetType(*B, mat_type));
45 if (call_setup) PetscCall(MatSetUp(*B));
46 PetscFunctionReturn(PETSC_SUCCESS);
47 }
48
TestUsability(Mat B,Vec pattern,PetscRandom rand,PetscBool test_pre)49 static PetscErrorCode TestUsability(Mat B, Vec pattern, PetscRandom rand, PetscBool test_pre)
50 {
51 Vec dx, df;
52
53 PetscFunctionBegin;
54 PetscCall(VecDuplicate(pattern, &dx));
55 PetscCall(VecDuplicate(pattern, &df));
56 PetscCall(VecSetRandom(dx, rand));
57
58 if (test_pre) {
59 // Mult and Solve should work before any updates
60 PetscCall(MatMult(B, dx, df));
61 PetscCall(MatSolve(B, dx, df));
62 }
63
64 for (PetscInt i = 0; i < 2; i++) {
65 PetscCall(VecCopy(dx, df));
66 PetscCall(MatLMVMUpdate(B, dx, df));
67 PetscCall(MatMult(B, dx, df));
68 PetscCall(MatSolve(B, dx, df));
69 }
70
71 PetscCall(VecDestroy(&df));
72 PetscCall(VecDestroy(&dx));
73 PetscFunctionReturn(PETSC_SUCCESS);
74 }
75
main(int argc,char ** argv)76 int main(int argc, char **argv)
77 {
78 MPI_Comm comm;
79 PetscRandom rand;
80 PetscInt N = 9;
81 PetscInt n = PETSC_DETERMINE;
82 Vec pattern;
83
84 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
85 PetscCall(KSPInitializePackage());
86 comm = PETSC_COMM_WORLD;
87 PetscCall(PetscRandomCreate(comm, &rand));
88 PetscCall(VecCreate(comm, &pattern));
89 PetscCall(PetscSplitOwnership(comm, &n, &N));
90 PetscCall(VecSetSizes(pattern, n, N));
91 PetscCall(VecSetType(pattern, VECSTANDARD));
92 for (TestInitType init = TEST_INIT_NONE; init < TEST_INIT_COUNT; init++) {
93 for (PetscInt call_setup = 0; call_setup < 2; call_setup++) {
94 for (TestSizeType size = (call_setup ? TEST_SIZE_NONE : TEST_SIZE_LOCAL); size < TEST_SIZE_COUNT; size++) {
95 Mat B;
96
97 // skip invalid case
98 if (init == TEST_INIT_NONE && size == TEST_SIZE_NONE) continue;
99 if (size == TEST_SIZE_NONE && call_setup) continue;
100
101 PetscCall(CreateMatWithTestSizes(comm, MATLMVMBFGS, n, N, size, call_setup ? PETSC_TRUE : PETSC_FALSE, &B));
102
103 switch (init) {
104 case TEST_INIT_NONE:
105 PetscCall(TestUsability(B, pattern, rand, call_setup ? PETSC_TRUE : PETSC_FALSE));
106 break;
107 case TEST_INIT_VECS: {
108 Vec x, f;
109
110 PetscCall(VecDuplicate(pattern, &x));
111 PetscCall(VecDuplicate(pattern, &f));
112 PetscCall(MatLMVMAllocate(B, x, f));
113 PetscCall(VecDestroy(&x));
114 PetscCall(VecDestroy(&f));
115 PetscCall(TestUsability(B, pattern, rand, PETSC_TRUE));
116 break;
117 }
118 case TEST_INIT_DIAG: {
119 Vec d;
120
121 PetscCall(VecDuplicate(pattern, &d));
122 PetscCall(VecSet(d, 1.0));
123 PetscCall(MatLMVMSetJ0Diag(B, d));
124 PetscCall(VecDestroy(&d));
125 PetscCall(TestUsability(B, pattern, rand, PETSC_TRUE));
126 break;
127 }
128 case TEST_INIT_MAT: {
129 for (PetscInt j0_call_setup = 0; j0_call_setup < 2; j0_call_setup++) {
130 for (TestSizeType j0_size = TEST_SIZE_LOCAL; j0_size < TEST_SIZE_COUNT; j0_size++) {
131 Mat J0;
132
133 PetscCall(CreateMatWithTestSizes(comm, MATDENSE, n, N, j0_size, j0_call_setup ? PETSC_TRUE : PETSC_FALSE, &J0));
134 PetscCall(MatLMVMSetJ0(B, J0));
135 PetscCall(MatZeroEntries(J0));
136 PetscCall(MatShift(J0, 1.0));
137 PetscCall(MatDestroy(&J0));
138 PetscCall(TestUsability(B, pattern, rand, PETSC_TRUE));
139 }
140 }
141 break;
142 }
143 case TEST_INIT_PC: {
144 for (PetscInt j0_call_setup = 0; j0_call_setup < 2; j0_call_setup++) {
145 for (TestSizeType j0_size = TEST_SIZE_LOCAL; j0_size < TEST_SIZE_COUNT; j0_size++) {
146 PC J0pc;
147 Mat J0;
148
149 PetscCall(CreateMatWithTestSizes(comm, MATCONSTANTDIAGONAL, n, N, j0_size, j0_call_setup ? PETSC_TRUE : PETSC_FALSE, &J0));
150 PetscCall(PCCreate(comm, &J0pc));
151 PetscCall(PCSetType(J0pc, PCMAT));
152 PetscCall(PCMatSetApplyOperation(J0pc, MATOP_SOLVE));
153 PetscCall(PCSetOperators(J0pc, J0, J0));
154 PetscCall(MatLMVMSetJ0PC(B, J0pc));
155 PetscCall(MatZeroEntries(J0));
156 PetscCall(MatShift(J0, 1.0));
157 PetscCall(MatDestroy(&J0));
158 PetscCall(PCDestroy(&J0pc));
159 PetscCall(TestUsability(B, pattern, rand, PETSC_TRUE));
160 }
161 }
162 break;
163 }
164 case TEST_INIT_KSP: {
165 for (PetscInt j0_call_setup = 0; j0_call_setup < 2; j0_call_setup++) {
166 for (TestSizeType j0_size = TEST_SIZE_LOCAL; j0_size < TEST_SIZE_COUNT; j0_size++) {
167 KSP J0ksp;
168 PC J0pc;
169 Mat J0;
170
171 PetscCall(CreateMatWithTestSizes(comm, MATCONSTANTDIAGONAL, n, N, j0_size, j0_call_setup ? PETSC_TRUE : PETSC_FALSE, &J0));
172 PetscCall(KSPCreate(comm, &J0ksp));
173 PetscCall(KSPSetOperators(J0ksp, J0, J0));
174 PetscCall(KSPGetPC(J0ksp, &J0pc));
175 PetscCall(PCSetType(J0pc, PCNONE));
176 PetscCall(MatLMVMSetJ0KSP(B, J0ksp));
177 PetscCall(MatZeroEntries(J0));
178 PetscCall(MatShift(J0, 1.0));
179 PetscCall(MatDestroy(&J0));
180 PetscCall(KSPDestroy(&J0ksp));
181 PetscCall(TestUsability(B, pattern, rand, PETSC_TRUE));
182 }
183 }
184 break;
185 }
186 default:
187 break;
188 }
189 PetscCall(MatDestroy(&B));
190 }
191 }
192 }
193
194 PetscCall(VecDestroy(&pattern));
195 PetscCall(PetscRandomDestroy(&rand));
196 PetscCall(PetscFinalize());
197 return 0;
198 }
199
200 /*TEST
201
202 test:
203 suffix: 0
204 nsize: 2
205 output_file: output/empty.out
206
207 TEST*/
208