Actual source code: svdimpl.h

slepc-3.13.3 2020-06-14
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2020, 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: */

 11: #if !defined(SLEPCSVDIMPL_H)
 12: #define SLEPCSVDIMPL_H

 14:  #include <slepcsvd.h>
 15:  #include <slepc/private/slepcimpl.h>

 17: SLEPC_EXTERN PetscBool SVDRegisterAllCalled;
 18: SLEPC_EXTERN PetscErrorCode SVDRegisterAll(void);
 19: SLEPC_EXTERN PetscLogEvent SVD_SetUp,SVD_Solve;

 21: typedef struct _SVDOps *SVDOps;

 23: struct _SVDOps {
 24:   PetscErrorCode (*solve)(SVD);
 25:   PetscErrorCode (*setup)(SVD);
 26:   PetscErrorCode (*setfromoptions)(PetscOptionItems*,SVD);
 27:   PetscErrorCode (*publishoptions)(SVD);
 28:   PetscErrorCode (*destroy)(SVD);
 29:   PetscErrorCode (*reset)(SVD);
 30:   PetscErrorCode (*view)(SVD,PetscViewer);
 31: };

 33: /*
 34:      Maximum number of monitors you can run with a single SVD
 35: */
 36: #define MAXSVDMONITORS 5

 38: typedef enum { SVD_STATE_INITIAL,
 39:                SVD_STATE_SETUP,
 40:                SVD_STATE_SOLVED,
 41:                SVD_STATE_VECTORS } SVDStateType;

 43: /*
 44:    Defines the SVD data structure.
 45: */
 46: struct _p_SVD {
 47:   PETSCHEADER(struct _SVDOps);
 48:   /*------------------------- User parameters ---------------------------*/
 49:   Mat            OP;               /* problem matrix */
 50:   PetscInt       max_it;           /* max iterations */
 51:   PetscInt       nsv;              /* number of requested values */
 52:   PetscInt       ncv;              /* basis size */
 53:   PetscInt       mpd;              /* maximum dimension of projected problem */
 54:   PetscInt       nini,ninil;       /* number of initial vecs (negative means not copied yet) */
 55:   PetscReal      tol;              /* tolerance */
 56:   SVDConv        conv;             /* convergence test */
 57:   SVDStop        stop;             /* stopping test */
 58:   SVDWhich       which;            /* which singular values are computed */
 59:   PetscBool      impltrans;        /* implicit transpose mode */
 60:   PetscBool      trackall;         /* whether all the residuals must be computed */

 62:   /*-------------- User-provided functions and contexts -----------------*/
 63:   PetscErrorCode (*converged)(SVD,PetscReal,PetscReal,PetscReal*,void*);
 64:   PetscErrorCode (*convergeduser)(SVD,PetscReal,PetscReal,PetscReal*,void*);
 65:   PetscErrorCode (*convergeddestroy)(void*);
 66:   PetscErrorCode (*stopping)(SVD,PetscInt,PetscInt,PetscInt,PetscInt,SVDConvergedReason*,void*);
 67:   PetscErrorCode (*stoppinguser)(SVD,PetscInt,PetscInt,PetscInt,PetscInt,SVDConvergedReason*,void*);
 68:   PetscErrorCode (*stoppingdestroy)(void*);
 69:   void           *convergedctx;
 70:   void           *stoppingctx;
 71:   PetscErrorCode (*monitor[MAXSVDMONITORS])(SVD,PetscInt,PetscInt,PetscReal*,PetscReal*,PetscInt,void*);
 72:   PetscErrorCode (*monitordestroy[MAXSVDMONITORS])(void**);
 73:   void           *monitorcontext[MAXSVDMONITORS];
 74:   PetscInt       numbermonitors;

 76:   /*----------------- Child objects and working data -------------------*/
 77:   DS             ds;               /* direct solver object */
 78:   BV             U,V;              /* left and right singular vectors */
 79:   SlepcSC        sc;               /* sorting criterion data */
 80:   Mat            A;                /* problem matrix (m>n) */
 81:   Mat            AT;               /* transposed matrix */
 82:   Vec            *IS,*ISL;         /* placeholder for references to user initial space */
 83:   PetscReal      *sigma;           /* singular values */
 84:   PetscInt       *perm;            /* permutation for singular value ordering */
 85:   PetscReal      *errest;          /* error estimates */
 86:   void           *data;            /* placeholder for solver-specific stuff */

 88:   /* ----------------------- Status variables -------------------------- */
 89:   SVDStateType   state;            /* initial -> setup -> solved -> vectors */
 90:   PetscInt       nconv;            /* number of converged values */
 91:   PetscInt       its;              /* iteration counter */
 92:   PetscBool      leftbasis;        /* if U is filled by the solver */
 93:   SVDConvergedReason reason;
 94: };

 96: /*
 97:     Macros to test valid SVD arguments
 98: */
 99: #if !defined(PETSC_USE_DEBUG)

101: #define SVDCheckSolved(h,arg) do {} while (0)

103: #else

105: #define SVDCheckSolved(h,arg) \
106:   do { \
107:     if ((h)->state<SVD_STATE_SOLVED) SETERRQ1(PetscObjectComm((PetscObject)(h)),PETSC_ERR_ARG_WRONGSTATE,"Must call SVDSolve() first: Parameter #%d",arg); \
108:   } while (0)

110: #endif

112: PETSC_STATIC_INLINE PetscErrorCode SVDMatMult(SVD svd,PetscBool trans,Vec x,Vec y)
113: {

117:   if (trans) {
118:     if (svd->AT) {
119:       MatMult(svd->AT,x,y);
120:     } else {
121: #if defined(PETSC_USE_COMPLEX)
122:       MatMultHermitianTranspose(svd->A,x,y);
123: #else
124:       MatMultTranspose(svd->A,x,y);
125: #endif
126:     }
127:   } else {
128:     if (svd->A) {
129:       MatMult(svd->A,x,y);
130:     } else {
131: #if defined(PETSC_USE_COMPLEX)
132:       MatMultHermitianTranspose(svd->AT,x,y);
133: #else
134:       MatMultTranspose(svd->AT,x,y);
135: #endif
136:     }
137:   }
138:   return(0);
139: }

141: PETSC_STATIC_INLINE PetscErrorCode SVDMatCreateVecs(SVD svd,Vec *x,Vec *y)
142: {

146:   if (svd->A) {
147:     MatCreateVecs(svd->A,x,y);
148:   } else {
149:     MatCreateVecs(svd->AT,y,x);
150:   }
151:   return(0);
152: }

154: PETSC_STATIC_INLINE PetscErrorCode SVDMatCreateVecsEmpty(SVD svd,Vec *x,Vec *y)
155: {

159:   if (svd->A) {
160:     MatCreateVecsEmpty(svd->A,x,y);
161:   } else {
162:     MatCreateVecsEmpty(svd->AT,y,x);
163:   }
164:   return(0);
165: }

167: PETSC_STATIC_INLINE PetscErrorCode SVDMatGetSize(SVD svd,PetscInt *m,PetscInt *n)
168: {

172:   if (svd->A) {
173:     MatGetSize(svd->A,m,n);
174:   } else {
175:     MatGetSize(svd->AT,n,m);
176:   }
177:   return(0);
178: }

180: PETSC_STATIC_INLINE PetscErrorCode SVDMatGetLocalSize(SVD svd,PetscInt *m,PetscInt *n)
181: {

185:   if (svd->A) {
186:     MatGetLocalSize(svd->A,m,n);
187:   } else {
188:     MatGetLocalSize(svd->AT,n,m);
189:   }
190:   return(0);
191: }

193: SLEPC_INTERN PetscErrorCode SVDTwoSideLanczos(SVD,PetscReal*,PetscReal*,BV,BV,PetscInt,PetscInt);
194: SLEPC_INTERN PetscErrorCode SVDSetDimensions_Default(SVD);
195: SLEPC_INTERN PetscErrorCode SVDComputeVectors(SVD);

197: #endif