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 - the 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: [](ch:lme), `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->E,"E","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 linear matrix equation solver 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 $AXE+DXB=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 `NULL` in the corresponding
143: argument.
145: Level: beginner
147: .seealso: [](ch:lme), `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 linear matrix equation solver 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: [](ch:lme), `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 linear matrix equation solver context
219: - C - the right-hand side matrix
221: Notes:
222: The matrix equation takes the general form $AXE+DXB=C$, where matrix $C$ is
223: given with this function. $C$ must be a low-rank matrix of type `MATLRC`, that is,
224: $C = UDV^*$ 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: [](ch:lme), `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 linear matrix equation solver context
262: Output Parameter:
263: . C - the right-hand side matrix
265: Level: intermediate
267: .seealso: [](ch:lme), `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 linear matrix equation solver context
286: - X - the solution matrix
288: Notes:
289: The matrix equation takes the general form $AXE+DXB=C$, where the solution
290: matrix is of low rank and is written in factored form $X = UDV^*$. 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: [](ch:lme), `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 linear matrix equation solver context
333: Output Parameter:
334: . X - the low-rank matrix
336: Notes:
337: If called after `LMESolve()`, `X` will contain the solution of the equation.
339: The matrix `X` may have been passed by the user via `LMESetSolution()`,
340: although this is not required.
342: Level: beginner
344: .seealso: [](ch:lme), `LMESolve()`, `LMESetSolution()`
345: @*/
346: PetscErrorCode LMEGetSolution(LME lme,Mat *X)
347: {
348: PetscFunctionBegin;
350: PetscAssertPointer(X,2);
351: *X = lme->X;
352: PetscFunctionReturn(PETSC_SUCCESS);
353: }
355: /*@
356: LMEAllocateSolution - Allocate memory storage for common variables such
357: as the basis vectors.
359: Collective
361: Input Parameters:
362: + lme - the linear matrix equation solver context
363: - extra - number of additional positions, used for methods that require a
364: working basis slightly larger than `ncv`
366: Developer Notes:
367: This is `SLEPC_EXTERN` because it may be required by user plugin `LME`
368: implementations.
370: Level: developer
372: .seealso: [](ch:lme), `LMESetUp()`, `LMESetDimensions()`
373: @*/
374: PetscErrorCode LMEAllocateSolution(LME lme,PetscInt extra)
375: {
376: PetscInt oldsize,requested;
377: Vec t;
379: PetscFunctionBegin;
380: requested = lme->ncv + extra;
382: /* oldsize is zero if this is the first time setup is called */
383: PetscCall(BVGetSizes(lme->V,NULL,NULL,&oldsize));
385: /* allocate basis vectors */
386: if (!lme->V) PetscCall(LMEGetBV(lme,&lme->V));
387: if (!oldsize) {
388: if (!((PetscObject)lme->V)->type_name) PetscCall(BVSetType(lme->V,BVMAT));
389: PetscCall(MatCreateVecsEmpty(lme->A,&t,NULL));
390: PetscCall(BVSetSizesFromVec(lme->V,t,requested));
391: PetscCall(VecDestroy(&t));
392: } else PetscCall(BVResize(lme->V,requested,PETSC_FALSE));
393: PetscFunctionReturn(PETSC_SUCCESS);
394: }