Actual source code: lmesolve.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  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: }