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: }