Actual source code: lmesolve.c
slepc-3.11.2 2019-07-30
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 the solution process
12: */
14: #include <slepc/private/lmeimpl.h> /*I "slepclme.h" I*/
15: #include <slepcblaslapack.h>
17: /*@
18: LMESolve - Solves the linear matrix equation.
20: Collective on LME
22: Input Parameters:
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: {
48: /* call setup */
49: LMESetUp(lme);
50: lme->its = 0;
51: lme->errest = 0.0;
53: LMEViewFromOptions(lme,NULL,"-lme_view_pre");
55: /* call solver */
56: if (!lme->ops->solve[lme->problem_type]) SETERRQ1(PetscObjectComm((PetscObject)lme),PETSC_ERR_SUP,"The specified solver does not support equation type %s",LMEProblemTypes[lme->problem_type]);
57: PetscLogEventBegin(LME_Solve,lme,0,0,0);
58: (*lme->ops->solve[lme->problem_type])(lme);
59: PetscLogEventEnd(LME_Solve,lme,0,0,0);
61: if (!lme->reason) SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
63: if (lme->errorifnotconverged && lme->reason < 0) SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_NOT_CONVERGED,"LMESolve has not converged");
65: /* various viewers */
66: LMEViewFromOptions(lme,NULL,"-lme_view");
67: LMEReasonViewFromOptions(lme);
68: MatViewFromOptions(lme->A,(PetscObject)lme,"-lme_view_mat");
69: MatViewFromOptions(lme->C,(PetscObject)lme,"-lme_view_rhs");
70: MatViewFromOptions(lme->X,(PetscObject)lme,"-lme_view_solution");
71: return(0);
72: }
74: /*@
75: LMEGetIterationNumber - Gets the current iteration number. If the
76: call to LMESolve() is complete, then it returns the number of iterations
77: carried out by the solution method.
79: Not Collective
81: Input Parameter:
82: . lme - the linear matrix equation solver context
84: Output Parameter:
85: . its - number of iterations
87: Note:
88: During the i-th iteration this call returns i-1. If LMESolve() is
89: complete, then parameter "its" contains either the iteration number at
90: which convergence was successfully reached, or failure was detected.
91: Call LMEGetConvergedReason() to determine if the solver converged or
92: failed and why.
94: Level: intermediate
96: .seealso: LMEGetConvergedReason(), LMESetTolerances()
97: @*/
98: PetscErrorCode LMEGetIterationNumber(LME lme,PetscInt *its)
99: {
103: *its = lme->its;
104: return(0);
105: }
107: /*@
108: LMEGetConvergedReason - Gets the reason why the LMESolve() iteration was
109: stopped.
111: Not Collective
113: Input Parameter:
114: . lme - the linear matrix equation solver context
116: Output Parameter:
117: . reason - negative value indicates diverged, positive value converged
119: Notes:
121: Possible values for reason are
122: + LME_CONVERGED_TOL - converged up to tolerance
123: . LME_DIVERGED_ITS - required more than max_it iterations to reach convergence
124: - LME_DIVERGED_BREAKDOWN - generic breakdown in method
126: Can only be called after the call to LMESolve() is complete.
128: Level: intermediate
130: .seealso: LMESetTolerances(), LMESolve(), LMEConvergedReason, LMESetErrorIfNotConverged()
131: @*/
132: PetscErrorCode LMEGetConvergedReason(LME lme,LMEConvergedReason *reason)
133: {
137: *reason = lme->reason;
138: return(0);
139: }
141: /*@
142: LMEGetErrorEstimate - Returns the error estimate obtained during solve.
144: Not Collective
146: Input Parameter:
147: . lme - linear matrix equation solver context
149: Output Parameter:
150: . errest - the error estimate
152: Notes:
153: This is the error estimated internally by the solver. The actual
154: error bound can be computed with LMEComputeError(). Note that some
155: solvers may not be able to provide an error estimate.
157: Level: advanced
159: .seealso: LMEComputeError()
160: @*/
161: PetscErrorCode LMEGetErrorEstimate(LME lme,PetscReal *errest)
162: {
166: *errest = lme->errest;
167: return(0);
168: }
170: /*
171: LMEComputeResidualNorm_Lyapunov - Computes the Frobenius norm of the residual matrix
172: associated with the Lyapunov equation.
173: */
174: PetscErrorCode LMEComputeResidualNorm_Lyapunov(LME lme,PetscReal *norm)
175: {
176: PetscErrorCode ierr;
177: PetscInt j,n,N,k,l;
178: PetscBLASInt n_,N_,k_,l_;
179: PetscScalar *Rarray,alpha=1.0,beta=0.0;
180: const PetscScalar *A,*B;
181: BV W,AX,X1,C1;
182: Mat R,X1m,C1m;
183: Vec v,w;
184: VecScatter vscat;
187: MatLRCGetMats(lme->C,NULL,&C1m,NULL,NULL);
188: MatLRCGetMats(lme->X,NULL,&X1m,NULL,NULL);
189: BVCreateFromMat(C1m,&C1);
190: BVSetFromOptions(C1);
191: BVCreateFromMat(X1m,&X1);
192: BVSetFromOptions(X1);
193: BVGetSizes(X1,&n,&N,&k);
194: BVGetSizes(C1,NULL,NULL,&l);
195: PetscBLASIntCast(n,&n_);
196: PetscBLASIntCast(N,&N_);
197: PetscBLASIntCast(k,&k_);
198: PetscBLASIntCast(l,&l_);
200: /* create W to store a redundant copy of a BV in each process */
201: BVCreate(PETSC_COMM_SELF,&W);
202: BVSetSizes(W,N,N,k);
203: BVSetFromOptions(W);
204: BVGetColumn(X1,0,&v);
205: VecScatterCreateToAll(v,&vscat,NULL);
206: BVRestoreColumn(X1,0,&v);
208: /* create AX to hold the product A*X1 */
209: BVDuplicate(X1,&AX);
210: BVMatMult(X1,lme->A,AX);
212: /* create dense matrix to hold the residual R=C1*C1'+AX*X1'+X1*AX' */
213: MatCreateDense(PetscObjectComm((PetscObject)lme),n,n,N,N,NULL,&R);
215: /* R=C1*C1' */
216: MatDenseGetArray(R,&Rarray);
217: for (j=0;j<l;j++) {
218: BVGetColumn(C1,j,&v);
219: BVGetColumn(W,j,&w);
220: VecScatterBegin(vscat,v,w,INSERT_VALUES,SCATTER_FORWARD);
221: VecScatterEnd(vscat,v,w,INSERT_VALUES,SCATTER_FORWARD);
222: BVRestoreColumn(C1,j,&v);
223: BVRestoreColumn(W,j,&w);
224: }
225: if (n) {
226: BVGetArrayRead(C1,&A);
227: BVGetArrayRead(W,&B);
228: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","C",&n_,&N_,&l_,&alpha,(PetscScalar*)A,&n_,(PetscScalar*)B,&N_,&beta,Rarray,&n_));
229: BVRestoreArrayRead(C1,&A);
230: BVRestoreArrayRead(W,&B);
231: }
232: beta = 1.0;
234: /* R+=AX*X1' */
235: for (j=0;j<k;j++) {
236: BVGetColumn(X1,j,&v);
237: BVGetColumn(W,j,&w);
238: VecScatterBegin(vscat,v,w,INSERT_VALUES,SCATTER_FORWARD);
239: VecScatterEnd(vscat,v,w,INSERT_VALUES,SCATTER_FORWARD);
240: BVRestoreColumn(X1,j,&v);
241: BVRestoreColumn(W,j,&w);
242: }
243: if (n) {
244: BVGetArrayRead(AX,&A);
245: BVGetArrayRead(W,&B);
246: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","C",&n_,&N_,&k_,&alpha,(PetscScalar*)A,&n_,(PetscScalar*)B,&N_,&beta,Rarray,&n_));
247: BVRestoreArrayRead(AX,&A);
248: BVRestoreArrayRead(W,&B);
249: }
251: /* R+=X1*AX' */
252: for (j=0;j<k;j++) {
253: BVGetColumn(AX,j,&v);
254: BVGetColumn(W,j,&w);
255: VecScatterBegin(vscat,v,w,INSERT_VALUES,SCATTER_FORWARD);
256: VecScatterEnd(vscat,v,w,INSERT_VALUES,SCATTER_FORWARD);
257: BVRestoreColumn(AX,j,&v);
258: BVRestoreColumn(W,j,&w);
259: }
260: if (n) {
261: BVGetArrayRead(X1,&A);
262: BVGetArrayRead(W,&B);
263: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","C",&n_,&N_,&k_,&alpha,(PetscScalar*)A,&n_,(PetscScalar*)B,&N_,&beta,Rarray,&n_));
264: BVRestoreArrayRead(X1,&A);
265: BVRestoreArrayRead(W,&B);
266: }
267: MatDenseRestoreArray(R,&Rarray);
269: /* compute ||R||_F */
270: MatAssemblyBegin(R,MAT_FINAL_ASSEMBLY);
271: MatAssemblyEnd(R,MAT_FINAL_ASSEMBLY);
272: MatNorm(R,NORM_FROBENIUS,norm);
274: BVDestroy(&W);
275: VecScatterDestroy(&vscat);
276: BVDestroy(&AX);
277: MatDestroy(&R);
278: BVDestroy(&C1);
279: BVDestroy(&X1);
280: return(0);
281: }
283: /*@
284: LMEComputeError - Computes the error (based on the residual norm) associated
285: with the last equation solved.
287: Collective on LME
289: Input Parameter:
290: . lme - the linear matrix equation solver context
292: Output Parameter:
293: . error - the error
295: Notes:
296: This function is not scalable (in terms of memory or parallel communication),
297: so it should not be called except in the case of small problem size. For
298: large equations, use LMEGetErrorEstimate().
300: Level: advanced
302: .seealso: LMESolve(), LMEGetErrorEstimate()
303: @*/
304: PetscErrorCode LMEComputeError(LME lme,PetscReal *error)
305: {
312: PetscLogEventBegin(LME_ComputeError,lme,0,0,0);
313: /* compute residual norm */
314: switch (lme->problem_type) {
315: case LME_LYAPUNOV:
316: LMEComputeResidualNorm_Lyapunov(lme,error);
317: break;
318: default:
319: SETERRQ1(PetscObjectComm((PetscObject)lme),PETSC_ERR_SUP,"Not implemented for equation type %s",LMEProblemTypes[lme->problem_type]);
320: }
322: /* compute error */
323: /* currently we only support absolute error, so just return the norm */
324: PetscLogEventEnd(LME_ComputeError,lme,0,0,0);
325: return(0);
326: }