1 // -*- Mode : c++ -*-
2 //
3 // SUMMARY :
4 // USAGE :
5 // ORG :
6 // AUTHOR : Frederic Hecht
7 // E-MAIL : hecht@ann.jussieu.fr
8 //
9
10 /*
11
12 This file is part of Freefem++
13
14 Freefem++ is free software; you can redistribute it and/or modify
15 it under the terms of the GNU Lesser General Public License as published by
16 the Free Software Foundation; either version 2.1 of the License, or
17 (at your option) any later version.
18
19 Freefem++ is distributed in the hope that it will be useful,
20 but WITHOUT ANY WARRANTY; without even the implied warranty of
21 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 GNU Lesser General Public License for more details.
23
24 You should have received a copy of the GNU Lesser General Public License
25 along with Freefem++; if not, write to the Free Software
26 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
27 */
28 #ifndef REMOVE_CODE
29 // to do F. HECHT
30 #ifdef __MWERKS__
31 #pragma optimization_level 0
32 #endif
33 #include "ff++.hpp"
34 #include "array_resize.hpp"
35
36 #include "AFunction_ext.hpp"
37 #include "CGNL.hpp"
38
39 namespace bamg { class Triangles; }
40 namespace Fem2D { void DrawIsoT(const R2 Pt[3],const R ff[3],const RN_ & Viso); }
41
42 extern const basicForEachType *aatypeknlongp; //// for compilation error with g++ 3.2.2
43 #include "BamgFreeFem.hpp"
44
45 #include "PlotStream.hpp"
46
47 extern FILE *ThePlotStream;
48
49
50 // Debut FH Houston -------- avril 2004 ---------
51 // class for operator on sparse matrix
52 // ------------------------------------
53 // pour le addition de matrice ----matrice creuse in french)
54 // MatriceCreuse real class for matrix sparse
55 // Matrice_Creuse class for matrix sparse + poiteur on FE space def
56 // to recompute matrice in case of mesh change
57 // list<tuple<R,MatriceCreuse<R> *,bool> > * is liste of
58 // \sum_i a_i*A_i where a_i is a scalare and A_i is a sparse matrix
59 //
60
61
62
63 template<class R>
to(Matrice_Creuse<R> * M)64 list<tuple<R,MatriceCreuse<R> *, bool > > * to(Matrice_Creuse<R> * M)
65 {
66 list<tuple<R,MatriceCreuse<R> *,bool> > * l=new list<tuple<R,MatriceCreuse<R> *,bool> >;
67 l ->push_back(make_tuple<R,MatriceCreuse<R> *,bool>(1,M->A,false));
68 return l;
69 }
70 template<class R>
to(Matrice_Creuse_Transpose<R> M)71 list<tuple<R,MatriceCreuse<R> *, bool > > * to(Matrice_Creuse_Transpose<R> M)
72 {
73 list<tuple<R,MatriceCreuse<R> *,bool> > * l=new list<tuple<R,MatriceCreuse<R> *,bool> >;
74 l ->push_back(make_tuple<R,MatriceCreuse<R> *,bool>(1,M.A->A,true));
75 return l;
76 }
77
78 template<class R>
M2L3(Stack,const AnyType & pp)79 AnyType M2L3 (Stack , const AnyType & pp)
80 {
81 return to(GetAny<Matrice_Creuse<R> *>(pp));
82 }
83
84
85 template<class R>
tM2L3(Stack,const AnyType & pp)86 AnyType tM2L3 (Stack , const AnyType & pp)
87 {
88 return to(GetAny<Matrice_Creuse_Transpose<R> >(pp));
89 }
90
91
92 template<class R>
93 struct Op2_ListCM: public binary_function<R,Matrice_Creuse<R> *,list<tuple<R,MatriceCreuse<R> *,bool> > *>
94 {
95 typedef tuple<R,MatriceCreuse<R> *,bool> P;
96 typedef list<P> L;
97 typedef L * RR;
98 typedef R AA;
99 typedef Matrice_Creuse<R> * BB;
100
fOp2_ListCM101 static RR f(const AA & a,const BB & b)
102 {
103 RR r= new list<P> ;
104 P p(a,b->pMC(),false);
105 r ->push_back(p);
106 return r;}
107 };
PrintL(const char * cc,list<tuple<R,VirtualMatrix<int,R> *,bool>> const * const lM)108 template<class R> void PrintL(const char* cc, list<tuple<R,VirtualMatrix<int,R>*,bool> > const * const lM)
109 {
110 if(verbosity<100) return;
111 typedef typename list<tuple<R,VirtualMatrix<int,R> *,bool> >::const_iterator lconst_iterator;
112
113 lconst_iterator begin=lM->begin();
114 lconst_iterator end=lM->end();
115 lconst_iterator i;
116 cout << cc << " (" ;
117 for(i=begin;i!=end;i++++)
118 {
119 if(std::get<1>(*i)) // M == 0 => zero matrix
120 {
121 VirtualMatrix<int,R>& M=*std::get<1>(*i);
122 bool transpose = std::get<2>(*i) ;
123 R coef=std::get<0>(*i);
124 cout << " + " << coef << "*" << &M <<"^" << transpose ;
125
126 }
127 }
128
129
130 cout << ") "<< endl;
131 }
132 template<class R>
133 struct Op2_ListMC: public binary_function<Matrice_Creuse<R> *,R,list<tuple<R,MatriceCreuse<R> *,bool> > *>
134 {
135 typedef tuple<R,MatriceCreuse<R> *,bool> P;
136 typedef list<P> L;
137 typedef L * RR;
138 typedef R AA;
139 typedef Matrice_Creuse<R> * BB;
140
fOp2_ListMC141 static RR f(const BB & b,const AA & a)
142 {
143 RR r= new list<P> ;
144 P p(a,b->pMC(),false);
145 r ->push_back(p);
146 return r;}
147 };
148 // ADD FH 16/02/2007
149
150 template<class R>
151 struct Op2_ListCMt: public binary_function<R,Matrice_Creuse_Transpose<R> ,list<tuple<R,MatriceCreuse<R> *,bool> > *>
152 {
153 typedef tuple<R,MatriceCreuse<R> *,bool> P;
154 typedef list<P> L;
155 typedef L * RR;
156 typedef R AA;
157 typedef Matrice_Creuse_Transpose<R> BB;
158
fOp2_ListCMt159 static RR f(const AA & a,const BB & b)
160 {
161 RR r= new list<P> ;
162 P p(a,b.A->pMC(),true);
163 r ->push_back(p);
164 return r;}
165 };
166
167 template<class R>
168 struct Op2_ListMtC: public binary_function<Matrice_Creuse_Transpose<R> ,R,list<tuple<R,MatriceCreuse<R> *,bool> > *>
169 {
170 typedef tuple<R,MatriceCreuse<R> *,bool> P;
171 typedef list<P> L;
172 typedef L * RR;
173 typedef R AA;
174 typedef Matrice_Creuse_Transpose<R> BB;
175
fOp2_ListMtC176 static RR f(const BB & b,const AA & a)
177 {
178 RR r= new list<P> ;
179 P p(a,b.A->pMC(),true);
180 r ->push_back(p);
181 return r;}
182 };
183 // FIN ADD 16/02/2007
184
185
186
187 template<class R,int c=-1>
188 struct Op1_LCMd: public unary_function<list<tuple<R,MatriceCreuse<R> *,bool> > *,
189 list<tuple<R,MatriceCreuse<R> *,bool> > * >
190 { // - ...
191 typedef tuple<R,MatriceCreuse<R> *,bool> P;
192 typedef list<P> L;
193 typedef L * RR;
194
fOp1_LCMd195 static RR f(const RR & l)
196 {
197 typedef typename list<tuple<R,MatriceCreuse<R> *,bool> >::iterator lci;
198 for(lci i= l->begin();i !=l->end();++i)
199 get<0>(*i) *= R(c);
200 PrintL(" - Op1_LCMd: ",l);
201 return l;
202 }
203
204 };
205
206 template<class R>
207 struct Op2_ListCMCMadd: public binary_function<list<tuple<R,MatriceCreuse<R> *,bool> > *,
208 list<tuple<R,MatriceCreuse<R> *,bool> > *,
209 list<tuple<R,MatriceCreuse<R> *,bool> > * >
210 { // ... + ...
211 typedef tuple<R,MatriceCreuse<R> *,bool> P;
212 typedef list<P> L;
213 typedef L * RR;
214
fOp2_ListCMCMadd215 static RR f(const RR & a,const RR & b)
216 {
217 a->insert(a->end(),b->begin(),b->end());
218
219 delete b;
220 return a;
221 }
222
223 };
224 template<class R>
225 struct Op2_ListCMCMsub: public binary_function<list<tuple<R,MatriceCreuse<R> *,bool> > *,
226 list<tuple<R,MatriceCreuse<R> *,bool> > *,
227 list<tuple<R,MatriceCreuse<R> *,bool> > * >
228 { // ... + ...
229 typedef tuple<R,MatriceCreuse<R> *,bool> P;
230 typedef list<P> L;
231 typedef L * RR;
232
fOp2_ListCMCMsub233 static RR f(const RR & a,const RR & b)
234 {
235 Op1_LCMd<R,-1>::f(b);
236 PrintL("Op2_ListCMCMsub +",a);
237 PrintL(" -(-) ",b);
238
239 a->insert(a->end(),b->begin(),b->end());
240 PrintL(" => ",a);
241 delete b;
242 return a;
243 }
244
245 };
246
247 template<class R,int cc=1>
248 struct Op2_ListMCMadd: public binary_function<Matrice_Creuse<R> *,
249 list<tuple<R,MatriceCreuse<R> *,bool> > *,
250 list<tuple<R,MatriceCreuse<R> *,bool> > * >
251 { // M + ....
252 typedef tuple<R,MatriceCreuse<R> *,bool> P;
253 typedef list<P> L;
254 typedef L * RR;
255 typedef Matrice_Creuse<R> * MM;
256
fOp2_ListMCMadd257 static RR f(const MM & a,const RR & b)
258 {
259 // M + c*L
260 Op1_LCMd<R,cc>::f(b);
261 PrintL("Op2_ListMCMadd M +",b);
262
263 b->push_front(make_tuple<R,MatriceCreuse<R> *,bool>(R(1.),a->A,false));
264 return b;
265 }
266
267
268 };
269
270 template<class R,int cc=1>
271 struct Op2_ListCMMadd: public binary_function< list<tuple<R,MatriceCreuse<R> *,bool> > *,
272 Matrice_Creuse<R> * ,
273 list<tuple<R,MatriceCreuse<R> *,bool> > *>
274 { // .... + M
275 typedef tuple<R,MatriceCreuse<R> *,bool> P;
276 typedef list<P> L;
277 typedef L * RR;
278 typedef Matrice_Creuse<R> * MM;
279
fOp2_ListCMMadd280 static RR f(const RR & a,const MM & b)
281 {
282 // L + c*M
283 a->push_back(make_tuple<R,MatriceCreuse<R> *,bool>(R(cc),b->A,false));
284 return a;
285 }
286
287
288 };
289
290 template<class R,int cc=1>
291 struct Op2_ListMMadd: public binary_function< Matrice_Creuse<R> *,
292 Matrice_Creuse<R> * ,
293 list<tuple<R,MatriceCreuse<R> *,bool> > *>
294 { // M + M
295 typedef tuple<R,MatriceCreuse<R> *,bool> P;
296 typedef list<P> L;
297 typedef L * RR;
298 typedef Matrice_Creuse<R> * MM;
299
fOp2_ListMMadd300 static RR f(const MM & a,const MM & b)
301 {
302 // M + c M
303 L * l=to(a);
304 l->push_back(make_tuple<R,MatriceCreuse<R> *>(R(cc),b->A,false));
305 return l;
306 }
307
308
309 };
310 // Fin Add FH Houston --------
311
312 // for Jolivet to build restriction jan 2014
313 // t[int] I= restrict(VCh,VGh,IPG); // ou
314 template<class pfes>
315 class RestrictArray : public OneOperator { public:
316 template< typename T > struct Base { typedef T B; };
317 template< typename T > struct Base< T* >{ typedef T B;};
318
319 typedef typename Base<pfes>::B::FESpace FESpace;
320 typedef typename FESpace::FElement FElement;
321
322 class Op : public E_F0info { public: // passe de l'info ..
323 typedef pfes * A;
324 Expression a,b,c,d;
325 static const int n_name_param =0;
326 Expression nargs[n_name_param];
arg(int i,Stack stack,bool a) const327 bool arg(int i,Stack stack,bool a) const{ return nargs[i] ? GetAny<bool>( (*nargs[i])(stack) ): a;}
arg(int i,Stack stack,long a) const328 long arg(int i,Stack stack,long a) const{ return nargs[i] ? GetAny<long>( (*nargs[i])(stack) ): a;}
arg(int i,Stack stack,KN_<long> a) const329 KN_<long> arg(int i,Stack stack,KN_<long> a ) const{ return nargs[i] ? GetAny<KN_<long> >( (*nargs[i])(stack) ): a;}
330
331 public:
Op(const basicAC_F0 & args,Expression aa,Expression bb,Expression cc)332 Op(const basicAC_F0 & args,Expression aa,Expression bb,Expression cc) : a(aa),b(bb),c(cc),d(0) {
333 args.SetNameParam(n_name_param,0,nargs);
334
335 }
336 };
RestrictArray()337 RestrictArray() : OneOperator( atype<const typename RestrictArray<pfes>::Op *>(),//atype<KN<long>* >(),
338 atype<pfes *>(),
339 atype<pfes *>(),
340 atype<KN_<long> >()) {}
341
code(const basicAC_F0 & args) const342 E_F0 * code(const basicAC_F0 & args) const
343 {
344 if(args.size()!=3)CompileError("Bug in RestrictArray code nb of args != 3 ????? bizarre" );
345 return new Op(args,t[0]->CastTo(args[0]),
346 t[1]->CastTo(args[1]),
347 t[2]->CastTo(args[2]));
348 }
349 };
350 // end restrict ...
351 template< typename T > struct Base { typedef T B; };
352 template< typename T > struct Base< T* >{ typedef T B;};
353
354 template<class pfes, int INIT>
SetRestrict(Stack stack,Expression einj,Expression erest)355 AnyType SetRestrict(Stack stack,Expression einj,Expression erest)
356 {
357
358 typedef typename Base<pfes>::B::FESpace FESpace;
359 typedef typename FESpace::FElement FElement;
360
361 KN<long> * pinj =GetAny<KN<long>*>((*einj)(stack));
362 const typename RestrictArray<pfes>::Op * ar(dynamic_cast<const typename RestrictArray<pfes>::Op *>(erest));
363 ffassert(ar);
364 if( verbosity>9) cout << " -- RestrictArray "<< endl;
365 pfes * pCUh = GetAny< pfes * >((* ar->a)(stack));
366 pfes * pFVh = GetAny< pfes * >((* ar->b)(stack));
367 // verif same FE.
368 KN_<long> ncf= GetAny< KN_<long> >((* ar->c)(stack));
369 FESpace * pVCh = **pCUh;
370 FESpace * pVFh = **pFVh;
371 FESpace & VCh = *pVCh;
372 FESpace & VFh = *pVFh;
373 long neC = VCh.NbOfElements ;
374 long neF = VFh.NbOfElements ;
375 long ndfC = VCh.NbOfDF ;
376 long ndfF = VFh.NbOfDF ;
377
378 KN_<long> nc2f= ncf;
379 if(INIT==0)
380 pinj->init(ndfC);
381 else pinj->resize(ndfC);
382 KN<long> & inj=*pinj;
383 inj = -1; // un set ..
384 if( verbosity>9) cout<< " ne =" << neC << " " << neF << endl;
385
386 for(int kc=0; kc <VCh.NbOfElements; kc++)
387 {
388
389 int kf = nc2f(kc);
390 FElement KC(pVCh,kc);
391 FElement KF(pVFh,kf);
392
393 int ndofKC = KC.NbDoF() ;
394 int ndofKF = KF.NbDoF() ;
395 if( verbosity>99) cout << kc << " " << kf << " : " <<ndofKC << " " << ndofKF << endl;
396 ffassert(ndofKC== ndofKF );
397 ffassert( kf >= 0 && kf < neF);
398 ffassert( kc >= 0 && kc< neC);
399
400 for(int df=0; df<ndofKC; df++)
401 {
402 int dfC =KC(df), dfF=KF(df);
403 if( verbosity>99) cout << dfC <<" -> "<< dfF << endl;
404 assert(dfC >= 0 && dfC < ndfC);
405 inj[dfC] = dfF;
406 }
407
408 }
409 if( verbosity>9) cout << " restrict: Inject= " << inj << endl;
410 ffassert(inj.min() != -1);
411
412 return pinj;
413 }
414
415 // Fin Add FH Houston --------
416 template<class pfes1,class pfes2>
417 class MatrixInterpolation : public OneOperator { public:
418
419 class Op : public E_F0info { public:
420 //typedef pfes * A;
421 Expression a,b,c,d;
422 static const int n_name_param =5;
423 static basicAC_F0::name_and_type name_param[] ;
424 Expression nargs[n_name_param];
arg(int i,Stack stack,bool a) const425 bool arg(int i,Stack stack,bool a) const{ return nargs[i] ? GetAny<bool>( (*nargs[i])(stack) ): a;}
arg(int i,Stack stack,long a) const426 long arg(int i,Stack stack,long a) const{ return nargs[i] ? GetAny<long>( (*nargs[i])(stack) ): a;}
arg(int i,Stack stack,KN_<long> a) const427 KN_<long> arg(int i,Stack stack,KN_<long> a ) const{ return nargs[i] ? GetAny<KN_<long> >( (*nargs[i])(stack) ): a;}
428
429 public:
Op(const basicAC_F0 & args,Expression aa,Expression bb)430 Op(const basicAC_F0 & args,Expression aa,Expression bb) : a(aa),b(bb),c(0),d(0) {
431 args.SetNameParam(n_name_param,name_param,nargs); }
Op(const basicAC_F0 & args,Expression aa,Expression bb,Expression cc)432 Op(const basicAC_F0 & args,Expression aa,Expression bb,Expression cc) : a(aa),b(bb),c(cc),d(0) {
433 args.SetNameParam(n_name_param,name_param,nargs); }
Op(const basicAC_F0 & args,Expression aa,Expression bb,Expression cc,Expression dd)434 Op(const basicAC_F0 & args,Expression aa,Expression bb,Expression cc,Expression dd) : a(aa),b(bb),c(cc),d(dd) {
435 args.SetNameParam(n_name_param,name_param,nargs); }
436
437
438 };
439 // interpolation(Vh,Vh)
MatrixInterpolation()440 MatrixInterpolation() : OneOperator(atype<const typename MatrixInterpolation<pfes1,pfes2>::Op *>(),
441 atype<pfes1 *>(),
442 atype<pfes2 *>()) {}
443 // interpolation(Vh,xx,yy) // 2d
MatrixInterpolation(int bidon)444 MatrixInterpolation(int bidon) : OneOperator(atype<const typename MatrixInterpolation<pfes1,pfes1>::Op *>(),
445 atype<pfes1 *>(),atype<KN_<double> >(),atype<KN_<double> >()) {}
446
447 // interpolation(Vh,xx,yy,zz) // 3d
MatrixInterpolation(int bidon,int bidon2)448 MatrixInterpolation(int bidon,int bidon2) : OneOperator(atype<const typename MatrixInterpolation<pfes1,pfes1>::Op *>(),
449 atype<pfes1 *>(),atype<KN_<double> >(),atype<KN_<double> >(),atype<KN_<double> >()) {}
450
451
code(const basicAC_F0 & args) const452 E_F0 * code(const basicAC_F0 & args) const
453 {
454 if(args.size()==2)
455 return new Op(args,t[0]->CastTo(args[0]),
456 t[1]->CastTo(args[1]));
457 else if(args.size()==3)
458 return new Op(args,t[0]->CastTo(args[0]),
459 t[1]->CastTo(args[1]),
460 t[2]->CastTo(args[2]));
461 else if(args.size()==4)
462 return new Op(args,t[0]->CastTo(args[0]),
463 t[1]->CastTo(args[1]),
464 t[2]->CastTo(args[2]),
465 t[2]->CastTo(args[3])
466 );
467 else CompileError("Bug in MatrixInterpolation code nb != 2 or 3 ????? bizarre" );
468 return 0;
469 }
470 };
471
472 template<class pfes1,class pfes2>
473 basicAC_F0::name_and_type MatrixInterpolation<pfes1,pfes2>::Op::name_param[]= {
474 { "t", &typeid(bool)},
475 { "op", &typeid(long)},
476 { "inside",&typeid(bool)},
477 { "composante",&typeid(long)},
478 { "U2Vc",&typeid(KN_<long>)}
479
480 };
481
482 /*template<>
483 basicAC_F0::name_and_type MatrixInterpolation<pfes3,pfes3>::Op::name_param[]= {
484 { "t", &typeid(bool)},
485 { "op", &typeid(long)},
486 { "inside",&typeid(bool)},
487 { "composante",&typeid(long)},
488 { "U2Vc",&typeid(KN_<long>)}
489
490 };
491
492 template<>
493 basicAC_F0::name_and_type MatrixInterpolation<pfesS,pfesS>::Op::name_param[]= {
494 { "t", &typeid(bool)},
495 { "op", &typeid(long)},
496 { "inside",&typeid(bool)},
497 { "composante",&typeid(long)},
498 { "U2Vc",&typeid(KN_<long>)}
499
500 };
501
502 template<>
503 basicAC_F0::name_and_type MatrixInterpolation<pfesL,pfesL>::Op::name_param[]= {
504 { "t", &typeid(bool)},
505 { "op", &typeid(long)},
506 { "inside",&typeid(bool)},
507 { "composante",&typeid(long)},
508 { "U2Vc",&typeid(KN_<long>)}
509
510 };
511
512 template<>
513 basicAC_F0::name_and_type MatrixInterpolation<pfes3,pfesS>::Op::name_param[]= {
514 { "t", &typeid(bool)},
515 { "op", &typeid(long)},
516 { "inside",&typeid(bool)},
517 { "composante",&typeid(long)},
518 { "U2Vc",&typeid(KN_<long>)}
519
520 };
521
522 template<>
523 basicAC_F0::name_and_type MatrixInterpolation<pfesL,pfesS>::Op::name_param[]= {
524 { "t", &typeid(bool)},
525 { "op", &typeid(long)},
526 { "inside",&typeid(bool)},
527 { "composante",&typeid(long)},
528 { "U2Vc",&typeid(KN_<long>)}
529
530 };
531
532 template<>
533 basicAC_F0::name_and_type MatrixInterpolation<pfesL,pfes>::Op::name_param[]= {
534 { "t", &typeid(bool)},
535 { "op", &typeid(long)},
536 { "inside",&typeid(bool)},
537 { "composante",&typeid(long)},
538 { "U2Vc",&typeid(KN_<long>)}
539
540 };
541 */
542 template<class R>
543 class SetMatrix_Op : public E_F0mps { public:
544 Expression a;
545
546 static aType btype;
547 static const int n_name_param =NB_NAME_PARM_MAT; // add nbiter FH 30/01/2007 11 -> 12 //add var MUMPS+autre
548 static basicAC_F0::name_and_type name_param[] ;
549 Expression nargs[n_name_param];
550 const OneOperator * precon;
arg(int i,Stack stack,bool a) const551 bool arg(int i,Stack stack,bool a) const{ return nargs[i] ? GetAny<bool>( (*nargs[i])(stack) ): a;}
arg(int i,Stack stack,long a) const552 long arg(int i,Stack stack,long a) const{ return nargs[i] ? GetAny<long>( (*nargs[i])(stack) ): a;}
553
554 public:
SetMatrix_Op(const basicAC_F0 & args,Expression aa)555 SetMatrix_Op(const basicAC_F0 & args,Expression aa) : a(aa) {
556 args.SetNameParam(n_name_param,name_param,nargs);
557 precon = 0; // a changer
558 if ( nargs[3])
559 {
560 const Polymorphic * op= dynamic_cast<const Polymorphic *>(nargs[3]);
561 assert(op);
562 precon = op->Find("(",ArrayOfaType(atype<KN<R>* >(),false)); // strange bug in g++ is R become a double
563 ffassert(precon);
564 }
565
566 }
567 AnyType operator()(Stack stack) const ;
568 };
569
570
571 template<class R>
572 class SetMatrix : public OneOperator { public:
573
SetMatrix()574 SetMatrix() : OneOperator(SetMatrix_Op<R>::btype,atype<Matrice_Creuse<R> *>() ) {}
575
code(const basicAC_F0 & args) const576 E_F0 * code(const basicAC_F0 & args) const
577 {
578 return new SetMatrix_Op<R>(args,t[0]->CastTo(args[0]));
579 }
580 };
581
582 template<class R>
583 aType SetMatrix_Op<R>::btype=0;
584
585 template <class R>
586 basicAC_F0::name_and_type SetMatrix_Op<R>::name_param[]= {
587 LIST_NAME_PARM_MAT
588
589 };
590
591 template<class R>
operator ()(Stack stack) const592 AnyType SetMatrix_Op<R>::operator()(Stack stack) const
593 {
594 Matrice_Creuse<R> * A= GetAny<Matrice_Creuse<R> *>((*a)(stack));
595
596 ffassert(A);
597 Data_Sparse_Solver ds;
598 bool VF=false;
599 ds.factorize=0;
600 // get previous value of sym ??? FH & PHT fev 2020
601 int syma=-1;
602 if( A->A)
603 { HashMatrix<int,R> * phm= A->pHM();
604 if(phm)
605 syma=phm->half;
606 }
607
608 if( !A->A) A->A.master(new MatriceMorse<R>(0,0,0,0));// set to empty matrix .. mars 2014 FH ..
609 SetEnd_Data_Sparse_Solver<R>(stack,ds,nargs,n_name_param,syma);
610 if( verbosity >4) cout <<" SetMatrix_Op " << ds.sym << " "<< ds.positive << " " << syma << endl;
611 VirtualMatrix<int,R> *pvm =A->pMC();
612 ffassert(pvm);
613 pvm->setsdp(ds.sym,ds.positive); // put the matrix in rigth format
614 SetSolver<R>(stack,VF,*A->pMC(),ds);
615
616 return Nothing;
617 }
618
619
620
621 template<int init>
622 AnyType SetMatrixInterpolation(Stack,Expression ,Expression);
623 template<int init>
624 AnyType SetMatrixInterpolation3(Stack,Expression ,Expression);
625 template<int init>
626 AnyType SetMatrixInterpolationS(Stack,Expression ,Expression);
627
628 template<class R>
BuildCombMat(MatriceMorse<R> & mij,const KNM_<R> & A,int ii00=0,int jj00=0,R coef=R (1.),bool cnj=false)629 void BuildCombMat(MatriceMorse<R> & mij,const KNM_<R> & A, int ii00=0,int jj00=0,R coef=R(1.),bool cnj=false)
630 {
631 double eps0=numeric_limits<double>::min();
632 int i,j;
633 int n = A.N(),m=A.M();
634 for ( i=0;i<n;i++)
635 for ( j=0;j<m;j++)
636 {
637 R cij=coef*A(i,j);
638 if (cnj) cij = RNM::conj(cij);
639 mij[ij_mat(false,ii00,jj00,i,j)] += cij;
640
641 }
642
643 }
644
buildInterpolationMatrix(const FESpace & Uh,const FESpace & Vh,void * data)645 MatriceMorse<R> * buildInterpolationMatrix(const FESpace & Uh,const FESpace & Vh,void *data)
646 { // Uh = Vh
647 MatriceMorse<R> * m=0;
648 int op=op_id; // value of the function
649 bool transpose=false;
650 bool inside=false;
651 int * iU2V=0;
652 if (data)
653 {
654 int * idata=static_cast<int*>(data);
655 transpose=idata[0];
656 op=idata[1];
657 inside=idata[2];
658 iU2V= idata + 4;
659 ffassert(op>=0 && op < 4);
660 }
661 if(verbosity>2)
662 {
663 cout << " -- buildInterpolationMatrix transpose =" << transpose << endl
664 << " value, dx , dy op = " << op << endl
665 << " just inside = " << inside << endl;
666 }
667 using namespace Fem2D;
668 int n=Uh.NbOfDF;
669 int mm=Vh.NbOfDF;
670 if(transpose) Exchange(n,mm);
671 m = new MatriceMorse<R>(n,mm,0,0);
672 const Mesh & ThU =Uh.Th; // line
673 const Mesh & ThV =Vh.Th; // colunm
674 bool samemesh = &Uh.Th == &Vh.Th; // same Mesh
675 int thecolor =0;
676
677 int nnz =0;
678
679 KN<int> color(ThV.nt);
680 KN<int> mark(n);
681 mark=0;
682
683 color=thecolor++;
684 FElement Uh0 = Uh[0];
685 FElement Vh0 = Vh[0];
686
687 FElement::aIPJ ipjU(Uh0.Pi_h_ipj());
688 FElement::aR2 PtHatU(Uh0.Pi_h_R2());
689
690 int nbdfVK= Vh0.NbDoF();
691 int NVh= Vh0.N;
692
693 int nbp= PtHatU.N();
694 KN<R2> PV(nbp); // the PtHat in ThV mesh
695 KN<int> itV(nbp); // the Triangle number
696 KN<bool> intV(nbp); // ouside or not
697 KN<R> AipjU(ipjU.N());
698
699 KNM<R> aaa(nbp,nbdfVK);
700
701
702 const R eps = 1.0e-10;
703 const int sfb1=Vh0.N*last_operatortype*Vh0.NbDoF();
704 KN<R> kv(sfb1*nbp);
705 R * v = kv;
706 KN<int> ik(nbp); // the Triangle number
707
708 bool whatd[last_operatortype];
709 for (int i=0;i<last_operatortype;i++)
710 whatd[i]=false;
711 whatd[op]=true; // the value of function
712 KN<bool> fait(Uh.NbOfDF);
713 fait=false;
714 R2 Gh(1./3,1./3);
715 {
716
717 for (int it=0;it<ThU.nt;it++)
718 {
719 thecolor++; // change the current color
720 const Triangle & TU(ThU[it]);
721 FElement KU(Uh[it]);
722 KU.Pi_h(AipjU);
723 if (samemesh)
724 {
725 PV = PtHatU;
726 itV = it;
727 intV= false;// add July 2009 (unset varaible FH)
728 }
729 else
730 {
731 const Triangle *ts=0,*ts0=0;
732 // Search barycenter
733 bool outside;
734 R2 G;// Add frev 2019 more cleaver F. Hecht
735 ts0=ThV.Find(TU(Gh),G,outside,ts0);// find barycenter good if imbricated mesh ....
736 if(outside) ts0=0; // not find
737 for (int i=0;i<nbp;i++)
738 {
739 ts=ThV.Find(TU(PtHatU[i]),PV[i],outside,ts0);
740 if(!outside && ts0==0) ts0=ts;
741 if(outside && verbosity>9 )
742 cout << it << " " << i << " :: " << TU(PtHatU[i]) << " -- "<< outside << PV[i] << " " << ThV(ts) << " -> " << (*ts)(PV[i]) <<endl;
743 itV[i]= ThV(ts);
744 intV[i]=outside && inside; // ouside and inside flag
745 }
746 }
747
748
749 for (int p=0;p<nbp;p++)
750 {
751
752 KNMK_<R> fb(v+p*sfb1,nbdfVK,NVh,last_operatortype); // valeur de fonction de base de Vh
753 // ou: fb(idf,j,0) valeur de la j composante de la fonction idf
754 Vh0.tfe->FB(whatd,ThV,ThV[itV[p]],PV[p],fb);
755 }
756
757 for (int i=0;i<ipjU.N();i++)
758 { // pour tous le terme
759 const FElement::IPJ &ipj_i(ipjU[i]);
760 int dfu = KU(ipj_i.i); // le numero de df global
761 if(fait[dfu]) continue;
762 int jU = ipj_i.j; // la composante dans U
763 int p=ipj_i.p; // le points
764 if (intV[p]) continue; // ouside and inside flag => next
765 R aipj = AipjU[i];
766 FElement KV(Vh[itV[p]]);
767 int jV=jU;
768 if(iU2V) jV=iU2V[jU];
769
770 if(jV>=0 && jV<NVh)
771 {
772 KNMK_<R> fb(v+p*sfb1,nbdfVK,NVh,last_operatortype);
773 KN_<R> fbj(fb('.',jV,op));
774
775 for (int idfv=0;idfv<nbdfVK;idfv++)
776 if (Abs(fbj[idfv])>eps)
777 {
778 int dfv=KV(idfv);
779 int ii=dfu, jj=dfv;
780 if(transpose) Exchange(ii,jj);
781 // le term dfu,dfv existe dans la matrice
782 R c= fbj[idfv]*aipj;
783 if(Abs(c)>eps)
784 (*m)(ii,jj) += c;
785 }
786 }
787
788 }
789
790 for (int df=0;df<KU.NbDoF();df++)
791 {
792 int dfu = KU(df); // le numero de df global
793 fait[dfu]=true;
794 }
795 }
796 }
797 // sij.clear();
798 return m;
799 }
800
801 template<typename RdHatT1,typename RdHatT2>
copyPt(KN<RdHatT1> & A,KN<RdHatT2> B)802 void copyPt( KN<RdHatT1> &A, KN<RdHatT2> B)
803 { A=B; };
804
805 template<>
copyPt(KN<R2> & A,KN<R1> B)806 void copyPt<R2,R1>( KN<R2> &A, KN<R1> B) {
807 for( int ik=0;ik<A.N();ik++)
808 A[ik].x = B[ik].x;
809 }
810 template<>
copyPt(KN<R3> & A,KN<R2> B)811 void copyPt<R3,R2>( KN<R3> &A, KN<R2> B) {
812 for( int ik=0;ik<A.N();ik++){
813 A[ik].x = B[ik].x;
814 A[ik].y = B[ik].y;
815 }
816 }
817
818 template<class FESpaceT1,class FESpaceT2>
buildInterpolationMatrixT(const FESpaceT1 & Uh,const FESpaceT2 & Vh,void * data)819 MatriceMorse<R> * buildInterpolationMatrixT(const FESpaceT1 & Uh,const FESpaceT2 & Vh,void *data)
820 {
821 MatriceMorse<R> * m=0;
822 typedef typename FESpaceT1::Mesh Mesh1;
823 typedef typename FESpaceT1::FElement FElement1;
824 typedef typename Mesh1::Element Element1;
825 typedef typename FESpaceT1::Rd Rd1;
826 typedef typename Element1::RdHat RdHat1;
827
828 typedef typename FESpaceT2::Mesh Mesh2;
829 typedef typename FESpaceT2::FElement FElement2;
830 typedef typename Mesh2::Element Element2;
831 typedef typename FESpaceT2::Rd Rd2;
832 typedef typename Element2::RdHat RdHat2;
833
834
835 int op=op_id; // value of the function
836 bool transpose=false;
837 bool inside=false;
838 int * iU2V=0;
839 if (data)
840 {
841 int * idata=static_cast<int*>(data);
842 transpose=idata[0];
843 op=idata[1];
844 inside=idata[2];
845 iU2V= idata + 4;
846 ffassert(op>=0 && op < 4);
847 }
848 if(verbosity>2)
849 {
850 cout << " -- buildInterpolationMatrix transpose =" << transpose << endl
851 << " value, dx , dy op = " << op << endl
852 << " just inside = " << inside << endl;
853 }
854 using namespace Fem2D;
855 int n=Uh.NbOfDF;
856 int mm=Vh.NbOfDF;
857 if(transpose) Exchange(n,mm);
858 m = new MatriceMorse<R>(n,mm,0,0);
859
860 RdHat1 Gh= RdHat1::diag(1./(RdHat1::d+1));
861 RdHat2 G;
862
863 int n1=n+1;
864 const Mesh1 & ThU =Uh.Th; // line
865 const Mesh2 & ThV =Vh.Th; // colunm
866 bool samemesh = (is_same< Mesh1, Mesh2 >::value) ? (void*)&Uh.Th == (void*)&Vh.Th : 0 ; // same Mesh
867 int thecolor =0;
868
869 int nnz =0;
870
871 KN<int> color(ThV.nt);
872 KN<int> mark(n);
873 mark=0;
874
875 int * cl = 0;
876 double *a=0;
877
878 color=thecolor++;
879 FElement1 Uh0 = Uh[0];
880 FElement2 Vh0 = Vh[0];
881
882
883 int nbdfVK= Vh0.NbDoF();
884 int NVh= Vh0.N;
885
886 InterpolationMatrix<RdHat1> ipmat(Uh);
887
888 int nbp=ipmat.np; //
889 KN<RdHat2> PV(nbp); // the PtHat in ThV mesh
890 KN<int> itV(nbp); // the Triangle number
891 KN<bool> intV(nbp); // ouside or not
892
893 KNM<R> aaa(nbp,nbdfVK);
894
895
896 const R eps = 1.0e-10;
897 const int sfb1=Vh0.N*last_operatortype*Vh0.NbDoF();
898 KN<R> kv(sfb1*nbp);
899 R * v = kv;
900 KN<int> ik(nbp); // the Triangle number
901 op= op==3 ? op_dz : op; // renumber op ???? dec 2010 FH.
902 What_d whatd= 1<< op;
903 KN<bool> fait(Uh.NbOfDF);
904 fait=false;
905 {
906
907 for (int it=0;it<ThU.nt;it++)
908 {
909 thecolor++; // change the current color
910 const Element1 & TU(ThU[it]);
911 FElement1 KU(Uh[it]);
912 ipmat.set(KU);
913 if (samemesh) {
914 copyPt<RdHat2,RdHat1>(PV,ipmat.P);
915 itV = it;
916 intV= false;// add July 2009 (unset varaible FH)
917 }
918 else
919 {
920 const Element2 *ts=0,*ts0=0;
921 bool outside;
922 ts0=ThV.Find(TU(Gh),G,outside,ts0);
923 if(outside) ts0=0; // bad starting tet
924 for (int i=0;i<nbp;i++)
925 {
926 ts=ThV.Find(TU(ipmat.P[i]),PV[i],outside,ts0);
927 if( ts0 ==0 && !outside) ts0=ts;
928 if(outside && verbosity>9 )
929 cout << it << " " << i << " :: " << TU(ipmat.P[i]) << " -- "<< outside << PV[i] << " " << ThV(ts) << " -> " << (*ts)(PV[i]) <<endl;
930 itV[i]= ThV(ts);
931 intV[i]=outside && inside; // ouside and inside flag
932 }
933 }
934
935 for (int p=0;p<nbp;p++)
936 {
937 KNMK_<R> fb(v+p*sfb1,nbdfVK,NVh,last_operatortype); // valeur de fonction de base de Vh
938 // ou: fb(idf,j,0) valeur de la j composante de la fonction idf
939 Vh0.tfe->FB(whatd,ThV,ThV[itV[p]],PV[p],fb);
940 }
941
942 for (int i=0;i<ipmat.ncoef;i++)
943 { // pour tous le terme
944 int dfu = KU(ipmat.dofe[i]); // le numero de df global
945 if(fait[dfu]) continue;
946 int jU = ipmat.comp[i]; // la composante dans U
947 int p=ipmat.p[i]; // le point
948 if (intV[p]) continue; // ouside and inside flag => next
949 R aipj = ipmat.coef[i];
950 FElement2 KV(Vh[itV[p]]);
951 int jV=jU;
952 if(iU2V) jV=iU2V[jU];
953
954 if(jV>=0 && jV<NVh)
955 {
956 KNMK_<R> fb(v+p*sfb1,nbdfVK,NVh,last_operatortype);
957 KN_<R> fbj(fb('.',jV,op));
958
959 for (int idfv=0;idfv<nbdfVK;idfv++)
960 if (Abs(fbj[idfv])>eps)
961 {
962 int dfv=KV(idfv);
963 int ii=dfu, jj=dfv;
964 if(transpose) Exchange(ii,jj);
965 // le term dfu,dfv existe dans la matrice
966 R c= fbj[idfv]*aipj;
967 if(Abs(c)>eps)
968 (*m)(ii,jj) += c;
969 }
970 }
971
972 }
973
974 for (int df=0;df<KU.NbDoF();df++)
975 {
976 int dfu = KU(df); // le numero de df global
977 fait[dfu]=true;
978 }
979
980
981 }
982
983 }
984 return m;
985 }
986
987
988 template< >
buildInterpolationMatrixT(const FESpaceL & Uh,const FESpace & Vh,void * data)989 MatriceMorse<R> * buildInterpolationMatrixT<FESpaceL,FESpace>(const FESpaceL & Uh,const FESpace & Vh,void *data)
990 {
991 MatriceMorse<R> * m=0;
992 typedef typename FESpaceL::Mesh Mesh1;
993 typedef typename FESpaceL::FElement FElement1;
994 typedef typename Mesh1::Element Element1;
995 typedef typename FESpaceL::Rd Rd1;
996 typedef typename Element1::RdHat RdHat1;
997
998 typedef typename FESpace::Mesh Mesh2;
999 typedef typename FESpace::FElement FElement2;
1000 typedef typename Mesh2::Element Element2;
1001 typedef typename FESpace::Rd Rd2;
1002 typedef typename Element2::RdHat RdHat2;
1003
1004
1005 int op=op_id; // value of the function
1006 bool transpose=false;
1007 bool inside=false;
1008 int * iU2V=0;
1009 if (data)
1010 {
1011 int * idata=static_cast<int*>(data);
1012 transpose=idata[0];
1013 op=idata[1];
1014 inside=idata[2];
1015 iU2V= idata + 4;
1016 ffassert(op>=0 && op < 4);
1017 }
1018 if(verbosity>2)
1019 {
1020 cout << " -- buildInterpolationMatrix transpose =" << transpose << endl
1021 << " value, dx , dy op = " << op << endl
1022 << " just inside = " << inside << endl;
1023 }
1024 using namespace Fem2D;
1025 int n=Uh.NbOfDF;
1026 int mm=Vh.NbOfDF;
1027 if(transpose) Exchange(n,mm);
1028 m = new MatriceMorse<R>(n,mm,0,0);
1029
1030 RdHat1 Gh= RdHat1::diag(1./(RdHat1::d+1));
1031 RdHat2 G;
1032
1033 int n1=n+1;
1034 const Mesh1 & ThU =Uh.Th; // line
1035 const Mesh2 & ThV =Vh.Th; // colunm
1036 bool samemesh = (is_same< Mesh1, Mesh2 >::value) ? (void*)&Uh.Th == (void*)&Vh.Th : 0 ; // same Mesh
1037 int thecolor =0;
1038
1039 int nnz =0;
1040
1041 KN<int> color(ThV.nt);
1042 KN<int> mark(n);
1043 mark=0;
1044
1045 int * cl = 0;
1046 double *a=0;
1047
1048 color=thecolor++;
1049 FElement1 Uh0 = Uh[0];
1050 FElement2 Vh0 = Vh[0];
1051
1052
1053 int nbdfVK= Vh0.NbDoF();
1054 int NVh= Vh0.N;
1055
1056 InterpolationMatrix<RdHat1> ipmat(Uh);
1057
1058 int nbp=ipmat.np; //
1059 KN<RdHat2> PV(nbp); // the PtHat in ThV mesh
1060 KN<int> itV(nbp); // the Triangle number
1061 KN<bool> intV(nbp); // ouside or not
1062
1063 KNM<R> aaa(nbp,nbdfVK);
1064
1065
1066 const R eps = 1.0e-10;
1067 const int sfb1=Vh0.N*last_operatortype*Vh0.NbDoF();
1068 KN<R> kv(sfb1*nbp);
1069 R * v = kv;
1070 KN<int> ik(nbp); // the Triangle number
1071
1072 bool whatd[last_operatortype];
1073 for (int i=0;i<last_operatortype;i++)
1074 whatd[i]=false;
1075 whatd[op]=true; // the value of function
1076
1077 KN<bool> fait(Uh.NbOfDF);
1078 fait=false;
1079 double epsP=0.; // must be choose
1080 {
1081
1082 for (int it=0;it<ThU.nt;it++)
1083 {
1084 thecolor++; // change the current color
1085 const Element1 & TU(ThU[it]);
1086 FElement1 KU(Uh[it]);
1087 ipmat.set(KU);
1088 if (samemesh) {
1089 copyPt<RdHat2,RdHat1>(PV,ipmat.P);
1090 itV = it;
1091 intV= false;// add July 2009 (unset varaible FH)
1092 }
1093 else {
1094 const Element2 *ts=0,*ts0=0;
1095 bool outside;
1096 R3 P1(TU(Gh));
1097 R2 P12(P1.p2());
1098 if(abs(P1.z)>epsP) {outside=true;ts0=0;}
1099 else
1100 ts0=ThV.Find(P12,G,outside,ts0);
1101 if(outside) ts0=0; // bad starting tet
1102 for (int i=0;i<nbp;i++) {
1103 R3 P1(TU(ipmat.P[i]));
1104 R2 P12(P1.p2());
1105 if(abs(P1.z)>epsP
1106
1107
1108
1109
1110
1111 ) {outside=true;ts=0;}
1112 else
1113 ts=ThV.Find(P12,PV[i],outside,ts0);
1114 if( ts0 ==0 && !outside) ts0=ts;
1115 if(outside && verbosity>9 )
1116 cout << it << " " << i << " :: " << TU(ipmat.P[i]) << " -- "<< outside << PV[i] << " " << ThV(ts) << " -> " << (*ts)(PV[i]) <<endl;
1117 itV[i]= ThV(ts);
1118 intV[i]=outside && inside; // ouside and inside flag
1119 }
1120 }
1121
1122 for (int p=0;p<nbp;p++)
1123 {
1124 KNMK_<R> fb(v+p*sfb1,nbdfVK,NVh,last_operatortype); // valeur de fonction de base de Vh
1125 // ou: fb(idf,j,0) valeur de la j composante de la fonction idf
1126 Vh0.tfe->FB(whatd,ThV,ThV[itV[p]],PV[p],fb);
1127 }
1128
1129 for (int i=0;i<ipmat.ncoef;i++)
1130 { // pour tous le terme
1131 int dfu = KU(ipmat.dofe[i]); // le numero de df global
1132 if(fait[dfu]) continue;
1133 int jU = ipmat.comp[i]; // la composante dans U
1134 int p=ipmat.p[i]; // le point
1135 if (intV[p]) continue; // ouside and inside flag => next
1136 R aipj = ipmat.coef[i];
1137 FElement2 KV(Vh[itV[p]]);
1138 int jV=jU;
1139 if(iU2V) jV=iU2V[jU];
1140
1141 if(jV>=0 && jV<NVh)
1142 {
1143 KNMK_<R> fb(v+p*sfb1,nbdfVK,NVh,last_operatortype);
1144 KN_<R> fbj(fb('.',jV,op));
1145
1146 for (int idfv=0;idfv<nbdfVK;idfv++)
1147 if (Abs(fbj[idfv])>eps)
1148 {
1149 int dfv=KV(idfv);
1150 int ii=dfu, jj=dfv;
1151 if(transpose) Exchange(ii,jj);
1152 // le term dfu,dfv existe dans la matrice
1153 R c= fbj[idfv]*aipj;
1154 if(Abs(c)>eps)
1155 (*m)(ii,jj) += c;
1156 }
1157 }
1158
1159 }
1160
1161 for (int df=0;df<KU.NbDoF();df++)
1162 {
1163 int dfu = KU(df); // le numero de df global
1164 fait[dfu]=true;
1165 }
1166
1167
1168 }
1169
1170 }
1171 return m;
1172 }
1173
1174
buildInterpolationMatrix1(const FESpace & Uh,const KN_<double> & xx,const KN_<double> & yy,int * data)1175 MatriceMorse<R> * buildInterpolationMatrix1(const FESpace & Uh,const KN_<double> & xx,const KN_<double> & yy ,int *data)
1176 {
1177 int op=op_id; // value of the function
1178 int icomp=0;
1179 bool transpose=false;
1180 bool inside=false;
1181 if (data)
1182 {
1183 transpose=data[0];
1184 op=data[1];
1185 inside=data[2];
1186 icomp = data[3];
1187 ffassert(op>=0 && op < 4);
1188
1189 }
1190 if(verbosity>2)
1191 {
1192 cout << " -- buildInterpolationMatrix transpose =" << transpose << endl
1193 << " value, dx , dy op = " << op << endl
1194 << " composante = " << icomp << endl
1195 << " just inside = " << inside << endl;
1196 }
1197 using namespace Fem2D;
1198 int n=Uh.NbOfDF;
1199 int mm=xx.N();
1200 int nbxx= mm;
1201 if(transpose) Exchange(n,mm);
1202 const Mesh & ThU =Uh.Th; // line
1203
1204
1205 FElement Uh0 = Uh[0];
1206
1207
1208
1209 int nbdfUK= Uh0.NbDoF();
1210 int NUh= Uh0.N;
1211
1212 ffassert(icomp < NUh && icomp >=0);
1213
1214
1215 const int sfb1=Uh0.N*last_operatortype*Uh0.NbDoF();
1216 KN<R> kv(sfb1);
1217 R * v = kv;
1218 const R eps = 1.0e-10;
1219
1220 bool whatd[last_operatortype];
1221 for (int i=0;i<last_operatortype;i++)
1222 whatd[i]=false;
1223 whatd[op]=true; // the value of function
1224 KN<bool> fait(Uh.NbOfDF);
1225 fait=false;
1226 R2 Phat;
1227 bool outside;
1228 MatriceMorse<R> * m = new MatriceMorse<R>(n,mm,0,0);
1229 for(int ii=0;ii<nbxx;ii++)
1230 {
1231 const Triangle *ts=ThU.Find(R2(xx[ii],yy[ii]),Phat,outside);
1232 if(outside && !inside) continue;
1233 int it = ThU(ts);
1234 FElement KU(Uh[it]);
1235 KNMK_<R> fb(v,nbdfUK,NUh,last_operatortype);
1236 Uh0.tfe->FB(whatd,ThU,ThU[it],Phat,fb);
1237 KN_<R> Fwi(fb('.',icomp,op));
1238 for (int idfu=0;idfu<nbdfUK;idfu++)
1239 {
1240 int j = ii;
1241 int i = KU(idfu);
1242 if(transpose) Exchange(i,j);
1243 R c = Fwi(idfu);
1244 if(Abs(c)>eps)
1245 (*m)(i,j) += c;
1246 }
1247 }
1248
1249
1250 // sij.clear();
1251 return m;
1252 }
1253
1254 template<class FESpaceT>
buildInterpolationMatrixT1(const FESpaceT & Uh,const KN_<double> & xx,const KN_<double> & yy,const KN_<double> & zz,int * data)1255 MatriceMorse<R> * buildInterpolationMatrixT1(const FESpaceT & Uh,const KN_<double> & xx,const KN_<double> & yy ,const KN_<double> & zz,int *data)
1256 {
1257 typedef typename FESpaceT::Mesh MeshT;
1258 typedef typename FESpaceT::FElement FElementT;
1259 typedef typename MeshT::Element ElementT;
1260 typedef typename FESpaceT::Rd RdT;
1261 typedef typename ElementT::RdHat RdHatT;
1262
1263 int op=op_id; // value of the function
1264 int icomp=0;
1265 bool transpose=false;
1266 bool inside=false;
1267 if (data){
1268 transpose=data[0];
1269 op=data[1];
1270 inside=data[2];
1271 icomp = data[3];
1272 ffassert(op>=0 && op < 4);
1273 if(op==3) op=op_dz;// correct missing
1274 }
1275 if(verbosity>2){
1276 cout << " -- buildInterpolationMatrix transpose =" << transpose << endl
1277 << " value, dx , dy op = " << op << endl
1278 << " composante = " << icomp << endl
1279 << " just inside = " << inside << endl;
1280 }
1281 using namespace Fem2D;
1282 int n=Uh.NbOfDF;
1283 int mm=xx.N();
1284 int nbxx= mm;
1285 if(transpose) Exchange(n,mm);
1286 const MeshT & ThU =Uh.Th; // line
1287
1288 FElementT Uh0 = Uh[0];
1289
1290 int nbdfUK= Uh0.NbDoF();
1291 int NUh= Uh0.N;
1292
1293 ffassert(icomp < NUh && icomp >=0);
1294
1295 const int sfb1=Uh0.N*last_operatortype*Uh0.NbDoF();
1296 KN<R> kv(sfb1);
1297 R * v = kv;
1298 const R eps = 1.0e-10;
1299
1300 What_d whatd= 1 <<op;
1301 KN<bool> fait(Uh.NbOfDF);
1302 fait=false;
1303 // map< pair<int,int> , double > sij;
1304 MatriceMorse<R> * m = new MatriceMorse<R>(n,mm,0,0);
1305 RdHatT Phat;
1306 bool outside;
1307
1308 for(int ii=0;ii<nbxx;ii++){
1309 if(verbosity>9) cout << " Find ThU " <<ii << ":" << RdT(xx[ii],yy[ii],zz[ii]) << endl;
1310 const ElementT *ts=ThU.Find(RdT(xx[ii],yy[ii],zz[ii]),Phat,outside);
1311 if(outside && !inside) continue;
1312 int it = ThU(ts);
1313 FElementT KU(Uh[it]);
1314 KNMK_<R> fb(v,nbdfUK,NUh,last_operatortype);
1315 Uh0.tfe->FB(whatd,ThU,ThU[it],Phat,fb);
1316 KN_<R> Fwi(fb('.',icomp,op));
1317 for (int idfu=0;idfu<nbdfUK;idfu++){
1318 int j = ii;
1319 int i = KU(idfu);
1320 if(transpose) Exchange(i,j);
1321 R c = Fwi(idfu);
1322 if(Abs(c)>eps)
1323 (*m)(i,j) += c;
1324 }
1325 }
1326
1327 return m;
1328 }
1329
SetMatrixInterpolation1(Stack stack,Expression emat,Expression einter,int init)1330 AnyType SetMatrixInterpolation1(Stack stack,Expression emat,Expression einter,int init)
1331 {
1332 using namespace Fem2D;
1333
1334 Matrice_Creuse<R> * sparse_mat =GetAny<Matrice_Creuse<R>* >((*emat)(stack));
1335 const MatrixInterpolation<pfes,pfes>::Op * mi(dynamic_cast<const MatrixInterpolation<pfes,pfes>::Op *>(einter));
1336 ffassert(einter);
1337 pfes * pUh = GetAny< pfes * >((* mi->a)(stack));
1338 FESpace * Uh = **pUh;
1339 int NUh =Uh->N;
1340 int* data = new int[4 + NUh];
1341 data[0]=mi->arg(0,stack,false); // transpose not
1342 data[1]=mi->arg(1,stack,(long) op_id); ; // get just value
1343 data[2]=mi->arg(2,stack,false); ; // get just value
1344 data[3]=mi->arg(3,stack,0L); ; // get just value
1345 KN<long> U2Vc;
1346 U2Vc= mi->arg(4,stack,U2Vc); ;
1347 if( mi->c==0)
1348 { // old cas
1349 pfes * pVh = GetAny< pfes * >((* mi->b)(stack));
1350 FESpace * Vh = **pVh;
1351 int NVh =Vh->N;
1352
1353 for(int i=0;i<NUh;++i)
1354 data[4+i]=i;//
1355 for(int i=0;i<min(NUh,(int) U2Vc.size());++i)
1356 data[4+i]= U2Vc[i];//
1357 if(verbosity>3)
1358 for(int i=0;i<NUh;++i)
1359 {
1360 cout << "The Uh componante " << i << " -> " << data[4+i] << " Componante of Vh " <<endl;
1361 }
1362 for(int i=0;i<NUh;++i)
1363 if(data[4+i]>=NVh)
1364 {
1365 cout << "The Uh componante " << i << " -> " << data[4+i] << " >= " << NVh << " number of Vh Componante " <<endl;
1366 ExecError("Interpolation incompability beetween componante ");
1367 }
1368
1369 ffassert(Vh);
1370 ffassert(Uh);
1371
1372 if(!init) sparse_mat->init();
1373 sparse_mat->typemat=0; //TypeSolveMat(TypeSolveMat::NONESQUARE); // none square matrice (morse)
1374 sparse_mat->A.master(buildInterpolationMatrix(*Uh,*Vh,data));
1375 }
1376 else
1377 { // new cas mars 2006
1378 KN_<double> xx = GetAny< KN_<double> >((* mi->b)(stack));
1379 KN_<double> yy = GetAny< KN_<double> >((* mi->c)(stack));
1380 ffassert( xx.N() == yy.N());
1381 ffassert(Uh);
1382
1383 if(!init) sparse_mat->init();
1384 sparse_mat->typemat=0;//TypeSolveMat(TypeSolveMat::NONESQUARE); // none square matrice (morse)
1385 sparse_mat->A.master(buildInterpolationMatrix1(*Uh,xx,yy,data));
1386 }
1387 delete [] data;
1388 return sparse_mat; // Warning .. no correct gestion of temp ptr ..
1389 }
1390
1391 template<class pfesT1, class FESpaceT1, class pfesT2, class FESpaceT2 >
SetMatrixInterpolationT1(Stack stack,Expression emat,Expression einter,int init)1392 AnyType SetMatrixInterpolationT1(Stack stack,Expression emat,Expression einter,int init)
1393 {
1394 using namespace Fem2D;
1395
1396 Matrice_Creuse<R> * sparse_mat =GetAny<Matrice_Creuse<R>* >((*emat)(stack));
1397 const typename MatrixInterpolation<pfesT1,pfesT2>::Op * mi(dynamic_cast<const typename MatrixInterpolation<pfesT1,pfesT2>::Op *>(einter));
1398 ffassert(einter);
1399 pfesT1 * pUh = GetAny< pfesT1 * >((* mi->a)(stack));
1400 FESpaceT1 * Uh = **pUh;
1401 int NUh =Uh->N;
1402
1403 int* data = new int[4 + NUh];
1404 data[0]=mi->arg(0,stack,false); // transpose not
1405 data[1]=mi->arg(1,stack,(long) op_id); ; // get just value
1406 data[2]=mi->arg(2,stack,false); ; // get just value
1407 data[3]=mi->arg(3,stack,0L); ; // get just value
1408 KN<long> U2Vc;
1409 U2Vc= mi->arg(4,stack,U2Vc); ;
1410 if( mi->c==0)
1411 { // old cas
1412 pfesT2 * pVh = GetAny< pfesT2 * >((* mi->b)(stack));
1413 FESpaceT2 * Vh = **pVh;
1414 int NVh =Vh->N;
1415
1416 for(int i=0;i<NUh;++i)
1417 data[4+i]=i;//
1418 for(int i=0;i<min(NUh,(int) U2Vc.size());++i)
1419 data[4+i]= U2Vc[i];//
1420 if(verbosity>3)
1421 for(int i=0;i<NUh;++i)
1422 {
1423 cout << "The Uh componante " << i << " -> " << data[4+i] << " Componante of Vh " <<endl;
1424 }
1425 for(int i=0;i<NUh;++i)
1426 if(data[4+i]>=NVh)
1427 {
1428 cout << "The Uh componante " << i << " -> " << data[4+i] << " >= " << NVh << " number of Vh Componante " <<endl;
1429 ExecError("Interpolation incompability beetween componante ");
1430 }
1431
1432 ffassert(Vh);
1433 ffassert(Uh);
1434 if(!init) sparse_mat->init();
1435 sparse_mat->typemat=0;//(TypeSolveMat::NONESQUARE); // none square matrice (morse)
1436 sparse_mat->A.master(buildInterpolationMatrixT<FESpaceT1,FESpaceT2>(*Uh,*Vh,data)); // sparse_mat->A.master(new MatriceMorse<R>(*Uh,*Vh,buildInterpolationMatrix,data));
1437 }
1438 else
1439 { // new cas mars 2006
1440 KN_<double> xx = GetAny< KN_<double> >((* mi->b)(stack));
1441 KN_<double> yy = GetAny< KN_<double> >((* mi->c)(stack));
1442 KN_<double> zz = GetAny< KN_<double> >((* mi->d)(stack));
1443 ffassert( xx.N() == yy.N());
1444 ffassert( xx.N() == zz.N());
1445 ffassert(Uh);
1446 if(!init) sparse_mat->init();
1447 sparse_mat->typemat=0;//(TypeSolveMat::NONESQUARE); // none square matrice (morse)
1448 sparse_mat->A.master(buildInterpolationMatrixT1<FESpaceT1>(*Uh,xx,yy,zz,data));
1449 }
1450 delete [] data;
1451 return sparse_mat;
1452 }
1453
1454
1455 template<int init>
SetMatrixInterpolation(Stack stack,Expression emat,Expression einter)1456 AnyType SetMatrixInterpolation(Stack stack,Expression emat,Expression einter)
1457 { return SetMatrixInterpolation1(stack,emat,einter,init);}
1458 template<int init>
SetMatrixInterpolation3(Stack stack,Expression emat,Expression einter)1459 AnyType SetMatrixInterpolation3(Stack stack,Expression emat,Expression einter)
1460 { return SetMatrixInterpolationT1<pfes3,FESpace3,pfes3,FESpace3>(stack,emat,einter,init);}
1461 template<int init>
SetMatrixInterpolationS(Stack stack,Expression emat,Expression einter)1462 AnyType SetMatrixInterpolationS(Stack stack,Expression emat,Expression einter)
1463 { return SetMatrixInterpolationT1<pfesS,FESpaceS,pfesS,FESpaceS>(stack,emat,einter,init);}
1464 template<int init>
SetMatrixInterpolationL(Stack stack,Expression emat,Expression einter)1465 AnyType SetMatrixInterpolationL(Stack stack,Expression emat,Expression einter)
1466 { return SetMatrixInterpolationT1<pfesL,FESpaceL,pfesL,FESpaceL>(stack,emat,einter,init);}
1467 template<int init>
SetMatrixInterpolationS3(Stack stack,Expression emat,Expression einter)1468 AnyType SetMatrixInterpolationS3(Stack stack,Expression emat,Expression einter)
1469 { return SetMatrixInterpolationT1<pfesS,FESpaceS,pfes3,FESpace3>(stack,emat,einter,init);}
1470 template<int init>
SetMatrixInterpolationL2(Stack stack,Expression emat,Expression einter)1471 AnyType SetMatrixInterpolationL2(Stack stack,Expression emat,Expression einter)
1472 { return SetMatrixInterpolationT1<pfesL,FESpaceL,pfes,FESpace>(stack,emat,einter,init);}
1473 template<int init>
SetMatrixInterpolationLS(Stack stack,Expression emat,Expression einter)1474 AnyType SetMatrixInterpolationLS(Stack stack,Expression emat,Expression einter)
1475 { return SetMatrixInterpolationT1<pfesL,FESpaceL,pfesS,FESpaceS>(stack,emat,einter,init);}
1476
1477
1478 template<class RA,class RB,class RAB,int init>
ProdMat(Stack stack,Expression emat,Expression prodmat)1479 AnyType ProdMat(Stack stack,Expression emat,Expression prodmat)
1480 {
1481 using namespace Fem2D;
1482
1483 Matrice_Creuse<RAB> * sparse_mat =GetAny<Matrice_Creuse<RA>* >((*emat)(stack));
1484 const Matrix_Prod<RA,RB> AB = GetAny<Matrix_Prod<RA,RB> >((*prodmat)(stack));
1485
1486 MatriceMorse<RA> *mA= AB.A->pHM();
1487 MatriceMorse<RB> *mB= AB.B->pHM();
1488 bool ta = AB.ta, tb = AB.tb;
1489 if( !mA || !mB)
1490 {
1491 ExecError(" numll matrix in pod mat ");
1492 }
1493 int An= mA->n, Am =mA->m;
1494 int Bn= mB->n, Bm =mB->m;
1495 if(ta) swap(An,Am);
1496 if(tb) swap(Bn,Bm);
1497
1498 if( Am!= Bn) {
1499 cerr << " -- Error dim ProdMat A*B : tA =" << AB.ta << " = tB " << AB.tb << endl;
1500 cerr << " --MatProd " << mA->n<< " "<< mA->m << " x " << mB->n<< " "<< mB->m << endl;
1501 ExecError(" Wrong mat dim in MatProd");
1502 }
1503 MatriceMorse<RAB> *mAB=new MatriceMorse<RAB>(An,Bm,0,0);
1504 AddMul(*mAB,*mA,*mB,ta,tb);
1505
1506 if(!init) sparse_mat->init();
1507 sparse_mat->typemat=0;
1508 sparse_mat->A.master(mAB);
1509 return sparse_mat;
1510
1511 }
1512
1513
1514 template<class R,int init>
CombMat(Stack stack,Expression emat,Expression combMat)1515 AnyType CombMat(Stack stack,Expression emat,Expression combMat)
1516 {
1517 using namespace Fem2D;
1518
1519 Matrice_Creuse<R> * sparse_mat =GetAny<Matrice_Creuse<R>* >((*emat)(stack));
1520 list<tuple<R,VirtualMatrix<int,R> *,bool> > * lcB = GetAny<list<tuple<R,VirtualMatrix<int,R> *,bool> >*>((*combMat)(stack));
1521 HashMatrix<int,R> * AA=BuildCombMat<R>(*lcB,false,0,0);
1522
1523 if(!init) sparse_mat->init();
1524 sparse_mat->A.master(AA);
1525 sparse_mat->typemat=0;
1526 delete lcB;
1527 return sparse_mat;
1528 }
1529 template<class R,int cc> // July 2019 FH A += c M + ...
AddCombMat(Stack stack,Expression emat,Expression combMat)1530 AnyType AddCombMat(Stack stack,Expression emat,Expression combMat)
1531 {
1532 using namespace Fem2D;
1533
1534 Matrice_Creuse<R> * pMCA =GetAny<Matrice_Creuse<R>* >((*emat)(stack));
1535 HashMatrix<int,R> * pA=pMCA->pHM();
1536 ffassert(pA);
1537 list<tuple<R,VirtualMatrix<int,R> *,bool> > * lcB = GetAny<list<tuple<R,VirtualMatrix<int,R> *,bool> >*>((*combMat)(stack));
1538 if( cc!=1) Op1_LCMd<R,cc>::f(lcB);
1539 BuildCombMat<R>(*pA,*lcB,false,0,0,false);
1540
1541 delete lcB;
1542 return pMCA;
1543 }
1544
1545 template<class R,int init>
DiagMat(Stack stack,Expression emat,Expression edia)1546 AnyType DiagMat(Stack stack,Expression emat,Expression edia)
1547 {
1548 using namespace Fem2D;
1549 KN<R> * diag=GetAny<KN<R>* >((*edia)(stack));
1550 Matrice_Creuse<R> * sparse_mat =GetAny<Matrice_Creuse<R>* >((*emat)(stack));
1551 if(!init) sparse_mat->init();
1552 sparse_mat->typemat=VirtualMatrix<int,R>::TS_SYM;//TypeSolveMat(TypeSolveMat::GC); // none square matrice (morse)
1553 sparse_mat->A.master(new MatriceMorse<R>((int) diag->N(),(const R*) *diag));
1554 return sparse_mat;
1555 }
1556
1557
1558
1559 template<class Rin,class Rout>
1560 struct ChangeMatriceMorse {
fChangeMatriceMorse1561 static MatriceMorse<Rout> *f(MatriceMorse<Rin> *mr)
1562 {
1563 MatriceMorse<Rout>* mrr=new MatriceMorse<Rout>(*mr);
1564 delete mr;
1565 return mrr;
1566 }
1567 };
1568
1569 template<class R>
1570 struct ChangeMatriceMorse<R,R> {
fChangeMatriceMorse1571 static MatriceMorse<R>* f(MatriceMorse<R>* mr)
1572 {
1573 return mr;
1574 }
1575 };
1576 template<class R,class RR,int init>
CopyMat_tt(Stack stack,Expression emat,Expression eA,bool transp)1577 AnyType CopyMat_tt(Stack stack,Expression emat,Expression eA,bool transp)
1578 {
1579 using namespace Fem2D;
1580 Matrice_Creuse<R> * Mat;
1581
1582 if(transp)
1583 {
1584 Matrice_Creuse_Transpose<R> tMat=GetAny<Matrice_Creuse_Transpose<R> >((*eA)(stack));
1585 Mat=tMat;
1586 }
1587 else Mat =GetAny<Matrice_Creuse<R>*>((*eA)(stack));
1588 MatriceMorse<R> * mr=Mat->pHM();
1589
1590 Matrice_Creuse<RR> * sparse_mat =GetAny<Matrice_Creuse<RR>* >((*emat)(stack));
1591 if(mr) {
1592 MatriceMorse<RR> * mrr = new MatriceMorse<RR>(mr->n,mr->m,0,0);
1593 *mrr = *mr;
1594 if(transp) mrr->dotranspose();
1595
1596
1597 if(!init) sparse_mat->init() ;
1598 sparse_mat->typemat=Mat->typemat; // none square matrice (morse)
1599 sparse_mat->A.master(mrr);
1600 VirtualMatrix<int,RR> *pvm = sparse_mat->pMC();
1601 pvm->SetSolver(); // copy solver ???
1602 }
1603 return sparse_mat;
1604 }
1605
1606 template<class R,class RR,int init>
CopyTrans(Stack stack,Expression emat,Expression eA)1607 AnyType CopyTrans(Stack stack,Expression emat,Expression eA)
1608 {
1609 return CopyMat_tt<R,RR,init>(stack,emat,eA,true);
1610 }
1611 template<class R,class RR,int init>
CopyMat(Stack stack,Expression emat,Expression eA)1612 AnyType CopyMat(Stack stack,Expression emat,Expression eA)
1613 {
1614 return CopyMat_tt<R,RR,init>(stack,emat,eA,false);
1615 }
1616
1617
1618 template<class R,int init>
MatFull2Sparse(Stack stack,Expression emat,Expression eA)1619 AnyType MatFull2Sparse(Stack stack,Expression emat,Expression eA)
1620 {
1621 KNM<R> * A=GetAny<KNM<R>* >((*eA)(stack));
1622 Matrice_Creuse<R> * sparse_mat =GetAny<Matrice_Creuse<R>* >((*emat)(stack));
1623 if(!init) sparse_mat->init() ;
1624 sparse_mat->typemat=0;//(TypeSolveMat::GMRES); // none square matrice (morse)
1625 sparse_mat->A.master(new MatriceMorse<R>((KNM_<R> &)*A,0.0));
1626
1627 return sparse_mat;
1628 }
1629
1630 template<class RR,class AA=RR,class BB=AA>
1631 struct Op2_pair: public binary_function<AA,BB,RR> {
fOp2_pair1632 static RR f(const AA & a,const BB & b)
1633 { return RR( a, b);}
1634 };
1635
1636
1637 template<class R>
get_mat_n(Matrice_Creuse<R> * p)1638 long get_mat_n(Matrice_Creuse<R> * p)
1639 { ffassert(p ) ; return p->A ?p->A->n: 0 ;}
1640
1641 template<class R>
set_mat_COO(Matrice_Creuse<R> * p)1642 bool set_mat_COO(Matrice_Creuse<R> * p)
1643 { ffassert(p ) ;
1644 HashMatrix<int,R> *phm=p->pHM();
1645 if(phm) phm->COO();
1646 return phm;
1647 }
1648 template<class R>
set_mat_CSR(Matrice_Creuse<R> * p)1649 bool set_mat_CSR(Matrice_Creuse<R> * p)
1650 { ffassert(p ) ;
1651 HashMatrix<int,R> *phm=p->pHM();
1652 if(phm) phm->CSR();
1653 return phm;
1654 }
1655 template<class R>
set_mat_CSC(Matrice_Creuse<R> * p)1656 bool set_mat_CSC(Matrice_Creuse<R> * p)
1657 { ffassert(p ) ;
1658 HashMatrix<int,R> *phm=p->pHM();
1659 if(phm) phm->CSC();
1660 return phm;
1661 }
1662
1663 template<class R>
get_mat_m(Matrice_Creuse<R> * p)1664 long get_mat_m(Matrice_Creuse<R> * p)
1665 { ffassert(p ) ; return p->A ?p->A->m: 0 ;}
1666
1667 template<class R>
get_mat_nbcoef(Matrice_Creuse<R> * p)1668 long get_mat_nbcoef(Matrice_Creuse<R> * p)
1669 { ffassert(p ) ; return p->A ?p->A->NbCoef(): 0 ;}
1670 template<class R>
get_mat_half(Matrice_Creuse<R> * p)1671 long get_mat_half(Matrice_Creuse<R> * p)
1672 { return p && p->pHM() ? p->pHM()->half: -1 ;}
1673
1674 template<class R>
get_NM(const list<tuple<R,MatriceCreuse<R> *,bool>> & lM)1675 pair<long,long> get_NM(const list<tuple<R,MatriceCreuse<R> *,bool> > & lM)
1676 {
1677 typedef typename list<tuple<R,MatriceCreuse<R> *,bool> >::const_iterator lconst_iterator;
1678
1679 lconst_iterator begin=lM.begin();
1680 lconst_iterator end=lM.end();
1681 lconst_iterator i;
1682
1683 long n=0,m=0;
1684 for(i=begin;i!=end;i++++)
1685 {
1686 ffassert(get<1>(*i));
1687 MatriceCreuse<R> * M=get<1>(*i);
1688 bool t=get<2>(*i);
1689 int nn= M->n,mm=M->m;
1690 if (t) swap(nn,mm);
1691 if ( n==0) n = nn;
1692 if ( m==0) m = mm;
1693 if (n != 0) ffassert(nn == n);
1694 if (m != 0) ffassert(mm == m);
1695 }
1696 return make_pair(n,m);
1697 }
1698
1699 template<class R>
get_diag(Matrice_Creuse<R> * p,KN<R> * x)1700 long get_diag(Matrice_Creuse<R> * p, KN<R> * x)
1701 { ffassert(p && x ) ; return p->A ?p->A->getdiag(*x): 0 ;}
1702 template<class R>
set_diag(Matrice_Creuse<R> * p,KN<R> * x)1703 long set_diag(Matrice_Creuse<R> * p, KN<R> * x)
1704 { ffassert(p && x ) ; return p->A ?p->A->setdiag(*x): 0 ;}
1705
1706
1707 template<class R>
get_elementp2mc(Matrice_Creuse<R> * const & ac,const long & b,const long & c)1708 R * get_elementp2mc(Matrice_Creuse<R> * const & ac,const long & b,const long & c){
1709 MatriceMorse<R> * a= ac ? ac->pHM() : 0 ;
1710 if( !a || a->n <= b || c<0 || a->m <= c )
1711 { cerr << " Out of bound 0 <=" << b << " < " << (a ? a->n : 0) << ", 0 <= " << c << " < " << (a ? a->m : 0)
1712 << " Matrix type = " << typeid(ac).name() << endl;
1713 cerr << ac << " " << a << endl;
1714 ExecError("Out of bound in operator Matrice_Creuse<R> (,)");}
1715 R * p =a->npij(b,c);
1716 if( !p) { if(verbosity) cerr << "Error: the coef a(" << b << "," << c << ") do'nt exist in sparse matrix "
1717 << " Matrix type = " << typeid(ac).name() << endl;
1718 ExecError("Use of unexisting coef in sparse matrix operator a(i,j) ");}
1719 return p;}
1720
1721 template<class R>
get_element2mc(Matrice_Creuse<R> * const & ac,const long & b,const long & c)1722 R get_element2mc(Matrice_Creuse<R> * const & ac,const long & b,const long & c){
1723 MatriceCreuse<R> * a= ac ? ac->A:0 ;
1724 R r=R();
1725 if( !a || a->n <= b || c<0 || a->m <= c )
1726 { cerr << " Out of bound 0 <=" << b << " < " << a->n << ", 0 <= " << c << " < " << a->m
1727 << " Matrix type = " << typeid(ac).name() << endl;
1728 cerr << ac << " " << a << endl;
1729 ExecError("Out of bound in operator Matrice_Creuse<R> (,)");}
1730 R * p =a->pij(b,c);
1731 if(p) r=*p;
1732 return r;}
1733
1734 template<class RR,class AA=RR,class BB=AA>
1735 struct Op2_mulAv: public binary_function<AA,BB,RR> {
fOp2_mulAv1736 static RR f(const AA & a,const BB & b)
1737 { return (*a->A * *b );}
1738 };
1739
1740 template<class RR,class AA=RR,class BB=AA>
1741 struct Op2_mulvirtAv: public binary_function<AA,BB,RR> {
fOp2_mulvirtAv1742 static RR f(const AA & a,const BB & b)
1743 { return RR( (*a).A, b );}
1744 };
1745
1746 class Matrice_Creuse_C2R { public:
1747 typedef Complex K;
1748 Matrice_Creuse<K> * A;
1749 int cas; // 0 re , 1 im
Matrice_Creuse_C2R(Matrice_Creuse<K> * AA,int cass)1750 Matrice_Creuse_C2R(Matrice_Creuse<K> * AA,int cass) : A(AA),cas(cass) {assert(A);}
operator MatriceCreuse<K>&() const1751 operator MatriceCreuse<K> & () const {return *A->A;}
operator Matrice_Creuse<K>*() const1752 operator Matrice_Creuse<K> * () const {return A;}
1753 };
1754
1755 // ZZZZZ
realC(Complex c)1756 R realC(Complex c) {return c.real();}
imagC(Complex c)1757 R imagC(Complex c) {return c.imag();}
1758
1759
1760 template<int cas>
Build_Matrice_Creuse_C2R(Stack stack,Matrice_Creuse<Complex> * const & Mat)1761 newpMatrice_Creuse<double> Build_Matrice_Creuse_C2R(Stack stack,Matrice_Creuse<Complex> * const & Mat)
1762 {
1763
1764 typedef Complex C;
1765 typedef double R;
1766 using namespace Fem2D;
1767 MatriceMorse<C> * mr= Mat->pHM();
1768 ffassert(mr);
1769 MatriceMorse<R> * mrr = 0;
1770 if(cas==0)
1771 mrr = new MatriceMorse<R>(*mr,realC);
1772 else if(cas==1)
1773 mrr = new MatriceMorse<R>(*mr,imagC);
1774 else {
1775 cout << " cas = " << cas <<endl;
1776 ffassert(0);
1777 }
1778 return newpMatrice_Creuse<double>(stack,mrr);
1779
1780 }
1781
1782 template<class K>
1783 class OneBinaryOperatorA_inv : public OneOperator { public:
OneBinaryOperatorA_inv()1784 OneBinaryOperatorA_inv() : OneOperator(atype<Matrice_Creuse_inv<K> >(),atype<Matrice_Creuse<K> *>(),atype<long>()) {}
code(const basicAC_F0 & args) const1785 E_F0 * code(const basicAC_F0 & args) const
1786 { Expression p=args[1];
1787 if ( ! p->EvaluableWithOutStack() )
1788 {
1789 bool bb=p->EvaluableWithOutStack();
1790 cout << bb << " " << * p << endl;
1791 CompileError(" A^p, The p must be a constant == -1, sorry");}
1792 long pv = GetAny<long>((*p)(NullStack));
1793 if (pv !=-1)
1794 { char buf[100];
1795 sprintf(buf," A^%ld, The pow must be == -1, sorry",pv);
1796 CompileError(buf);}
1797 return new E_F_F0<Matrice_Creuse_inv<K>,Matrice_Creuse<K> *>(Build<Matrice_Creuse_inv<K>,Matrice_Creuse<K> *>,t[0]->CastTo(args[0]));
1798 }
1799 };
1800 template<class K>
1801 class OneBinaryOperatorAt_inv : public OneOperator { public:
OneBinaryOperatorAt_inv()1802 OneBinaryOperatorAt_inv() : OneOperator(atype<Matrice_Creuse_inv_trans<K> >(),atype<Matrice_Creuse_Transpose<K> >(),atype<long>()) {}
code(const basicAC_F0 & args) const1803 E_F0 * code(const basicAC_F0 & args) const
1804 { Expression p=args[1];
1805 if ( ! p->EvaluableWithOutStack() )
1806 {
1807 bool bb=p->EvaluableWithOutStack();
1808 cout << bb << " " << * p << endl;
1809 CompileError(" A^p, The p must be a constant == -1, sorry");}
1810 long pv = GetAny<long>((*p)(NullStack));
1811 if (pv !=-1)
1812 { char buf[100];
1813 sprintf(buf," A^%ld, The pow must be == -1, sorry",pv);
1814 CompileError(buf);}
1815 return new E_F_F0<Matrice_Creuse_inv_trans<K>,Matrice_Creuse_Transpose<K> >(Build<Matrice_Creuse_inv_trans<K>,Matrice_Creuse_Transpose<K> >,t[0]->CastTo(args[0]));
1816 }
1817 };
1818
1819
1820
1821
1822 template<class K>
1823 class Psor : public E_F0 { public:
1824
1825 typedef double Result;
1826 Expression mat;
1827 Expression xx,gmn,gmx,oomega;
Psor(const basicAC_F0 & args)1828 Psor(const basicAC_F0 & args)
1829 {
1830 args.SetNameParam();
1831 mat=to<Matrice_Creuse<K> *>(args[0]);
1832 gmn=to<KN<K>*>(args[1]);
1833 gmx=to<KN<K>*>(args[2]);
1834 xx=to<KN<K>*>(args[3]);
1835 oomega=to<double>(args[4]);
1836
1837 }
typeargs()1838 static ArrayOfaType typeargs() {
1839 return ArrayOfaType( atype<double>(),
1840 atype<Matrice_Creuse<K> *>(),
1841 atype<KN<K>*>(),
1842 atype<KN<K>*>(),
1843 atype<KN<K>*>(),
1844 atype<double>(),false);}
1845
f(const basicAC_F0 & args)1846 static E_F0 * f(const basicAC_F0 & args){ return new Psor(args);}
1847
operator ()(Stack s) const1848 AnyType operator()(Stack s) const {
1849 Matrice_Creuse<K>* A= GetAny<Matrice_Creuse<K>* >( (*mat)(s) );
1850 KN<K>* gmin = GetAny<KN<K>* >( (*gmn)(s) );
1851 KN<K>* gmax = GetAny<KN<K>* >( (*gmx)(s) );
1852 KN<K>* x = GetAny<KN<K>* >( (*xx)(s) );
1853 double omega = GetAny<double>((*oomega)(s));
1854 return A->A->psor(*gmin,*gmax,*x,omega);
1855 }
1856
1857 };
1858 template <class R>
1859 struct TheDiagMat {
1860 Matrice_Creuse<R> * A;
TheDiagMatTheDiagMat1861 TheDiagMat(Matrice_Creuse<R> * AA) :A(AA) {ffassert(A);}
get_mat_daigTheDiagMat1862 void get_mat_daig( KN_<R> & x) { ffassert(A && A->A && x.N() == A->A->n && A->A->n == A->A->m );
1863 A->A->getdiag(x);}
init_get_mat_daigTheDiagMat1864 void init_get_mat_daig( KN<R> & x) {
1865 ffassert(A && A->A && A->A->n == A->A->m );
1866 x.init(A->A->n);
1867 A->A->getdiag(x);}
set_mat_daigTheDiagMat1868 void set_mat_daig(const KN_<R> & x) { ffassert(A && A->A && x.N() == A->A->n && A->A->n == A->A->m );
1869 A->A->setdiag(x);}
1870 };
1871
1872 template <class R>
1873 struct TheCoefMat {
1874 Matrice_Creuse<R> * A;
TheCoefMatTheCoefMat1875 TheCoefMat(Matrice_Creuse<R> * AA) :A(AA) {ffassert(A);}
get_mat_coefTheCoefMat1876 void get_mat_coef( KN_<R> & x) { ffassert(A && A->A && x.N() == A->A->NbCoef() );
1877 A->A->getcoef(x);}
set_mat_coefTheCoefMat1878 void set_mat_coef(const KN_<R> & x) { ffassert(A && A->A && x.N() == A->A->NbCoef() );
1879 A->A->setcoef(x);}
set_mat_coefTheCoefMat1880 void set_mat_coef(const R & v) { ffassert(A && A->A );
1881 KN<R> x(1,0,v);// vertor constant to v
1882 A->A->setcoef(x);}
1883
traceTheCoefMat1884 R trace() { return A->A->trace(); }
1885 };
1886
1887 template<class R>
get_trace_mat(Matrice_Creuse<R> * p)1888 R get_trace_mat(Matrice_Creuse<R> * p)
1889 {
1890 return p ? p->A->trace():0.;
1891
1892 }
1893 template<class R>
clear_mat(Matrice_Creuse<R> * p)1894 bool clear_mat(Matrice_Creuse<R> * p)
1895 {
1896 if(p) p->A->clear();
1897 return true;
1898 }
1899
1900
1901 template<class R>
thediag(Matrice_Creuse<R> * p)1902 TheDiagMat<R> thediag(Matrice_Creuse<R> * p)
1903 { return TheDiagMat<R>(p);}
1904
1905 template<class R>
thecoef(Matrice_Creuse<R> * p)1906 TheCoefMat<R> thecoef(Matrice_Creuse<R> * p)
1907 { return TheCoefMat<R>(p);}
1908
1909 template<class R>
set_mat_daig(TheDiagMat<R> dm,KN<R> * x)1910 TheDiagMat<R> set_mat_daig(TheDiagMat<R> dm,KN<R> * x)
1911 {
1912 dm.set_mat_daig(*x);
1913 return dm;
1914 }
1915 template<class R>
get_mat_daig(KN<R> * x,TheDiagMat<R> dm)1916 KN<R> * get_mat_daig(KN<R> * x,TheDiagMat<R> dm)
1917 {
1918 dm.get_mat_daig(*x);
1919 return x;
1920 }
1921 template<class R>
init_get_mat_daig(KN<R> * x,TheDiagMat<R> dm)1922 KN<R> * init_get_mat_daig(KN<R> * x,TheDiagMat<R> dm)
1923 {
1924 dm.init_get_mat_daig(*x);
1925 return x;
1926 }
1927
1928
1929 template<class R>
set_mat_coef(TheCoefMat<R> dm,KN<R> * x)1930 TheCoefMat<R> set_mat_coef(TheCoefMat<R> dm,KN<R> * x)
1931 {
1932 dm.set_mat_coef(*x);
1933 return dm;
1934 }
1935 template<class R>
set_mat_coef(TheCoefMat<R> dm,R x)1936 TheCoefMat<R> set_mat_coef(TheCoefMat<R> dm,R x)
1937 {
1938 dm.set_mat_coef(x);
1939 return dm;
1940 }
1941 template<class R>
get_mat_coef(KN<R> * x,TheCoefMat<R> dm)1942 KN<R> * get_mat_coef(KN<R> * x,TheCoefMat<R> dm)
1943 {
1944 dm.get_mat_coef(*x);
1945 return x;
1946 }
1947
1948 template<class T> struct Thresholding {
1949 Matrice_Creuse<T> *v;
ThresholdingThresholding1950 Thresholding (Matrice_Creuse<T> *vv): v(vv) {}
1951 };
1952
1953 template<class R>
thresholding2(const Thresholding<R> & t,const double & threshold)1954 Matrice_Creuse<R>*thresholding2 (const Thresholding<R> &t, const double &threshold) {
1955 typedef HashMatrix<int,R> HMat;
1956 Matrice_Creuse<R> *sparse_mat = t.v;
1957 if (sparse_mat) {
1958 HMat *phm=sparse_mat->pHM() ;
1959 if( phm)
1960 {
1961 int n = phm->n, m = phm->m;
1962 int nnzo = phm->nnz;
1963 phm->resize(n,m,0,threshold);
1964 if (verbosity) {cout << " thresholding : remove " << nnzo-phm->nnz << " them in the matrix " << sparse_mat << " " << threshold << endl;}
1965 } else if (verbosity) {cout << " empty matrix " << sparse_mat << endl;}
1966 }
1967
1968 return t.v;
1969 }
1970
1971 template<class T>
symmetrizeCSR(Matrice_Creuse<T> * const & sparse_mat)1972 long symmetrizeCSR (Matrice_Creuse<T> *const &sparse_mat)
1973 {
1974
1975 typedef HashMatrix<int,T> HMat;
1976 if (sparse_mat) {
1977 HMat *phm=sparse_mat->pHM() ;
1978 if( phm)
1979 {
1980 int n = phm->n, m = phm->m;
1981 int nnzo = phm->nnz;
1982 phm->resize(n,m,0,-1,true);
1983 if (verbosity) {cout << " symmetrizeCSR remove " << (long) nnzo-(long) phm->nnz << " them in the matrix " << sparse_mat << endl;}
1984 } else if (verbosity) {cout << " empty matrix " << sparse_mat << endl;}
1985 }
1986
1987 return 1L;
1988 }
1989
1990 template<class T>
to_Thresholding(Matrice_Creuse<T> * v)1991 Thresholding<T> to_Thresholding (Matrice_Creuse<T> *v) {return Thresholding<T>(v);}
1992
1993 template<class R>
IsRawMat(const basicAC_F0 & args)1994 bool IsRawMat(const basicAC_F0 & args)
1995 {
1996
1997 const E_Array * pee= dynamic_cast<const E_Array*>((Expression) args[1]);
1998 if (!pee) return 0;
1999 const E_Array &ee=*pee;
2000 int N=ee.size();
2001 if (N==1)
2002 {
2003 C_F0 c0(ee[0]);
2004 return
2005 atype<KN_<R> >()->CastingFrom(ee[0].left());
2006
2007 }
2008 else if (N==3)
2009 {
2010 C_F0 c0(ee[0]),c1(ee[1]),c2(ee[2]);
2011 return
2012 atype<KN_<long> >()->CastingFrom(ee[0].left())
2013 && atype<KN_<long> >()->CastingFrom(ee[1].left())
2014 && atype<KN_<R> >()->CastingFrom(ee[2].left());
2015
2016 }
2017 return 0;
2018 }
2019
2020
2021 template<typename R>
2022 class RawMatrix : public E_F0 { public:
2023 int init;
2024 typedef Matrice_Creuse<R> * Result;
2025 Expression emat;
2026 Expression coef,col,lig;
2027 RawMatrix(const basicAC_F0 & args,int initt) ;
typeargs()2028 static ArrayOfaType typeargs() { return ArrayOfaType(atype<Matrice_Creuse<R>*>(),atype<E_Array>());}
2029 AnyType operator()(Stack s) const ;
2030 };
2031
2032 template<typename R>
2033 class BlockMatrix : public E_F0 { public:
2034 typedef Matrice_Creuse<R> * Result;
2035 int N,M;
2036 int init;
2037 Expression emat;
2038 Expression ** e_Mij;
2039 int ** t_Mij;
2040 BlockMatrix(const basicAC_F0 & args,int iinit=0) ;
2041 ~BlockMatrix() ;
2042
typeargs()2043 static ArrayOfaType typeargs() { return ArrayOfaType(atype<Matrice_Creuse<R>*>(),atype<E_Array>());}
f(const basicAC_F0 & args)2044 static E_F0 * f(const basicAC_F0 & args){
2045 if(IsRawMat<R>(args)) return new RawMatrix<R>(args,0);
2046 else return new BlockMatrix(args,0);
2047 }
2048 AnyType operator()(Stack s) const ;
2049
2050 };
2051 template<typename R>
2052 class BlockMatrix1 : public BlockMatrix<R> { public:
BlockMatrix1(const basicAC_F0 & args)2053 BlockMatrix1(const basicAC_F0 & args): BlockMatrix<R>(args,1) {}
f(const basicAC_F0 & args)2054 static E_F0 * f(const basicAC_F0 & args){
2055 if(IsRawMat<R>(args)) return new RawMatrix<R>(args,1);
2056 else return new BlockMatrix<R>(args,1);
2057 }
2058
2059 };
2060
2061 template<typename R>
Matrixfull2mapIJ_inv(Stack s,KNM<R> * const & pa,const Inv_KN_long & iii,const Inv_KN_long & jjj)2062 newpMatrice_Creuse<R> Matrixfull2mapIJ_inv (Stack s,KNM<R> * const & pa,const Inv_KN_long & iii,const Inv_KN_long & jjj)
2063 {
2064
2065 const KN_<long> &ii(iii), &jj(jjj);
2066 const KNM<R> & a(*pa);
2067 int N=a.N(),M=a.M();
2068 long n = ii(SubArray(N)).max()+1;
2069 long m= jj(SubArray(M)).max()+1;
2070
2071
2072 HashMatrix<int,R> *pA= new HashMatrix<int,R>((int) n,(int) m,0,0);
2073 HashMatrix<int,R> & A =*pA;
2074
2075 for (long i=0;i<N;++i)
2076 for (long j=0;j<M;++j)
2077 { R aij=a(i,j);
2078 if(ii[i]>=0 && jj[j]>=0 && std::norm(aij)>1e-40)
2079 A(ii[i],jj[j]) += aij;
2080 }
2081
2082 return newpMatrice_Creuse<R>(s,pA);
2083 }
2084
2085 template<typename R>
Matrixfull2mapIJ(Stack s,KNM<R> * const & pa,const KN_<long> & ii,const KN_<long> & jj)2086 newpMatrice_Creuse<R> Matrixfull2mapIJ (Stack s, KNM<R> * const & pa,const KN_<long> & ii,const KN_<long> & jj)
2087 {
2088 const KNM<R> & a(*pa);
2089 int N=a.N(),M=a.M();
2090 long n = ii.N();
2091 long m= jj.N();
2092 HashMatrix<int,R> *pA= new HashMatrix<int,R>((int)n,(int)m,0,0);
2093 HashMatrix<int,R> & A =*pA;
2094
2095 for (long il=0;il<n;++il)// correct juil 2017 FH N --> n
2096 for (long jl=0;jl<m;++jl)// correct juil 2017 FH M --> m
2097 {
2098 long i = ii[il];
2099 long j = jj[jl];
2100 if( i>=0 && j >=0) {
2101 if ( !(0 <= i && i < N && 0 <= j && j < M ) )
2102 {
2103 cerr << " Out of Bound in A(I,J) : " << i << " " << j << " not in " << "[0,"<<N<<"[x[0," << M << "[ \n";
2104 ExecError("Out of Bound Error");
2105 }
2106
2107 R aij=a(i,j);
2108 if (std::norm(aij)>1e-40)
2109 A(il,jl) += aij;
2110 }
2111 }
2112
2113 return newpMatrice_Creuse<R> (s,pA);//;pA;
2114 }
2115
2116 template<class R>
Matrixfull2map(Stack s,const AnyType & pp)2117 AnyType Matrixfull2map (Stack s , const AnyType & pp)
2118 {
2119 const KNM<R> & a(*GetAny<KNM<R> *>(pp));
2120 int N=a.N(),M=a.M();
2121 int n = N;
2122 int m= M;
2123 HashMatrix<int,R> *pA= new HashMatrix<int,R>(n,m,0,0);
2124 HashMatrix<int,R> & A =*pA;
2125
2126 A(n-1,m-1) = R(); // Hack to be sure that the last term existe
2127
2128 for (int i=0;i<N;++i)
2129 for (int j=0;j<M;++j)
2130 { R aij=a(i,j);
2131 if (std::norm(aij)>1e-40)
2132 A(i,j) += aij;
2133 }
2134
2135 return SetAny<newpMatrice_Creuse<R> >(newpMatrice_Creuse<R> (s,pA));//
2136 }
2137
2138
2139 template<class R>
Matrixoutp2mapIJ_inv(Stack s,outProduct_KN_<R> * const & pop,const Inv_KN_long & iii,const Inv_KN_long & jjj)2140 newpMatrice_Creuse<R> Matrixoutp2mapIJ_inv (Stack s,outProduct_KN_<R> * const & pop,const Inv_KN_long & iii,const Inv_KN_long & jjj)
2141 {
2142 const KN_<long> &ii(iii), &jj(jjj);
2143 const outProduct_KN_<R> & op(*pop);
2144 long N=op.a.N(),M=op.b.N();
2145 long n = ii(SubArray(N)).max()+1;
2146 long m= jj(SubArray(M)).max()+1;
2147 HashMatrix<int,R> *pA= new HashMatrix<int,R>((int)n,(int)m,0,0);
2148 HashMatrix<int,R> & A =*pA;
2149
2150 for (int i=0;i<N;++i)
2151 for (int j=0;j<M;++j)
2152 {
2153 R aij=op.a[i]*RNM::conj(op.b[j]);
2154 if(ii[i]>=0 && jj[j]>=0 && std::norm(aij)>1e-40)
2155 A(ii[i],jj[j]) += aij;
2156 }
2157 delete pop;
2158
2159 return newpMatrice_Creuse<R> (s,pA);
2160 }
2161
2162
2163 template<class R>
2164 newpMatrice_Creuse<R>
Matrixmapp2mapIJ1(Stack s,Matrice_Creuse<R> * const & mcB,const Inv_KN_long & iii,const Inv_KN_long & jjj)2165 Matrixmapp2mapIJ1 (Stack s,Matrice_Creuse<R> *const & mcB,const Inv_KN_long & iii,const Inv_KN_long & jjj)
2166 {
2167 const KN_<long> &ii(iii), &jj(jjj);
2168 typedef typename map< pair<int,int>, R>::const_iterator It;
2169
2170 int n=0,m=0;
2171 HashMatrix<int,R> *B = dynamic_cast<HashMatrix<int,R> *>(mcB->pMC());
2172 ffassert(B);
2173
2174 int N=ii.N(),M=jj.N();
2175 int nn = ii.max(), mm= jj.max();
2176 HashMatrix<int,R> *pA= new HashMatrix<int,R>(nn,mm,0,0);
2177 HashMatrix<int,R> & A =*pA;
2178
2179 for (int k=0;k<B->nnz;++k)
2180 {
2181 int il = B->i[k];
2182 int jl = B->j[k];
2183 if ( !( 0 <= il && il < N && 0 <= jl && jl < M ) )
2184 {
2185 cerr << " Out of Bound in (Map)(I,J) : " << il << " " << jl << " not in " << "[0,"<<N<<"[x[0," << M << "[ \n";
2186 ExecError("Out of Bound Error");
2187 }
2188 int i=ii(il);
2189 int j=jj(jl);
2190 n=max(i,n);
2191 m=max(j,m);
2192 R aij = B->aij[k];
2193 if(i >=0 && j>=0)
2194 A(i,j) += aij;
2195 }
2196 // A[make_pair(n,m)] += R(); // Hack to be sure that the last term existe
2197 // delete B;
2198 // FH question resize or not pA to (n,m);
2199 return newpMatrice_Creuse<R> (s,pA);//
2200 }
2201
2202 template<class R>
2203 // map< pair<int,int>, R> *
Matrixmapp2mapIJ(Stack s,Matrice_Creuse<R> * const & mcB,const KN_<long> & ii,const KN_<long> & jj)2204 newpMatrice_Creuse<R> Matrixmapp2mapIJ (Stack s,Matrice_Creuse<R> *const & mcB,const KN_<long> & ii,const KN_<long> & jj)
2205 {
2206
2207
2208 typedef typename map< pair<int,int>, R>::const_iterator It;
2209 typedef typename multimap< int,int>::iterator MI;
2210 HashMatrix<int,R> *B = dynamic_cast<HashMatrix<int,R> *>(mcB->pMC());
2211 ffassert(B);
2212 multimap< int,int > I,J;
2213 int N=ii.N(),M=jj.N();
2214 for (int i=0;i<N;++i)
2215 if(ii[i]>=0)
2216 I.insert(make_pair(ii[i],i));
2217 for (int j=0;j<M;++j)
2218 if(jj[j]>=0)
2219 J.insert(make_pair(jj[j],j));
2220 int n=N-1,m=M-1;// change FH sep 2009 to have the correct size..
2221 HashMatrix<int,R> *pA= new HashMatrix<int,R>(N,M,0,0);
2222 HashMatrix<int,R> & A =*pA;
2223
2224 for (int k=0;k!=B->nnz;++k)
2225 {
2226 int il = B->i[k];
2227 int jl = B->j[k];
2228 R aij = B->aij[k];
2229 pair<MI,MI> PPi=I.equal_range(il);
2230 pair<MI,MI> PPj=J.equal_range(jl);
2231 for(MI pi=PPi.first ; pi != PPi.second; ++pi)
2232 {
2233 int i=pi->second;
2234 for(MI pj=PPj.first ; pj != PPj.second; ++pj)
2235 {
2236 int j=pj->second;
2237 n=max(i,n);
2238 m=max(j,m);
2239 if(i >=0 && j>=0)
2240 A(i,j) += aij;
2241 }
2242 }
2243 }
2244 // resier A to n,m ?????
2245 // A[make_pair(n,m)] += R(); // Hack to be sure that the last term existe
2246 // delete B;
2247
2248 // return pA;
2249 return newpMatrice_Creuse<R> (s,pA);//
2250 }
2251
2252 template<class R>
2253 newpMatrice_Creuse<R>
Matrixoutp2mapIJ(Stack s,outProduct_KN_<R> * const & pop,const KN_<long> & ii,const KN_<long> & jj)2254 Matrixoutp2mapIJ (Stack s,outProduct_KN_<R> * const & pop,const KN_<long> & ii,const KN_<long> & jj)
2255 {
2256 const outProduct_KN_<R> & op(*pop);
2257 long N=op.a.N(),M=op.b.N();
2258 long n=ii.N(),m=jj.N();
2259 HashMatrix<int,R> *pA= new HashMatrix<int,R>((int)n,(int)m,0,0);
2260 HashMatrix<int,R> & A =*pA;
2261
2262 for (long il=0;il<n;++il)
2263 for (long jl=0;jl<m;++jl)
2264 {
2265 long i = ii[il];
2266 long j = jj[jl];
2267 if(i>=0 && j >=0)
2268 {
2269 if ( !( 0 <= i && i < N && 0 <= j && j < M ) )
2270 {
2271 cerr << " Out of Bound in (a*b')(I,J) : " << i << " " << j << " not in " << "[0,"<<N<<"[x[0," << M << "[ \n";
2272 ExecError("Out of Bound Error");
2273 }
2274 R aij=op.a[i]*RNM::conj(op.b[j]);
2275 if (std::norm(aij)>1e-40)
2276 A(il,jl) += aij;
2277 }
2278 }
2279 delete pop;
2280 return newpMatrice_Creuse<R> (s,pA);
2281 }
2282
2283
2284 template<class R>
Matrixoutp2map(Stack s,const AnyType & pp)2285 AnyType Matrixoutp2map (Stack s, const AnyType & pp)
2286 {
2287 const outProduct_KN_<R> & op(*GetAny<outProduct_KN_<R> *>(pp));
2288 long N=op.a.N(),M=op.b.N();
2289 long n = N;
2290 long m= M;
2291 // A[make_pair(n-1,m-1)] = R(); // Hack to be sure that the last term existe
2292 HashMatrix<int,R> *pA= new HashMatrix<int,R>((int)n,(int)m,0,0);
2293 HashMatrix<int,R> & A =*pA;
2294
2295 for (long i=0;i<N;++i)
2296 for (long j=0;j<M;++j)
2297 {
2298 R aij=op.a[i]*RNM::conj(op.b[j]);
2299 if (std::norm(aij)>1e-40)
2300 A(i,j) += aij;
2301 }
2302 delete &op;
2303 return SetAny<newpMatrice_Creuse<R>>(newpMatrice_Creuse<R> (s,pA));//
2304 }
2305
2306
~BlockMatrix()2307 template<typename R> BlockMatrix<R>::~BlockMatrix()
2308 {
2309 if (e_Mij)
2310 { if(verbosity>9999) cout << " del Block matrix "<< this << " " << e_Mij <<" N = " << N << " M = " << M << endl;
2311 for (int i=0;i<N;i++)
2312 { delete [] e_Mij[i];
2313 delete [] t_Mij[i];
2314 }
2315 delete [] e_Mij;
2316 delete [] t_Mij;
2317 N=0;
2318 M=0;
2319 e_Mij=0;
2320 t_Mij=0; }
2321 }
2322
RawMatrix(const basicAC_F0 & args,int iinit)2323 template<typename R> RawMatrix<R>::RawMatrix(const basicAC_F0 & args,int iinit)
2324 : init(iinit)
2325 {
2326 args.SetNameParam();
2327 emat = args[0];
2328
2329 const E_Array & ee= *dynamic_cast<const E_Array*>((Expression) args[1]);
2330
2331 int N=ee.size();
2332 if (N==1)
2333 {
2334 C_F0 c0(ee[0]);
2335 coef=to<KN_<R> >(ee[0]);
2336 lig=0;
2337 col=0;
2338 }
2339
2340 else if (N==3)
2341 {
2342 C_F0 c0(ee[0]),c1(ee[1]),c2(ee[2]);
2343 coef=to<KN_<R> >(ee[2]);
2344 lig=to<KN_<long> >(ee[0]);
2345 col=to<KN_<long> >(ee[1]);
2346
2347 }
2348
2349
2350 }
BlockMatrix(const basicAC_F0 & args,int iinit)2351 template<typename R> BlockMatrix<R>::BlockMatrix(const basicAC_F0 & args,int iinit)
2352 : init(iinit)
2353 {
2354 N=0;
2355 M=0;
2356 args.SetNameParam();
2357 emat = args[0];
2358 const E_Array & eMij= *dynamic_cast<const E_Array*>((Expression) args[1]);
2359 N=eMij.size();
2360 int err =0;
2361 for (int i=0;i<N;i++)
2362 {
2363 const E_Array* emi= dynamic_cast<const E_Array*>((Expression) eMij[i]);
2364 if (!emi) err++;
2365 else
2366 {
2367 if ( i==0)
2368 M = emi->size();
2369 else
2370 if(M != emi->size()) err++;
2371 }
2372 }
2373 if (err) {
2374 CompileError(" Block matrix : [[ a, b, c], [ a,b,c ]] or Raw Matrix [a] or [ l, c, a ] ");
2375 }
2376 assert(N && M);
2377 e_Mij = new Expression * [N];
2378 t_Mij = new int * [N];
2379 for (int i=0;i<N;i++)
2380 {
2381 const E_Array li= *dynamic_cast<const E_Array*>((Expression) eMij[i]);
2382
2383 e_Mij[i] = new Expression [M];
2384 t_Mij[i] = new int [M];
2385 for (int j=0; j<M;j++)
2386 {
2387 C_F0 c_Mij(li[j]);
2388 Expression eij=c_Mij.LeftValue();
2389 aType rij = c_Mij.left();
2390 if ( rij == atype<long>() && eij->EvaluableWithOutStack() )
2391 {
2392 long contm = GetAny<long>((*eij)(NullStack));
2393 if(contm==0)
2394 {
2395 e_Mij[i][j]=0;
2396 t_Mij[i][j]=0;
2397 }
2398 else if ( atype<R >()->CastingFrom(rij) )
2399 { // frev 2007
2400 e_Mij[i][j]=to<R>(c_Mij);
2401 t_Mij[i][j]=7; // just un scalaire
2402 }
2403 else CompileError(" Block matrix , Just 0 matrix");
2404 }
2405 else if ( rij == atype<Matrice_Creuse<R> *>())
2406 {
2407 e_Mij[i][j]=eij;
2408 t_Mij[i][j]=1;
2409 }
2410 else if ( rij == atype<Matrice_Creuse_Transpose<R> >())
2411 {
2412 e_Mij[i][j]=eij;
2413 t_Mij[i][j]=2;
2414 }
2415 else if ( atype<KNM<R> * >()->CastingFrom(rij) )
2416 { // before KN_ because KNM can be cast in KN_
2417
2418 e_Mij[i][j]=to<KNM<R> * >(c_Mij);
2419 t_Mij[i][j]=5;
2420 }
2421 else if ( atype<KN_<R> >()->CastingFrom(rij) )
2422 {
2423 e_Mij[i][j]=to<KN_<R> >(c_Mij);
2424 t_Mij[i][j]=3;
2425
2426 }
2427 else if ( atype<Transpose<KN_<R> > >()->CastingFrom(rij) )
2428 {
2429
2430 e_Mij[i][j]=to<Transpose<KN_<R> > >(c_Mij);
2431 t_Mij[i][j]=4;
2432 }
2433 else if ( atype<Transpose< KNM<R> * > >()->CastingFrom(rij) )
2434 {
2435
2436 e_Mij[i][j]=to<Transpose<KNM<R> *> >(c_Mij);
2437 t_Mij[i][j]=6;
2438 }
2439 else if ( atype<R >()->CastingFrom(rij) )
2440 { // frev 2007
2441 e_Mij[i][j]=to<R>(c_Mij);
2442 t_Mij[i][j]=7; // just un scalaire
2443 }
2444
2445 else {
2446
2447 CompileError(" Block matrix , bad type in block matrix");
2448 }
2449 }
2450
2451 }
2452 }
2453
2454 template<typename RR>
2455 class SetRawMatformMat : public OneOperator {
2456 public:
2457 typedef Matrice_Creuse<RR> * A; // Warning B type of 2 parameter
2458 typedef Matrice_Creuse<RR> * R;
2459 typedef E_Array B; // A type of 1 parameter
2460
2461 class CODE : public E_F0 { public:
2462 Expression Mat;
2463 Expression lig;
2464 Expression col;
2465 Expression coef;
2466 bool mi;
CODE(Expression a,const E_Array & tt)2467 CODE(Expression a,const E_Array & tt)
2468 : Mat(a),
2469 mi(tt.MeshIndependent())
2470 {
2471
2472 assert(&tt);
2473 if(tt.size()!=3)
2474 CompileError("Set raw matrix: [ lg,col, a] = A (size !=3) ");
2475 if ( aatypeknlongp->CastingFrom(tt[0].left() ) //// for compilation error with g++ 3.2.2 (4 times)
2476 && aatypeknlongp->CastingFrom(tt[1].left() )
2477 && atype<KN<RR>* >()->CastingFrom(tt[2].left() ) )
2478 {
2479 lig = aatypeknlongp->CastTo(tt[0]);
2480 col = aatypeknlongp->CastTo(tt[1]);
2481 coef = atype<KN<RR>* >()->CastTo(tt[2]);
2482 }
2483 else
2484 CompileError(" we are waiting for [ lg,col,a] = A");
2485 }
2486
operator ()(Stack stack) const2487 AnyType operator()(Stack stack) const
2488 {
2489
2490 //V4
2491 A a=GetAny<A>((*Mat)(stack));
2492
2493 KN<long> *lg,*cl;
2494 KN<RR> *cc;
2495 lg = GetAny<KN<long>*>((*lig)(stack));
2496 cl = GetAny<KN<long>*>((*col)(stack));
2497 cc = GetAny<KN<RR>*>((*coef)(stack));
2498 int n=a->N(),m=a->M();
2499 HashMatrix<int,RR> *mh = a->pHM();
2500 mh->COO();
2501
2502 int kk = mh->nnz,k1=0;
2503 if( mh->pij(n-1,m-1)==0 ) k1=1;
2504 lg->resize(kk+k1);
2505 cc->resize(kk+k1);
2506 cl->resize(kk+k1);
2507 int k=0;
2508
2509 for ( k=0 ; k < kk;++k)
2510 {
2511 (*lg)[k]= mh->i[k];
2512 (*cl)[k]= mh->j[k];
2513 (*cc)[k]= mh->aij[k];
2514 }
2515
2516 if(k1)
2517 {
2518 (*lg)[kk]= n;
2519 (*cl)[kk]= m;
2520 (*cc)[kk]= 0;
2521 }
2522
2523 return SetAny<R>(a);
2524 }
MeshIndependent() const2525 bool MeshIndependent() const {return mi;} //
~CODE()2526 ~CODE() {}
operator aType() const2527 operator aType () const { return atype<R>();}
2528 }; // end sub class CODE
2529
2530
2531 public: // warning hack A and B
code(const basicAC_F0 & args) const2532 E_F0 * code(const basicAC_F0 & args) const
2533 { return new CODE(t[1]->CastTo(args[1]),*dynamic_cast<const E_Array*>( t[0]->CastTo(args[0]).RightValue()));}
SetRawMatformMat()2534 SetRawMatformMat(): OneOperator(atype<R>(),atype<B>(),atype<A>()) {} // warning with A and B
2535
2536 };
2537
operator ()(Stack stack) const2538 template<typename R> AnyType RawMatrix<R>::operator()(Stack stack) const
2539 {
2540 MatriceMorse<R> * amorse =0;
2541 KN_<R> cc(GetAny< KN_<R> >((*coef)(stack)));
2542 int k= cc.N();
2543 int n= k;
2544 int m=n;
2545 bool sym=false;
2546 if( lig && col)
2547 {
2548 KN_<long> lg(GetAny< KN_<long> >((*lig)(stack)));
2549 KN_<long> cl=(GetAny< KN_<long> >((*col)(stack)));
2550 n = lg.max()+1;
2551 m = cl.max()+1;
2552 ffassert( lg.N()==k && cl.N()==k && lg.min()>=0 && lg.max()>=0);
2553 amorse = new MatriceMorse<R>(n,m,k,0);
2554 sym=false;
2555 for(int i=0;i<k;++i)
2556 (*amorse)[make_pair<int,int>((int)lg[i],(int)cl[i])]+=cc[i];
2557 }
2558 else
2559 {
2560 amorse = new MatriceMorse<R>(n,cc);
2561 }
2562
2563 if(verbosity)
2564 cout << " -- Raw Matrix nxm =" <<n<< "x" << m << " nb none zero coef. " << amorse->nnz << endl;
2565
2566 Matrice_Creuse<R> * sparse_mat =GetAny<Matrice_Creuse<R>* >((*emat)(stack));
2567 if( !init) sparse_mat->init();
2568
2569 sparse_mat->A.master(amorse);
2570 sparse_mat->typemat=0; //(amorse->n == amorse->m) ? TypeSolveMat(TypeSolveMat::GMRES) : TypeSolveMat(TypeSolveMat::NONESQUARE); // none square matrice (morse)
2571
2572 if(verbosity>3) { cout << " End Raw Matrix : " << endl;}
2573
2574 return sparse_mat;
2575 }
operator ()(Stack s) const2576 template<typename R> AnyType BlockMatrix<R>::operator()(Stack s) const
2577 {
2578 typedef list<tuple<R,MatriceCreuse<R> *,bool> > * L;
2579 KNM<L> Bij(N,M);
2580 KNM<KNM_<R> * > Fij(N,M);
2581 KNM<bool> cnjij(N,M);
2582 KNM<R> Rij(N,M); // to sto
2583
2584 cnjij = false;
2585 KN<long> Oi(N+1), Oj(M+1);
2586 if(verbosity>9) { cout << " Build Block Matrix : " << N << " x " << M << endl;}
2587 Bij = (L) 0;
2588 Oi = (long) 0;
2589 Oj = (long)0;
2590 for (int i=0;i<N;++i)
2591 for (int j=0;j<M;++j)
2592 {
2593 Fij(i,j)=0;
2594 Expression eij = e_Mij[i][j];
2595 int tij=t_Mij[i][j];
2596 if (eij)
2597 {
2598 cnjij(i,j) = tij%2 == 0;
2599 AnyType e=(*eij)(s);
2600 if (tij==1) Bij(i,j) = to( GetAny< Matrice_Creuse<R>* >( e)) ;
2601 else if (tij==2) Bij(i,j) = to( GetAny<Matrice_Creuse_Transpose<R> >(e));
2602 else if (tij==3) { KN_<R> x=GetAny< KN_<R> >( e); Fij(i,j) = new KNM_<R>(x,x.N(),1);}
2603 else if (tij==4) { KN_<R> x=GetAny< Transpose< KN_<R> > >( e).t ; Fij(i,j) = new KNM_<R>(x,1,x.N());}
2604 else if (tij==5) { KNM<R> * m= GetAny< KNM<R>* >( e); Fij(i,j) = new KNM_<R>(*m);}
2605 else if (tij==6) { KNM<R> * m= GetAny< Transpose< KNM<R>* > >( e).t; Fij(i,j) = new KNM_<R>(m->t()); }
2606 else if (tij==7) { Rij(i,j)=GetAny< R >( e); Fij(i,j) = new KNM_<R>(&(Rij(i,j)),1,1);}
2607
2608 else {
2609 cout << " Bug " << tij << endl;
2610 ExecError(" Type sub matrix block unknown ");
2611 }
2612 }
2613 }
2614 // compute size of matrix
2615 int err=0;
2616 for (int i=0;i<N;++i)
2617 for (int j=0;j<M;++j)
2618 {
2619 pair<long,long> nm(0,0);
2620
2621 if (Bij(i,j))
2622 nm = get_NM( *Bij(i,j));
2623 else if(Fij(i,j))
2624 nm = make_pair<long,long>(Fij(i,j)->N(), Fij(i,j)->M());
2625
2626 if (( nm.first || nm.second) && verbosity>3)
2627 cout << " Block [ " << i << "," << j << " ] = " << nm.first << " x " << nm.second << " cnj = " << cnjij(i,j) << endl;
2628 if (nm.first)
2629 {
2630 if ( Oi(i+1) ==0 ) Oi(i+1)=nm.first;
2631 else if(Oi(i+1) != nm.first)
2632 {
2633 err++;
2634 cerr <<"Error Block Matrix, size sub matrix" << i << ","<< j << " n (old) " << Oi(i+1)
2635 << " n (new) " << nm.first << endl;
2636
2637 }
2638 }
2639 if(nm.second)
2640 {
2641 if ( Oj(j+1) ==0) Oj(j+1)=nm.second;
2642 else if(Oj(j+1) != nm.second)
2643 {
2644 cerr <<"Error Block Matrix, size sub matrix" << i << ","<< j << " m (old) " << Oj(j+1)
2645 << " m (new) " << nm.second << endl;
2646 err++;}
2647 }
2648 }
2649
2650 if (err) ExecError("Error Block Matrix, size sub matrix");
2651 // gestion of zero block ????
2652
2653 for (int j=0;j<M;++j)
2654 { if(verbosity>99) cout << j << " colum size" << Oj(j+1) << endl;
2655 if ( Oj(j+1) ==0) {
2656 Oj(j+1)=1;
2657 if( Oj(j+1) !=1) err++;}
2658 }
2659 for (int i=0;i<N;++i)
2660 {
2661 if(verbosity>99) cout << i << " row size" << Oi(i+1) << endl;
2662 if ( Oi(i+1) ==0) {
2663 Oi(i+1)=1;
2664 if( Oi(i+1) !=1) err++;}
2665 }
2666 if (err) ExecError("Error Block Matrix with 0 line or 0 colomn..");
2667
2668 for (int i=0;i<N;++i)
2669 Oi(i+1) += Oi(i);
2670 for (int j=0;j<M;++j) // correct 10/01/2007 FH
2671 Oj(j+1) += Oj(j);// correct 07/03/2010 FH
2672 long n=Oi(N),m=Oj(M);
2673 if(verbosity>99)
2674 {
2675 cout << " Oi = " << Oi << endl;
2676 cout << " Oj = " << Oj << endl;
2677 }
2678 MatriceMorse<R> * amorse =0;
2679 {
2680 HashMatrix<int,R> *Aij = new HashMatrix<int,R>( n, m,0,0);
2681 for (int i=0;i<N;++i)
2682 for (int j=0;j<M;++j)
2683 if (Bij(i,j))
2684 {
2685 if(verbosity>99)
2686 cout << " Add Block S " << i << "," << j << " = at " << Oi(i) << " x " << Oj(j) << " conj = " << cnjij(i,j) << endl;
2687 HashMatrix<int,R> & mmij=*Aij;
2688 const list<tuple<R,MatriceCreuse<R>*,bool> > &lM=*Bij(i,j);
2689 bool ttrans=false;
2690 int ii00=Oi(i);
2691 int jj00=Oj(j);
2692 bool cnj=cnjij(i,j);
2693 BuildCombMat(mmij,lM,ttrans,ii00,jj00,cnj);
2694
2695
2696 }
2697 else if (Fij(i,j))
2698 {
2699 if(verbosity>99)
2700 cout << " Add Block F " << i << "," << j << " = at " << Oi(i) << " x " << Oj(j) << endl;
2701 BuildCombMat(*Aij,*Fij(i,j),Oi(i),Oj(j),R(1.),cnjij(i,j));// BuildCombMat
2702 }
2703
2704
2705 amorse= Aij;
2706 }
2707 if(verbosity>9)
2708 cout << " -- Block Matrix NxM = " << N << "x" << M << " nxm =" <<n<< "x" << m << " nb none zero coef. " << amorse->nnz << endl;
2709
2710 Matrice_Creuse<R> * sparse_mat =GetAny<Matrice_Creuse<R>* >((*emat)(s));
2711 if(!init) sparse_mat->init();
2712 sparse_mat->A.master(amorse);
2713 sparse_mat->typemat=0;//(amorse->n == amorse->m) ? TypeSolveMat(TypeSolveMat::GMRES) : TypeSolveMat(TypeSolveMat::NONESQUARE); // none square matrice (morse)
2714
2715
2716 // cleanning
2717 for (int i=0;i<N;++i)
2718 for (int j=0;j<M;++j)
2719 if(Bij(i,j)) delete Bij(i,j);
2720 else if(Fij(i,j)) delete Fij(i,j);
2721 if(verbosity>9) { cout << " End Build Blok Matrix : " << endl;}
2722
2723 return sparse_mat;
2724
2725 }
2726
2727 template<class R>
2728 class minusMat { public:
2729 list<tuple<R,MatriceCreuse<R> *,bool> > *l;
minusMat(list<tuple<R,MatriceCreuse<R> *,bool>> * ll)2730 minusMat(list<tuple<R,MatriceCreuse<R> *,bool> > *ll):
2731 l(new list<tuple<R,MatriceCreuse<R> *,bool> >(*ll) )
2732 {
2733 typedef typename list<tuple<R,MatriceCreuse<R> *,bool> >::iterator lci;
2734 for (lci i= l->begin();i !=l->end();++i)
2735 get<0>(*i) *= R(-1);
2736 }
2737 };
2738
2739 template<class R>
mM2L3(Stack,const AnyType & pp)2740 AnyType mM2L3 (Stack , const AnyType & pp)
2741 {
2742 minusMat<R> mpp(to(GetAny<Matrice_Creuse<R> *>(pp)));
2743 return SetAny<minusMat<R> >(mpp);
2744 }
2745
2746 template<class R>
2747 class E_ForAllLoopMatrix
2748 { public:
2749
2750 typedef R *VV;
2751 typedef long KK;
2752
2753 typedef Matrice_Creuse<R> * Tab;
2754
2755 typedef ForAllLoopOpBase DataL;
2756 const DataL *data;
E_ForAllLoopMatrix(const DataL * t)2757 E_ForAllLoopMatrix(const DataL *t): data(t){}
f(Stack s) const2758 AnyType f(Stack s) const {
2759 Tab t= GetAny<Tab >(data->tab(s));
2760 KK * i = GetAny<KK*>(data->i(s));
2761 KK * j = GetAny<KK*>(data->j(s));
2762 VV v = GetAny<VV >(data->v(s));
2763 if(verbosity>1000) {
2764 cout << " i " << (char*) (void *) i - (char*)(void*) s ;
2765 cout << " j " << (char*) (void *) j - (char*)(void*) s ;
2766 cout << " vi " << (char*) (void *) v - (char*)(void*) s ;
2767 cout << endl;
2768 }
2769
2770 ffassert(i && v);
2771 MatriceCreuse<R> *m=t->A;
2772 MatriceMorse<R> *mm = dynamic_cast<MatriceMorse<R>*>(m);
2773 if(!mm) ExecError(" Matrix sparse of bad type ( not HMatrix ) , sorry.. ");
2774 if(mm)
2775 for (long kk=0;kk< mm->nnz; ++kk)
2776 {
2777 *i=mm->i[kk];
2778 *j= mm->j[kk];
2779 *v = mm->aij[kk];
2780 data->code(s);
2781 mm->aij[kk] = *v;
2782 }
2783 return Nothing ;
2784 }
2785
2786 };
2787
2788 // Mat real -> mat complex .... ??? FH. april 2016 ....
2789 struct VirtualMatCR :public RNM_VirtualMatrix<Complex>{ public:
2790 RNM_VirtualMatrix<double>& VM;
2791 typedef Complex R;
VirtualMatCRVirtualMatCR2792 VirtualMatCR( RNM_VirtualMatrix<double> & MM): RNM_VirtualMatrix<Complex>(MM.N,MM.M), VM(MM) {}
addMatMulVirtualMatCR2793 void addMatMul(const KN_<R> & cx, KN_<R> & cy) const {
2794 double *px = static_cast<double*>(static_cast<void*>(cx));
2795 double *py = static_cast<double*>(static_cast<void*>(cy));
2796 KN_<double> rx(px+0,cx.N(),cx.step*2);
2797 KN_<double> ix(px+1,cx.N(),cx.step*2);
2798 KN_<double> ry(py+0,cy.N(),cy.step*2);
2799 KN_<double> iy(py+1,cy.N(),cy.step*2);
2800 VM.addMatMul(rx,ry);
2801 VM.addMatMul(ix,iy);
2802 }
addMatTransMulVirtualMatCR2803 void addMatTransMul(const KN_<R> & cx , KN_<R> & cy ) const {
2804 double *px = static_cast<double*>(static_cast<void*>(cx));
2805 double *py = static_cast<double*>(static_cast<void*>(cy));
2806 KN_<double> rx(px+0,cx.N(),cx.step*2);
2807 KN_<double> ix(px+1,cx.N(),cx.step*2);
2808 KN_<double> ry(py+0,cy.N(),cy.step*2);
2809 KN_<double> iy(py+1,cy.N(),cy.step*2);
2810 VM.addMatTransMul(rx,ry);
2811 VM.addMatTransMul(ix,iy);
2812
2813 }
WithSolverVirtualMatCR2814 bool WithSolver() const {return VM.WithSolver();} // by default no solver
SolveVirtualMatCR2815 virtual void Solve( KN_<R> & cx ,const KN_<R> & cy) const
2816 { if( !VM.WithSolver()) InternalError("RNM_VirtualMatrix::solve not implemented ");
2817 double *px = static_cast<double*>(static_cast<void*>(cx));
2818 double *py = static_cast<double*>(static_cast<void*>(cy));
2819 KN_<double> rx(px+0,cx.N(),cx.step*2);
2820 KN_<double> ix(px+1,cx.N(),cx.step*2);
2821 KN_<double> ry(py+0,cy.N(),cy.step*2);
2822 KN_<double> iy(py+1,cy.N(),cy.step*2);
2823 VM.Solve(rx,ry);
2824 VM.Solve(ix,iy);
2825 }
2826
ChecknbLineVirtualMatCR2827 bool ChecknbLine (int n) const { return VM.ChecknbLine(n); }
ChecknbColumnVirtualMatCR2828 bool ChecknbColumn (int m) const { return VM.ChecknbColumn(m); }
2829 };
2830
2831 template<class R,class A,class B> // extend (4th arg.)
2832 class Op2_mulvirtAvCR : public OneOperator { //
2833 aType r; // return type
2834
2835
2836 public:
2837 class CODE :public E_F0 { public: // extend
2838 Expression a0,a1; // extend
CODE(Expression aa0,Expression aa1)2839 CODE( Expression aa0,Expression aa1) : a0(aa0), a1(aa1) {} // extend (2th arg.)
operator ()(Stack s) const2840 AnyType operator()(Stack s) const
2841 {
2842
2843 RNM_VirtualMatrix<Complex> *pv = new VirtualMatCR ((*GetAny<A>((*a0)(s))).A);
2844 Add2StackOfPtr2Free(s,pv);
2845 return SetAny<R>(R(pv,GetAny<B>((*a1)(s))));
2846 }
nbitem() const2847 virtual size_t nbitem() const {return a1->nbitem(); } // modif ???
MeshIndependent() const2848 bool MeshIndependent() const {return a0->MeshIndependent() && a1->MeshIndependent() ;}
2849 };
2850
code(const basicAC_F0 & args) const2851 E_F0 * code(const basicAC_F0 & args) const
2852 { if ( args.named_parameter && !args.named_parameter->empty() )
2853 CompileError( " They are used Named parameter ");
2854
2855 return new CODE( t[0]->CastTo(args[0]),
2856 t[1]->CastTo(args[1]));} // extend
Op2_mulvirtAvCR(int preff=0)2857 Op2_mulvirtAvCR(int preff=0): // 3->4
2858 OneOperator(map_type[typeid(R).name()],
2859 map_type[typeid(A).name()],
2860 map_type[typeid(B).name()]) {pref=preff;}
2861
2862 };
2863 // Norme Hashmatrix
2864 template<class K>
get_norme_linfty(Matrice_Creuse<K> * p)2865 double get_norme_linfty(Matrice_Creuse<K> * p){
2866 if(p==0) return 0.;
2867 HashMatrix<int,K>* ph=p->pHM();
2868 ffassert(ph);
2869 return ph->norminfty();}
2870 template<class K>
get_norme_l2(Matrice_Creuse<K> * p)2871 double get_norme_l2(Matrice_Creuse<K> * p){
2872 if(p==0) return 0.;
2873 HashMatrix<int,K>* ph=p->pHM();
2874 ffassert(ph);
2875 return ph->FrobeniusNorm();}
2876
2877 template<class K>
get_norme_l1(Matrice_Creuse<K> * p)2878 double get_norme_l1(Matrice_Creuse<K> * p){
2879 if(p==0) return 0.;
2880 HashMatrix<int,K>* ph=p->pHM();
2881 ffassert(ph);
2882 return ph->norm1();}
2883
2884 template<class R,int Init>
SetMatrice_Creuse(Matrice_Creuse<R> * p,newpMatrice_Creuse<R> np)2885 Matrice_Creuse<R> * SetMatrice_Creuse(Matrice_Creuse<R> * p,newpMatrice_Creuse<R> np)
2886 {
2887 return np.set(p,Init);
2888 }
2889 // Add F.H July 2019
2890 template<class R>
InitMatrice_Creuse_nm(Matrice_Creuse<R> * const & p,const long & n,const long & m)2891 Matrice_Creuse<R> * InitMatrice_Creuse_nm(Matrice_Creuse<R> * const & p,const long &n,const long &m)
2892 {
2893 p->init() ;
2894 HashMatrix<int,R> *phm= new HashMatrix<int,R>((int) n,(int) m,0,0);
2895 MatriceCreuse<R> *pmc(phm);
2896 p->A.master(pmc);
2897 return p;
2898 }
2899 template<class R>
InitMatrice_Creuse_n(Matrice_Creuse<R> * const & p,const long & n)2900 Matrice_Creuse<R> * InitMatrice_Creuse_n(Matrice_Creuse<R> * const & p,const long &n)
2901 {
2902 p->init() ;
2903 HashMatrix<int,R> *phm= new HashMatrix<int,R>((int)n,(int) n,0,0);
2904 MatriceCreuse<R> *pmc(phm);
2905 p->A.master(pmc);
2906 return p;
2907 }
2908
2909 template<class R,int c>
AddtoMatrice_Creuse(Matrice_Creuse<R> * p,newpMatrice_Creuse<R> np)2910 Matrice_Creuse<R> * AddtoMatrice_Creuse(Matrice_Creuse<R> * p,newpMatrice_Creuse<R> np)
2911 {
2912 return np.add(p,double(c));
2913 }
2914
2915 template<class K, bool init>
set_H_Eye(Matrice_Creuse<K> * pA,const Eye eye)2916 Matrice_Creuse<K>* set_H_Eye(Matrice_Creuse<K> *pA,const Eye eye)
2917 {
2918 int n = eye.n, m=eye.m, nn= min(n,m);
2919 if( init) pA->init();
2920 pA->resize(n,m);
2921 HashMatrix<int,K> * pH= pA->pHM();
2922 ffassert(pH);
2923 pH->clear();
2924 pH->resize(n,m,nn);
2925 for(int i=0; i< n; ++i)
2926 (*pH)(i,i)=1.;
2927 return pA;
2928 }
2929 template <class R>
AddSparseMat()2930 void AddSparseMat()
2931 {
2932 SetMatrix_Op<R>::btype = Dcl_Type<const SetMatrix_Op<R> * >();
2933 Dcl_Type<TheDiagMat<R> >();
2934 Dcl_Type<TheCoefMat<R> >(); // Add FH oct 2005
2935 Dcl_Type< map< pair<int,int>, R> * >(); // Add FH mars 2005
2936 Dcl_Type< minusMat<R> >(); // Add FJH mars 2007
2937
2938 basicForEachType * t_MC=atype< Matrice_Creuse<R>* >();
2939
2940 t_MC->SetTypeLoop(atype< R* >(),atype< long* >(),atype< long* >());
2941
2942 basicForEachType * t_MM=atype<map< pair<int,int>, R> * >();
2943
2944 TheOperators->Add("*",
2945 new OneBinaryOperator<Op2_mulvirtAv<typename RNM_VirtualMatrix<R>::plusAx,Matrice_Creuse<R>*,KN_<R> > >,
2946 new OneBinaryOperator<Op2_mulvirtAv<typename RNM_VirtualMatrix<R>::plusAtx,Matrice_Creuse_Transpose<R>,KN_<R> > >,
2947 new OneBinaryOperator<Op2_mulvirtAv<typename RNM_VirtualMatrix<R>::solveAxeqb,Matrice_Creuse_inv<R>,KN_<R> > > ,
2948 new OneBinaryOperator<Op2_mulvirtAv<typename RNM_VirtualMatrix<R>::solveAtxeqb,Matrice_Creuse_inv_trans<R>,KN_<R> > >
2949 );
2950
2951 TheOperators->Add("^", new OneBinaryOperatorA_inv<R>());
2952 TheOperators->Add("^", new OneBinaryOperatorAt_inv<R>());
2953
2954 // matrix new code FH (Houston 2004)
2955 TheOperators->Add("=",
2956 new OneOperator2<Matrice_Creuse<R>*,Matrice_Creuse<R>*,newpMatrice_Creuse<R> > (SetMatrice_Creuse<R,0> ),
2957 new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,const Matrix_Prod<R,R>,E_F_StackF0F0>(ProdMat<R,R,R,1>),
2958 new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,KN<R> *,E_F_StackF0F0>(DiagMat<R,1>),
2959 new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,Matrice_Creuse_Transpose<R>,E_F_StackF0F0>(CopyTrans<R,R,1>),
2960 new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,Matrice_Creuse<R>*,E_F_StackF0F0>(CopyMat<R,R,1>) ,
2961 new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,KNM<R>*,E_F_StackF0F0>(MatFull2Sparse<R,1>) ,
2962 new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,list<tuple<R,MatriceCreuse<R> *,bool> > *,E_F_StackF0F0>(CombMat<R,1>) ,
2963 new OneOperatorCode<BlockMatrix1<R> >(),
2964 new OneOperator2<Matrice_Creuse<R>*,Matrice_Creuse<R>*,Eye>(set_H_Eye<R,false> )
2965
2966 );
2967 TheOperators->Add("+=",
2968 new OneOperator2<Matrice_Creuse<R>*,Matrice_Creuse<R>*,newpMatrice_Creuse<R> > (AddtoMatrice_Creuse<R, 1> ),
2969 new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,list<tuple<R,MatriceCreuse<R> *,bool> > *,E_F_StackF0F0>(AddCombMat<R,1>));
2970
2971 TheOperators->Add("-=",
2972 new OneOperator2<Matrice_Creuse<R>*,Matrice_Creuse<R>*,newpMatrice_Creuse<R> > (AddtoMatrice_Creuse<R, -1> ),
2973 new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,list<tuple<R,MatriceCreuse<R> *,bool> > *,E_F_StackF0F0>(AddCombMat<R,-1>));
2974
2975
2976
2977 TheOperators->Add("<-",
2978 new OneOperatorCode<BlockMatrix<R> >(),
2979 new OneOperator2<Matrice_Creuse<R>*,Matrice_Creuse<R>*,newpMatrice_Creuse<R> > (SetMatrice_Creuse<R,1> ),
2980 new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,long > (InitMatrice_Creuse_n<R> ),
2981 new OneOperator3_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,long,long > (InitMatrice_Creuse_nm<R> ),
2982 new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,const Matrix_Prod<R,R>,E_F_StackF0F0>(ProdMat<R,R,R,0>),
2983 new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,KN<R> *,E_F_StackF0F0>(DiagMat<R,0>) ,
2984 new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,Matrice_Creuse_Transpose<R>,E_F_StackF0F0>(CopyTrans<R,R,0>),
2985 new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,Matrice_Creuse<R>*,E_F_StackF0F0>(CopyMat<R,R,0>) ,
2986 new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,KNM<R>*,E_F_StackF0F0>(MatFull2Sparse<R,0>) ,
2987 new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,list<tuple<R,MatriceCreuse<R> *,bool> > *,E_F_StackF0F0>(CombMat<R,0>),
2988 new OneOperator2<Matrice_Creuse<R>*,Matrice_Creuse<R>*,Eye>(set_H_Eye<R,true> )
2989
2990
2991 );
2992 TheOperators->Add("*",
2993 new OneBinaryOperator<Op2_pair<Matrix_Prod<R,R>,Matrice_Creuse<R>*,Matrice_Creuse<R>*> >,
2994 new OneBinaryOperator<Op2_pair<Matrix_Prod<R,R>,Matrice_Creuse_Transpose<R>,Matrice_Creuse<R>* > >,
2995 new OneBinaryOperator<Op2_pair<Matrix_Prod<R,R>,Matrice_Creuse_Transpose<R>,Matrice_Creuse_Transpose<R> > >,
2996 new OneBinaryOperator<Op2_pair<Matrix_Prod<R,R>,Matrice_Creuse<R>*,Matrice_Creuse_Transpose<R> > > ,
2997 new OneBinaryOperator<Op2_ListCM<R> > ,
2998 new OneBinaryOperator<Op2_ListMC<R> > ,
2999 new OneBinaryOperator<Op2_ListCMt<R> > ,
3000 new OneBinaryOperator<Op2_ListMtC<R> >
3001
3002 );
3003 TheOperators->Add("+",
3004 new OneBinaryOperator<Op2_ListCMCMadd<R> >,
3005 new OneBinaryOperator<Op2_ListCMMadd<R> >,
3006 new OneBinaryOperator<Op2_ListMCMadd<R> >,
3007 new OneBinaryOperator<Op2_ListMMadd<R> >
3008
3009 );
3010 TheOperators->Add("-",
3011 new OneBinaryOperator<Op2_ListCMCMsub<R> >, // (L) - (L)
3012 new OneBinaryOperator<Op2_ListCMMadd<R,-1> >, // L - M
3013 new OneBinaryOperator<Op2_ListMCMadd<R,-1> >, // M - L
3014 new OneBinaryOperator<Op2_ListMMadd<R,-1> > // M - M
3015
3016 );
3017 TheOperators->Add("-",
3018 new OneUnaryOperator<Op1_LCMd<R,-1> >
3019 );
3020 TheOperators->Add("+",
3021 new OneUnaryOperator<Op1_LCMd<R,1> >
3022 );
3023
3024 Add<Matrice_Creuse<R> *>("COO",".",new OneOperator1<bool,Matrice_Creuse<R> *>(set_mat_COO<R>) );
3025 Add<Matrice_Creuse<R> *>("CSR",".",new OneOperator1<bool,Matrice_Creuse<R> *>(set_mat_CSR<R>) );
3026 Add<Matrice_Creuse<R> *>("CSC",".",new OneOperator1<bool,Matrice_Creuse<R> *>(set_mat_CSC<R>) );
3027 Add<Matrice_Creuse<R> *>("n",".",new OneOperator1<long,Matrice_Creuse<R> *>(get_mat_n<R>) );
3028 Add<Matrice_Creuse<R> *>("m",".",new OneOperator1<long,Matrice_Creuse<R> *>(get_mat_m<R>) );
3029 Add<Matrice_Creuse<R> *>("nbcoef",".",new OneOperator1<long,Matrice_Creuse<R> *>(get_mat_nbcoef<R>) );
3030 Add<Matrice_Creuse<R> *>("nnz",".",new OneOperator1<long,Matrice_Creuse<R> *>(get_mat_nbcoef<R>) );
3031 Add<Matrice_Creuse<R> *>("half",".",new OneOperator1<long,Matrice_Creuse<R> *>(get_mat_half<R>) );
3032 Add<Matrice_Creuse<R> *>("size",".",new OneOperator1<long,Matrice_Creuse<R> *>(get_mat_nbcoef<R>) );
3033 Add<Matrice_Creuse<R> *>("trace",".",new OneOperator1<R,Matrice_Creuse<R>* >(get_trace_mat<R>) );
3034 Add<Matrice_Creuse<R> *>("clear",".",new OneOperator1<bool,Matrice_Creuse<R>* >(clear_mat<R>) );
3035
3036 Add<Matrice_Creuse<R> *>("diag",".",new OneOperator1<TheDiagMat<R> ,Matrice_Creuse<R> *>(thediag<R>) );
3037 Add<Matrice_Creuse<R> *>("coef",".",new OneOperator1<TheCoefMat<R> ,Matrice_Creuse<R> *>(thecoef<R>) );
3038
3039 TheOperators->Add("=", new OneOperator2<KN<R>*,KN<R>*,TheDiagMat<R> >(get_mat_daig<R>) );
3040
3041 TheOperators->Add("<-", new OneOperator2<KN<R>*,KN<R>*,TheDiagMat<R> >(init_get_mat_daig<R>) );
3042 TheOperators->Add("=", new OneOperator2<TheDiagMat<R>,TheDiagMat<R>,KN<R>*>(set_mat_daig<R>) );
3043
3044 // ADD oct 2005
3045 TheOperators->Add("=", new OneOperator2<KN<R>*,KN<R>*,TheCoefMat<R> >(get_mat_coef<R>) );
3046 TheOperators->Add("=", new OneOperator2<TheCoefMat<R>,TheCoefMat<R>,KN<R>*>(set_mat_coef<R>) );
3047 TheOperators->Add("=", new OneOperator2<TheCoefMat<R>,TheCoefMat<R>,R>(set_mat_coef<R>) );
3048
3049 Global.Add("set","(",new SetMatrix<R>);
3050 Add<Matrice_Creuse<R> *>("linfty",".",new OneOperator1<double,Matrice_Creuse<R> *>(get_norme_linfty));
3051 Add<Matrice_Creuse<R> *>("l2",".",new OneOperator1<double,Matrice_Creuse<R> *>(get_norme_l2));
3052 Add<Matrice_Creuse<R> *>("l1",".",new OneOperator1<double,Matrice_Creuse<R> *>(get_norme_l1));
3053 atype<Matrice_Creuse<R> * >()->Add("(","",new OneOperator3_<R*,Matrice_Creuse<R> *,long,long >(1,get_elementp2mc<R>));
3054
3055 atype<Matrice_Creuse<R> * >()->Add("[","",new OneOperator3_<R,Matrice_Creuse<R> *,long,long >(10,get_element2mc<R>));
3056
3057 // typedef map< pair<int,int>, R> MAPMAT;
3058 typedef Matrice_Creuse<R> * MAPMATC;
3059 typedef newpMatrice_Creuse<R> MAPMATN;
3060 atype<KNM<R>*>()->Add("(","",new OneOperator3s_<MAPMATN ,KNM<R>*,Inv_KN_long,Inv_KN_long >(Matrixfull2mapIJ_inv<R>));
3061 atype<KNM<R>*>()->Add("(","",new OneOperator3s_<MAPMATN ,KNM<R>*,KN_<long>,KN_<long> >(Matrixfull2mapIJ<R>));
3062
3063 atype<outProduct_KN_<R>*>()->Add("(","",new OneOperator3s_<MAPMATN ,outProduct_KN_<R>*,Inv_KN_long,Inv_KN_long >(Matrixoutp2mapIJ_inv<R>));
3064 atype<outProduct_KN_<R>*>()->Add("(","",new OneOperator3s_<MAPMATN ,outProduct_KN_<R>*,KN_<long>,KN_<long> >(Matrixoutp2mapIJ<R>));
3065
3066
3067 TheOperators->Add("=", new SetRawMatformMat<R>);
3068
3069
3070
3071 t_MM->Add("(","",new OneOperator3s_<MAPMATN ,MAPMATC ,Inv_KN_long,Inv_KN_long >(Matrixmapp2mapIJ1<R>));
3072 t_MM->Add("(","",new OneOperator3s_<MAPMATN ,MAPMATC ,KN_<long>,KN_<long> >(Matrixmapp2mapIJ<R>));
3073
3074 t_MC->Add("(","",new OneOperator3s_<MAPMATN ,MAPMATC ,Inv_KN_long,Inv_KN_long >(Matrixmapp2mapIJ1<R>,t_MC));
3075 t_MC->Add("(","",new OneOperator3s_<MAPMATN ,MAPMATC ,KN_<long>,KN_<long> >(Matrixmapp2mapIJ<R>,t_MC));
3076
3077 map_type[typeid(MAPMATN ).name()]->AddCast(
3078 new E_F1_funcT<MAPMATN ,KNM<R>* >(Matrixfull2map<R>),
3079 new E_F1_funcT<MAPMATN ,outProduct_KN_<R>* >(Matrixoutp2map<R>)
3080
3081 );
3082
3083 map_type[typeid(list<tuple<R,MatriceCreuse<R> *,bool> > *).name()]->AddCast(
3084 new E_F1_funcT<list<tuple<R,MatriceCreuse<R> *,bool> > *,Matrice_Creuse<R>* >(M2L3<R>),
3085 new E_F1_funcT<list<tuple<R,MatriceCreuse<R> *,bool> > *,Matrice_Creuse_Transpose<R> >(tM2L3<R>),
3086 new E_F1_funcT<list<tuple<R,MatriceCreuse<R> *,bool> > *,minusMat<R> >(mM2L3<R> )
3087 );
3088
3089
3090 TheOperators->Add("{}",new ForAllLoop<E_ForAllLoopMatrix<R> >);
3091 // remove 2 plugin thresholding and symmetrizeCSR
3092 typedef Thresholding<R> TMR;
3093 typedef Matrice_Creuse<R> MR;
3094 Dcl_Type<TMR>();
3095 TMR t(0);
3096 //thresholding2(t, 0.);
3097 Add<MR *>("thresholding", ".", new OneOperator1<TMR, MR *>(to_Thresholding));
3098 Add<TMR>("(", "", new OneOperator2_<MR *, TMR, double>(thresholding2));
3099 Global.Add("symmetrizeCSR", "(", new OneOperator1_<long, Matrice_Creuse<R> *>(symmetrizeCSR<R> ));
3100 // --- end
3101 }
3102
3103 extern int lineno();
3104 class PrintErrorCompile : public OneOperator {
3105 public:
3106 const char * cmm;
code(const basicAC_F0 &) const3107 E_F0 * code(const basicAC_F0 & ) const
3108 { ErrorCompile(cmm,lineno());
3109 return 0;}
PrintErrorCompile(const char * cc)3110 PrintErrorCompile(const char * cc): OneOperator(map_type[typeid(R).name()]),cmm(cc){}
3111
3112 };
3113
3114 class PrintErrorCompileIM : public E_F0info { public:
3115 typedef double Result;
f(const basicAC_F0 & args)3116 static E_F0 * f(const basicAC_F0 & args)
3117 {
3118 lgerror("\n\n *** change interplotematrix in interpole.\n *** Bad name in previous version,\n *** sorry FH.\n\n");
3119 return 0; }
typeargs()3120 static ArrayOfaType typeargs() {return ArrayOfaType(true);}
operator aType() const3121 operator aType () const { return atype<double>();}
3122
3123 };
3124
3125 template<class T>
3126 class removeDOF_Op : public E_F0mps {
3127 public:
3128 Expression eA;
3129 Expression eR;
3130 Expression ex;
3131 Expression eout;
3132 bool transpose;
3133 static const int n_name_param = 4;
3134 static basicAC_F0::name_and_type name_param[];
3135 Expression nargs[n_name_param];
removeDOF_Op(const basicAC_F0 & args,Expression param1,Expression param2,Expression param3,Expression param4)3136 removeDOF_Op(const basicAC_F0& args, Expression param1, Expression param2, Expression param3, Expression param4)
3137 : eA(param1), eR(param2), ex(param3), eout(param4), transpose(false) {
3138 args.SetNameParam(n_name_param, name_param, nargs);
3139 }
removeDOF_Op(const basicAC_F0 & args,Expression param2,Expression param3,Expression param4,bool t,int)3140 removeDOF_Op(const basicAC_F0& args, Expression param2, Expression param3, Expression param4, bool t, int)
3141 : eA(0), eR(param2), ex(param3), eout(param4), transpose(t) {
3142 args.SetNameParam(n_name_param, name_param, nargs);
3143 }
removeDOF_Op(const basicAC_F0 & args,Expression param1,Expression param2)3144 removeDOF_Op(const basicAC_F0& args, Expression param1, Expression param2)
3145 : eA(param1), eR(param2), ex(0), eout(0), transpose(false) {
3146 args.SetNameParam(n_name_param, name_param, nargs);
3147 }
3148
3149 AnyType operator()(Stack stack) const;
3150 };
3151
3152 template<class T>
3153 basicAC_F0::name_and_type removeDOF_Op<T>::name_param[] = {
3154 {"symmetrize", &typeid(bool)},
3155 {"condensation", &typeid(KN<long>*)},
3156 {"R", &typeid(Matrice_Creuse<double>*)},
3157 {"eps",&typeid(double)}
3158 };
3159
3160 template<class T>
3161 class removeDOF : public OneOperator {
3162 const unsigned short withA;
3163 public:
removeDOF()3164 removeDOF() : OneOperator(atype<long>(), atype<Matrice_Creuse<T>*>(), atype<Matrice_Creuse<double>*>(), atype<KN<T>*>(), atype<KN<T>*>()),withA(1) {}
removeDOF(int)3165 removeDOF(int) : OneOperator(atype<long>(), atype<Matrice_Creuse<double>*>(), atype<KN<T>*>(), atype<KN<T>*>()),withA(0) {}
removeDOF(int,int)3166 removeDOF(int, int) : OneOperator(atype<long>(), atype<Matrice_Creuse<T>*>(), atype<Matrice_Creuse<double>*>()),withA(2) {}
removeDOF(int,int,int)3167 removeDOF(int, int, int) : OneOperator(atype<long>(), atype<Matrice_Creuse_Transpose<double>>(), atype<KN<T>*>(), atype<KN<T>*>()),withA(3) {}
3168
code(const basicAC_F0 & args) const3169 E_F0* code(const basicAC_F0& args) const {
3170 if(withA == 1)
3171 return new removeDOF_Op<T>(args, t[0]->CastTo(args[0]), t[1]->CastTo(args[1]), t[2]->CastTo(args[2]), t[3]->CastTo(args[3]));
3172 else if(withA == 0 || withA == 3)
3173 return new removeDOF_Op<T>(args, t[0]->CastTo(args[0]), t[1]->CastTo(args[1]), t[2]->CastTo(args[2]), withA == 3, 1);
3174 else
3175 return new removeDOF_Op<T>(args, t[0]->CastTo(args[0]), t[1]->CastTo(args[1]));
3176 }
3177 };
cmp(const std::pair<unsigned int,T> & lhs,const std::pair<unsigned int,T> & rhs)3178 template<class T> bool cmp(const std::pair<unsigned int, T>& lhs, const std::pair<unsigned int, T>& rhs) { return lhs.first < rhs.first; }
3179 template<class T>
operator ()(Stack stack) const3180 AnyType removeDOF_Op<T>::operator()(Stack stack) const {
3181
3182 static const double defEPS=1e-12;
3183 typedef double R;
3184 // code wri-ng no ...
3185 Matrice_Creuse<T>* pA = eA ? GetAny<Matrice_Creuse<T>* >((*eA)(stack)):0;
3186 Matrice_Creuse<R>* pR;
3187 if(transpose) {
3188 Matrice_Creuse_Transpose<R> tpR = GetAny<Matrice_Creuse_Transpose<R>>((*eR)(stack));
3189 pR = tpR;
3190 }
3191 else
3192 pR = GetAny<Matrice_Creuse<R>*>((*eR)(stack));
3193 KN<T>* pX = ex ? GetAny<KN<T>* >((*ex)(stack)) : 0;
3194 KN<T>* pOut = eout ? GetAny<KN<T>* >((*eout)(stack)) : 0;
3195 Matrice_Creuse<R>* pC = nargs[2] ? GetAny<Matrice_Creuse<R>*>((*nargs[2])(stack)) : 0;
3196 MatriceMorse<R> *mC = 0;
3197 if(pC && pC->A) {
3198 mC = static_cast<MatriceMorse<R>*>(&(*pC->A));
3199 }
3200 ffassert(pR);
3201 bool rhs = (pX && pOut) && (pOut->n > 0 || pX->n > 0);
3202 if(pA)
3203 {
3204 if(!pC)
3205 pC = pR;
3206 MatriceMorse<T> *mA = pA->pHM();
3207 MatriceMorse<R> *mR = pR->pHM();
3208 if(!mC)
3209 mC = mR;
3210 pA->Uh = pR->Uh;
3211 pA->Vh = pC->Vh;
3212
3213 bool symmetrize = nargs[0] ? GetAny<bool>((*nargs[0])(stack)) : false;
3214 double EPS=nargs[3] ? GetAny<double>((*nargs[3])(stack)) :defEPS ;
3215 KN<long>* condensation = nargs[1] ? GetAny<KN<long>* >((*nargs[1])(stack)) : (KN<long>*) 0;
3216 ffassert(condensation ||( mR && mC) );
3217 unsigned int n = condensation ? condensation->n : mR->nnz;
3218 unsigned int m = condensation ? condensation->n : mC->nnz;
3219 KN<int> lg(n+1,0);
3220
3221 if(rhs && pOut->n != n) pOut->resize(n);
3222 mC->COO();
3223 mR->COO();
3224 mA->CSR();
3225 std::vector<signed int> tmpVec;
3226 if(!condensation)
3227 {
3228 tmpVec.resize(mA->m);
3229 for(long i = 0; i < m; ++i)
3230 tmpVec[mC->j[i]] = i + 1;
3231 if(!mA->half) {
3232 std::vector<std::pair<int, T> > tmp;
3233 tmp.reserve(mA->nnz);
3234
3235 lg[0] = 0;
3236 for(long i = 0; i < n; ++i) {
3237 for(long j = mA->p[mR->j[i]]; j < mA->p[mR->j[i] + 1]; ++j) {
3238 long col = tmpVec[mA->j[j]];
3239 if(col != 0 && abs(mA->aij[j]) > EPS) {
3240 if(symmetrize) {
3241 if(col - 1 <= i)
3242 tmp.push_back(std::make_pair(col - 1, mA->aij[j]));
3243 }
3244 else
3245 tmp.push_back(std::make_pair(col - 1, mA->aij[j]));
3246
3247 }
3248 }
3249 std::sort(tmp.begin() + lg[i], tmp.end(),cmp<T> );
3250 // c++11 , [](const std::pair<unsigned int, T>& lhs, const std::pair<unsigned int, T>& rhs) { return lhs.first < rhs.first; });
3251 if(rhs)
3252 *(*pOut + i) = *(*pX + mC->j[i]);
3253 lg[i + 1] = tmp.size();
3254 }
3255 mA->clear();
3256 mA->resize(n,m);
3257 MatriceMorse<T> &A = *mA;
3258 A.half = symmetrize;
3259 for(int i=0; i<n; ++i)
3260 {
3261 for(int k= lg[i]; k < lg[i+1]; ++k)
3262 {
3263 int j= tmp[k].first;
3264 T aij = tmp[k].second;
3265 A(i,j) =aij;
3266 }
3267 }
3268
3269
3270 }// !mA->Half
3271 else {
3272 std::vector<std::vector<std::pair<unsigned int, T> > > tmp(n);
3273 for(unsigned int i = 0; i < n; ++i)
3274 tmp[i].reserve(mA->p[mR->j[i] + 1] - mA->p[mR->j[i]]);
3275
3276 unsigned int nnz = 0;
3277 for(unsigned int i = 0; i < n; ++i) {
3278 for(unsigned int j = mA->p[mR->j[i]]; j < mA->p[mR->j[i] + 1]; ++j) {
3279 unsigned int col = tmpVec[mA->j[j]];
3280 if(col != 0 && abs(mA->aij[j]) > EPS) {
3281 if(i < col - 1)
3282 tmp[col - 1].push_back(make_pair(i, mA->aij[j]));
3283 else
3284 tmp[i].push_back(make_pair(col - 1, mA->aij[j]));
3285 ++nnz;
3286 }
3287 }
3288 if(rhs)
3289 *(*pOut + i) = *(*pX + mC->j[i]);
3290 }
3291 int Half = mA->half;
3292 mA->clear();
3293 mA->resize(n,m,nnz);
3294 MatriceMorse<T> &A = *mA;
3295 A.half = Half;
3296 for(unsigned int i = 0; i < n; ++i) {
3297 std::sort(tmp[i].begin(), tmp[i].end(),cmp<T>);
3298 // c++11, [](const std::pair<unsigned int, T>& lhs, const std::pair<unsigned int, T>& rhs) { return lhs.first < rhs.first; });
3299 for(typename std::vector<std::pair<unsigned int, T> >::const_iterator it = tmp[i].begin(); it != tmp[i].end(); ++it)
3300 A(i,it->first) =it->second;
3301
3302 }
3303
3304 }
3305 }
3306
3307 else
3308 {
3309 tmpVec.reserve(mA->n);
3310 unsigned int i = 0, j = 1;
3311 for(unsigned int k = 0; k < mA->n; ++k) {
3312 if(k == *(*condensation + i)) {
3313 ++i;
3314 tmpVec.push_back(i);
3315 }
3316 else {
3317 tmpVec.push_back(-j);
3318 ++j;
3319 }
3320
3321 }
3322
3323 std::vector<std::pair<int, T> > tmpInterior;
3324 std::vector<std::pair<int, T> > tmpBoundary;
3325 std::vector<std::pair<int, T> > tmpInteraction;
3326 tmpInterior.reserve(mA->nnz);
3327 tmpBoundary.reserve(mA->nnz);
3328 tmpInteraction.reserve(mA->nnz);
3329
3330 lg[0] = 0;
3331 for(long i = 0; i < mA->n; ++i) {
3332 int row = tmpVec[i];
3333 if(row < 0) {
3334 for(unsigned int j = mA->p[i]; j < mA->p[i + 1]; ++j) {
3335 int col = tmpVec[mA->j[j]];
3336 if(col < 0)
3337 tmpInterior.push_back(make_pair(-col - 1, mA->aij[j]));
3338 else
3339 tmpInteraction.push_back(make_pair(col - 1, mA->aij[j]));
3340 }
3341
3342 }
3343 else {
3344 for(unsigned int j = mA->p[i]; j < mA->p[i + 1]; ++j) {
3345 int col = tmpVec[mA->j[j]];
3346 if(col > 0)
3347 tmpBoundary.push_back(make_pair(col - 1, mA->aij[j]));
3348 }
3349 if(rhs)
3350 *(*pOut + i) = *(*pX + *(*condensation + i));
3351 lg[i + 1] = tmpBoundary.size();
3352 }
3353 }
3354
3355 mA->clear();
3356 mA->resize(n,n);
3357 MatriceMorse<T> &mR = *new MatriceMorse<T>(n,m,tmpBoundary.size(),0);
3358 for(unsigned int i = 0; i < tmpBoundary.size(); ++i) {
3359 mR(i, tmpBoundary[i].first)= tmpBoundary[i].second;
3360 }
3361 ffassert(0);
3362 pR->typemat = 0; //TypeSolveMat(TypeSolveMat::GMRES);
3363 // bug ici ::: FH.. pR->A.master(&mR);
3364 }
3365
3366 }
3367 else if(rhs)
3368 {
3369
3370
3371 MatriceMorse<R> *mR = pR->pHM();
3372 if(mR && mR->n && mR->m) {
3373 mR->COO();
3374 unsigned int n = mR->nnz;
3375 if(transpose) {
3376 if(pOut->n != mR->m) pOut->resize(mR->m);
3377 *pOut = T();
3378 for(unsigned int i = 0; i < n; ++i) {
3379 *(*pOut + mR->j[i]) = *(*pX + mR->i[i]);
3380 }
3381 }
3382 else {
3383 if(pOut->n != n) pOut->resize(n);
3384 for(unsigned int i = 0; i < n; ++i) {
3385 *(*pOut + i) = *(*pX + mR->j[i]);
3386 }
3387 }
3388 }
3389 else
3390 *pOut = T();
3391 }
3392 return 0L;
3393 }
3394
3395
SparseDefault()3396 bool SparseDefault()
3397 {
3398 return 1;//TypeSolveMat::SparseSolver== TypeSolveMat::defaultvalue;
3399 }
3400
3401 bool Have_UMFPACK_=false;
Have_UMFPACK()3402 bool Have_UMFPACK() { return Have_UMFPACK_;}
3403
removeHalf(MatriceMorse<R> & A,long half,double tol)3404 MatriceMorse<R> * removeHalf(MatriceMorse<R> & A,long half,double tol)
3405 {
3406
3407 // half < 0 => L
3408 // half > 0 => U
3409 // half = 0 => L and the result will be sym
3410 int nnz =0;
3411
3412 if( A.half )
3413 return &A;// copy
3414 // do alloc
3415 MatriceMorse<R> *r=new MatriceMorse<R>(A);
3416 r->RemoveHalf(half,tol);
3417 if(verbosity )
3418 cout << " removeHalf: new nnz = "<< r->nnz << " "<< r->half << endl;
3419
3420 return r;
3421 }
3422
removeHalf(Stack stack,Matrice_Creuse<R> * const & pA,long const & half,const double & tol)3423 newpMatrice_Creuse<R> removeHalf(Stack stack,Matrice_Creuse<R> *const & pA,long const & half,const double & tol)
3424 {
3425 MatriceCreuse<R> * pa=pA->A;
3426 MatriceMorse<R> *pma= dynamic_cast<MatriceMorse<R>* > (pa);
3427 ffassert(pma);
3428 return newpMatrice_Creuse<R>(stack,removeHalf(*pma,half,tol));
3429 }
3430
removeHalf(Stack stack,Matrice_Creuse<R> * const & pR,Matrice_Creuse<R> * const & pA,long const & half,const double & tol)3431 bool removeHalf(Stack stack,Matrice_Creuse<R> *const & pR,Matrice_Creuse<R> *const & pA,long const & half,const double & tol)
3432 {
3433 MatriceCreuse<R> * pa=pA->A;
3434 MatriceMorse<R> *pma= dynamic_cast<MatriceMorse<R>* > (pa);
3435 MatriceCreuse<R> * pr= removeHalf(*pma,half,tol);
3436
3437 pR->A.master(pr);
3438 return true;
3439 }
removeHalf(Stack stack,Matrice_Creuse<R> * const & pA,long const & half)3440 newpMatrice_Creuse<R> removeHalf(Stack stack,Matrice_Creuse<R> *const & pA,long const & half)
3441 {
3442 return removeHalf(stack,pA,half,-1.);
3443 }
3444
3445 template<class K>
3446 class plotMatrix : public OneOperator {
3447 public:
3448
3449 class Op : public E_F0info {
3450 public:
3451 Expression a;
3452
3453 static const int n_name_param = 1;
3454 static basicAC_F0::name_and_type name_param[] ;
3455 Expression nargs[n_name_param];
arg(int i,Stack stack,bool a) const3456 bool arg(int i,Stack stack,bool a) const{ return nargs[i] ? GetAny<bool>( (*nargs[i])(stack) ): a;}
arg(int i,Stack stack,long a) const3457 long arg(int i,Stack stack,long a) const{ return nargs[i] ? GetAny<long>( (*nargs[i])(stack) ): a;}
3458
3459 public:
Op(const basicAC_F0 & args,Expression aa)3460 Op(const basicAC_F0 & args,Expression aa) : a(aa) {
3461 args.SetNameParam(n_name_param,name_param,nargs);
3462 }
3463
operator ()(Stack stack) const3464 AnyType operator()(Stack stack) const{
3465
3466 if (mpirank == 0 && ThePlotStream) {
3467 Matrice_Creuse<K>* A =GetAny<Matrice_Creuse<K>* >((*a)(stack));
3468 bool wait = arg(0,stack,false);
3469
3470 PlotStream theplot(ThePlotStream);
3471 theplot.SendNewPlot();
3472 theplot << 3L;
3473 theplot <= wait;
3474 theplot.SendEndArgPlot();
3475 theplot.SendMeshes();
3476 theplot << 0L;
3477 theplot.SendPlots();
3478 theplot << 1L;
3479 theplot << 31L;
3480
3481 HashMatrix<int,K>* ph=A->pHM();
3482
3483 if (!ph) {
3484 theplot << 0;
3485 theplot << 0;
3486 theplot << 0L;
3487 theplot << 0L;
3488 }
3489 else {
3490 theplot << (int)ph->n;
3491 theplot << (int)ph->m;
3492 theplot << 0L;
3493 theplot << (long)ph->nnz;
3494
3495 for (int i=0;i<ph->nnz;i++) {
3496 theplot << ph->i[i];
3497 theplot << ph->j[i];
3498 theplot << 1;
3499 theplot << 1;
3500 theplot << 1;
3501 theplot << abs(ph->aij[i]);
3502 }
3503 }
3504
3505 theplot.SendEndPlot();
3506
3507 }
3508
3509 return 0L;
3510 }
3511 };
3512
plotMatrix()3513 plotMatrix() : OneOperator(atype<long>(),atype<Matrice_Creuse<K>*>()) {}
3514
code(const basicAC_F0 & args) const3515 E_F0 * code(const basicAC_F0 & args) const
3516 {
3517 return new Op(args,t[0]->CastTo(args[0]));
3518 }
3519 };
3520
3521 template<class K>
3522 basicAC_F0::name_and_type plotMatrix<K>::Op::name_param[]= {
3523 { "wait", &typeid(bool)}
3524 };
3525
init_lgmat()3526 void init_lgmat()
3527
3528 {
3529
3530 Dcl_Type<const MatrixInterpolation<pfes,pfes>::Op *>();
3531 Dcl_Type<const MatrixInterpolation<pfes3,pfes3>::Op *>();
3532 Dcl_Type<const MatrixInterpolation<pfesS,pfesS>::Op *>();
3533 Dcl_Type<const MatrixInterpolation<pfesL,pfesL>::Op *>();
3534 Dcl_Type<const MatrixInterpolation<pfesS,pfes3>::Op *>();
3535 Dcl_Type<const MatrixInterpolation<pfesL,pfes>::Op *>();
3536 Dcl_Type<const MatrixInterpolation<pfesL,pfesS>::Op *>();
3537
3538 map_type_of_map[make_pair(atype<Matrice_Creuse<double>* >(),atype<double*>())]=atype<Matrice_Creuse<double> *>();
3539 map_type_of_map[make_pair(atype<Matrice_Creuse<double>* >(),atype<Complex*>())]=atype<Matrice_Creuse<Complex> *>();
3540 AddSparseMat<double>();
3541 AddSparseMat<Complex>();
3542
3543 Add<const MatrixInterpolation<pfes,pfes>::Op *>("<-","(", new MatrixInterpolation<pfes,pfes>);
3544 Add<const MatrixInterpolation<pfes,pfes>::Op *>("<-","(", new MatrixInterpolation<pfes,pfes>(1));
3545 Add<const MatrixInterpolation<pfes,pfes>::Op *>("<-","(", new MatrixInterpolation<pfes3,pfes3>);
3546 Add<const MatrixInterpolation<pfes,pfes>::Op *>("<-","(", new MatrixInterpolation<pfes3,pfes3>(1,1));
3547 Add<const MatrixInterpolation<pfes,pfes>::Op *>("<-","(", new MatrixInterpolation<pfesS,pfesS>);
3548 Add<const MatrixInterpolation<pfes,pfes>::Op *>("<-","(", new MatrixInterpolation<pfesS,pfesS>(1,1));
3549 Add<const MatrixInterpolation<pfes,pfes>::Op *>("<-","(", new MatrixInterpolation<pfesL,pfesL>);
3550 Add<const MatrixInterpolation<pfes,pfes>::Op *>("<-","(", new MatrixInterpolation<pfesL,pfesL>(1,1));
3551 Add<const MatrixInterpolation<pfes,pfes>::Op *>("<-","(", new MatrixInterpolation<pfesS,pfes3>);
3552 //Add<const MatrixInterpolation<pfes,pfes>::Op *>("<-","(", new MatrixInterpolation<pfesS,pfes3>(1,1));
3553 Add<const MatrixInterpolation<pfes,pfes>::Op *>("<-","(", new MatrixInterpolation<pfesL,pfes>);
3554 // Add<const MatrixInterpolation<pfes,pfes>::Op *>("<-","(", new MatrixInterpolation<pfesL,pfes>(1,1));
3555 Add<const MatrixInterpolation<pfes,pfes>::Op *>("<-","(", new MatrixInterpolation<pfesL,pfesS>);
3556 // Add<const MatrixInterpolation<pfes,pfes>::Op *>("<-","(", new MatrixInterpolation<pfesL,pfesS>(1,1));
3557
3558 Dcl_Type<const RestrictArray<pfes>::Op *>();
3559 Dcl_Type<const RestrictArray<pfes3>::Op *>();
3560 Dcl_Type<const RestrictArray<pfesS>::Op *>();
3561 Dcl_Type<const RestrictArray<pfesL>::Op *>();
3562
3563 Global.Add("restrict","(",new RestrictArray<pfes>);// FH Jan 2014
3564 Global.Add("restrict","(",new RestrictArray<pfes3>);// FH Jan 2014
3565 Global.Add("restrict","(",new RestrictArray<pfesS>);// PHT Apr 2019
3566 Global.Add("restrict","(",new RestrictArray<pfesL>);
3567
3568 TheOperators->Add("=",
3569 new OneOperator2_<KN<long>*,KN<long>*,const RestrictArray<pfes>::Op*,E_F_StackF0F0>(SetRestrict<pfes,1>),
3570 new OneOperator2_<KN<long>*,KN<long>*,const RestrictArray<pfes3>::Op*,E_F_StackF0F0>(SetRestrict<pfes3,1>),
3571 new OneOperator2_<KN<long>*,KN<long>*,const RestrictArray<pfesS>::Op*,E_F_StackF0F0>(SetRestrict<pfesS,1>),
3572 new OneOperator2_<KN<long>*,KN<long>*,const RestrictArray<pfesL>::Op*,E_F_StackF0F0>(SetRestrict<pfesL,1>)
3573 );
3574 TheOperators->Add("<-",
3575 new OneOperator2_<KN<long>*,KN<long>*,const RestrictArray<pfes>::Op*,E_F_StackF0F0>(SetRestrict<pfes,0>),
3576 new OneOperator2_<KN<long>*,KN<long>*,const RestrictArray<pfes3>::Op*,E_F_StackF0F0>(SetRestrict<pfes3,0>),
3577 new OneOperator2_<KN<long>*,KN<long>*,const RestrictArray<pfesS>::Op*,E_F_StackF0F0>(SetRestrict<pfesS,0>),
3578 new OneOperator2_<KN<long>*,KN<long>*,const RestrictArray<pfesL>::Op*,E_F_StackF0F0>(SetRestrict<pfesL,0>)
3579 );
3580
3581
3582 Global.Add("interpolate","(",new MatrixInterpolation<pfes,pfes>);
3583 Global.Add("interpolate","(",new MatrixInterpolation<pfes,pfes>(1));
3584 Global.Add("interpolate","(",new MatrixInterpolation<pfes3,pfes3>);
3585 Global.Add("interpolate","(",new MatrixInterpolation<pfes3,pfes3>(1,1));
3586 Global.Add("interpolate","(",new MatrixInterpolation<pfesS,pfesS>);
3587 Global.Add("interpolate","(",new MatrixInterpolation<pfesS,pfesS>(1,1));
3588 Global.Add("interpolate","(",new MatrixInterpolation<pfesL,pfesL>);
3589 Global.Add("interpolate","(",new MatrixInterpolation<pfesL,pfesL>(1,1));
3590
3591 Global.Add("interpolate","(",new MatrixInterpolation<pfesS,pfes3>);
3592 // Global.Add("interpolate","(",new MatrixInterpolation<pfesS,pfes3>(1,1));
3593 Global.Add("interpolate","(",new MatrixInterpolation<pfesL,pfes>);
3594 // Global.Add("interpolate","(",new MatrixInterpolation<pfesL,pfes>(1,1));
3595 Global.Add("interpolate","(",new MatrixInterpolation<pfesL,pfesS>);
3596 // Global.Add("interpolate","(",new MatrixInterpolation<pfesL,pfesS>(1,1));
3597
3598 Global.Add("interplotematrix","(",new OneOperatorCode<PrintErrorCompileIM>);
3599 zzzfff->Add("mapmatrix",atype<map< pair<int,int>, double> *>());
3600 zzzfff->Add("Cmapmatrix",atype<map< pair<int,int>, Complex> *>()); // a voir
3601
3602 Dcl_Type< Resize<Matrice_Creuse<double> > > ();
3603
3604 Add<Matrice_Creuse<double> *>("resize",".",new OneOperator1< Resize<Matrice_Creuse<double> >,Matrice_Creuse<double> *>(to_Resize));
3605 Add<Resize<Matrice_Creuse<double> > >("(","",new OneOperator3_<Matrice_Creuse<double> *,Resize<Matrice_Creuse<double> > , long, long >(resize2));
3606 // add missing in
3607 Dcl_Type< Resize<Matrice_Creuse<Complex> > > ();
3608 Add<Matrice_Creuse<Complex> *>("resize",".",new OneOperator1< Resize<Matrice_Creuse<Complex> >,Matrice_Creuse<Complex> *>(to_Resize));
3609 Add<Resize<Matrice_Creuse<Complex> > >("(","",new OneOperator3_<Matrice_Creuse<Complex> *,Resize<Matrice_Creuse<Complex> > , long, long >(resize2));
3610
3611
3612 // pour compatibiliter
3613
3614 TheOperators->Add("=",
3615 new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,const MatrixInterpolation<pfes,pfes>::Op*,E_F_StackF0F0>(SetMatrixInterpolation<1>),
3616 new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,const MatrixInterpolation<pfes3,pfes3>::Op*,E_F_StackF0F0>(SetMatrixInterpolation3<1>),
3617 new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,const MatrixInterpolation<pfesS,pfesS>::Op*,E_F_StackF0F0>(SetMatrixInterpolationS<1>),
3618 new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,const MatrixInterpolation<pfesL,pfesL>::Op*,E_F_StackF0F0>(SetMatrixInterpolationL<1>),
3619 new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,const MatrixInterpolation<pfesL,pfes>::Op*,E_F_StackF0F0>(SetMatrixInterpolationL2<1>),
3620 new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,const MatrixInterpolation<pfesL,pfesS>::Op*,E_F_StackF0F0>(SetMatrixInterpolationLS<1>),
3621 new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,const MatrixInterpolation<pfesS,pfes3>::Op*,E_F_StackF0F0>(SetMatrixInterpolationS3<1>)
3622
3623 );
3624
3625
3626 TheOperators->Add("<-",
3627 new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,const MatrixInterpolation<pfes,pfes>::Op*,E_F_StackF0F0>(SetMatrixInterpolation<0>),
3628 new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,const MatrixInterpolation<pfes3,pfes3>::Op*,E_F_StackF0F0>(SetMatrixInterpolation3<0>),
3629 new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,const MatrixInterpolation<pfesS,pfesS>::Op*,E_F_StackF0F0>(SetMatrixInterpolationS<0>),
3630 new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,const MatrixInterpolation<pfesL,pfesL>::Op*,E_F_StackF0F0>(SetMatrixInterpolationL<0>),
3631 new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,const MatrixInterpolation<pfesL,pfesS>::Op*,E_F_StackF0F0>(SetMatrixInterpolationLS<0>),
3632 new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,const MatrixInterpolation<pfesL,pfes>::Op*,E_F_StackF0F0>(SetMatrixInterpolationL2<0>),
3633 new OneOperator2_<Matrice_Creuse<R>*,Matrice_Creuse<R>*,const MatrixInterpolation<pfesS,pfes3>::Op*,E_F_StackF0F0>(SetMatrixInterpolationS3<0>)
3634 );
3635 // construction of complex matrix form a double matrix
3636 TheOperators->Add("=", new OneOperator2_<Matrice_Creuse<Complex>*,Matrice_Creuse<Complex>*,Matrice_Creuse<double>*,E_F_StackF0F0>(CopyMat<R,Complex,1>)
3637 );
3638
3639 TheOperators->Add("<-", new OneOperator2_<Matrice_Creuse<Complex>*,Matrice_Creuse<Complex>*,Matrice_Creuse<double>*,E_F_StackF0F0>(CopyMat<R,Complex,0>)
3640 );
3641 Dcl_Type<Matrice_Creuse_C2R>();
3642 Add<Matrice_Creuse<Complex>*>("re",".",new OneOperator1s_<newpMatrice_Creuse<double> ,Matrice_Creuse<Complex>* >(Build_Matrice_Creuse_C2R<0> ));
3643 Add<Matrice_Creuse<Complex>*>("im",".",new OneOperator1s_<newpMatrice_Creuse<double> ,Matrice_Creuse<Complex>* >(Build_Matrice_Creuse_C2R<1> ));
3644 // construction of complex matrix form a double matrix
3645 init_UMFPack_solver();
3646 init_HashMatrix ();
3647
3648 Global.Add("renumbering", "(", new removeDOF<double>);
3649 Global.Add("renumbering", "(", new removeDOF<Complex>);
3650 Global.Add("renumbering", "(", new removeDOF<double>(1));
3651 Global.Add("renumbering", "(", new removeDOF<Complex>(1));
3652 Global.Add("renumbering", "(", new removeDOF<double>(1, 1));
3653 Global.Add("renumbering", "(", new removeDOF<Complex>(1, 1));
3654 Global.Add("renumbering", "(", new removeDOF<double>(1, 1, 1));
3655 Global.Add("renumbering", "(", new removeDOF<Complex>(1, 1, 1));
3656
3657 Global.Add("display", "(", new plotMatrix<double>);
3658 Global.Add("display", "(", new plotMatrix<Complex>);
3659
3660 // ZZZZZZ ne marche pas FH....
3661 TheOperators->Add("*",
3662 new Op2_mulvirtAvCR< RNM_VirtualMatrix<Complex>::plusAx,Matrice_Creuse<double>*,KN_<Complex> > ,
3663 new Op2_mulvirtAvCR< RNM_VirtualMatrix<Complex>::plusAtx,Matrice_Creuse_Transpose<double>,KN_<Complex> > ,
3664 new Op2_mulvirtAvCR< RNM_VirtualMatrix<Complex>::solveAxeqb,Matrice_Creuse_inv<R>,KN_<Complex> >,
3665 new Op2_mulvirtAvCR< RNM_VirtualMatrix<Complex>::solveAtxeqb,Matrice_Creuse_inv<R>,KN_<Complex> >
3666 );
3667 init_SparseLinearSolver();
3668
3669
3670 Global.New("DefaultSolver",CPValue<string*>(def_solver));
3671 Global.New("DefaultSolverSym",CPValue<string*>(def_solver_sym));
3672 Global.New("DefaultSolverSDP",CPValue<string*>(def_solver_sym_dp));
3673
3674 Global.Add("removeHalf", "(", new OneOperator2s_<newpMatrice_Creuse<R> ,Matrice_Creuse<R> * ,long>(removeHalf));
3675 Global.Add("removeHalf", "(", new OneOperator3s_<newpMatrice_Creuse<R> ,Matrice_Creuse<R> * ,long,double>(removeHalf));
3676 Global.Add("removeHalf", "(", new OneOperator4s_<bool,Matrice_Creuse<R> * ,Matrice_Creuse<R> * ,long,double>(removeHalf));
3677
3678 }
3679
Data_Sparse_Solver_version()3680 int Data_Sparse_Solver_version() { return VDATASPARSESOLVER;}
3681 #else
3682 #include <iostream>
3683
init_lgmat()3684 void init_lgmat()
3685 { std::cout << "\n warning init_lgmat EMPTY\n"<< std::endl;
3686
3687 }
3688 #endif
3689