Actual source code: lmesetup.c

slepc-3.11.2 2019-07-30
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 requires 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 assumes a positive-definite solution X");
 30:   }
 31:   return(0);
 32: }

 34: PETSC_STATIC_INLINE PetscErrorCode LMESetUp_Sylvester(LME lme)
 35: {
 37:   if (!lme->B) SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONGSTATE,"Sylvester matrix equation requires coefficient matrix B");
 38:   SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_SUP,"There is no solver yet for this matrix equation type");
 39:   return(0);
 40: }

 42: PETSC_STATIC_INLINE PetscErrorCode LMESetUp_Gen_Lyapunov(LME lme)
 43: {
 45:   if (!lme->D) SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONGSTATE,"Generalized Lyapunov matrix equation requires coefficient matrix D");
 46:   SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_SUP,"There is no solver yet for this matrix equation type");
 47:   return(0);
 48: }

 50: PETSC_STATIC_INLINE PetscErrorCode LMESetUp_Gen_Sylvester(LME lme)
 51: {
 53:   if (!lme->B) SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONGSTATE,"Generalized Sylvester matrix equation requires coefficient matrix B");
 54:   if (!lme->D) SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONGSTATE,"Generalized Sylvester matrix equation requires coefficient matrix D");
 55:   if (!lme->E) SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONGSTATE,"Generalized Sylvester matrix equation requires coefficient matrix E");
 56:   SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_SUP,"There is no solver yet for this matrix equation type");
 57:   return(0);
 58: }

 60: PETSC_STATIC_INLINE PetscErrorCode LMESetUp_Stein(LME lme)
 61: {
 63:   if (!lme->D) SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONGSTATE,"Stein matrix equation requires coefficient matrix D");
 64:   SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_SUP,"There is no solver yet for this matrix equation type");
 65:   return(0);
 66: }

 68: PETSC_STATIC_INLINE PetscErrorCode LMESetUp_DT_Lyapunov(LME lme)
 69: {
 71:   SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_SUP,"There is no solver yet for this matrix equation type");
 72:   return(0);
 73: }

 75: /*@
 76:    LMESetUp - Sets up all the internal data structures necessary for the
 77:    execution of the linear matrix equation solver.

 79:    Collective on LME

 81:    Input Parameter:
 82: .  lme   - linear matrix equation solver context

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

 89:    Level: developer

 91: .seealso: LMECreate(), LMESolve(), LMEDestroy()
 92: @*/
 93: PetscErrorCode LMESetUp(LME lme)
 94: {
 96:   PetscInt       N;


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

104:   if (lme->setupcalled) return(0);
105:   PetscLogEventBegin(LME_SetUp,lme,0,0,0);

107:   /* Set default solver type (LMESetFromOptions was not called) */
108:   if (!((PetscObject)lme)->type_name) {
109:     LMESetType(lme,LMEKRYLOV);
110:   }

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

117:   /* setup options for the particular equation type */
118:   switch (lme->problem_type) {
119:     case LME_LYAPUNOV:
120:       LMESetUp_Lyapunov(lme);
121:       break;
122:     case LME_SYLVESTER:
123:       LMESetUp_Sylvester(lme);
124:       break;
125:     case LME_GEN_LYAPUNOV:
126:       LMESetUp_Gen_Lyapunov(lme);
127:       break;
128:     case LME_GEN_SYLVESTER:
129:       LMESetUp_Gen_Sylvester(lme);
130:       break;
131:     case LME_DT_LYAPUNOV:
132:       LMESetUp_DT_Lyapunov(lme);
133:       break;
134:     case LME_STEIN:
135:       LMESetUp_Stein(lme);
136:       break;
137:   }

139:   /* call specific solver setup */
140:   (*lme->ops->setup)(lme);

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

145:   PetscLogEventEnd(LME_SetUp,lme,0,0,0);
146:   lme->setupcalled = 1;
147:   return(0);
148: }

150: /*@
151:    LMESetCoefficients - Sets the coefficient matrices that define the linear matrix
152:    equation to be solved.

154:    Collective on LME and Mat

156:    Input Parameters:
157: +  lme - the matrix function context
158: .  A   - first coefficient matrix
159: .  B   - second coefficient matrix
160: .  D   - third coefficient matrix
161: -  E   - fourth coefficient matrix

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

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

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

174:    Level: beginner

176: .seealso: LMESolve(), LMESetUp(), LMESetRHS(), LMESetProblemType()
177: @*/
178: PetscErrorCode LMESetCoefficients(LME lme,Mat A,Mat B,Mat D,Mat E)
179: {
181:   PetscInt       m,n;

187:   if (B) {
190:   }
191:   if (D) {
194:   }
195:   if (E) {
198:   }

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

202:   MatGetSize(A,&m,&n);
203:   if (m!=n) SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONG,"A is a non-square matrix");
204:   if (!lme->setupcalled) { MatDestroy(&lme->A); }
205:   PetscObjectReference((PetscObject)A);
206:   lme->A = A;
207:   if (B) {
208:     MatGetSize(B,&m,&n);
209:     if (m!=n) SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONG,"B is a non-square matrix");
210:     if (!lme->setupcalled) { MatDestroy(&lme->B); }
211:     PetscObjectReference((PetscObject)B);
212:     lme->B = B;
213:   } else if (!lme->setupcalled) { MatDestroy(&lme->B); }
214:   if (D) {
215:     MatGetSize(D,&m,&n);
216:     if (m!=n) SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONG,"D is a non-square matrix");
217:     if (!lme->setupcalled) { MatDestroy(&lme->D); }
218:     PetscObjectReference((PetscObject)D);
219:     lme->D = D;
220:   } else if (!lme->setupcalled) { MatDestroy(&lme->D); }
221:   if (E) {
222:     MatGetSize(E,&m,&n);
223:     if (m!=n) SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONG,"E is a non-square matrix");
224:     if (!lme->setupcalled) { MatDestroy(&lme->E); }
225:     PetscObjectReference((PetscObject)E);
226:     lme->E = E;
227:   } else if (!lme->setupcalled) { MatDestroy(&lme->E); }

229:   lme->setupcalled = 0;
230:   return(0);
231: }

233: /*@
234:    LMEGetCoefficients - Gets the coefficient matrices of the matrix equation.

236:    Collective on LME and Mat

238:    Input Parameter:
239: .  lme - the LME context

241:    Output Parameters:
242: +  A   - first coefficient matrix
243: .  B   - second coefficient matrix
244: .  D   - third coefficient matrix
245: -  E   - fourth coefficient matrix

247:    Level: intermediate

249: .seealso: LMESolve(), LMESetCoefficients()
250: @*/
251: PetscErrorCode LMEGetCoefficients(LME lme,Mat *A,Mat *B,Mat *D,Mat *E)
252: {
255:   if (A) *A = lme->A;
256:   if (B) *B = lme->B;
257:   if (D) *D = lme->D;
258:   if (E) *E = lme->E;
259:   return(0);
260: }

262: /*@
263:    LMESetRHS - Sets the right-hand side of the matrix equation, as a low-rank
264:    matrix.

266:    Collective on LME and Mat

268:    Input Parameters:
269: +  lme - the matrix function context
270: -  C   - the right-hand side matrix

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

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

282:    Level: beginner

284: .seealso: LMESetSolution(), LMESetProblemType()
285: @*/
286: PetscErrorCode LMESetRHS(LME lme,Mat C)
287: {
289:   PetscBool      match;
290:   Mat            A;


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

302:   PetscObjectReference((PetscObject)C);
303:   MatDestroy(&lme->C);
304:   lme->C = C;
305:   return(0);
306: }

308: /*@
309:    LMEGetRHS - Gets the right-hand side of the matrix equation.

311:    Collective on LME and Mat

313:    Input Parameter:
314: .  lme - the LME context

316:    Output Parameters:
317: .  C   - the low-rank matrix

319:    Level: intermediate

321: .seealso: LMESolve(), LMESetRHS()
322: @*/
323: PetscErrorCode LMEGetRHS(LME lme,Mat *C)
324: {
328:   *C = lme->C;
329:   return(0);
330: }

332: /*@
333:    LMESetSolution - Sets the placeholder for the solution of the matrix
334:    equation, as a low-rank matrix.

336:    Collective on LME and Mat

338:    Input Parameters:
339: +  lme - the matrix function context
340: -  X   - the solution matrix

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

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

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

356:    Level: intermediate

358: .seealso: LMEGetSolution(), LMESetRHS(), LMESetProblemType(), LMESolve()
359: @*/
360: PetscErrorCode LMESetSolution(LME lme,Mat X)
361: {
363:   PetscBool      match;
364:   Mat            A;

368:   if (X) {
371:     PetscObjectTypeCompare((PetscObject)X,MATLRC,&match);
372:     if (!match) SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_SUP,"Mat argument must have been created with MatCreateLRC");
373:     MatLRCGetMats(X,&A,NULL,NULL,NULL);
374:     if (A) SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_SUP,"The MatLRC must not have a sparse matrix term");
375:     PetscObjectReference((PetscObject)X);
376:   }
377:   MatDestroy(&lme->X);
378:   lme->X = X;
379:   return(0);
380: }

382: /*@
383:    LMEGetSolution - Gets the solution of the matrix equation.

385:    Collective on LME and Mat

387:    Input Parameter:
388: .  lme - the LME context

390:    Output Parameters:
391: .  X   - the low-rank matrix

393:    Level: intermediate

395: .seealso: LMESolve(), LMESetSolution()
396: @*/
397: PetscErrorCode LMEGetSolution(LME lme,Mat *X)
398: {
402:   *X = lme->X;
403:   return(0);
404: }

406: /*@
407:    LMEAllocateSolution - Allocate memory storage for common variables such
408:    as the basis vectors.

410:    Collective on LME

412:    Input Parameters:
413: +  lme   - linear matrix equation solver context
414: -  extra - number of additional positions, used for methods that require a
415:            working basis slightly larger than ncv

417:    Developers Note:
418:    This is SLEPC_EXTERN because it may be required by user plugin LME
419:    implementations.

421:    Level: developer
422: @*/
423: PetscErrorCode LMEAllocateSolution(LME lme,PetscInt extra)
424: {
426:   PetscInt       oldsize,requested;
427:   Vec            t;

430:   requested = lme->ncv + extra;

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

435:   /* allocate basis vectors */
436:   if (!lme->V) { LMEGetBV(lme,&lme->V); }
437:   if (!oldsize) {
438:     if (!((PetscObject)(lme->V))->type_name) {
439:       BVSetType(lme->V,BVSVEC);
440:     }
441:     MatCreateVecsEmpty(lme->A,&t,NULL);
442:     BVSetSizesFromVec(lme->V,t,requested);
443:     VecDestroy(&t);
444:   } else {
445:     BVResize(lme->V,requested,PETSC_FALSE);
446:   }
447:   return(0);
448: }