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 : /*
18 : Special STApply() for the BSE structured matrix
19 :
20 : H = [ R C; -C^H -R^T ].
21 :
22 : Assumes that H is a MATNEST and x,y are VECNEST.
23 :
24 : Only the upper part of the product y=H*x is computed.
25 :
26 : y1 = R*x1+C*x2
27 :
28 : The bottom part of the input vector x2 should normally contain
29 : either conj(x1) or -conj(x1).
30 : The bottom part of the output vector y2 is not referenced.
31 : */
32 8423 : static PetscErrorCode STApply_Shift_BSE(ST st,Vec x,Vec y)
33 : {
34 8423 : Mat H,R,C;
35 8423 : Vec x1,x2,y1;
36 :
37 8423 : PetscFunctionBegin;
38 8423 : H = st->T[0];
39 8423 : PetscCall(MatNestGetSubMat(H,0,0,&R));
40 8423 : PetscCall(MatNestGetSubMat(H,0,1,&C));
41 8423 : PetscCall(VecNestGetSubVec(x,0,&x1));
42 8423 : PetscCall(VecNestGetSubVec(x,1,&x2));
43 8423 : PetscCall(VecNestGetSubVec(y,0,&y1));
44 8423 : PetscCall(MatMult(C,x2,y1));
45 8423 : PetscCall(MatMultAdd(R,x1,y1,y1));
46 8423 : PetscFunctionReturn(PETSC_SUCCESS);
47 : }
48 :
49 853324 : static PetscErrorCode STBackTransform_Shift(ST st,PetscInt n,PetscScalar *eigr,PetscScalar *eigi)
50 : {
51 853324 : PetscInt j;
52 :
53 853324 : PetscFunctionBegin;
54 2548173 : for (j=0;j<n;j++) {
55 1694849 : eigr[j] += st->sigma;
56 : }
57 853324 : PetscFunctionReturn(PETSC_SUCCESS);
58 : }
59 :
60 573 : static PetscErrorCode STPostSolve_Shift(ST st)
61 : {
62 573 : PetscFunctionBegin;
63 573 : if (st->matmode == ST_MATMODE_INPLACE) {
64 4 : if (st->nmat>1) PetscCall(MatAXPY(st->A[0],st->sigma,st->A[1],st->str));
65 3 : else PetscCall(MatShift(st->A[0],st->sigma));
66 4 : st->Astate[0] = ((PetscObject)st->A[0])->state;
67 4 : st->state = ST_STATE_INITIAL;
68 4 : st->opready = PETSC_FALSE;
69 : }
70 573 : PetscFunctionReturn(PETSC_SUCCESS);
71 : }
72 :
73 : /*
74 : Operator (shift):
75 : Op P M
76 : if nmat=1: A-sI NULL A-sI
77 : if nmat=2: B^-1 (A-sB) B A-sB
78 : */
79 424 : static PetscErrorCode STComputeOperator_Shift(ST st)
80 : {
81 424 : PetscFunctionBegin;
82 424 : st->usesksp = (st->nmat>1)? PETSC_TRUE: PETSC_FALSE;
83 424 : PetscCall(PetscObjectReference((PetscObject)st->A[1]));
84 424 : PetscCall(MatDestroy(&st->T[1]));
85 424 : st->T[1] = st->A[1];
86 424 : PetscCall(STMatMAXPY_Private(st,-st->sigma,0.0,0,NULL,PetscNot(st->state==ST_STATE_UPDATED),PETSC_FALSE,&st->T[0]));
87 424 : if (st->nmat>1) PetscCall(PetscObjectReference((PetscObject)st->T[1]));
88 424 : PetscCall(MatDestroy(&st->P));
89 424 : st->P = (st->nmat>1)? st->T[1]: NULL;
90 424 : st->M = st->T[0];
91 424 : if (st->nmat>1 && st->Psplit) { /* build custom preconditioner from the split matrices */
92 1 : PetscCall(MatDestroy(&st->Pmat));
93 1 : PetscCall(PetscObjectReference((PetscObject)st->Psplit[1]));
94 1 : st->Pmat = st->Psplit[1];
95 : }
96 424 : PetscFunctionReturn(PETSC_SUCCESS);
97 : }
98 :
99 501 : static PetscErrorCode STSetUp_Shift(ST st)
100 : {
101 501 : PetscInt k,nc,nmat=st->nmat;
102 501 : PetscScalar *coeffs=NULL;
103 :
104 501 : PetscFunctionBegin;
105 501 : if (nmat>1) PetscCall(STSetWorkVecs(st,1));
106 501 : if (nmat>2) { /* set-up matrices for polynomial eigenproblems */
107 77 : if (st->transform) {
108 14 : nc = (nmat*(nmat+1))/2;
109 14 : PetscCall(PetscMalloc1(nc,&coeffs));
110 : /* Compute coeffs */
111 14 : PetscCall(STCoeffs_Monomial(st,coeffs));
112 : /* T[n] = A_n */
113 14 : k = nmat-1;
114 14 : PetscCall(PetscObjectReference((PetscObject)st->A[k]));
115 14 : PetscCall(MatDestroy(&st->T[k]));
116 14 : st->T[k] = st->A[k];
117 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]));
118 14 : PetscCall(PetscFree(coeffs));
119 14 : PetscCall(PetscObjectReference((PetscObject)st->T[nmat-1]));
120 14 : PetscCall(MatDestroy(&st->P));
121 14 : st->P = st->T[nmat-1];
122 14 : if (st->Psplit) { /* build custom preconditioner from the split matrices */
123 1 : PetscCall(STMatMAXPY_Private(st,st->sigma,0.0,nmat-1,coeffs?coeffs:NULL,PETSC_TRUE,PETSC_TRUE,&st->Pmat));
124 : }
125 14 : PetscCall(ST_KSPSetOperators(st,st->P,st->Pmat?st->Pmat:st->P));
126 : } else {
127 276 : for (k=0;k<nmat;k++) {
128 213 : PetscCall(PetscObjectReference((PetscObject)st->A[k]));
129 213 : PetscCall(MatDestroy(&st->T[k]));
130 213 : st->T[k] = st->A[k];
131 : }
132 : }
133 : }
134 501 : if (st->P) PetscCall(KSPSetUp(st->ksp));
135 501 : if (st->structured) st->ops->apply = STApply_Shift_BSE;
136 501 : PetscFunctionReturn(PETSC_SUCCESS);
137 : }
138 :
139 18 : static PetscErrorCode STSetShift_Shift(ST st,PetscScalar newshift)
140 : {
141 18 : PetscInt k,nc,nmat=PetscMax(st->nmat,2);
142 18 : PetscScalar *coeffs=NULL;
143 :
144 18 : PetscFunctionBegin;
145 18 : if (st->transform) {
146 16 : if (st->matmode == ST_MATMODE_COPY && nmat>2) {
147 2 : nc = (nmat*(nmat+1))/2;
148 2 : PetscCall(PetscMalloc1(nc,&coeffs));
149 : /* Compute coeffs */
150 2 : PetscCall(STCoeffs_Monomial(st,coeffs));
151 : }
152 38 : 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]));
153 16 : if (st->matmode == ST_MATMODE_COPY && nmat>2) PetscCall(PetscFree(coeffs));
154 16 : if (st->nmat<=2) st->M = st->T[0];
155 : }
156 18 : PetscFunctionReturn(PETSC_SUCCESS);
157 : }
158 :
159 576 : SLEPC_EXTERN PetscErrorCode STCreate_Shift(ST st)
160 : {
161 576 : PetscFunctionBegin;
162 576 : st->usesksp = PETSC_TRUE;
163 :
164 576 : st->ops->apply = STApply_Generic;
165 576 : st->ops->applytrans = STApplyTranspose_Generic;
166 576 : st->ops->applyhermtrans = STApplyHermitianTranspose_Generic;
167 576 : st->ops->backtransform = STBackTransform_Shift;
168 576 : st->ops->setshift = STSetShift_Shift;
169 576 : st->ops->getbilinearform = STGetBilinearForm_Default;
170 576 : st->ops->setup = STSetUp_Shift;
171 576 : st->ops->computeoperator = STComputeOperator_Shift;
172 576 : st->ops->postsolve = STPostSolve_Shift;
173 576 : st->ops->setdefaultksp = STSetDefaultKSP_Default;
174 576 : PetscFunctionReturn(PETSC_SUCCESS);
175 : }
|
