Actual source code: svdimpl.h
slepc-3.6.3 2016-03-29
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2015, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
8: SLEPc is free software: you can redistribute it and/or modify it under the
9: terms of version 3 of the GNU Lesser General Public License as published by
10: the Free Software Foundation.
12: SLEPc is distributed in the hope that it will be useful, but WITHOUT ANY
13: WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14: FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
15: more details.
17: You should have received a copy of the GNU Lesser General Public License
18: along with SLEPc. If not, see <http://www.gnu.org/licenses/>.
19: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
20: */
22: #if !defined(_SVDIMPL)
23: #define _SVDIMPL
25: #include <slepcsvd.h>
26: #include <slepc/private/slepcimpl.h>
28: PETSC_EXTERN PetscBool SVDRegisterAllCalled;
29: PETSC_EXTERN PetscErrorCode SVDRegisterAll(void);
30: PETSC_EXTERN PetscLogEvent SVD_SetUp,SVD_Solve;
32: typedef struct _SVDOps *SVDOps;
34: struct _SVDOps {
35: PetscErrorCode (*solve)(SVD);
36: PetscErrorCode (*setup)(SVD);
37: PetscErrorCode (*setfromoptions)(PetscOptions*,SVD);
38: PetscErrorCode (*publishoptions)(SVD);
39: PetscErrorCode (*destroy)(SVD);
40: PetscErrorCode (*reset)(SVD);
41: PetscErrorCode (*view)(SVD,PetscViewer);
42: };
44: /*
45: Maximum number of monitors you can run with a single SVD
46: */
47: #define MAXSVDMONITORS 5
49: typedef enum { SVD_STATE_INITIAL,
50: SVD_STATE_SETUP,
51: SVD_STATE_SOLVED,
52: SVD_STATE_VECTORS } SVDStateType;
54: /*
55: Defines the SVD data structure.
56: */
57: struct _p_SVD {
58: PETSCHEADER(struct _SVDOps);
59: /*------------------------- User parameters ---------------------------*/
60: Mat OP; /* problem matrix */
61: PetscInt max_it; /* max iterations */
62: PetscInt nsv; /* number of requested values */
63: PetscInt ncv; /* basis size */
64: PetscInt mpd; /* maximum dimension of projected problem */
65: PetscInt nini,ninil; /* number of initial vectors (negative means not copied yet) */
66: PetscReal tol; /* tolerance */
67: SVDWhich which; /* which singular values are computed */
68: PetscBool impltrans; /* implicit transpose mode */
69: PetscBool trackall; /* whether all the residuals must be computed */
71: /*-------------- User-provided functions and contexts -----------------*/
72: PetscErrorCode (*monitor[MAXSVDMONITORS])(SVD,PetscInt,PetscInt,PetscReal*,PetscReal*,PetscInt,void*);
73: PetscErrorCode (*monitordestroy[MAXSVDMONITORS])(void**);
74: void *monitorcontext[MAXSVDMONITORS];
75: PetscInt numbermonitors;
77: /*----------------- Child objects and working data -------------------*/
78: DS ds; /* direct solver object */
79: BV U,V; /* left and right singular vectors */
80: PetscRandom rand; /* random number generator */
81: SlepcSC sc; /* sorting criterion data */
82: Mat A; /* problem matrix (m>n) */
83: Mat AT; /* transposed matrix */
84: Vec *IS,*ISL; /* placeholder for references to user-provided initial space */
85: PetscReal *sigma; /* singular values */
86: PetscInt *perm; /* permutation for singular value ordering */
87: PetscReal *errest; /* error estimates */
88: void *data; /* placeholder for solver-specific stuff */
90: /* ----------------------- Status variables -------------------------- */
91: SVDStateType state; /* initial -> setup -> solved -> vectors */
92: PetscInt nconv; /* number of converged values */
93: PetscInt its; /* iteration counter */
94: PetscBool leftbasis; /* if U is filled by the solver */
95: SVDConvergedReason reason;
96: };
98: /*
99: Macros to test valid SVD arguments
100: */
101: #if !defined(PETSC_USE_DEBUG)
103: #define SVDCheckSolved(h,arg) do {} while (0)
105: #else
107: #define SVDCheckSolved(h,arg) \
108: do { \
109: if (h->state<SVD_STATE_SOLVED) SETERRQ1(PetscObjectComm((PetscObject)h),PETSC_ERR_ARG_WRONGSTATE,"Must call SVDSolve() first: Parameter #%d",arg); \
110: } while (0)
112: #endif
116: PETSC_STATIC_INLINE PetscErrorCode SVDMatMult(SVD svd,PetscBool trans,Vec x,Vec y)
117: {
121: if (trans) {
122: if (svd->AT) {
123: MatMult(svd->AT,x,y);
124: } else {
125: #if defined(PETSC_USE_COMPLEX)
126: MatMultHermitianTranspose(svd->A,x,y);
127: #else
128: MatMultTranspose(svd->A,x,y);
129: #endif
130: }
131: } else {
132: if (svd->A) {
133: MatMult(svd->A,x,y);
134: } else {
135: #if defined(PETSC_USE_COMPLEX)
136: MatMultHermitianTranspose(svd->AT,x,y);
137: #else
138: MatMultTranspose(svd->AT,x,y);
139: #endif
140: }
141: }
142: return(0);
143: }
147: PETSC_STATIC_INLINE PetscErrorCode SVDMatCreateVecs(SVD svd,Vec *x,Vec *y)
148: {
152: if (svd->A) {
153: MatCreateVecs(svd->A,x,y);
154: } else {
155: MatCreateVecs(svd->AT,y,x);
156: }
157: return(0);
158: }
162: PETSC_STATIC_INLINE PetscErrorCode SVDMatGetSize(SVD svd,PetscInt *m,PetscInt *n)
163: {
167: if (svd->A) {
168: MatGetSize(svd->A,m,n);
169: } else {
170: MatGetSize(svd->AT,n,m);
171: }
172: return(0);
173: }
177: PETSC_STATIC_INLINE PetscErrorCode SVDMatGetLocalSize(SVD svd,PetscInt *m,PetscInt *n)
178: {
182: if (svd->A) {
183: MatGetLocalSize(svd->A,m,n);
184: } else {
185: MatGetLocalSize(svd->AT,n,m);
186: }
187: return(0);
188: }
190: PETSC_INTERN PetscErrorCode SVDTwoSideLanczos(SVD,PetscReal*,PetscReal*,BV,BV,PetscInt,PetscInt);
191: PETSC_INTERN PetscErrorCode SVDSetDimensions_Default(SVD);
192: PETSC_INTERN PetscErrorCode SVDComputeVectors(SVD);
194: #endif