1: #ifndef __LMVMMAT_H


  5: #include <petsc/private/matimpl.h>
  6: #include <petscksp.h>

  8: #define MatLMVM_Scale_None              0
  9: #define MatLMVM_Scale_Scalar            1
 10: #define MatLMVM_Scale_Broyden           2
 11: #define MatLMVM_Scale_H0                3
 12: #define MatLMVM_Scale_Types             4

 14: #define MatLMVM_Rescale_None            0
 15: #define MatLMVM_Rescale_Scalar          1
 16: #define MatLMVM_Rescale_GL              2
 17: #define MatLMVM_Rescale_Types           3

 19: #define MatLMVM_Limit_None              0
 20: #define MatLMVM_Limit_Average           1
 21: #define MatLMVM_Limit_Relative          2
 22: #define MatLMVM_Limit_Absolute          3
 23: #define MatLMVM_Limit_Types             4

 25: #define TAO_ZERO_SAFEGUARD      1e-8
 26: #define TAO_INF_SAFEGUARD       1e+8

 28: typedef struct{
 29:     PetscBool allocated;
 30:     PetscInt lm;
 31:     PetscReal eps;
 32:     PetscInt limitType;
 33:     PetscInt scaleType;
 34:     PetscInt rScaleType;

 36:     PetscReal s_alpha;  /*  Factor for scalar scaling */
 37:     PetscReal r_alpha;  /*  Factor on scalar for rescaling diagonal matrix */
 38:     PetscReal r_beta;   /*  Factor on diagonal for rescaling diagonal matrix */
 39:     PetscReal mu;               /*  Factor for using historical information */
 40:     PetscReal nu;               /*  Factor for using historical information */
 41:     PetscReal phi;              /*  Factor for Broyden scaling */

 43:   PetscInt scalar_history;      /*  Amount of history to keep for scalar scaling */
 44:   PetscReal *yy_history;        /*  Past information for scalar scaling */
 45:   PetscReal *ys_history;        /*  Past information for scalar scaling */
 46:   PetscReal *ss_history;        /*  Past information for scalar scaling */

 48:   PetscInt rescale_history;  /*  Amount of history to keep for rescaling diagonal */
 49:   PetscReal *yy_rhistory;       /*  Past information for scalar rescaling */
 50:   PetscReal *ys_rhistory;       /*  Past information for scalar rescaling */
 51:   PetscReal *ss_rhistory;       /*  Past information for scalar rescaling */

 53:   PetscReal delta_max;  /*  Maximum value for delta */
 54:   PetscReal delta_min;  /*  Minimum value for delta */

 56:   PetscInt lmnow;
 57:   PetscInt iter;
 58:   PetscInt nupdates;
 59:   PetscInt nrejects;

 61:   Vec *S;
 62:   Vec *Y;
 63:   Vec Gprev;
 64:   Vec Xprev;

 66:   Vec D;
 67:   Vec U;
 68:   Vec V;
 69:   Vec W;
 70:   Vec P;
 71:   Vec Q;

 73:   PetscReal delta;
 74:   PetscReal sigma;
 75:   PetscReal theta;
 76:   PetscReal *rho;
 77:   PetscReal *beta;

 79:   PetscBool useDefaultH0;
 80:   Mat H0_mat;
 81:   KSP H0_ksp;
 82:   Vec H0_norm;

 84:   PetscBool useScale;
 85:   Vec scale;


 88: } MatLMVMCtx;


 91: extern PetscErrorCode MatCreateLMVM(MPI_Comm,PetscInt,PetscInt,Mat*);


 94: /* PETSc Mat overrides */
 95: extern PetscErrorCode MatView_LMVM(Mat,PetscViewer);
 96: extern PetscErrorCode MatDestroy_LMVM(Mat);

 98: /*
 99: int MatMultTranspose_LMVM(Mat,Vec,Vec);
100: int MatDiagonalShift_LMVM(Vec,Mat);
101: int MatDestroy_LMVM(Mat);
102: int MatShift_LMVM(Mat,PetscReal);
103: int MatDuplicate_LMVM(Mat,MatDuplicateOption,Mat*);
104: int MatEqual_LMVM(Mat,Mat,PetscBool*);
105: int MatScale_LMVM(Mat,PetscReal);
106: int MatGetSubMatrix_LMVM(Mat,IS,IS,int,MatReuse,Mat *);
107: int MatGetSubMatrices_LMVM(Mat,int,IS*,IS*,MatReuse,Mat**);
108: int MatTranspose_LMVM(Mat,Mat*);
109: int MatGetDiagonal_LMVM(Mat,Vec);
110: int MatGetColumnVector_LMVM(Mat,Vec, int);
111: int MatNorm_LMVM(Mat,NormType,PetscReal *);
112: */

114: /* Functions used by TAO */
115: PetscErrorCode MatLMVMReset(Mat);
116: PetscErrorCode MatLMVMUpdate(Mat,Vec, Vec);
117: PetscErrorCode MatLMVMSetDelta(Mat,PetscReal);
118: PetscErrorCode MatLMVMSetScale(Mat,Vec);
119: PetscErrorCode MatLMVMGetRejects(Mat,PetscInt*);
120: PetscErrorCode MatLMVMSetH0(Mat,Mat);
121: PetscErrorCode MatLMVMGetH0(Mat,Mat*);
122: PetscErrorCode MatLMVMGetH0KSP(Mat,KSP*);
123: PetscErrorCode MatLMVMSetPrev(Mat,Vec,Vec);
124: PetscErrorCode MatLMVMGetX0(Mat,Vec);
125: PetscErrorCode MatLMVMRefine(Mat, Mat, Vec, Vec);
126: PetscErrorCode MatLMVMAllocateVectors(Mat m, Vec v);
127: PetscErrorCode MatLMVMSolve(Mat, Vec, Vec);


130: #endif