Actual source code: shift.c

slepc-3.23.1 2025-05-01
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain

  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: */

 15: #include <slepc/private/stimpl.h>
 16: #include <petsc/private/matimpl.h>  /* Mat_MPIDense */

 18: /*
 19:    Special STApply() for the BSE structured matrix

 21:        H = [ R  C; -C^H -R^T ].

 23:    Assumes that H is a MATNEST and x,y are VECNEST.

 25:    Only the upper part of the product y=H*x is computed.

 27:        y1 = R*x1+C*x2

 29:    The bottom part of the input vector x2 is computed as
 30:    either conj(x1) or -conj(x1), where the sign is given by
 31:    s in the context SlepcMatStruct.
 32:    The bottom part of the output vector y2 is not referenced.
 33: */
 34: static PetscErrorCode STApply_Shift_BSE(ST st,Vec x,Vec y)
 35: {
 36:   Mat            H,R,C;
 37:   Vec            x1,x2,y1;
 38:   PetscContainer container;
 39:   SlepcMatStruct mctx;

 41:   PetscFunctionBegin;
 42:   H = st->T[0];
 43:   PetscCall(PetscObjectQuery((PetscObject)H,"SlepcMatStruct",(PetscObject*)&container));
 44:   PetscCall(PetscContainerGetPointer(container,(void**)&mctx));
 45:   PetscCall(MatNestGetSubMat(H,0,0,&R));
 46:   PetscCall(MatNestGetSubMat(H,0,1,&C));
 47:   PetscCall(VecNestGetSubVec(x,0,&x1));
 48:   PetscCall(VecNestGetSubVec(x,1,&x2));
 49:   PetscCall(VecNestGetSubVec(y,0,&y1));

 51:   /* x2 = +/-conj(x1) */
 52:   PetscCall(VecCopy(x1,x2));
 53:   PetscCall(VecConjugate(x2));
 54:   if (mctx->s==-1.0) PetscCall(VecScale(x2,-1.0));

 56:   PetscCall(MatMult(C,x2,y1));
 57:   PetscCall(MatMultAdd(R,x1,y1,y1));
 58:   PetscFunctionReturn(PETSC_SUCCESS);
 59: }

 61: /*
 62:    Specialized version that avoids communication for multiplication by R.
 63:    It needs to access internal data structures of MATMPIDENSE.
 64: */
 65: static PetscErrorCode STApply_Shift_BSE_Dense(ST st,Vec x,Vec y)
 66: {
 67:   Mat            H,R,C;
 68:   Vec            x1,x2,y1;
 69:   Mat_MPIDense   *mdnR,*mdnC;
 70:   PetscContainer container;
 71:   SlepcMatStruct mctx;

 73:   PetscFunctionBegin;
 74:   H = st->T[0];
 75:   PetscCall(PetscObjectQuery((PetscObject)H,"SlepcMatStruct",(PetscObject*)&container));
 76:   PetscCall(PetscContainerGetPointer(container,(void**)&mctx));
 77:   PetscCall(MatNestGetSubMat(H,0,0,&R));
 78:   PetscCall(MatNestGetSubMat(H,0,1,&C));
 79:   mdnR = (Mat_MPIDense*)R->data;
 80:   mdnC = (Mat_MPIDense*)C->data;
 81:   PetscCall(VecNestGetSubVec(x,0,&x1));
 82:   PetscCall(VecNestGetSubVec(x,1,&x2));
 83:   PetscCall(VecNestGetSubVec(y,0,&y1));

 85:   /* x2 = +/-conj(x1) */
 86:   PetscCall(VecCopy(x1,x2));
 87:   PetscCall(VecConjugate(x2));
 88:   if (mctx->s==-1.0) PetscCall(VecScale(x2,-1.0));

 90:   PetscCall(MatMult(C,x2,y1));
 91:   /* PetscCall(MatMultAdd(R,x1,y1,y1)); */
 92:   PetscCall(VecConjugate(mdnC->lvec));
 93:   if (mctx->s==-1.0) PetscCall(VecScale(mdnC->lvec,-1.0));
 94:   PetscCall(PetscLogEventBegin(MAT_MultAdd,R,mdnC->lvec,y1,y1));
 95:   PetscCall(VecLockReadPush(mdnC->lvec));
 96:   PetscCall((*mdnR->A->ops->multadd)(mdnR->A,mdnC->lvec,y1,y1));
 97:   PetscCall(VecLockReadPop(mdnC->lvec));
 98:   PetscCall(PetscLogEventEnd(MAT_MultAdd,R,mdnC->lvec,y1,y1));
 99:   PetscCall(PetscObjectStateIncrease((PetscObject)y1));
100:   PetscFunctionReturn(PETSC_SUCCESS);
101: }

103: static PetscErrorCode STBackTransform_Shift(ST st,PetscInt n,PetscScalar *eigr,PetscScalar *eigi)
104: {
105:   PetscInt j;

107:   PetscFunctionBegin;
108:   for (j=0;j<n;j++) {
109:     eigr[j] += st->sigma;
110:   }
111:   PetscFunctionReturn(PETSC_SUCCESS);
112: }

114: static PetscErrorCode STPostSolve_Shift(ST st)
115: {
116:   PetscFunctionBegin;
117:   if (st->matmode == ST_MATMODE_INPLACE) {
118:     if (st->nmat>1) PetscCall(MatAXPY(st->A[0],st->sigma,st->A[1],st->str));
119:     else PetscCall(MatShift(st->A[0],st->sigma));
120:     st->Astate[0] = ((PetscObject)st->A[0])->state;
121:     st->state   = ST_STATE_INITIAL;
122:     st->opready = PETSC_FALSE;
123:   }
124:   PetscFunctionReturn(PETSC_SUCCESS);
125: }

127: /*
128:    Operator (shift):
129:                Op               P         M
130:    if nmat=1:  A-sI             NULL      A-sI
131:    if nmat=2:  B^-1 (A-sB)      B         A-sB
132: */
133: static PetscErrorCode STComputeOperator_Shift(ST st)
134: {
135:   PetscFunctionBegin;
136:   st->usesksp = (st->nmat>1)? PETSC_TRUE: PETSC_FALSE;
137:   PetscCall(PetscObjectReference((PetscObject)st->A[1]));
138:   PetscCall(MatDestroy(&st->T[1]));
139:   st->T[1] = st->A[1];
140:   PetscCall(STMatMAXPY_Private(st,-st->sigma,0.0,0,NULL,PetscNot(st->state==ST_STATE_UPDATED),PETSC_FALSE,&st->T[0]));
141:   if (st->nmat>1) PetscCall(PetscObjectReference((PetscObject)st->T[1]));
142:   PetscCall(MatDestroy(&st->P));
143:   st->P = (st->nmat>1)? st->T[1]: NULL;
144:   st->M = st->T[0];
145:   if (st->nmat>1 && st->Psplit) {  /* build custom preconditioner from the split matrices */
146:     PetscCall(MatDestroy(&st->Pmat));
147:     PetscCall(PetscObjectReference((PetscObject)st->Psplit[1]));
148:     st->Pmat = st->Psplit[1];
149:   }
150:   PetscFunctionReturn(PETSC_SUCCESS);
151: }

153: static PetscErrorCode STSetUp_Shift(ST st)
154: {
155:   PetscInt       k,nc,nmat=st->nmat;
156:   PetscScalar    *coeffs=NULL;
157:   PetscBool      denseR,denseC;
158:   Mat            H,R,C;

160:   PetscFunctionBegin;
161:   if (nmat>1) PetscCall(STSetWorkVecs(st,1));
162:   if (nmat>2) {  /* set-up matrices for polynomial eigenproblems */
163:     if (st->transform) {
164:       nc = (nmat*(nmat+1))/2;
165:       PetscCall(PetscMalloc1(nc,&coeffs));
166:       /* Compute coeffs */
167:       PetscCall(STCoeffs_Monomial(st,coeffs));
168:       /* T[n] = A_n */
169:       k = nmat-1;
170:       PetscCall(PetscObjectReference((PetscObject)st->A[k]));
171:       PetscCall(MatDestroy(&st->T[k]));
172:       st->T[k] = st->A[k];
173:       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]));
174:       PetscCall(PetscFree(coeffs));
175:       PetscCall(PetscObjectReference((PetscObject)st->T[nmat-1]));
176:       PetscCall(MatDestroy(&st->P));
177:       st->P = st->T[nmat-1];
178:       if (st->Psplit) {  /* build custom preconditioner from the split matrices */
179:         PetscCall(STMatMAXPY_Private(st,st->sigma,0.0,nmat-1,coeffs?coeffs:NULL,PETSC_TRUE,PETSC_TRUE,&st->Pmat));
180:       }
181:       PetscCall(ST_KSPSetOperators(st,st->P,st->Pmat?st->Pmat:st->P));
182:     } else {
183:       for (k=0;k<nmat;k++) {
184:         PetscCall(PetscObjectReference((PetscObject)st->A[k]));
185:         PetscCall(MatDestroy(&st->T[k]));
186:         st->T[k] = st->A[k];
187:       }
188:     }
189:   }
190:   if (st->P) PetscCall(KSPSetUp(st->ksp));
191:   if (st->structured) {
192:     H = st->T[0];
193:     PetscCall(MatNestGetSubMat(H,0,0,&R));
194:     PetscCall(MatNestGetSubMat(H,0,1,&C));
195:     PetscCall(PetscObjectTypeCompareAny((PetscObject)R,&denseR,MATMPIDENSE,MATMPIDENSECUDA,MATMPIDENSEHIP,""));
196:     PetscCall(PetscObjectTypeCompareAny((PetscObject)C,&denseC,MATMPIDENSE,MATMPIDENSECUDA,MATMPIDENSEHIP,""));
197:     st->ops->apply = (denseR && denseC)? STApply_Shift_BSE_Dense: STApply_Shift_BSE;
198:   }
199:   PetscFunctionReturn(PETSC_SUCCESS);
200: }

202: static PetscErrorCode STSetShift_Shift(ST st,PetscScalar newshift)
203: {
204:   PetscInt       k,nc,nmat=PetscMax(st->nmat,2);
205:   PetscScalar    *coeffs=NULL;

207:   PetscFunctionBegin;
208:   if (st->transform) {
209:     if (st->matmode == ST_MATMODE_COPY && nmat>2) {
210:       nc = (nmat*(nmat+1))/2;
211:       PetscCall(PetscMalloc1(nc,&coeffs));
212:       /* Compute coeffs */
213:       PetscCall(STCoeffs_Monomial(st,coeffs));
214:     }
215:     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]));
216:     if (st->matmode == ST_MATMODE_COPY && nmat>2) PetscCall(PetscFree(coeffs));
217:     if (st->nmat<=2) st->M = st->T[0];
218:   }
219:   PetscFunctionReturn(PETSC_SUCCESS);
220: }

222: SLEPC_EXTERN PetscErrorCode STCreate_Shift(ST st)
223: {
224:   PetscFunctionBegin;
225:   st->usesksp = PETSC_TRUE;

227:   st->ops->apply           = STApply_Generic;
228:   st->ops->applytrans      = STApplyTranspose_Generic;
229:   st->ops->applyhermtrans  = STApplyHermitianTranspose_Generic;
230:   st->ops->backtransform   = STBackTransform_Shift;
231:   st->ops->setshift        = STSetShift_Shift;
232:   st->ops->getbilinearform = STGetBilinearForm_Default;
233:   st->ops->setup           = STSetUp_Shift;
234:   st->ops->computeoperator = STComputeOperator_Shift;
235:   st->ops->postsolve       = STPostSolve_Shift;
236:   st->ops->setdefaultksp   = STSetDefaultKSP_Default;
237:   PetscFunctionReturn(PETSC_SUCCESS);
238: }