Actual source code: lmesetup.c

slepc-3.12.1 2019-11-08
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2019, 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>       /*I "slepclme.h" I*/

 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:   PetscBool      match;
245:   Mat            A;


252:   PetscObjectTypeCompare((PetscObject)C,MATLRC,&match);
253:   if (!match) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_SUP,"Mat argument must have been created with MatCreateLRC");
254:   MatLRCGetMats(C,&A,NULL,NULL,NULL);
255:   if (A) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_SUP,"The MatLRC must not have a sparse matrix term");

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

263: /*@
264:    LMEGetRHS - Gets the right-hand side of the matrix equation.

266:    Collective on lme

268:    Input Parameter:
269: .  lme - the LME context

271:    Output Parameters:
272: .  C   - the low-rank matrix

274:    Level: intermediate

276: .seealso: LMESolve(), LMESetRHS()
277: @*/
278: PetscErrorCode LMEGetRHS(LME lme,Mat *C)
279: {
283:   *C = lme->C;
284:   return(0);
285: }

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

291:    Collective on lme

293:    Input Parameters:
294: +  lme - the matrix function context
295: -  X   - the solution matrix

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

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

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

311:    Level: intermediate

313: .seealso: LMEGetSolution(), LMESetRHS(), LMESetProblemType(), LMESolve()
314: @*/
315: PetscErrorCode LMESetSolution(LME lme,Mat X)
316: {
318:   PetscBool      match;
319:   Mat            A;

323:   if (X) {
326:     PetscObjectTypeCompare((PetscObject)X,MATLRC,&match);
327:     if (!match) SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_SUP,"Mat argument must have been created with MatCreateLRC");
328:     MatLRCGetMats(X,&A,NULL,NULL,NULL);
329:     if (A) SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_SUP,"The MatLRC must not have a sparse matrix term");
330:     PetscObjectReference((PetscObject)X);
331:   }
332:   MatDestroy(&lme->X);
333:   lme->X = X;
334:   return(0);
335: }

337: /*@
338:    LMEGetSolution - Gets the solution of the matrix equation.

340:    Collective on lme

342:    Input Parameter:
343: .  lme - the LME context

345:    Output Parameters:
346: .  X   - the low-rank matrix

348:    Level: intermediate

350: .seealso: LMESolve(), LMESetSolution()
351: @*/
352: PetscErrorCode LMEGetSolution(LME lme,Mat *X)
353: {
357:   *X = lme->X;
358:   return(0);
359: }

361: /*@
362:    LMEAllocateSolution - Allocate memory storage for common variables such
363:    as the basis vectors.

365:    Collective on lme

367:    Input Parameters:
368: +  lme   - linear matrix equation solver context
369: -  extra - number of additional positions, used for methods that require a
370:            working basis slightly larger than ncv

372:    Developers Note:
373:    This is SLEPC_EXTERN because it may be required by user plugin LME
374:    implementations.

376:    Level: developer
377: @*/
378: PetscErrorCode LMEAllocateSolution(LME lme,PetscInt extra)
379: {
381:   PetscInt       oldsize,requested;
382:   Vec            t;

385:   requested = lme->ncv + extra;

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

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