Actual source code: lmesetup.c

slepc-3.20.1 2023-11-27
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:    LME routines related to problem setup
 12: */

 14: #include <slepc/private/lmeimpl.h>

 16: static inline PetscErrorCode LMESetUp_Lyapunov(LME lme)
 17: {
 18:   Mat            C1,C2,X1,X2;
 19:   Vec            dc,dx;

 21:   PetscFunctionBegin;
 22:   PetscCall(MatLRCGetMats(lme->C,NULL,&C1,&dc,&C2));
 23:   PetscCheck(C1==C2,PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONGSTATE,"Lyapunov matrix equation requires symmetric right-hand side C");
 24:   PetscCheck(!dc,PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONGSTATE,"Lyapunov solvers currently require positive-definite right-hand side C");
 25:   if (lme->X) {
 26:     PetscCall(MatLRCGetMats(lme->X,NULL,&X1,&dx,&X2));
 27:     PetscCheck(X1==X2,PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONGSTATE,"Lyapunov matrix equation requires symmetric solution X");
 28:     PetscCheck(!dx,PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONGSTATE,"Lyapunov solvers currently assume a positive-definite solution X");
 29:   }
 30:   PetscFunctionReturn(PETSC_SUCCESS);
 31: }

 33: /*@
 34:    LMESetUp - Sets up all the internal data structures necessary for the
 35:    execution of the linear matrix equation solver.

 37:    Collective

 39:    Input Parameter:
 40: .  lme   - linear matrix equation solver context

 42:    Notes:
 43:    This function need not be called explicitly in most cases, since LMESolve()
 44:    calls it. It can be useful when one wants to measure the set-up time
 45:    separately from the solve time.

 47:    Level: developer

 49: .seealso: LMECreate(), LMESolve(), LMEDestroy()
 50: @*/
 51: PetscErrorCode LMESetUp(LME lme)
 52: {
 53:   PetscInt       N;

 55:   PetscFunctionBegin;

 58:   /* reset the convergence flag from the previous solves */
 59:   lme->reason = LME_CONVERGED_ITERATING;

 61:   if (lme->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
 62:   PetscCall(PetscLogEventBegin(LME_SetUp,lme,0,0,0));

 64:   /* Set default solver type (LMESetFromOptions was not called) */
 65:   if (!((PetscObject)lme)->type_name) PetscCall(LMESetType(lme,LMEKRYLOV));

 67:   /* Check problem dimensions */
 68:   PetscCheck(lme->A,PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONGSTATE,"LMESetCoefficients must be called first");
 69:   PetscCall(MatGetSize(lme->A,&N,NULL));
 70:   if (lme->ncv > N) lme->ncv = N;

 72:   /* setup options for the particular equation type */
 73:   switch (lme->problem_type) {
 74:     case LME_LYAPUNOV:
 75:       PetscCall(LMESetUp_Lyapunov(lme));
 76:       break;
 77:     case LME_SYLVESTER:
 78:       LMECheckCoeff(lme,lme->B,"B","Sylvester");
 79:       break;
 80:     case LME_GEN_LYAPUNOV:
 81:       LMECheckCoeff(lme,lme->D,"D","Generalized Lyapunov");
 82:       break;
 83:     case LME_GEN_SYLVESTER:
 84:       LMECheckCoeff(lme,lme->B,"B","Generalized Sylvester");
 85:       LMECheckCoeff(lme,lme->D,"D","Generalized Sylvester");
 86:       LMECheckCoeff(lme,lme->E,"E","Generalized Sylvester");
 87:       break;
 88:     case LME_DT_LYAPUNOV:
 89:       break;
 90:     case LME_STEIN:
 91:       LMECheckCoeff(lme,lme->D,"D","Stein");
 92:       break;
 93:   }
 94:   PetscCheck(lme->problem_type==LME_LYAPUNOV,PetscObjectComm((PetscObject)lme),PETSC_ERR_SUP,"There is no solver yet for this matrix equation type");

 96:   /* call specific solver setup */
 97:   PetscUseTypeMethod(lme,setup);

 99:   /* set tolerance if not yet set */
100:   if (lme->tol==(PetscReal)PETSC_DEFAULT) lme->tol = SLEPC_DEFAULT_TOL;

102:   PetscCall(PetscLogEventEnd(LME_SetUp,lme,0,0,0));
103:   lme->setupcalled = 1;
104:   PetscFunctionReturn(PETSC_SUCCESS);
105: }

107: static inline PetscErrorCode LMESetCoefficients_Private(LME lme,Mat A,Mat *lmeA)
108: {
109:   PetscInt       m,n;

111:   PetscFunctionBegin;
112:   PetscCall(MatGetSize(A,&m,&n));
113:   PetscCheck(m==n,PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONG,"Matrix is non-square");
114:   if (!lme->setupcalled) PetscCall(MatDestroy(lmeA));
115:   PetscCall(PetscObjectReference((PetscObject)A));
116:   *lmeA = A;
117:   PetscFunctionReturn(PETSC_SUCCESS);
118: }

120: /*@
121:    LMESetCoefficients - Sets the coefficient matrices that define the linear matrix
122:    equation to be solved.

124:    Collective

126:    Input Parameters:
127: +  lme - the matrix function context
128: .  A   - first coefficient matrix
129: .  B   - second coefficient matrix
130: .  D   - third coefficient matrix
131: -  E   - fourth coefficient matrix

133:    Notes:
134:    The matrix equation takes the general form A*X*E+D*X*B=C, where matrix C is not
135:    provided here but with LMESetRHS(). Not all four matrices must be passed, some
136:    can be NULL instead, see LMESetProblemType() for details.

138:    It must be called before LMESetUp(). If it is called again after LMESetUp() then
139:    the LME object is reset.

141:    In order to delete a previously set matrix, pass a NULL in the corresponding
142:    argument.

144:    Level: beginner

146: .seealso: LMESolve(), LMESetUp(), LMESetRHS(), LMESetProblemType()
147: @*/
148: PetscErrorCode LMESetCoefficients(LME lme,Mat A,Mat B,Mat D,Mat E)
149: {
150:   PetscFunctionBegin;
153:   PetscCheckSameComm(lme,1,A,2);
154:   if (B) {
156:     PetscCheckSameComm(lme,1,B,3);
157:   }
158:   if (D) {
160:     PetscCheckSameComm(lme,1,D,4);
161:   }
162:   if (E) {
164:     PetscCheckSameComm(lme,1,E,5);
165:   }

167:   if (lme->setupcalled) PetscCall(LMEReset(lme));

169:   PetscCall(LMESetCoefficients_Private(lme,A,&lme->A));
170:   if (B) PetscCall(LMESetCoefficients_Private(lme,B,&lme->B));
171:   else if (!lme->setupcalled) PetscCall(MatDestroy(&lme->B));
172:   if (D) PetscCall(LMESetCoefficients_Private(lme,D,&lme->D));
173:   else if (!lme->setupcalled) PetscCall(MatDestroy(&lme->D));
174:   if (E) PetscCall(LMESetCoefficients_Private(lme,E,&lme->E));
175:   else if (!lme->setupcalled) PetscCall(MatDestroy(&lme->E));

177:   lme->setupcalled = 0;
178:   PetscFunctionReturn(PETSC_SUCCESS);
179: }

181: /*@
182:    LMEGetCoefficients - Gets the coefficient matrices of the matrix equation.

184:    Not Collective

186:    Input Parameter:
187: .  lme - the LME context

189:    Output Parameters:
190: +  A   - first coefficient matrix
191: .  B   - second coefficient matrix
192: .  D   - third coefficient matrix
193: -  E   - fourth coefficient matrix

195:    Level: intermediate

197: .seealso: LMESolve(), LMESetCoefficients()
198: @*/
199: PetscErrorCode LMEGetCoefficients(LME lme,Mat *A,Mat *B,Mat *D,Mat *E)
200: {
201:   PetscFunctionBegin;
203:   if (A) *A = lme->A;
204:   if (B) *B = lme->B;
205:   if (D) *D = lme->D;
206:   if (E) *E = lme->E;
207:   PetscFunctionReturn(PETSC_SUCCESS);
208: }

210: /*@
211:    LMESetRHS - Sets the right-hand side of the matrix equation, as a low-rank
212:    matrix.

214:    Collective

216:    Input Parameters:
217: +  lme - the matrix function context
218: -  C   - the right-hand side matrix

220:    Notes:
221:    The matrix equation takes the general form A*X*E+D*X*B=C, where matrix C is
222:    given with this function. C must be a low-rank matrix of type MATLRC, that is,
223:    C = U*D*V' where D is diagonal of order k, and U, V are dense tall-skinny
224:    matrices with k columns. No sparse matrix must be provided when creating the
225:    MATLRC matrix.

227:    In equation types that require C to be symmetric, such as Lyapunov, C must be
228:    created with V=U (or V=NULL).

230:    Level: beginner

232: .seealso: LMESetSolution(), LMESetProblemType()
233: @*/
234: PetscErrorCode LMESetRHS(LME lme,Mat C)
235: {
236:   Mat            A;

238:   PetscFunctionBegin;
241:   PetscCheckSameComm(lme,1,C,2);
242:   PetscCheckTypeName(C,MATLRC);

244:   PetscCall(MatLRCGetMats(C,&A,NULL,NULL,NULL));
245:   PetscCheck(!A,PetscObjectComm((PetscObject)C),PETSC_ERR_SUP,"The MatLRC must not have a sparse matrix term");

247:   PetscCall(PetscObjectReference((PetscObject)C));
248:   PetscCall(MatDestroy(&lme->C));
249:   lme->C = C;
250:   PetscFunctionReturn(PETSC_SUCCESS);
251: }

253: /*@
254:    LMEGetRHS - Gets the right-hand side of the matrix equation.

256:    Not Collective

258:    Input Parameter:
259: .  lme - the LME context

261:    Output Parameters:
262: .  C   - the low-rank matrix

264:    Level: intermediate

266: .seealso: LMESolve(), LMESetRHS()
267: @*/
268: PetscErrorCode LMEGetRHS(LME lme,Mat *C)
269: {
270:   PetscFunctionBegin;
272:   PetscAssertPointer(C,2);
273:   *C = lme->C;
274:   PetscFunctionReturn(PETSC_SUCCESS);
275: }

277: /*@
278:    LMESetSolution - Sets the placeholder for the solution of the matrix
279:    equation, as a low-rank matrix.

281:    Collective

283:    Input Parameters:
284: +  lme - the matrix function context
285: -  X   - the solution matrix

287:    Notes:
288:    The matrix equation takes the general form A*X*E+D*X*B=C, where the solution
289:    matrix is of low rank and is written in factored form X = U*D*V'. This function
290:    provides a Mat object of type MATLRC that stores U, V and (optionally) D.
291:    These factors will be computed during LMESolve().

293:    In equation types whose solution X is symmetric, such as Lyapunov, X must be
294:    created with V=U (or V=NULL).

296:    If the user provides X with this function, then the solver will
297:    return a solution with rank at most the number of columns of U. Alternatively,
298:    it is possible to let the solver choose the rank of the solution, by
299:    setting X to NULL and then calling LMEGetSolution() after LMESolve().

301:    Level: intermediate

303: .seealso: LMEGetSolution(), LMESetRHS(), LMESetProblemType(), LMESolve()
304: @*/
305: PetscErrorCode LMESetSolution(LME lme,Mat X)
306: {
307:   Mat            A;

309:   PetscFunctionBegin;
311:   if (X) {
313:     PetscCheckSameComm(lme,1,X,2);
314:     PetscCheckTypeName(X,MATLRC);
315:     PetscCall(MatLRCGetMats(X,&A,NULL,NULL,NULL));
316:     PetscCheck(!A,PetscObjectComm((PetscObject)X),PETSC_ERR_SUP,"The MatLRC must not have a sparse matrix term");
317:     PetscCall(PetscObjectReference((PetscObject)X));
318:   }
319:   PetscCall(MatDestroy(&lme->X));
320:   lme->X = X;
321:   PetscFunctionReturn(PETSC_SUCCESS);
322: }

324: /*@
325:    LMEGetSolution - Gets the solution of the matrix equation.

327:    Not Collective

329:    Input Parameter:
330: .  lme - the LME context

332:    Output Parameters:
333: .  X   - the low-rank matrix

335:    Level: intermediate

337: .seealso: LMESolve(), LMESetSolution()
338: @*/
339: PetscErrorCode LMEGetSolution(LME lme,Mat *X)
340: {
341:   PetscFunctionBegin;
343:   PetscAssertPointer(X,2);
344:   *X = lme->X;
345:   PetscFunctionReturn(PETSC_SUCCESS);
346: }

348: /*@
349:    LMEAllocateSolution - Allocate memory storage for common variables such
350:    as the basis vectors.

352:    Collective

354:    Input Parameters:
355: +  lme   - linear matrix equation solver context
356: -  extra - number of additional positions, used for methods that require a
357:            working basis slightly larger than ncv

359:    Developer Notes:
360:    This is SLEPC_EXTERN because it may be required by user plugin LME
361:    implementations.

363:    Level: developer

365: .seealso: LMESetUp()
366: @*/
367: PetscErrorCode LMEAllocateSolution(LME lme,PetscInt extra)
368: {
369:   PetscInt       oldsize,requested;
370:   Vec            t;

372:   PetscFunctionBegin;
373:   requested = lme->ncv + extra;

375:   /* oldsize is zero if this is the first time setup is called */
376:   PetscCall(BVGetSizes(lme->V,NULL,NULL,&oldsize));

378:   /* allocate basis vectors */
379:   if (!lme->V) PetscCall(LMEGetBV(lme,&lme->V));
380:   if (!oldsize) {
381:     if (!((PetscObject)(lme->V))->type_name) PetscCall(BVSetType(lme->V,BVMAT));
382:     PetscCall(MatCreateVecsEmpty(lme->A,&t,NULL));
383:     PetscCall(BVSetSizesFromVec(lme->V,t,requested));
384:     PetscCall(VecDestroy(&t));
385:   } else PetscCall(BVResize(lme->V,requested,PETSC_FALSE));
386:   PetscFunctionReturn(PETSC_SUCCESS);
387: }