Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
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 : */
13 :
14 : #include <slepc/private/lmeimpl.h> /*I "slepclme.h" I*/
15 :
16 7 : static inline PetscErrorCode LMESetUp_Lyapunov(LME lme)
17 : {
18 7 : Mat C1,C2,X1,X2;
19 7 : Vec dc,dx;
20 :
21 7 : PetscFunctionBegin;
22 7 : PetscCall(MatLRCGetMats(lme->C,NULL,&C1,&dc,&C2));
23 7 : PetscCheck(C1==C2,PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONGSTATE,"Lyapunov matrix equation requires symmetric right-hand side C");
24 7 : PetscCheck(!dc,PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONGSTATE,"Lyapunov solvers currently require positive-definite right-hand side C");
25 7 : if (lme->X) {
26 4 : PetscCall(MatLRCGetMats(lme->X,NULL,&X1,&dx,&X2));
27 4 : PetscCheck(X1==X2,PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONGSTATE,"Lyapunov matrix equation requires symmetric solution X");
28 4 : PetscCheck(!dx,PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONGSTATE,"Lyapunov solvers currently assume a positive-definite solution X");
29 : }
30 7 : PetscFunctionReturn(PETSC_SUCCESS);
31 : }
32 :
33 : /*@
34 : LMESetUp - Sets up all the internal data structures necessary for the
35 : execution of the linear matrix equation solver.
36 :
37 : Collective
38 :
39 : Input Parameter:
40 : . lme - linear matrix equation solver context
41 :
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.
46 :
47 : Level: developer
48 :
49 : .seealso: LMECreate(), LMESolve(), LMEDestroy()
50 : @*/
51 30 : PetscErrorCode LMESetUp(LME lme)
52 : {
53 30 : PetscInt N;
54 :
55 30 : PetscFunctionBegin;
56 30 : PetscValidHeaderSpecific(lme,LME_CLASSID,1);
57 :
58 : /* reset the convergence flag from the previous solves */
59 30 : lme->reason = LME_CONVERGED_ITERATING;
60 :
61 30 : if (lme->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
62 7 : PetscCall(PetscLogEventBegin(LME_SetUp,lme,0,0,0));
63 :
64 : /* Set default solver type (LMESetFromOptions was not called) */
65 7 : if (!((PetscObject)lme)->type_name) PetscCall(LMESetType(lme,LMEKRYLOV));
66 :
67 : /* Check problem dimensions */
68 7 : PetscCheck(lme->A,PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONGSTATE,"LMESetCoefficients must be called first");
69 7 : PetscCall(MatGetSize(lme->A,&N,NULL));
70 7 : if (lme->ncv > N) lme->ncv = N;
71 :
72 : /* setup options for the particular equation type */
73 7 : switch (lme->problem_type) {
74 7 : case LME_LYAPUNOV:
75 7 : PetscCall(LMESetUp_Lyapunov(lme));
76 : break;
77 0 : case LME_SYLVESTER:
78 0 : LMECheckCoeff(lme,lme->B,"B","Sylvester");
79 : break;
80 0 : case LME_GEN_LYAPUNOV:
81 0 : LMECheckCoeff(lme,lme->D,"D","Generalized Lyapunov");
82 : break;
83 0 : case LME_GEN_SYLVESTER:
84 0 : LMECheckCoeff(lme,lme->B,"B","Generalized Sylvester");
85 0 : LMECheckCoeff(lme,lme->D,"D","Generalized Sylvester");
86 0 : LMECheckCoeff(lme,lme->E,"E","Generalized Sylvester");
87 : break;
88 : case LME_DT_LYAPUNOV:
89 : break;
90 0 : case LME_STEIN:
91 0 : LMECheckCoeff(lme,lme->D,"D","Stein");
92 : break;
93 : }
94 7 : PetscCheck(lme->problem_type==LME_LYAPUNOV,PetscObjectComm((PetscObject)lme),PETSC_ERR_SUP,"There is no solver yet for this matrix equation type");
95 :
96 : /* call specific solver setup */
97 7 : PetscUseTypeMethod(lme,setup);
98 :
99 : /* set tolerance if not yet set */
100 7 : if (lme->tol==(PetscReal)PETSC_DETERMINE) lme->tol = SLEPC_DEFAULT_TOL;
101 :
102 7 : PetscCall(PetscLogEventEnd(LME_SetUp,lme,0,0,0));
103 7 : lme->setupcalled = 1;
104 7 : PetscFunctionReturn(PETSC_SUCCESS);
105 : }
106 :
107 7 : static inline PetscErrorCode LMESetCoefficients_Private(LME lme,Mat A,Mat *lmeA)
108 : {
109 7 : PetscInt m,n;
110 :
111 7 : PetscFunctionBegin;
112 7 : PetscCall(MatGetSize(A,&m,&n));
113 7 : PetscCheck(m==n,PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONG,"Matrix is non-square");
114 7 : if (!lme->setupcalled) PetscCall(MatDestroy(lmeA));
115 7 : PetscCall(PetscObjectReference((PetscObject)A));
116 7 : *lmeA = A;
117 7 : PetscFunctionReturn(PETSC_SUCCESS);
118 : }
119 :
120 : /*@
121 : LMESetCoefficients - Sets the coefficient matrices that define the linear matrix
122 : equation to be solved.
123 :
124 : Collective
125 :
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
132 :
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.
137 :
138 : It must be called before LMESetUp(). If it is called again after LMESetUp() then
139 : the LME object is reset.
140 :
141 : In order to delete a previously set matrix, pass a NULL in the corresponding
142 : argument.
143 :
144 : Level: beginner
145 :
146 : .seealso: LMESolve(), LMESetUp(), LMESetRHS(), LMESetProblemType()
147 : @*/
148 7 : PetscErrorCode LMESetCoefficients(LME lme,Mat A,Mat B,Mat D,Mat E)
149 : {
150 7 : PetscFunctionBegin;
151 7 : PetscValidHeaderSpecific(lme,LME_CLASSID,1);
152 7 : PetscValidHeaderSpecific(A,MAT_CLASSID,2);
153 7 : PetscCheckSameComm(lme,1,A,2);
154 7 : if (B) {
155 0 : PetscValidHeaderSpecific(B,MAT_CLASSID,3);
156 0 : PetscCheckSameComm(lme,1,B,3);
157 : }
158 7 : if (D) {
159 0 : PetscValidHeaderSpecific(D,MAT_CLASSID,4);
160 0 : PetscCheckSameComm(lme,1,D,4);
161 : }
162 7 : if (E) {
163 0 : PetscValidHeaderSpecific(E,MAT_CLASSID,5);
164 0 : PetscCheckSameComm(lme,1,E,5);
165 : }
166 :
167 7 : if (lme->setupcalled) PetscCall(LMEReset(lme));
168 :
169 7 : PetscCall(LMESetCoefficients_Private(lme,A,&lme->A));
170 7 : if (B) PetscCall(LMESetCoefficients_Private(lme,B,&lme->B));
171 7 : else if (!lme->setupcalled) PetscCall(MatDestroy(&lme->B));
172 7 : if (D) PetscCall(LMESetCoefficients_Private(lme,D,&lme->D));
173 7 : else if (!lme->setupcalled) PetscCall(MatDestroy(&lme->D));
174 7 : if (E) PetscCall(LMESetCoefficients_Private(lme,E,&lme->E));
175 7 : else if (!lme->setupcalled) PetscCall(MatDestroy(&lme->E));
176 :
177 7 : lme->setupcalled = 0;
178 7 : PetscFunctionReturn(PETSC_SUCCESS);
179 : }
180 :
181 : /*@
182 : LMEGetCoefficients - Gets the coefficient matrices of the matrix equation.
183 :
184 : Not Collective
185 :
186 : Input Parameter:
187 : . lme - the LME context
188 :
189 : Output Parameters:
190 : + A - first coefficient matrix
191 : . B - second coefficient matrix
192 : . D - third coefficient matrix
193 : - E - fourth coefficient matrix
194 :
195 : Level: intermediate
196 :
197 : .seealso: LMESolve(), LMESetCoefficients()
198 : @*/
199 2 : PetscErrorCode LMEGetCoefficients(LME lme,Mat *A,Mat *B,Mat *D,Mat *E)
200 : {
201 2 : PetscFunctionBegin;
202 2 : PetscValidHeaderSpecific(lme,LME_CLASSID,1);
203 2 : if (A) *A = lme->A;
204 2 : if (B) *B = lme->B;
205 2 : if (D) *D = lme->D;
206 2 : if (E) *E = lme->E;
207 2 : PetscFunctionReturn(PETSC_SUCCESS);
208 : }
209 :
210 : /*@
211 : LMESetRHS - Sets the right-hand side of the matrix equation, as a low-rank
212 : matrix.
213 :
214 : Collective
215 :
216 : Input Parameters:
217 : + lme - the matrix function context
218 : - C - the right-hand side matrix
219 :
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.
226 :
227 : In equation types that require C to be symmetric, such as Lyapunov, C must be
228 : created with V=U (or V=NULL).
229 :
230 : Level: beginner
231 :
232 : .seealso: LMESetSolution(), LMESetProblemType()
233 : @*/
234 30 : PetscErrorCode LMESetRHS(LME lme,Mat C)
235 : {
236 30 : Mat A;
237 :
238 30 : PetscFunctionBegin;
239 30 : PetscValidHeaderSpecific(lme,LME_CLASSID,1);
240 30 : PetscValidHeaderSpecific(C,MAT_CLASSID,2);
241 30 : PetscCheckSameComm(lme,1,C,2);
242 30 : PetscCheckTypeName(C,MATLRC);
243 :
244 30 : PetscCall(MatLRCGetMats(C,&A,NULL,NULL,NULL));
245 30 : PetscCheck(!A,PetscObjectComm((PetscObject)C),PETSC_ERR_SUP,"The MatLRC must not have a sparse matrix term");
246 :
247 30 : PetscCall(PetscObjectReference((PetscObject)C));
248 30 : PetscCall(MatDestroy(&lme->C));
249 30 : lme->C = C;
250 30 : PetscFunctionReturn(PETSC_SUCCESS);
251 : }
252 :
253 : /*@
254 : LMEGetRHS - Gets the right-hand side of the matrix equation.
255 :
256 : Not Collective
257 :
258 : Input Parameter:
259 : . lme - the LME context
260 :
261 : Output Parameters:
262 : . C - the low-rank matrix
263 :
264 : Level: intermediate
265 :
266 : .seealso: LMESolve(), LMESetRHS()
267 : @*/
268 2 : PetscErrorCode LMEGetRHS(LME lme,Mat *C)
269 : {
270 2 : PetscFunctionBegin;
271 2 : PetscValidHeaderSpecific(lme,LME_CLASSID,1);
272 2 : PetscAssertPointer(C,2);
273 2 : *C = lme->C;
274 2 : PetscFunctionReturn(PETSC_SUCCESS);
275 : }
276 :
277 : /*@
278 : LMESetSolution - Sets the placeholder for the solution of the matrix
279 : equation, as a low-rank matrix.
280 :
281 : Collective
282 :
283 : Input Parameters:
284 : + lme - the matrix function context
285 : - X - the solution matrix
286 :
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().
292 :
293 : In equation types whose solution X is symmetric, such as Lyapunov, X must be
294 : created with V=U (or V=NULL).
295 :
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().
300 :
301 : Level: intermediate
302 :
303 : .seealso: LMEGetSolution(), LMESetRHS(), LMESetProblemType(), LMESolve()
304 : @*/
305 27 : PetscErrorCode LMESetSolution(LME lme,Mat X)
306 : {
307 27 : Mat A;
308 :
309 27 : PetscFunctionBegin;
310 27 : PetscValidHeaderSpecific(lme,LME_CLASSID,1);
311 27 : if (X) {
312 27 : PetscValidHeaderSpecific(X,MAT_CLASSID,2);
313 27 : PetscCheckSameComm(lme,1,X,2);
314 27 : PetscCheckTypeName(X,MATLRC);
315 27 : PetscCall(MatLRCGetMats(X,&A,NULL,NULL,NULL));
316 27 : PetscCheck(!A,PetscObjectComm((PetscObject)X),PETSC_ERR_SUP,"The MatLRC must not have a sparse matrix term");
317 27 : PetscCall(PetscObjectReference((PetscObject)X));
318 : }
319 27 : PetscCall(MatDestroy(&lme->X));
320 27 : lme->X = X;
321 27 : PetscFunctionReturn(PETSC_SUCCESS);
322 : }
323 :
324 : /*@
325 : LMEGetSolution - Gets the solution of the matrix equation.
326 :
327 : Not Collective
328 :
329 : Input Parameter:
330 : . lme - the LME context
331 :
332 : Output Parameters:
333 : . X - the low-rank matrix
334 :
335 : Level: intermediate
336 :
337 : .seealso: LMESolve(), LMESetSolution()
338 : @*/
339 1 : PetscErrorCode LMEGetSolution(LME lme,Mat *X)
340 : {
341 1 : PetscFunctionBegin;
342 1 : PetscValidHeaderSpecific(lme,LME_CLASSID,1);
343 1 : PetscAssertPointer(X,2);
344 1 : *X = lme->X;
345 1 : PetscFunctionReturn(PETSC_SUCCESS);
346 : }
347 :
348 : /*@
349 : LMEAllocateSolution - Allocate memory storage for common variables such
350 : as the basis vectors.
351 :
352 : Collective
353 :
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
358 :
359 : Developer Notes:
360 : This is SLEPC_EXTERN because it may be required by user plugin LME
361 : implementations.
362 :
363 : Level: developer
364 :
365 : .seealso: LMESetUp()
366 : @*/
367 7 : PetscErrorCode LMEAllocateSolution(LME lme,PetscInt extra)
368 : {
369 7 : PetscInt oldsize,requested;
370 7 : Vec t;
371 :
372 7 : PetscFunctionBegin;
373 7 : requested = lme->ncv + extra;
374 :
375 : /* oldsize is zero if this is the first time setup is called */
376 7 : PetscCall(BVGetSizes(lme->V,NULL,NULL,&oldsize));
377 :
378 : /* allocate basis vectors */
379 7 : if (!lme->V) PetscCall(LMEGetBV(lme,&lme->V));
380 7 : if (!oldsize) {
381 7 : if (!((PetscObject)lme->V)->type_name) PetscCall(BVSetType(lme->V,BVMAT));
382 7 : PetscCall(MatCreateVecsEmpty(lme->A,&t,NULL));
383 7 : PetscCall(BVSetSizesFromVec(lme->V,t,requested));
384 7 : PetscCall(VecDestroy(&t));
385 0 : } else PetscCall(BVResize(lme->V,requested,PETSC_FALSE));
386 7 : PetscFunctionReturn(PETSC_SUCCESS);
387 : }
|
