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 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 - Registers an SVD monitor routine that may be accessed with SVDMonitorSetFromOptions().

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

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

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

333:    Level: advanced

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

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

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

353:    Collective

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

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

364:    Level: advanced

366: .seealso: `SVDGetBV()`
367: @*/
368: PetscErrorCode SVDSetBV(SVD svd,BV V,BV U)
369: {
370:   PetscFunctionBegin;
372:   if (V) {
374:     PetscCheckSameComm(svd,1,V,2);
375:     PetscCall(PetscObjectReference((PetscObject)V));
376:     PetscCall(BVDestroy(&svd->V));
377:     svd->V = V;
378:   }
379:   if (U) {
381:     PetscCheckSameComm(svd,1,U,3);
382:     PetscCall(PetscObjectReference((PetscObject)U));
383:     PetscCall(BVDestroy(&svd->U));
384:     svd->U = U;
385:   }
386:   PetscFunctionReturn(PETSC_SUCCESS);
387: }

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

393:    Not Collective

395:    Input Parameter:
396: .  svd - singular value solver context obtained from SVDCreate()

398:    Output Parameters:
399: +  V - basis vectors context for right singular vectors
400: -  U - basis vectors context for left singular vectors

402:    Level: advanced

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

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

432:    Collective

434:    Input Parameters:
435: +  svd - singular value solver context obtained from SVDCreate()
436: -  ds  - the direct solver object

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

442:    Level: advanced

444: .seealso: `SVDGetDS()`
445: @*/
446: PetscErrorCode SVDSetDS(SVD svd,DS ds)
447: {
448:   PetscFunctionBegin;
451:   PetscCheckSameComm(svd,1,ds,2);
452:   PetscCall(PetscObjectReference((PetscObject)ds));
453:   PetscCall(DSDestroy(&svd->ds));
454:   svd->ds = ds;
455:   PetscFunctionReturn(PETSC_SUCCESS);
456: }

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

462:    Not Collective

464:    Input Parameters:
465: .  svd - singular value solver context obtained from SVDCreate()

467:    Output Parameter:
468: .  ds - direct solver context

470:    Level: advanced

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