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 97661 : PetscErrorCode STApply_Generic(ST st,Vec x,Vec y)
17 : {
18 97661 : PetscFunctionBegin;
19 97661 : if (st->M && st->P) {
20 18320 : PetscCall(MatMult(st->M,x,st->work[0]));
21 18320 : PetscCall(STMatSolve(st,st->work[0],y));
22 79341 : } else if (st->M) PetscCall(MatMult(st->M,x,y));
23 16636 : else PetscCall(STMatSolve(st,x,y));
24 97661 : 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 45331 : PetscErrorCode STApply(ST st,Vec x,Vec y)
46 : {
47 45331 : Mat Op;
48 :
49 45331 : PetscFunctionBegin;
50 45331 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
51 45331 : PetscValidHeaderSpecific(x,VEC_CLASSID,2);
52 45331 : PetscValidHeaderSpecific(y,VEC_CLASSID,3);
53 45331 : PetscValidType(st,1);
54 45331 : STCheckMatrices(st,1);
55 45331 : PetscCheck(x!=y,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
56 45331 : PetscCall(VecSetErrorIfLocked(y,3));
57 45331 : PetscCall(STGetOperator_Private(st,&Op));
58 45331 : PetscCall(MatMult(Op,x,y));
59 45331 : PetscFunctionReturn(PETSC_SUCCESS);
60 : }
61 :
62 5141 : PetscErrorCode STApplyMat_Generic(ST st,Mat B,Mat C)
63 : {
64 5141 : Mat work;
65 :
66 5141 : PetscFunctionBegin;
67 5141 : if (st->M && st->P) {
68 1998 : PetscCall(MatMatMult(st->M,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&work));
69 1998 : PetscCall(STMatMatSolve(st,work,C));
70 1998 : PetscCall(MatDestroy(&work));
71 3143 : } else if (st->M) PetscCall(MatMatMult(st->M,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C));
72 1456 : else PetscCall(STMatMatSolve(st,B,C));
73 5141 : 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 1412 : PetscErrorCode STApplyMat(ST st,Mat X,Mat Y)
95 : {
96 1412 : PetscFunctionBegin;
97 1412 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
98 1412 : PetscValidHeaderSpecific(X,MAT_CLASSID,2);
99 1412 : PetscValidHeaderSpecific(Y,MAT_CLASSID,3);
100 1412 : PetscValidType(st,1);
101 1412 : STCheckMatrices(st,1);
102 1412 : PetscCheck(X!=Y,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"X and Y must be different matrices");
103 1412 : PetscUseTypeMethod(st,applymat,X,Y);
104 1412 : 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 1314 : PetscErrorCode STApplyHermitianTranspose_Generic(ST st,Vec x,Vec y)
154 : {
155 1314 : PetscFunctionBegin;
156 1314 : 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 1257 : } else if (st->M) PetscCall(MatMultHermitianTranspose(st->M,x,y));
160 527 : else PetscCall(STMatSolveHermitianTranspose(st,x,y));
161 1314 : 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 474 : PetscErrorCode STApplyHermitianTranspose(ST st,Vec x,Vec y)
186 : {
187 474 : Mat Op;
188 :
189 474 : PetscFunctionBegin;
190 474 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
191 474 : PetscValidHeaderSpecific(x,VEC_CLASSID,2);
192 474 : PetscValidHeaderSpecific(y,VEC_CLASSID,3);
193 474 : PetscValidType(st,1);
194 474 : STCheckMatrices(st,1);
195 474 : PetscCheck(x!=y,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
196 474 : PetscCall(VecSetErrorIfLocked(y,3));
197 474 : PetscCall(STGetOperator_Private(st,&Op));
198 474 : PetscCall(MatMultHermitianTranspose(Op,x,y));
199 474 : 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 100789 : static PetscErrorCode MatMult_STOperator(Mat Op,Vec x,Vec y)
245 : {
246 100789 : ST st;
247 :
248 100789 : PetscFunctionBegin;
249 100789 : PetscCall(MatShellGetContext(Op,&st));
250 100789 : PetscCall(STSetUp(st));
251 100789 : PetscCall(PetscLogEventBegin(ST_Apply,st,x,y,0));
252 100789 : 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 100077 : } else PetscUseTypeMethod(st,apply,x,y);
257 100789 : PetscCall(PetscLogEventEnd(ST_Apply,st,x,y,0));
258 100789 : 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 1390 : static PetscErrorCode MatMultHermitianTranspose_STOperator(Mat Op,Vec x,Vec y)
280 : {
281 1390 : ST st;
282 :
283 1390 : PetscFunctionBegin;
284 1390 : PetscCall(MatShellGetContext(Op,&st));
285 1390 : PetscCall(STSetUp(st));
286 1390 : PetscCall(PetscLogEventBegin(ST_ApplyHermitianTranspose,st,x,y,0));
287 1390 : if (st->ops->applyhermtrans) {
288 1390 : 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 1340 : } 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 1390 : PetscCall(PetscLogEventEnd(ST_ApplyHermitianTranspose,st,x,y,0));
307 1390 : PetscFunctionReturn(PETSC_SUCCESS);
308 : }
309 : #endif
310 :
311 3729 : static PetscErrorCode MatMatMult_STOperator(Mat Op,Mat B,Mat C,void *ctx)
312 : {
313 3729 : ST st;
314 :
315 3729 : PetscFunctionBegin;
316 3729 : PetscCall(MatShellGetContext(Op,&st));
317 3729 : PetscCall(STSetUp(st));
318 3729 : PetscCall(PetscLogEventBegin(ST_Apply,st,B,C,0));
319 3729 : PetscCall(STApplyMat_Generic(st,B,C));
320 3729 : PetscCall(PetscLogEventEnd(ST_Apply,st,B,C,0));
321 3729 : PetscFunctionReturn(PETSC_SUCCESS);
322 : }
323 :
324 49824 : PetscErrorCode STGetOperator_Private(ST st,Mat *Op)
325 : {
326 49824 : PetscInt m,n,M,N;
327 49824 : Vec v;
328 49824 : VecType vtype;
329 :
330 49824 : PetscFunctionBegin;
331 49824 : if (!st->Op) {
332 469 : if (Op) *Op = NULL;
333 : /* create the shell matrix */
334 469 : PetscCall(MatGetLocalSize(st->A[0],&m,&n));
335 469 : PetscCall(MatGetSize(st->A[0],&M,&N));
336 469 : PetscCall(MatCreateShell(PetscObjectComm((PetscObject)st),m,n,M,N,st,&st->Op));
337 469 : PetscCall(MatShellSetOperation(st->Op,MATOP_MULT,(void(*)(void))MatMult_STOperator));
338 469 : PetscCall(MatShellSetOperation(st->Op,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_STOperator));
339 : #if defined(PETSC_USE_COMPLEX)
340 469 : 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 469 : if (!st->D && st->ops->apply==STApply_Generic) {
345 437 : PetscCall(MatShellSetMatProductOperation(st->Op,MATPRODUCT_AB,NULL,MatMatMult_STOperator,NULL,MATDENSE,MATDENSE));
346 437 : PetscCall(MatShellSetMatProductOperation(st->Op,MATPRODUCT_AB,NULL,MatMatMult_STOperator,NULL,MATDENSECUDA,MATDENSECUDA));
347 : }
348 : /* make sure the shell matrix generates a vector of the same type as the problem matrices */
349 469 : PetscCall(MatCreateVecs(st->A[0],&v,NULL));
350 469 : PetscCall(VecGetType(v,&vtype));
351 469 : PetscCall(MatShellSetVecType(st->Op,vtype));
352 469 : PetscCall(VecDestroy(&v));
353 : /* build the operator matrices */
354 469 : PetscCall(STComputeOperator(st));
355 : }
356 49824 : if (Op) *Op = st->Op;
357 49824 : PetscFunctionReturn(PETSC_SUCCESS);
358 : }
359 :
360 : /*@
361 : STGetOperator - Returns a shell matrix that represents the operator of the
362 : spectral transformation.
363 :
364 : Collective
365 :
366 : Input Parameter:
367 : . st - the spectral transformation context
368 :
369 : Output Parameter:
370 : . Op - operator matrix
371 :
372 : Notes:
373 : The operator is defined in linear eigenproblems only, not in polynomial ones,
374 : so the call will fail if more than 2 matrices were passed in STSetMatrices().
375 :
376 : The returned shell matrix is essentially a wrapper to the STApply() and
377 : STApplyTranspose() operations. The operator can often be expressed as
378 :
379 : $ Op = D*inv(K)*M*inv(D)
380 :
381 : where D is the balancing matrix, and M and K are two matrices corresponding
382 : to the numerator and denominator for spectral transformations that represent
383 : a rational matrix function. In the case of STSHELL, the inner part inv(K)*M
384 : is replaced by the user-provided operation from STShellSetApply().
385 :
386 : The preconditioner matrix K typically depends on the value of the shift, and
387 : its inverse is handled via an internal KSP object. Normal usage does not
388 : require explicitly calling STGetOperator(), but it can be used to force the
389 : creation of K and M, and then K is passed to the KSP. This is useful for
390 : setting options associated with the PCFactor (to set MUMPS options, for instance).
391 :
392 : The returned matrix must NOT be destroyed by the user. Instead, when no
393 : longer needed it must be returned with STRestoreOperator(). In particular,
394 : this is required before modifying the ST matrices or the shift.
395 :
396 : A NULL pointer can be passed in Op in case the matrix is not required but we
397 : want to force its creation. In this case, STRestoreOperator() should not be
398 : called.
399 :
400 : Level: advanced
401 :
402 : .seealso: STApply(), STApplyTranspose(), STSetBalanceMatrix(), STShellSetApply(),
403 : STGetKSP(), STSetShift(), STRestoreOperator(), STSetMatrices()
404 : @*/
405 4007 : PetscErrorCode STGetOperator(ST st,Mat *Op)
406 : {
407 4007 : PetscFunctionBegin;
408 4007 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
409 4007 : PetscValidType(st,1);
410 4007 : STCheckMatrices(st,1);
411 4007 : STCheckNotSeized(st,1);
412 4007 : PetscCheck(st->nmat<=2,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_WRONGSTATE,"The operator is not defined in polynomial eigenproblems");
413 4007 : PetscCall(STGetOperator_Private(st,Op));
414 4007 : if (Op) st->opseized = PETSC_TRUE;
415 4007 : PetscFunctionReturn(PETSC_SUCCESS);
416 : }
417 :
418 : /*@
419 : STRestoreOperator - Restore the previously seized operator matrix.
420 :
421 : Logically Collective
422 :
423 : Input Parameters:
424 : + st - the spectral transformation context
425 : - Op - operator matrix
426 :
427 : Notes:
428 : The arguments must match the corresponding call to STGetOperator().
429 :
430 : Level: advanced
431 :
432 : .seealso: STGetOperator()
433 : @*/
434 3977 : PetscErrorCode STRestoreOperator(ST st,Mat *Op)
435 : {
436 3977 : PetscFunctionBegin;
437 3977 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
438 3977 : PetscAssertPointer(Op,2);
439 3977 : PetscValidHeaderSpecific(*Op,MAT_CLASSID,2);
440 3977 : PetscCheck(st->opseized,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_WRONGSTATE,"Must be called after STGetOperator()");
441 3977 : *Op = NULL;
442 3977 : st->opseized = PETSC_FALSE;
443 3977 : PetscFunctionReturn(PETSC_SUCCESS);
444 : }
445 :
446 : /*
447 : STComputeOperator - Computes the matrices that constitute the operator
448 :
449 : Op = D*inv(K)*M*inv(D).
450 :
451 : K and M are computed here (D is user-provided) from the system matrices
452 : and the shift sigma (whenever these are changed, this function recomputes
453 : K and M). This is used only in linear eigenproblems (nmat<3).
454 :
455 : K is the "preconditioner matrix": it is the denominator in rational operators,
456 : e.g. (A-sigma*B) in shift-and-invert. In non-rational transformations such
457 : as STFILTER, K=NULL which means identity. After computing K, it is passed to
458 : the internal KSP object via KSPSetOperators.
459 :
460 : M is the numerator in rational operators. If unused it is set to NULL (e.g.
461 : in STPRECOND).
462 :
463 : STSHELL does not compute anything here, but sets the flag as if it was ready.
464 : */
465 1237 : PetscErrorCode STComputeOperator(ST st)
466 : {
467 1237 : PC pc;
468 :
469 1237 : PetscFunctionBegin;
470 1237 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
471 1237 : PetscValidType(st,1);
472 1237 : if (!st->opready && st->ops->computeoperator) {
473 748 : PetscCall(PetscInfo(st,"Building the operator matrices\n"));
474 748 : STCheckMatrices(st,1);
475 748 : if (!st->T) PetscCall(PetscCalloc1(PetscMax(2,st->nmat),&st->T));
476 748 : PetscCall(PetscLogEventBegin(ST_ComputeOperator,st,0,0,0));
477 748 : PetscUseTypeMethod(st,computeoperator);
478 748 : PetscCall(PetscLogEventEnd(ST_ComputeOperator,st,0,0,0));
479 748 : if (st->usesksp) {
480 411 : if (!st->ksp) PetscCall(STGetKSP(st,&st->ksp));
481 411 : if (st->P) {
482 372 : PetscCall(STSetDefaultKSP(st));
483 372 : PetscCall(ST_KSPSetOperators(st,st->P,st->Pmat?st->Pmat:st->P));
484 : } else {
485 : /* STPRECOND defaults to PCNONE if st->P is empty */
486 39 : PetscCall(KSPGetPC(st->ksp,&pc));
487 39 : PetscCall(PCSetType(pc,PCNONE));
488 : }
489 : }
490 : }
491 1237 : st->opready = PETSC_TRUE;
492 1237 : PetscFunctionReturn(PETSC_SUCCESS);
493 : }
494 :
495 : /*@
496 : STSetUp - Prepares for the use of a spectral transformation.
497 :
498 : Collective
499 :
500 : Input Parameter:
501 : . st - the spectral transformation context
502 :
503 : Level: advanced
504 :
505 : .seealso: STCreate(), STApply(), STDestroy()
506 : @*/
507 107175 : PetscErrorCode STSetUp(ST st)
508 : {
509 107175 : PetscInt i,n,k;
510 :
511 107175 : PetscFunctionBegin;
512 107175 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
513 107175 : PetscValidType(st,1);
514 107175 : STCheckMatrices(st,1);
515 107175 : switch (st->state) {
516 952 : case ST_STATE_INITIAL:
517 952 : PetscCall(PetscInfo(st,"Setting up new ST\n"));
518 952 : if (!((PetscObject)st)->type_name) PetscCall(STSetType(st,STSHIFT));
519 : break;
520 106185 : case ST_STATE_SETUP:
521 106185 : PetscFunctionReturn(PETSC_SUCCESS);
522 38 : case ST_STATE_UPDATED:
523 38 : PetscCall(PetscInfo(st,"Setting up updated ST\n"));
524 : break;
525 : }
526 990 : PetscCall(PetscLogEventBegin(ST_SetUp,st,0,0,0));
527 990 : if (st->state!=ST_STATE_UPDATED) {
528 952 : if (!(st->nmat<3 && st->opready)) {
529 924 : if (st->T) {
530 276 : for (i=0;i<PetscMax(2,st->nmat);i++) PetscCall(MatDestroy(&st->T[i]));
531 : }
532 924 : PetscCall(MatDestroy(&st->P));
533 : }
534 : }
535 990 : if (st->D) {
536 13 : PetscCall(MatGetLocalSize(st->A[0],NULL,&n));
537 13 : PetscCall(VecGetLocalSize(st->D,&k));
538 13 : PetscCheck(n==k,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_SIZ,"Balance matrix has wrong dimension %" PetscInt_FMT " (should be %" PetscInt_FMT ")",k,n);
539 13 : if (!st->wb) PetscCall(VecDuplicate(st->D,&st->wb));
540 : }
541 990 : if (st->nmat<3 && st->transform) PetscCall(STComputeOperator(st));
542 : else {
543 222 : if (!st->T) PetscCall(PetscCalloc1(PetscMax(2,st->nmat),&st->T));
544 : }
545 990 : PetscTryTypeMethod(st,setup);
546 990 : st->state = ST_STATE_SETUP;
547 990 : PetscCall(PetscLogEventEnd(ST_SetUp,st,0,0,0));
548 990 : PetscFunctionReturn(PETSC_SUCCESS);
549 : }
550 :
551 : /*
552 : Computes coefficients for the transformed polynomial,
553 : and stores the result in argument S.
554 :
555 : alpha - value of the parameter of the transformed polynomial
556 : beta - value of the previous shift (only used in inplace mode)
557 : k - index of first matrix included in the computation
558 : coeffs - coefficients of the expansion
559 : initial - true if this is the first time
560 : precond - whether the preconditioner matrix must be computed
561 : */
562 1617 : PetscErrorCode STMatMAXPY_Private(ST st,PetscScalar alpha,PetscScalar beta,PetscInt k,PetscScalar *coeffs,PetscBool initial,PetscBool precond,Mat *S)
563 : {
564 1617 : PetscInt *matIdx=NULL,nmat,i,ini=-1;
565 1617 : PetscScalar t=1.0,ta,gamma;
566 1617 : PetscBool nz=PETSC_FALSE;
567 1617 : Mat *A=precond?st->Psplit:st->A;
568 1617 : MatStructure str=precond?st->strp:st->str;
569 :
570 1617 : PetscFunctionBegin;
571 1617 : nmat = st->nmat-k;
572 1617 : switch (st->matmode) {
573 28 : case ST_MATMODE_INPLACE:
574 28 : PetscCheck(st->nmat<=2,PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"ST_MATMODE_INPLACE not supported for polynomial eigenproblems");
575 28 : PetscCheck(!precond,PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"ST_MATMODE_INPLACE not supported for split preconditioner");
576 28 : if (initial) {
577 16 : PetscCall(PetscObjectReference((PetscObject)A[0]));
578 16 : *S = A[0];
579 16 : gamma = alpha;
580 12 : } else gamma = alpha-beta;
581 28 : if (gamma != 0.0) {
582 23 : if (st->nmat>1) PetscCall(MatAXPY(*S,gamma,A[1],str));
583 18 : else PetscCall(MatShift(*S,gamma));
584 : }
585 : break;
586 109 : case ST_MATMODE_SHELL:
587 109 : PetscCheck(!precond,PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"ST_MATMODE_SHELL not supported for split preconditioner");
588 109 : if (initial) {
589 85 : if (st->nmat>2) {
590 10 : PetscCall(PetscMalloc1(nmat,&matIdx));
591 39 : for (i=0;i<nmat;i++) matIdx[i] = k+i;
592 : }
593 85 : PetscCall(STMatShellCreate(st,alpha,nmat,matIdx,coeffs,S));
594 85 : if (st->nmat>2) PetscCall(PetscFree(matIdx));
595 24 : } else PetscCall(STMatShellShift(*S,alpha));
596 : break;
597 1480 : case ST_MATMODE_COPY:
598 1480 : if (coeffs) {
599 564 : for (i=0;i<nmat && ini==-1;i++) {
600 328 : if (coeffs[i]!=0.0) ini = i;
601 92 : else t *= alpha;
602 : }
603 236 : if (coeffs[ini] != 1.0) nz = PETSC_TRUE;
604 445 : for (i=ini+1;i<nmat&&!nz;i++) if (coeffs[i]!=0.0) nz = PETSC_TRUE;
605 : } else { nz = PETSC_TRUE; ini = 0; }
606 1480 : if ((alpha == 0.0 || !nz) && t==1.0) {
607 610 : PetscCall(PetscObjectReference((PetscObject)A[k+ini]));
608 610 : PetscCall(MatDestroy(S));
609 610 : *S = A[k+ini];
610 : } else {
611 870 : if (*S && *S!=A[k+ini]) {
612 637 : PetscCall(MatSetOption(*S,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE));
613 637 : PetscCall(MatCopy(A[k+ini],*S,DIFFERENT_NONZERO_PATTERN));
614 : } else {
615 233 : PetscCall(MatDestroy(S));
616 233 : PetscCall(MatDuplicate(A[k+ini],MAT_COPY_VALUES,S));
617 233 : PetscCall(MatSetOption(*S,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE));
618 : }
619 870 : if (coeffs && coeffs[ini]!=1.0) PetscCall(MatScale(*S,coeffs[ini]));
620 1994 : for (i=ini+k+1;i<PetscMax(2,st->nmat);i++) {
621 1124 : t *= alpha;
622 1124 : ta = t;
623 1124 : if (coeffs) ta *= coeffs[i-k];
624 1124 : if (ta!=0.0) {
625 1123 : if (st->nmat>1) PetscCall(MatAXPY(*S,ta,A[i],str));
626 1124 : else PetscCall(MatShift(*S,ta));
627 : }
628 : }
629 : }
630 : }
631 1617 : PetscCall(MatSetOption(*S,MAT_SYMMETRIC,st->asymm));
632 1617 : PetscCall(MatSetOption(*S,MAT_HERMITIAN,(PetscImaginaryPart(st->sigma)==0.0)?st->aherm:PETSC_FALSE));
633 1617 : PetscFunctionReturn(PETSC_SUCCESS);
634 : }
635 :
636 : /*
637 : Computes the values of the coefficients required by STMatMAXPY_Private
638 : for the case of monomial basis.
639 : */
640 41 : PetscErrorCode STCoeffs_Monomial(ST st, PetscScalar *coeffs)
641 : {
642 41 : PetscInt k,i,ini,inip;
643 :
644 41 : PetscFunctionBegin;
645 : /* Compute binomial coefficients */
646 41 : ini = (st->nmat*(st->nmat-1))/2;
647 174 : for (i=0;i<st->nmat;i++) coeffs[ini+i]=1.0;
648 133 : for (k=st->nmat-1;k>=1;k--) {
649 92 : inip = ini+1;
650 92 : ini = (k*(k-1))/2;
651 92 : coeffs[ini] = 1.0;
652 153 : for (i=1;i<k;i++) coeffs[ini+i] = coeffs[ini+i-1]+coeffs[inip+i-1];
653 : }
654 41 : PetscFunctionReturn(PETSC_SUCCESS);
655 : }
656 :
657 : /*@
658 : STPostSolve - Optional post-solve phase, intended for any actions that must
659 : be performed on the ST object after the eigensolver has finished.
660 :
661 : Collective
662 :
663 : Input Parameters:
664 : . st - the spectral transformation context
665 :
666 : Level: developer
667 :
668 : .seealso: EPSSolve()
669 : @*/
670 1021 : PetscErrorCode STPostSolve(ST st)
671 : {
672 1021 : PetscFunctionBegin;
673 1021 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
674 1021 : PetscValidType(st,1);
675 1021 : PetscTryTypeMethod(st,postsolve);
676 1021 : PetscFunctionReturn(PETSC_SUCCESS);
677 : }
678 :
679 : /*@
680 : STBackTransform - Back-transformation phase, intended for
681 : spectral transformations which require to transform the computed
682 : eigenvalues back to the original eigenvalue problem.
683 :
684 : Not Collective
685 :
686 : Input Parameters:
687 : + st - the spectral transformation context
688 : . n - number of eigenvalues
689 : . eigr - real part of a computed eigenvalues
690 : - eigi - imaginary part of a computed eigenvalues
691 :
692 : Level: developer
693 :
694 : .seealso: STIsInjective()
695 : @*/
696 1751132 : PetscErrorCode STBackTransform(ST st,PetscInt n,PetscScalar* eigr,PetscScalar* eigi)
697 : {
698 1751132 : PetscFunctionBegin;
699 1751132 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
700 1751132 : PetscValidType(st,1);
701 1751132 : PetscTryTypeMethod(st,backtransform,n,eigr,eigi);
702 1751132 : PetscFunctionReturn(PETSC_SUCCESS);
703 : }
704 :
705 : /*@
706 : STIsInjective - Ask if this spectral transformation is injective or not
707 : (that is, if it corresponds to a one-to-one mapping). If not, then it
708 : does not make sense to call STBackTransform().
709 :
710 : Not Collective
711 :
712 : Input Parameter:
713 : . st - the spectral transformation context
714 :
715 : Output Parameter:
716 : . is - the answer
717 :
718 : Level: developer
719 :
720 : .seealso: STBackTransform()
721 : @*/
722 608 : PetscErrorCode STIsInjective(ST st,PetscBool* is)
723 : {
724 608 : PetscBool shell;
725 :
726 608 : PetscFunctionBegin;
727 608 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
728 608 : PetscValidType(st,1);
729 608 : PetscAssertPointer(is,2);
730 :
731 608 : PetscCall(PetscObjectTypeCompare((PetscObject)st,STSHELL,&shell));
732 608 : if (shell) PetscCall(STIsInjective_Shell(st,is));
733 576 : else *is = st->ops->backtransform? PETSC_TRUE: PETSC_FALSE;
734 608 : PetscFunctionReturn(PETSC_SUCCESS);
735 : }
736 :
737 : /*@
738 : STMatSetUp - Build the preconditioner matrix used in STMatSolve().
739 :
740 : Collective
741 :
742 : Input Parameters:
743 : + st - the spectral transformation context
744 : . sigma - the shift
745 : - coeffs - the coefficients (may be NULL)
746 :
747 : Note:
748 : This function is not intended to be called by end users, but by SLEPc
749 : solvers that use ST. It builds matrix st->P as follows, then calls KSPSetUp().
750 : .vb
751 : If (coeffs) st->P = Sum_{i=0..nmat-1} coeffs[i]*sigma^i*A_i
752 : else st->P = Sum_{i=0..nmat-1} sigma^i*A_i
753 : .ve
754 :
755 : Level: developer
756 :
757 : .seealso: STMatSolve()
758 : @*/
759 151 : PetscErrorCode STMatSetUp(ST st,PetscScalar sigma,PetscScalar *coeffs)
760 : {
761 151 : PetscFunctionBegin;
762 151 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
763 604 : PetscValidLogicalCollectiveScalar(st,sigma,2);
764 151 : STCheckMatrices(st,1);
765 :
766 151 : PetscCall(PetscLogEventBegin(ST_MatSetUp,st,0,0,0));
767 151 : PetscCall(STMatMAXPY_Private(st,sigma,0.0,0,coeffs,PETSC_TRUE,PETSC_FALSE,&st->P));
768 151 : if (st->Psplit) PetscCall(STMatMAXPY_Private(st,sigma,0.0,0,coeffs,PETSC_TRUE,PETSC_TRUE,&st->Pmat));
769 151 : PetscCall(ST_KSPSetOperators(st,st->P,st->Pmat?st->Pmat:st->P));
770 151 : PetscCall(KSPSetUp(st->ksp));
771 151 : PetscCall(PetscLogEventEnd(ST_MatSetUp,st,0,0,0));
772 151 : PetscFunctionReturn(PETSC_SUCCESS);
773 : }
774 :
775 : /*@
776 : STSetWorkVecs - Sets a number of work vectors into the ST object.
777 :
778 : Collective
779 :
780 : Input Parameters:
781 : + st - the spectral transformation context
782 : - nw - number of work vectors to allocate
783 :
784 : Developer Notes:
785 : This is SLEPC_EXTERN because it may be required by shell STs.
786 :
787 : Level: developer
788 :
789 : .seealso: STMatCreateVecs()
790 : @*/
791 363 : PetscErrorCode STSetWorkVecs(ST st,PetscInt nw)
792 : {
793 363 : PetscInt i;
794 :
795 363 : PetscFunctionBegin;
796 363 : PetscValidHeaderSpecific(st,ST_CLASSID,1);
797 1452 : PetscValidLogicalCollectiveInt(st,nw,2);
798 363 : PetscCheck(nw>0,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_OUTOFRANGE,"nw must be > 0: nw = %" PetscInt_FMT,nw);
799 363 : if (st->nwork < nw) {
800 343 : PetscCall(VecDestroyVecs(st->nwork,&st->work));
801 343 : st->nwork = nw;
802 343 : PetscCall(PetscMalloc1(nw,&st->work));
803 723 : for (i=0;i<nw;i++) PetscCall(STMatCreateVecs(st,&st->work[i],NULL));
804 : }
805 363 : PetscFunctionReturn(PETSC_SUCCESS);
806 : }
|
