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 5505 : static PetscErrorCode STApply_Shift_BSE(ST st,Vec x,Vec y)
33 : {
34 5505 : Mat H,R,C;
35 5505 : Vec x1,x2,y1;
36 :
37 5505 : PetscFunctionBegin;
38 5505 : H = st->T[0];
39 5505 : PetscCall(MatNestGetSubMat(H,0,0,&R));
40 5505 : PetscCall(MatNestGetSubMat(H,0,1,&C));
41 5505 : PetscCall(VecNestGetSubVec(x,0,&x1));
42 5505 : PetscCall(VecNestGetSubVec(x,1,&x2));
43 5505 : PetscCall(VecNestGetSubVec(y,0,&y1));
44 5505 : PetscCall(MatMult(C,x2,y1));
45 5505 : PetscCall(MatMultAdd(R,x1,y1,y1));
46 5505 : PetscFunctionReturn(PETSC_SUCCESS);
47 : }
48 :
49 1495919 : static PetscErrorCode STBackTransform_Shift(ST st,PetscInt n,PetscScalar *eigr,PetscScalar *eigi)
50 : {
51 1495919 : PetscInt j;
52 :
53 1495919 : PetscFunctionBegin;
54 4376837 : for (j=0;j<n;j++) {
55 2880918 : eigr[j] += st->sigma;
56 : }
57 1495919 : PetscFunctionReturn(PETSC_SUCCESS);
58 : }
59 :
60 590 : static PetscErrorCode STPostSolve_Shift(ST st)
61 : {
62 590 : PetscFunctionBegin;
63 590 : if (st->matmode == ST_MATMODE_INPLACE) {
64 3 : if (st->nmat>1) PetscCall(MatAXPY(st->A[0],st->sigma,st->A[1],st->str));
65 2 : else PetscCall(MatShift(st->A[0],st->sigma));
66 3 : st->Astate[0] = ((PetscObject)st->A[0])->state;
67 3 : st->state = ST_STATE_INITIAL;
68 3 : st->opready = PETSC_FALSE;
69 : }
70 590 : 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 447 : static PetscErrorCode STComputeOperator_Shift(ST st)
80 : {
81 447 : PetscFunctionBegin;
82 447 : st->usesksp = (st->nmat>1)? PETSC_TRUE: PETSC_FALSE;
83 447 : PetscCall(PetscObjectReference((PetscObject)st->A[1]));
84 447 : PetscCall(MatDestroy(&st->T[1]));
85 447 : st->T[1] = st->A[1];
86 447 : PetscCall(STMatMAXPY_Private(st,-st->sigma,0.0,0,NULL,PetscNot(st->state==ST_STATE_UPDATED),PETSC_FALSE,&st->T[0]));
87 447 : if (st->nmat>1) PetscCall(PetscObjectReference((PetscObject)st->T[1]));
88 447 : PetscCall(MatDestroy(&st->P));
89 447 : st->P = (st->nmat>1)? st->T[1]: NULL;
90 447 : st->M = st->T[0];
91 447 : 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 447 : PetscFunctionReturn(PETSC_SUCCESS);
97 : }
98 :
99 515 : static PetscErrorCode STSetUp_Shift(ST st)
100 : {
101 515 : PetscInt k,nc,nmat=st->nmat;
102 515 : PetscScalar *coeffs=NULL;
103 :
104 515 : PetscFunctionBegin;
105 515 : if (nmat>1) PetscCall(STSetWorkVecs(st,1));
106 515 : if (nmat>2) { /* set-up matrices for polynomial eigenproblems */
107 68 : 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 224 : for (k=0;k<nmat;k++) {
128 170 : PetscCall(PetscObjectReference((PetscObject)st->A[k]));
129 170 : PetscCall(MatDestroy(&st->T[k]));
130 170 : st->T[k] = st->A[k];
131 : }
132 : }
133 : }
134 515 : if (st->P) PetscCall(KSPSetUp(st->ksp));
135 515 : if (st->structured) st->ops->apply = STApply_Shift_BSE;
136 515 : PetscFunctionReturn(PETSC_SUCCESS);
137 : }
138 :
139 15 : static PetscErrorCode STSetShift_Shift(ST st,PetscScalar newshift)
140 : {
141 15 : PetscInt k,nc,nmat=PetscMax(st->nmat,2);
142 15 : PetscScalar *coeffs=NULL;
143 :
144 15 : PetscFunctionBegin;
145 15 : if (st->transform) {
146 13 : 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 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]));
153 13 : if (st->matmode == ST_MATMODE_COPY && nmat>2) PetscCall(PetscFree(coeffs));
154 13 : if (st->nmat<=2) st->M = st->T[0];
155 : }
156 15 : PetscFunctionReturn(PETSC_SUCCESS);
157 : }
158 :
159 588 : SLEPC_EXTERN PetscErrorCode STCreate_Shift(ST st)
160 : {
161 588 : PetscFunctionBegin;
162 588 : st->usesksp = PETSC_TRUE;
163 :
164 588 : st->ops->apply = STApply_Generic;
165 588 : st->ops->applytrans = STApplyTranspose_Generic;
166 588 : st->ops->applyhermtrans = STApplyHermitianTranspose_Generic;
167 588 : st->ops->backtransform = STBackTransform_Shift;
168 588 : st->ops->setshift = STSetShift_Shift;
169 588 : st->ops->getbilinearform = STGetBilinearForm_Default;
170 588 : st->ops->setup = STSetUp_Shift;
171 588 : st->ops->computeoperator = STComputeOperator_Shift;
172 588 : st->ops->postsolve = STPostSolve_Shift;
173 588 : st->ops->setdefaultksp = STSetDefaultKSP_Default;
174 588 : PetscFunctionReturn(PETSC_SUCCESS);
175 : }
|
