1 #ifndef MatriceCreuse_h_
2 #define MatriceCreuse_h_
3
4 #define VERSION_MATRICE_CREUSE 4
5
6 // add Sep 2007 for generic Space solver
7 #define DCL_ARG_SPARSE_SOLVER(T,A) Stack stack, HashMatrix<int,T> *A, Data_Sparse_Solver & ds
8 #define ARG_SPARSE_SOLVER(A) stack,A, ds
9
10 #ifdef HAVE_LIBUMFPACK_XXXXXXXXXXXXX
11 extern "C" {
12 #ifdef HAVE_UMFPACK_H
13 #include <umfpack.h>
14 #else
15 #ifdef HAVE_UMFPACK_UMFPACK_H
16 #include <umfpack/umfpack.h>
17 #else
18 #ifdef HAVE_BIG_UMFPACK_UMFPACK_H
19 #include <UMFPACK/umfpack.h>
20 #else
21 #ifdef HAVE_UFSPARSE_UMFPACK_H
22 #include <ufsparse/umfpack.h>
23 #else
24 #ifdef HAVE_SUITESPARSE_UMFPACK_H
25 #include <suitesparse/umfpack.h>
26 #else
27
28 // Defaults to a local version of the UMFPACK headers
29 #include "../../3rdparty/include/umfpack.h"
30
31 #endif // HAVE_SUITESPARSE_UMFPACK_H
32 #endif // HAVE_UFSPARSE_UMFPACK_H
33 #endif // HAVE_BIG_UMFPACK_UMFPACK_H
34 #endif // HAVE_UMFPACK_UMFPACK_H
35 #endif // HAVE_UMFPACK_H
36 }
37 #endif
38
39 #include "RNM.hpp"
40 #include "fem.hpp"
41 #include "FESpace.hpp"
42 #include "DOperator.hpp"
43 #include "QuadratureFormular.hpp"
44 #include "HashMatrix.hpp"
45
46 extern double ff_tgv;
47 using Fem2D::Mesh;
48 using Fem2D::FESpace;
49 using Fem2D::FElement;
50 using Fem2D::baseFElement;
51 using Fem2D::FMortar;
52 using Fem2D::TypeOfMortar;
53 using Fem2D::QuadratureFormular;
54 using Fem2D::QuadratureFormular1d;
55 using Fem2D::QuadratureFormular_T_5;
56 using Fem2D::QF_GaussLegendre3;
57 const double EPSILON=1e-20;
58 using Fem2D::onWhatIsEdge;
59
60 //#define APROGRAMMER(a) {cerr << "A PROGRAMMER " #a << endl; exit (1) ;}
61 #define ERREUR(a,b) {cerr << "ERROR " #a<< b <<endl; throw(ErrorExec("FATAL ERROR in " __FILE__ "\n" #a " line: ",__LINE__)) ;}
62
63 template<class TypeIndex=int,class TypeScalar=double> class VirtualMatrix;
64 template<class TypeIndex,class TypeScalaire> class HashMatrix ;
65 template<class K> using MatriceCreuse=VirtualMatrix<int, K> ;
66 template<class K> using MatriceMorse=HashMatrix<int,K>;
67
68 //template <class R> class MatriceCreuse;
69 template <class R> class MatriceElementaire;
70 template <class R,class FESpace> class MatriceElementaireSymetrique;
71 template <class R,class FESpace> class MatriceElementairePleine;
72 //template <class R> class MatriceMorseOld;
73 template <class R> class MatriceProdTensoriel;
74
75 //template <class R> R Square(R x){ return x*x;}
76
docpyornot(bool nocpy,T * p,int n)77 template <class T> T* docpyornot(bool nocpy,T* p,int n)
78 {
79 T * r=p;
80 if( !nocpy) { // do copy
81 r= new T[n]; ffassert(r);
82 if(p)
83 for(int i=0;i<n;i++)
84 r[i]=p[i];
85 }
86 else if( r==0) // always do allocation July 2009 for mpi matrice
87 { r= new T[n]; ffassert(r);}
88 return r;
89 }
90 template <class T,class TT> T* docpy(TT* p,int n)
91 {
92 T * r=0;
93 if(p && n) { // do copy
94 r= new T[n]; ffassert(r);
95 for(int i=0;i<n;i++)
96 r[i]=T(p[i]); // pour cadna ???? FH
97 }
98 return r;
99 }
100
101
102
103
104 template <typename Z,typename R> class VirtualMatrix;
105 template <typename Z,typename R> class HashMatrix;
106
107
108
Stop(int iter,R *,R *)109 template<class R> class StopGC { public: virtual bool Stop(int iter, R *, R * ){cout << " Stop !!!!!\n"; return false;} };
110 template<class R,class M,class P,class S >// S=StopGC<Real>
ConjuguedGradient(const M & A,const P & C,const KN_<R> & b,KN_<R> & x,const int nbitermax,double & eps,long kprint=1000000000,S * Stop=0)111 int ConjuguedGradient(const M & A,const P & C,const KN_<R> &b,KN_<R> &x,const int nbitermax, double &eps,long kprint=1000000000,S *Stop=0)
112 {
113
114 // ConjuguedGradient lineare A*x est appele avec des conditions au limites
115 // non-homogene puis homogene pour calculer le gradient
116 if (verbosity>50)
117 kprint=2;
118 if (verbosity>99) cout << A << endl;
119 // throwassert(&x && &b && &A && &C);
120 typedef KN<R> Rn;
121 int n=b.N();
122 throwassert(n==x.N());
123 Rn g(n), h(n), Ah(n), & Cg(Ah); // on utilise Ah pour stocke Cg
124 g = A*x;
125 double xx= std::real((x,conj(x)));
126 double epsold=eps;
127 g -= b;// g = Ax-b
128 Cg = C*g; // gradient preconditionne
129 h =-Cg;
130 double g2 = std::real((Cg,conj(g)));
131 if (g2 < 1e-30)
132 { if(verbosity>1 || (kprint<100000))
133 cout << "GC g^2 =" << g2 << " < 1.e-30 Nothing to do " << endl;
134 return 2; }
135 double reps2 =eps >0 ? eps*eps*g2 : -eps; // epsilon relatif
136 eps = reps2;
137 for (int iter=0;iter<=nbitermax;iter++)
138 {
139 Ah = A*h;
140 double hAh =std::real((h,conj(Ah)));
141 // if (Abs(hAh)<1e-30) ExecError("CG2: Matrix non defined, sorry ");
142 R ro = - std::real((g,conj(h)))/ hAh; // ro optimal (produit scalaire usuel)
143 x += ro *h;
144 g += ro *Ah; // plus besoin de Ah, on utilise avec Cg optimisation
145 Cg = C*g;
146 double g2p=g2;
147 g2 = std::real((Cg,conj(g)));
148 bool stop = Stop && Stop->Stop(iter,x,g);
149
150 if ( !(iter%kprint) && iter && (verbosity>3) )
151 cout << "CG:" <<iter << " ro = " << ro << " ||g||^2 = " << g2 << " " << stop << endl;
152
153 if (g2 < reps2 || stop) {
154 if ( !iter && !xx && g2 && epsold >0 ) {
155 // change fo eps converge to fast due to the
156 // penalization of boundary condition.
157 eps = epsold*epsold*g2;
158 if (verbosity>3 || (kprint<3))
159 cout << "CG converge to fast (pb of BC) restart: " << iter << " ro = "
160 << ro << " ||g||^2 = " << g2 << " <= " << reps2 << " new eps2 =" << eps <<endl;
161 reps2=eps;
162 }
163 else
164 {
165 if (verbosity>1 || (kprint<100000) )
166 cout << "CG converge: " << iter << " ro = " << ro << " ||g||^2 = " << g2 << endl;
167 return 1;// ok
168 }
169 }
170 double gamma = g2/g2p;
171 h *= gamma;
172 h -= Cg; // h = -Cg * gamma* h
173 }
174 if(verbosity)
175 cout << " GC: method doesn't converge in " << nbitermax
176 << " iteration , xx= " << xx<< endl;
177 return 0;
178 }
179
180 template<class R,class M,class P,class S>// S=StopGC<Real>
ConjuguedGradient2(const M & A,const P & C,KN_<R> & x,const KN_<R> & b,const int nbitermax,double & eps,long kprint=1000000000,S * Stop=0)181 int ConjuguedGradient2(const M & A,const P & C,KN_<R> &x,const KN_<R> &b,const int nbitermax, double &eps,long kprint=1000000000,S *Stop=0)
182 {
183 // ConjuguedGradient2 affine A*x = 0 est toujours appele avec les condition aux limites
184 // -------------
185
186 // throwassert(&x && &A && &C);
187 typedef KN<R> Rn;
188 int n=x.N();
189 // if (verbosity>99) kprint=1;
190 R ro=1;
191 Rn g(n),h(n),Ah(n), & Cg(Ah); // on utilise Ah pour stocke Cg
192 g = A*x;
193 g -= b;
194 Cg = C*g; // gradient preconditionne
195 h =-Cg;
196 R g2 = (Cg,g);
197 if (g2 < 1e-30)
198 { if(verbosity>1 || kprint< 1000000)
199 cout << "GC g^2 =" << g2 << " < 1.e-30 Nothing to do " << endl;
200 return 2; }
201 if (verbosity>5 || (kprint<2))
202 cout << " 0 GC g^2 =" << g2 << endl;
203 R reps2 =eps >0 ? eps*eps*g2 : -eps; // epsilon relatif
204 eps = reps2;
205 for (int iter=0;iter<=nbitermax;iter++)
206 {
207 R rop = ro;
208 x += rop*h; // x+ rop*h , g=Ax (x old)
209 // ((Ah = A*x - b) - g);
210 // Ah -= b; // Ax + rop*Ah = rop*Ah + g =
211 // Ah -= g; // Ah*rop
212 Ah = A*x;
213 Ah -= b; // Ax + rop*Ah = rop*Ah + g =
214 Ah -= g; // Ah*rop
215 R hAh =(h,Ah);
216 if (std::norm(hAh)<1e-60) ExecError("CG2: Matrix is not defined (/0), sorry ");
217 ro = - (g,h)*rop/hAh ; // ro optimal (produit scalaire usuel)
218 if ( ro != ro) ExecError("CG2: Bug : ro is NaN ??? ");
219 x += (ro-rop) *h;
220 g += (ro/rop) *Ah; // plus besoin de Ah, on utilise avec Cg optimisation
221 Cg = C*g;
222 R g2p=g2;
223 g2 = (Cg,g);
224 bool stop = Stop && Stop->Stop(iter,x,g);
225 if ( ( (iter%kprint) == kprint-1) /*&& verbosity >1*/ )
226 cout << "CG:" <<iter << " ro = " << ro << " ||g||^2 = " << g2 << " " << stop << endl;
227 if (stop || g2 < reps2 ) {
228 if (kprint <= nbitermax )
229 cout << "CG converges " << iter << " ro = " << ro << " ||g||^2 = " << g2
230 << " stop=" << stop << " /Stop= " << Stop << endl;
231 return 1;// ok
232 }
233 R gamma = g2/g2p;
234 h *= gamma;
235 h -= Cg; // h = -Cg * gamma* h
236 }
237 //if (nbitermax <= nbitermax )
238 cout << "CG does'nt converge: " << nbitermax << " ||g||^2 = " << g2 << " reps2= " << reps2 << endl;
239 return 0;
240 }
241
242 template <class R>
243 class MatriceIdentite:public RNM_VirtualMatrix<R> { public:
244 typedef typename RNM_VirtualMatrix<R>::plusAx plusAx;
MatriceIdentite(int n)245 MatriceIdentite(int n) :RNM_VirtualMatrix<R>(n) {};
addMatMul(const KN_<R> & x,KN_<R> & Ax) const246 void addMatMul(const KN_<R> & x, KN_<R> & Ax) const {
247 ffassert(x.N()==Ax.N());
248 Ax+=x; }
addMatTransMul(const KN_<R> & x,KN_<R> & Ax) const249 void addMatTransMul(const KN_<R> & x, KN_<R> & Ax) const {
250 ffassert(x.N()==Ax.N());
251 Ax+=x; }
Solve(KN_<R> & y,const KN_<R> & x) const252 void Solve( KN_<R> & y ,const KN_<R> & x) const { y=x; }
WithSolver() const253 bool WithSolver() const {return true;}
operator *(const KN<R> & x) const254 typename RNM_VirtualMatrix<R>::plusAx operator*(const KN<R> & x) const {return typename RNM_VirtualMatrix<R>::plusAx(this,x);}
ChecknbLine(int n) const255 bool ChecknbLine(int n) const { return true;}
ChecknbColumn(int m) const256 bool ChecknbColumn(int m) const { return true;}
257
258 };
259
260 template<class R>
261 class SolveGCDiag : public MatriceMorse<R>::VirtualSolver , public RNM_VirtualMatrix<R>{
262 int n;
263 int nbitermax;
264 double eps;
265 mutable double epsr;
266 KN<R> D1;
267 public:
268 typedef typename RNM_VirtualMatrix<R>::plusAx plusAx;
SolveGCDiag(const MatriceMorse<R> & A,int itmax,double epsilon=1e-6)269 SolveGCDiag(const MatriceMorse<R> &A,int itmax,double epsilon=1e-6) :
270 RNM_VirtualMatrix<R>(A.n),
271 n(A.n),nbitermax(itmax?itmax: Max(100,n)),eps(epsilon),epsr(0),D1(n)
272 { //throwassert(A.sym());
273 for (int i=0;i<n;i++)
274 D1[i] = 1./A(i,i);}
Solver(const MatriceMorse<R> & a,KN_<R> & x,const KN_<R> & b) const275 void Solver(const MatriceMorse<R> &a,KN_<R> &x,const KN_<R> &b) const {
276 epsr = (eps < 0) ? (epsr >0 ? -epsr : -eps ) : eps ;
277 // cout << " epsr = " << epsr << endl;
278 ConjuguedGradient<R,MatriceMorse<R>,SolveGCDiag<R>,StopGC<R> >(a,*this,b,x,nbitermax,epsr );
279 }
operator *(const KN_<R> & x) const280 plusAx operator*(const KN_<R> & x) const {return plusAx(this,x);}
281
282
addMatMul(const KN_<R> & x,KN_<R> & Ax) const283 void addMatMul(const KN_<R> & x, KN_<R> & Ax) const
284 { ffassert(x.N()==Ax.N());
285 for (int i=0;i<n;i++)
286 Ax[i]+= D1[i]*x[i];}
287
ChecknbLine(int nn) const288 bool ChecknbLine(int nn) const { return n==nn;}
ChecknbColumn(int mm) const289 bool ChecknbColumn(int mm) const { return n==mm;}
290
291 };
292
293
294
295
296 // alloc solver and solevr def
297 //template<class R,int sympos> typename DefSparseSolverNew<R,sympos>::SparseMatSolver DefSparseSolverNew<R,sympos>::solver=0;
298 //template<class R,int sympos> typename DefSparseSolverNew<R,sympos>::SparseMatSolver DefSparseSolverNew<R,sympos>::solverdef=0;
299
300 // End Sep 2007 for generic Space solver
301
302
303
C2RR(int n,Complex * c,double * cr,double * ci)304 inline void C2RR(int n,Complex *c,double *cr,double *ci)
305 {
306 for (int i=0;i<n;i++)
307 {
308 cr[i]=real(c[i]);
309 ci[i]=imag(c[i]);
310 }
311 }
312
RR2C(int n,double * cr,double * ci,Complex * c)313 inline void RR2C(int n,double *cr,double *ci,Complex *c)
314 {
315 for (int i=0;i<n;i++)
316 {
317 c[i]=Complex(cr[i],ci[i]);
318 }
319 }
320
321
322
323 #endif
324