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 the solution process
12 : */
13 :
14 : #include <slepc/private/lmeimpl.h> /*I "slepclme.h" I*/
15 : #include <slepcblaslapack.h>
16 :
17 : /*@
18 : LMESolve - Solves the linear matrix equation.
19 :
20 : Collective
21 :
22 : Input Parameter:
23 : . lme - linear matrix equation solver context obtained from LMECreate()
24 :
25 : Options Database Keys:
26 : + -lme_view - print information about the solver used
27 : . -lme_view_mat binary - save the matrix to the default binary viewer
28 : . -lme_view_rhs binary - save right hand side to the default binary viewer
29 : . -lme_view_solution binary - save computed solution to the default binary viewer
30 : - -lme_converged_reason - print reason for convergence, and number of iterations
31 :
32 : Notes:
33 : The matrix coefficients are specified with LMESetCoefficients().
34 : The right-hand side is specified with LMESetRHS().
35 : The placeholder for the solution is specified with LMESetSolution().
36 :
37 : Level: beginner
38 :
39 : .seealso: LMECreate(), LMESetUp(), LMEDestroy(), LMESetTolerances(), LMESetCoefficients(), LMESetRHS(), LMESetSolution()
40 : @*/
41 30 : PetscErrorCode LMESolve(LME lme)
42 : {
43 30 : PetscFunctionBegin;
44 30 : PetscValidHeaderSpecific(lme,LME_CLASSID,1);
45 :
46 : /* call setup */
47 30 : PetscCall(LMESetUp(lme));
48 30 : lme->its = 0;
49 30 : lme->errest = 0.0;
50 :
51 30 : PetscCall(LMEViewFromOptions(lme,NULL,"-lme_view_pre"));
52 :
53 : /* call solver */
54 30 : PetscCheck(lme->ops->solve[lme->problem_type],PetscObjectComm((PetscObject)lme),PETSC_ERR_SUP,"The specified solver does not support equation type %s",LMEProblemTypes[lme->problem_type]);
55 30 : PetscCall(PetscLogEventBegin(LME_Solve,lme,0,0,0));
56 30 : PetscUseTypeMethod(lme,solve[lme->problem_type]);
57 30 : PetscCall(PetscLogEventEnd(LME_Solve,lme,0,0,0));
58 :
59 30 : PetscCheck(lme->reason,PetscObjectComm((PetscObject)lme),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
60 :
61 30 : PetscCheck(!lme->errorifnotconverged || lme->reason>=0,PetscObjectComm((PetscObject)lme),PETSC_ERR_NOT_CONVERGED,"LMESolve has not converged");
62 :
63 : /* various viewers */
64 30 : PetscCall(LMEViewFromOptions(lme,NULL,"-lme_view"));
65 30 : PetscCall(LMEConvergedReasonViewFromOptions(lme));
66 30 : PetscCall(MatViewFromOptions(lme->A,(PetscObject)lme,"-lme_view_mat"));
67 30 : PetscCall(MatViewFromOptions(lme->C,(PetscObject)lme,"-lme_view_rhs"));
68 30 : PetscCall(MatViewFromOptions(lme->X,(PetscObject)lme,"-lme_view_solution"));
69 30 : PetscFunctionReturn(PETSC_SUCCESS);
70 : }
71 :
72 : /*@
73 : LMEGetIterationNumber - Gets the current iteration number. If the
74 : call to LMESolve() is complete, then it returns the number of iterations
75 : carried out by the solution method.
76 :
77 : Not Collective
78 :
79 : Input Parameter:
80 : . lme - the linear matrix equation solver context
81 :
82 : Output Parameter:
83 : . its - number of iterations
84 :
85 : Note:
86 : During the i-th iteration this call returns i-1. If LMESolve() is
87 : complete, then parameter "its" contains either the iteration number at
88 : which convergence was successfully reached, or failure was detected.
89 : Call LMEGetConvergedReason() to determine if the solver converged or
90 : failed and why.
91 :
92 : Level: intermediate
93 :
94 : .seealso: LMEGetConvergedReason(), LMESetTolerances()
95 : @*/
96 2 : PetscErrorCode LMEGetIterationNumber(LME lme,PetscInt *its)
97 : {
98 2 : PetscFunctionBegin;
99 2 : PetscValidHeaderSpecific(lme,LME_CLASSID,1);
100 2 : PetscAssertPointer(its,2);
101 2 : *its = lme->its;
102 2 : PetscFunctionReturn(PETSC_SUCCESS);
103 : }
104 :
105 : /*@
106 : LMEGetConvergedReason - Gets the reason why the LMESolve() iteration was
107 : stopped.
108 :
109 : Not Collective
110 :
111 : Input Parameter:
112 : . lme - the linear matrix equation solver context
113 :
114 : Output Parameter:
115 : . reason - negative value indicates diverged, positive value converged
116 :
117 : Notes:
118 :
119 : Possible values for reason are
120 : + LME_CONVERGED_TOL - converged up to tolerance
121 : . LME_DIVERGED_ITS - required more than max_it iterations to reach convergence
122 : - LME_DIVERGED_BREAKDOWN - generic breakdown in method
123 :
124 : Can only be called after the call to LMESolve() is complete.
125 :
126 : Level: intermediate
127 :
128 : .seealso: LMESetTolerances(), LMESolve(), LMEConvergedReason, LMESetErrorIfNotConverged()
129 : @*/
130 2 : PetscErrorCode LMEGetConvergedReason(LME lme,LMEConvergedReason *reason)
131 : {
132 2 : PetscFunctionBegin;
133 2 : PetscValidHeaderSpecific(lme,LME_CLASSID,1);
134 2 : PetscAssertPointer(reason,2);
135 2 : *reason = lme->reason;
136 2 : PetscFunctionReturn(PETSC_SUCCESS);
137 : }
138 :
139 : /*@
140 : LMEGetErrorEstimate - Returns the error estimate obtained during solve.
141 :
142 : Not Collective
143 :
144 : Input Parameter:
145 : . lme - linear matrix equation solver context
146 :
147 : Output Parameter:
148 : . errest - the error estimate
149 :
150 : Notes:
151 : This is the error estimated internally by the solver. The actual
152 : error bound can be computed with LMEComputeError(). Note that some
153 : solvers may not be able to provide an error estimate.
154 :
155 : Level: advanced
156 :
157 : .seealso: LMEComputeError()
158 : @*/
159 4 : PetscErrorCode LMEGetErrorEstimate(LME lme,PetscReal *errest)
160 : {
161 4 : PetscFunctionBegin;
162 4 : PetscValidHeaderSpecific(lme,LME_CLASSID,1);
163 4 : PetscAssertPointer(errest,2);
164 4 : *errest = lme->errest;
165 4 : PetscFunctionReturn(PETSC_SUCCESS);
166 : }
167 :
168 : /*
169 : LMEComputeResidualNorm_Lyapunov - Computes the Frobenius norm of the residual matrix
170 : associated with the Lyapunov equation.
171 : */
172 4 : static PetscErrorCode LMEComputeResidualNorm_Lyapunov(LME lme,PetscReal *norm)
173 : {
174 4 : PetscInt j,n,N,k,l,lda,ldb;
175 4 : PetscBLASInt n_,N_,k_,l_,lda_,ldb_;
176 4 : PetscScalar *Rarray,alpha=1.0,beta=0.0;
177 4 : const PetscScalar *A,*B;
178 4 : BV W,AX,X1,C1;
179 4 : Mat R,X1m,C1m;
180 4 : Vec v,w;
181 4 : VecScatter vscat;
182 :
183 4 : PetscFunctionBegin;
184 4 : PetscCall(MatLRCGetMats(lme->C,NULL,&C1m,NULL,NULL));
185 4 : PetscCall(MatLRCGetMats(lme->X,NULL,&X1m,NULL,NULL));
186 4 : PetscCall(BVCreateFromMat(C1m,&C1));
187 4 : PetscCall(BVSetFromOptions(C1));
188 4 : PetscCall(BVCreateFromMat(X1m,&X1));
189 4 : PetscCall(BVSetFromOptions(X1));
190 4 : PetscCall(BVGetSizes(X1,&n,&N,&k));
191 4 : PetscCall(BVGetSizes(C1,NULL,NULL,&l));
192 4 : PetscCall(PetscBLASIntCast(n,&n_));
193 4 : PetscCall(PetscBLASIntCast(N,&N_));
194 4 : PetscCall(PetscBLASIntCast(k,&k_));
195 4 : PetscCall(PetscBLASIntCast(l,&l_));
196 :
197 : /* create W to store a redundant copy of a BV in each process */
198 4 : PetscCall(BVCreate(PETSC_COMM_SELF,&W));
199 4 : PetscCall(BVSetSizes(W,N,N,k));
200 4 : PetscCall(BVSetFromOptions(W));
201 4 : PetscCall(BVGetColumn(X1,0,&v));
202 4 : PetscCall(VecScatterCreateToAll(v,&vscat,NULL));
203 4 : PetscCall(BVRestoreColumn(X1,0,&v));
204 :
205 : /* create AX to hold the product A*X1 */
206 4 : PetscCall(BVDuplicate(X1,&AX));
207 4 : PetscCall(BVMatMult(X1,lme->A,AX));
208 :
209 : /* create dense matrix to hold the residual R=C1*C1'+AX*X1'+X1*AX' */
210 4 : PetscCall(MatCreateDense(PetscObjectComm((PetscObject)lme),n,n,N,N,NULL,&R));
211 :
212 : /* R=C1*C1' */
213 4 : PetscCall(MatDenseGetArrayWrite(R,&Rarray));
214 12 : for (j=0;j<l;j++) {
215 8 : PetscCall(BVGetColumn(C1,j,&v));
216 8 : PetscCall(BVGetColumn(W,j,&w));
217 8 : PetscCall(VecScatterBegin(vscat,v,w,INSERT_VALUES,SCATTER_FORWARD));
218 8 : PetscCall(VecScatterEnd(vscat,v,w,INSERT_VALUES,SCATTER_FORWARD));
219 8 : PetscCall(BVRestoreColumn(C1,j,&v));
220 8 : PetscCall(BVRestoreColumn(W,j,&w));
221 : }
222 4 : if (n) {
223 4 : PetscCall(BVGetArrayRead(C1,&A));
224 4 : PetscCall(BVGetLeadingDimension(C1,&lda));
225 4 : PetscCall(PetscBLASIntCast(lda,&lda_));
226 4 : PetscCall(BVGetArrayRead(W,&B));
227 4 : PetscCall(BVGetLeadingDimension(W,&ldb));
228 4 : PetscCall(PetscBLASIntCast(ldb,&ldb_));
229 4 : PetscCallBLAS("BLASgemm",BLASgemm_("N","C",&n_,&N_,&l_,&alpha,(PetscScalar*)A,&lda_,(PetscScalar*)B,&ldb_,&beta,Rarray,&n_));
230 4 : PetscCall(BVRestoreArrayRead(C1,&A));
231 4 : PetscCall(BVRestoreArrayRead(W,&B));
232 : }
233 4 : beta = 1.0;
234 :
235 : /* R+=AX*X1' */
236 272 : for (j=0;j<k;j++) {
237 268 : PetscCall(BVGetColumn(X1,j,&v));
238 268 : PetscCall(BVGetColumn(W,j,&w));
239 268 : PetscCall(VecScatterBegin(vscat,v,w,INSERT_VALUES,SCATTER_FORWARD));
240 268 : PetscCall(VecScatterEnd(vscat,v,w,INSERT_VALUES,SCATTER_FORWARD));
241 268 : PetscCall(BVRestoreColumn(X1,j,&v));
242 268 : PetscCall(BVRestoreColumn(W,j,&w));
243 : }
244 4 : if (n) {
245 4 : PetscCall(BVGetArrayRead(AX,&A));
246 4 : PetscCall(BVGetLeadingDimension(AX,&lda));
247 4 : PetscCall(PetscBLASIntCast(lda,&lda_));
248 4 : PetscCall(BVGetArrayRead(W,&B));
249 4 : PetscCallBLAS("BLASgemm",BLASgemm_("N","C",&n_,&N_,&k_,&alpha,(PetscScalar*)A,&lda_,(PetscScalar*)B,&ldb_,&beta,Rarray,&n_));
250 4 : PetscCall(BVRestoreArrayRead(AX,&A));
251 4 : PetscCall(BVRestoreArrayRead(W,&B));
252 : }
253 :
254 : /* R+=X1*AX' */
255 272 : for (j=0;j<k;j++) {
256 268 : PetscCall(BVGetColumn(AX,j,&v));
257 268 : PetscCall(BVGetColumn(W,j,&w));
258 268 : PetscCall(VecScatterBegin(vscat,v,w,INSERT_VALUES,SCATTER_FORWARD));
259 268 : PetscCall(VecScatterEnd(vscat,v,w,INSERT_VALUES,SCATTER_FORWARD));
260 268 : PetscCall(BVRestoreColumn(AX,j,&v));
261 268 : PetscCall(BVRestoreColumn(W,j,&w));
262 : }
263 4 : if (n) {
264 4 : PetscCall(BVGetArrayRead(X1,&A));
265 4 : PetscCall(BVGetLeadingDimension(AX,&lda));
266 4 : PetscCall(PetscBLASIntCast(lda,&lda_));
267 4 : PetscCall(BVGetArrayRead(W,&B));
268 4 : PetscCallBLAS("BLASgemm",BLASgemm_("N","C",&n_,&N_,&k_,&alpha,(PetscScalar*)A,&lda_,(PetscScalar*)B,&ldb_,&beta,Rarray,&n_));
269 4 : PetscCall(BVRestoreArrayRead(X1,&A));
270 4 : PetscCall(BVRestoreArrayRead(W,&B));
271 : }
272 4 : PetscCall(MatDenseRestoreArrayWrite(R,&Rarray));
273 :
274 : /* compute ||R||_F */
275 4 : PetscCall(MatNorm(R,NORM_FROBENIUS,norm));
276 :
277 4 : PetscCall(BVDestroy(&W));
278 4 : PetscCall(VecScatterDestroy(&vscat));
279 4 : PetscCall(BVDestroy(&AX));
280 4 : PetscCall(MatDestroy(&R));
281 4 : PetscCall(BVDestroy(&C1));
282 4 : PetscCall(BVDestroy(&X1));
283 4 : PetscFunctionReturn(PETSC_SUCCESS);
284 : }
285 :
286 : /*@
287 : LMEComputeError - Computes the error (based on the residual norm) associated
288 : with the last equation solved.
289 :
290 : Collective
291 :
292 : Input Parameter:
293 : . lme - the linear matrix equation solver context
294 :
295 : Output Parameter:
296 : . error - the error
297 :
298 : Notes:
299 : This function is not scalable (in terms of memory or parallel communication),
300 : so it should not be called except in the case of small problem size. For
301 : large equations, use LMEGetErrorEstimate().
302 :
303 : Level: advanced
304 :
305 : .seealso: LMESolve(), LMEGetErrorEstimate()
306 : @*/
307 4 : PetscErrorCode LMEComputeError(LME lme,PetscReal *error)
308 : {
309 4 : PetscFunctionBegin;
310 4 : PetscValidHeaderSpecific(lme,LME_CLASSID,1);
311 4 : PetscAssertPointer(error,2);
312 :
313 4 : PetscCall(PetscLogEventBegin(LME_ComputeError,lme,0,0,0));
314 : /* compute residual norm */
315 4 : switch (lme->problem_type) {
316 4 : case LME_LYAPUNOV:
317 4 : PetscCall(LMEComputeResidualNorm_Lyapunov(lme,error));
318 4 : break;
319 0 : default:
320 0 : SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_SUP,"Not implemented for equation type %s",LMEProblemTypes[lme->problem_type]);
321 : }
322 :
323 : /* compute error */
324 : /* currently we only support absolute error, so just return the norm */
325 4 : PetscCall(PetscLogEventEnd(LME_ComputeError,lme,0,0,0));
326 4 : PetscFunctionReturn(PETSC_SUCCESS);
327 : }
|
