1 #ifndef __HashMatrix_HPP__
2 #define __HashMatrix_HPP__
3
4 #include <cstddef>
5 #include <algorithm>
6 #include <map>
7 #include <list>
8
9 #include <tuple>
10 using std::tuple;
11 #include <iostream>
12 #include <cassert>
13 #include <complex>
14 #include <cstdint>
15 #include "RNM.hpp"
16 #include "RefCounter.hpp"
17 #include "VirtualMatrix.hpp"
18
19 using std::max;
20 using std::min;
21 using std::cout;
22 using std::endl;
23 using std::pair;
24 using std::swap;
25 using std::complex;
26 using std::list;
27 extern long verbosity ;
28
29 void init_HashMatrix ();
30 template<class Z,class ZZ>
roll64(Z y,ZZ r)31 inline uint64_t roll64(Z y,ZZ r){uint64_t x=y; r %= 64; return (x<<r) | (x << (64-r)) ;}
32
33
34
35 int WhichMatrix(istream & f);
36 template<class TypeIndex,class TypeScalaire> class HashMatrix;
37
38 template<class I,class R,class K>
39 void Addto(HashMatrix<I,R> *P0,const HashMatrix<I,K> *PA,R (*f)(K) ,bool trans=false, I ii00=0,I jj00=0);
40
41
42
43 template<class TypeIndex,class TypeScalaire>
44 class HashMatrix : public VirtualMatrix<TypeIndex,TypeScalaire>
45 {
46 static void HeapSort(TypeIndex *ii,TypeIndex *jj, TypeScalaire *aij,long n);
47
48
49 public:
50
conj(R * x,TypeIndex nnz)51 template<class R> static void conj(R *x,TypeIndex nnz){}
conj(complex<double> * x,TypeIndex nnz)52 static void conj(complex<double> *x,TypeIndex nnz){ for(int k=0; k<nnz; ++k) x[k]=std::conj(x[k]) ; }
conj(complex<float> * x,TypeIndex nnz)53 static void conj(complex<float> *x,TypeIndex nnz){for(int k=0; k<nnz; ++k) x[k]=std::conj(x[k]) ; }
conj(double x)54 static double conj(double x){return x;}
conj(float x)55 static float conj(float x){return x;}
conj(const complex<double> & x)56 static complex<double> conj(const complex<double> &x){return std::conj(x);}
conj(const complex<float> & x)57 static complex<float> conj(const complex<float> &x){return std::conj(x);}
58
59 typedef TypeIndex I;
60 typedef TypeScalaire R;
61 typedef uint64_t uniquecodeInt;
62 static const int type_isdeleted=-1,type_HM=0,type_COO=1, type_CSR=2,type_CSC=3;
63 static const int unsorted=0, sorted_ij=type_CSR,sorted_ji=type_CSC;
64 typedef size_t iterator;
65 typedef size_t Hash;
66 typedef pair<size_t,size_t> Pair;
67 size_t nnz,nnzmax,nhash;
68 mutable size_t nbcollision,nbfind;
69 mutable double matmulcpu;
70 I * i,*j;
71 I *p;
72 R * aij;
73 size_t * head;
74 size_t * next;
75 int half;
76 int state,type_state;
77 size_t nbsort;
78 I sizep;
79 int lock;
80 int fortran; // index start a one ..
81 mutable int re_do_numerics,re_do_symbolic;
82 static const size_t empty= (I) -1;
83 // for Dirichlet BC
84 double tgv;
85 I ntgv;
86
NbCoef() const87 I NbCoef() const {return (I) nnz;}
setcoef(const KN_<R> & x)88 void setcoef(const KN_<R> & x){KN_<R>c(this->aij,nnz); ffassert(x.SameShape(c));
89 if( x.constant()) c=x[0]; else c = x;}
getcoef(KN_<R> & x) const90 void getcoef( KN_<R> & x) const {ffassert(x.N()==(I) nnz);x =KN_<R>(this->aij,nnz);}
91
92 void setdiag(const KN_<R> & d);
93 void getdiag( KN_<R> & d) const;
94 R pscal(R *x,R *y,I sx=1,I sy=1);
pscal(const KN_<R> & x,const KN_<R> & y)95 R pscal(const KN_<R> & x,const KN_<R> & y) { return pscal(x,y,(I) x.step,(I) y.step);}
96 void SetMorse();
97 void UnSetMorse();
98 uniquecodeInt CodeIJ() const ;
99 void init(I nn,I mm,size_t nnnz,int halff);
100 HashMatrix(I nn,I mm,I nnnz,int halff);
101 HashMatrix(istream & f,int cas=-1);
102 HashMatrix(KNM_<R> F,double threshold=1e-30);
103 HashMatrix(int Half,I nn);
104 HashMatrix(const HashMatrix& A);
105 HashMatrix(I nn,const R *diag);
106 void RenumberingInv(KN_<I> II,KN_<I> JJ);
107 void Renumbering(I nn,I mm,KN_<I> II,KN_<I> JJ);
108 void RemoveDoubleij(int kk);// remove
cast_funct(K x)109 template<class R,class K> static R cast_funct(K x) { return (R) x;}
110
111 template<class J,class K> HashMatrix(const HashMatrix<J,K> &A , R (*ff)(K) );
112 template<class J> HashMatrix(const HashMatrix<J,R> &A );
113
operator =(const HashMatrix<II,RR> & A)114 template<class II,class RR> HashMatrix & operator=(const HashMatrix<II,RR>& A )
115 {
116 if( (const void*) this == (const void*) & A) return *this;
117
118 set(A.n,A.m,A.half,A.nnz,A.i,A.j,A.aij,A.fortran,cast_funct<R,RR>);
119 return *this;
120 }
121
122 template<class II,class RR> HashMatrix & operator+=(const HashMatrix<II,RR>& A );
123
124
125 int IsTrianglulare() const ;
126
CheckUnLock(const char * cmm)127 void CheckUnLock(const char * cmm)
128 {
129 if( lock)
130 {
131 std::cerr << " Sorry Forbidder operation ( "
132 << cmm<< " ) matix is lock " << endl;
133 assert(0);
134 }
135 }
136 void setp(I sp);
resize(I nn,I mm=0)137 void resize(I nn, I mm=0) {resize(nn,mm,nnz); }
138
139 void resize(I nn, I mm,size_t nnnz, double tol = -1., int sym=0 );
140 void SymmetrizePattern(); // To do for Suzuki , Paradiso
141 void clear();
hash(size_t ii,size_t jj) const142 Hash hash(size_t ii,size_t jj) const{ return ( (ii-fortran)+ (jj-fortran)*this->n )%nhash; }
143
144 void setfortran(int yes);
145
find(size_t ii,size_t jj) const146 size_t find(size_t ii,size_t jj) const { return find(ii,jj, hash(ii,jj)); }
147
find(I ii,I jj,Hash h) const148 size_t find(I ii,I jj,Hash h) const
149 {
150 nbfind++;
151 for (size_t k=head[h];k!=empty;k=next[k])
152 {
153 ++nbcollision;
154 if( ii== i[k] && jj==j[k] ) return k;
155 }
156 return empty;
157 }
158
insert(I ii,I jj,const R & aa)159 size_t insert(I ii, I jj,const R & aa)
160 {
161 state=unsorted;
162 Hash h= hash(ii,jj);
163 size_t k=find(ii,jj,h);
164 if(k==empty)
165 k =simpleinsert(ii,jj,h);
166 aij[k] += aa;
167 return k;
168 }
169
170 template<typename T> static void HMresize(T *&t,size_t no,size_t nn);
171 template<typename T,typename TT> static void HMcopy( T *dst,const TT *from, size_t nn);
172 template<typename T,typename TT> static void HMcopy( T *dst,const TT *from, size_t nn, T (*ff)(TT) );
173
174 bool do2Triangular(bool lower) ; // put half tp lower or upper
175 void dotranspose();
176 void Increaze(size_t nnznew=0,size_t newnnz=0);// newnnz<0 => newnnz is set to nnz (change value of nnz)
177 void ReHash();
size() const178 size_t size() const { return nnz;}
179 size_t simpleinsert(I ii, I jj,Hash &h);
180 R *pij(I ii,I jj) const;
181 R *npij(I ii,I jj); // with add if no term ii,jj
182
diag(I ii)183 R & diag(I ii) { return operator()(ii,ii);}
diag(I ii) const184 R diag(I ii) const { return operator()(ii,ii);}
operator ()(I ii,I jj) const185 R operator()(I ii,I jj) const
186 {
187 Hash h = hash(ii,jj);
188 size_t k = find(ii,jj,h);
189 if(k==empty) return R();
190 else return aij[k];
191 }
192
operator ()(I ii,I jj)193 R & operator()(I ii,I jj)
194 {
195 return *npij(ii,jj);
196 }
197
operator ()(pair<I,I> ij) const198 R operator()(pair<I,I> ij) const {return operator()(ij.first,ij.second);}
operator ()(pair<I,I> ij)199 R & operator()(pair<I,I> ij) {return operator()(ij.first,ij.second);}
operator [](pair<I,I> ij) const200 R operator[](pair<I,I> ij) const {return operator()(ij.first,ij.second);}
operator [](pair<I,I> ij)201 R & operator[](pair<I,I> ij) {return operator()(ij.first,ij.second);}
202 ~HashMatrix();
203
204
205 void Sortij();
206 void Sortji();
207
208 void set(I nn,I mm,int hhalf,size_t nnnz, I *ii, I *jj, R *aa,int f77=0,int tcsr=0);
209 template<class II> void set(II nn,II mm,int hhalf,size_t nnnz, II *ii, II*jj, R *aa,int f77);
210 template<class II,class RR> void set(II nn,II mm,int hhalf,size_t nnnz, II *ii, II*jj, RR *aa,int f77,R (*ff)(RR));
211
212 void Add(const HashMatrix<I,R> *PA,R coef=R(1),bool trans=false, I ii00=0,I jj00=0);
213
214 HashMatrix &operator=(const HashMatrix &A) ;
215 HashMatrix &operator+=(const HashMatrix &A) ;// {Add(&A); return *this;};
216 HashMatrix &operator-=(const HashMatrix &A) ; //{Add(&A,R(-1.)); return *this;}
217
218
219
220 void operator*=(R v);
221 void operator=(const R & v);
222 void HM();// un sorted ... Default type ..
223 void COO();
224 void COO(I *& IA, I *& IJ, R *& A);
225 void CSR(I *& IA, I *& JA, R *& A);
226 void CSR();
227 void CSC();
228 void CSC(I *& JA, I *& IA, R *& A);
addstateLU(int U)229 static int addstateLU(int U) { return U>0 ? 4 : 5; };
230
231 size_t SortLU(int U);
232 size_t CSC_U(I *& JA, I *& IA, R *& A);
233 size_t CSR_L(I *& IA, I *& JA, R *& A);
234 void Buildp(I nn,I * IA,int type_m,size_t nnzz=0);
Row(I ii)235 Pair Row(I ii) { return RoworCol(ii, true /*!trans*/);}
Col(I jj)236 Pair Col(I jj) { return RoworCol(jj, false /*trans*/);}
237 Pair RoworCol(I ii,bool row);
checksize(size_t nn,size_t mm=0) const238 void checksize(size_t nn,size_t mm=0) const
239 { mm= mm ? mm : nn; assert( (nn ==this->n) && (mm==this->m));}
240 R* addMatMul(R *x,R*Ax,bool Transpose,I sx=1,I sAx=1) const;
addMatTransMul(R * x,R * Ax) const241 R* addMatTransMul(R *x,R*Ax) const { return addMatMul(x,Ax,true); }
addMatMul(R * x,R * Ax) const242 R* addMatMul(R *x,R*Ax) const { return addMatMul(x,Ax,false);}
243 R trace () const;
244 double FrobeniusNorm() const;
245 double norm1() const;
246 double norminfty() const;
sym() const247 bool sym() const {return (half > 0);}
248
typemat() const249 int typemat() const { return int(half > 0)*VirtualMatrix<int,R>::TS_SYM ;}
250 void SetBC(char *wbc,double ttgv);
251
252
253 void addMap(R coef,std::map< pair<I,I>, R> &mij,bool trans=false,I ii00=0,I jj00=0,bool cnj=false,double threshold=0.);
254 bool addMatTo(R coef,HashMatrix<I,R> & mij,bool trans=false,I ii00=0,I jj00=0,bool cnj=false,double threshold=0.,const bool keepSym=false) ;
255
256
257 VirtualMatrix<I,R> & operator +=(MatriceElementaire<R> & me) ;
dump(ostream & f) const258 ostream& dump (ostream&f) const { return f<<*this;}
SetBC(I ii,double ttgv)259 void SetBC(I ii,double ttgv) { diag(ii)=ttgv;};
260
261 double gettgv(I * pntgv=0,double ratio=1e6) const ;
GetReDoNumerics() const262 bool GetReDoNumerics() const { bool b=re_do_numerics; re_do_numerics=0;return b;}
GetReDoSymbolic() const263 bool GetReDoSymbolic() const { bool b=re_do_symbolic; re_do_symbolic=0;return b;}
264
265
toMatriceMorse(bool transpose=false,bool copy=false) const266 HashMatrix<I, R> *toMatriceMorse(bool transpose=false,bool copy=false) const {ffassert(0); return 0;}
psor(KN_<R> & x,const KN_<R> & gmin,const KN_<R> & gmax,double omega)267 double psor(KN_<R> & x,const KN_<R> & gmin,const KN_<R> & gmax , double omega) {ffassert(0); };
268
269 void UnHalf();
Half(int wsym)270 void Half(int wsym) {resize(this->n,this->m,nnz,-1,wsym);}
271 void RemoveHalf(int cas,double tol=-1) ;
272
273 void setsdp(int sym,bool dp); // set of unset to sym / defpos or not
274
ChecknbLine(I n) const275 virtual bool ChecknbLine (I n) const {return this->n==n;}
ChecknbColumn(I m) const276 virtual bool ChecknbColumn (I m) const {return this->m==m;}
CPUsecond()277 static double CPUsecond() {
278 return (double)clock()/CLOCKS_PER_SEC;
279 }
280 };
281 // 0 good , -1 delete, ...
GoodPtrHashMatrix(const HashMatrix<I,R> * p)282 template<class I,class R> int GoodPtrHashMatrix(const HashMatrix<I,R> *p ) {
283 if( p==0) return 1;
284 if( p->N != p->n) return -2;
285 if( p->M != p->m) return -3;
286 if (p->nnz ==-1234567802) return -4;
287 if( p->i && p->j && p->aij ) return 0;
288 return -5;
289 }
CheckPtrHashMatrix(const HashMatrix<I,R> * p,const char * where)290 template<class I,class R> void CheckPtrHashMatrix(const HashMatrix<I,R> *p,const char * where )
291 {
292 int gm=GoodPtrHashMatrix(p);
293 if( gm !=0)
294 {
295 if(gm <0)
296 cout << " n = " << p->n << " == " << p->N
297 << " , m= " << p->m << " "<< p->M
298 << " nzz "<< p->nnz << endl;
299 cerr << " Fatal Error " << where << " invalide HashMatrix Ptr "<< gm << " "<< p << endl;
300 ffassert(0);
301 }
302 }
303
304 // END OF CLASS HashMatrix
305
306 template<class I,class R>
simpleinsert(I ii,I jj,Hash & h)307 inline size_t HashMatrix<I,R>::simpleinsert(I ii, I jj,Hash &h)
308 {
309 state=unsorted;
310 re_do_numerics=1;
311 re_do_symbolic=1;
312
313 type_state=type_HM;
314 if(nnz==nnzmax) {
315 Increaze();
316 h = hash(ii,jj);
317 }
318 i[nnz] = ii;
319 j[nnz] = jj;
320 aij[nnz] = R();
321 next[nnz]=head[h];
322 head[h]=nnz;
323 return nnz++;
324 }
325 template<class I,class R>
pij(I ii,I jj) const326 inline R * HashMatrix<I,R>::pij(I ii,I jj) const
327 {
328 re_do_numerics=1;
329 Hash h = hash(ii,jj);
330 size_t k = find(ii,jj,h);
331 return k==empty ? 0 : aij+k;
332 }
333
334 template<class I,class R>
npij(I ii,I jj)335 inline R *HashMatrix<I,R>::npij(I ii,I jj) // with add if no term ii,jj
336 {
337 re_do_numerics=1;
338 Hash h = hash(ii,jj);
339 size_t k = find(ii,jj,h);
340 if(k==empty)
341 {
342 k =simpleinsert(ii,jj,h);
343 aij[k]=0;
344 }
345 return aij+k;
346 }
347
348 //
349
350
351 template<class I,class RA,class RB=RA,class RAB=RA>
AddMul(HashMatrix<I,RAB> & AB,HashMatrix<I,RA> & A,HashMatrix<I,RB> & B,bool ta=false,bool tb=false,R c=R (1))352 void AddMul(HashMatrix<I,RAB> &AB,HashMatrix<I,RA> &A, HashMatrix<I,RB> &B,bool ta=false,bool tb=false,R c=R(1))
353 {
354 AB.half=false;
355 int An= A.n, Am =A.m;
356 int Bn= B.n, Bm =B.m;
357 int halfA = A.half;
358 int halfB = B.half;
359
360 if(ta) swap(An,Am);
361 if(tb) swap(Bn,Bm);
362 bool tcb = (std::is_same<RB,complex<double> >::value|| std::is_same<RB,complex<float> >::value ) && tb;
363 bool tca = (std::is_same<RA,complex<double> >::value|| std::is_same<RA,complex<float> >::value ) && ta;
364 AB.checksize(An,Bm);
365 ffassert(Am == Bn);
366 // need A col sort , b row sort
367 if( tb)
368 B.CSC(); // sort by COL nd build p.
369 else
370 B.CSR(); // sort by row... and build p.
371 int * Bj = tb ? B.i : B.j;
372 int * Bi = tb ? B.j : B.i;
373 // first Half
374 for(size_t l=0; l< A.nnz;++l)
375 {
376 I i=A.i[l],j=A.j[l];
377 RA aij=A.aij[l];
378 if(ta) swap(i,j);
379 if(tca) aij=HashMatrix<I,RA>::conj(aij);
380
381 for(size_t ll=B.p[j]; ll< B.p[j+1] ;++ll)
382 {
383 I k = Bj[ll];
384 if(verbosity>1000000000) cout << " *** " << i<< " " << " " << k << " : " << j << " : "
385 << ll << " " << B.i[ll] <<" " << B.j[ll]<< " :: " << A.aij[l]*B.aij[ll] <<endl;
386 assert(j == Bi[ll]);
387 RB bjk = tcb ? HashMatrix<I,RB>::conj(B.aij[ll]) : B.aij[ll];
388
389 AB(i,k) += c* aij*bjk;
390 }
391 }
392 // second Half of A
393 if( halfA )
394 {
395 if(verbosity>10) cout<< " ** halfA "<< endl;
396 for(size_t l=0; l< A.nnz;++l)
397 {
398 I i=A.i[l],j=A.j[l];
399 if( i==j) continue;
400 RA aij=A.aij[l];
401 if(!ta) swap(i,j);
402 //if(!tca) aij=HashMatrix<I,RA>::conj(aij);
403
404 for(size_t ll=B.p[j]; ll< B.p[j+1] ;++ll)
405 {
406 I k = Bj[ll];
407 if(verbosity>1000000000) cout << " *** " << i<< " " << " " << k << " : " << j << " : "
408 << ll << " " << B.i[ll] <<" " << B.j[ll]<< " :: " << A.aij[l]*B.aij[ll] <<endl;
409 assert(j == Bi[ll]);
410 RB bjk = tcb ? HashMatrix<I,RB>::conj(B.aij[ll]) : B.aij[ll];
411
412 AB(i,k) += c* aij*bjk;
413 }
414 }
415 }
416 if( halfB)
417 {
418 if( !tb)
419 B.CSC(); // sort by COL nd build p.
420 else
421 B.CSR(); // sort by row... and build p.
422 // first Half miss ..
423 for(size_t l=0; l< A.nnz;++l)
424 {
425 I i=A.i[l],j=A.j[l];
426 RA aij=A.aij[l];
427 if(ta) swap(i,j);
428 if(tca) aij=HashMatrix<I,RA>::conj(aij);
429
430 for(size_t ll=B.p[j]; ll< B.p[j+1] ;++ll)
431 {
432 if( Bi[ll]== Bj[ll]) continue;
433 I k = Bi[ll];
434
435 assert(j == Bj[ll]);
436 RB bjk = B.aij[ll];
437 if(tb) bjk = HashMatrix<I,RB>::conj(bjk);
438 AB(i,k) += c* aij*bjk;
439 }
440 }
441 // second Half of A
442 if( halfA )
443 for(size_t l=0; l< A.nnz;++l)
444 {
445 I i=A.i[l],j=A.j[l];
446 if( i==j) continue;
447 RA aij=A.aij[l];
448 if(!ta) swap(i,j);
449 //if(!tca) aij=HashMatrix<I,RA>::conj(aij);
450
451 for(size_t ll=B.p[j]; ll< B.p[j+1] ;++ll)
452 {
453 if( Bi[ll]== Bj[ll]) continue;
454 I k = Bi[ll];
455 if(verbosity>1000000000) cout << " *** " << i<< " " << " " << k << " : " << j << " : "
456 << ll << " " << B.i[ll] <<" " << B.j[ll]<< " :: " << A.aij[l]*B.aij[ll] <<endl;
457 assert(j == Bj[ll]);
458 RB bjk = tcb ? HashMatrix<I,RB>::conj(B.aij[ll]) : B.aij[ll];
459
460 AB(i,k) += c* aij*bjk;
461 }
462 }
463
464
465 }
466
467 }
468
469 template<class I,class R>
operator <<(std::ostream & f,const HashMatrix<I,R> & A)470 std::ostream & operator<<(std::ostream & f, const HashMatrix<I,R> &A)
471 {
472 int p20=20;
473 long pold= f.precision();
474 if( pold > 20) p20= (int) pold;
475 if(A.type_state==HashMatrix<I,R>::type_CSR)
476 {
477 using namespace std;
478 f << "# Sparse Matrix (Morse) " << &A << endl;
479 f << "# first line: n m (is symmetic) nnz \n";
480 f << "# after for each nonzero coefficient: i j a_ij where (i,j) \\in {1,...,n}x{1,...,m} \n";
481
482 f << A.n << " " << A.m << " " << A.half << " " << A.nnz <<endl;
483 I k=A.p[0];
484
485 for (I i=0;i<A.n;i++)
486 {
487 I ke=A.p[i+1];
488 for (;k<ke;k++)
489 f << setw(9) << i+1 << ' ' << setw(9) << A.j[k]+1 << ' ' << setprecision( p20) << RNM::removeeps(A.aij[k])<< '\n' ;
490
491 }
492 }
493 else
494 {
495 f << "# HashMatrix Matrix (COO) "<< &A << endl;
496 f << "# n m nnz half fortran state \n";
497 f << A.n << " " << A.m << " " << A.nnz << " "<< A.half << " " << A.fortran
498 << " " << A.state<< " " << A.type_state<< " " << endl;
499 for(size_t k=0; k < A.nnz; ++k)
500 f << setw(10) << A.i[k] << setw(10) << A.j[k] << ' '<< setprecision( p20) << A.aij[k] << endl;
501 }
502 f.precision(pold);
503 return f;
504 }
505
506 template<class R>
BuildCombMat(HashMatrix<int,R> & mij,const list<tuple<R,VirtualMatrix<int,R> *,bool>> & lM,bool trans,int ii00,int jj00,bool cnj)507 tuple<int,int,bool> BuildCombMat(HashMatrix<int,R> & mij,const list<tuple<R,VirtualMatrix<int,R>*,bool> > &lM,bool trans,int ii00,int jj00,bool cnj)
508 {
509
510 typedef typename list<tuple<R,VirtualMatrix<int,R> *,bool> >::const_iterator lconst_iterator;
511
512 lconst_iterator begin=lM.begin();
513 lconst_iterator end=lM.end();
514 lconst_iterator i;
515
516
517 int n=0,m=0;
518 bool sym=true;
519 for(i=begin;i!=end&&sym;i++)
520 {
521 if(get<1>(*i))// M == 0 => zero matrix
522 {
523 VirtualMatrix<int,R>& M=*get<1>(*i);
524 if(!M.sym())
525 sym = false;
526 }
527 }
528 int iter=0;
529 for(i=begin;i!=end;i++)
530 {
531 if(get<1>(*i)) // M == 0 => zero matrix
532 {
533 VirtualMatrix<int,R> & M=*get<1>(*i);
534 bool transpose = get<2>(*i) != trans;
535 bool conjuge = get<2>(*i) != cnj;
536
537 ffassert( &M);
538 R coef= get<0>(*i);
539 if(verbosity>99)
540 cout << " "<< iter++<< " BuildCombMat + " << coef << "*" << &M << " " << sym << " t = " << transpose << " " << get<2>(*i) << endl;
541 { if(transpose) {m=max(m,M.n); n=max(n,M.m);} else{n=max(M.n,n); m=max(M.m,m);}}
542
543 M.addMatTo(coef,mij,transpose,ii00,jj00,conjuge,0.0,sym);
544 }
545 }
546
547 //V4 return new MatriceMorseOld<R>(n,m,mij,sym);
548 return make_tuple(n,m,sym);
549 }
550
551 template<class R>
nmCombMat(const list<tuple<R,VirtualMatrix<int,R> *,bool>> & lM,bool trans,int ii00,int jj00,bool cnj=false)552 tuple<int,int,bool> nmCombMat(const list<tuple<R,VirtualMatrix<int,R>*,bool> > &lM,bool trans,int ii00,int jj00,bool cnj=false)
553 {
554
555 typedef typename list<tuple<R,VirtualMatrix<int,R> *,bool> >::const_iterator lconst_iterator;
556
557 lconst_iterator begin=lM.begin();
558 lconst_iterator end=lM.end();
559 lconst_iterator i;
560
561
562 int n=0,m=0;
563 bool sym=true;
564 for(i=begin;i!=end&&sym;i++++)
565 {
566 if(std::get<1>(*i)) // M == 0 => zero matrix
567 {
568 VirtualMatrix<int,R>& M=*std::get<1>(*i);
569 if(!M.sym())
570 sym = false;
571 }
572 }
573
574 for(i=begin;i!=end;i++++)
575 {
576 if(std::get<1>(*i)) // M == 0 => zero matrix
577 {
578 VirtualMatrix<int,R>& M=*std::get<1>(*i);
579 bool transpose = std::get<2>(*i) != trans;
580 ffassert( &M);
581 R coef=std::get<0>(*i);
582 if(verbosity>99)
583 cout << " BuildCombMat + " << coef << "*" << &M << " " << sym << " t = " << transpose << " " << std::get<2>(*i) << endl;
584 { if(transpose) {m=max(m,M.n); n=max(n,M.m);} else{n=max(M.n,n); m=max(M.m,m);}}
585
586 }
587 }
588
589 return make_tuple(n,m,sym);
590 }
591
592 template<class R>
BuildCombMat(const list<tuple<R,VirtualMatrix<int,R> *,bool>> & lM,bool trans=false,int ii00=0,int jj00=0)593 HashMatrix<int,R>* BuildCombMat(const list<tuple<R,VirtualMatrix<int,R>*,bool> > &lM,bool trans=false,int ii00=0,int jj00=0)
594 {
595
596 auto nmsym=nmCombMat(lM,trans,ii00,jj00);
597 int n = std::get<0>(nmsym), m =std::get<1>(nmsym);
598 bool half= std::get<2>(nmsym);
599 HashMatrix<int,R> * mij= new HashMatrix<int,R>(n,m,0,int(half));
600 nmsym=BuildCombMat(*mij,lM,trans,ii00,jj00,trans);// remember trans => conj
601
602 return mij; // V4 mij;
603
604 }
605
606 template<class I>
ij_mat(bool trans,I ii00,I jj00,I i,I j)607 static inline pair<I,I> ij_mat(bool trans,I ii00,I jj00,I i,I j) {
608 // warning trans sub matrix and not the block.
609 return trans ? make_pair<I,I>(j+ii00,i+jj00)
610 : make_pair<I,I>(i+ii00,j+jj00) ; }
611
612
613
614 template<class I,class R> template<typename T,typename TT>
HMcopy(T * dst,const TT * from,size_t nn,T (* ff)(TT))615 void HashMatrix<I,R>::HMcopy( T *dst,const TT *from, size_t nn, T(*ff)(TT))
616 {
617 for(size_t i=0; i< nn; ++i)
618 dst[i]= ff(from[i]);
619 }
620 template<class I,class R> template<typename T,typename TT>
HMcopy(T * dst,const TT * from,size_t nn)621 void HashMatrix<I,R>::HMcopy( T *dst,const TT *from, size_t nn)
622 {
623 for(size_t i=0; i< nn; ++i)
624 dst[i]= (T) from[i];
625 }
626 template<class I,class R> template<class II,class RR> HashMatrix<I,R> &
operator +=(const HashMatrix<II,RR> & A)627 HashMatrix<I,R>::operator+=(const HashMatrix<II,RR>& A )
628 {
629 Addto(this,&A,cast_funct);
630 }
631
632
633 template<class I,class R> template<class II,class RR>
set(II nn,II mm,int hhalf,size_t nnnz,II * ii,II * jj,RR * aa,int f77,R (* ff)(RR))634 void HashMatrix<I,R>::set(II nn,II mm,int hhalf,size_t nnnz, II *ii, II*jj, RR *aa,int f77,R(*ff)(RR))
635 {
636 clear();
637 this->n=nn;
638 this->m=mm;
639 this->N=nn;
640 this->M=mm;
641 fortran=f77;
642 half=hhalf;
643 Increaze(nnnz);
644 nnz=nnnz;
645
646 HMcopy(i,ii,nnnz);
647 HMcopy(j,jj,nnnz);
648 HMcopy(aij,aa,nnnz,ff);
649 ReHash();
650 }
651 template<class I,class R> template<class II>
set(II nn,II mm,int hhalf,size_t nnnz,II * ii,II * jj,R * aa,int f77)652 void HashMatrix<I,R>::set(II nn,II mm,int hhalf,size_t nnnz, II *ii, II*jj, R *aa,int f77)
653 {
654 clear();
655 this->n=nn;
656 this->m=mm;
657 this->N=nn;
658 this->M=mm;
659 fortran=f77;
660 half=hhalf;
661 Increaze(nnnz);
662 nnz=nnnz;
663
664 HMcopy(i,ii,nnnz);
665 HMcopy(j,jj,nnnz);
666 HMcopy(aij,aa,nnnz);
667 ReHash();
668 }
669 #endif
670