1 #ifndef VIRTUALMATRIX_HPP__
2 #define VIRTUALMATRIX_HPP__
3 #include "MatriceElementaire.hpp"
4
MATERROR(int i,const char * cmm)5 inline void MATERROR(int i,const char *cmm)
6 {
7 std::cerr << " MATERROR " << i << " : " << cmm << std::endl;
8 ErrorExec("MATERROR",1);
9 std::abort();
10 }
11 template<class Z,class R> class HashMatrix ;
12 template<class TypeIndex,class TypeScalar>
13 class VirtualMatrix: public RefCounter, public RNM_VirtualMatrix<TypeScalar,TypeIndex> {
14 public:
15 typedef TypeIndex I;
16 typedef TypeScalar R;
17 // 1 unsym , 2 herm, 4 sym, 8 pos , 16 nopos, 32 seq, 64 ompi, 128 mpi
18 static const int TS_unsym=1, TS_herm=2, TS_sym=4, TS_def_positif=8, TS_not_def_positif=16, TS_sequental = 32, TS_mpi = 64;// for verification
19 static const int TS_SYM=1,TS_DEF_POS=2,TS_PARA=4;
20 typedef VirtualMatrix<I,R> VMat ;
21 typedef void (*ERRORFunc)(int i,const char *cmm);
22 ERRORFunc ERRORHandle;
23 class VSolver {public:
24 virtual R* solve(R *x,R*b,int N=0,int trans=0) =0;
factorize(int st)25 virtual void factorize(int st){ if(verbosity>4) cout << " ** warning no factorisation" << st << endl;};
~VSolver()26 virtual ~VSolver(){}
27 int count;
VSolver()28 VSolver(): count(0) {}
copyptr()29 VSolver * copyptr(){ count++; return this;}
destroy()30 void destroy() { if(count==0) delete this; else count --;}
31 // virtual VSolver * clone() =0; to hard
32 };
33
ERROR(int i,const char * cmm) const34 void ERROR(int i,const char *cmm) const
35 {
36 if(ERRORHandle) (*ERRORHandle)(i,cmm);
37 else MATERROR(i,cmm);
38 }
39 I n,m; // size of matrix
40 bool symetric,positive_definite; // for cholesky or CG
41
42 VSolver *vsolver;
43 bool delvsolver;
44
45 // long state,codeini,codesym,codenum;
VirtualMatrix(I NN,I MM=-1,bool sym=false,bool dp=false)46 VirtualMatrix(I NN,I MM=-1,bool sym=false,bool dp=false) : RNM_VirtualMatrix<R,I> (NN,MM),
47 n(NN),m(this->M),symetric(sym),positive_definite(dp),vsolver(0),delvsolver(false)
48 {}
setsdp(int sym,bool dp)49 virtual void setsdp(int sym,bool dp) { symetric=(sym>0); positive_definite=dp;}
Set2Const(I n,R * x,R c=R ())50 static R *Set2Const(I n,R *x,R c=R()) { std::fill(x,x+n,c); return x;}
51
solve(R * x,R * b,int N=1,int transpo=0) const52 R* solve(R *x,R*b,int N=1,int transpo=0) const
53 { if(vsolver) return vsolver->solve(x,b,N,transpo);
54 else
55 ERROR(1,"VirtualMatrix:: no solver ?????");
56 return 0;
57 }
addMatMul(R * x,R * Ax,bool Transpose,I sx=1,I sAx=1) const58 virtual R* addMatMul(R *x,R*Ax,bool Transpose,I sx=1,I sAx=1) const { ERROR(1,"VirtualMatrix:: no AddMatMul ?????"); return 0;}
addMatMul(R * x,R * Ax) const59 R* addMatMul(R *x,R*Ax) const { return addMatMul(x,Ax,0); }
addMatTransMul(R * x,R * Atx) const60 R* addMatTransMul(R *x,R*Atx) const { return addMatMul(x,Atx,1); }
MatMul(R * x,R * Ax) const61 R* MatMul(R *x,R*Ax) const { return addMatMul(x, Set2Const(n,Ax)); }
MatTransMul(R * x,R * Atx) const62 R* MatTransMul(R *x,R*Atx) const { return addMatMul(x, Set2Const(m,Atx));}
WithSolver() const63 bool WithSolver() const {return vsolver;} // by default no solver
64 // void CloneSolver() { return vsolver ? vsolver->clone(): 0;}// for copy matrix to hard ...
SetSolver(const VirtualMatrix<I,R> & A)65 void SetSolver(const VirtualMatrix<I,R> & A)
66 {
67
68 if(A.vsolver) SetSolver(vsolver->copyptr(),A.delvsolver);
69 else SetSolver(0);
70 }
SetSolver(VSolver * f=0,bool del=true)71 void SetSolver(VSolver *f=0, bool del = true) {
72 if(verbosity>4) cout<< " ## SetSolver "<< this << " " << vsolver <<" "<< f << endl;
73 if(vsolver && delvsolver)
74 vsolver->destroy();
75 vsolver=f;
76 delvsolver=del;
77 // cout << "\n *** type SetSolver = " << typeid(f).name() << endl;
78 }
79
MatMul(KN_<R> & ax,const KN_<R> & x) const80 KN_<R> & MatMul(KN_<R> &ax,const KN_<R> &x) const {
81 MatMul((R*)x,(R*)ax);
82 return ax;}
addMatMul(const KN_<R> & x,KN_<R> & y) const83 void addMatMul(const KN_<R> & x, KN_<R> & y) const { ffassert(Checknm(y.N(),x.N())); addMatMul(x,y,0,x.step,y.step);}
addMatTransMul(const KN_<R> & x,KN_<R> & y) const84 void addMatTransMul(const KN_<R> & x , KN_<R> & y ) const {ffassert(Checknm(x.N(),y.N())); addMatMul(x,y,1,x.step,y.step);}
Solve(KN_<R> & x,const KN_<R> & b) const85 void Solve(KN_<R> & x,const KN_<R> & b) const {solve(x,b);}
SolveT(KN_<R> & x,const KN_<R> & b) const86 void SolveT(KN_<R> & x,const KN_<R> & b) const {solve(x,b,1,1);}
87
ChecknbLine(int nn) const88 bool ChecknbLine(int nn) const { return this->n==nn;}
ChecknbColumn(int mm) const89 bool ChecknbColumn(int mm) const { return this->m==mm;}
Checknm(int nn,int mm) const90 bool Checknm(int nn,int mm) const { return nn==this->n && mm==this->m;}
~VirtualMatrix()91 virtual ~VirtualMatrix(){
92 if(verbosity>99999) cout << " ** ~VirtualMatrix " << this << endl;
93 if(vsolver && delvsolver) vsolver->destroy();// Warning the solver is del afer the Matrix
94 } // clean solver
95
size() const96 virtual size_t size() const {return 0; };
operator +=(MatriceElementaire<R> &)97 virtual VirtualMatrix & operator +=(MatriceElementaire<R> & ){AFAIRE("VirtualMatrix::+=");}
operator =(const R & v)98 virtual void operator=(const R & v){AFAIRE("VirtualMatrix::=v");} // Mise a zero
dump(ostream &) const99 virtual ostream& dump (ostream&) const {cout << mpirank << " BUG virtualmatrix " << this << endl; AFAIRE("VirtualMatrix::dump");}
diag(I i)100 virtual R & diag(I i){AFAIRE("VirtualMatrix::diab");}
SetBC(I i,double tgv)101 virtual void SetBC(I i,double tgv){AFAIRE("VirtualMatrix::setbc");}
operator ()(I i,I j)102 virtual R & operator()(I i,I j){AFAIRE("VirtualMatrix::()(i,j)");}
pij(I i,I j) const103 virtual R * pij(I i,I j) const {AFAIRE("VirtualMatrix::pij");} // Add FH
104
toMatriceMorse(bool transpose=false,bool copy=false) const105 virtual HashMatrix<I, R> *toMatriceMorse(bool transpose=false,bool copy=false) const {return 0;} // not
addMatTo(R coef,HashMatrix<I,R> & mij,bool trans=false,int ii00=0,int jj00=0,bool cnj=false,double threshold=0.,const bool keepSym=false)106 virtual bool addMatTo(R coef,HashMatrix<I,R> &mij,bool trans=false,int ii00=0,int jj00=0,bool cnj=false,double threshold=0.,const bool keepSym=false)
107 {AFAIRE("VirtualMatrix::addMatTo");};
pscal(const KN_<R> & x,const KN_<R> & y)108 virtual R pscal(const KN_<R> & x,const KN_<R> & y) {AFAIRE("VirtualMatrix::pscal");} ; // produit scalaire
psor(KN_<R> & x,const KN_<R> & gmin,const KN_<R> & gmax,double omega)109 virtual double psor(KN_<R> & x,const KN_<R> & gmin,const KN_<R> & gmax , double omega) {AFAIRE("VirtualMatrix::psor");}
setdiag(const KN_<R> & x)110 virtual void setdiag(const KN_<R> & x){AFAIRE("VirtualMatrix::setdiag");} ;
getdiag(KN_<R> & x) const111 virtual void getdiag( KN_<R> & x) const {AFAIRE("VirtualMatrix::getdiag");}
NbCoef() const112 virtual I NbCoef() const {return 0;};
setcoef(const KN_<R> & x)113 virtual void setcoef(const KN_<R> & x){AFAIRE("VirtualMatrix::setcoef");}
getcoef(KN_<R> & x) const114 virtual void getcoef( KN_<R> & x) const {AFAIRE("VirtualMatrix::getcoef");}
sym() const115 virtual bool sym() const {return false;}
116
resize(I n,I m)117 virtual void resize(I n,I m) {AFAIRE("VirtualMatrix::resize");}
clear()118 virtual void clear()
119 { I NN=this->N, MM=this->M;
120 resize(0,0);
121 resize(NN,MM);
122 }
trace() const123 virtual R trace() const {ffassert(this->n==this->m); R t=R(), *p; for(int i=0; i<this->n; ++i) { p=pij(i,i); if(p) t+= *p;} return t; }
SetBC(char * wbc,double tgv)124 virtual void SetBC(char *wbc,double tgv) { for (int i=0; i<this->n; ++i) if(wbc[i]) SetBC(i,tgv);}
125 // void init(int nn=0,int mm=0) { VMat *p=new VMat(nn,mm); }
126
127 struct plusAx { const VirtualMatrix * A; const KN_<R> x;
plusAxVirtualMatrix::plusAx128 plusAx( const VirtualMatrix * B,const KN_<R> & y) :A(B),x(y)
129 { if(B) { ffassert(B->ChecknbColumn(y.N())); } }
callVirtualMatrix::plusAx130 void call(KN_<R> &ax,int init=0) const { if(init) ax=R(); A->addMatMul(x,ax,0,x.step,ax.step); }
131 };
132
133
134 struct plusAtx { const VirtualMatrix * A; const KN_<R> x;
plusAtxVirtualMatrix::plusAtx135 plusAtx( const VirtualMatrix * B,const KN_<R> & y) :A(B),x(y)
136 { if(B) { ffassert(B->ChecknbLine(y.N())); } }
callVirtualMatrix::plusAtx137 void call(KN_<R> &ax,int init=0) const {if(init) ax=R(); A->addMatMul(x,ax,1,x.step,ax.step); }
138 };
139
140 struct solveAxeqb { const VirtualMatrix * A; const KN_<R> b;
solveAxeqbVirtualMatrix::solveAxeqb141 solveAxeqb( const VirtualMatrix * B,const KN_<R> & y) :A(B),b(y)
142 { if(B) { ffassert(B->ChecknbColumn(y.N())); } }
callVirtualMatrix::solveAxeqb143 void call(KN_<R> &ax) const {ffassert(ax.contiguous() &&b.contiguous()); A->Solve(ax,b); }
144 };
145
146 struct solveAtxeqb { const VirtualMatrix * A; const KN_<R> b;
solveAtxeqbVirtualMatrix::solveAtxeqb147 solveAtxeqb( const VirtualMatrix * B,const KN_<R> & y) :A(B),b(y)
148 { if(B) { ffassert(B->ChecknbColumn(y.N())); } }
callVirtualMatrix::solveAtxeqb149 void call(KN_<R> &ax) const {ffassert(ax.contiguous() &&b.contiguous()); A->SolveT(ax,b); }
150 };
151
VirtualMatrix(const VirtualMatrix<I,R> & A)152 VirtualMatrix(const VirtualMatrix<I,R> &A)
153 : RNM_VirtualMatrix<TypeScalar,TypeIndex>(A) { operator=(A); }
operator =(const VirtualMatrix<I,R> & A)154 void operator=(const VirtualMatrix<I,R> &A)
155 {
156 n=A.n;
157 m=A.m; // size of matrix
158 symetric=A.symetric;
159 positive_definite=A.positive_definite; // for cholesky or CG
160 delvsolver = A.delvsolver;
161 vsolver = A.vsolver && delvsolver ? A.vsolver->copyptr() : 0;
162 }
163
164 };
165
166
167 template<class TypeIndex=int,class TypeScalar=double>
ProduitMatVec(const VirtualMatrix<TypeIndex,TypeScalar> * A,TypeScalar * x,TypeScalar * Ax)168 inline double * ProduitMatVec(const VirtualMatrix<TypeIndex,TypeScalar> *A,TypeScalar *x, TypeScalar *Ax) { return A->MatMul(x,Ax);}
169 template<class TypeIndex=int,class TypeScalar=double>
ProduitMatVec(const VirtualMatrix<TypeIndex,TypeScalar> & A,TypeScalar * x,TypeScalar * Ax)170 inline double * ProduitMatVec(const VirtualMatrix<TypeIndex,TypeScalar> &A,TypeScalar *x, TypeScalar *Ax) { return A.MatMul(x,Ax);}
171
172
173 #endif
174