Actual source code: svdbasic.c

slepc-3.16.1 2021-11-17
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2021, 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:    Basic SVD routines
 12: */

 14: #include <slepc/private/svdimpl.h>

 16: /* Logging support */
 17: PetscClassId      SVD_CLASSID = 0;
 18: PetscLogEvent     SVD_SetUp = 0,SVD_Solve = 0;

 20: /* List of registered SVD routines */
 21: PetscFunctionList SVDList = 0;
 22: PetscBool         SVDRegisterAllCalled = PETSC_FALSE;

 24: /* List of registered SVD monitors */
 25: PetscFunctionList SVDMonitorList              = NULL;
 26: PetscFunctionList SVDMonitorCreateList        = NULL;
 27: PetscFunctionList SVDMonitorDestroyList       = NULL;
 28: PetscBool         SVDMonitorRegisterAllCalled = PETSC_FALSE;

 30: /*@
 31:    SVDCreate - Creates the default SVD context.

 33:    Collective

 35:    Input Parameter:
 36: .  comm - MPI communicator

 38:    Output Parameter:
 39: .  svd - location to put the SVD context

 41:    Note:
 42:    The default SVD type is SVDCROSS

 44:    Level: beginner

 46: .seealso: SVDSetUp(), SVDSolve(), SVDDestroy(), SVD
 47: @*/
 48: PetscErrorCode SVDCreate(MPI_Comm comm,SVD *outsvd)
 49: {
 51:   SVD            svd;

 55:   *outsvd = 0;
 56:   SVDInitializePackage();
 57:   SlepcHeaderCreate(svd,SVD_CLASSID,"SVD","Singular Value Decomposition","SVD",comm,SVDDestroy,SVDView);

 59:   svd->OP               = NULL;
 60:   svd->OPb              = NULL;
 61:   svd->max_it           = PETSC_DEFAULT;
 62:   svd->nsv              = 1;
 63:   svd->ncv              = PETSC_DEFAULT;
 64:   svd->mpd              = PETSC_DEFAULT;
 65:   svd->nini             = 0;
 66:   svd->ninil            = 0;
 67:   svd->tol              = PETSC_DEFAULT;
 68:   svd->conv             = SVD_CONV_REL;
 69:   svd->stop             = SVD_STOP_BASIC;
 70:   svd->which            = SVD_LARGEST;
 71:   svd->problem_type     = (SVDProblemType)0;
 72:   svd->impltrans        = PETSC_FALSE;
 73:   svd->trackall         = PETSC_FALSE;

 75:   svd->converged        = SVDConvergedRelative;
 76:   svd->convergeduser    = NULL;
 77:   svd->convergeddestroy = NULL;
 78:   svd->stopping         = SVDStoppingBasic;
 79:   svd->stoppinguser     = NULL;
 80:   svd->stoppingdestroy  = NULL;
 81:   svd->convergedctx     = NULL;
 82:   svd->stoppingctx      = NULL;
 83:   svd->numbermonitors   = 0;

 85:   svd->ds               = NULL;
 86:   svd->U                = NULL;
 87:   svd->V                = NULL;
 88:   svd->A                = NULL;
 89:   svd->B                = NULL;
 90:   svd->AT               = NULL;
 91:   svd->BT               = NULL;
 92:   svd->IS               = NULL;
 93:   svd->ISL              = NULL;
 94:   svd->sigma            = NULL;
 95:   svd->errest           = NULL;
 96:   svd->perm             = NULL;
 97:   svd->nworkl           = 0;
 98:   svd->nworkr           = 0;
 99:   svd->workl            = NULL;
100:   svd->workr            = NULL;
101:   svd->data             = NULL;

103:   svd->state            = SVD_STATE_INITIAL;
104:   svd->nconv            = 0;
105:   svd->its              = 0;
106:   svd->leftbasis        = PETSC_FALSE;
107:   svd->swapped          = PETSC_FALSE;
108:   svd->expltrans        = PETSC_FALSE;
109:   svd->nrma             = 0.0;
110:   svd->nrmb             = 0.0;
111:   svd->isgeneralized    = PETSC_FALSE;
112:   svd->reason           = SVD_CONVERGED_ITERATING;

114:   PetscNewLog(svd,&svd->sc);
115:   *outsvd = svd;
116:   return(0);
117: }

119: /*@
120:    SVDReset - Resets the SVD context to the initial state (prior to setup)
121:    and destroys any allocated Vecs and Mats.

123:    Collective on svd

125:    Input Parameter:
126: .  svd - singular value solver context obtained from SVDCreate()

128:    Level: advanced

130: .seealso: SVDDestroy()
131: @*/
132: PetscErrorCode SVDReset(SVD svd)
133: {

138:   if (!svd) return(0);
139:   if (svd->ops->reset) { (svd->ops->reset)(svd); }
140:   MatDestroy(&svd->OP);
141:   MatDestroy(&svd->OPb);
142:   MatDestroy(&svd->A);
143:   MatDestroy(&svd->B);
144:   MatDestroy(&svd->AT);
145:   MatDestroy(&svd->BT);
146:   BVDestroy(&svd->U);
147:   BVDestroy(&svd->V);
148:   VecDestroyVecs(svd->nworkl,&svd->workl);
149:   svd->nworkl = 0;
150:   VecDestroyVecs(svd->nworkr,&svd->workr);
151:   svd->nworkr = 0;
152:   svd->state = SVD_STATE_INITIAL;
153:   return(0);
154: }

156: /*@C
157:    SVDDestroy - Destroys the SVD context.

159:    Collective on svd

161:    Input Parameter:
162: .  svd - singular value solver context obtained from SVDCreate()

164:    Level: beginner

166: .seealso: SVDCreate(), SVDSetUp(), SVDSolve()
167: @*/
168: PetscErrorCode SVDDestroy(SVD *svd)
169: {

173:   if (!*svd) return(0);
175:   if (--((PetscObject)(*svd))->refct > 0) { *svd = 0; return(0); }
176:   SVDReset(*svd);
177:   if ((*svd)->ops->destroy) { (*(*svd)->ops->destroy)(*svd); }
178:   if ((*svd)->sigma) {
179:     PetscFree3((*svd)->sigma,(*svd)->perm,(*svd)->errest);
180:   }
181:   DSDestroy(&(*svd)->ds);
182:   PetscFree((*svd)->sc);
183:   /* just in case the initial vectors have not been used */
184:   SlepcBasisDestroy_Private(&(*svd)->nini,&(*svd)->IS);
185:   SlepcBasisDestroy_Private(&(*svd)->ninil,&(*svd)->ISL);
186:   SVDMonitorCancel(*svd);
187:   PetscHeaderDestroy(svd);
188:   return(0);
189: }

191: /*@C
192:    SVDSetType - Selects the particular solver to be used in the SVD object.

194:    Logically Collective on svd

196:    Input Parameters:
197: +  svd      - the singular value solver context
198: -  type     - a known method

200:    Options Database Key:
201: .  -svd_type <method> - Sets the method; use -help for a list
202:     of available methods

204:    Notes:
205:    See "slepc/include/slepcsvd.h" for available methods. The default
206:    is SVDCROSS.

208:    Normally, it is best to use the SVDSetFromOptions() command and
209:    then set the SVD type from the options database rather than by using
210:    this routine.  Using the options database provides the user with
211:    maximum flexibility in evaluating the different available methods.
212:    The SVDSetType() routine is provided for those situations where it
213:    is necessary to set the iterative solver independently of the command
214:    line or options database.

216:    Level: intermediate

218: .seealso: SVDType
219: @*/
220: PetscErrorCode SVDSetType(SVD svd,SVDType type)
221: {
222:   PetscErrorCode ierr,(*r)(SVD);
223:   PetscBool      match;


229:   PetscObjectTypeCompare((PetscObject)svd,type,&match);
230:   if (match) return(0);

232:   PetscFunctionListFind(SVDList,type,&r);
233:   if (!r) SETERRQ1(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown SVD type given: %s",type);

235:   if (svd->ops->destroy) { (*svd->ops->destroy)(svd); }
236:   PetscMemzero(svd->ops,sizeof(struct _SVDOps));

238:   svd->state = SVD_STATE_INITIAL;
239:   PetscObjectChangeTypeName((PetscObject)svd,type);
240:   (*r)(svd);
241:   return(0);
242: }

244: /*@C
245:    SVDGetType - Gets the SVD type as a string from the SVD object.

247:    Not Collective

249:    Input Parameter:
250: .  svd - the singular value solver context

252:    Output Parameter:
253: .  name - name of SVD method

255:    Level: intermediate

257: .seealso: SVDSetType()
258: @*/
259: PetscErrorCode SVDGetType(SVD svd,SVDType *type)
260: {
264:   *type = ((PetscObject)svd)->type_name;
265:   return(0);
266: }

268: /*@C
269:    SVDRegister - Adds a method to the singular value solver package.

271:    Not Collective

273:    Input Parameters:
274: +  name - name of a new user-defined solver
275: -  function - routine to create the solver context

277:    Notes:
278:    SVDRegister() may be called multiple times to add several user-defined solvers.

280:    Sample usage:
281: .vb
282:     SVDRegister("my_solver",MySolverCreate);
283: .ve

285:    Then, your solver can be chosen with the procedural interface via
286: $     SVDSetType(svd,"my_solver")
287:    or at runtime via the option
288: $     -svd_type my_solver

290:    Level: advanced

292: .seealso: SVDRegisterAll()
293: @*/
294: PetscErrorCode SVDRegister(const char *name,PetscErrorCode (*function)(SVD))
295: {

299:   SVDInitializePackage();
300:   PetscFunctionListAdd(&SVDList,name,function);
301:   return(0);
302: }

304: /*@C
305:    SVDMonitorRegister - Adds SVD monitor routine.

307:    Not Collective

309:    Input Parameters:
310: +  name    - name of a new monitor routine
311: .  vtype   - a PetscViewerType for the output
312: .  format  - a PetscViewerFormat for the output
313: .  monitor - monitor routine
314: .  create  - creation routine, or NULL
315: -  destroy - destruction routine, or NULL

317:    Notes:
318:    SVDMonitorRegister() may be called multiple times to add several user-defined monitors.

320:    Sample usage:
321: .vb
322:    SVDMonitorRegister("my_monitor",PETSCVIEWERASCII,PETSC_VIEWER_ASCII_INFO_DETAIL,MyMonitor,NULL,NULL);
323: .ve

325:    Then, your monitor can be chosen with the procedural interface via
326: $      SVDMonitorSetFromOptions(svd,"-svd_monitor_my_monitor","my_monitor",NULL)
327:    or at runtime via the option
328: $      -svd_monitor_my_monitor

330:    Level: advanced

332: .seealso: SVDMonitorRegisterAll()
333: @*/
334: PetscErrorCode SVDMonitorRegister(const char name[],PetscViewerType vtype,PetscViewerFormat format,PetscErrorCode (*monitor)(SVD,PetscInt,PetscInt,PetscReal*,PetscReal*,PetscInt,PetscViewerAndFormat*),PetscErrorCode (*create)(PetscViewer,PetscViewerFormat,void*,PetscViewerAndFormat**),PetscErrorCode (*destroy)(PetscViewerAndFormat**))
335: {
336:   char           key[PETSC_MAX_PATH_LEN];

340:   SVDInitializePackage();
341:   SlepcMonitorMakeKey_Internal(name,vtype,format,key);
342:   PetscFunctionListAdd(&SVDMonitorList,key,monitor);
343:   if (create)  { PetscFunctionListAdd(&SVDMonitorCreateList,key,create); }
344:   if (destroy) { PetscFunctionListAdd(&SVDMonitorDestroyList,key,destroy); }
345:   return(0);
346: }

348: /*@
349:    SVDSetBV - Associates basis vectors objects to the singular value solver.

351:    Collective on svd

353:    Input Parameters:
354: +  svd - singular value solver context obtained from SVDCreate()
355: .  V   - the basis vectors object for right singular vectors
356: -  U   - the basis vectors object for left singular vectors

358:    Note:
359:    Use SVDGetBV() to retrieve the basis vectors contexts (for example,
360:    to free them at the end of the computations).

362:    Level: advanced

364: .seealso: SVDGetBV()
365: @*/
366: PetscErrorCode SVDSetBV(SVD svd,BV V,BV U)
367: {

372:   if (V) {
375:     PetscObjectReference((PetscObject)V);
376:     BVDestroy(&svd->V);
377:     svd->V = V;
378:     PetscLogObjectParent((PetscObject)svd,(PetscObject)svd->V);
379:   }
380:   if (U) {
383:     PetscObjectReference((PetscObject)U);
384:     BVDestroy(&svd->U);
385:     svd->U = U;
386:     PetscLogObjectParent((PetscObject)svd,(PetscObject)svd->U);
387:   }
388:   return(0);
389: }

391: /*@
392:    SVDGetBV - Obtain the basis vectors objects associated to the singular
393:    value solver object.

395:    Not Collective

397:    Input Parameter:
398: .  svd - singular value solver context obtained from SVDCreate()

400:    Output Parameters:
401: +  V - basis vectors context for right singular vectors
402: -  U - basis vectors context for left singular vectors

404:    Level: advanced

406: .seealso: SVDSetBV()
407: @*/
408: PetscErrorCode SVDGetBV(SVD svd,BV *V,BV *U)
409: {

414:   if (V) {
415:     if (!svd->V) {
416:       BVCreate(PetscObjectComm((PetscObject)svd),&svd->V);
417:       PetscObjectIncrementTabLevel((PetscObject)svd->V,(PetscObject)svd,0);
418:       PetscLogObjectParent((PetscObject)svd,(PetscObject)svd->V);
419:       PetscObjectSetOptions((PetscObject)svd->V,((PetscObject)svd)->options);
420:     }
421:     *V = svd->V;
422:   }
423:   if (U) {
424:     if (!svd->U) {
425:       BVCreate(PetscObjectComm((PetscObject)svd),&svd->U);
426:       PetscObjectIncrementTabLevel((PetscObject)svd->U,(PetscObject)svd,0);
427:       PetscLogObjectParent((PetscObject)svd,(PetscObject)svd->U);
428:       PetscObjectSetOptions((PetscObject)svd->U,((PetscObject)svd)->options);
429:     }
430:     *U = svd->U;
431:   }
432:   return(0);
433: }

435: /*@
436:    SVDSetDS - Associates a direct solver object to the singular value solver.

438:    Collective on svd

440:    Input Parameters:
441: +  svd - singular value solver context obtained from SVDCreate()
442: -  ds  - the direct solver object

444:    Note:
445:    Use SVDGetDS() to retrieve the direct solver context (for example,
446:    to free it at the end of the computations).

448:    Level: advanced

450: .seealso: SVDGetDS()
451: @*/
452: PetscErrorCode SVDSetDS(SVD svd,DS ds)
453: {

460:   PetscObjectReference((PetscObject)ds);
461:   DSDestroy(&svd->ds);
462:   svd->ds = ds;
463:   PetscLogObjectParent((PetscObject)svd,(PetscObject)svd->ds);
464:   return(0);
465: }

467: /*@
468:    SVDGetDS - Obtain the direct solver object associated to the singular value
469:    solver object.

471:    Not Collective

473:    Input Parameters:
474: .  svd - singular value solver context obtained from SVDCreate()

476:    Output Parameter:
477: .  ds - direct solver context

479:    Level: advanced

481: .seealso: SVDSetDS()
482: @*/
483: PetscErrorCode SVDGetDS(SVD svd,DS *ds)
484: {

490:   if (!svd->ds) {
491:     DSCreate(PetscObjectComm((PetscObject)svd),&svd->ds);
492:     PetscObjectIncrementTabLevel((PetscObject)svd->ds,(PetscObject)svd,0);
493:     PetscLogObjectParent((PetscObject)svd,(PetscObject)svd->ds);
494:     PetscObjectSetOptions((PetscObject)svd->ds,((PetscObject)svd)->options);
495:   }
496:   *ds = svd->ds;
497:   return(0);
498: }