Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
6 : This file is part of SLEPc.
7 : SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9 : */
10 : /*
11 : ST interface routines, callable by users
12 : */
13 :
14 : #include <slepc/private/stimpl.h> /*I "slepcst.h" I*/
15 :
16 144069 : PetscErrorCode STApply_Generic(ST st,Vec x,Vec y)
17 : {
18 144069 : PetscFunctionBegin;
19 144069 : if (st->M && st->P) {
20 57499 : PetscCall(MatMult(st->M,x,st->work[0]));
21 57499 : PetscCall(STMatSolve(st,st->work[0],y));
22 86570 : } else if (st->M) PetscCall(MatMult(st->M,x,y));
23 15704 : else PetscCall(STMatSolve(st,x,y));
24 144069 : PetscFunctionReturn(PETSC_SUCCESS);
25 : }
26 :
27 : /*@
28 : STApply - Applies the spectral transformation operator to a vector, for
29 : instance (A - sB)^-1 B in the case of the shift-and-invert transformation
30 : and generalized eigenproblem.
31 :
32 : Collective
33 :
34 : Input Parameters:
35 : + st - the spectral transformation context
36 : - x - input vector
37 :
38 : Output Parameter:
39 : . y - output vector
40 :
41 : Level: developer
42 :
43 : .seealso: STApplyTranspose(), STApplyHermitianTranspose()
44 : @*/
45 84725 : PetscErrorCode STApply(ST st,Vec x,Vec y)
46 : {
47 84725 : Mat Op;
48 :
49 84725 : PetscFunctionBegin;
50 84725 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
51 84725 : PetscValidHeaderSpecific(x,VEC_CLASSID,2);
52 84725 : PetscValidHeaderSpecific(y,VEC_CLASSID,3);
53 84725 : PetscValidType(st,1);
54 84725 : STCheckMatrices(st,1);
55 84725 : PetscCheck(x!=y,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
56 84725 : PetscCall(VecSetErrorIfLocked(y,3));
57 84725 : PetscCall(STGetOperator_Private(st,&Op));
58 84725 : PetscCall(MatMult(Op,x,y));
59 84725 : PetscFunctionReturn(PETSC_SUCCESS);
60 : }
61 :
62 4986 : PetscErrorCode STApplyMat_Generic(ST st,Mat B,Mat C)
63 : {
64 4986 : Mat work;
65 :
66 4986 : PetscFunctionBegin;
67 4986 : if (st->M && st->P) {
68 1936 : PetscCall(MatMatMult(st->M,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&work));
69 1936 : PetscCall(STMatMatSolve(st,work,C));
70 1936 : PetscCall(MatDestroy(&work));
71 3050 : } else if (st->M) PetscCall(MatMatMult(st->M,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C));
72 1704 : else PetscCall(STMatMatSolve(st,B,C));
73 4986 : PetscFunctionReturn(PETSC_SUCCESS);
74 : }
75 :
76 : /*@
77 : STApplyMat - Applies the spectral transformation operator to a matrix, for
78 : instance (A - sB)^-1 B in the case of the shift-and-invert transformation
79 : and generalized eigenproblem.
80 :
81 : Collective
82 :
83 : Input Parameters:
84 : + st - the spectral transformation context
85 : - X - input matrix
86 :
87 : Output Parameter:
88 : . Y - output matrix
89 :
90 : Level: developer
91 :
92 : .seealso: STApply()
93 : @*/
94 1658 : PetscErrorCode STApplyMat(ST st,Mat X,Mat Y)
95 : {
96 1658 : PetscFunctionBegin;
97 1658 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
98 1658 : PetscValidHeaderSpecific(X,MAT_CLASSID,2);
99 1658 : PetscValidHeaderSpecific(Y,MAT_CLASSID,3);
100 1658 : PetscValidType(st,1);
101 1658 : STCheckMatrices(st,1);
102 1658 : PetscCheck(X!=Y,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"X and Y must be different matrices");
103 1658 : PetscUseTypeMethod(st,applymat,X,Y);
104 1658 : PetscFunctionReturn(PETSC_SUCCESS);
105 : }
106 :
107 1244 : PetscErrorCode STApplyTranspose_Generic(ST st,Vec x,Vec y)
108 : {
109 1244 : PetscFunctionBegin;
110 1244 : if (st->M && st->P) {
111 37 : PetscCall(STMatSolveTranspose(st,x,st->work[0]));
112 37 : PetscCall(MatMultTranspose(st->M,st->work[0],y));
113 1207 : } else if (st->M) PetscCall(MatMultTranspose(st->M,x,y));
114 510 : else PetscCall(STMatSolveTranspose(st,x,y));
115 1244 : PetscFunctionReturn(PETSC_SUCCESS);
116 : }
117 :
118 : /*@
119 : STApplyTranspose - Applies the transpose of the operator to a vector, for
120 : instance B^T(A - sB)^-T in the case of the shift-and-invert transformation
121 : and generalized eigenproblem.
122 :
123 : Collective
124 :
125 : Input Parameters:
126 : + st - the spectral transformation context
127 : - x - input vector
128 :
129 : Output Parameter:
130 : . y - output vector
131 :
132 : Level: developer
133 :
134 : .seealso: STApply(), STApplyHermitianTranspose()
135 : @*/
136 6 : PetscErrorCode STApplyTranspose(ST st,Vec x,Vec y)
137 : {
138 6 : Mat Op;
139 :
140 6 : PetscFunctionBegin;
141 6 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
142 6 : PetscValidHeaderSpecific(x,VEC_CLASSID,2);
143 6 : PetscValidHeaderSpecific(y,VEC_CLASSID,3);
144 6 : PetscValidType(st,1);
145 6 : STCheckMatrices(st,1);
146 6 : PetscCheck(x!=y,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
147 6 : PetscCall(VecSetErrorIfLocked(y,3));
148 6 : PetscCall(STGetOperator_Private(st,&Op));
149 6 : PetscCall(MatMultTranspose(Op,x,y));
150 6 : PetscFunctionReturn(PETSC_SUCCESS);
151 : }
152 :
153 0 : PetscErrorCode STApplyHermitianTranspose_Generic(ST st,Vec x,Vec y)
154 : {
155 0 : PetscFunctionBegin;
156 0 : if (st->M && st->P) {
157 0 : PetscCall(STMatSolveHermitianTranspose(st,x,st->work[0]));
158 0 : PetscCall(MatMultHermitianTranspose(st->M,st->work[0],y));
159 0 : } else if (st->M) PetscCall(MatMultHermitianTranspose(st->M,x,y));
160 0 : else PetscCall(STMatSolveHermitianTranspose(st,x,y));
161 0 : PetscFunctionReturn(PETSC_SUCCESS);
162 : }
163 :
164 : /*@
165 : STApplyHermitianTranspose - Applies the hermitian-transpose of the operator
166 : to a vector, for instance B^H(A - sB)^-H in the case of the shift-and-invert
167 : transformation and generalized eigenproblem.
168 :
169 : Collective
170 :
171 : Input Parameters:
172 : + st - the spectral transformation context
173 : - x - input vector
174 :
175 : Output Parameter:
176 : . y - output vector
177 :
178 : Note:
179 : Currently implemented via STApplyTranspose() with appropriate conjugation.
180 :
181 : Level: developer
182 :
183 : .seealso: STApply(), STApplyTranspose()
184 : @*/
185 463 : PetscErrorCode STApplyHermitianTranspose(ST st,Vec x,Vec y)
186 : {
187 463 : Mat Op;
188 :
189 463 : PetscFunctionBegin;
190 463 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
191 463 : PetscValidHeaderSpecific(x,VEC_CLASSID,2);
192 463 : PetscValidHeaderSpecific(y,VEC_CLASSID,3);
193 463 : PetscValidType(st,1);
194 463 : STCheckMatrices(st,1);
195 463 : PetscCheck(x!=y,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
196 463 : PetscCall(VecSetErrorIfLocked(y,3));
197 463 : PetscCall(STGetOperator_Private(st,&Op));
198 463 : PetscCall(MatMultHermitianTranspose(Op,x,y));
199 463 : PetscFunctionReturn(PETSC_SUCCESS);
200 : }
201 :
202 : /*@
203 : STGetBilinearForm - Returns the matrix used in the bilinear form with a
204 : generalized problem with semi-definite B.
205 :
206 : Logically Collective
207 :
208 : Input Parameters:
209 : . st - the spectral transformation context
210 :
211 : Output Parameter:
212 : . B - output matrix
213 :
214 : Notes:
215 : The output matrix B must be destroyed after use. It will be NULL in
216 : case of standard eigenproblems.
217 :
218 : Level: developer
219 :
220 : .seealso: BVSetMatrix()
221 : @*/
222 154 : PetscErrorCode STGetBilinearForm(ST st,Mat *B)
223 : {
224 154 : PetscFunctionBegin;
225 154 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
226 154 : PetscValidType(st,1);
227 154 : PetscAssertPointer(B,2);
228 154 : STCheckMatrices(st,1);
229 154 : PetscUseTypeMethod(st,getbilinearform,B);
230 154 : PetscFunctionReturn(PETSC_SUCCESS);
231 : }
232 :
233 147 : PetscErrorCode STGetBilinearForm_Default(ST st,Mat *B)
234 : {
235 147 : PetscFunctionBegin;
236 147 : if (st->nmat==1) *B = NULL;
237 : else {
238 147 : *B = st->A[1];
239 147 : PetscCall(PetscObjectReference((PetscObject)*B));
240 : }
241 147 : PetscFunctionReturn(PETSC_SUCCESS);
242 : }
243 :
244 147634 : static PetscErrorCode MatMult_STOperator(Mat Op,Vec x,Vec y)
245 : {
246 147634 : ST st;
247 :
248 147634 : PetscFunctionBegin;
249 147634 : PetscCall(MatShellGetContext(Op,&st));
250 147634 : PetscCall(STSetUp(st));
251 147634 : PetscCall(PetscLogEventBegin(ST_Apply,st,x,y,0));
252 147634 : if (st->D) { /* with balancing */
253 697 : PetscCall(VecPointwiseDivide(st->wb,x,st->D));
254 697 : PetscUseTypeMethod(st,apply,st->wb,y);
255 697 : PetscCall(VecPointwiseMult(y,y,st->D));
256 146937 : } else PetscUseTypeMethod(st,apply,x,y);
257 147634 : PetscCall(PetscLogEventEnd(ST_Apply,st,x,y,0));
258 147634 : PetscFunctionReturn(PETSC_SUCCESS);
259 : }
260 :
261 1320 : static PetscErrorCode MatMultTranspose_STOperator(Mat Op,Vec x,Vec y)
262 : {
263 1320 : ST st;
264 :
265 1320 : PetscFunctionBegin;
266 1320 : PetscCall(MatShellGetContext(Op,&st));
267 1320 : PetscCall(STSetUp(st));
268 1320 : PetscCall(PetscLogEventBegin(ST_ApplyTranspose,st,x,y,0));
269 1320 : if (st->D) { /* with balancing */
270 50 : PetscCall(VecPointwiseMult(st->wb,x,st->D));
271 50 : PetscUseTypeMethod(st,applytrans,st->wb,y);
272 50 : PetscCall(VecPointwiseDivide(y,y,st->D));
273 1270 : } else PetscUseTypeMethod(st,applytrans,x,y);
274 1320 : PetscCall(PetscLogEventEnd(ST_ApplyTranspose,st,x,y,0));
275 1320 : PetscFunctionReturn(PETSC_SUCCESS);
276 : }
277 :
278 : #if defined(PETSC_USE_COMPLEX)
279 : static PetscErrorCode MatMultHermitianTranspose_STOperator(Mat Op,Vec x,Vec y)
280 : {
281 : ST st;
282 :
283 : PetscFunctionBegin;
284 : PetscCall(MatShellGetContext(Op,&st));
285 : PetscCall(STSetUp(st));
286 : PetscCall(PetscLogEventBegin(ST_ApplyHermitianTranspose,st,x,y,0));
287 : if (st->ops->applyhermtrans) {
288 : if (st->D) { /* with balancing */
289 : PetscCall(VecConjugate(st->D));
290 : PetscCall(VecPointwiseMult(st->wb,x,st->D));
291 : PetscUseTypeMethod(st,applyhermtrans,st->wb,y);
292 : PetscCall(VecPointwiseDivide(y,y,st->D));
293 : PetscCall(VecConjugate(st->D));
294 : } else PetscUseTypeMethod(st,applyhermtrans,x,y);
295 : } else {
296 : if (!st->wht) PetscCall(MatCreateVecs(st->A[0],&st->wht,NULL));
297 : PetscCall(VecCopy(x,st->wht));
298 : PetscCall(VecConjugate(st->wht));
299 : if (st->D) { /* with balancing */
300 : PetscCall(VecPointwiseMult(st->wb,st->wht,st->D));
301 : PetscUseTypeMethod(st,applytrans,st->wb,y);
302 : PetscCall(VecPointwiseDivide(y,y,st->D));
303 : } else PetscUseTypeMethod(st,applytrans,st->wht,y);
304 : PetscCall(VecConjugate(y));
305 : }
306 : PetscCall(PetscLogEventEnd(ST_ApplyHermitianTranspose,st,x,y,0));
307 : PetscFunctionReturn(PETSC_SUCCESS);
308 : }
309 : #endif
310 :
311 3328 : static PetscErrorCode MatMatMult_STOperator(Mat Op,Mat B,Mat C,void *ctx)
312 : {
313 3328 : ST st;
314 :
315 3328 : PetscFunctionBegin;
316 3328 : PetscCall(MatShellGetContext(Op,&st));
317 3328 : PetscCall(STSetUp(st));
318 3328 : PetscCall(PetscLogEventBegin(ST_Apply,st,B,C,0));
319 3328 : PetscCall(STApplyMat_Generic(st,B,C));
320 3328 : PetscCall(PetscLogEventEnd(ST_Apply,st,B,C,0));
321 3328 : PetscFunctionReturn(PETSC_SUCCESS);
322 : }
323 :
324 90039 : PetscErrorCode STGetOperator_Private(ST st,Mat *Op)
325 : {
326 90039 : PetscInt m,n,M,N;
327 90039 : Vec v;
328 90039 : VecType vtype;
329 :
330 90039 : PetscFunctionBegin;
331 90039 : if (!st->Op) {
332 508 : if (Op) *Op = NULL;
333 : /* create the shell matrix */
334 508 : PetscCall(MatGetLocalSize(st->A[0],&m,&n));
335 508 : PetscCall(MatGetSize(st->A[0],&M,&N));
336 508 : PetscCall(MatCreateShell(PetscObjectComm((PetscObject)st),m,n,M,N,st,&st->Op));
337 508 : PetscCall(MatShellSetOperation(st->Op,MATOP_MULT,(void(*)(void))MatMult_STOperator));
338 508 : PetscCall(MatShellSetOperation(st->Op,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_STOperator));
339 : #if defined(PETSC_USE_COMPLEX)
340 : PetscCall(MatShellSetOperation(st->Op,MATOP_MULT_HERMITIAN_TRANSPOSE,(void(*)(void))MatMultHermitianTranspose_STOperator));
341 : #else
342 508 : PetscCall(MatShellSetOperation(st->Op,MATOP_MULT_HERMITIAN_TRANSPOSE,(void(*)(void))MatMultTranspose_STOperator));
343 : #endif
344 508 : if (!st->D && st->ops->apply==STApply_Generic) {
345 466 : PetscCall(MatShellSetMatProductOperation(st->Op,MATPRODUCT_AB,NULL,MatMatMult_STOperator,NULL,MATDENSE,MATDENSE));
346 466 : PetscCall(MatShellSetMatProductOperation(st->Op,MATPRODUCT_AB,NULL,MatMatMult_STOperator,NULL,MATDENSECUDA,MATDENSECUDA));
347 466 : PetscCall(MatShellSetMatProductOperation(st->Op,MATPRODUCT_AB,NULL,MatMatMult_STOperator,NULL,MATDENSEHIP,MATDENSEHIP));
348 : }
349 : /* make sure the shell matrix generates a vector of the same type as the problem matrices */
350 508 : PetscCall(MatCreateVecs(st->A[0],&v,NULL));
351 508 : PetscCall(VecGetType(v,&vtype));
352 508 : PetscCall(MatShellSetVecType(st->Op,vtype));
353 508 : PetscCall(VecDestroy(&v));
354 : /* build the operator matrices */
355 508 : PetscCall(STComputeOperator(st));
356 : }
357 90039 : if (Op) *Op = st->Op;
358 90039 : PetscFunctionReturn(PETSC_SUCCESS);
359 : }
360 :
361 : /*@
362 : STGetOperator - Returns a shell matrix that represents the operator of the
363 : spectral transformation.
364 :
365 : Collective
366 :
367 : Input Parameter:
368 : . st - the spectral transformation context
369 :
370 : Output Parameter:
371 : . Op - operator matrix
372 :
373 : Notes:
374 : The operator is defined in linear eigenproblems only, not in polynomial ones,
375 : so the call will fail if more than 2 matrices were passed in STSetMatrices().
376 :
377 : The returned shell matrix is essentially a wrapper to the STApply() and
378 : STApplyTranspose() operations. The operator can often be expressed as
379 :
380 : $ Op = D*inv(K)*M*inv(D)
381 :
382 : where D is the balancing matrix, and M and K are two matrices corresponding
383 : to the numerator and denominator for spectral transformations that represent
384 : a rational matrix function. In the case of STSHELL, the inner part inv(K)*M
385 : is replaced by the user-provided operation from STShellSetApply().
386 :
387 : The preconditioner matrix K typically depends on the value of the shift, and
388 : its inverse is handled via an internal KSP object. Normal usage does not
389 : require explicitly calling STGetOperator(), but it can be used to force the
390 : creation of K and M, and then K is passed to the KSP. This is useful for
391 : setting options associated with the PCFactor (to set MUMPS options, for instance).
392 :
393 : The returned matrix must NOT be destroyed by the user. Instead, when no
394 : longer needed it must be returned with STRestoreOperator(). In particular,
395 : this is required before modifying the ST matrices or the shift.
396 :
397 : A NULL pointer can be passed in Op in case the matrix is not required but we
398 : want to force its creation. In this case, STRestoreOperator() should not be
399 : called.
400 :
401 : Level: advanced
402 :
403 : .seealso: STApply(), STApplyTranspose(), STSetBalanceMatrix(), STShellSetApply(),
404 : STGetKSP(), STSetShift(), STRestoreOperator(), STSetMatrices()
405 : @*/
406 4845 : PetscErrorCode STGetOperator(ST st,Mat *Op)
407 : {
408 4845 : PetscFunctionBegin;
409 4845 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
410 4845 : PetscValidType(st,1);
411 4845 : STCheckMatrices(st,1);
412 4845 : STCheckNotSeized(st,1);
413 4845 : PetscCheck(st->nmat<=2,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_WRONGSTATE,"The operator is not defined in polynomial eigenproblems");
414 4845 : PetscCall(STGetOperator_Private(st,Op));
415 4845 : if (Op) st->opseized = PETSC_TRUE;
416 4845 : PetscFunctionReturn(PETSC_SUCCESS);
417 : }
418 :
419 : /*@
420 : STRestoreOperator - Restore the previously seized operator matrix.
421 :
422 : Logically Collective
423 :
424 : Input Parameters:
425 : + st - the spectral transformation context
426 : - Op - operator matrix
427 :
428 : Notes:
429 : The arguments must match the corresponding call to STGetOperator().
430 :
431 : Level: advanced
432 :
433 : .seealso: STGetOperator()
434 : @*/
435 4815 : PetscErrorCode STRestoreOperator(ST st,Mat *Op)
436 : {
437 4815 : PetscFunctionBegin;
438 4815 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
439 4815 : PetscAssertPointer(Op,2);
440 4815 : PetscValidHeaderSpecific(*Op,MAT_CLASSID,2);
441 4815 : PetscCheck(st->opseized,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_WRONGSTATE,"Must be called after STGetOperator()");
442 4815 : *Op = NULL;
443 4815 : st->opseized = PETSC_FALSE;
444 4815 : PetscFunctionReturn(PETSC_SUCCESS);
445 : }
446 :
447 : /*
448 : STComputeOperator - Computes the matrices that constitute the operator
449 :
450 : Op = D*inv(K)*M*inv(D).
451 :
452 : K and M are computed here (D is user-provided) from the system matrices
453 : and the shift sigma (whenever these are changed, this function recomputes
454 : K and M). This is used only in linear eigenproblems (nmat<3).
455 :
456 : K is the "preconditioner matrix": it is the denominator in rational operators,
457 : e.g. (A-sigma*B) in shift-and-invert. In non-rational transformations such
458 : as STFILTER, K=NULL which means identity. After computing K, it is passed to
459 : the internal KSP object via KSPSetOperators.
460 :
461 : M is the numerator in rational operators. If unused it is set to NULL (e.g.
462 : in STPRECOND).
463 :
464 : STSHELL does not compute anything here, but sets the flag as if it was ready.
465 : */
466 1307 : PetscErrorCode STComputeOperator(ST st)
467 : {
468 1307 : PC pc;
469 :
470 1307 : PetscFunctionBegin;
471 1307 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
472 1307 : PetscValidType(st,1);
473 1307 : if (!st->opready && st->ops->computeoperator) {
474 780 : PetscCall(PetscInfo(st,"Building the operator matrices\n"));
475 780 : STCheckMatrices(st,1);
476 780 : if (!st->T) PetscCall(PetscCalloc1(PetscMax(2,st->nmat),&st->T));
477 780 : PetscCall(PetscLogEventBegin(ST_ComputeOperator,st,0,0,0));
478 780 : PetscUseTypeMethod(st,computeoperator);
479 780 : PetscCall(PetscLogEventEnd(ST_ComputeOperator,st,0,0,0));
480 780 : if (st->usesksp) {
481 421 : if (!st->ksp) PetscCall(STGetKSP(st,&st->ksp));
482 421 : if (st->P) {
483 384 : PetscCall(STSetDefaultKSP(st));
484 384 : PetscCall(ST_KSPSetOperators(st,st->P,st->Pmat?st->Pmat:st->P));
485 : } else {
486 : /* STPRECOND defaults to PCNONE if st->P is empty */
487 37 : PetscCall(KSPGetPC(st->ksp,&pc));
488 37 : PetscCall(PCSetType(pc,PCNONE));
489 : }
490 : }
491 : }
492 1307 : st->opready = PETSC_TRUE;
493 1307 : PetscFunctionReturn(PETSC_SUCCESS);
494 : }
495 :
496 : /*@
497 : STSetUp - Prepares for the use of a spectral transformation.
498 :
499 : Collective
500 :
501 : Input Parameter:
502 : . st - the spectral transformation context
503 :
504 : Level: advanced
505 :
506 : .seealso: STCreate(), STApply(), STDestroy()
507 : @*/
508 153538 : PetscErrorCode STSetUp(ST st)
509 : {
510 153538 : PetscInt i,n,k;
511 :
512 153538 : PetscFunctionBegin;
513 153538 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
514 153538 : PetscValidType(st,1);
515 153538 : STCheckMatrices(st,1);
516 153538 : switch (st->state) {
517 967 : case ST_STATE_INITIAL:
518 967 : PetscCall(PetscInfo(st,"Setting up new ST\n"));
519 967 : if (!((PetscObject)st)->type_name) PetscCall(STSetType(st,STSHIFT));
520 : break;
521 152537 : case ST_STATE_SETUP:
522 152537 : PetscFunctionReturn(PETSC_SUCCESS);
523 34 : case ST_STATE_UPDATED:
524 34 : PetscCall(PetscInfo(st,"Setting up updated ST\n"));
525 : break;
526 : }
527 1001 : PetscCall(PetscLogEventBegin(ST_SetUp,st,0,0,0));
528 1001 : if (st->state!=ST_STATE_UPDATED) {
529 967 : if (!(st->nmat<3 && st->opready)) {
530 942 : if (st->T) {
531 249 : for (i=0;i<PetscMax(2,st->nmat);i++) PetscCall(MatDestroy(&st->T[i]));
532 : }
533 942 : PetscCall(MatDestroy(&st->P));
534 : }
535 : }
536 1001 : if (st->D) {
537 13 : PetscCall(MatGetLocalSize(st->A[0],NULL,&n));
538 13 : PetscCall(VecGetLocalSize(st->D,&k));
539 13 : PetscCheck(n==k,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_SIZ,"Balance matrix has wrong dimension %" PetscInt_FMT " (should be %" PetscInt_FMT ")",k,n);
540 13 : if (!st->wb) PetscCall(VecDuplicate(st->D,&st->wb));
541 : }
542 1001 : if (st->nmat<3 && st->transform) PetscCall(STComputeOperator(st));
543 : else {
544 202 : if (!st->T) PetscCall(PetscCalloc1(PetscMax(2,st->nmat),&st->T));
545 : }
546 1001 : PetscTryTypeMethod(st,setup);
547 1001 : st->state = ST_STATE_SETUP;
548 1001 : PetscCall(PetscLogEventEnd(ST_SetUp,st,0,0,0));
549 1001 : PetscFunctionReturn(PETSC_SUCCESS);
550 : }
551 :
552 : /*
553 : Computes coefficients for the transformed polynomial,
554 : and stores the result in argument S.
555 :
556 : alpha - value of the parameter of the transformed polynomial
557 : beta - value of the previous shift (only used in inplace mode)
558 : k - index of first matrix included in the computation
559 : coeffs - coefficients of the expansion
560 : initial - true if this is the first time
561 : precond - whether the preconditioner matrix must be computed
562 : */
563 1415 : PetscErrorCode STMatMAXPY_Private(ST st,PetscScalar alpha,PetscScalar beta,PetscInt k,PetscScalar *coeffs,PetscBool initial,PetscBool precond,Mat *S)
564 : {
565 1415 : PetscInt *matIdx=NULL,nmat,i,ini=-1;
566 1415 : PetscScalar t=1.0,ta,gamma;
567 1415 : PetscBool nz=PETSC_FALSE;
568 1415 : Mat *A=precond?st->Psplit:st->A;
569 1415 : MatStructure str=precond?st->strp:st->str;
570 :
571 1415 : PetscFunctionBegin;
572 1415 : nmat = st->nmat-k;
573 1415 : switch (st->matmode) {
574 20 : case ST_MATMODE_INPLACE:
575 20 : PetscCheck(st->nmat<=2,PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"ST_MATMODE_INPLACE not supported for polynomial eigenproblems");
576 20 : PetscCheck(!precond,PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"ST_MATMODE_INPLACE not supported for split preconditioner");
577 20 : if (initial) {
578 11 : PetscCall(PetscObjectReference((PetscObject)A[0]));
579 11 : *S = A[0];
580 11 : gamma = alpha;
581 9 : } else gamma = alpha-beta;
582 20 : if (gamma != 0.0) {
583 16 : if (st->nmat>1) PetscCall(MatAXPY(*S,gamma,A[1],str));
584 11 : else PetscCall(MatShift(*S,gamma));
585 : }
586 : break;
587 101 : case ST_MATMODE_SHELL:
588 101 : PetscCheck(!precond,PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"ST_MATMODE_SHELL not supported for split preconditioner");
589 101 : if (initial) {
590 82 : if (st->nmat>2) {
591 10 : PetscCall(PetscMalloc1(nmat,&matIdx));
592 39 : for (i=0;i<nmat;i++) matIdx[i] = k+i;
593 : }
594 82 : PetscCall(STMatShellCreate(st,alpha,nmat,matIdx,coeffs,S));
595 82 : if (st->nmat>2) PetscCall(PetscFree(matIdx));
596 19 : } else PetscCall(STMatShellShift(*S,alpha));
597 : break;
598 1294 : case ST_MATMODE_COPY:
599 1294 : if (coeffs) {
600 546 : for (i=0;i<nmat && ini==-1;i++) {
601 318 : if (coeffs[i]!=0.0) ini = i;
602 90 : else t *= alpha;
603 : }
604 228 : if (coeffs[ini] != 1.0) nz = PETSC_TRUE;
605 430 : for (i=ini+1;i<nmat&&!nz;i++) if (coeffs[i]!=0.0) nz = PETSC_TRUE;
606 : } else { nz = PETSC_TRUE; ini = 0; }
607 1294 : if ((alpha == 0.0 || !nz) && t==1.0) {
608 650 : PetscCall(PetscObjectReference((PetscObject)A[k+ini]));
609 650 : PetscCall(MatDestroy(S));
610 650 : *S = A[k+ini];
611 : } else {
612 644 : if (*S && *S!=A[k+ini]) {
613 427 : PetscCall(MatSetOption(*S,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE));
614 427 : PetscCall(MatCopy(A[k+ini],*S,DIFFERENT_NONZERO_PATTERN));
615 : } else {
616 217 : PetscCall(MatDestroy(S));
617 217 : PetscCall(MatDuplicate(A[k+ini],MAT_COPY_VALUES,S));
618 217 : PetscCall(MatSetOption(*S,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE));
619 : }
620 644 : if (coeffs && coeffs[ini]!=1.0) PetscCall(MatScale(*S,coeffs[ini]));
621 1535 : for (i=ini+k+1;i<PetscMax(2,st->nmat);i++) {
622 891 : t *= alpha;
623 891 : ta = t;
624 891 : if (coeffs) ta *= coeffs[i-k];
625 891 : if (ta!=0.0) {
626 890 : if (st->nmat>1) PetscCall(MatAXPY(*S,ta,A[i],str));
627 891 : else PetscCall(MatShift(*S,ta));
628 : }
629 : }
630 : }
631 : }
632 1415 : PetscCall(MatSetOption(*S,MAT_SYMMETRIC,st->asymm));
633 1415 : PetscCall(MatSetOption(*S,MAT_HERMITIAN,(PetscImaginaryPart(st->sigma)==0.0)?st->aherm:PETSC_FALSE));
634 1415 : PetscFunctionReturn(PETSC_SUCCESS);
635 : }
636 :
637 : /*
638 : Computes the values of the coefficients required by STMatMAXPY_Private
639 : for the case of monomial basis.
640 : */
641 40 : PetscErrorCode STCoeffs_Monomial(ST st, PetscScalar *coeffs)
642 : {
643 40 : PetscInt k,i,ini,inip;
644 :
645 40 : PetscFunctionBegin;
646 : /* Compute binomial coefficients */
647 40 : ini = (st->nmat*(st->nmat-1))/2;
648 170 : for (i=0;i<st->nmat;i++) coeffs[ini+i]=1.0;
649 130 : for (k=st->nmat-1;k>=1;k--) {
650 90 : inip = ini+1;
651 90 : ini = (k*(k-1))/2;
652 90 : coeffs[ini] = 1.0;
653 150 : for (i=1;i<k;i++) coeffs[ini+i] = coeffs[ini+i-1]+coeffs[inip+i-1];
654 : }
655 40 : PetscFunctionReturn(PETSC_SUCCESS);
656 : }
657 :
658 : /*@
659 : STPostSolve - Optional post-solve phase, intended for any actions that must
660 : be performed on the ST object after the eigensolver has finished.
661 :
662 : Collective
663 :
664 : Input Parameters:
665 : . st - the spectral transformation context
666 :
667 : Level: developer
668 :
669 : .seealso: EPSSolve()
670 : @*/
671 1032 : PetscErrorCode STPostSolve(ST st)
672 : {
673 1032 : PetscFunctionBegin;
674 1032 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
675 1032 : PetscValidType(st,1);
676 1032 : PetscTryTypeMethod(st,postsolve);
677 1032 : PetscFunctionReturn(PETSC_SUCCESS);
678 : }
679 :
680 : /*@
681 : STBackTransform - Back-transformation phase, intended for
682 : spectral transformations which require to transform the computed
683 : eigenvalues back to the original eigenvalue problem.
684 :
685 : Not Collective
686 :
687 : Input Parameters:
688 : + st - the spectral transformation context
689 : . n - number of eigenvalues
690 : . eigr - real part of a computed eigenvalues
691 : - eigi - imaginary part of a computed eigenvalues
692 :
693 : Level: developer
694 :
695 : .seealso: STIsInjective()
696 : @*/
697 2170843 : PetscErrorCode STBackTransform(ST st,PetscInt n,PetscScalar* eigr,PetscScalar* eigi)
698 : {
699 2170843 : PetscFunctionBegin;
700 2170843 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
701 2170843 : PetscValidType(st,1);
702 2170843 : PetscTryTypeMethod(st,backtransform,n,eigr,eigi);
703 2170843 : PetscFunctionReturn(PETSC_SUCCESS);
704 : }
705 :
706 : /*@
707 : STIsInjective - Ask if this spectral transformation is injective or not
708 : (that is, if it corresponds to a one-to-one mapping). If not, then it
709 : does not make sense to call STBackTransform().
710 :
711 : Not Collective
712 :
713 : Input Parameter:
714 : . st - the spectral transformation context
715 :
716 : Output Parameter:
717 : . is - the answer
718 :
719 : Level: developer
720 :
721 : .seealso: STBackTransform()
722 : @*/
723 653 : PetscErrorCode STIsInjective(ST st,PetscBool* is)
724 : {
725 653 : PetscBool shell;
726 :
727 653 : PetscFunctionBegin;
728 653 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
729 653 : PetscValidType(st,1);
730 653 : PetscAssertPointer(is,2);
731 :
732 653 : PetscCall(PetscObjectTypeCompare((PetscObject)st,STSHELL,&shell));
733 653 : if (shell) PetscCall(STIsInjective_Shell(st,is));
734 622 : else *is = st->ops->backtransform? PETSC_TRUE: PETSC_FALSE;
735 653 : PetscFunctionReturn(PETSC_SUCCESS);
736 : }
737 :
738 : /*@
739 : STMatSetUp - Build the preconditioner matrix used in STMatSolve().
740 :
741 : Collective
742 :
743 : Input Parameters:
744 : + st - the spectral transformation context
745 : . sigma - the shift
746 : - coeffs - the coefficients (may be NULL)
747 :
748 : Note:
749 : This function is not intended to be called by end users, but by SLEPc
750 : solvers that use ST. It builds matrix st->P as follows, then calls KSPSetUp().
751 : .vb
752 : If (coeffs) st->P = Sum_{i=0..nmat-1} coeffs[i]*sigma^i*A_i
753 : else st->P = Sum_{i=0..nmat-1} sigma^i*A_i
754 : .ve
755 :
756 : Level: developer
757 :
758 : .seealso: STMatSolve()
759 : @*/
760 145 : PetscErrorCode STMatSetUp(ST st,PetscScalar sigma,PetscScalar *coeffs)
761 : {
762 145 : PetscFunctionBegin;
763 145 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
764 435 : PetscValidLogicalCollectiveScalar(st,sigma,2);
765 145 : STCheckMatrices(st,1);
766 :
767 145 : PetscCall(PetscLogEventBegin(ST_MatSetUp,st,0,0,0));
768 145 : PetscCall(STMatMAXPY_Private(st,sigma,0.0,0,coeffs,PETSC_TRUE,PETSC_FALSE,&st->P));
769 145 : if (st->Psplit) PetscCall(STMatMAXPY_Private(st,sigma,0.0,0,coeffs,PETSC_TRUE,PETSC_TRUE,&st->Pmat));
770 145 : PetscCall(ST_KSPSetOperators(st,st->P,st->Pmat?st->Pmat:st->P));
771 145 : PetscCall(KSPSetUp(st->ksp));
772 145 : PetscCall(PetscLogEventEnd(ST_MatSetUp,st,0,0,0));
773 145 : PetscFunctionReturn(PETSC_SUCCESS);
774 : }
775 :
776 : /*@
777 : STSetWorkVecs - Sets a number of work vectors into the ST object.
778 :
779 : Collective
780 :
781 : Input Parameters:
782 : + st - the spectral transformation context
783 : - nw - number of work vectors to allocate
784 :
785 : Developer Notes:
786 : This is SLEPC_EXTERN because it may be required by shell STs.
787 :
788 : Level: developer
789 :
790 : .seealso: STMatCreateVecs()
791 : @*/
792 362 : PetscErrorCode STSetWorkVecs(ST st,PetscInt nw)
793 : {
794 362 : PetscInt i;
795 :
796 362 : PetscFunctionBegin;
797 362 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
798 1086 : PetscValidLogicalCollectiveInt(st,nw,2);
799 362 : PetscCheck(nw>0,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_OUTOFRANGE,"nw must be > 0: nw = %" PetscInt_FMT,nw);
800 362 : if (st->nwork < nw) {
801 342 : PetscCall(VecDestroyVecs(st->nwork,&st->work));
802 342 : st->nwork = nw;
803 342 : PetscCall(PetscMalloc1(nw,&st->work));
804 715 : for (i=0;i<nw;i++) PetscCall(STMatCreateVecs(st,&st->work[i],NULL));
805 : }
806 362 : PetscFunctionReturn(PETSC_SUCCESS);
807 : }
|
