Actual source code: lmesetup.c
slepc-3.22.1 2024-10-28
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: PetscCall(MatGetSize(lme->A,&N,NULL));
70: if (lme->ncv > N) lme->ncv = N;
72: /* setup options for the particular equation type */
73: switch (lme->problem_type) {
74: case LME_LYAPUNOV:
75: PetscCall(LMESetUp_Lyapunov(lme));
76: break;
77: case LME_SYLVESTER:
78: LMECheckCoeff(lme,lme->B,"B","Sylvester");
79: break;
80: case LME_GEN_LYAPUNOV:
81: LMECheckCoeff(lme,lme->D,"D","Generalized Lyapunov");
82: break;
83: case LME_GEN_SYLVESTER:
84: LMECheckCoeff(lme,lme->B,"B","Generalized Sylvester");
85: LMECheckCoeff(lme,lme->D,"D","Generalized Sylvester");
86: LMECheckCoeff(lme,lme->E,"E","Generalized Sylvester");
87: break;
88: case LME_DT_LYAPUNOV:
89: break;
90: case LME_STEIN:
91: LMECheckCoeff(lme,lme->D,"D","Stein");
92: break;
93: }
94: PetscCheck(lme->problem_type==LME_LYAPUNOV,PetscObjectComm((PetscObject)lme),PETSC_ERR_SUP,"There is no solver yet for this matrix equation type");
96: /* call specific solver setup */
97: PetscUseTypeMethod(lme,setup);
99: /* set tolerance if not yet set */
100: if (lme->tol==(PetscReal)PETSC_DETERMINE) lme->tol = SLEPC_DEFAULT_TOL;
102: PetscCall(PetscLogEventEnd(LME_SetUp,lme,0,0,0));
103: lme->setupcalled = 1;
104: PetscFunctionReturn(PETSC_SUCCESS);
105: }
107: static inline PetscErrorCode LMESetCoefficients_Private(LME lme,Mat A,Mat *lmeA)
108: {
109: PetscInt m,n;
111: PetscFunctionBegin;
112: PetscCall(MatGetSize(A,&m,&n));
113: PetscCheck(m==n,PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONG,"Matrix is non-square");
114: if (!lme->setupcalled) PetscCall(MatDestroy(lmeA));
115: PetscCall(PetscObjectReference((PetscObject)A));
116: *lmeA = A;
117: PetscFunctionReturn(PETSC_SUCCESS);
118: }
120: /*@
121: LMESetCoefficients - Sets the coefficient matrices that define the linear matrix
122: equation to be solved.
124: Collective
126: Input Parameters:
127: + lme - the matrix function context
128: . A - first coefficient matrix
129: . B - second coefficient matrix
130: . D - third coefficient matrix
131: - E - fourth coefficient matrix
133: Notes:
134: The matrix equation takes the general form A*X*E+D*X*B=C, where matrix C is not
135: provided here but with LMESetRHS(). Not all four matrices must be passed, some
136: can be NULL instead, see LMESetProblemType() for details.
138: It must be called before LMESetUp(). If it is called again after LMESetUp() then
139: the LME object is reset.
141: In order to delete a previously set matrix, pass a NULL in the corresponding
142: argument.
144: Level: beginner
146: .seealso: LMESolve(), LMESetUp(), LMESetRHS(), LMESetProblemType()
147: @*/
148: PetscErrorCode LMESetCoefficients(LME lme,Mat A,Mat B,Mat D,Mat E)
149: {
150: PetscFunctionBegin;
153: PetscCheckSameComm(lme,1,A,2);
154: if (B) {
156: PetscCheckSameComm(lme,1,B,3);
157: }
158: if (D) {
160: PetscCheckSameComm(lme,1,D,4);
161: }
162: if (E) {
164: PetscCheckSameComm(lme,1,E,5);
165: }
167: if (lme->setupcalled) PetscCall(LMEReset(lme));
169: PetscCall(LMESetCoefficients_Private(lme,A,&lme->A));
170: if (B) PetscCall(LMESetCoefficients_Private(lme,B,&lme->B));
171: else if (!lme->setupcalled) PetscCall(MatDestroy(&lme->B));
172: if (D) PetscCall(LMESetCoefficients_Private(lme,D,&lme->D));
173: else if (!lme->setupcalled) PetscCall(MatDestroy(&lme->D));
174: if (E) PetscCall(LMESetCoefficients_Private(lme,E,&lme->E));
175: else if (!lme->setupcalled) PetscCall(MatDestroy(&lme->E));
177: lme->setupcalled = 0;
178: PetscFunctionReturn(PETSC_SUCCESS);
179: }
181: /*@
182: LMEGetCoefficients - Gets the coefficient matrices of the matrix equation.
184: Not Collective
186: Input Parameter:
187: . lme - the LME context
189: Output Parameters:
190: + A - first coefficient matrix
191: . B - second coefficient matrix
192: . D - third coefficient matrix
193: - E - fourth coefficient matrix
195: Level: intermediate
197: .seealso: LMESolve(), LMESetCoefficients()
198: @*/
199: PetscErrorCode LMEGetCoefficients(LME lme,Mat *A,Mat *B,Mat *D,Mat *E)
200: {
201: PetscFunctionBegin;
203: if (A) *A = lme->A;
204: if (B) *B = lme->B;
205: if (D) *D = lme->D;
206: if (E) *E = lme->E;
207: PetscFunctionReturn(PETSC_SUCCESS);
208: }
210: /*@
211: LMESetRHS - Sets the right-hand side of the matrix equation, as a low-rank
212: matrix.
214: Collective
216: Input Parameters:
217: + lme - the matrix function context
218: - C - the right-hand side matrix
220: Notes:
221: The matrix equation takes the general form A*X*E+D*X*B=C, where matrix C is
222: given with this function. C must be a low-rank matrix of type MATLRC, that is,
223: C = U*D*V' where D is diagonal of order k, and U, V are dense tall-skinny
224: matrices with k columns. No sparse matrix must be provided when creating the
225: MATLRC matrix.
227: In equation types that require C to be symmetric, such as Lyapunov, C must be
228: created with V=U (or V=NULL).
230: Level: beginner
232: .seealso: LMESetSolution(), LMESetProblemType()
233: @*/
234: PetscErrorCode LMESetRHS(LME lme,Mat C)
235: {
236: Mat A;
238: PetscFunctionBegin;
241: PetscCheckSameComm(lme,1,C,2);
242: PetscCheckTypeName(C,MATLRC);
244: PetscCall(MatLRCGetMats(C,&A,NULL,NULL,NULL));
245: PetscCheck(!A,PetscObjectComm((PetscObject)C),PETSC_ERR_SUP,"The MatLRC must not have a sparse matrix term");
247: PetscCall(PetscObjectReference((PetscObject)C));
248: PetscCall(MatDestroy(&lme->C));
249: lme->C = C;
250: PetscFunctionReturn(PETSC_SUCCESS);
251: }
253: /*@
254: LMEGetRHS - Gets the right-hand side of the matrix equation.
256: Not Collective
258: Input Parameter:
259: . lme - the LME context
261: Output Parameters:
262: . C - the low-rank matrix
264: Level: intermediate
266: .seealso: LMESolve(), LMESetRHS()
267: @*/
268: PetscErrorCode LMEGetRHS(LME lme,Mat *C)
269: {
270: PetscFunctionBegin;
272: PetscAssertPointer(C,2);
273: *C = lme->C;
274: PetscFunctionReturn(PETSC_SUCCESS);
275: }
277: /*@
278: LMESetSolution - Sets the placeholder for the solution of the matrix
279: equation, as a low-rank matrix.
281: Collective
283: Input Parameters:
284: + lme - the matrix function context
285: - X - the solution matrix
287: Notes:
288: The matrix equation takes the general form A*X*E+D*X*B=C, where the solution
289: matrix is of low rank and is written in factored form X = U*D*V'. This function
290: provides a Mat object of type MATLRC that stores U, V and (optionally) D.
291: These factors will be computed during LMESolve().
293: In equation types whose solution X is symmetric, such as Lyapunov, X must be
294: created with V=U (or V=NULL).
296: If the user provides X with this function, then the solver will
297: return a solution with rank at most the number of columns of U. Alternatively,
298: it is possible to let the solver choose the rank of the solution, by
299: setting X to NULL and then calling LMEGetSolution() after LMESolve().
301: Level: intermediate
303: .seealso: LMEGetSolution(), LMESetRHS(), LMESetProblemType(), LMESolve()
304: @*/
305: PetscErrorCode LMESetSolution(LME lme,Mat X)
306: {
307: Mat A;
309: PetscFunctionBegin;
311: if (X) {
313: PetscCheckSameComm(lme,1,X,2);
314: PetscCheckTypeName(X,MATLRC);
315: PetscCall(MatLRCGetMats(X,&A,NULL,NULL,NULL));
316: PetscCheck(!A,PetscObjectComm((PetscObject)X),PETSC_ERR_SUP,"The MatLRC must not have a sparse matrix term");
317: PetscCall(PetscObjectReference((PetscObject)X));
318: }
319: PetscCall(MatDestroy(&lme->X));
320: lme->X = X;
321: PetscFunctionReturn(PETSC_SUCCESS);
322: }
324: /*@
325: LMEGetSolution - Gets the solution of the matrix equation.
327: Not Collective
329: Input Parameter:
330: . lme - the LME context
332: Output Parameters:
333: . X - the low-rank matrix
335: Level: intermediate
337: .seealso: LMESolve(), LMESetSolution()
338: @*/
339: PetscErrorCode LMEGetSolution(LME lme,Mat *X)
340: {
341: PetscFunctionBegin;
343: PetscAssertPointer(X,2);
344: *X = lme->X;
345: PetscFunctionReturn(PETSC_SUCCESS);
346: }
348: /*@
349: LMEAllocateSolution - Allocate memory storage for common variables such
350: as the basis vectors.
352: Collective
354: Input Parameters:
355: + lme - linear matrix equation solver context
356: - extra - number of additional positions, used for methods that require a
357: working basis slightly larger than ncv
359: Developer Notes:
360: This is SLEPC_EXTERN because it may be required by user plugin LME
361: implementations.
363: Level: developer
365: .seealso: LMESetUp()
366: @*/
367: PetscErrorCode LMEAllocateSolution(LME lme,PetscInt extra)
368: {
369: PetscInt oldsize,requested;
370: Vec t;
372: PetscFunctionBegin;
373: requested = lme->ncv + extra;
375: /* oldsize is zero if this is the first time setup is called */
376: PetscCall(BVGetSizes(lme->V,NULL,NULL,&oldsize));
378: /* allocate basis vectors */
379: if (!lme->V) PetscCall(LMEGetBV(lme,&lme->V));
380: if (!oldsize) {
381: if (!((PetscObject)lme->V)->type_name) PetscCall(BVSetType(lme->V,BVMAT));
382: PetscCall(MatCreateVecsEmpty(lme->A,&t,NULL));
383: PetscCall(BVSetSizesFromVec(lme->V,t,requested));
384: PetscCall(VecDestroy(&t));
385: } else PetscCall(BVResize(lme->V,requested,PETSC_FALSE));
386: PetscFunctionReturn(PETSC_SUCCESS);
387: }