1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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> /*I "slepclme.h" I*/
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 requires 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 assumes a positive-definite solution X");
30: }
31: return(0);
32: }
34: PETSC_STATIC_INLINE PetscErrorCode LMESetUp_Sylvester(LME lme) 35: {
37: if (!lme->B) SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONGSTATE,"Sylvester matrix equation requires coefficient matrix B");
38: SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_SUP,"There is no solver yet for this matrix equation type");
39: return(0);
40: }
42: PETSC_STATIC_INLINE PetscErrorCode LMESetUp_Gen_Lyapunov(LME lme) 43: {
45: if (!lme->D) SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONGSTATE,"Generalized Lyapunov matrix equation requires coefficient matrix D");
46: SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_SUP,"There is no solver yet for this matrix equation type");
47: return(0);
48: }
50: PETSC_STATIC_INLINE PetscErrorCode LMESetUp_Gen_Sylvester(LME lme) 51: {
53: if (!lme->B) SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONGSTATE,"Generalized Sylvester matrix equation requires coefficient matrix B");
54: if (!lme->D) SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONGSTATE,"Generalized Sylvester matrix equation requires coefficient matrix D");
55: if (!lme->E) SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONGSTATE,"Generalized Sylvester matrix equation requires coefficient matrix E");
56: SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_SUP,"There is no solver yet for this matrix equation type");
57: return(0);
58: }
60: PETSC_STATIC_INLINE PetscErrorCode LMESetUp_Stein(LME lme) 61: {
63: if (!lme->D) SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONGSTATE,"Stein matrix equation requires coefficient matrix D");
64: SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_SUP,"There is no solver yet for this matrix equation type");
65: return(0);
66: }
68: PETSC_STATIC_INLINE PetscErrorCode LMESetUp_DT_Lyapunov(LME lme) 69: {
71: SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_SUP,"There is no solver yet for this matrix equation type");
72: return(0);
73: }
75: /*@
76: LMESetUp - Sets up all the internal data structures necessary for the
77: execution of the linear matrix equation solver.
79: Collective on LME 81: Input Parameter:
82: . lme - linear matrix equation solver context
84: Notes:
85: This function need not be called explicitly in most cases, since LMESolve()
86: calls it. It can be useful when one wants to measure the set-up time
87: separately from the solve time.
89: Level: developer
91: .seealso: LMECreate(), LMESolve(), LMEDestroy()
92: @*/
93: PetscErrorCode LMESetUp(LME lme) 94: {
96: PetscInt N;
101: /* reset the convergence flag from the previous solves */
102: lme->reason = LME_CONVERGED_ITERATING;
104: if (lme->setupcalled) return(0);
105: PetscLogEventBegin(LME_SetUp,lme,0,0,0);
107: /* Set default solver type (LMESetFromOptions was not called) */
108: if (!((PetscObject)lme)->type_name) {
109: LMESetType(lme,LMEKRYLOV);
110: }
112: /* Check problem dimensions */
113: if (!lme->A) SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONGSTATE,"LMESetCoefficients must be called first");
114: MatGetSize(lme->A,&N,NULL);
115: if (lme->ncv > N) lme->ncv = N;
117: /* setup options for the particular equation type */
118: switch (lme->problem_type) {
119: case LME_LYAPUNOV:
120: LMESetUp_Lyapunov(lme);
121: break;
122: case LME_SYLVESTER:
123: LMESetUp_Sylvester(lme);
124: break;
125: case LME_GEN_LYAPUNOV:
126: LMESetUp_Gen_Lyapunov(lme);
127: break;
128: case LME_GEN_SYLVESTER:
129: LMESetUp_Gen_Sylvester(lme);
130: break;
131: case LME_DT_LYAPUNOV:
132: LMESetUp_DT_Lyapunov(lme);
133: break;
134: case LME_STEIN:
135: LMESetUp_Stein(lme);
136: break;
137: }
139: /* call specific solver setup */
140: (*lme->ops->setup)(lme);
142: /* set tolerance if not yet set */
143: if (lme->tol==PETSC_DEFAULT) lme->tol = SLEPC_DEFAULT_TOL;
145: PetscLogEventEnd(LME_SetUp,lme,0,0,0);
146: lme->setupcalled = 1;
147: return(0);
148: }
150: /*@
151: LMESetCoefficients - Sets the coefficient matrices that define the linear matrix
152: equation to be solved.
154: Collective on LME and Mat
156: Input Parameters:
157: + lme - the matrix function context
158: . A - first coefficient matrix
159: . B - second coefficient matrix
160: . D - third coefficient matrix
161: - E - fourth coefficient matrix
163: Notes:
164: The matrix equation takes the general form A*X*E+D*X*B=C, where matrix C is not
165: provided here but with LMESetRHS(). Not all four matrices must be passed, some
166: can be NULL instead, see LMESetProblemType() for details.
168: It must be called before LMESetUp(). If it is called again after LMESetUp() then
169: the LME object is reset.
171: In order to delete a previously set matrix, pass a NULL in the corresponding
172: argument.
174: Level: beginner
176: .seealso: LMESolve(), LMESetUp(), LMESetRHS(), LMESetProblemType()
177: @*/
178: PetscErrorCode LMESetCoefficients(LME lme,Mat A,Mat B,Mat D,Mat E)179: {
181: PetscInt m,n;
187: if (B) {
190: }
191: if (D) {
194: }
195: if (E) {
198: }
200: if (lme->setupcalled) { LMEReset(lme); }
202: MatGetSize(A,&m,&n);
203: if (m!=n) SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONG,"A is a non-square matrix");
204: if (!lme->setupcalled) { MatDestroy(&lme->A); }
205: PetscObjectReference((PetscObject)A);
206: lme->A = A;
207: if (B) {
208: MatGetSize(B,&m,&n);
209: if (m!=n) SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONG,"B is a non-square matrix");
210: if (!lme->setupcalled) { MatDestroy(&lme->B); }
211: PetscObjectReference((PetscObject)B);
212: lme->B = B;
213: } else if (!lme->setupcalled) { MatDestroy(&lme->B); }
214: if (D) {
215: MatGetSize(D,&m,&n);
216: if (m!=n) SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONG,"D is a non-square matrix");
217: if (!lme->setupcalled) { MatDestroy(&lme->D); }
218: PetscObjectReference((PetscObject)D);
219: lme->D = D;
220: } else if (!lme->setupcalled) { MatDestroy(&lme->D); }
221: if (E) {
222: MatGetSize(E,&m,&n);
223: if (m!=n) SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONG,"E is a non-square matrix");
224: if (!lme->setupcalled) { MatDestroy(&lme->E); }
225: PetscObjectReference((PetscObject)E);
226: lme->E = E;
227: } else if (!lme->setupcalled) { MatDestroy(&lme->E); }
229: lme->setupcalled = 0;
230: return(0);
231: }
233: /*@
234: LMEGetCoefficients - Gets the coefficient matrices of the matrix equation.
236: Collective on LME and Mat
238: Input Parameter:
239: . lme - the LME context
241: Output Parameters:
242: + A - first coefficient matrix
243: . B - second coefficient matrix
244: . D - third coefficient matrix
245: - E - fourth coefficient matrix
247: Level: intermediate
249: .seealso: LMESolve(), LMESetCoefficients()
250: @*/
251: PetscErrorCode LMEGetCoefficients(LME lme,Mat *A,Mat *B,Mat *D,Mat *E)252: {
255: if (A) *A = lme->A;
256: if (B) *B = lme->B;
257: if (D) *D = lme->D;
258: if (E) *E = lme->E;
259: return(0);
260: }
262: /*@
263: LMESetRHS - Sets the right-hand side of the matrix equation, as a low-rank
264: matrix.
266: Collective on LME and Mat
268: Input Parameters:
269: + lme - the matrix function context
270: - C - the right-hand side matrix
272: Notes:
273: The matrix equation takes the general form A*X*E+D*X*B=C, where matrix C is
274: given with this function. C must be a low-rank matrix of type MATLRC, that is,
275: C = U*D*V' where D is diagonal of order k, and U, V are dense tall-skinny
276: matrices with k columns. No sparse matrix must be provided when creating the
277: MATLRC matrix.
279: In equation types that require C to be symmetric, such as Lyapunov, C must be
280: created with V=U (or V=NULL).
282: Level: beginner
284: .seealso: LMESetSolution(), LMESetProblemType()
285: @*/
286: PetscErrorCode LMESetRHS(LME lme,Mat C)287: {
289: PetscBool match;
290: Mat A;
297: PetscObjectTypeCompare((PetscObject)C,MATLRC,&match);
298: if (!match) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_SUP,"Mat argument must have been created with MatCreateLRC");
299: MatLRCGetMats(C,&A,NULL,NULL,NULL);
300: if (A) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_SUP,"The MatLRC must not have a sparse matrix term");
302: PetscObjectReference((PetscObject)C);
303: MatDestroy(&lme->C);
304: lme->C = C;
305: return(0);
306: }
308: /*@
309: LMEGetRHS - Gets the right-hand side of the matrix equation.
311: Collective on LME and Mat
313: Input Parameter:
314: . lme - the LME context
316: Output Parameters:
317: . C - the low-rank matrix
319: Level: intermediate
321: .seealso: LMESolve(), LMESetRHS()
322: @*/
323: PetscErrorCode LMEGetRHS(LME lme,Mat *C)324: {
328: *C = lme->C;
329: return(0);
330: }
332: /*@
333: LMESetSolution - Sets the placeholder for the solution of the matrix
334: equation, as a low-rank matrix.
336: Collective on LME and Mat
338: Input Parameters:
339: + lme - the matrix function context
340: - X - the solution matrix
342: Notes:
343: The matrix equation takes the general form A*X*E+D*X*B=C, where the solution
344: matrix is of low rank and is written in factored form X = U*D*V'. This function
345: provides a Mat object of type MATLRC that stores U, V and (optionally) D.
346: These factors will be computed during LMESolve().
348: In equation types whose solution X is symmetric, such as Lyapunov, X must be
349: created with V=U (or V=NULL).
351: If the user provides X with this function, then the solver will
352: return a solution with rank at most the number of columns of U. Alternatively,
353: it is possible to let the solver choose the rank of the solution, by
354: setting X to NULL and then calling LMEGetSolution() after LMESolve().
356: Level: intermediate
358: .seealso: LMEGetSolution(), LMESetRHS(), LMESetProblemType(), LMESolve()
359: @*/
360: PetscErrorCode LMESetSolution(LME lme,Mat X)361: {
363: PetscBool match;
364: Mat A;
368: if (X) {
371: PetscObjectTypeCompare((PetscObject)X,MATLRC,&match);
372: if (!match) SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_SUP,"Mat argument must have been created with MatCreateLRC");
373: MatLRCGetMats(X,&A,NULL,NULL,NULL);
374: if (A) SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_SUP,"The MatLRC must not have a sparse matrix term");
375: PetscObjectReference((PetscObject)X);
376: }
377: MatDestroy(&lme->X);
378: lme->X = X;
379: return(0);
380: }
382: /*@
383: LMEGetSolution - Gets the solution of the matrix equation.
385: Collective on LME and Mat
387: Input Parameter:
388: . lme - the LME context
390: Output Parameters:
391: . X - the low-rank matrix
393: Level: intermediate
395: .seealso: LMESolve(), LMESetSolution()
396: @*/
397: PetscErrorCode LMEGetSolution(LME lme,Mat *X)398: {
402: *X = lme->X;
403: return(0);
404: }
406: /*@
407: LMEAllocateSolution - Allocate memory storage for common variables such
408: as the basis vectors.
410: Collective on LME412: Input Parameters:
413: + lme - linear matrix equation solver context
414: - extra - number of additional positions, used for methods that require a
415: working basis slightly larger than ncv
417: Developers Note:
418: This is SLEPC_EXTERN because it may be required by user plugin LME419: implementations.
421: Level: developer
422: @*/
423: PetscErrorCode LMEAllocateSolution(LME lme,PetscInt extra)424: {
426: PetscInt oldsize,requested;
427: Vec t;
430: requested = lme->ncv + extra;
432: /* oldsize is zero if this is the first time setup is called */
433: BVGetSizes(lme->V,NULL,NULL,&oldsize);
435: /* allocate basis vectors */
436: if (!lme->V) { LMEGetBV(lme,&lme->V); }
437: if (!oldsize) {
438: if (!((PetscObject)(lme->V))->type_name) {
439: BVSetType(lme->V,BVSVEC);
440: }
441: MatCreateVecsEmpty(lme->A,&t,NULL);
442: BVSetSizesFromVec(lme->V,t,requested);
443: VecDestroy(&t);
444: } else {
445: BVResize(lme->V,requested,PETSC_FALSE);
446: }
447: return(0);
448: }