Actual source code: lmesetup.c

slepc-3.15.1 2021-05-28
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2021, 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: PETSC_STATIC_INLINE PetscErrorCode LMESetUp_Lyapunov(LME lme)
 17: {
 19:   Mat            C1,C2,X1,X2;
 20:   Vec            dc,dx;

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

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

 38:    Collective on lme

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

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

 48:    Level: developer

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


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

 63:   if (lme->setupcalled) return(0);
 64:   PetscLogEventBegin(LME_SetUp,lme,0,0,0);

 66:   /* Set default solver type (LMESetFromOptions was not called) */
 67:   if (!((PetscObject)lme)->type_name) {
 68:     LMESetType(lme,LMEKRYLOV);
 69:   }

 71:   /* Check problem dimensions */
 72:   if (!lme->A) SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONGSTATE,"LMESetCoefficients must be called first");
 73:   MatGetSize(lme->A,&N,NULL);
 74:   if (lme->ncv > N) lme->ncv = N;

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

100:   /* call specific solver setup */
101:   (*lme->ops->setup)(lme);

103:   /* set tolerance if not yet set */
104:   if (lme->tol==PETSC_DEFAULT) lme->tol = SLEPC_DEFAULT_TOL;

106:   PetscLogEventEnd(LME_SetUp,lme,0,0,0);
107:   lme->setupcalled = 1;
108:   return(0);
109: }

111: PETSC_STATIC_INLINE PetscErrorCode LMESetCoefficients_Private(LME lme,Mat A,Mat *lmeA)
112: {
114:   PetscInt       m,n;

117:   MatGetSize(A,&m,&n);
118:   if (m!=n) SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONG,"Matrix is non-square");
119:   if (!lme->setupcalled) { MatDestroy(lmeA); }
120:   PetscObjectReference((PetscObject)A);
121:   *lmeA = A;
122:   return(0);
123: }

125: /*@
126:    LMESetCoefficients - Sets the coefficient matrices that define the linear matrix
127:    equation to be solved.

129:    Collective on lme

131:    Input Parameters:
132: +  lme - the matrix function context
133: .  A   - first coefficient matrix
134: .  B   - second coefficient matrix
135: .  D   - third coefficient matrix
136: -  E   - fourth coefficient matrix

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

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

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

149:    Level: beginner

151: .seealso: LMESolve(), LMESetUp(), LMESetRHS(), LMESetProblemType()
152: @*/
153: PetscErrorCode LMESetCoefficients(LME lme,Mat A,Mat B,Mat D,Mat E)
154: {

161:   if (B) {
164:   }
165:   if (D) {
168:   }
169:   if (E) {
172:   }

174:   if (lme->setupcalled) { LMEReset(lme); }

176:   LMESetCoefficients_Private(lme,A,&lme->A);
177:   if (B) { LMESetCoefficients_Private(lme,B,&lme->B); }
178:   else if (!lme->setupcalled) { MatDestroy(&lme->B); }
179:   if (D) { LMESetCoefficients_Private(lme,D,&lme->D); }
180:   else if (!lme->setupcalled) { MatDestroy(&lme->D); }
181:   if (E) { LMESetCoefficients_Private(lme,E,&lme->E); }
182:   else if (!lme->setupcalled) { MatDestroy(&lme->E); }

184:   lme->setupcalled = 0;
185:   return(0);
186: }

188: /*@
189:    LMEGetCoefficients - Gets the coefficient matrices of the matrix equation.

191:    Collective on lme

193:    Input Parameter:
194: .  lme - the LME context

196:    Output Parameters:
197: +  A   - first coefficient matrix
198: .  B   - second coefficient matrix
199: .  D   - third coefficient matrix
200: -  E   - fourth coefficient matrix

202:    Level: intermediate

204: .seealso: LMESolve(), LMESetCoefficients()
205: @*/
206: PetscErrorCode LMEGetCoefficients(LME lme,Mat *A,Mat *B,Mat *D,Mat *E)
207: {
210:   if (A) *A = lme->A;
211:   if (B) *B = lme->B;
212:   if (D) *D = lme->D;
213:   if (E) *E = lme->E;
214:   return(0);
215: }

217: /*@
218:    LMESetRHS - Sets the right-hand side of the matrix equation, as a low-rank
219:    matrix.

221:    Collective on lme

223:    Input Parameters:
224: +  lme - the matrix function context
225: -  C   - the right-hand side matrix

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

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

237:    Level: beginner

239: .seealso: LMESetSolution(), LMESetProblemType()
240: @*/
241: PetscErrorCode LMESetRHS(LME lme,Mat C)
242: {
244:   Mat            A;


252:   MatLRCGetMats(C,&A,NULL,NULL,NULL);
253:   if (A) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_SUP,"The MatLRC must not have a sparse matrix term");

255:   PetscObjectReference((PetscObject)C);
256:   MatDestroy(&lme->C);
257:   lme->C = C;
258:   return(0);
259: }

261: /*@
262:    LMEGetRHS - Gets the right-hand side of the matrix equation.

264:    Collective on lme

266:    Input Parameter:
267: .  lme - the LME context

269:    Output Parameters:
270: .  C   - the low-rank matrix

272:    Level: intermediate

274: .seealso: LMESolve(), LMESetRHS()
275: @*/
276: PetscErrorCode LMEGetRHS(LME lme,Mat *C)
277: {
281:   *C = lme->C;
282:   return(0);
283: }

285: /*@
286:    LMESetSolution - Sets the placeholder for the solution of the matrix
287:    equation, as a low-rank matrix.

289:    Collective on lme

291:    Input Parameters:
292: +  lme - the matrix function context
293: -  X   - the solution matrix

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

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

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

309:    Level: intermediate

311: .seealso: LMEGetSolution(), LMESetRHS(), LMESetProblemType(), LMESolve()
312: @*/
313: PetscErrorCode LMESetSolution(LME lme,Mat X)
314: {
316:   Mat            A;

320:   if (X) {
324:     MatLRCGetMats(X,&A,NULL,NULL,NULL);
325:     if (A) SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_SUP,"The MatLRC must not have a sparse matrix term");
326:     PetscObjectReference((PetscObject)X);
327:   }
328:   MatDestroy(&lme->X);
329:   lme->X = X;
330:   return(0);
331: }

333: /*@
334:    LMEGetSolution - Gets the solution of the matrix equation.

336:    Collective on lme

338:    Input Parameter:
339: .  lme - the LME context

341:    Output Parameters:
342: .  X   - the low-rank matrix

344:    Level: intermediate

346: .seealso: LMESolve(), LMESetSolution()
347: @*/
348: PetscErrorCode LMEGetSolution(LME lme,Mat *X)
349: {
353:   *X = lme->X;
354:   return(0);
355: }

357: /*@
358:    LMEAllocateSolution - Allocate memory storage for common variables such
359:    as the basis vectors.

361:    Collective on lme

363:    Input Parameters:
364: +  lme   - linear matrix equation solver context
365: -  extra - number of additional positions, used for methods that require a
366:            working basis slightly larger than ncv

368:    Developers Note:
369:    This is SLEPC_EXTERN because it may be required by user plugin LME
370:    implementations.

372:    Level: developer
373: @*/
374: PetscErrorCode LMEAllocateSolution(LME lme,PetscInt extra)
375: {
377:   PetscInt       oldsize,requested;
378:   Vec            t;

381:   requested = lme->ncv + extra;

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

386:   /* allocate basis vectors */
387:   if (!lme->V) { LMEGetBV(lme,&lme->V); }
388:   if (!oldsize) {
389:     if (!((PetscObject)(lme->V))->type_name) {
390:       BVSetType(lme->V,BVSVEC);
391:     }
392:     MatCreateVecsEmpty(lme->A,&t,NULL);
393:     BVSetSizesFromVec(lme->V,t,requested);
394:     VecDestroy(&t);
395:   } else {
396:     BVResize(lme->V,requested,PETSC_FALSE);
397:   }
398:   return(0);
399: }