1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2020, 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: 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 LME374: 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: }