Actual source code: lmesetup.c

  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:   PetscCheck(lme->C,PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONGSTATE,"LMESetRHS must be called first");
 70:   PetscCall(MatGetSize(lme->A,&N,NULL));
 71:   if (lme->ncv > N) lme->ncv = N;

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

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

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

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

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

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

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

125:    Collective

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

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

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

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

145:    Level: beginner

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

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

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

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

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

185:    Collective

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

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

196:    Level: intermediate

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

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

215:    Collective

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

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

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

231:    Level: beginner

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

239:   PetscFunctionBegin;
242:   PetscCheckSameComm(lme,1,C,2);
243:   PetscCheckTypeName(C,MATLRC);

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

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

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

257:    Collective

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

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

265:    Level: intermediate

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

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

282:    Collective

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

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

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

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

302:    Level: intermediate

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

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

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

328:    Collective

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

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

336:    Level: intermediate

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

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

353:    Collective

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

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

364:    Level: developer

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

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

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

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