Actual source code: svdbasic.c

  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 `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: [](ch:svd), `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 - the singular value solver context

130:    Level: advanced

132: .seealso: [](ch:svd), `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 - the singular value solver context

166:    Level: beginner

168: .seealso: [](ch:svd), `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 `SVDType` for available methods. The default 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: [](ch:svd), `SVDType`
219: @*/
220: PetscErrorCode SVDSetType(SVD svd,SVDType type)
221: {
222:   PetscErrorCode (*r)(SVD);
223:   PetscBool      match;

225:   PetscFunctionBegin;
227:   PetscAssertPointer(type,2);

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

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

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

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

244: /*@
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: .  type - name of `SVD` method

255:    Level: intermediate

257: .seealso: [](ch:svd), `SVDSetType()`
258: @*/
259: PetscErrorCode SVDGetType(SVD svd,SVDType *type)
260: {
261:   PetscFunctionBegin;
263:   PetscAssertPointer(type,2);
264:   *type = ((PetscObject)svd)->type_name;
265:   PetscFunctionReturn(PETSC_SUCCESS);
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:    Note:
278:    `SVDRegister()` may be called multiple times to add several user-defined solvers.

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

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

291:    Level: advanced

293: .seealso: [](ch:svd), `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 - Registers an `SVD` monitor routine that may be accessed with
305:    `SVDMonitorSetFromOptions()`.

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, see `SVDMonitorRegisterFn`
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:    The calling sequence for the given function matches the calling sequence of `SVDMonitorFn`
321:    functions passed to `SVDMonitorSet()` with the additional requirement that its final argument
322:    be a `PetscViewerAndFormat`.

324:    Example Usage:
325: .vb
326:    SVDMonitorRegister("my_monitor",PETSCVIEWERASCII,PETSC_VIEWER_ASCII_INFO_DETAIL,MyMonitor,NULL,NULL);
327: .ve

329:    Then, your monitor can be chosen with the procedural interface via
330: .vb
331:    SVDMonitorSetFromOptions(svd,"-svd_monitor_my_monitor","my_monitor",NULL);
332: .ve
333:    or at runtime via the option `-svd_monitor_my_monitor`.

335:    Level: advanced

337: .seealso: [](ch:svd), `SVDMonitorSet()`, `SVDMonitorRegisterAll()`, `SVDMonitorSetFromOptions()`
338: @*/
339: PetscErrorCode SVDMonitorRegister(const char name[],PetscViewerType vtype,PetscViewerFormat format,SVDMonitorRegisterFn *monitor,SVDMonitorRegisterCreateFn *create,SVDMonitorRegisterDestroyFn *destroy)
340: {
341:   char           key[PETSC_MAX_PATH_LEN];

343:   PetscFunctionBegin;
344:   PetscCall(SVDInitializePackage());
345:   PetscCall(SlepcMonitorMakeKey_Internal(name,vtype,format,key));
346:   PetscCall(PetscFunctionListAdd(&SVDMonitorList,key,monitor));
347:   if (create)  PetscCall(PetscFunctionListAdd(&SVDMonitorCreateList,key,create));
348:   if (destroy) PetscCall(PetscFunctionListAdd(&SVDMonitorDestroyList,key,destroy));
349:   PetscFunctionReturn(PETSC_SUCCESS);
350: }

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

355:    Collective

357:    Input Parameters:
358: +  svd - the singular value solver context
359: .  V   - the basis vectors object for right singular vectors
360: -  U   - the basis vectors object for left singular vectors

362:    Note:
363:    Use SVDGetBV() to retrieve the basis vectors contexts (for example,
364:    to free them at the end of the computations).

366:    Level: advanced

368: .seealso: [](ch:svd), `SVDGetBV()`
369: @*/
370: PetscErrorCode SVDSetBV(SVD svd,BV V,BV U)
371: {
372:   PetscFunctionBegin;
374:   if (V) {
376:     PetscCheckSameComm(svd,1,V,2);
377:     PetscCall(PetscObjectReference((PetscObject)V));
378:     PetscCall(BVDestroy(&svd->V));
379:     svd->V = V;
380:   }
381:   if (U) {
383:     PetscCheckSameComm(svd,1,U,3);
384:     PetscCall(PetscObjectReference((PetscObject)U));
385:     PetscCall(BVDestroy(&svd->U));
386:     svd->U = U;
387:   }
388:   PetscFunctionReturn(PETSC_SUCCESS);
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 - the singular value solver context

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: [](ch:svd), `SVDSetBV()`
407: @*/
408: PetscErrorCode SVDGetBV(SVD svd,BV *V,BV *U)
409: {
410:   PetscFunctionBegin;
412:   if (V) {
413:     if (!svd->V) {
414:       PetscCall(BVCreate(PetscObjectComm((PetscObject)svd),&svd->V));
415:       PetscCall(PetscObjectIncrementTabLevel((PetscObject)svd->V,(PetscObject)svd,0));
416:       PetscCall(PetscObjectSetOptions((PetscObject)svd->V,((PetscObject)svd)->options));
417:     }
418:     *V = svd->V;
419:   }
420:   if (U) {
421:     if (!svd->U) {
422:       PetscCall(BVCreate(PetscObjectComm((PetscObject)svd),&svd->U));
423:       PetscCall(PetscObjectIncrementTabLevel((PetscObject)svd->U,(PetscObject)svd,0));
424:       PetscCall(PetscObjectSetOptions((PetscObject)svd->U,((PetscObject)svd)->options));
425:     }
426:     *U = svd->U;
427:   }
428:   PetscFunctionReturn(PETSC_SUCCESS);
429: }

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

434:    Collective

436:    Input Parameters:
437: +  svd - the singular value solver context
438: -  ds  - the direct solver object

440:    Note:
441:    Use SVDGetDS() to retrieve the direct solver context (for example,
442:    to free it at the end of the computations).

444:    Level: advanced

446: .seealso: [](ch:svd), `SVDGetDS()`
447: @*/
448: PetscErrorCode SVDSetDS(SVD svd,DS ds)
449: {
450:   PetscFunctionBegin;
453:   PetscCheckSameComm(svd,1,ds,2);
454:   PetscCall(PetscObjectReference((PetscObject)ds));
455:   PetscCall(DSDestroy(&svd->ds));
456:   svd->ds = ds;
457:   PetscFunctionReturn(PETSC_SUCCESS);
458: }

460: /*@
461:    SVDGetDS - Obtain the direct solver object associated to the singular value
462:    solver object.

464:    Not Collective

466:    Input Parameter:
467: .  svd - the singular value solver context

469:    Output Parameter:
470: .  ds - direct solver context

472:    Level: advanced

474: .seealso: [](ch:svd), `SVDSetDS()`
475: @*/
476: PetscErrorCode SVDGetDS(SVD svd,DS *ds)
477: {
478:   PetscFunctionBegin;
480:   PetscAssertPointer(ds,2);
481:   if (!svd->ds) {
482:     PetscCall(DSCreate(PetscObjectComm((PetscObject)svd),&svd->ds));
483:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)svd->ds,(PetscObject)svd,0));
484:     PetscCall(PetscObjectSetOptions((PetscObject)svd->ds,((PetscObject)svd)->options));
485:   }
486:   *ds = svd->ds;
487:   PetscFunctionReturn(PETSC_SUCCESS);
488: }