1 
2 #include "HashMatrix.hpp"
3 //   Warning the instation will be don a the add Of the File
4 // F. hecht
5 // version v0 ...
6 
WhichMatrix(istream & f)7 int WhichMatrix(istream & f)
8 {
9     string line;
10     while ( isspace(f.peek()))
11         f.get();
12     if  ( f.peek() =='#' )
13     {
14         line="";
15         while ( f.good()  )
16         {
17             char c=f.get();
18             if(c=='\n' || c=='\r') { break;}
19             line += c;
20         }
21         if( line.find("(Morse)") != std::string::npos)
22             return 2; // morse
23 
24         else  if( line.find("(COO)") != std::string::npos)
25             return 3; // morse
26 
27         return 0;
28     }
29     return 0;
30 }
31 
32 
33 template<class I,class R>
HeapSort(I * ii,I * jj,R * aij,long n)34 void  HashMatrix<I,R>::HeapSort(I *ii,I *jj, R *aij,long n)
35 {
36     long l,j,r,i;
37     I criti,critj;
38     R crita;
39 #define HM__criteq(i) criti=ii[i],critj=jj[i], crita=aij[i]
40 #define HM__eqcrit(i) ii[i]=criti,jj[i]=critj, aij[i]=crita
41 #define HM__eqij(i,j) ii[i]=ii[j],jj[i]=jj[j], aij[i]=aij[j]
42 #define HM__cmpij(i,j) ( ii[i] != ii[j] ? ii[i] < ii[j] : jj[i] < jj[j])
43 #define HM__cmpcrit(j) ( criti != ii[j] ? criti < ii[j] : critj < jj[j])
44 
45     if( n <= 1) return;
46     l = n/2 + 1;
47     r = n;
48     while (1) { // label 2
49         if(l <= 1 ) { // label 20
50             HM__criteq(r-1);//crit = c[r];
51             HM__eqij(r-1,0);
52             --r; //c[r--] = c[1];
53             if ( r == 1 ) { HM__eqcrit(0); /*c[1]=crit;*/ return;}
54         } else  {--l;  HM__criteq(l-1);}// crit = c[--l];
55         j=l;
56         while (1) {// label 4
57             i=j;
58             j=2*j;
59             if  (j>r) {HM__eqcrit(i-1);/*c[i]=crit;*/break;} // L8 -> G2
60             if ((j<r) && (HM__cmpij(j-1,j) )) j++; // L5
61             if (HM__cmpcrit(j-1)) HM__eqij(i-1,j-1);//c[i]=c[j]; // L6+1 G4
62             else {HM__eqcrit(i-1); /*c[i]=crit;*/break;} //L8 -> G2
63         }
64     }
65 #undef HM__criteq
66 #undef HM__eqcrit
67 #undef HM__eqij
68 #undef HM__cmpij
69 #undef HM__cmpcrit
70 }
71 
72 template<class I,class R>
setdiag(const KN_<R> & d)73 void HashMatrix<I,R>::setdiag(const KN_<R> & d)
74 {
75     for(int ii=0; ii<this->n; ++ii)
76         diag(ii)=d[ii];
77 }
78 
79 template<class I,class R>
getdiag(KN_<R> & d) const80 void HashMatrix<I,R>::getdiag( KN_<R> & d) const
81 {
82     for(int ii=0; ii<this->n; ++ii)
83         d[ii]=diag(ii);
84 }
85 
86 template<class I,class R>
pscal(R * x,R * y,I sx,I sy)87 R HashMatrix<I,R>::pscal(R *x,R *y,I sx,I sy)
88 {
89     R s=0;
90     if( fortran)
91     {
92         x -= fortran;
93         y -= fortran;
94     }
95     if(half > 0)
96         for (size_t k=0; k<nnz;++i)
97         {
98             if(i[k] != j[k] )
99                 s += conj(x[i[k]*sx])*aij[k]*y[j[k]*sy] + conj(x[j[k]*sx])* (half == 1?conj(aij[k]):aij[k]) *y[i[k]*sy] ;
100             else
101                 s += conj(x[i[k]*sx])*aij[k]*y[j[k]*sy];
102         }
103     else
104         for (size_t k=0; k<nnz;++i)
105             s += conj(x[i[k]*sx])*aij[k]*y[j[k]*sy];
106     return s;
107 }
108 
109 
110 template<class I,class R>
CodeIJ() const111 typename   HashMatrix<I,R>::uniquecodeInt HashMatrix<I,R>::CodeIJ() const  {
112     uniquecodeInt code=this->n;
113     code ^=roll64(this->n,0);
114     code ^=roll64(this->m,24);
115     code ^=roll64(nnz, 48);
116     for(size_t k=0,kk=0; k< nnz;++k)
117     {
118         code^=roll64(i[k],++kk);
119         code^=roll64(j[k],++kk);
120     }
121     return code;
122 }
123 
124 template<class I,class R>
HashMatrix(I nn,I mm,I nnnz,int halff)125 HashMatrix<I,R>::HashMatrix(I nn,I mm,I nnnz,int halff)
126 :  VirtualMatrix<I,R>(nn,mm),  nnz(0),nnzmax(0),nhash(0),nbcollision(0),nbfind(0),matmulcpu(0.),i(0),j(0),p(0),aij(0),
127 head(0), next(0),
128 // trans(false),
129 half(halff), state(unsorted),type_state(type_HM),
130 nbsort(0),sizep(0),lock(0), fortran(0) ,
131 re_do_numerics(0),re_do_symbolic(0)
132 {
133     Increaze(nnnz);
134 }
135 
136 template<class I,class R>
HashMatrix(istream & f,int cas)137 HashMatrix<I,R>::HashMatrix(istream & f,int cas):
138 VirtualMatrix<I,R>(0,0),  nnz(0),nnzmax(0),nhash(0),nbcollision(0),nbfind(0),matmulcpu(0.),i(0),j(0),p(0),aij(0),
139 head(0), next(0),
140 // trans(false),
141 half(0), state(unsorted),type_state(type_HM),
142 nbsort(0),sizep(0),lock(0), fortran(0) ,
143 re_do_numerics(0),re_do_symbolic(0)
144 {
145     if(cas ==-1) cas=  WhichMatrix(f)  ;
146     // eat lines with #
147     string line;
148     int k=0;
149     while ( isspace(f.peek()))
150         f.get();
151     while ( f.peek() =='#' )
152     {
153         line="";
154         while ( f.good()  )
155         {
156             char c=f.get();
157             if(c=='\n' || c=='\r') { break;}
158             line += c;
159         }
160         if( f.peek()=='\n' || f.peek()=='\r') f.get();
161         if(verbosity>9)
162             cout << "   --  Matrx: "<< k << " :"   << line << " cas " << cas << endl;
163         k++;
164     }
165     if(cas== 2)
166     {
167         I rn,rm;
168         int rsymetrique;
169         size_t rnbcoef;
170         ffassert(f.good() );
171         f >> rn >> rm >> rsymetrique >>rnbcoef;
172 
173         if(verbosity>3)
174             cout << "     -- Read Mat: " <<  this->n << " x " <<  this->m << " sym : " << rsymetrique << " nnz=" << rnbcoef <<endl;
175         ffassert(f.good() && rn>=0 && rm>=0 && rnbcoef>=0 );
176         resize(rn,rm,rnbcoef);
177         half =rsymetrique;
178         I ii,jj;
179 
180         R aaij;
181 
182         for (size_t kk =0;kk<rnbcoef; ++kk)
183         {
184             f >> ii >> jj >> aaij;
185             ffassert(f.good() );
186             operator()(ii-1,jj-1) = aaij;
187         }
188     }
189     else if(cas== 3)
190     {
191         I nn,mm,f77, sstate, tstate;
192         int hhalf;
193         size_t nnnz;
194         f >> nn >> mm >> nnnz >> hhalf >> f77 >> sstate >> tstate ;
195         ffassert(f.good() && nn>0 && mm>0 && nnnz>0  );
196         resize(nn,mm,nnnz);
197         half =hhalf;
198         for(size_t kk=0; kk < nnnz; ++kk)
199         {
200             I ii,jj;
201             R aaij;
202             f >>  ii >> jj >> aaij ;
203             ffassert(f.good() );
204             operator()(ii-f77,jj-f77) = aaij;
205         }
206 
207     }
208     else {
209         cerr << " Unkown matrix type" << endl << endl;
210         ffassert(0);
211     }
212 
213 
214 }
215 
216 template<class I,class R>
HashMatrix(KNM_<R> F,double threshold)217 HashMatrix<I,R>::HashMatrix(KNM_<R> F,double  threshold)
218 : VirtualMatrix<I,R>(F.N(),F.M()),  nnz(0),nnzmax(0),nhash(0),nbcollision(0),nbfind(0),matmulcpu(0.),i(0),j(0),p(0),aij(0),
219 head(0), next(0),
220 half(0), state(unsorted),type_state(type_HM),
221 nbsort(0),sizep(0),lock(0), fortran(0) ,
222 re_do_numerics(0),re_do_symbolic(0),tgv(0), ntgv(0)
223 
224 {
225     Increaze();
226     for(int ii=0; ii <F.N(); ++ii)
227         for(int jj=0; jj <F.M(); ++jj)
228         {
229             R Fiijj = F(ii,jj);
230             if(abs(Fiijj) > threshold)
231                 operator()(ii,jj)=Fiijj;
232         }
233 }
234 
235 template<class I,class R>
HashMatrix(int Half,I nn)236 HashMatrix<I,R>::HashMatrix(int Half,I nn)
237 : VirtualMatrix<I,R>(nn,nn),  nnz(0),nnzmax(0),nhash(0),nbcollision(0),nbfind(0),matmulcpu(0.),i(0),j(0),p(0),aij(0),
238 head(0), next(0),
239 half(Half), state(unsorted),type_state(type_HM),
240 nbsort(0),sizep(0),lock(0), fortran(0) ,
241 re_do_numerics(0),re_do_symbolic(0)
242 
243 {
244     Increaze(nn*4);
245 }
246 
247 template<class I,class R>
HashMatrix(const HashMatrix & A)248 HashMatrix<I,R>::HashMatrix(const HashMatrix& A)
249 :  VirtualMatrix<I,R> (A),  nnz(0),nnzmax(0),nhash(0),nbcollision(0),nbfind(0),matmulcpu(0.),i(0),j(0),p(0),aij(0),
250 head(0), next(0),
251 half(A.half), state(unsorted),type_state(type_HM),
252 nbsort(0),sizep(0),lock(0), fortran(0) ,
253 re_do_numerics(0),re_do_symbolic(0)
254 {
255     Increaze(A.nnz);
256     operator=(A);
257 }
258 
259 template<class I,class R> template<class II>
HashMatrix(const HashMatrix<II,R> & A)260 HashMatrix<I,R>::HashMatrix(const HashMatrix<II,R>& A)
261 :  VirtualMatrix<I,R> (A.n,A.m),  nnz(0),nnzmax(0),nhash(0),nbcollision(0),nbfind(0),matmulcpu(0.),i(0),j(0),p(0),aij(0),
262 head(0), next(0),
263 half(A.half), state(unsorted),type_state(type_HM),
264 nbsort(0),sizep(0),lock(0), fortran(0) ,
265 re_do_numerics(0),re_do_symbolic(0)
266 {
267     HashMatrix<I,R>::set<II>(A.n,A.m,A.half,A.nnz,A.i,A.j,A.aij,fortran);
268 }
269 
270 template<class I,class R> template<class II,class RR>
HashMatrix(const HashMatrix<II,RR> & A,R (* ff)(RR))271 HashMatrix<I,R>::HashMatrix(const HashMatrix<II,RR>& A, R (*ff)(RR) )
272 :  VirtualMatrix<I,R> (A.n,A.m),  nnz(0),nnzmax(0),nhash(0),nbcollision(0),nbfind(0),matmulcpu(0.),i(0),j(0),p(0),aij(0),
273 head(0), next(0),
274 half(A.half), state(unsorted),type_state(type_HM),
275 nbsort(0),sizep(0),lock(0), fortran(0) ,
276 re_do_numerics(0),re_do_symbolic(0)
277 {
278      HashMatrix<I,R>::set<II,RR>(A.n,A.m,A.half,A.nnz,A.i,A.j,A.aij,fortran,ff);
279 }
280 
281 template<class I,class R>
HashMatrix(I nn,const R * diag)282 HashMatrix<I,R>::HashMatrix(I nn,const R *diag)
283 :  VirtualMatrix<I,R> (nn,nn),  nnz(0),nnzmax(0),nhash(0),nbcollision(0),nbfind(0),matmulcpu(0.),i(0),j(0),p(0),aij(0),
284 head(0), next(0),
285 half(0), state(unsorted),type_state(type_HM),
286 nbsort(0),sizep(0),lock(0), fortran(0) ,
287 re_do_numerics(0),re_do_symbolic(0)
288 {
289     Increaze(nn);
290     for(int k=0;k<nn;++k)
291         operator()(k,k) = diag[k];
292 
293 }
294 
295 
296 template<class I,class R>
setp(I sp)297 void HashMatrix<I,R>::setp(I sp)
298 {
299 
300     if(sp == 0)
301     {
302         delete [] p;
303         p=0;
304     }
305 
306     else if( sp != sizep)
307     {
308         if(p) delete [] p;
309         p = new I[sp];
310 
311     }
312     if(p )
313         for(I ii=0; ii<sp; ++ii)
314             p[ii]=-1;
315      if(verbosity>999) cout << "  HashMatrix:: setp "<< this << " sp= " << sp << " p  "<< p << " / " << (sp != sizep) << endl;
316     ffassert( (sp==0) ==  (p==0)  );
317     sizep=sp;
318 
319 }
320 template<class I,class R>
RemoveHalf(int cas,double tol)321 void  HashMatrix<I,R>::RemoveHalf(int cas,double tol)
322 {
323 
324     size_t kk=0;
325     if( cas <=0) // L cas
326     {
327     for(size_t k=0; k <nnz ;++k)
328         if(  abs(aij[k])> tol && ( i[k] >= j[k] ) )
329         {
330             i[kk] = i[k];
331             j[kk] = j[k];
332             aij[kk] = aij[k];
333             ++kk;
334         }
335     }
336     else
337     {
338      for(size_t k=0; k <nnz ;++k)
339         if(  abs(aij[k])> tol && ( i[k] <= j[k] ) )
340         {
341             i[kk] = i[k];
342             j[kk] = j[k];
343             aij[kk] = aij[k];
344             ++kk;
345         }
346     }
347     half= cas==0;
348     size_t nrt = nnz - kk;
349     if (verbosity>2)
350         cout << " RemoveHalf " << cas << " tol "  << tol << " , remove term : " << nrt << " half =" << half << endl;
351     nnz = kk;
352     if(nrt)  Increaze(kk);
353     else ReHash();
354     state= unsorted;
355     type_state=type_COO;
356 
357 }
358 
359 template<class I,class R>
resize(I nn,I mm,size_t nnnz,double tol,int sym)360 void  HashMatrix<I,R>::resize(I nn, I mm,size_t nnnz, double tol , int sym )
361 {
362 
363     mm= mm ? mm : nn;
364     if( nn == this->n && mm == this->m && nnz == nnnz && sym == half && tol <0) return ;
365     this->m=mm;
366     this->n=nn;
367     this->N=nn;
368     this->M=mm;
369     R mxt =0;
370     size_t kk=0;
371     if( (nn>0) && (mm >0))
372     for(size_t k=0; k <nnz ;++k)
373     {
374         double t =abs(aij[k]);
375         if( i[k] < this->n && j[k] < this->m && t > tol && (!sym || i[k] <= j[k] ) )
376         {
377             i[kk] = i[k];
378             j[kk] = j[k];
379             aij[kk] = aij[k];
380             ++kk;
381         }
382     }
383     if(verbosity>9 ) cout << "HashMatrix<I,R>::resize: new nnz  = " << kk << " " << nnz << " " << " sym " << sym << " "<< tol << " " << nn << " " << mm << endl;
384     half = (half ? half : sym);
385     nnnz = max(nnnz,kk);
386     bool rresize = (nnnz < 0.8*nnz) || ((nnnz > 1.2*nnz)) ;
387     nnz = kk; // forget after value ...
388     if(rresize) Increaze(nnnz);
389     else ReHash();
390     state= unsorted;
391     type_state=type_COO;
392 }
393 template<class I,class R>
SymmetrizePattern()394 void HashMatrix<I,R>::SymmetrizePattern()
395 {
396     ffassert( this->n==this->m);
397     for(size_t k=0; k <nnz ;++k)
398         npij( j[k], i[k]);
399 
400 }
401 
402 template<class I,class R>
clear()403 void HashMatrix<I,R>::clear()
404 {
405     nbsort=0;
406     re_do_numerics=1;
407     re_do_symbolic=1;
408     half=0;
409     lock=0;
410     fortran = 0;
411     state=unsorted;
412     type_state = type_HM;
413     setp(0); // unset
414     if(nnz)
415     {
416         nnz=0;
417         nbcollision=0;
418         for (size_t h=0;h<nhash;++h)
419             head[h]=empty;
420     }
421 
422 }
423 
424 template<class I,class R>
setfortran(int yes)425 void HashMatrix<I,R>::setfortran(int yes)
426 {
427 // trop casse geule
428     int shift =  yes-fortran;
429     if( shift == 0) return ;
430     CheckUnLock("setfortran");
431     for( int k = 0; k<nnz ; ++ k)
432     {
433         i[k] += shift;
434         j[k] += shift;
435     }
436     if (  type_CSC == type_state)
437         for (int k=0; k<= this->m; ++k)
438             p[k] += shift;
439     else  if (  type_CSR == type_state)
440         for (int k=0; k<= this->n; ++k)
441             p[k] += shift;
442     fortran = yes ;
443     // do the shift
444     return ;
445 }
446 
447 
448 
449 
450 template<class I,class R>  template<typename T>
HMresize(T * & t,size_t no,size_t nn)451 void HashMatrix<I,R>::HMresize(T *&t,size_t no,size_t nn)
452 {
453     if( no != nn || t==0)
454     {
455     T * tt= new T[nn];
456     if(t)
457         for(size_t i=0; i< no; ++i)
458            tt[i]= t[i];
459     if(t) delete [] t;
460     t=tt;
461     }
462 }
463 
464 
do2Triangular(bool lower)465 template<class I,class R>   bool  HashMatrix<I,R>::do2Triangular(bool lower)
466 {
467     ffassert(half); //
468     size_t nc=0;
469     if( lower )
470        for(size_t k=0; k< nnz; ++k)
471        {
472            if(  (i[k] < j[k]) ) // Upper
473                ++nc,swap(i[k],j[k]);
474        }
475     else
476         for(size_t k=0; k< nnz; ++k)
477         {
478             if(  (i[k] > j[k]) ) // Lower
479                 ++nc,swap(i[k],j[k]);
480         }
481     if( nc)
482     {
483     ReHash();
484     state=unsorted;
485     }
486     return  nc;
487 }
IsTrianglulare() const488 template<class I,class R>   int HashMatrix<I,R>::IsTrianglulare() const
489 {
490     size_t nU, nL=0,nD=0;
491     for(size_t k=0; k< nnz; ++k)
492     {
493         if( i[k] < j[k] ) ++nU;// 0,10
494         else if( i[k] == j[k] ) ++nD;// 10,10
495         else ++nL;// 10,0
496     }
497     // 3 => no sym
498     // 1 nU==0 =>   Trian Lower (no upper term)
499     // 2 nL==0 =>   Trian Upper (no upper term)
500     return 2*!nL +  !nU ;
501 }
502 template<class I,class R>
dotranspose()503 void HashMatrix<I,R>::dotranspose()
504 {
505     swap(i,j);
506     swap(this->n,this->m);
507     swap(this->N,this->M);
508     conj(aij,this->nnz);
509     ReHash();
510     state=unsorted;
511 }
512 template<class I,class R>
Renumbering(I nn,I mm,KN_<I> ii,KN_<I> jj)513 void HashMatrix<I,R>::Renumbering(I nn,I mm,KN_<I> ii,KN_<I> jj)
514 {
515     //Do  B_ii(i),jj(j) : A_i,j
516     size_t kk=0;
517     for(size_t k=0; k<nnz; ++k)
518     {
519         I i1=ii[i[k]],j1=jj[i[k]];
520         if( i1 >=0 && j1 >=0 && i1 <nn && j1 < mm) // coef existe
521         {
522             i[kk] =i1;
523             j[kk] =j1;
524             aij[kk++]=aij[k];
525         }
526     }
527     nnz = kk;
528     this->n=nn;
529     this->m=mm;
530     this->N=nn;
531     this->M=mm;
532     state=unsorted;
533 
534     RemoveDoubleij(kk); // remove double term
535 
536 
537 
538 }
539 template<class I,class R>
RemoveDoubleij(int kkk)540 void HashMatrix<I,R>::RemoveDoubleij(int kkk)
541 {
542     nnz=kkk;
543     COO();
544     I ip=-1,jp=-1;
545     long kk=-1;
546     for(size_t k=0; k<nnz; ++k)
547     {
548         if( ip != i[k] && jp != j[k] )
549         {  // new term
550             ++kk;
551             i[kk] =ip=i[k];
552             j[kk] =jp=j[k];
553             aij[kk]=aij[k];
554         }
555         else
556             aij[kk]=+aij[k];
557     }
558     if(verbosity>4)
559         cout << " HashMatrix::RemoveDoubleij  remove "<< nnz - kk << " coef "<< endl;
560     Increaze(kk,kk);
561 }
562 template<class I,class R>
RenumberingInv(KN_<I> II,KN_<I> JJ)563 void HashMatrix<I,R>::RenumberingInv(KN_<I> II,KN_<I> JJ)
564 {
565        //Do  B_(i),(j) : A_ii(i),jj(j)
566     I n = this->n, m = this->m;
567     I nn= II.N(), mm= JJ.N();
568     const I minus1 =-1;
569     KN<I> ii(n,minus1), jj(m,minus1);
570     // build inversion
571     int notinjection=0;
572     // build invertion ..
573     for( I l=0; l<nn; ++l)
574     {
575         I IIl = II[l];
576         if( IIl >= 0 && IIl < n)
577         {
578             if(ii[IIl]>=0 ) notinjection =1;
579           ii[IIl]=l;
580         }
581     }
582     for( int l=0; l<mm; ++l)
583     {
584         I IIl = JJ[l];
585         if( IIl >= 0 && IIl < m)
586         {
587             if(jj[IIl]>=0 ) notinjection |=2 ;
588             jj[IIl]=l;
589         }
590     }
591     if( !notinjection)
592     {
593         cerr << " HashMatrix<I,R>::Renumbering not injection " <<notinjection <<endl;
594         ffassert(0); // to do
595     }
596     Renumbering(nn,mm,ii,jj);
597 }
598 
599 template<class I,class R>
Increaze(size_t nnzxnew,size_t newnnz)600 void HashMatrix<I,R>::Increaze(size_t nnzxnew,size_t newnnz)
601 {
602     size_t mnx =(size_t)max(this->n,this->m);
603     if(newnnz==0) newnnz = nnz;
604     else nnz=newnnz;
605     if( nnzxnew==0) {
606         nnzxnew = max(max(mnx*2,newnnz),(size_t) (nnzmax*1.2) );
607     }
608 
609     size_t nzzx =  nnzxnew;
610     double nnzl = min(max(1.,double(nzzx)/mnx),50.);// bof bof !!!! FH
611     size_t nh = mnx*nnzl;
612     if(verbosity>999) cout << "HMresize "<< nzzx << " " <<nnzxnew << " mpirank" << mpirank << " " << this << endl;
613     HMresize(i,nnz,nzzx);
614     HMresize(j,nnz,nzzx);
615     HMresize(aij,nnz,nzzx);
616     HMresize(next,0,nzzx);
617     HMresize(head,0,nh);
618     nnzmax=nzzx;
619     nhash = nh;
620     ReHash();
621     state=unsorted;
622 }
623 
624 template<class I,class R>
ReHash()625 void HashMatrix<I,R>::ReHash()
626 {
627 
628     for(size_t h=0;h<nhash;++h)
629         head[h]= empty;
630     for(size_t k=0;k<nnz;++k)
631     {
632         Hash h= hash(i[k],j[k]);
633         next[k]=head[h];
634         head[h]=k;
635     }
636 }
637 
638 
639 template<class I,class R>
~HashMatrix()640 HashMatrix<I,R>::~HashMatrix()
641 {
642     if(nbfind && verbosity>4)
643         cout << "    ~HashMatrix:   Mean collision in hash: " << (double) nbcollision/ nbfind << " "
644         << this << " rank: "<< mpirank << " matmul " << matmulcpu << "s" << endl;
645 
646     delete [] i;
647     delete [] j;
648     delete [] aij;
649     delete [] next;
650     delete [] head;
651     if(p) delete [] p;
652     type_state=type_isdeleted;// Mark Matrix is type_isdeleted
653     i =0;
654     j=0;
655     aij=0;
656     next=0;
657     head =0;
658     this->n=-1234567890;//  stupide number
659     this->m=-1234567801;//  stupide number
660     this->nnz=-1234567802;
661     // because the solver is some time deleted after ...
662 }
663 
664 template<class I,class R>
Sortij()665 void HashMatrix<I,R>::Sortij()
666 {
667 
668     if( state != sorted_ij)
669     {
670         CheckUnLock("Sortij");
671         HeapSort(i,j,aij,nnz);
672         ReHash();
673         nbsort++;
674         state = sorted_ij;
675     }
676 }
677 
678 template<class I,class R>
Sortji()679 void HashMatrix<I,R>::Sortji()
680 {
681     if( state != sorted_ji)
682     {
683         CheckUnLock("Sortji");
684         HeapSort(j,i,aij,nnz);
685         ReHash();
686         state = sorted_ji;
687         nbsort++;
688     }
689 }
690 
691 template<class I,class R>
set(I nn,I mm,int hhalf,size_t nnnz,I * ii,I * jj,R * aa,int f77,int tcsr)692 void HashMatrix<I,R>::set(I nn,I mm,int hhalf,size_t nnnz, I *ii, I*jj, R *aa,int f77,int tcsr)
693 {
694 //    tcsr >0 => CSR ii pointer size nn+1, jj col size nnnz
695 //    tcrs <0 = CSC  jj pointer size mm+1, ii row size nnnz
696 //    tcsr == 0 => COO
697 
698     clear();
699     this->n=nn;
700     this->m=mm;
701     this->N=nn;
702     this->M=mm;
703     fortran=f77;
704     half=hhalf;
705     Increaze(nnnz);
706     nnz=nnnz;
707     if(tcsr >=0) //  input CSR
708       HMcopy(j,jj,nnnz);// copy of col
709     if(tcsr <=0) //  input CSC
710       HMcopy(i,ii,nnnz);// copy of  row
711     if( tcsr>0)//  input CSR
712         for(I ip=0; ip< nn; ++ip)
713             for(I k=ii[ip]; k<ii[ip+1]; ++k)
714                 i[k-f77]=ip+f77;
715     else if( tcsr<0) //  input CSC
716         for(I jp=0; jp< mm; ++jp)
717             for(I k=jj[jp]; k<jj[jp+1]; ++k)
718                 j[k-f77]=jp+f77;
719 
720     HMcopy(aij,aa,nnnz);
721     ReHash();
722 }
723 
724 
725 template<class I,class R>
operator =(const HashMatrix & A)726 HashMatrix<I,R> & HashMatrix<I,R>::operator=(const HashMatrix &A)
727 {
728     if( (const void*)  this == (const void*) & A) return *this;
729     set(A.n,A.m,A.half,A.nnz,A.i,A.j,A.aij,A.fortran);
730     return *this;
731 }
732 
733 template<class I,class R>
operator +=(const HashMatrix & A)734 HashMatrix<I,R> & HashMatrix<I,R>::operator+=(const HashMatrix &A)
735 {
736     Add(&A);
737     return *this;
738 }
739 
740 template<class I,class R>
operator -=(const HashMatrix & A)741 HashMatrix<I,R> & HashMatrix<I,R>::operator-=(const HashMatrix &A)
742 {
743     Add(&A,R(-1.));
744     return *this;
745 }
746 
747 template<class I,class R>
Add(const HashMatrix<I,R> * PA,R coef,bool trans,I ii00,I jj00)748 void HashMatrix<I,R>::Add(const HashMatrix<I,R> *PA,R coef,bool trans, I ii00,I jj00)
749 {
750     const HashMatrix<I,R> &A=*PA;
751     I nn=A.n,mm=A.m;
752     I *ii=A.i;
753     I *jj=A.j;
754     size_t Annz= A.nnz;  ;
755     if( trans ) swap(nn,mm),swap(ii,jj);
756     nn = max(this->n,nn);
757     mm = max(this->m,mm);
758     resize(nn,mm);
759     if( (const void*) this == (const void*) & A && ii00==0 && jj00==0 && !trans )
760         for(I k=0; k < nnz; ++k)
761             aij[k] += coef*aij[k];
762     else
763     {   ffassert((const void*) this != (const void*) & A); // not the same matrix
764         if ( this->half!= A.half ) MATERROR(1,"+= on diff type of mat of type HashMatrix ");
765 
766         for(size_t k=0;k<Annz;++k)
767             operator()(ii[k]+ii00,jj[k]+jj00) += coef*A.aij[k];
768     }
769 }
770 
771 
772 template<class I,class R,class K>
Addto(HashMatrix<I,R> * P0,const HashMatrix<I,K> * PA,R (* f)(K),bool trans,I ii00,I jj00)773 void Addto(HashMatrix<I,R> *P0, const HashMatrix<I,K> *PA,R (*f)(K) ,bool trans, I ii00,I jj00)
774 {
775     const HashMatrix<I,K> &A=*PA;
776     I nn=A.n,mm=A.m;
777     I *ii=A.i;
778     I *jj=A.j;
779     size_t Annz= A.nnz;  ;
780     if( trans ) swap(nn,mm),swap(ii,jj);
781     nn = max(P0->n,nn);
782     mm = max(P0->m,mm);
783     P0->resize(nn,mm);
784 
785     {
786         if( (const void*) P0 == (const void*) & A && ii00==0 && jj00==0 && !trans )
787             for(I k=0; k < P0->nnz; ++k)
788                 P0->aij[k] += f(A.aij[k]);
789         else
790         {   ffassert((const void*) P0 != (const void*) & A); // not the same matrix
791             if ( P0->half!= A.half ) MATERROR(1,"+= on diff type of mat of type HashMatrix ");
792 
793             for(size_t k=0;k<Annz;++k)
794                 P0->operator()(ii[k]+ii00,jj[k]+jj00) += f(A.aij[k]);
795         }
796     }
797 }
798 
799 
800 
801 
802 template<class I,class R>
operator *=(R v)803 void HashMatrix<I,R>::operator*=(R v)
804 {
805     re_do_numerics=1;
806     for(int k=0; k < nnz; ++k)
807         aij[k] *= v;
808 
809 }
810 
811 template<class I,class R>
operator =(const R & v)812 void HashMatrix<I,R>::operator=(const R & v)
813 {
814     re_do_numerics=1;
815     for(int k=0; k < nnz; ++k)
816         aij[k] = v;
817 
818 }
819 
820 template<class I,class R>
HM()821 void HashMatrix<I,R>::HM()
822 {
823     re_do_numerics++;
824     re_do_symbolic++;
825     setp(0);
826     type_state=type_HM;
827 
828 }
829 
830 template<class I,class R>
COO()831 void HashMatrix<I,R>::COO()
832 {
833      if(verbosity>999) cout << " HashMatrix:: COO "<< this->n << " x " << this->m << " "<< p << " / "<< sizep << endl;
834     Sortij();
835     setp(0);
836     type_state=type_COO;
837 }
838 
839 
840 template<class I,class R>
COO(I * & IA,I * & IJ,R * & A)841 void HashMatrix<I,R>::COO(I *& IA, I *& IJ, R *& A)
842 {
843     Sortij();
844     IA=i;
845     IJ=j;
846     A=aij;
847     setp(0);
848     type_state=type_COO;
849 }
850 
851 
852 template<class I,class R>
CSR(I * & IA,I * & JA,R * & A)853 void HashMatrix<I,R>::CSR(I *& IA, I *& JA, R *& A)
854 {
855     Sortij();
856     Buildp(this->n,i,type_CSR);
857     IA=p;
858     JA=j;
859     A=aij;
860     type_state=type_CSR;
861 }
862 
863 template<class I,class R>
CSR()864 void HashMatrix<I,R>::CSR()
865 {
866      if(verbosity>999) cout << " HashMatrix:: csr()  " << state << " " << sorted_ij<< endl;
867     Sortij();
868     Buildp(this->n,i,type_CSR);
869     type_state=type_CSR;
870 }
871 
872 
873 template<class I,class R>
CSC()874 void HashMatrix<I,R>::CSC()
875 {
876     Sortji();
877     Buildp(this->m,j,type_CSC);
878     type_state=type_CSC;
879 }
880 
881 
882 template<class I,class R>
CSC(I * & JA,I * & IA,R * & A)883 void HashMatrix<I,R>::CSC(I *& JA, I *& IA, R *& A)
884 {
885 
886     Sortji();
887     Buildp(this->m,j,type_CSC);
888     JA=p;
889     IA=i;
890     A=aij;
891     type_state=type_CSC;
892 }
893 
894 
895 
896 template<class I,class R>
SortLU(int U)897 size_t HashMatrix<I,R>::SortLU(int U)
898 {   // put U or L at the top and return the nnz term associated
899     if( half) {
900 
901     }
902     size_t nnzu =0,nnzd=0;
903     for(int k=0; k<nnz; ++ k)
904         if( i[k]<j[k]) ++nnzu; // 1 2 in U => i <= j
905         else if (i[k]==j[k]) ++nnzd;
906     size_t nnzl=nnz-nnzu-nnzd;
907     cout << "SortLU " << U << ":s " << nnzl << " "<< nnzd << " " << nnzu << endl;
908     size_t kR=0;
909     if(U>0)
910     {
911         size_t kU=0L,k;
912 
913         for( k=0; k<nnz; ++ k)
914             if( i[k]<=j[k])
915             {
916                 std::swap(aij[k],aij[kU]);
917                 std::swap(i[k],i[kU]);
918                 std::swap(j[k],j[kU++]);
919             } // 1 2 in U => i <= j
920         // verif
921         cout << " SortLU "<< kU << " " << nnzu+nnzd << " " << k << endl;
922 
923         for( k=0; k<kU ; ++k)
924             assert(i[k]<=j[k]);
925         for( k=kU; k<nnz ; ++k)
926             assert(i[k]>j[k]);
927 
928         kR = kU;
929         assert(kU == nnzu+nnzd);
930     }
931     else
932     {
933         size_t kL=0L,k;
934 
935         for( k=0; k<nnz; ++ k)
936             if( i[k]>=j[k])
937             {
938                 std::swap(aij[k],aij[kL]);
939                 std::swap(i[k],i[kL]);
940                 std::swap(j[k],j[kL++]);
941             }
942 
943         kR = kL;
944 
945 
946     }
947     ReHash();
948     nbsort++;
949     state += addstateLU(U) ;
950     return kR;
951 
952 }
953 
954 template<class I,class R>
CSC_U(I * & JA,I * & IA,R * & A)955 size_t HashMatrix<I,R>::CSC_U(I *& JA, I *& IA, R *& A)
956 {
957 
958     if(type_CSC!=type_state)
959         HeapSort(j,i,aij,nnz);
960     state=type_CSC;
961     size_t nnzu=SortLU(1);
962     Buildp(this->m,j,type_CSC+addstateLU(1),nnzu);
963 
964     JA=p;
965     IA=i;
966     A=aij;
967     return nnzu;
968 }
969 
970 template<class I,class R>
CSR_L(I * & IA,I * & JA,R * & A)971 size_t HashMatrix<I,R>::CSR_L(I *& IA, I *& JA, R *& A)
972 {
973     if(type_CSR!=type_state)
974         HeapSort(i,j,aij,nnz);
975     state=type_CSR;
976     size_t nnzl=SortLU(-1);
977     Buildp(this->m,i,type_CSR+addstateLU(-1),nnzl);
978 
979     IA=p;
980     JA=j;
981     A=aij;
982     return nnzl;
983 }
984 
985 
986 template<class I,class R>
Buildp(I nn,I * IA,int type_m,size_t nnzz)987 void HashMatrix<I,R>::Buildp(I nn,I * IA,int type_m,size_t nnzz)
988 {
989    if(verbosity>999) cout << " HashMatrix:: Buildp"<< this->n<< " x " << this->m << " " << nn << " " << IA<< " " << type_m << " " << nnzz << " / " << nnz << " / p=" << p <<endl;
990     if(nnzz==0) nnzz=nnz;
991     if(type_m != type_state)
992     {
993 
994         assert( state==type_m);
995 
996         setp(nn+1);
997          //int shift =  fortran;
998         std::fill(p,p+nn+1,-1);
999         p[nn] = fortran+ (I) nnzz;
1000         for( I k=I(nnzz)-1; k>=0 ; --k )
1001         {
1002             p[IA[k]-fortran] = k+fortran;
1003             ffassert( (IA[k]-fortran>=0 ) && (IA[k]-fortran<=nn));
1004         }
1005         p[nn]=nnzz+fortran;
1006         // remove empty row
1007         for(I ii=nn-1;ii>=0;--ii)
1008             if(p[ii]<0)//  empty row
1009                 p[ii]=p[ii+1];
1010 
1011 
1012 
1013         type_state=type_m;
1014         if(verbosity>10)
1015         {
1016             cout << "Buildp  p=" << nn<< " type_m =" << type_m << endl;
1017             for(int ii=0;ii<=nn;++ii)
1018                 cout << ii << " " << p[ii] << endl;
1019             cout << " end Buildp\n";
1020         }
1021 
1022     }
1023     assert(p);
1024 }
1025 
1026 
1027 
1028 template<class I,class R>
RoworCol(I ii,bool row)1029 typename HashMatrix<I,R>::Pair HashMatrix<I,R>::RoworCol(I ii,bool row)
1030 {
1031     size_t k0=0,k1=nnz-1,k=nnz-1;
1032     I * pp=0;
1033     if(row)
1034     {
1035         if(state != sorted_ij) Sortij();
1036         pp= i;
1037     }
1038     else
1039     {
1040         if(state != sorted_ji) Sortji();
1041         pp = j;
1042     }
1043     assert(nbsort < (size_t) this->n+100 );
1044 
1045     while(pp[k0] < ii)
1046     {
1047         if( k-k0 <=1) break;
1048         size_t kk = (k0+k)/2;
1049         //cout << kk << " " << k0 << endl;
1050         if(pp[kk]<ii)
1051             k0=kk;
1052         else
1053         {
1054             if(pp[kk]>ii) k1=kk;
1055             k=kk;
1056         }
1057     }
1058     if(pp[k0]<ii)
1059         ++k0;
1060     k=k0;
1061     while(ii<pp[k1])
1062     {
1063         if( ( k1-k ) <= 1) break;
1064         size_t kk = (k1+k)/2;
1065         //cout << kk << " " << k1 << endl;
1066         if(pp[kk]>ii)
1067             k1=kk;
1068         else
1069             k=kk;
1070     }
1071     if(pp[k1]>ii)
1072         --k1;
1073     //cout << ii << "k0 =" << k0 << " k1 " << k1 << endl;
1074     Pair  ret(k0,k1);
1075     return ret;
1076 }
1077 
1078 
1079 
1080 template<class I,class R>
addMatMul(R * x,R * Ax,bool Transpose,I sx,I sAx) const1081 R* HashMatrix<I,R>::addMatMul(R *x,R*Ax,bool Transpose,I sx,I sAx) const {
1082     I *ii=i,*jj=j;
1083     R *aa=aij;
1084     double t0= CPUsecond();
1085     //   if(Transpose != trans) {std::swap(ii,jj);}
1086     if(Transpose ) {std::swap(ii,jj);}
1087     if(fortran) {aa++;}
1088     if(sx==1 && sAx==1)
1089     {
1090         if( half)
1091             for(int k=0; k<nnz;++k)
1092             {
1093                 Ax[ii[k]] += aa[k]*x[jj[k]];
1094                 if( ii[k] != jj[k]) Ax[jj[k]] += (half == 1?conj(aa[k]):aa[k]) *x[ii[k]];
1095             }
1096         else
1097             for(int k=0; k<nnz;++k)
1098                 Ax[ii[k]] += aa[k]*x[jj[k]];
1099     }
1100     else
1101     if( half)
1102         for(int k=0; k<nnz;++k)
1103         {
1104             Ax[ii[k]*sAx] += aa[k]*x[jj[k]*sx];
1105             if( ii[k] != jj[k]) Ax[jj[k]*sAx] += (half == 1?conj(aa[k]):aa[k])*x[ii[k]*sx];
1106         }
1107     else
1108         for(int k=0; k<nnz;++k)
1109             Ax[ii[k]*sAx] += aa[k]*x[jj[k]*sx];
1110     matmulcpu+=  CPUsecond()-t0;
1111     return Ax;}
1112 
1113 
1114 template<class I,class R>
trace() const1115 R HashMatrix<I,R>::trace () const {
1116     R t=R();
1117     for(int ii=0; ii<this->n; ++ii)
1118         t+= diag(ii);
1119     return t;
1120 }
1121 
1122 template<class I,class R>
FrobeniusNorm() const1123 double HashMatrix<I,R>::FrobeniusNorm() const
1124 {
1125     double s=0;
1126     for(size_t k=0; k<nnz;++k)
1127         s += norm(aij[k]);
1128     s = sqrt(s);
1129     return  s;
1130 }
1131 
1132 template<class I,class R>
norm1() const1133 double HashMatrix<I,R>::norm1() const
1134 {
1135     double  s=0;
1136     for(size_t k=0; k<nnz;++k)
1137         s += abs(aij[k]);
1138     return  s;
1139 }
1140 
1141 template<class I,class R>
norminfty() const1142 double HashMatrix<I,R>::norminfty() const
1143 {
1144     double s=abs(aij[0]);
1145     for(size_t k=0; k<nnz;++k)
1146         s = std::max(abs(aij[k]),s);
1147     return  s;
1148 }
1149 
1150 
1151 template<class I,class R>
SetBC(char * wbc,double ttgv)1152 void HashMatrix<I,R>::SetBC(char *wbc,double ttgv)
1153 {
1154     tgv = ttgv;
1155     ntgv =0;
1156     if(ttgv<0)
1157         CSR();
1158     if ( this->n != this->m) MATERROR(1,"SetBC on none square matrix  ?????");
1159     for(I ii=0; ii< this->n; ++ii)
1160         if( wbc[ii] )
1161         {
1162             ntgv++;
1163             if(ttgv<0)
1164             {
1165 
1166                 if( wbc[ii] )
1167                 {
1168                     for (I k=p[ii];k<p[ii+1]; ++k)
1169                         if( j[k]==ii )
1170                             aij[k] = (std::abs(ttgv+10.0) < 1.0e-10 ? 0.0 : 1.0);
1171                         else
1172                             aij[k]=0;// put the line to Zero.
1173                 }
1174 
1175             }
1176             else
1177                 operator()(ii,ii)=ttgv;
1178         }
1179     if(std::abs(ttgv+2.0) < 1.0e-10 || std::abs(ttgv+20.0) < 1.0e-10) //  remove also columm tgv == -2 .....
1180     {
1181         CSC();
1182         for(I jj=0; jj< this->n; ++jj)
1183             if( wbc[jj] ) {
1184                 for (I k=p[jj];k<p[jj+1]; ++k)
1185                     if( i[k]!=jj || std::abs(ttgv+20.0) < 1.0e-10)
1186                         aij[k]=0;//
1187             }
1188     }
1189 
1190 
1191 }
1192 
1193 
1194 
1195 template<class I,class R>
addMap(R coef,std::map<pair<I,I>,R> & mij,bool trans,I ii00,I jj00,bool cnj,double threshold)1196 void HashMatrix<I,R>::addMap(R coef,std::map< pair<I,I>, R> &mij,bool trans,I ii00,I jj00,bool cnj,double threshold)
1197 {
1198     for (auto pm=mij.begin();  pm != mij.end(); ++pm)
1199     {
1200         I ii = pm->first.first+ii00, jj= pm->first.second+jj00;
1201         R cmij = coef* pm->second;
1202 
1203         if(trans) swap(ii,jj);
1204         if( abs(cmij)  > threshold )
1205         {
1206           if(cnj) cmij = conj(cmij);
1207           operator()(ii,jj) += cmij;
1208         }
1209     }
1210 }
1211 
1212 template<class I,class R>
addMatTo(R coef,HashMatrix<I,R> & mij,bool trans,I ii00,I jj00,bool cnj,double threshold,const bool keepSym)1213 bool HashMatrix<I,R>::addMatTo(R coef,HashMatrix<I,R> & mij,bool trans,I ii00,I jj00,bool cnj,double threshold,const bool keepSym)
1214 {
1215     //  add a mij + = coef * [(this)^trans^cnj ,
1216     double eps0=max(numeric_limits<double>::min(),threshold);
1217 
1218     if (half)
1219     {
1220         for( size_t kk= 0; kk<nnz ;++kk)
1221         {
1222             I ii=i[kk], jj=j[kk];
1223             R cij =  coef* ( cnj ? RNM::conj(aij[kk]) : aij[kk]);
1224             if(threshold==0 || std::norm(cij)>eps0) // remove for IPOPT april 2019 FH only if threshold >0
1225             {
1226                 mij[ij_mat(trans,ii00,jj00,ii,jj)] += cij ;
1227                 if (ii!=jj&&!keepSym)
1228                     mij[ij_mat(trans,ii00,jj00,jj,ii)] += cij;
1229             }
1230         }
1231     }
1232     else
1233     {
1234         for(size_t  kk= 0; kk<nnz ;++kk)
1235         {
1236             I ii=i[kk], jj=j[kk];
1237             R cij =  coef* ( cnj ? RNM::conj(aij[kk]) : aij[kk]);
1238             if(threshold==0 || std::norm(cij)>eps0) // / remove for IPOPT april 2019 FH
1239                 mij[ij_mat(trans,ii00,jj00,ii,jj)] += cij;
1240         }
1241     }
1242 
1243     return keepSym;
1244 }
1245 
1246 
1247 template<class I,class R>
operator +=(MatriceElementaire<R> & me)1248 VirtualMatrix<I,R>  & HashMatrix<I,R>::operator +=(MatriceElementaire<R> & me) {
1249     //  R zero=R();
1250     int il,jl,i,j;
1251     int * mi=me.ni, *mj=me.nj;
1252     if ((this->n==0) && (this->m==0))
1253     {
1254 
1255 
1256         cout << "  -- Bug: HashMat  is empty let's build it" << endl;
1257         ffassert(0);
1258 
1259     }
1260     R * al = me.a;
1261     R * aij;
1262     switch (me.mtype) {
1263         case MatriceElementaire<R>::Full : ffassert(!half);
1264             for (il=0; il<me.n; ++il)  { i=mi[il];
1265                 for ( jl=0; jl< me.m ; ++jl,++al)  {j=mj[jl];
1266                     aij = npij(i,j);
1267                     {
1268                         throwassert(aij);
1269                         *aij += *al;}}}
1270             break;
1271 
1272         case MatriceElementaire<R>::Symmetric : ffassert(half);
1273             for (il=0; il<me.n; ++il) {  i=mi[il] ;
1274                 for (jl=0;jl< il+1 ; ++jl) { j=mj[jl];
1275                     aij =    (j<i) ? npij(i,j) : npij(j,i);
1276                     throwassert(aij);
1277                     *aij += *al++;}}
1278             break;
1279         default:
1280             cerr << "Big bug type MatriceElementaire unknown" << (int) me.mtype << endl;
1281             ErrorExec("Bi-ug in  Hashmat += MatriceElementaire",1);
1282             break;
1283     }
1284     return *this;
1285 }
1286 
1287 
1288 template<class I,class R>
SetBC(I ii,double ttgv)1289 void SetBC(I ii,double ttgv) { diag(ii)=ttgv;};
1290 
1291 
1292 
1293 template<class I,class R>
gettgv(I * pntgv,double ratio) const1294 double HashMatrix<I,R>::gettgv(I * pntgv,double ratio) const
1295 {
1296     if( this->n != this->m) return 0; // no ttgv
1297     double ttgv =0, max1=0;
1298     I ntgv=0;
1299     for (I ii=0; ii<this->n;++ii)
1300     {
1301         R * p=pij(ii,ii);
1302         if (p)
1303         {
1304 
1305             double a=real(*p);
1306             if( a> ttgv )
1307             {
1308                 max1=ttgv;
1309                 ttgv=a;
1310                 ntgv =1;
1311             }
1312             else if (a== ttgv)
1313                 ++ntgv;
1314             else if( a >max1)
1315                 max1 = a;
1316         }
1317 
1318     }
1319     if( max1*ratio> ttgv) // ttgv to small => no ttgv ....
1320         ttgv=0,ntgv=0; // no ttgv
1321     if(pntgv) *pntgv = ntgv;
1322     return ttgv;
1323 
1324 }
1325 
1326 template<class I,class R>
setsdp(int sym,bool dp)1327 void HashMatrix<I,R>::setsdp(int sym,bool dp) // sym dpos para
1328 {
1329     this->symetric=(sym > 0);
1330     this->positive_definite=dp;
1331     if( (half>0) != sym)
1332     {
1333         if( sym)
1334             Half(sym);//  remove half part
1335         else
1336             UnHalf();
1337     }
1338     else {
1339       half = sym;
1340     }
1341 }
1342 
1343 template<class I,class R>
UnHalf()1344 void HashMatrix<I,R>::UnHalf()
1345 {
1346 
1347     if (half==0) return;
1348     HM();
1349     size_t nnz0=nnz,err=0;
1350     for(int k=0; k<nnz0; ++k)
1351       if( i[k] > j[k] )
1352        {
1353         size_t ki=insert(j[k],i[k],(half == 1?conj(aij[k]):aij[k]));
1354         err += ki < nnz0;
1355         }
1356     else if ( i[k] < j[k] )
1357         err++;
1358     half = 0;
1359     if( err )
1360         cerr << " Try of unsymmetrize no half matrix (Bug) ... ????" <<endl;
1361     ffassert(err==0);
1362 
1363 }
1364 
1365 typedef double R;
1366 typedef complex<R> C;
1367 //  because for UMFPACK 64  because long are only 32 bits under windows
1368 #ifdef _WIN32
1369 typedef long long  int64;
1370 #else
1371 typedef long int64;
1372 #endif
1373 
1374 template class HashMatrix<int,R>;
1375 template class HashMatrix<int,C >;
1376 template class HashMatrix<int64,R>;
1377 template class HashMatrix<int64,C >;
1378 
1379 template HashMatrix<int64,R>::HashMatrix(const HashMatrix<int,R> & );
1380 template HashMatrix<int64,C>::HashMatrix(const HashMatrix<int,C> & );
1381 template HashMatrix<int,C>::HashMatrix(const HashMatrix<int64,C> & );
1382 template HashMatrix<int,R>::HashMatrix(const HashMatrix<int,C> & , R(*ff)(C));
1383 //template HashMatrix<int,C>::HashMatrix(const HashMatrix<int,R> & );
1384 
1385 //template HashMatrix<int,R> & HashMatrix<int,R>::operator=(const HashMatrix<int,C> & );
1386 //template HashMatrix<int,C> & HashMatrix<int,C>::operator=(const HashMatrix<int,R> & );
1387 //template HashMatrix<int,R> & HashMatrix<int,R>::operator+=(const HashMatrix<int,C> & );
1388 //template HashMatrix<int,C> & HashMatrix<int,C>::operator+=(const HashMatrix<int,R> & );
1389 
1390 template  void Addto<int,R,C>(HashMatrix<int,R> *P0, const HashMatrix<int,C> *PA,R (*f)(C) ,bool trans, int ii00,int jj00);
1391 template  void Addto<int,C,R>(HashMatrix<int,C> *P0, const HashMatrix<int,R> *PA,C (*f)(R) ,bool trans, int ii00,int jj00);
1392 
1393 template void HashMatrix<int,R>::set<int64>(int64 nn,int64 mm,int hhalf,size_t nnnz, int64 *ii, int64*jj, R *aa,int f77);
1394 //template void HashMatrix<int,R>::set<int,R>(int nn,int mm,bool hhalf,size_t nnnz, int *ii, int*jj, R *aa,int f77,R(*ff)(R));
1395 //template void HashMatrix<int,R>::set<long,R>(long nn,long mm,bool hhalf,size_t nnnz, long *ii, long *jj, R *aa,int f77,R(*ff)(R));
1396 template void HashMatrix<int,C>::set<int,C>(int nn,int mm,int hhalf,size_t nnnz, int *ii, int*jj, C *aa,int f77,C(*ff)(C));
1397 template void HashMatrix<int,R>::set<int,C>(int nn,int mm,int hhalf,size_t nnnz, int *ii, int*jj, C *aa,int f77,R(*ff)(C));
1398 
1399 //template void HashMatrix<int,C>::set<int,R>(int nn,int mm,bool hhalf,size_t nnnz, int *ii, int*jj, R *aa,int f77,C(*ff)(R));
1400 //template void HashMatrix<int,C>::set<long,C>(long nn,long mm,bool hhalf,size_t nnnz, long *ii, long *jj, C *aa,int f77,C(*ff)(C));
1401 
1402 
init_HashMatrix()1403  void init_HashMatrix ()
1404 {
1405 
1406 }
1407 // just to test
1408 /*
1409 static void tttt() {
1410     HashMatrix<int,R> AiR(10);
1411     HashMatrix<long,R> AlR(AiR);//,HashMatrix<long,R>::cast_funct);
1412 
1413     HashMatrix<long,C> AlC(10);
1414     HashMatrix<int,C> AiC(10);
1415 
1416 }
1417 
1418 */
1419