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