Actual source code: lmvm_copy_test.c
1: const char help[] = "Test that MatCopy() does not affect the copied LMVM matrix";
3: #include <petsc.h>
5: static PetscErrorCode positiveVectorUpdate(PetscRandom rand, Vec x, Vec f)
6: {
7: Vec _x, _f;
8: PetscScalar dot;
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: }
23: static PetscErrorCode VecEqualToTolerance(Vec a, Vec b, NormType norm_type, PetscReal tol, PetscBool *flg)
24: {
25: Vec diff;
26: PetscReal diff_norm;
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: }
38: // unlike MatTestEqual(), this test tests MatMult() and MatSolve()
39: static PetscErrorCode testMatEqual(PetscRandom rand, Mat A, Mat B, PetscBool *flg)
40: {
41: Vec x, y_A, y_B;
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;
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;
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: }
75: static PetscErrorCode testUnchangedBegin(PetscRandom rand, Mat A, Vec *x, Vec *y, Vec *z)
76: {
77: Vec _x, _y, _z;
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: }
91: static PetscErrorCode testUnchangedEnd(PetscRandom rand, Mat A, Vec *x, Vec *y, Vec *z, PetscBool *unchanged)
92: {
93: Vec _x, _y, _z, _y2, _z2;
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: }
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;
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");
149: PetscCall(VecDuplicate(x, &x_copy));
150: PetscCall(VecCopy(x, x_copy));
151: PetscCall(VecDuplicate(f, &f_copy));
152: PetscCall(VecCopy(f, f_copy));
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");
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");
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: }
176: int main(int argc, char **argv)
177: {
178: MPI_Comm comm;
179: PetscRandom rand;
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: }
192: /*TEST
194: test:
195: suffix: 0
196: args: -mat_type {{lmvmbfgs lmvmdfp lmvmsr1 lmvmbroyden lmvmbadbroyden lmvmsymbroyden lmvmsymbadbroyden lmvmdiagbroyden lmvmdbfgs lmvmddfp lmvmdqn}}
198: TEST*/