Actual source code: lmesolve.c

slepc-3.18.3 2023-03-24
Report Typos and Errors
  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: }