1 const char help[] = "Test MATLMVMDIAGBROYDEN by comparing it to MATLMVMSYMBROYDEN for a scalar problem";
2
3 #include <petscksp.h>
4
main(int argc,char ** argv)5 int main(int argc, char **argv)
6 {
7 MPI_Comm comm;
8 Vec x, g, s, y, u, v_diag, v_sym;
9 Mat sym, diag;
10 PetscInt m = 5;
11 PetscReal phi = 0.618;
12 PetscRandom rand;
13 PetscBool is_sym, is_symbad;
14
15 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
16 comm = PETSC_COMM_SELF;
17 PetscCall(VecCreateSeq(comm, 1, &s));
18 PetscCall(VecDuplicate(s, &y));
19 PetscCall(VecDuplicate(s, &x));
20 PetscCall(VecDuplicate(s, &g));
21 PetscCall(VecDuplicate(s, &u));
22 PetscCall(VecDuplicate(s, &v_diag));
23 PetscCall(VecDuplicate(s, &v_sym));
24
25 PetscCall(MatCreateLMVMDiagBroyden(comm, 1, 1, &diag));
26 PetscCall(MatCreateLMVMSymBroyden(comm, 1, 1, &sym));
27
28 PetscCall(MatLMVMSetHistorySize(diag, m));
29 PetscCall(MatLMVMSetHistorySize(sym, m));
30
31 PetscCall(MatSetOptionsPrefix(sym, "sym_"));
32 PetscCall(MatSetFromOptions(sym));
33
34 PetscCall(PetscObjectTypeCompare((PetscObject)sym, MATLMVMSYMBROYDEN, &is_sym));
35 PetscCall(PetscObjectTypeCompare((PetscObject)sym, MATLMVMSYMBADBROYDEN, &is_symbad));
36 if (is_sym) PetscCall(MatLMVMSymBroydenSetPhi(sym, phi));
37 if (is_symbad) PetscCall(MatLMVMSymBadBroydenSetPsi(sym, phi));
38 PetscCall(MatLMVMSymBroydenSetScaleType(sym, MAT_LMVM_SYMBROYDEN_SCALE_NONE));
39
40 PetscCall(MatSetOptionsPrefix(diag, "diag_"));
41 PetscCall(MatSetFromOptions(diag));
42
43 PetscCall(MatSetUp(sym));
44 PetscCall(MatSetUp(diag));
45
46 PetscCall(MatViewFromOptions(sym, NULL, "-view"));
47 PetscCall(MatViewFromOptions(diag, NULL, "-view"));
48
49 PetscCall(PetscRandomCreate(comm, &rand));
50 if (PetscDefined(USE_COMPLEX)) {
51 PetscScalar i = PetscSqrtScalar(-1.0);
52
53 PetscCall(PetscRandomSetInterval(rand, -1.0 - i, 1.0 + i));
54 } else {
55 PetscCall(PetscRandomSetInterval(rand, -1.0, 1.0));
56 }
57
58 PetscCall(VecSetRandom(x, rand));
59 PetscCall(VecSetRandom(g, rand));
60
61 PetscCall(MatLMVMUpdate(sym, x, g));
62 PetscCall(MatLMVMUpdate(diag, x, g));
63
64 for (PetscInt i = 0; i < m; i++) {
65 PetscScalar dot;
66 PetscReal err, scale;
67
68 PetscCall(VecSetRandom(s, rand));
69 PetscCall(VecSetRandom(y, rand));
70
71 PetscCall(VecDot(s, y, &dot));
72 PetscCall(VecScale(s, PetscAbsScalar(dot) / dot));
73
74 PetscCall(VecAXPY(x, 1.0, s));
75 PetscCall(VecAXPY(g, 1.0, y));
76
77 PetscCall(MatLMVMUpdate(sym, x, g));
78 PetscCall(MatLMVMUpdate(diag, x, g));
79
80 PetscCall(VecSet(u, 1.0));
81
82 PetscCall(MatMult(diag, u, v_diag));
83 PetscCall(MatMult(sym, u, v_sym));
84
85 PetscCall(VecAXPY(v_diag, -1.0, v_sym));
86 PetscCall(VecNorm(v_sym, NORM_2, &scale));
87 PetscCall(VecNorm(v_diag, NORM_2, &err));
88 PetscCall(PetscInfo(diag, "Diagonal Broyden error %g, relative error %g\n", (double)err, (double)(err / scale)));
89 PetscCheck(err <= PETSC_SMALL * scale, comm, PETSC_ERR_PLIB, "Diagonal Broyden error %g, relative error %g", (double)err, (double)(err / scale));
90 }
91
92 PetscCall(PetscRandomDestroy(&rand));
93
94 PetscCall(MatDestroy(&sym));
95 PetscCall(MatDestroy(&diag));
96
97 PetscCall(VecDestroy(&v_sym));
98 PetscCall(VecDestroy(&v_diag));
99 PetscCall(VecDestroy(&u));
100 PetscCall(VecDestroy(&g));
101 PetscCall(VecDestroy(&x));
102 PetscCall(VecDestroy(&y));
103 PetscCall(VecDestroy(&s));
104
105 PetscCall(PetscFinalize());
106 return 0;
107 }
108
109 /*TEST
110
111 test:
112 suffix: 0
113 output_file: output/empty.out
114 args: -diag_mat_lmvm_theta 0.618 -diag_mat_lmvm_sigma_hist 0 -diag_mat_lmvm_forward
115
116 test:
117 suffix: 1
118 output_file: output/empty.out
119 args: -diag_mat_lmvm_theta 0.618 -diag_mat_lmvm_sigma_hist 0 -diag_mat_lmvm_forward
120 args: -diag_mat_lmvm_scale_type scalar -diag_mat_lmvm_sigma_hist 1
121
122 test:
123 suffix: 2
124 output_file: output/empty.out
125 args: -sym_mat_type lmvmbfgs -diag_mat_lmvm_theta 0.0 -diag_mat_lmvm_sigma_hist 0
126
127 test:
128 suffix: 3
129 output_file: output/empty.out
130 args: -sym_mat_type lmvmbfgs -diag_mat_lmvm_theta 0.0
131 args: -diag_mat_lmvm_scale_type scalar -diag_mat_lmvm_sigma_hist 1
132
133 test:
134 suffix: 4
135 output_file: output/empty.out
136 args: -sym_mat_type lmvmdfp -diag_mat_lmvm_theta 0.0 -diag_mat_lmvm_sigma_hist 0
137
138 test:
139 suffix: 5
140 output_file: output/empty.out
141 args: -sym_mat_type lmvmdfp -diag_mat_lmvm_theta 0.0
142 args: -diag_mat_lmvm_scale_type scalar -diag_mat_lmvm_sigma_hist 1
143
144 test:
145 suffix: 6
146 output_file: output/empty.out
147 args: -sym_mat_type lmvmsymbadbroyden -diag_mat_lmvm_theta 0.618 -diag_mat_lmvm_sigma_hist 0 -sym_mat_lmvm_scale_type none
148
149 TEST*/
150