1: #ifndef __LMVMMAT_H

  4:  #include <petscksp.h>

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

 12: #define MatLMVM_Rescale_None            0
 13: #define MatLMVM_Rescale_Scalar          1
 14: #define MatLMVM_Rescale_GL              2
 15: #define MatLMVM_Rescale_Types           3

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

 23: #define TAO_ZERO_SAFEGUARD      1e-8
 24: #define TAO_INF_SAFEGUARD       1e+8

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

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

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

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

 51:   PetscReal delta_max;  /*  Maximum value for delta */
 52:   PetscReal delta_min;  /*  Minimum value for delta */

 54:   PetscInt lmnow;
 55:   PetscInt iter;
 56:   PetscInt nupdates;
 57:   PetscInt nrejects;

 59:   Vec *S;
 60:   Vec *Y;
 61:   Vec Gprev;
 62:   Vec Xprev;

 64:   Vec D;
 65:   Vec U;
 66:   Vec V;
 67:   Vec W;
 68:   Vec P;
 69:   Vec Q;

 71:   PetscReal delta;
 72:   PetscReal sigma;
 73:   PetscReal theta;
 74:   PetscReal *rho;
 75:   PetscReal *beta;

 77:   PetscBool useDefaultH0;
 78:   Mat H0_mat;
 79:   KSP H0_ksp;
 80:   Vec H0_norm;

 82:   PetscBool useScale;
 83:   Vec scale;


 86: } MatLMVMCtx;


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


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

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

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


128: #endif