Actual source code: lmesolve.c
slepc-3.18.2 2023-01-26
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 the solution process
12: */
14: #include <slepc/private/lmeimpl.h>
15: #include <slepcblaslapack.h>
17: /*@
18: LMESolve - Solves the linear matrix equation.
20: Collective on lme
22: Input Parameter:
23: . lme - linear matrix equation solver context obtained from LMECreate()
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
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().
37: Level: beginner
39: .seealso: LMECreate(), LMESetUp(), LMEDestroy(), LMESetTolerances(), LMESetCoefficients(), LMESetRHS(), LMESetSolution()
40: @*/
41: PetscErrorCode LMESolve(LME lme)
42: {
45: /* call setup */
46: LMESetUp(lme);
47: lme->its = 0;
48: lme->errest = 0.0;
50: LMEViewFromOptions(lme,NULL,"-lme_view_pre");
52: /* call solver */
54: PetscLogEventBegin(LME_Solve,lme,0,0,0);
55: PetscUseTypeMethod(lme,solve[lme->problem_type]);
56: PetscLogEventEnd(LME_Solve,lme,0,0,0);
62: /* various viewers */
63: LMEViewFromOptions(lme,NULL,"-lme_view");
64: LMEConvergedReasonViewFromOptions(lme);
65: MatViewFromOptions(lme->A,(PetscObject)lme,"-lme_view_mat");
66: MatViewFromOptions(lme->C,(PetscObject)lme,"-lme_view_rhs");
67: MatViewFromOptions(lme->X,(PetscObject)lme,"-lme_view_solution");
68: return 0;
69: }
71: /*@
72: LMEGetIterationNumber - Gets the current iteration number. If the
73: call to LMESolve() is complete, then it returns the number of iterations
74: carried out by the solution method.
76: Not Collective
78: Input Parameter:
79: . lme - the linear matrix equation solver context
81: Output Parameter:
82: . its - number of iterations
84: Note:
85: During the i-th iteration this call returns i-1. If LMESolve() is
86: complete, then parameter "its" contains either the iteration number at
87: which convergence was successfully reached, or failure was detected.
88: Call LMEGetConvergedReason() to determine if the solver converged or
89: failed and why.
91: Level: intermediate
93: .seealso: LMEGetConvergedReason(), LMESetTolerances()
94: @*/
95: PetscErrorCode LMEGetIterationNumber(LME lme,PetscInt *its)
96: {
99: *its = lme->its;
100: return 0;
101: }
103: /*@
104: LMEGetConvergedReason - Gets the reason why the LMESolve() iteration was
105: stopped.
107: Not Collective
109: Input Parameter:
110: . lme - the linear matrix equation solver context
112: Output Parameter:
113: . reason - negative value indicates diverged, positive value converged
115: Notes:
117: Possible values for reason are
118: + LME_CONVERGED_TOL - converged up to tolerance
119: . LME_DIVERGED_ITS - required more than max_it iterations to reach convergence
120: - LME_DIVERGED_BREAKDOWN - generic breakdown in method
122: Can only be called after the call to LMESolve() is complete.
124: Level: intermediate
126: .seealso: LMESetTolerances(), LMESolve(), LMEConvergedReason, LMESetErrorIfNotConverged()
127: @*/
128: PetscErrorCode LMEGetConvergedReason(LME lme,LMEConvergedReason *reason)
129: {
132: *reason = lme->reason;
133: return 0;
134: }
136: /*@
137: LMEGetErrorEstimate - Returns the error estimate obtained during solve.
139: Not Collective
141: Input Parameter:
142: . lme - linear matrix equation solver context
144: Output Parameter:
145: . errest - the error estimate
147: Notes:
148: This is the error estimated internally by the solver. The actual
149: error bound can be computed with LMEComputeError(). Note that some
150: solvers may not be able to provide an error estimate.
152: Level: advanced
154: .seealso: LMEComputeError()
155: @*/
156: PetscErrorCode LMEGetErrorEstimate(LME lme,PetscReal *errest)
157: {
160: *errest = lme->errest;
161: return 0;
162: }
164: /*
165: LMEComputeResidualNorm_Lyapunov - Computes the Frobenius norm of the residual matrix
166: associated with the Lyapunov equation.
167: */
168: PetscErrorCode LMEComputeResidualNorm_Lyapunov(LME lme,PetscReal *norm)
169: {
170: PetscInt j,n,N,k,l;
171: PetscBLASInt n_,N_,k_,l_;
172: PetscScalar *Rarray,alpha=1.0,beta=0.0;
173: const PetscScalar *A,*B;
174: BV W,AX,X1,C1;
175: Mat R,X1m,C1m;
176: Vec v,w;
177: VecScatter vscat;
179: MatLRCGetMats(lme->C,NULL,&C1m,NULL,NULL);
180: MatLRCGetMats(lme->X,NULL,&X1m,NULL,NULL);
181: BVCreateFromMat(C1m,&C1);
182: BVSetFromOptions(C1);
183: BVCreateFromMat(X1m,&X1);
184: BVSetFromOptions(X1);
185: BVGetSizes(X1,&n,&N,&k);
186: BVGetSizes(C1,NULL,NULL,&l);
187: PetscBLASIntCast(n,&n_);
188: PetscBLASIntCast(N,&N_);
189: PetscBLASIntCast(k,&k_);
190: PetscBLASIntCast(l,&l_);
192: /* create W to store a redundant copy of a BV in each process */
193: BVCreate(PETSC_COMM_SELF,&W);
194: BVSetSizes(W,N,N,k);
195: BVSetFromOptions(W);
196: BVGetColumn(X1,0,&v);
197: VecScatterCreateToAll(v,&vscat,NULL);
198: BVRestoreColumn(X1,0,&v);
200: /* create AX to hold the product A*X1 */
201: BVDuplicate(X1,&AX);
202: BVMatMult(X1,lme->A,AX);
204: /* create dense matrix to hold the residual R=C1*C1'+AX*X1'+X1*AX' */
205: MatCreateDense(PetscObjectComm((PetscObject)lme),n,n,N,N,NULL,&R);
207: /* R=C1*C1' */
208: MatDenseGetArrayWrite(R,&Rarray);
209: for (j=0;j<l;j++) {
210: BVGetColumn(C1,j,&v);
211: BVGetColumn(W,j,&w);
212: VecScatterBegin(vscat,v,w,INSERT_VALUES,SCATTER_FORWARD);
213: VecScatterEnd(vscat,v,w,INSERT_VALUES,SCATTER_FORWARD);
214: BVRestoreColumn(C1,j,&v);
215: BVRestoreColumn(W,j,&w);
216: }
217: if (n) {
218: BVGetArrayRead(C1,&A);
219: BVGetArrayRead(W,&B);
220: PetscCallBLAS("BLASgemm",BLASgemm_("N","C",&n_,&N_,&l_,&alpha,(PetscScalar*)A,&n_,(PetscScalar*)B,&N_,&beta,Rarray,&n_));
221: BVRestoreArrayRead(C1,&A);
222: BVRestoreArrayRead(W,&B);
223: }
224: beta = 1.0;
226: /* R+=AX*X1' */
227: for (j=0;j<k;j++) {
228: BVGetColumn(X1,j,&v);
229: BVGetColumn(W,j,&w);
230: VecScatterBegin(vscat,v,w,INSERT_VALUES,SCATTER_FORWARD);
231: VecScatterEnd(vscat,v,w,INSERT_VALUES,SCATTER_FORWARD);
232: BVRestoreColumn(X1,j,&v);
233: BVRestoreColumn(W,j,&w);
234: }
235: if (n) {
236: BVGetArrayRead(AX,&A);
237: BVGetArrayRead(W,&B);
238: PetscCallBLAS("BLASgemm",BLASgemm_("N","C",&n_,&N_,&k_,&alpha,(PetscScalar*)A,&n_,(PetscScalar*)B,&N_,&beta,Rarray,&n_));
239: BVRestoreArrayRead(AX,&A);
240: BVRestoreArrayRead(W,&B);
241: }
243: /* R+=X1*AX' */
244: for (j=0;j<k;j++) {
245: BVGetColumn(AX,j,&v);
246: BVGetColumn(W,j,&w);
247: VecScatterBegin(vscat,v,w,INSERT_VALUES,SCATTER_FORWARD);
248: VecScatterEnd(vscat,v,w,INSERT_VALUES,SCATTER_FORWARD);
249: BVRestoreColumn(AX,j,&v);
250: BVRestoreColumn(W,j,&w);
251: }
252: if (n) {
253: BVGetArrayRead(X1,&A);
254: BVGetArrayRead(W,&B);
255: PetscCallBLAS("BLASgemm",BLASgemm_("N","C",&n_,&N_,&k_,&alpha,(PetscScalar*)A,&n_,(PetscScalar*)B,&N_,&beta,Rarray,&n_));
256: BVRestoreArrayRead(X1,&A);
257: BVRestoreArrayRead(W,&B);
258: }
259: MatDenseRestoreArrayWrite(R,&Rarray);
261: /* compute ||R||_F */
262: MatNorm(R,NORM_FROBENIUS,norm);
264: BVDestroy(&W);
265: VecScatterDestroy(&vscat);
266: BVDestroy(&AX);
267: MatDestroy(&R);
268: BVDestroy(&C1);
269: BVDestroy(&X1);
270: return 0;
271: }
273: /*@
274: LMEComputeError - Computes the error (based on the residual norm) associated
275: with the last equation solved.
277: Collective on lme
279: Input Parameter:
280: . lme - the linear matrix equation solver context
282: Output Parameter:
283: . error - the error
285: Notes:
286: This function is not scalable (in terms of memory or parallel communication),
287: so it should not be called except in the case of small problem size. For
288: large equations, use LMEGetErrorEstimate().
290: Level: advanced
292: .seealso: LMESolve(), LMEGetErrorEstimate()
293: @*/
294: PetscErrorCode LMEComputeError(LME lme,PetscReal *error)
295: {
299: PetscLogEventBegin(LME_ComputeError,lme,0,0,0);
300: /* compute residual norm */
301: switch (lme->problem_type) {
302: case LME_LYAPUNOV:
303: LMEComputeResidualNorm_Lyapunov(lme,error);
304: break;
305: default:
306: SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_SUP,"Not implemented for equation type %s",LMEProblemTypes[lme->problem_type]);
307: }
309: /* compute error */
310: /* currently we only support absolute error, so just return the norm */
311: PetscLogEventEnd(LME_ComputeError,lme,0,0,0);
312: return 0;
313: }