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 : Shift spectral transformation, applies (A + sigma I) as operator, or
12 : inv(B)(A + sigma B) for generalized problems
13 : */
14 :
15 : #include <slepc/private/stimpl.h>
16 :
17 1592369 : static PetscErrorCode STBackTransform_Shift(ST st,PetscInt n,PetscScalar *eigr,PetscScalar *eigi)
18 : {
19 1592369 : PetscInt j;
20 :
21 1592369 : PetscFunctionBegin;
22 4657740 : for (j=0;j<n;j++) {
23 3065371 : eigr[j] += st->sigma;
24 : }
25 1592369 : PetscFunctionReturn(PETSC_SUCCESS);
26 : }
27 :
28 567 : static PetscErrorCode STPostSolve_Shift(ST st)
29 : {
30 567 : PetscFunctionBegin;
31 567 : if (st->matmode == ST_MATMODE_INPLACE) {
32 3 : if (st->nmat>1) PetscCall(MatAXPY(st->A[0],st->sigma,st->A[1],st->str));
33 2 : else PetscCall(MatShift(st->A[0],st->sigma));
34 3 : st->Astate[0] = ((PetscObject)st->A[0])->state;
35 3 : st->state = ST_STATE_INITIAL;
36 3 : st->opready = PETSC_FALSE;
37 : }
38 567 : PetscFunctionReturn(PETSC_SUCCESS);
39 : }
40 :
41 : /*
42 : Operator (shift):
43 : Op P M
44 : if nmat=1: A-sI NULL A-sI
45 : if nmat=2: B^-1 (A-sB) B A-sB
46 : */
47 424 : static PetscErrorCode STComputeOperator_Shift(ST st)
48 : {
49 424 : PetscFunctionBegin;
50 424 : st->usesksp = (st->nmat>1)? PETSC_TRUE: PETSC_FALSE;
51 424 : PetscCall(PetscObjectReference((PetscObject)st->A[1]));
52 424 : PetscCall(MatDestroy(&st->T[1]));
53 424 : st->T[1] = st->A[1];
54 424 : PetscCall(STMatMAXPY_Private(st,-st->sigma,0.0,0,NULL,PetscNot(st->state==ST_STATE_UPDATED),PETSC_FALSE,&st->T[0]));
55 424 : if (st->nmat>1) PetscCall(PetscObjectReference((PetscObject)st->T[1]));
56 424 : PetscCall(MatDestroy(&st->P));
57 424 : st->P = (st->nmat>1)? st->T[1]: NULL;
58 424 : st->M = st->T[0];
59 424 : if (st->nmat>1 && st->Psplit) { /* build custom preconditioner from the split matrices */
60 1 : PetscCall(MatDestroy(&st->Pmat));
61 1 : PetscCall(PetscObjectReference((PetscObject)st->Psplit[1]));
62 1 : st->Pmat = st->Psplit[1];
63 : }
64 424 : PetscFunctionReturn(PETSC_SUCCESS);
65 : }
66 :
67 492 : static PetscErrorCode STSetUp_Shift(ST st)
68 : {
69 492 : PetscInt k,nc,nmat=st->nmat;
70 492 : PetscScalar *coeffs=NULL;
71 :
72 492 : PetscFunctionBegin;
73 492 : if (nmat>1) PetscCall(STSetWorkVecs(st,1));
74 492 : if (nmat>2) { /* set-up matrices for polynomial eigenproblems */
75 68 : if (st->transform) {
76 14 : nc = (nmat*(nmat+1))/2;
77 14 : PetscCall(PetscMalloc1(nc,&coeffs));
78 : /* Compute coeffs */
79 14 : PetscCall(STCoeffs_Monomial(st,coeffs));
80 : /* T[n] = A_n */
81 14 : k = nmat-1;
82 14 : PetscCall(PetscObjectReference((PetscObject)st->A[k]));
83 14 : PetscCall(MatDestroy(&st->T[k]));
84 14 : st->T[k] = st->A[k];
85 45 : for (k=0;k<nmat-1;k++) PetscCall(STMatMAXPY_Private(st,nmat>2?st->sigma:-st->sigma,0.0,k,coeffs?coeffs+((nmat-k)*(nmat-k-1))/2:NULL,PetscNot(st->state==ST_STATE_UPDATED),PETSC_FALSE,&st->T[k]));
86 14 : PetscCall(PetscFree(coeffs));
87 14 : PetscCall(PetscObjectReference((PetscObject)st->T[nmat-1]));
88 14 : PetscCall(MatDestroy(&st->P));
89 14 : st->P = st->T[nmat-1];
90 14 : if (st->Psplit) { /* build custom preconditioner from the split matrices */
91 1 : PetscCall(STMatMAXPY_Private(st,st->sigma,0.0,nmat-1,coeffs?coeffs:NULL,PETSC_TRUE,PETSC_TRUE,&st->Pmat));
92 : }
93 14 : PetscCall(ST_KSPSetOperators(st,st->P,st->Pmat?st->Pmat:st->P));
94 : } else {
95 224 : for (k=0;k<nmat;k++) {
96 170 : PetscCall(PetscObjectReference((PetscObject)st->A[k]));
97 170 : PetscCall(MatDestroy(&st->T[k]));
98 170 : st->T[k] = st->A[k];
99 : }
100 : }
101 : }
102 492 : if (st->P) PetscCall(KSPSetUp(st->ksp));
103 492 : PetscFunctionReturn(PETSC_SUCCESS);
104 : }
105 :
106 15 : static PetscErrorCode STSetShift_Shift(ST st,PetscScalar newshift)
107 : {
108 15 : PetscInt k,nc,nmat=PetscMax(st->nmat,2);
109 15 : PetscScalar *coeffs=NULL;
110 :
111 15 : PetscFunctionBegin;
112 15 : if (st->transform) {
113 13 : if (st->matmode == ST_MATMODE_COPY && nmat>2) {
114 2 : nc = (nmat*(nmat+1))/2;
115 2 : PetscCall(PetscMalloc1(nc,&coeffs));
116 : /* Compute coeffs */
117 2 : PetscCall(STCoeffs_Monomial(st,coeffs));
118 : }
119 32 : for (k=0;k<nmat-1;k++) PetscCall(STMatMAXPY_Private(st,nmat>2?newshift:-newshift,nmat>2?st->sigma:-st->sigma,k,coeffs?coeffs+((nmat-k)*(nmat-k-1))/2:NULL,PETSC_FALSE,PETSC_FALSE,&st->T[k]));
120 13 : if (st->matmode == ST_MATMODE_COPY && nmat>2) PetscCall(PetscFree(coeffs));
121 13 : if (st->nmat<=2) st->M = st->T[0];
122 : }
123 15 : PetscFunctionReturn(PETSC_SUCCESS);
124 : }
125 :
126 559 : SLEPC_EXTERN PetscErrorCode STCreate_Shift(ST st)
127 : {
128 559 : PetscFunctionBegin;
129 559 : st->usesksp = PETSC_TRUE;
130 :
131 559 : st->ops->apply = STApply_Generic;
132 559 : st->ops->applytrans = STApplyTranspose_Generic;
133 559 : st->ops->applyhermtrans = STApplyHermitianTranspose_Generic;
134 559 : st->ops->backtransform = STBackTransform_Shift;
135 559 : st->ops->setshift = STSetShift_Shift;
136 559 : st->ops->getbilinearform = STGetBilinearForm_Default;
137 559 : st->ops->setup = STSetUp_Shift;
138 559 : st->ops->computeoperator = STComputeOperator_Shift;
139 559 : st->ops->postsolve = STPostSolve_Shift;
140 559 : st->ops->setdefaultksp = STSetDefaultKSP_Default;
141 559 : PetscFunctionReturn(PETSC_SUCCESS);
142 : }
|
