1 const char help[] = "Test that MatCopy() does not affect the copied LMVM matrix";
2
3 #include <petsc.h>
4
positiveVectorUpdate(PetscRandom rand,Vec x,Vec f)5 static PetscErrorCode positiveVectorUpdate(PetscRandom rand, Vec x, Vec f)
6 {
7 Vec x_, f_;
8 PetscScalar dot;
9
10 PetscFunctionBegin;
11 PetscCall(VecDuplicate(x, &x_));
12 PetscCall(VecDuplicate(f, &f_));
13 PetscCall(VecSetRandom(x_, rand));
14 PetscCall(VecSetRandom(f_, rand));
15 PetscCall(VecDot(x_, f_, &dot));
16 PetscCall(VecAXPY(x, PetscAbsScalar(dot) / dot, x_));
17 PetscCall(VecAXPY(f, 1.0, f_));
18 PetscCall(VecDestroy(&f_));
19 PetscCall(VecDestroy(&x_));
20 PetscFunctionReturn(PETSC_SUCCESS);
21 }
22
VecEqualToTolerance(Vec a,Vec b,NormType norm_type,PetscReal tol,PetscBool * flg)23 static PetscErrorCode VecEqualToTolerance(Vec a, Vec b, NormType norm_type, PetscReal tol, PetscBool *flg)
24 {
25 Vec diff;
26 PetscReal diff_norm;
27
28 PetscFunctionBegin;
29 PetscCall(VecDuplicate(a, &diff));
30 PetscCall(VecCopy(a, diff));
31 PetscCall(VecAXPY(diff, -1.0, b));
32 PetscCall(VecNorm(diff, norm_type, &diff_norm));
33 PetscCall(VecDestroy(&diff));
34 *flg = (diff_norm <= tol) ? PETSC_TRUE : PETSC_FALSE;
35 PetscFunctionReturn(PETSC_SUCCESS);
36 }
37
38 // unlike MatTestEqual(), this test tests MatMult() and MatSolve()
testMatEqual(PetscRandom rand,Mat A,Mat B,PetscBool * flg)39 static PetscErrorCode testMatEqual(PetscRandom rand, Mat A, Mat B, PetscBool *flg)
40 {
41 Vec x, y_A, y_B;
42
43 PetscFunctionBegin;
44 *flg = PETSC_TRUE;
45 PetscCall(MatCreateVecs(A, &x, &y_A));
46 PetscCall(MatCreateVecs(B, NULL, &y_B));
47 PetscCall(VecSetRandom(x, rand));
48 PetscCall(MatMult(A, x, y_A));
49 PetscCall(MatMult(B, x, y_B));
50 PetscCall(VecEqualToTolerance(y_A, y_B, NORM_2, PETSC_SMALL, flg));
51 if (*flg == PETSC_TRUE) {
52 PetscCall(MatSolve(A, x, y_A));
53 PetscCall(MatSolve(B, x, y_B));
54 PetscCall(VecEqualToTolerance(y_A, y_B, NORM_2, PETSC_SMALL, flg));
55 if (*flg == PETSC_FALSE) {
56 PetscReal norm;
57
58 PetscCall(VecAXPY(y_A, -1.0, y_B));
59 PetscCall(VecNorm(y_A, NORM_INFINITY, &norm));
60 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)A), "MatSolve() norm error %g\n", (double)norm));
61 }
62 } else {
63 PetscReal norm;
64
65 PetscCall(VecAXPY(y_A, -1.0, y_B));
66 PetscCall(VecNorm(y_A, NORM_INFINITY, &norm));
67 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)A), "MatMult() norm error %g\n", (double)norm));
68 }
69 PetscCall(VecDestroy(&y_B));
70 PetscCall(VecDestroy(&y_A));
71 PetscCall(VecDestroy(&x));
72 PetscFunctionReturn(PETSC_SUCCESS);
73 }
74
testUnchangedBegin(PetscRandom rand,Mat A,Vec * x,Vec * y,Vec * z)75 static PetscErrorCode testUnchangedBegin(PetscRandom rand, Mat A, Vec *x, Vec *y, Vec *z)
76 {
77 Vec x_, y_, z_;
78
79 PetscFunctionBegin;
80 PetscCall(MatCreateVecs(A, &x_, &y_));
81 PetscCall(MatCreateVecs(A, NULL, &z_));
82 PetscCall(VecSetRandom(x_, rand));
83 PetscCall(MatMult(A, x_, y_));
84 PetscCall(MatSolve(A, x_, z_));
85 *x = x_;
86 *y = y_;
87 *z = z_;
88 PetscFunctionReturn(PETSC_SUCCESS);
89 }
90
testUnchangedEnd(PetscRandom rand,Mat A,Vec * x,Vec * y,Vec * z,PetscBool * unchanged)91 static PetscErrorCode testUnchangedEnd(PetscRandom rand, Mat A, Vec *x, Vec *y, Vec *z, PetscBool *unchanged)
92 {
93 Vec x_, y_, z_, y2_, z2_;
94
95 PetscFunctionBegin;
96 *unchanged = PETSC_TRUE;
97 x_ = *x;
98 y_ = *y;
99 z_ = *z;
100 *x = NULL;
101 *y = NULL;
102 *z = NULL;
103 PetscCall(MatCreateVecs(A, NULL, &y2_));
104 PetscCall(MatCreateVecs(A, NULL, &z2_));
105 PetscCall(MatMult(A, x_, y2_));
106 PetscCall(MatSolve(A, x_, z2_));
107 PetscCall(VecEqual(y_, y2_, unchanged));
108 if (*unchanged == PETSC_TRUE) PetscCall(VecEqual(z_, z2_, unchanged));
109 PetscCall(VecDestroy(&z2_));
110 PetscCall(VecDestroy(&y2_));
111 PetscCall(VecDestroy(&z_));
112 PetscCall(VecDestroy(&y_));
113 PetscCall(VecDestroy(&x_));
114 PetscFunctionReturn(PETSC_SUCCESS);
115 }
116
testMatLMVMCopy(PetscRandom rand)117 static PetscErrorCode testMatLMVMCopy(PetscRandom rand)
118 {
119 PetscInt N = 10;
120 MPI_Comm comm = PetscObjectComm((PetscObject)rand);
121 PetscInt k_pre = 2; // number of updates before copy
122 Mat A, A_copy;
123 Vec u, x, f, x_copy, f_copy, v1, v2, v3;
124 PetscBool equal;
125 PetscLayout layout;
126
127 PetscFunctionBegin;
128 PetscCall(VecCreateMPI(comm, PETSC_DECIDE, N, &u));
129 PetscCall(VecSetFromOptions(u));
130 PetscCall(VecSetUp(u));
131 PetscCall(VecDuplicate(u, &x));
132 PetscCall(VecDuplicate(u, &f));
133 PetscCall(VecGetLayout(u, &layout));
134 PetscCall(MatCreate(comm, &A));
135 PetscCall(MatSetLayouts(A, layout, layout));
136 PetscCall(MatSetType(A, MATLMVMBFGS));
137 PetscCall(MatSetFromOptions(A));
138 PetscCall(positiveVectorUpdate(rand, x, f));
139 PetscCall(MatLMVMAllocate(A, x, f));
140 PetscCall(MatSetUp(A));
141 for (PetscInt k = 0; k <= k_pre; k++) {
142 PetscCall(positiveVectorUpdate(rand, x, f));
143 PetscCall(MatLMVMUpdate(A, x, f));
144 }
145 PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A_copy));
146 PetscCall(testMatEqual(rand, A, A_copy, &equal));
147 PetscCheck(equal, comm, PETSC_ERR_PLIB, "MatCopy() not the same after initial copy");
148
149 PetscCall(VecDuplicate(x, &x_copy));
150 PetscCall(VecCopy(x, x_copy));
151 PetscCall(VecDuplicate(f, &f_copy));
152 PetscCall(VecCopy(f, f_copy));
153
154 PetscCall(testUnchangedBegin(rand, A_copy, &v1, &v2, &v3));
155 PetscCall(positiveVectorUpdate(rand, x, f));
156 PetscCall(MatLMVMUpdate(A, x, f));
157 PetscCall(testUnchangedEnd(rand, A_copy, &v1, &v2, &v3, &equal));
158 PetscCheck(equal, comm, PETSC_ERR_PLIB, "MatLMVMUpdate() to original matrix affects copy");
159
160 PetscCall(testUnchangedBegin(rand, A, &v1, &v2, &v3));
161 PetscCall(positiveVectorUpdate(rand, x_copy, f_copy));
162 PetscCall(MatLMVMUpdate(A_copy, x_copy, f_copy));
163 PetscCall(testUnchangedEnd(rand, A, &v1, &v2, &v3, &equal));
164 PetscCheck(equal, comm, PETSC_ERR_PLIB, "MatLMVMUpdate() to copy matrix affects original");
165
166 PetscCall(VecDestroy(&f_copy));
167 PetscCall(VecDestroy(&x_copy));
168 PetscCall(MatDestroy(&A_copy));
169 PetscCall(MatDestroy(&A));
170 PetscCall(VecDestroy(&f));
171 PetscCall(VecDestroy(&x));
172 PetscCall(VecDestroy(&u));
173 PetscFunctionReturn(PETSC_SUCCESS);
174 }
175
main(int argc,char ** argv)176 int main(int argc, char **argv)
177 {
178 MPI_Comm comm;
179 PetscRandom rand;
180
181 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
182 comm = PETSC_COMM_WORLD;
183 PetscCall(PetscRandomCreate(comm, &rand));
184 PetscCall(PetscRandomSetFromOptions(rand));
185 PetscCall(KSPInitializePackage());
186 PetscCall(testMatLMVMCopy(rand));
187 PetscCall(PetscRandomDestroy(&rand));
188 PetscCall(PetscFinalize());
189 return 0;
190 }
191
192 /*TEST
193
194 # dense != compact_dense
195 test:
196 suffix: 0
197 output_file: output/empty.out
198 args: -mat_lmvm_mult_algorithm {{recursive dense compact_dense}} -mat_type {{lmvmbfgs lmvmdfp lmvmbroyden lmvmbadbroyden}}
199
200 # dense == compact_dense
201 test:
202 suffix: 1
203 output_file: output/empty.out
204 args: -mat_lmvm_mult_algorithm {{recursive dense}} -mat_type {{lmvmsr1 lmvmsymbroyden lmvmsymbadbroyden}}
205
206 test:
207 suffix: 2
208 output_file: output/empty.out
209 args: -mat_type {{lmvmdiagbroyden lmvmdqn}}
210
211 test:
212 suffix: 3
213 output_file: output/empty.out
214 args: -mat_type lmvmdbfgs -mat_lbfgs_recursive {{0 1}}
215
216 test:
217 suffix: 4
218 output_file: output/empty.out
219 args: -mat_type lmvmddfp -mat_ldfp_recursive {{0 1}}
220
221 TEST*/
222