Actual source code: svdbasic.c
slepc-3.18.2 2023-01-26
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;
53: *outsvd = NULL;
54: SVDInitializePackage();
55: 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_DEFAULT;
61: svd->nsv = 1;
62: svd->ncv = PETSC_DEFAULT;
63: svd->mpd = PETSC_DEFAULT;
64: svd->nini = 0;
65: svd->ninil = 0;
66: svd->tol = PETSC_DEFAULT;
67: svd->conv = (SVDConv)-1;
68: svd->stop = SVD_STOP_BASIC;
69: svd->which = SVD_LARGEST;
70: svd->problem_type = (SVDProblemType)0;
71: svd->impltrans = PETSC_FALSE;
72: svd->trackall = PETSC_FALSE;
74: svd->converged = NULL;
75: svd->convergeduser = NULL;
76: svd->convergeddestroy = NULL;
77: svd->stopping = SVDStoppingBasic;
78: svd->stoppinguser = NULL;
79: svd->stoppingdestroy = NULL;
80: svd->convergedctx = NULL;
81: svd->stoppingctx = NULL;
82: svd->numbermonitors = 0;
84: svd->ds = NULL;
85: svd->U = NULL;
86: svd->V = NULL;
87: svd->A = NULL;
88: svd->B = NULL;
89: svd->AT = NULL;
90: svd->BT = NULL;
91: svd->IS = NULL;
92: svd->ISL = NULL;
93: svd->sigma = NULL;
94: svd->errest = NULL;
95: svd->sign = 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: PetscNew(&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: {
135: if (!svd) return 0;
136: PetscTryTypeMethod(svd,reset);
137: MatDestroy(&svd->OP);
138: MatDestroy(&svd->OPb);
139: VecDestroy(&svd->omega);
140: MatDestroy(&svd->A);
141: MatDestroy(&svd->B);
142: MatDestroy(&svd->AT);
143: MatDestroy(&svd->BT);
144: BVDestroy(&svd->U);
145: BVDestroy(&svd->V);
146: VecDestroyVecs(svd->nworkl,&svd->workl);
147: svd->nworkl = 0;
148: VecDestroyVecs(svd->nworkr,&svd->workr);
149: svd->nworkr = 0;
150: svd->swapped = PETSC_FALSE;
151: svd->state = SVD_STATE_INITIAL;
152: return 0;
153: }
155: /*@C
156: SVDDestroy - Destroys the SVD context.
158: Collective on svd
160: Input Parameter:
161: . svd - singular value solver context obtained from SVDCreate()
163: Level: beginner
165: .seealso: SVDCreate(), SVDSetUp(), SVDSolve()
166: @*/
167: PetscErrorCode SVDDestroy(SVD *svd)
168: {
169: if (!*svd) return 0;
171: if (--((PetscObject)(*svd))->refct > 0) { *svd = NULL; return 0; }
172: SVDReset(*svd);
173: PetscTryTypeMethod(*svd,destroy);
174: if ((*svd)->sigma) PetscFree3((*svd)->sigma,(*svd)->perm,(*svd)->errest);
175: if ((*svd)->sign) PetscFree((*svd)->sign);
176: DSDestroy(&(*svd)->ds);
177: PetscFree((*svd)->sc);
178: /* just in case the initial vectors have not been used */
179: SlepcBasisDestroy_Private(&(*svd)->nini,&(*svd)->IS);
180: SlepcBasisDestroy_Private(&(*svd)->ninil,&(*svd)->ISL);
181: SVDMonitorCancel(*svd);
182: PetscHeaderDestroy(svd);
183: return 0;
184: }
186: /*@C
187: SVDSetType - Selects the particular solver to be used in the SVD object.
189: Logically Collective on svd
191: Input Parameters:
192: + svd - the singular value solver context
193: - type - a known method
195: Options Database Key:
196: . -svd_type <method> - Sets the method; use -help for a list
197: of available methods
199: Notes:
200: See "slepc/include/slepcsvd.h" for available methods. The default
201: is SVDCROSS.
203: Normally, it is best to use the SVDSetFromOptions() command and
204: then set the SVD type from the options database rather than by using
205: this routine. Using the options database provides the user with
206: maximum flexibility in evaluating the different available methods.
207: The SVDSetType() routine is provided for those situations where it
208: is necessary to set the iterative solver independently of the command
209: line or options database.
211: Level: intermediate
213: .seealso: SVDType
214: @*/
215: PetscErrorCode SVDSetType(SVD svd,SVDType type)
216: {
217: PetscErrorCode (*r)(SVD);
218: PetscBool match;
223: PetscObjectTypeCompare((PetscObject)svd,type,&match);
224: if (match) return 0;
226: PetscFunctionListFind(SVDList,type,&r);
229: PetscTryTypeMethod(svd,destroy);
230: PetscMemzero(svd->ops,sizeof(struct _SVDOps));
232: svd->state = SVD_STATE_INITIAL;
233: PetscObjectChangeTypeName((PetscObject)svd,type);
234: (*r)(svd);
235: return 0;
236: }
238: /*@C
239: SVDGetType - Gets the SVD type as a string from the SVD object.
241: Not Collective
243: Input Parameter:
244: . svd - the singular value solver context
246: Output Parameter:
247: . type - name of SVD method
249: Level: intermediate
251: .seealso: SVDSetType()
252: @*/
253: PetscErrorCode SVDGetType(SVD svd,SVDType *type)
254: {
257: *type = ((PetscObject)svd)->type_name;
258: return 0;
259: }
261: /*@C
262: SVDRegister - Adds a method to the singular value solver package.
264: Not Collective
266: Input Parameters:
267: + name - name of a new user-defined solver
268: - function - routine to create the solver context
270: Notes:
271: SVDRegister() may be called multiple times to add several user-defined solvers.
273: Sample usage:
274: .vb
275: SVDRegister("my_solver",MySolverCreate);
276: .ve
278: Then, your solver can be chosen with the procedural interface via
279: $ SVDSetType(svd,"my_solver")
280: or at runtime via the option
281: $ -svd_type my_solver
283: Level: advanced
285: .seealso: SVDRegisterAll()
286: @*/
287: PetscErrorCode SVDRegister(const char *name,PetscErrorCode (*function)(SVD))
288: {
289: SVDInitializePackage();
290: PetscFunctionListAdd(&SVDList,name,function);
291: return 0;
292: }
294: /*@C
295: SVDMonitorRegister - Adds SVD monitor routine.
297: Not Collective
299: Input Parameters:
300: + name - name of a new monitor routine
301: . vtype - a PetscViewerType for the output
302: . format - a PetscViewerFormat for the output
303: . monitor - monitor routine
304: . create - creation routine, or NULL
305: - destroy - destruction routine, or NULL
307: Notes:
308: SVDMonitorRegister() may be called multiple times to add several user-defined monitors.
310: Sample usage:
311: .vb
312: SVDMonitorRegister("my_monitor",PETSCVIEWERASCII,PETSC_VIEWER_ASCII_INFO_DETAIL,MyMonitor,NULL,NULL);
313: .ve
315: Then, your monitor can be chosen with the procedural interface via
316: $ SVDMonitorSetFromOptions(svd,"-svd_monitor_my_monitor","my_monitor",NULL)
317: or at runtime via the option
318: $ -svd_monitor_my_monitor
320: Level: advanced
322: .seealso: SVDMonitorRegisterAll()
323: @*/
324: 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**))
325: {
326: char key[PETSC_MAX_PATH_LEN];
328: SVDInitializePackage();
329: SlepcMonitorMakeKey_Internal(name,vtype,format,key);
330: PetscFunctionListAdd(&SVDMonitorList,key,monitor);
331: if (create) PetscFunctionListAdd(&SVDMonitorCreateList,key,create);
332: if (destroy) PetscFunctionListAdd(&SVDMonitorDestroyList,key,destroy);
333: return 0;
334: }
336: /*@
337: SVDSetBV - Associates basis vectors objects to the singular value solver.
339: Collective on svd
341: Input Parameters:
342: + svd - singular value solver context obtained from SVDCreate()
343: . V - the basis vectors object for right singular vectors
344: - U - the basis vectors object for left singular vectors
346: Note:
347: Use SVDGetBV() to retrieve the basis vectors contexts (for example,
348: to free them at the end of the computations).
350: Level: advanced
352: .seealso: SVDGetBV()
353: @*/
354: PetscErrorCode SVDSetBV(SVD svd,BV V,BV U)
355: {
357: if (V) {
360: PetscObjectReference((PetscObject)V);
361: BVDestroy(&svd->V);
362: svd->V = V;
363: }
364: if (U) {
367: PetscObjectReference((PetscObject)U);
368: BVDestroy(&svd->U);
369: svd->U = U;
370: }
371: return 0;
372: }
374: /*@
375: SVDGetBV - Obtain the basis vectors objects associated to the singular
376: value solver object.
378: Not Collective
380: Input Parameter:
381: . svd - singular value solver context obtained from SVDCreate()
383: Output Parameters:
384: + V - basis vectors context for right singular vectors
385: - U - basis vectors context for left singular vectors
387: Level: advanced
389: .seealso: SVDSetBV()
390: @*/
391: PetscErrorCode SVDGetBV(SVD svd,BV *V,BV *U)
392: {
394: if (V) {
395: if (!svd->V) {
396: BVCreate(PetscObjectComm((PetscObject)svd),&svd->V);
397: PetscObjectIncrementTabLevel((PetscObject)svd->V,(PetscObject)svd,0);
398: PetscObjectSetOptions((PetscObject)svd->V,((PetscObject)svd)->options);
399: }
400: *V = svd->V;
401: }
402: if (U) {
403: if (!svd->U) {
404: BVCreate(PetscObjectComm((PetscObject)svd),&svd->U);
405: PetscObjectIncrementTabLevel((PetscObject)svd->U,(PetscObject)svd,0);
406: PetscObjectSetOptions((PetscObject)svd->U,((PetscObject)svd)->options);
407: }
408: *U = svd->U;
409: }
410: return 0;
411: }
413: /*@
414: SVDSetDS - Associates a direct solver object to the singular value solver.
416: Collective on svd
418: Input Parameters:
419: + svd - singular value solver context obtained from SVDCreate()
420: - ds - the direct solver object
422: Note:
423: Use SVDGetDS() to retrieve the direct solver context (for example,
424: to free it at the end of the computations).
426: Level: advanced
428: .seealso: SVDGetDS()
429: @*/
430: PetscErrorCode SVDSetDS(SVD svd,DS ds)
431: {
435: PetscObjectReference((PetscObject)ds);
436: DSDestroy(&svd->ds);
437: svd->ds = ds;
438: return 0;
439: }
441: /*@
442: SVDGetDS - Obtain the direct solver object associated to the singular value
443: solver object.
445: Not Collective
447: Input Parameters:
448: . svd - singular value solver context obtained from SVDCreate()
450: Output Parameter:
451: . ds - direct solver context
453: Level: advanced
455: .seealso: SVDSetDS()
456: @*/
457: PetscErrorCode SVDGetDS(SVD svd,DS *ds)
458: {
461: if (!svd->ds) {
462: DSCreate(PetscObjectComm((PetscObject)svd),&svd->ds);
463: PetscObjectIncrementTabLevel((PetscObject)svd->ds,(PetscObject)svd,0);
464: PetscObjectSetOptions((PetscObject)svd->ds,((PetscObject)svd)->options);
465: }
466: *ds = svd->ds;
467: return 0;
468: }