1 #include <../src/ksp/ksp/utils/lmvm/lmvm.h> /*I "petscksp.h" I*/
2 
3 /*------------------------------------------------------------*/
4 
MatReset_LMVM(Mat B,PetscBool destructive)5 PetscErrorCode MatReset_LMVM(Mat B, PetscBool destructive)
6 {
7   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
8   PetscErrorCode    ierr;
9 
10   PetscFunctionBegin;
11   lmvm->k = -1;
12   lmvm->prev_set = PETSC_FALSE;
13   lmvm->shift = 0.0;
14   if (destructive && lmvm->allocated) {
15     ierr = MatLMVMClearJ0(B);CHKERRQ(ierr);
16     B->rmap->n = B->rmap->N = B->cmap->n = B->cmap->N = 0;
17     ierr = VecDestroyVecs(lmvm->m, &lmvm->S);CHKERRQ(ierr);
18     ierr = VecDestroyVecs(lmvm->m, &lmvm->Y);CHKERRQ(ierr);
19     ierr = VecDestroy(&lmvm->Xprev);CHKERRQ(ierr);
20     ierr = VecDestroy(&lmvm->Fprev);CHKERRQ(ierr);
21     lmvm->nupdates = 0;
22     lmvm->nrejects = 0;
23     lmvm->m_old = 0;
24     lmvm->allocated = PETSC_FALSE;
25     B->preallocated = PETSC_FALSE;
26     B->assembled = PETSC_FALSE;
27   }
28   ++lmvm->nresets;
29   PetscFunctionReturn(0);
30 }
31 
32 /*------------------------------------------------------------*/
33 
MatAllocate_LMVM(Mat B,Vec X,Vec F)34 PetscErrorCode MatAllocate_LMVM(Mat B, Vec X, Vec F)
35 {
36   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
37   PetscErrorCode    ierr;
38   PetscBool         same, allocate = PETSC_FALSE;
39   PetscInt          m, n, M, N;
40   VecType           type;
41 
42   PetscFunctionBegin;
43   if (lmvm->allocated) {
44     VecCheckMatCompatible(B, X, 2, F, 3);
45     ierr = VecGetType(X, &type);CHKERRQ(ierr);
46     ierr = PetscObjectTypeCompare((PetscObject)lmvm->Xprev, type, &same);CHKERRQ(ierr);
47     if (!same) {
48       /* Given X vector has a different type than allocated X-type data structures.
49          We need to destroy all of this and duplicate again out of the given vector. */
50       allocate = PETSC_TRUE;
51       ierr = MatLMVMReset(B, PETSC_TRUE);CHKERRQ(ierr);
52     }
53   } else {
54     allocate = PETSC_TRUE;
55   }
56   if (allocate) {
57     ierr = VecGetLocalSize(X, &n);CHKERRQ(ierr);
58     ierr = VecGetSize(X, &N);CHKERRQ(ierr);
59     ierr = VecGetLocalSize(F, &m);CHKERRQ(ierr);
60     ierr = VecGetSize(F, &M);CHKERRQ(ierr);
61     B->rmap->n = m;
62     B->cmap->n = n;
63     B->rmap->N = M > -1 ? M : B->rmap->N;
64     B->cmap->N = N > -1 ? N : B->cmap->N;
65     ierr = VecDuplicate(X, &lmvm->Xprev);CHKERRQ(ierr);
66     ierr = VecDuplicate(F, &lmvm->Fprev);CHKERRQ(ierr);
67     if (lmvm->m > 0) {
68       ierr = VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lmvm->S);CHKERRQ(ierr);
69       ierr = VecDuplicateVecs(lmvm->Fprev, lmvm->m, &lmvm->Y);CHKERRQ(ierr);
70     }
71     lmvm->m_old = lmvm->m;
72     lmvm->allocated = PETSC_TRUE;
73     B->preallocated = PETSC_TRUE;
74     B->assembled = PETSC_TRUE;
75   }
76   PetscFunctionReturn(0);
77 }
78 
79 /*------------------------------------------------------------*/
80 
MatUpdateKernel_LMVM(Mat B,Vec S,Vec Y)81 PetscErrorCode MatUpdateKernel_LMVM(Mat B, Vec S, Vec Y)
82 {
83   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
84   PetscErrorCode    ierr;
85   PetscInt          i;
86   Vec               Stmp, Ytmp;
87 
88   PetscFunctionBegin;
89   if (lmvm->k == lmvm->m-1) {
90     /* We hit the memory limit, so shift all the vectors back one spot
91        and shift the oldest to the front to receive the latest update. */
92     Stmp = lmvm->S[0];
93     Ytmp = lmvm->Y[0];
94     for (i = 0; i < lmvm->k; ++i) {
95       lmvm->S[i] = lmvm->S[i+1];
96       lmvm->Y[i] = lmvm->Y[i+1];
97     }
98     lmvm->S[lmvm->k] = Stmp;
99     lmvm->Y[lmvm->k] = Ytmp;
100   } else {
101     ++lmvm->k;
102   }
103   /* Put the precomputed update into the last vector */
104   ierr = VecCopy(S, lmvm->S[lmvm->k]);CHKERRQ(ierr);
105   ierr = VecCopy(Y, lmvm->Y[lmvm->k]);CHKERRQ(ierr);
106   ++lmvm->nupdates;
107   PetscFunctionReturn(0);
108 }
109 
110 /*------------------------------------------------------------*/
111 
MatUpdate_LMVM(Mat B,Vec X,Vec F)112 PetscErrorCode MatUpdate_LMVM(Mat B, Vec X, Vec F)
113 {
114   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
115   PetscErrorCode    ierr;
116 
117   PetscFunctionBegin;
118   if (!lmvm->m) PetscFunctionReturn(0);
119   if (lmvm->prev_set) {
120     /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
121     ierr = VecAXPBY(lmvm->Xprev, 1.0, -1.0, X);CHKERRQ(ierr);
122     ierr = VecAXPBY(lmvm->Fprev, 1.0, -1.0, F);CHKERRQ(ierr);
123     /* Update S and Y */
124     ierr = MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev);CHKERRQ(ierr);
125   }
126 
127   /* Save the solution and function to be used in the next update */
128   ierr = VecCopy(X, lmvm->Xprev);CHKERRQ(ierr);
129   ierr = VecCopy(F, lmvm->Fprev);CHKERRQ(ierr);
130   lmvm->prev_set = PETSC_TRUE;
131   PetscFunctionReturn(0);
132 }
133 
134 /*------------------------------------------------------------*/
135 
MatMultAdd_LMVM(Mat B,Vec X,Vec Y,Vec Z)136 static PetscErrorCode MatMultAdd_LMVM(Mat B, Vec X, Vec Y, Vec Z)
137 {
138   PetscErrorCode    ierr;
139 
140   PetscFunctionBegin;
141   ierr = MatMult(B, X, Z);CHKERRQ(ierr);
142   ierr = VecAXPY(Z, 1.0, Y);CHKERRQ(ierr);
143   PetscFunctionReturn(0);
144 }
145 
146 /*------------------------------------------------------------*/
147 
MatMult_LMVM(Mat B,Vec X,Vec Y)148 static PetscErrorCode MatMult_LMVM(Mat B, Vec X, Vec Y)
149 {
150   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
151   PetscErrorCode    ierr;
152 
153   PetscFunctionBegin;
154   VecCheckSameSize(X, 2, Y, 3);
155   VecCheckMatCompatible(B, X, 2, Y, 3);
156   if (!lmvm->allocated) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ORDER, "LMVM matrix must be allocated first");
157   ierr = (*lmvm->ops->mult)(B, X, Y);CHKERRQ(ierr);
158   if (lmvm->shift != 0.0) {
159     ierr = VecAXPY(Y, lmvm->shift, X);CHKERRQ(ierr);
160   }
161   PetscFunctionReturn(0);
162 }
163 
164 /*------------------------------------------------------------*/
165 
MatCopy_LMVM(Mat B,Mat M,MatStructure str)166 static PetscErrorCode MatCopy_LMVM(Mat B, Mat M, MatStructure str)
167 {
168   Mat_LMVM          *bctx = (Mat_LMVM*)B->data;
169   Mat_LMVM          *mctx;
170   PetscErrorCode    ierr;
171   PetscInt          i;
172   PetscBool         allocatedM;
173 
174   PetscFunctionBegin;
175   if (str == DIFFERENT_NONZERO_PATTERN) {
176     ierr = MatLMVMReset(M, PETSC_TRUE);CHKERRQ(ierr);
177     ierr = MatLMVMAllocate(M, bctx->Xprev, bctx->Fprev);CHKERRQ(ierr);
178   } else {
179     ierr = MatLMVMIsAllocated(M, &allocatedM);CHKERRQ(ierr);
180     if (!allocatedM) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONGSTATE, "Target matrix must be allocated first");
181     MatCheckSameSize(B, 1, M, 2);
182   }
183 
184   mctx = (Mat_LMVM*)M->data;
185   if (bctx->user_pc) {
186     ierr = MatLMVMSetJ0PC(M, bctx->J0pc);CHKERRQ(ierr);
187   } else if (bctx->user_ksp) {
188     ierr = MatLMVMSetJ0KSP(M, bctx->J0ksp);CHKERRQ(ierr);
189   } else if (bctx->J0) {
190     ierr = MatLMVMSetJ0(M, bctx->J0);CHKERRQ(ierr);
191   } else if (bctx->user_scale) {
192     if (bctx->J0diag) {
193       ierr = MatLMVMSetJ0Diag(M, bctx->J0diag);CHKERRQ(ierr);
194     } else {
195       ierr = MatLMVMSetJ0Scale(M, bctx->J0scalar);CHKERRQ(ierr);
196     }
197   }
198   mctx->nupdates = bctx->nupdates;
199   mctx->nrejects = bctx->nrejects;
200   mctx->k = bctx->k;
201   for (i=0; i<=bctx->k; ++i) {
202     ierr = VecCopy(bctx->S[i], mctx->S[i]);CHKERRQ(ierr);
203     ierr = VecCopy(bctx->Y[i], mctx->Y[i]);CHKERRQ(ierr);
204     ierr = VecCopy(bctx->Xprev, mctx->Xprev);CHKERRQ(ierr);
205     ierr = VecCopy(bctx->Fprev, mctx->Fprev);CHKERRQ(ierr);
206   }
207   if (bctx->ops->copy) {
208     ierr = (*bctx->ops->copy)(B, M, str);CHKERRQ(ierr);
209   }
210   PetscFunctionReturn(0);
211 }
212 
213 /*------------------------------------------------------------*/
214 
MatDuplicate_LMVM(Mat B,MatDuplicateOption op,Mat * mat)215 static PetscErrorCode MatDuplicate_LMVM(Mat B, MatDuplicateOption op, Mat *mat)
216 {
217   Mat_LMVM          *bctx = (Mat_LMVM*)B->data;
218   Mat_LMVM          *mctx;
219   PetscErrorCode    ierr;
220   MatType           lmvmType;
221   Mat               A;
222 
223   PetscFunctionBegin;
224   ierr = MatGetType(B, &lmvmType);CHKERRQ(ierr);
225   ierr = MatCreate(PetscObjectComm((PetscObject)B), mat);CHKERRQ(ierr);
226   ierr = MatSetType(*mat, lmvmType);CHKERRQ(ierr);
227 
228   A = *mat;
229   mctx = (Mat_LMVM*)A->data;
230   mctx->m = bctx->m;
231   mctx->ksp_max_it = bctx->ksp_max_it;
232   mctx->ksp_rtol = bctx->ksp_rtol;
233   mctx->ksp_atol = bctx->ksp_atol;
234   mctx->shift = bctx->shift;
235   ierr = KSPSetTolerances(mctx->J0ksp, mctx->ksp_rtol, mctx->ksp_atol, PETSC_DEFAULT, mctx->ksp_max_it);CHKERRQ(ierr);
236 
237   ierr = MatLMVMAllocate(*mat, bctx->Xprev, bctx->Fprev);CHKERRQ(ierr);
238   if (op == MAT_COPY_VALUES) {
239     ierr = MatCopy(B, *mat, SAME_NONZERO_PATTERN);CHKERRQ(ierr);
240   }
241   PetscFunctionReturn(0);
242 }
243 
244 /*------------------------------------------------------------*/
245 
MatShift_LMVM(Mat B,PetscScalar a)246 static PetscErrorCode MatShift_LMVM(Mat B, PetscScalar a)
247 {
248   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
249 
250   PetscFunctionBegin;
251   if (!lmvm->allocated) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ORDER, "LMVM matrix must be allocated first");
252   lmvm->shift += PetscRealPart(a);
253   PetscFunctionReturn(0);
254 }
255 
256 /*------------------------------------------------------------*/
257 
MatGetVecs_LMVM(Mat B,Vec * L,Vec * R)258 static PetscErrorCode MatGetVecs_LMVM(Mat B, Vec *L, Vec *R)
259 {
260   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
261   PetscErrorCode    ierr;
262 
263   PetscFunctionBegin;
264   if (!lmvm->allocated) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ORDER, "LMVM matrix must be allocated first");
265   ierr = VecDuplicate(lmvm->Xprev, L);CHKERRQ(ierr);
266   ierr = VecDuplicate(lmvm->Fprev, R);CHKERRQ(ierr);
267   PetscFunctionReturn(0);
268 }
269 
270 /*------------------------------------------------------------*/
271 
MatView_LMVM(Mat B,PetscViewer pv)272 PetscErrorCode MatView_LMVM(Mat B, PetscViewer pv)
273 {
274   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
275   PetscErrorCode    ierr;
276   PetscBool         isascii;
277   MatType           type;
278 
279   PetscFunctionBegin;
280   ierr = PetscObjectTypeCompare((PetscObject)pv,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr);
281   if (isascii) {
282     ierr = MatGetType(B, &type);CHKERRQ(ierr);
283     ierr = PetscViewerASCIIPrintf(pv,"Max. storage: %D\n",lmvm->m);CHKERRQ(ierr);
284     ierr = PetscViewerASCIIPrintf(pv,"Used storage: %D\n",lmvm->k+1);CHKERRQ(ierr);
285     ierr = PetscViewerASCIIPrintf(pv,"Number of updates: %D\n",lmvm->nupdates);CHKERRQ(ierr);
286     ierr = PetscViewerASCIIPrintf(pv,"Number of rejects: %D\n",lmvm->nrejects);CHKERRQ(ierr);
287     ierr = PetscViewerASCIIPrintf(pv,"Number of resets: %D\n",lmvm->nresets);CHKERRQ(ierr);
288     if (lmvm->J0) {
289       ierr = PetscViewerASCIIPrintf(pv,"J0 Matrix:\n");CHKERRQ(ierr);
290       ierr = PetscViewerPushFormat(pv, PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr);
291       ierr = MatView(lmvm->J0, pv);CHKERRQ(ierr);
292       ierr = PetscViewerPopFormat(pv);CHKERRQ(ierr);
293     }
294   }
295   PetscFunctionReturn(0);
296 }
297 
298 /*------------------------------------------------------------*/
299 
MatSetFromOptions_LMVM(PetscOptionItems * PetscOptionsObject,Mat B)300 PetscErrorCode MatSetFromOptions_LMVM(PetscOptionItems *PetscOptionsObject, Mat B)
301 {
302   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
303   PetscErrorCode    ierr;
304 
305   PetscFunctionBegin;
306   ierr = PetscOptionsHead(PetscOptionsObject,"Limited-memory Variable Metric matrix for approximating Jacobians");CHKERRQ(ierr);
307   ierr = PetscOptionsInt("-mat_lmvm_hist_size","number of past updates kept in memory for the approximation","",lmvm->m,&lmvm->m,NULL);CHKERRQ(ierr);
308   ierr = PetscOptionsInt("-mat_lmvm_ksp_its","(developer) fixed number of KSP iterations to take when inverting J0","",lmvm->ksp_max_it,&lmvm->ksp_max_it,NULL);CHKERRQ(ierr);
309   ierr = PetscOptionsReal("-mat_lmvm_eps","(developer) machine zero definition","",lmvm->eps,&lmvm->eps,NULL);CHKERRQ(ierr);
310   ierr = PetscOptionsTail();CHKERRQ(ierr);
311   ierr = KSPSetFromOptions(lmvm->J0ksp);CHKERRQ(ierr);
312   PetscFunctionReturn(0);
313 }
314 
315 /*------------------------------------------------------------*/
316 
MatSetUp_LMVM(Mat B)317 PetscErrorCode MatSetUp_LMVM(Mat B)
318 {
319   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
320   PetscErrorCode    ierr;
321   PetscInt          m, n, M, N;
322   PetscMPIInt       size;
323   MPI_Comm          comm = PetscObjectComm((PetscObject)B);
324 
325   PetscFunctionBegin;
326   ierr = MatGetSize(B, &M, &N);CHKERRQ(ierr);
327   if (M == 0 && N == 0) SETERRQ(comm, PETSC_ERR_ORDER, "MatSetSizes() must be called before MatSetUp()");
328   if (!lmvm->allocated) {
329     ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr);
330     if (size == 1) {
331       ierr = VecCreateSeq(comm, N, &lmvm->Xprev);CHKERRQ(ierr);
332       ierr = VecCreateSeq(comm, M, &lmvm->Fprev);CHKERRQ(ierr);
333     } else {
334       ierr = MatGetLocalSize(B, &m, &n);CHKERRQ(ierr);
335       ierr = VecCreateMPI(comm, n, N, &lmvm->Xprev);CHKERRQ(ierr);
336       ierr = VecCreateMPI(comm, m, M, &lmvm->Fprev);CHKERRQ(ierr);
337     }
338     if (lmvm->m > 0) {
339       ierr = VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lmvm->S);CHKERRQ(ierr);
340       ierr = VecDuplicateVecs(lmvm->Fprev, lmvm->m, &lmvm->Y);CHKERRQ(ierr);
341     }
342     lmvm->m_old = lmvm->m;
343     lmvm->allocated = PETSC_TRUE;
344     B->preallocated = PETSC_TRUE;
345     B->assembled = PETSC_TRUE;
346   }
347   PetscFunctionReturn(0);
348 }
349 
350 /*------------------------------------------------------------*/
351 
MatDestroy_LMVM(Mat B)352 PetscErrorCode MatDestroy_LMVM(Mat B)
353 {
354   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
355   PetscErrorCode    ierr;
356 
357   PetscFunctionBegin;
358   if (lmvm->allocated) {
359     ierr = VecDestroyVecs(lmvm->m, &lmvm->S);CHKERRQ(ierr);
360     ierr = VecDestroyVecs(lmvm->m, &lmvm->Y);CHKERRQ(ierr);
361     ierr = VecDestroy(&lmvm->Xprev);CHKERRQ(ierr);
362     ierr = VecDestroy(&lmvm->Fprev);CHKERRQ(ierr);
363   }
364   ierr = KSPDestroy(&lmvm->J0ksp);CHKERRQ(ierr);
365   ierr = MatLMVMClearJ0(B);CHKERRQ(ierr);
366   ierr = PetscFree(B->data);CHKERRQ(ierr);
367   PetscFunctionReturn(0);
368 }
369 
370 /*------------------------------------------------------------*/
371 
MatCreate_LMVM(Mat B)372 PetscErrorCode MatCreate_LMVM(Mat B)
373 {
374   Mat_LMVM          *lmvm;
375   PetscErrorCode    ierr;
376 
377   PetscFunctionBegin;
378   ierr = PetscNewLog(B, &lmvm);CHKERRQ(ierr);
379   B->data = (void*)lmvm;
380 
381   lmvm->m_old = 0;
382   lmvm->m = 5;
383   lmvm->k = -1;
384   lmvm->nupdates = 0;
385   lmvm->nrejects = 0;
386   lmvm->nresets = 0;
387 
388   lmvm->ksp_max_it = 20;
389   lmvm->ksp_rtol = 0.0;
390   lmvm->ksp_atol = 0.0;
391 
392   lmvm->shift = 0.0;
393 
394   lmvm->eps = PetscPowReal(PETSC_MACHINE_EPSILON, 2.0/3.0);
395   lmvm->allocated = PETSC_FALSE;
396   lmvm->prev_set = PETSC_FALSE;
397   lmvm->user_scale = PETSC_FALSE;
398   lmvm->user_pc = PETSC_FALSE;
399   lmvm->user_ksp = PETSC_FALSE;
400   lmvm->square = PETSC_FALSE;
401 
402   B->ops->destroy = MatDestroy_LMVM;
403   B->ops->setfromoptions = MatSetFromOptions_LMVM;
404   B->ops->view = MatView_LMVM;
405   B->ops->setup = MatSetUp_LMVM;
406   B->ops->getvecs = MatGetVecs_LMVM;
407   B->ops->shift = MatShift_LMVM;
408   B->ops->duplicate = MatDuplicate_LMVM;
409   B->ops->mult = MatMult_LMVM;
410   B->ops->multadd = MatMultAdd_LMVM;
411   B->ops->copy = MatCopy_LMVM;
412 
413   lmvm->ops->update = MatUpdate_LMVM;
414   lmvm->ops->allocate = MatAllocate_LMVM;
415   lmvm->ops->reset = MatReset_LMVM;
416 
417   ierr = KSPCreate(PetscObjectComm((PetscObject)B), &lmvm->J0ksp);CHKERRQ(ierr);
418   ierr = PetscObjectIncrementTabLevel((PetscObject)lmvm->J0ksp, (PetscObject)B, 1);CHKERRQ(ierr);
419   ierr = KSPSetOptionsPrefix(lmvm->J0ksp, "mat_lmvm_");CHKERRQ(ierr);
420   ierr = KSPSetType(lmvm->J0ksp, KSPGMRES);CHKERRQ(ierr);
421   ierr = KSPSetTolerances(lmvm->J0ksp, lmvm->ksp_rtol, lmvm->ksp_atol, PETSC_DEFAULT, lmvm->ksp_max_it);CHKERRQ(ierr);
422   PetscFunctionReturn(0);
423 }
424