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 97206 : PetscErrorCode STApply_Generic(ST st,Vec x,Vec y)
17 : {
18 97206 : PetscFunctionBegin;
19 97206 : if (st->M && st->P) {
20 18338 : PetscCall(MatMult(st->M,x,st->work[0]));
21 18338 : PetscCall(STMatSolve(st,st->work[0],y));
22 78868 : } else if (st->M) PetscCall(MatMult(st->M,x,y));
23 16832 : else PetscCall(STMatSolve(st,x,y));
24 97206 : 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 46551 : PetscErrorCode STApply(ST st,Vec x,Vec y)
46 : {
47 46551 : Mat Op;
48 :
49 46551 : PetscFunctionBegin;
50 46551 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
51 46551 : PetscValidHeaderSpecific(x,VEC_CLASSID,2);
52 46551 : PetscValidHeaderSpecific(y,VEC_CLASSID,3);
53 46551 : PetscValidType(st,1);
54 46551 : STCheckMatrices(st,1);
55 46551 : PetscCheck(x!=y,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
56 46551 : PetscCall(VecSetErrorIfLocked(y,3));
57 46551 : PetscCall(STGetOperator_Private(st,&Op));
58 46551 : PetscCall(MatMult(Op,x,y));
59 46551 : PetscFunctionReturn(PETSC_SUCCESS);
60 : }
61 :
62 5053 : PetscErrorCode STApplyMat_Generic(ST st,Mat B,Mat C)
63 : {
64 5053 : Mat work;
65 :
66 5053 : PetscFunctionBegin;
67 5053 : if (st->M && st->P) {
68 1999 : PetscCall(MatMatMult(st->M,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&work));
69 1999 : PetscCall(STMatMatSolve(st,work,C));
70 1999 : PetscCall(MatDestroy(&work));
71 3054 : } else if (st->M) PetscCall(MatMatMult(st->M,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C));
72 1363 : else PetscCall(STMatMatSolve(st,B,C));
73 5053 : 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 1319 : PetscErrorCode STApplyMat(ST st,Mat X,Mat Y)
95 : {
96 1319 : PetscFunctionBegin;
97 1319 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
98 1319 : PetscValidHeaderSpecific(X,MAT_CLASSID,2);
99 1319 : PetscValidHeaderSpecific(Y,MAT_CLASSID,3);
100 1319 : PetscValidType(st,1);
101 1319 : STCheckMatrices(st,1);
102 1319 : PetscCheck(X!=Y,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"X and Y must be different matrices");
103 1319 : PetscUseTypeMethod(st,applymat,X,Y);
104 1319 : PetscFunctionReturn(PETSC_SUCCESS);
105 : }
106 :
107 18 : PetscErrorCode STApplyTranspose_Generic(ST st,Vec x,Vec y)
108 : {
109 18 : PetscFunctionBegin;
110 18 : if (st->M && st->P) {
111 12 : PetscCall(STMatSolveTranspose(st,x,st->work[0]));
112 12 : PetscCall(MatMultTranspose(st->M,st->work[0],y));
113 6 : } else if (st->M) PetscCall(MatMultTranspose(st->M,x,y));
114 0 : else PetscCall(STMatSolveTranspose(st,x,y));
115 18 : 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 12 : PetscErrorCode STApplyTranspose(ST st,Vec x,Vec y)
137 : {
138 12 : Mat Op;
139 :
140 12 : PetscFunctionBegin;
141 12 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
142 12 : PetscValidHeaderSpecific(x,VEC_CLASSID,2);
143 12 : PetscValidHeaderSpecific(y,VEC_CLASSID,3);
144 12 : PetscValidType(st,1);
145 12 : STCheckMatrices(st,1);
146 12 : PetscCheck(x!=y,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
147 12 : PetscCall(VecSetErrorIfLocked(y,3));
148 12 : PetscCall(STGetOperator_Private(st,&Op));
149 12 : PetscCall(MatMultTranspose(Op,x,y));
150 12 : PetscFunctionReturn(PETSC_SUCCESS);
151 : }
152 :
153 1285 : PetscErrorCode STApplyHermitianTranspose_Generic(ST st,Vec x,Vec y)
154 : {
155 1285 : PetscFunctionBegin;
156 1285 : if (st->M && st->P) {
157 57 : PetscCall(STMatSolveHermitianTranspose(st,x,st->work[0]));
158 57 : PetscCall(MatMultHermitianTranspose(st->M,st->work[0],y));
159 1228 : } else if (st->M) PetscCall(MatMultHermitianTranspose(st->M,x,y));
160 498 : else PetscCall(STMatSolveHermitianTranspose(st,x,y));
161 1285 : 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 445 : PetscErrorCode STApplyHermitianTranspose(ST st,Vec x,Vec y)
186 : {
187 445 : Mat Op;
188 :
189 445 : PetscFunctionBegin;
190 445 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
191 445 : PetscValidHeaderSpecific(x,VEC_CLASSID,2);
192 445 : PetscValidHeaderSpecific(y,VEC_CLASSID,3);
193 445 : PetscValidType(st,1);
194 445 : STCheckMatrices(st,1);
195 445 : PetscCheck(x!=y,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
196 445 : PetscCall(VecSetErrorIfLocked(y,3));
197 445 : PetscCall(STGetOperator_Private(st,&Op));
198 445 : PetscCall(MatMultHermitianTranspose(Op,x,y));
199 445 : 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 136 : PetscErrorCode STGetBilinearForm(ST st,Mat *B)
223 : {
224 136 : PetscFunctionBegin;
225 136 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
226 136 : PetscValidType(st,1);
227 136 : PetscAssertPointer(B,2);
228 136 : STCheckMatrices(st,1);
229 136 : PetscUseTypeMethod(st,getbilinearform,B);
230 136 : PetscFunctionReturn(PETSC_SUCCESS);
231 : }
232 :
233 126 : PetscErrorCode STGetBilinearForm_Default(ST st,Mat *B)
234 : {
235 126 : PetscFunctionBegin;
236 126 : if (st->nmat==1) *B = NULL;
237 : else {
238 126 : *B = st->A[1];
239 126 : PetscCall(PetscObjectReference((PetscObject)*B));
240 : }
241 126 : PetscFunctionReturn(PETSC_SUCCESS);
242 : }
243 :
244 101417 : static PetscErrorCode MatMult_STOperator(Mat Op,Vec x,Vec y)
245 : {
246 101417 : ST st;
247 :
248 101417 : PetscFunctionBegin;
249 101417 : PetscCall(MatShellGetContext(Op,&st));
250 101417 : PetscCall(STSetUp(st));
251 101417 : PetscCall(PetscLogEventBegin(ST_Apply,st,x,y,0));
252 101417 : if (st->D) { /* with balancing */
253 712 : PetscCall(VecPointwiseDivide(st->wb,x,st->D));
254 712 : PetscUseTypeMethod(st,apply,st->wb,y);
255 712 : PetscCall(VecPointwiseMult(y,y,st->D));
256 100705 : } else PetscUseTypeMethod(st,apply,x,y);
257 101417 : PetscCall(PetscLogEventEnd(ST_Apply,st,x,y,0));
258 101417 : PetscFunctionReturn(PETSC_SUCCESS);
259 : }
260 :
261 18 : static PetscErrorCode MatMultTranspose_STOperator(Mat Op,Vec x,Vec y)
262 : {
263 18 : ST st;
264 :
265 18 : PetscFunctionBegin;
266 18 : PetscCall(MatShellGetContext(Op,&st));
267 18 : PetscCall(STSetUp(st));
268 18 : PetscCall(PetscLogEventBegin(ST_ApplyTranspose,st,x,y,0));
269 18 : if (st->D) { /* with balancing */
270 0 : PetscCall(VecPointwiseMult(st->wb,x,st->D));
271 0 : PetscUseTypeMethod(st,applytrans,st->wb,y);
272 0 : PetscCall(VecPointwiseDivide(y,y,st->D));
273 18 : } else PetscUseTypeMethod(st,applytrans,x,y);
274 18 : PetscCall(PetscLogEventEnd(ST_ApplyTranspose,st,x,y,0));
275 18 : PetscFunctionReturn(PETSC_SUCCESS);
276 : }
277 :
278 : #if defined(PETSC_USE_COMPLEX)
279 1361 : static PetscErrorCode MatMultHermitianTranspose_STOperator(Mat Op,Vec x,Vec y)
280 : {
281 1361 : ST st;
282 :
283 1361 : PetscFunctionBegin;
284 1361 : PetscCall(MatShellGetContext(Op,&st));
285 1361 : PetscCall(STSetUp(st));
286 1361 : PetscCall(PetscLogEventBegin(ST_ApplyHermitianTranspose,st,x,y,0));
287 1361 : if (st->ops->applyhermtrans) {
288 1361 : if (st->D) { /* with balancing */
289 50 : PetscCall(VecConjugate(st->D));
290 50 : PetscCall(VecPointwiseMult(st->wb,x,st->D));
291 50 : PetscUseTypeMethod(st,applyhermtrans,st->wb,y);
292 50 : PetscCall(VecPointwiseDivide(y,y,st->D));
293 50 : PetscCall(VecConjugate(st->D));
294 1311 : } else PetscUseTypeMethod(st,applyhermtrans,x,y);
295 : } else {
296 0 : if (!st->wht) PetscCall(MatCreateVecs(st->A[0],&st->wht,NULL));
297 0 : PetscCall(VecCopy(x,st->wht));
298 0 : PetscCall(VecConjugate(st->wht));
299 0 : if (st->D) { /* with balancing */
300 0 : PetscCall(VecPointwiseMult(st->wb,st->wht,st->D));
301 0 : PetscUseTypeMethod(st,applytrans,st->wb,y);
302 0 : PetscCall(VecPointwiseDivide(y,y,st->D));
303 0 : } else PetscUseTypeMethod(st,applytrans,st->wht,y);
304 0 : PetscCall(VecConjugate(y));
305 : }
306 1361 : PetscCall(PetscLogEventEnd(ST_ApplyHermitianTranspose,st,x,y,0));
307 1361 : PetscFunctionReturn(PETSC_SUCCESS);
308 : }
309 : #endif
310 :
311 3734 : static PetscErrorCode MatMatMult_STOperator(Mat Op,Mat B,Mat C,void *ctx)
312 : {
313 3734 : ST st;
314 :
315 3734 : PetscFunctionBegin;
316 3734 : PetscCall(MatShellGetContext(Op,&st));
317 3734 : PetscCall(STSetUp(st));
318 3734 : PetscCall(PetscLogEventBegin(ST_Apply,st,B,C,0));
319 3734 : PetscCall(STApplyMat_Generic(st,B,C));
320 3734 : PetscCall(PetscLogEventEnd(ST_Apply,st,B,C,0));
321 3734 : PetscFunctionReturn(PETSC_SUCCESS);
322 : }
323 :
324 50960 : PetscErrorCode STGetOperator_Private(ST st,Mat *Op)
325 : {
326 50960 : PetscInt m,n,M,N;
327 50960 : Vec v;
328 50960 : VecType vtype;
329 :
330 50960 : PetscFunctionBegin;
331 50960 : if (!st->Op) {
332 478 : if (Op) *Op = NULL;
333 : /* create the shell matrix */
334 478 : PetscCall(MatGetLocalSize(st->A[0],&m,&n));
335 478 : PetscCall(MatGetSize(st->A[0],&M,&N));
336 478 : PetscCall(MatCreateShell(PetscObjectComm((PetscObject)st),m,n,M,N,st,&st->Op));
337 478 : PetscCall(MatShellSetOperation(st->Op,MATOP_MULT,(void(*)(void))MatMult_STOperator));
338 478 : PetscCall(MatShellSetOperation(st->Op,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_STOperator));
339 : #if defined(PETSC_USE_COMPLEX)
340 478 : PetscCall(MatShellSetOperation(st->Op,MATOP_MULT_HERMITIAN_TRANSPOSE,(void(*)(void))MatMultHermitianTranspose_STOperator));
341 : #else
342 : PetscCall(MatShellSetOperation(st->Op,MATOP_MULT_HERMITIAN_TRANSPOSE,(void(*)(void))MatMultTranspose_STOperator));
343 : #endif
344 478 : if (!st->D && st->ops->apply==STApply_Generic) {
345 435 : PetscCall(MatShellSetMatProductOperation(st->Op,MATPRODUCT_AB,NULL,MatMatMult_STOperator,NULL,MATDENSE,MATDENSE));
346 435 : PetscCall(MatShellSetMatProductOperation(st->Op,MATPRODUCT_AB,NULL,MatMatMult_STOperator,NULL,MATDENSECUDA,MATDENSECUDA));
347 435 : 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 478 : PetscCall(MatCreateVecs(st->A[0],&v,NULL));
351 478 : PetscCall(VecGetType(v,&vtype));
352 478 : PetscCall(MatShellSetVecType(st->Op,vtype));
353 478 : PetscCall(VecDestroy(&v));
354 : /* build the operator matrices */
355 478 : PetscCall(STComputeOperator(st));
356 : }
357 50960 : if (Op) *Op = st->Op;
358 50960 : 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 3952 : PetscErrorCode STGetOperator(ST st,Mat *Op)
407 : {
408 3952 : PetscFunctionBegin;
409 3952 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
410 3952 : PetscValidType(st,1);
411 3952 : STCheckMatrices(st,1);
412 3952 : STCheckNotSeized(st,1);
413 3952 : PetscCheck(st->nmat<=2,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_WRONGSTATE,"The operator is not defined in polynomial eigenproblems");
414 3952 : PetscCall(STGetOperator_Private(st,Op));
415 3952 : if (Op) st->opseized = PETSC_TRUE;
416 3952 : 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 3922 : PetscErrorCode STRestoreOperator(ST st,Mat *Op)
436 : {
437 3922 : PetscFunctionBegin;
438 3922 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
439 3922 : PetscAssertPointer(Op,2);
440 3922 : PetscValidHeaderSpecific(*Op,MAT_CLASSID,2);
441 3922 : PetscCheck(st->opseized,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_WRONGSTATE,"Must be called after STGetOperator()");
442 3922 : *Op = NULL;
443 3922 : st->opseized = PETSC_FALSE;
444 3922 : 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 1255 : PetscErrorCode STComputeOperator(ST st)
467 : {
468 1255 : PC pc;
469 :
470 1255 : PetscFunctionBegin;
471 1255 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
472 1255 : PetscValidType(st,1);
473 1255 : if (!st->opready && st->ops->computeoperator) {
474 757 : PetscCall(PetscInfo(st,"Building the operator matrices\n"));
475 757 : STCheckMatrices(st,1);
476 757 : if (!st->T) PetscCall(PetscCalloc1(PetscMax(2,st->nmat),&st->T));
477 757 : PetscCall(PetscLogEventBegin(ST_ComputeOperator,st,0,0,0));
478 757 : PetscUseTypeMethod(st,computeoperator);
479 757 : PetscCall(PetscLogEventEnd(ST_ComputeOperator,st,0,0,0));
480 757 : if (st->usesksp) {
481 414 : if (!st->ksp) PetscCall(STGetKSP(st,&st->ksp));
482 414 : if (st->P) {
483 375 : PetscCall(STSetDefaultKSP(st));
484 375 : 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 39 : PetscCall(KSPGetPC(st->ksp,&pc));
488 39 : PetscCall(PCSetType(pc,PCNONE));
489 : }
490 : }
491 : }
492 1255 : st->opready = PETSC_TRUE;
493 1255 : 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 107774 : PetscErrorCode STSetUp(ST st)
509 : {
510 107774 : PetscInt i,n,k;
511 :
512 107774 : PetscFunctionBegin;
513 107774 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
514 107774 : PetscValidType(st,1);
515 107774 : STCheckMatrices(st,1);
516 107774 : switch (st->state) {
517 957 : case ST_STATE_INITIAL:
518 957 : PetscCall(PetscInfo(st,"Setting up new ST\n"));
519 957 : if (!((PetscObject)st)->type_name) PetscCall(STSetType(st,STSHIFT));
520 : break;
521 106779 : case ST_STATE_SETUP:
522 106779 : PetscFunctionReturn(PETSC_SUCCESS);
523 38 : case ST_STATE_UPDATED:
524 38 : PetscCall(PetscInfo(st,"Setting up updated ST\n"));
525 : break;
526 : }
527 995 : PetscCall(PetscLogEventBegin(ST_SetUp,st,0,0,0));
528 995 : if (st->state!=ST_STATE_UPDATED) {
529 957 : if (!(st->nmat<3 && st->opready)) {
530 929 : if (st->T) {
531 276 : for (i=0;i<PetscMax(2,st->nmat);i++) PetscCall(MatDestroy(&st->T[i]));
532 : }
533 929 : PetscCall(MatDestroy(&st->P));
534 : }
535 : }
536 995 : 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 995 : if (st->nmat<3 && st->transform) PetscCall(STComputeOperator(st));
543 : else {
544 218 : if (!st->T) PetscCall(PetscCalloc1(PetscMax(2,st->nmat),&st->T));
545 : }
546 995 : PetscTryTypeMethod(st,setup);
547 995 : st->state = ST_STATE_SETUP;
548 995 : PetscCall(PetscLogEventEnd(ST_SetUp,st,0,0,0));
549 995 : 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 1621 : PetscErrorCode STMatMAXPY_Private(ST st,PetscScalar alpha,PetscScalar beta,PetscInt k,PetscScalar *coeffs,PetscBool initial,PetscBool precond,Mat *S)
564 : {
565 1621 : PetscInt *matIdx=NULL,nmat,i,ini=-1;
566 1621 : PetscScalar t=1.0,ta,gamma;
567 1621 : PetscBool nz=PETSC_FALSE;
568 1621 : Mat *A=precond?st->Psplit:st->A;
569 1621 : MatStructure str=precond?st->strp:st->str;
570 :
571 1621 : PetscFunctionBegin;
572 1621 : nmat = st->nmat-k;
573 1621 : switch (st->matmode) {
574 28 : case ST_MATMODE_INPLACE:
575 28 : PetscCheck(st->nmat<=2,PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"ST_MATMODE_INPLACE not supported for polynomial eigenproblems");
576 28 : PetscCheck(!precond,PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"ST_MATMODE_INPLACE not supported for split preconditioner");
577 28 : if (initial) {
578 16 : PetscCall(PetscObjectReference((PetscObject)A[0]));
579 16 : *S = A[0];
580 16 : gamma = alpha;
581 12 : } else gamma = alpha-beta;
582 28 : if (gamma != 0.0) {
583 23 : if (st->nmat>1) PetscCall(MatAXPY(*S,gamma,A[1],str));
584 18 : else PetscCall(MatShift(*S,gamma));
585 : }
586 : break;
587 107 : case ST_MATMODE_SHELL:
588 107 : PetscCheck(!precond,PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"ST_MATMODE_SHELL not supported for split preconditioner");
589 107 : if (initial) {
590 83 : if (st->nmat>2) {
591 10 : PetscCall(PetscMalloc1(nmat,&matIdx));
592 39 : for (i=0;i<nmat;i++) matIdx[i] = k+i;
593 : }
594 83 : PetscCall(STMatShellCreate(st,alpha,nmat,matIdx,coeffs,S));
595 83 : if (st->nmat>2) PetscCall(PetscFree(matIdx));
596 24 : } else PetscCall(STMatShellShift(*S,alpha));
597 : break;
598 1486 : case ST_MATMODE_COPY:
599 1486 : if (coeffs) {
600 552 : for (i=0;i<nmat && ini==-1;i++) {
601 320 : if (coeffs[i]!=0.0) ini = i;
602 88 : else t *= alpha;
603 : }
604 232 : if (coeffs[ini] != 1.0) nz = PETSC_TRUE;
605 439 : for (i=ini+1;i<nmat&&!nz;i++) if (coeffs[i]!=0.0) nz = PETSC_TRUE;
606 : } else { nz = PETSC_TRUE; ini = 0; }
607 1486 : if ((alpha == 0.0 || !nz) && t==1.0) {
608 619 : PetscCall(PetscObjectReference((PetscObject)A[k+ini]));
609 619 : PetscCall(MatDestroy(S));
610 619 : *S = A[k+ini];
611 : } else {
612 867 : if (*S && *S!=A[k+ini]) {
613 636 : PetscCall(MatSetOption(*S,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE));
614 636 : PetscCall(MatCopy(A[k+ini],*S,DIFFERENT_NONZERO_PATTERN));
615 : } else {
616 231 : PetscCall(MatDestroy(S));
617 231 : PetscCall(MatDuplicate(A[k+ini],MAT_COPY_VALUES,S));
618 231 : PetscCall(MatSetOption(*S,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE));
619 : }
620 867 : if (coeffs && coeffs[ini]!=1.0) PetscCall(MatScale(*S,coeffs[ini]));
621 1985 : for (i=ini+k+1;i<PetscMax(2,st->nmat);i++) {
622 1118 : t *= alpha;
623 1118 : ta = t;
624 1118 : if (coeffs) ta *= coeffs[i-k];
625 1118 : if (ta!=0.0) {
626 1117 : if (st->nmat>1) PetscCall(MatAXPY(*S,ta,A[i],str));
627 1118 : else PetscCall(MatShift(*S,ta));
628 : }
629 : }
630 : }
631 : }
632 1621 : PetscCall(MatSetOption(*S,MAT_SYMMETRIC,st->asymm));
633 1621 : PetscCall(MatSetOption(*S,MAT_HERMITIAN,(PetscImaginaryPart(st->sigma)==0.0)?st->aherm:PETSC_FALSE));
634 1621 : 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 41 : PetscErrorCode STCoeffs_Monomial(ST st, PetscScalar *coeffs)
642 : {
643 41 : PetscInt k,i,ini,inip;
644 :
645 41 : PetscFunctionBegin;
646 : /* Compute binomial coefficients */
647 41 : ini = (st->nmat*(st->nmat-1))/2;
648 174 : for (i=0;i<st->nmat;i++) coeffs[ini+i]=1.0;
649 133 : for (k=st->nmat-1;k>=1;k--) {
650 92 : inip = ini+1;
651 92 : ini = (k*(k-1))/2;
652 92 : coeffs[ini] = 1.0;
653 153 : for (i=1;i<k;i++) coeffs[ini+i] = coeffs[ini+i-1]+coeffs[inip+i-1];
654 : }
655 41 : 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 1016 : PetscErrorCode STPostSolve(ST st)
672 : {
673 1016 : PetscFunctionBegin;
674 1016 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
675 1016 : PetscValidType(st,1);
676 1016 : PetscTryTypeMethod(st,postsolve);
677 1016 : 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 1733786 : PetscErrorCode STBackTransform(ST st,PetscInt n,PetscScalar* eigr,PetscScalar* eigi)
698 : {
699 1733786 : PetscFunctionBegin;
700 1733786 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
701 1733786 : PetscValidType(st,1);
702 1733786 : PetscTryTypeMethod(st,backtransform,n,eigr,eigi);
703 1733786 : 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 607 : PetscErrorCode STIsInjective(ST st,PetscBool* is)
724 : {
725 607 : PetscBool shell;
726 :
727 607 : PetscFunctionBegin;
728 607 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
729 607 : PetscValidType(st,1);
730 607 : PetscAssertPointer(is,2);
731 :
732 607 : PetscCall(PetscObjectTypeCompare((PetscObject)st,STSHELL,&shell));
733 607 : if (shell) PetscCall(STIsInjective_Shell(st,is));
734 575 : else *is = st->ops->backtransform? PETSC_TRUE: PETSC_FALSE;
735 607 : 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 147 : PetscErrorCode STMatSetUp(ST st,PetscScalar sigma,PetscScalar *coeffs)
761 : {
762 147 : PetscFunctionBegin;
763 147 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
764 441 : PetscValidLogicalCollectiveScalar(st,sigma,2);
765 147 : STCheckMatrices(st,1);
766 :
767 147 : PetscCall(PetscLogEventBegin(ST_MatSetUp,st,0,0,0));
768 147 : PetscCall(STMatMAXPY_Private(st,sigma,0.0,0,coeffs,PETSC_TRUE,PETSC_FALSE,&st->P));
769 147 : if (st->Psplit) PetscCall(STMatMAXPY_Private(st,sigma,0.0,0,coeffs,PETSC_TRUE,PETSC_TRUE,&st->Pmat));
770 147 : PetscCall(ST_KSPSetOperators(st,st->P,st->Pmat?st->Pmat:st->P));
771 147 : PetscCall(KSPSetUp(st->ksp));
772 147 : PetscCall(PetscLogEventEnd(ST_MatSetUp,st,0,0,0));
773 147 : 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 359 : PetscErrorCode STSetWorkVecs(ST st,PetscInt nw)
793 : {
794 359 : PetscInt i;
795 :
796 359 : PetscFunctionBegin;
797 359 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
798 1077 : PetscValidLogicalCollectiveInt(st,nw,2);
799 359 : PetscCheck(nw>0,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_OUTOFRANGE,"nw must be > 0: nw = %" PetscInt_FMT,nw);
800 359 : if (st->nwork < nw) {
801 339 : PetscCall(VecDestroyVecs(st->nwork,&st->work));
802 339 : st->nwork = nw;
803 339 : PetscCall(PetscMalloc1(nw,&st->work));
804 715 : for (i=0;i<nw;i++) PetscCall(STMatCreateVecs(st,&st->work[i],NULL));
805 : }
806 359 : PetscFunctionReturn(PETSC_SUCCESS);
807 : }
|
