Actual source code: svdbasic.c

slepc-main 2024-12-17
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:    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 = NULL;
 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: .  outsvd - 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: {
 50:   SVD            svd;

 52:   PetscFunctionBegin;
 53:   PetscAssertPointer(outsvd,2);
 54:   PetscCall(SVDInitializePackage());
 55:   PetscCall(SlepcHeaderCreate(svd,SVD_CLASSID,"SVD","Singular Value Decomposition","SVD",comm,SVDDestroy,SVDView));

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

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

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

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

116:   PetscCall(PetscNew(&svd->sc));
117:   *outsvd = svd;
118:   PetscFunctionReturn(PETSC_SUCCESS);
119: }

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

125:    Collective

127:    Input Parameter:
128: .  svd - singular value solver context obtained from SVDCreate()

130:    Level: advanced

132: .seealso: SVDDestroy()
133: @*/
134: PetscErrorCode SVDReset(SVD svd)
135: {
136:   PetscFunctionBegin;
138:   if (!svd) PetscFunctionReturn(PETSC_SUCCESS);
139:   PetscTryTypeMethod(svd,reset);
140:   PetscCall(MatDestroy(&svd->OP));
141:   PetscCall(MatDestroy(&svd->OPb));
142:   PetscCall(VecDestroy(&svd->omega));
143:   PetscCall(MatDestroy(&svd->A));
144:   PetscCall(MatDestroy(&svd->B));
145:   PetscCall(MatDestroy(&svd->AT));
146:   PetscCall(MatDestroy(&svd->BT));
147:   PetscCall(BVDestroy(&svd->U));
148:   PetscCall(BVDestroy(&svd->V));
149:   PetscCall(VecDestroyVecs(svd->nworkl,&svd->workl));
150:   svd->nworkl = 0;
151:   PetscCall(VecDestroyVecs(svd->nworkr,&svd->workr));
152:   svd->nworkr = 0;
153:   svd->swapped = PETSC_FALSE;
154:   svd->state = SVD_STATE_INITIAL;
155:   PetscFunctionReturn(PETSC_SUCCESS);
156: }

158: /*@
159:    SVDDestroy - Destroys the SVD context.

161:    Collective

163:    Input Parameter:
164: .  svd - singular value solver context obtained from SVDCreate()

166:    Level: beginner

168: .seealso: SVDCreate(), SVDSetUp(), SVDSolve()
169: @*/
170: PetscErrorCode SVDDestroy(SVD *svd)
171: {
172:   PetscFunctionBegin;
173:   if (!*svd) PetscFunctionReturn(PETSC_SUCCESS);
175:   if (--((PetscObject)*svd)->refct > 0) { *svd = NULL; PetscFunctionReturn(PETSC_SUCCESS); }
176:   PetscCall(SVDReset(*svd));
177:   PetscTryTypeMethod(*svd,destroy);
178:   if ((*svd)->sigma) PetscCall(PetscFree3((*svd)->sigma,(*svd)->perm,(*svd)->errest));
179:   if ((*svd)->sign) PetscCall(PetscFree((*svd)->sign));
180:   PetscCall(DSDestroy(&(*svd)->ds));
181:   PetscCall(PetscFree((*svd)->sc));
182:   /* just in case the initial vectors have not been used */
183:   PetscCall(SlepcBasisDestroy_Private(&(*svd)->nini,&(*svd)->IS));
184:   PetscCall(SlepcBasisDestroy_Private(&(*svd)->ninil,&(*svd)->ISL));
185:   if ((*svd)->convergeddestroy) PetscCall((*(*svd)->convergeddestroy)(&(*svd)->convergedctx));
186:   if ((*svd)->stoppingdestroy) PetscCall((*(*svd)->stoppingdestroy)(&(*svd)->stoppingctx));
187:   PetscCall(SVDMonitorCancel(*svd));
188:   PetscCall(PetscHeaderDestroy(svd));
189:   PetscFunctionReturn(PETSC_SUCCESS);
190: }

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

195:    Logically Collective

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

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

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

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

217:    Level: intermediate

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

226:   PetscFunctionBegin;
228:   PetscAssertPointer(type,2);

230:   PetscCall(PetscObjectTypeCompare((PetscObject)svd,type,&match));
231:   if (match) PetscFunctionReturn(PETSC_SUCCESS);

233:   PetscCall(PetscFunctionListFind(SVDList,type,&r));
234:   PetscCheck(r,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown SVD type given: %s",type);

236:   PetscTryTypeMethod(svd,destroy);
237:   PetscCall(PetscMemzero(svd->ops,sizeof(struct _SVDOps)));

239:   svd->state = SVD_STATE_INITIAL;
240:   PetscCall(PetscObjectChangeTypeName((PetscObject)svd,type));
241:   PetscCall((*r)(svd));
242:   PetscFunctionReturn(PETSC_SUCCESS);
243: }

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

248:    Not Collective

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

253:    Output Parameter:
254: .  type - name of SVD method

256:    Level: intermediate

258: .seealso: SVDSetType()
259: @*/
260: PetscErrorCode SVDGetType(SVD svd,SVDType *type)
261: {
262:   PetscFunctionBegin;
264:   PetscAssertPointer(type,2);
265:   *type = ((PetscObject)svd)->type_name;
266:   PetscFunctionReturn(PETSC_SUCCESS);
267: }

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

272:    Not Collective

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

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

281:    Example Usage:
282: .vb
283:     SVDRegister("my_solver",MySolverCreate);
284: .ve

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

291:    Level: advanced

293: .seealso: SVDRegisterAll()
294: @*/
295: PetscErrorCode SVDRegister(const char *name,PetscErrorCode (*function)(SVD))
296: {
297:   PetscFunctionBegin;
298:   PetscCall(SVDInitializePackage());
299:   PetscCall(PetscFunctionListAdd(&SVDList,name,function));
300:   PetscFunctionReturn(PETSC_SUCCESS);
301: }

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

306:    Not Collective

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

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

319:    Example Usage:
320: .vb
321:    SVDMonitorRegister("my_monitor",PETSCVIEWERASCII,PETSC_VIEWER_ASCII_INFO_DETAIL,MyMonitor,NULL,NULL);
322: .ve

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

329:    Level: advanced

331: .seealso: SVDMonitorRegisterAll()
332: @*/
333: 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**))
334: {
335:   char           key[PETSC_MAX_PATH_LEN];

337:   PetscFunctionBegin;
338:   PetscCall(SVDInitializePackage());
339:   PetscCall(SlepcMonitorMakeKey_Internal(name,vtype,format,key));
340:   PetscCall(PetscFunctionListAdd(&SVDMonitorList,key,monitor));
341:   if (create)  PetscCall(PetscFunctionListAdd(&SVDMonitorCreateList,key,create));
342:   if (destroy) PetscCall(PetscFunctionListAdd(&SVDMonitorDestroyList,key,destroy));
343:   PetscFunctionReturn(PETSC_SUCCESS);
344: }

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

349:    Collective

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

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

360:    Level: advanced

362: .seealso: SVDGetBV()
363: @*/
364: PetscErrorCode SVDSetBV(SVD svd,BV V,BV U)
365: {
366:   PetscFunctionBegin;
368:   if (V) {
370:     PetscCheckSameComm(svd,1,V,2);
371:     PetscCall(PetscObjectReference((PetscObject)V));
372:     PetscCall(BVDestroy(&svd->V));
373:     svd->V = V;
374:   }
375:   if (U) {
377:     PetscCheckSameComm(svd,1,U,3);
378:     PetscCall(PetscObjectReference((PetscObject)U));
379:     PetscCall(BVDestroy(&svd->U));
380:     svd->U = U;
381:   }
382:   PetscFunctionReturn(PETSC_SUCCESS);
383: }

385: /*@
386:    SVDGetBV - Obtain the basis vectors objects associated to the singular
387:    value solver object.

389:    Not Collective

391:    Input Parameter:
392: .  svd - singular value solver context obtained from SVDCreate()

394:    Output Parameters:
395: +  V - basis vectors context for right singular vectors
396: -  U - basis vectors context for left singular vectors

398:    Level: advanced

400: .seealso: SVDSetBV()
401: @*/
402: PetscErrorCode SVDGetBV(SVD svd,BV *V,BV *U)
403: {
404:   PetscFunctionBegin;
406:   if (V) {
407:     if (!svd->V) {
408:       PetscCall(BVCreate(PetscObjectComm((PetscObject)svd),&svd->V));
409:       PetscCall(PetscObjectIncrementTabLevel((PetscObject)svd->V,(PetscObject)svd,0));
410:       PetscCall(PetscObjectSetOptions((PetscObject)svd->V,((PetscObject)svd)->options));
411:     }
412:     *V = svd->V;
413:   }
414:   if (U) {
415:     if (!svd->U) {
416:       PetscCall(BVCreate(PetscObjectComm((PetscObject)svd),&svd->U));
417:       PetscCall(PetscObjectIncrementTabLevel((PetscObject)svd->U,(PetscObject)svd,0));
418:       PetscCall(PetscObjectSetOptions((PetscObject)svd->U,((PetscObject)svd)->options));
419:     }
420:     *U = svd->U;
421:   }
422:   PetscFunctionReturn(PETSC_SUCCESS);
423: }

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

428:    Collective

430:    Input Parameters:
431: +  svd - singular value solver context obtained from SVDCreate()
432: -  ds  - the direct solver object

434:    Note:
435:    Use SVDGetDS() to retrieve the direct solver context (for example,
436:    to free it at the end of the computations).

438:    Level: advanced

440: .seealso: SVDGetDS()
441: @*/
442: PetscErrorCode SVDSetDS(SVD svd,DS ds)
443: {
444:   PetscFunctionBegin;
447:   PetscCheckSameComm(svd,1,ds,2);
448:   PetscCall(PetscObjectReference((PetscObject)ds));
449:   PetscCall(DSDestroy(&svd->ds));
450:   svd->ds = ds;
451:   PetscFunctionReturn(PETSC_SUCCESS);
452: }

454: /*@
455:    SVDGetDS - Obtain the direct solver object associated to the singular value
456:    solver object.

458:    Not Collective

460:    Input Parameters:
461: .  svd - singular value solver context obtained from SVDCreate()

463:    Output Parameter:
464: .  ds - direct solver context

466:    Level: advanced

468: .seealso: SVDSetDS()
469: @*/
470: PetscErrorCode SVDGetDS(SVD svd,DS *ds)
471: {
472:   PetscFunctionBegin;
474:   PetscAssertPointer(ds,2);
475:   if (!svd->ds) {
476:     PetscCall(DSCreate(PetscObjectComm((PetscObject)svd),&svd->ds));
477:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)svd->ds,(PetscObject)svd,0));
478:     PetscCall(PetscObjectSetOptions((PetscObject)svd->ds,((PetscObject)svd)->options));
479:   }
480:   *ds = svd->ds;
481:   PetscFunctionReturn(PETSC_SUCCESS);
482: }