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 `Vec`s and `Mat`s.
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: 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 type - sets the method; use `-help` for a list of available methods
204: Notes:
205: See `SVDType` for available methods. The default is `SVDCROSS`.
207: Normally, it is best to use the `SVDSetFromOptions()` command and
208: then set the `SVD` type from the options database rather than by using
209: this routine. Using the options database provides the user with
210: maximum flexibility in evaluating the different available methods.
211: The `SVDSetType()` routine is provided for those situations where it
212: is necessary to set the iterative solver independently of the command
213: line or options database.
215: Level: intermediate
217: .seealso: [](ch:svd), `SVDType`
218: @*/
219: PetscErrorCode SVDSetType(SVD svd,SVDType type)
220: {
221: PetscErrorCode (*r)(SVD);
222: PetscBool match;
224: PetscFunctionBegin;
226: PetscAssertPointer(type,2);
228: PetscCall(PetscObjectTypeCompare((PetscObject)svd,type,&match));
229: if (match) PetscFunctionReturn(PETSC_SUCCESS);
231: PetscCall(PetscFunctionListFind(SVDList,type,&r));
232: PetscCheck(r,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown SVD type given: %s",type);
234: PetscTryTypeMethod(svd,destroy);
235: PetscCall(PetscMemzero(svd->ops,sizeof(struct _SVDOps)));
237: svd->state = SVD_STATE_INITIAL;
238: PetscCall(PetscObjectChangeTypeName((PetscObject)svd,type));
239: PetscCall((*r)(svd));
240: PetscFunctionReturn(PETSC_SUCCESS);
241: }
243: /*@
244: SVDGetType - Gets the `SVD` type as a string from the `SVD` object.
246: Not Collective
248: Input Parameter:
249: . svd - the singular value solver context
251: Output Parameter:
252: . type - name of `SVD` method
254: Note:
255: `type` should not be retained for later use as it will be an invalid pointer
256: if the `SVDType` of `svd` is changed.
258: Level: intermediate
260: .seealso: [](ch:svd), `SVDSetType()`, `PetscObjectTypeCompare()`, `PetscObjectTypeCompareAny()`
261: @*/
262: PetscErrorCode SVDGetType(SVD svd,SVDType *type)
263: {
264: PetscFunctionBegin;
266: PetscAssertPointer(type,2);
267: *type = ((PetscObject)svd)->type_name;
268: PetscFunctionReturn(PETSC_SUCCESS);
269: }
271: /*@C
272: SVDRegister - Adds a method to the singular value solver package.
274: Not Collective
276: Input Parameters:
277: + name - name of a new user-defined solver
278: - function - routine to create the solver context
280: Note:
281: `SVDRegister()` may be called multiple times to add several user-defined solvers.
283: Example Usage:
284: .vb
285: SVDRegister("my_solver",MySolverCreate);
286: .ve
288: Then, your solver can be chosen with the procedural interface via
289: .vb
290: SVDSetType(svd,"my_solver")
291: .ve
292: or at runtime via the option `-svd_type my_solver`.
294: Level: advanced
296: .seealso: [](ch:svd), `SVDRegisterAll()`
297: @*/
298: PetscErrorCode SVDRegister(const char *name,PetscErrorCode (*function)(SVD))
299: {
300: PetscFunctionBegin;
301: PetscCall(SVDInitializePackage());
302: PetscCall(PetscFunctionListAdd(&SVDList,name,function));
303: PetscFunctionReturn(PETSC_SUCCESS);
304: }
306: /*@C
307: SVDMonitorRegister - Registers an `SVD` monitor routine that may be accessed with
308: `SVDMonitorSetFromOptions()`.
310: Not Collective
312: Input Parameters:
313: + name - name of a new monitor routine
314: . vtype - a `PetscViewerType` for the output
315: . format - a `PetscViewerFormat` for the output
316: . monitor - monitor routine, see `SVDMonitorRegisterFn`
317: . create - creation routine, or `NULL`
318: - destroy - destruction routine, or `NULL`
320: Notes:
321: `SVDMonitorRegister()` may be called multiple times to add several user-defined monitors.
323: The calling sequence for the given function matches the calling sequence of `SVDMonitorFn`
324: functions passed to `SVDMonitorSet()` with the additional requirement that its final argument
325: be a `PetscViewerAndFormat`.
327: Example Usage:
328: .vb
329: SVDMonitorRegister("my_monitor",PETSCVIEWERASCII,PETSC_VIEWER_ASCII_INFO_DETAIL,MyMonitor,NULL,NULL);
330: .ve
332: Then, your monitor can be chosen with the procedural interface via
333: .vb
334: SVDMonitorSetFromOptions(svd,"-svd_monitor_my_monitor","my_monitor",NULL);
335: .ve
336: or at runtime via the option `-svd_monitor_my_monitor`.
338: Level: advanced
340: .seealso: [](ch:svd), `SVDMonitorSet()`, `SVDMonitorRegisterAll()`, `SVDMonitorSetFromOptions()`
341: @*/
342: PetscErrorCode SVDMonitorRegister(const char name[],PetscViewerType vtype,PetscViewerFormat format,SVDMonitorRegisterFn *monitor,SVDMonitorRegisterCreateFn *create,SVDMonitorRegisterDestroyFn *destroy)
343: {
344: char key[PETSC_MAX_PATH_LEN];
346: PetscFunctionBegin;
347: PetscCall(SVDInitializePackage());
348: PetscCall(SlepcMonitorMakeKey_Internal(name,vtype,format,key));
349: PetscCall(PetscFunctionListAdd(&SVDMonitorList,key,monitor));
350: if (create) PetscCall(PetscFunctionListAdd(&SVDMonitorCreateList,key,create));
351: if (destroy) PetscCall(PetscFunctionListAdd(&SVDMonitorDestroyList,key,destroy));
352: PetscFunctionReturn(PETSC_SUCCESS);
353: }
355: /*@
356: SVDSetBV - Associates basis vectors objects to the singular value solver.
358: Collective
360: Input Parameters:
361: + svd - the singular value solver context
362: . V - the basis vectors object for right singular vectors
363: - U - the basis vectors object for left singular vectors
365: Note:
366: Use `SVDGetBV()` to retrieve the basis vectors contexts (for example,
367: to free them at the end of the computations).
369: Level: advanced
371: .seealso: [](ch:svd), `SVDGetBV()`
372: @*/
373: PetscErrorCode SVDSetBV(SVD svd,BV V,BV U)
374: {
375: PetscFunctionBegin;
377: if (V) {
379: PetscCheckSameComm(svd,1,V,2);
380: PetscCall(PetscObjectReference((PetscObject)V));
381: PetscCall(BVDestroy(&svd->V));
382: svd->V = V;
383: }
384: if (U) {
386: PetscCheckSameComm(svd,1,U,3);
387: PetscCall(PetscObjectReference((PetscObject)U));
388: PetscCall(BVDestroy(&svd->U));
389: svd->U = U;
390: }
391: PetscFunctionReturn(PETSC_SUCCESS);
392: }
394: /*@
395: SVDGetBV - Obtain the basis vectors objects associated to the singular
396: value solver object.
398: Not Collective
400: Input Parameter:
401: . svd - the singular value solver context
403: Output Parameters:
404: + V - basis vectors context for right singular vectors
405: - U - basis vectors context for left singular vectors
407: Level: advanced
409: .seealso: [](ch:svd), `SVDSetBV()`
410: @*/
411: PetscErrorCode SVDGetBV(SVD svd,BV *V,BV *U)
412: {
413: PetscFunctionBegin;
415: if (V) {
416: if (!svd->V) {
417: PetscCall(BVCreate(PetscObjectComm((PetscObject)svd),&svd->V));
418: PetscCall(PetscObjectIncrementTabLevel((PetscObject)svd->V,(PetscObject)svd,0));
419: PetscCall(PetscObjectSetOptions((PetscObject)svd->V,((PetscObject)svd)->options));
420: }
421: *V = svd->V;
422: }
423: if (U) {
424: if (!svd->U) {
425: PetscCall(BVCreate(PetscObjectComm((PetscObject)svd),&svd->U));
426: PetscCall(PetscObjectIncrementTabLevel((PetscObject)svd->U,(PetscObject)svd,0));
427: PetscCall(PetscObjectSetOptions((PetscObject)svd->U,((PetscObject)svd)->options));
428: }
429: *U = svd->U;
430: }
431: PetscFunctionReturn(PETSC_SUCCESS);
432: }
434: /*@
435: SVDSetDS - Associates a direct solver object to the singular value solver.
437: Collective
439: Input Parameters:
440: + svd - the singular value solver context
441: - ds - the direct solver object
443: Note:
444: Use `SVDGetDS()` to retrieve the direct solver context (for example,
445: to free it at the end of the computations).
447: Level: advanced
449: .seealso: [](ch:svd), `SVDGetDS()`
450: @*/
451: PetscErrorCode SVDSetDS(SVD svd,DS ds)
452: {
453: PetscFunctionBegin;
456: PetscCheckSameComm(svd,1,ds,2);
457: PetscCall(PetscObjectReference((PetscObject)ds));
458: PetscCall(DSDestroy(&svd->ds));
459: svd->ds = ds;
460: PetscFunctionReturn(PETSC_SUCCESS);
461: }
463: /*@
464: SVDGetDS - Obtain the direct solver object associated to the singular value
465: solver object.
467: Not Collective
469: Input Parameter:
470: . svd - the singular value solver context
472: Output Parameter:
473: . ds - direct solver context
475: Level: advanced
477: .seealso: [](ch:svd), `SVDSetDS()`
478: @*/
479: PetscErrorCode SVDGetDS(SVD svd,DS *ds)
480: {
481: PetscFunctionBegin;
483: PetscAssertPointer(ds,2);
484: if (!svd->ds) {
485: PetscCall(DSCreate(PetscObjectComm((PetscObject)svd),&svd->ds));
486: PetscCall(PetscObjectIncrementTabLevel((PetscObject)svd->ds,(PetscObject)svd,0));
487: PetscCall(PetscObjectSetOptions((PetscObject)svd->ds,((PetscObject)svd)->options));
488: }
489: *ds = svd->ds;
490: PetscFunctionReturn(PETSC_SUCCESS);
491: }