Actual source code: lmesolve.c
slepc-3.21.0 2024-03-30
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
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: {
43: PetscFunctionBegin;
46: /* call setup */
47: PetscCall(LMESetUp(lme));
48: lme->its = 0;
49: lme->errest = 0.0;
51: PetscCall(LMEViewFromOptions(lme,NULL,"-lme_view_pre"));
53: /* call solver */
54: 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: PetscCall(PetscLogEventBegin(LME_Solve,lme,0,0,0));
56: PetscUseTypeMethod(lme,solve[lme->problem_type]);
57: PetscCall(PetscLogEventEnd(LME_Solve,lme,0,0,0));
59: PetscCheck(lme->reason,PetscObjectComm((PetscObject)lme),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
61: PetscCheck(!lme->errorifnotconverged || lme->reason>=0,PetscObjectComm((PetscObject)lme),PETSC_ERR_NOT_CONVERGED,"LMESolve has not converged");
63: /* various viewers */
64: PetscCall(LMEViewFromOptions(lme,NULL,"-lme_view"));
65: PetscCall(LMEConvergedReasonViewFromOptions(lme));
66: PetscCall(MatViewFromOptions(lme->A,(PetscObject)lme,"-lme_view_mat"));
67: PetscCall(MatViewFromOptions(lme->C,(PetscObject)lme,"-lme_view_rhs"));
68: PetscCall(MatViewFromOptions(lme->X,(PetscObject)lme,"-lme_view_solution"));
69: PetscFunctionReturn(PETSC_SUCCESS);
70: }
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.
77: Not Collective
79: Input Parameter:
80: . lme - the linear matrix equation solver context
82: Output Parameter:
83: . its - number of iterations
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.
92: Level: intermediate
94: .seealso: LMEGetConvergedReason(), LMESetTolerances()
95: @*/
96: PetscErrorCode LMEGetIterationNumber(LME lme,PetscInt *its)
97: {
98: PetscFunctionBegin;
100: PetscAssertPointer(its,2);
101: *its = lme->its;
102: PetscFunctionReturn(PETSC_SUCCESS);
103: }
105: /*@
106: LMEGetConvergedReason - Gets the reason why the LMESolve() iteration was
107: stopped.
109: Not Collective
111: Input Parameter:
112: . lme - the linear matrix equation solver context
114: Output Parameter:
115: . reason - negative value indicates diverged, positive value converged
117: Notes:
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
124: Can only be called after the call to LMESolve() is complete.
126: Level: intermediate
128: .seealso: LMESetTolerances(), LMESolve(), LMEConvergedReason, LMESetErrorIfNotConverged()
129: @*/
130: PetscErrorCode LMEGetConvergedReason(LME lme,LMEConvergedReason *reason)
131: {
132: PetscFunctionBegin;
134: PetscAssertPointer(reason,2);
135: *reason = lme->reason;
136: PetscFunctionReturn(PETSC_SUCCESS);
137: }
139: /*@
140: LMEGetErrorEstimate - Returns the error estimate obtained during solve.
142: Not Collective
144: Input Parameter:
145: . lme - linear matrix equation solver context
147: Output Parameter:
148: . errest - the error estimate
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.
155: Level: advanced
157: .seealso: LMEComputeError()
158: @*/
159: PetscErrorCode LMEGetErrorEstimate(LME lme,PetscReal *errest)
160: {
161: PetscFunctionBegin;
163: PetscAssertPointer(errest,2);
164: *errest = lme->errest;
165: PetscFunctionReturn(PETSC_SUCCESS);
166: }
168: /*
169: LMEComputeResidualNorm_Lyapunov - Computes the Frobenius norm of the residual matrix
170: associated with the Lyapunov equation.
171: */
172: static PetscErrorCode LMEComputeResidualNorm_Lyapunov(LME lme,PetscReal *norm)
173: {
174: PetscInt j,n,N,k,l,lda,ldb;
175: PetscBLASInt n_,N_,k_,l_,lda_,ldb_;
176: PetscScalar *Rarray,alpha=1.0,beta=0.0;
177: const PetscScalar *A,*B;
178: BV W,AX,X1,C1;
179: Mat R,X1m,C1m;
180: Vec v,w;
181: VecScatter vscat;
183: PetscFunctionBegin;
184: PetscCall(MatLRCGetMats(lme->C,NULL,&C1m,NULL,NULL));
185: PetscCall(MatLRCGetMats(lme->X,NULL,&X1m,NULL,NULL));
186: PetscCall(BVCreateFromMat(C1m,&C1));
187: PetscCall(BVSetFromOptions(C1));
188: PetscCall(BVCreateFromMat(X1m,&X1));
189: PetscCall(BVSetFromOptions(X1));
190: PetscCall(BVGetSizes(X1,&n,&N,&k));
191: PetscCall(BVGetSizes(C1,NULL,NULL,&l));
192: PetscCall(PetscBLASIntCast(n,&n_));
193: PetscCall(PetscBLASIntCast(N,&N_));
194: PetscCall(PetscBLASIntCast(k,&k_));
195: PetscCall(PetscBLASIntCast(l,&l_));
197: /* create W to store a redundant copy of a BV in each process */
198: PetscCall(BVCreate(PETSC_COMM_SELF,&W));
199: PetscCall(BVSetSizes(W,N,N,k));
200: PetscCall(BVSetFromOptions(W));
201: PetscCall(BVGetColumn(X1,0,&v));
202: PetscCall(VecScatterCreateToAll(v,&vscat,NULL));
203: PetscCall(BVRestoreColumn(X1,0,&v));
205: /* create AX to hold the product A*X1 */
206: PetscCall(BVDuplicate(X1,&AX));
207: PetscCall(BVMatMult(X1,lme->A,AX));
209: /* create dense matrix to hold the residual R=C1*C1'+AX*X1'+X1*AX' */
210: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)lme),n,n,N,N,NULL,&R));
212: /* R=C1*C1' */
213: PetscCall(MatDenseGetArrayWrite(R,&Rarray));
214: for (j=0;j<l;j++) {
215: PetscCall(BVGetColumn(C1,j,&v));
216: PetscCall(BVGetColumn(W,j,&w));
217: PetscCall(VecScatterBegin(vscat,v,w,INSERT_VALUES,SCATTER_FORWARD));
218: PetscCall(VecScatterEnd(vscat,v,w,INSERT_VALUES,SCATTER_FORWARD));
219: PetscCall(BVRestoreColumn(C1,j,&v));
220: PetscCall(BVRestoreColumn(W,j,&w));
221: }
222: if (n) {
223: PetscCall(BVGetArrayRead(C1,&A));
224: PetscCall(BVGetLeadingDimension(C1,&lda));
225: PetscCall(PetscBLASIntCast(lda,&lda_));
226: PetscCall(BVGetArrayRead(W,&B));
227: PetscCall(BVGetLeadingDimension(W,&ldb));
228: PetscCall(PetscBLASIntCast(ldb,&ldb_));
229: PetscCallBLAS("BLASgemm",BLASgemm_("N","C",&n_,&N_,&l_,&alpha,(PetscScalar*)A,&lda_,(PetscScalar*)B,&ldb_,&beta,Rarray,&n_));
230: PetscCall(BVRestoreArrayRead(C1,&A));
231: PetscCall(BVRestoreArrayRead(W,&B));
232: }
233: beta = 1.0;
235: /* R+=AX*X1' */
236: for (j=0;j<k;j++) {
237: PetscCall(BVGetColumn(X1,j,&v));
238: PetscCall(BVGetColumn(W,j,&w));
239: PetscCall(VecScatterBegin(vscat,v,w,INSERT_VALUES,SCATTER_FORWARD));
240: PetscCall(VecScatterEnd(vscat,v,w,INSERT_VALUES,SCATTER_FORWARD));
241: PetscCall(BVRestoreColumn(X1,j,&v));
242: PetscCall(BVRestoreColumn(W,j,&w));
243: }
244: if (n) {
245: PetscCall(BVGetArrayRead(AX,&A));
246: PetscCall(BVGetLeadingDimension(AX,&lda));
247: PetscCall(PetscBLASIntCast(lda,&lda_));
248: PetscCall(BVGetArrayRead(W,&B));
249: PetscCallBLAS("BLASgemm",BLASgemm_("N","C",&n_,&N_,&k_,&alpha,(PetscScalar*)A,&lda_,(PetscScalar*)B,&ldb_,&beta,Rarray,&n_));
250: PetscCall(BVRestoreArrayRead(AX,&A));
251: PetscCall(BVRestoreArrayRead(W,&B));
252: }
254: /* R+=X1*AX' */
255: for (j=0;j<k;j++) {
256: PetscCall(BVGetColumn(AX,j,&v));
257: PetscCall(BVGetColumn(W,j,&w));
258: PetscCall(VecScatterBegin(vscat,v,w,INSERT_VALUES,SCATTER_FORWARD));
259: PetscCall(VecScatterEnd(vscat,v,w,INSERT_VALUES,SCATTER_FORWARD));
260: PetscCall(BVRestoreColumn(AX,j,&v));
261: PetscCall(BVRestoreColumn(W,j,&w));
262: }
263: if (n) {
264: PetscCall(BVGetArrayRead(X1,&A));
265: PetscCall(BVGetLeadingDimension(AX,&lda));
266: PetscCall(PetscBLASIntCast(lda,&lda_));
267: PetscCall(BVGetArrayRead(W,&B));
268: PetscCallBLAS("BLASgemm",BLASgemm_("N","C",&n_,&N_,&k_,&alpha,(PetscScalar*)A,&lda_,(PetscScalar*)B,&ldb_,&beta,Rarray,&n_));
269: PetscCall(BVRestoreArrayRead(X1,&A));
270: PetscCall(BVRestoreArrayRead(W,&B));
271: }
272: PetscCall(MatDenseRestoreArrayWrite(R,&Rarray));
274: /* compute ||R||_F */
275: PetscCall(MatNorm(R,NORM_FROBENIUS,norm));
277: PetscCall(BVDestroy(&W));
278: PetscCall(VecScatterDestroy(&vscat));
279: PetscCall(BVDestroy(&AX));
280: PetscCall(MatDestroy(&R));
281: PetscCall(BVDestroy(&C1));
282: PetscCall(BVDestroy(&X1));
283: PetscFunctionReturn(PETSC_SUCCESS);
284: }
286: /*@
287: LMEComputeError - Computes the error (based on the residual norm) associated
288: with the last equation solved.
290: Collective
292: Input Parameter:
293: . lme - the linear matrix equation solver context
295: Output Parameter:
296: . error - the error
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().
303: Level: advanced
305: .seealso: LMESolve(), LMEGetErrorEstimate()
306: @*/
307: PetscErrorCode LMEComputeError(LME lme,PetscReal *error)
308: {
309: PetscFunctionBegin;
311: PetscAssertPointer(error,2);
313: PetscCall(PetscLogEventBegin(LME_ComputeError,lme,0,0,0));
314: /* compute residual norm */
315: switch (lme->problem_type) {
316: case LME_LYAPUNOV:
317: PetscCall(LMEComputeResidualNorm_Lyapunov(lme,error));
318: break;
319: default:
320: SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_SUP,"Not implemented for equation type %s",LMEProblemTypes[lme->problem_type]);
321: }
323: /* compute error */
324: /* currently we only support absolute error, so just return the norm */
325: PetscCall(PetscLogEventEnd(LME_ComputeError,lme,0,0,0));
326: PetscFunctionReturn(PETSC_SUCCESS);
327: }