Actual source code: svdbasic.c
slepc-3.22.1 2024-10-28
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 = 1;
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->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: PetscCall(PetscNew(&svd->sc));
115: *outsvd = svd;
116: PetscFunctionReturn(PETSC_SUCCESS);
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
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: {
134: PetscFunctionBegin;
136: if (!svd) PetscFunctionReturn(PETSC_SUCCESS);
137: PetscTryTypeMethod(svd,reset);
138: PetscCall(MatDestroy(&svd->OP));
139: PetscCall(MatDestroy(&svd->OPb));
140: PetscCall(VecDestroy(&svd->omega));
141: PetscCall(MatDestroy(&svd->A));
142: PetscCall(MatDestroy(&svd->B));
143: PetscCall(MatDestroy(&svd->AT));
144: PetscCall(MatDestroy(&svd->BT));
145: PetscCall(BVDestroy(&svd->U));
146: PetscCall(BVDestroy(&svd->V));
147: PetscCall(VecDestroyVecs(svd->nworkl,&svd->workl));
148: svd->nworkl = 0;
149: PetscCall(VecDestroyVecs(svd->nworkr,&svd->workr));
150: svd->nworkr = 0;
151: svd->swapped = PETSC_FALSE;
152: svd->state = SVD_STATE_INITIAL;
153: PetscFunctionReturn(PETSC_SUCCESS);
154: }
156: /*@
157: SVDDestroy - Destroys the SVD context.
159: Collective
161: Input Parameter:
162: . svd - singular value solver context obtained from SVDCreate()
164: Level: beginner
166: .seealso: SVDCreate(), SVDSetUp(), SVDSolve()
167: @*/
168: PetscErrorCode SVDDestroy(SVD *svd)
169: {
170: PetscFunctionBegin;
171: if (!*svd) PetscFunctionReturn(PETSC_SUCCESS);
173: if (--((PetscObject)*svd)->refct > 0) { *svd = NULL; PetscFunctionReturn(PETSC_SUCCESS); }
174: PetscCall(SVDReset(*svd));
175: PetscTryTypeMethod(*svd,destroy);
176: if ((*svd)->sigma) PetscCall(PetscFree3((*svd)->sigma,(*svd)->perm,(*svd)->errest));
177: if ((*svd)->sign) PetscCall(PetscFree((*svd)->sign));
178: PetscCall(DSDestroy(&(*svd)->ds));
179: PetscCall(PetscFree((*svd)->sc));
180: /* just in case the initial vectors have not been used */
181: PetscCall(SlepcBasisDestroy_Private(&(*svd)->nini,&(*svd)->IS));
182: PetscCall(SlepcBasisDestroy_Private(&(*svd)->ninil,&(*svd)->ISL));
183: PetscCall(SVDMonitorCancel(*svd));
184: PetscCall(PetscHeaderDestroy(svd));
185: PetscFunctionReturn(PETSC_SUCCESS);
186: }
188: /*@
189: SVDSetType - Selects the particular solver to be used in the SVD object.
191: Logically Collective
193: Input Parameters:
194: + svd - the singular value solver context
195: - type - a known method
197: Options Database Key:
198: . -svd_type <method> - Sets the method; use -help for a list
199: of available methods
201: Notes:
202: See "slepc/include/slepcsvd.h" for available methods. The default
203: is SVDCROSS.
205: Normally, it is best to use the SVDSetFromOptions() command and
206: then set the SVD type from the options database rather than by using
207: this routine. Using the options database provides the user with
208: maximum flexibility in evaluating the different available methods.
209: The SVDSetType() routine is provided for those situations where it
210: is necessary to set the iterative solver independently of the command
211: line or options database.
213: Level: intermediate
215: .seealso: SVDType
216: @*/
217: PetscErrorCode SVDSetType(SVD svd,SVDType type)
218: {
219: PetscErrorCode (*r)(SVD);
220: PetscBool match;
222: PetscFunctionBegin;
224: PetscAssertPointer(type,2);
226: PetscCall(PetscObjectTypeCompare((PetscObject)svd,type,&match));
227: if (match) PetscFunctionReturn(PETSC_SUCCESS);
229: PetscCall(PetscFunctionListFind(SVDList,type,&r));
230: PetscCheck(r,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown SVD type given: %s",type);
232: PetscTryTypeMethod(svd,destroy);
233: PetscCall(PetscMemzero(svd->ops,sizeof(struct _SVDOps)));
235: svd->state = SVD_STATE_INITIAL;
236: PetscCall(PetscObjectChangeTypeName((PetscObject)svd,type));
237: PetscCall((*r)(svd));
238: PetscFunctionReturn(PETSC_SUCCESS);
239: }
241: /*@
242: SVDGetType - Gets the SVD type as a string from the SVD object.
244: Not Collective
246: Input Parameter:
247: . svd - the singular value solver context
249: Output Parameter:
250: . type - name of SVD method
252: Level: intermediate
254: .seealso: SVDSetType()
255: @*/
256: PetscErrorCode SVDGetType(SVD svd,SVDType *type)
257: {
258: PetscFunctionBegin;
260: PetscAssertPointer(type,2);
261: *type = ((PetscObject)svd)->type_name;
262: PetscFunctionReturn(PETSC_SUCCESS);
263: }
265: /*@C
266: SVDRegister - Adds a method to the singular value solver package.
268: Not Collective
270: Input Parameters:
271: + name - name of a new user-defined solver
272: - function - routine to create the solver context
274: Notes:
275: SVDRegister() may be called multiple times to add several user-defined solvers.
277: Example Usage:
278: .vb
279: SVDRegister("my_solver",MySolverCreate);
280: .ve
282: Then, your solver can be chosen with the procedural interface via
283: $ SVDSetType(svd,"my_solver")
284: or at runtime via the option
285: $ -svd_type my_solver
287: Level: advanced
289: .seealso: SVDRegisterAll()
290: @*/
291: PetscErrorCode SVDRegister(const char *name,PetscErrorCode (*function)(SVD))
292: {
293: PetscFunctionBegin;
294: PetscCall(SVDInitializePackage());
295: PetscCall(PetscFunctionListAdd(&SVDList,name,function));
296: PetscFunctionReturn(PETSC_SUCCESS);
297: }
299: /*@C
300: SVDMonitorRegister - Adds SVD monitor routine.
302: Not Collective
304: Input Parameters:
305: + name - name of a new monitor routine
306: . vtype - a PetscViewerType for the output
307: . format - a PetscViewerFormat for the output
308: . monitor - monitor routine
309: . create - creation routine, or NULL
310: - destroy - destruction routine, or NULL
312: Notes:
313: SVDMonitorRegister() may be called multiple times to add several user-defined monitors.
315: Example Usage:
316: .vb
317: SVDMonitorRegister("my_monitor",PETSCVIEWERASCII,PETSC_VIEWER_ASCII_INFO_DETAIL,MyMonitor,NULL,NULL);
318: .ve
320: Then, your monitor can be chosen with the procedural interface via
321: $ SVDMonitorSetFromOptions(svd,"-svd_monitor_my_monitor","my_monitor",NULL)
322: or at runtime via the option
323: $ -svd_monitor_my_monitor
325: Level: advanced
327: .seealso: SVDMonitorRegisterAll()
328: @*/
329: 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**))
330: {
331: char key[PETSC_MAX_PATH_LEN];
333: PetscFunctionBegin;
334: PetscCall(SVDInitializePackage());
335: PetscCall(SlepcMonitorMakeKey_Internal(name,vtype,format,key));
336: PetscCall(PetscFunctionListAdd(&SVDMonitorList,key,monitor));
337: if (create) PetscCall(PetscFunctionListAdd(&SVDMonitorCreateList,key,create));
338: if (destroy) PetscCall(PetscFunctionListAdd(&SVDMonitorDestroyList,key,destroy));
339: PetscFunctionReturn(PETSC_SUCCESS);
340: }
342: /*@
343: SVDSetBV - Associates basis vectors objects to the singular value solver.
345: Collective
347: Input Parameters:
348: + svd - singular value solver context obtained from SVDCreate()
349: . V - the basis vectors object for right singular vectors
350: - U - the basis vectors object for left singular vectors
352: Note:
353: Use SVDGetBV() to retrieve the basis vectors contexts (for example,
354: to free them at the end of the computations).
356: Level: advanced
358: .seealso: SVDGetBV()
359: @*/
360: PetscErrorCode SVDSetBV(SVD svd,BV V,BV U)
361: {
362: PetscFunctionBegin;
364: if (V) {
366: PetscCheckSameComm(svd,1,V,2);
367: PetscCall(PetscObjectReference((PetscObject)V));
368: PetscCall(BVDestroy(&svd->V));
369: svd->V = V;
370: }
371: if (U) {
373: PetscCheckSameComm(svd,1,U,3);
374: PetscCall(PetscObjectReference((PetscObject)U));
375: PetscCall(BVDestroy(&svd->U));
376: svd->U = U;
377: }
378: PetscFunctionReturn(PETSC_SUCCESS);
379: }
381: /*@
382: SVDGetBV - Obtain the basis vectors objects associated to the singular
383: value solver object.
385: Not Collective
387: Input Parameter:
388: . svd - singular value solver context obtained from SVDCreate()
390: Output Parameters:
391: + V - basis vectors context for right singular vectors
392: - U - basis vectors context for left singular vectors
394: Level: advanced
396: .seealso: SVDSetBV()
397: @*/
398: PetscErrorCode SVDGetBV(SVD svd,BV *V,BV *U)
399: {
400: PetscFunctionBegin;
402: if (V) {
403: if (!svd->V) {
404: PetscCall(BVCreate(PetscObjectComm((PetscObject)svd),&svd->V));
405: PetscCall(PetscObjectIncrementTabLevel((PetscObject)svd->V,(PetscObject)svd,0));
406: PetscCall(PetscObjectSetOptions((PetscObject)svd->V,((PetscObject)svd)->options));
407: }
408: *V = svd->V;
409: }
410: if (U) {
411: if (!svd->U) {
412: PetscCall(BVCreate(PetscObjectComm((PetscObject)svd),&svd->U));
413: PetscCall(PetscObjectIncrementTabLevel((PetscObject)svd->U,(PetscObject)svd,0));
414: PetscCall(PetscObjectSetOptions((PetscObject)svd->U,((PetscObject)svd)->options));
415: }
416: *U = svd->U;
417: }
418: PetscFunctionReturn(PETSC_SUCCESS);
419: }
421: /*@
422: SVDSetDS - Associates a direct solver object to the singular value solver.
424: Collective
426: Input Parameters:
427: + svd - singular value solver context obtained from SVDCreate()
428: - ds - the direct solver object
430: Note:
431: Use SVDGetDS() to retrieve the direct solver context (for example,
432: to free it at the end of the computations).
434: Level: advanced
436: .seealso: SVDGetDS()
437: @*/
438: PetscErrorCode SVDSetDS(SVD svd,DS ds)
439: {
440: PetscFunctionBegin;
443: PetscCheckSameComm(svd,1,ds,2);
444: PetscCall(PetscObjectReference((PetscObject)ds));
445: PetscCall(DSDestroy(&svd->ds));
446: svd->ds = ds;
447: PetscFunctionReturn(PETSC_SUCCESS);
448: }
450: /*@
451: SVDGetDS - Obtain the direct solver object associated to the singular value
452: solver object.
454: Not Collective
456: Input Parameters:
457: . svd - singular value solver context obtained from SVDCreate()
459: Output Parameter:
460: . ds - direct solver context
462: Level: advanced
464: .seealso: SVDSetDS()
465: @*/
466: PetscErrorCode SVDGetDS(SVD svd,DS *ds)
467: {
468: PetscFunctionBegin;
470: PetscAssertPointer(ds,2);
471: if (!svd->ds) {
472: PetscCall(DSCreate(PetscObjectComm((PetscObject)svd),&svd->ds));
473: PetscCall(PetscObjectIncrementTabLevel((PetscObject)svd->ds,(PetscObject)svd,0));
474: PetscCall(PetscObjectSetOptions((PetscObject)svd->ds,((PetscObject)svd)->options));
475: }
476: *ds = svd->ds;
477: PetscFunctionReturn(PETSC_SUCCESS);
478: }