Actual source code: svdimpl.h

slepc-3.6.3 2016-03-29
Report Typos and Errors
  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