1 #ifndef VIRTUALSOLVERSKYLINE_HPP_
2 #define VIRTUALSOLVERSKYLINE_HPP_
3
4 #include <iostream>
5 #include <cmath>
6 #include "HashMatrix.hpp"
7 #include <vector>
8 #include "VirtualSolver.hpp"
9
10 #include <complex>
11 #include "SkyLineSolver.hpp"
12 #include <queue>
13 #include <vector>
14 #include <utility>
15 //#include <climits>
16
17 template <class Z>
order_CutHill_McKee(Z n,Z * Ap,Z * Ai,Z * p)18 Z order_CutHill_McKee(Z n, Z *Ap, Z* Ai,Z*p)
19 {// CutHill-McKee
20 std::vector<Z> mark(n,-1),r(n);
21 Z nr =0;
22 unsigned long long pft=0,pft0=0;
23 std::multimap<Z,Z> q,qn;
24 bool debug =0;
25 unsigned long long pfso=0 ,pfs=0,pfss=0,pfs0=0;
26 for(Z id=0; id<n; ++id)
27 if(mark[id]==-1)
28 { // start point
29 Z nri=nr;
30 Z s=id;
31 Z si=s;
32 for(int step=0; step < 5; ++step)
33 {
34
35 for (int i= nri; i< nr;++i)
36 mark[r[i]]=-1; // clean
37 nr=nri; //
38 s=si; // last point ..
39 mark[s]=nr;
40 r[nr++] = s;
41 Z sc = Ap[s+1]-Ap[s];
42 q.clear();
43 q.insert(std::make_pair(-sc,s));
44 if(debug) cout << " s= "<< s << " "<< nr << endl;
45 int lvl=0;
46 while (q.size()>0)
47 {
48 if(debug) cout << lvl++ <<" " << nr << " : ";
49
50 qn.clear();
51 for(auto iq=q.begin(); iq != q.end(); ++iq )
52 {
53 int s = iq->second;
54 for(int k= Ap[s]; k<Ap[s+1];++k)
55 {
56 Z ss= Ai[k];
57 if(mark[ss]<0)
58 {
59 Z ssc = Ap[ss+1]-Ap[ss];
60 qn.insert(std::make_pair(-ssc,ss));
61 mark[ss]=nr;
62 r[nr++] = ss;
63 if(debug) cout << ss <<" ";
64 }
65
66 }
67 }
68 if(debug) cout << endl;
69 swap(q,qn);
70 }
71
72 // compute SkyLine
73
74 pfs =0;
75 pfs0= 0;
76 for(int i= nri; i<nr; ++i)
77 {
78 Z s = r[i];
79 Z imn = s,i0mn=i;
80 Z imx = s,i0mx=i;
81 for(int k= Ap[s]; k<Ap[s+1];++k)
82 {
83 Z ss0=Ai[k];
84 Z ss= mark[ss0];
85 imn = min(imn,ss);
86 imx = max(imx,ss);
87 i0mn = min(i0mn,ss0);
88 i0mx = max(i0mx,ss0);
89 }
90 pfs += s-imn;
91 pfs0 += i-i0mn;
92 }
93 if(debug) cout << " ** " << step << " " << pfso << " " << pfs << " " << pfs0 << endl;
94 if(!pfso) pfso= pfs0;
95 if(pfs <pfso )
96 {
97 for(int i= nri; i<nr; ++i)
98 p[r[i]]=mark[r[i]];
99 if(debug)
100 cout << step<< " pfs " << pfs << " " << pfso << endl;
101 pfss=pfs;
102 }
103 else if(step==0)
104 {
105 if(debug) cout << step<< " ID pfs " << pfs << " " << pfs0 << endl;
106 for(int i= nri; i<nr; ++i)
107 p[r[i]]=r[i];
108 }
109 if((pfso <pfs) && (step > 5) ) break;
110 pfso=min(pfs,pfso);
111 }
112
113
114 pft0 += pfs0;
115 pft += pfss;
116 if(verbosity>1) cout << "order_CutHill_McKee: pft = " << pft << " pft id = " << pft0 << endl;
117 }
118 if(verbosity>9)
119 {
120 cout << " " << n << endl;
121 for(int i=0; i< n; ++i)
122 { if(i%10==0) cout << endl;
123 cout << p[i] << " ";
124 }
125 cout << endl<< endl;
126 }
127 return (Z) pft;
128 }
129
130 template<class Z=int,class K=double>
131 class VirtualSolverSkyLine: public VirtualSolver<Z,K> {
132 public:
133 // 1 unsym , 2 sym, 4 pos , 8 nopos, 16 seq, 32 ompi, 64 mpi ,
134 static const int orTypeSol = 1|2|4|8|16;
135
136
137 typedef HashMatrix<Z,K> HMat;
138 typedef SkyLineMatrix<Z,K> SLMat;
139 int typesolver;
140 HMat *A;
141 SLMat *SL;
142 int cs,cn;
143 Z *p; // permuation old -> new
144 K *v;
145 double tol_pivot;
146 int verb;
147 mutable int status;
ttypesolver(const string & name)148 static int ttypesolver(const string &name) {
149 if(name.length() <=2) return 3;
150 else if( name[0] =='C') return 1 + ( (name[1]) != 'H');
151 else return 3; }
VirtualSolverSkyLine(HMat & AA,const Data_Sparse_Solver & ds,Stack stack)152 VirtualSolverSkyLine(HMat &AA, const Data_Sparse_Solver & ds,Stack stack )
153 :typesolver(ttypesolver(ds.solver)),A(&AA),SL(0),cs(0),cn(0),p(0),v(0),
154 tol_pivot(ds.tol_pivot<0 ? 1e-15 :ds.tol_pivot) , verb(ds.verb) {
155 if( verb>3) cout << " -- SkyLineMatrix: "<<typesolver << " " <<ds.solver<< endl;
156 }
157
158
UpdateState()159 void UpdateState(){
160 if( A->GetReDoNumerics() ) cn++;
161 if( A->GetReDoSymbolic() ) cs++;
162 this->ChangeCodeState(A->n,cs,cn);
163 }
dosolver(K * x,K * b,int N,int trans)164 void dosolver(K *x,K*b,int N,int trans) {
165 if(verb>2|| verbosity>9)
166 std::cout <<" dosolver::SkyLine" << N<< " "<< trans << " WF " << SL->whichfac
167 << " size " << SL->size() << " typesolver: " <<typesolver
168 << std::endl;
169
170 for(int k=0,oo=0; k<N;++k, oo+= A->n)
171 {
172 K * xx = x+oo, *bb = b+oo;
173 for(Z i=0; i<A->n; ++i)
174 v[i] = bb[p[i]];
175 SL->solve(v,trans);
176 for(Z i=0; i<A->n; ++i)
177 xx[p[i]] = v[i];
178
179 }
180 }
fac_init()181 void fac_init(){
182
183 if(p) delete p;
184 p = new Z[A->n];
185 if(v) delete v;
186 v = new K[A->n];
187 }
fac_symbolic()188 void fac_symbolic(){
189 // build permuation
190
191 if( A->n !=A->m) MATERROR(1,"VirtualSolverSkyLine non squre mat");
192 Z *Ap,*Ai;
193 K * Ax;
194
195 A->CSC(Ap,Ai,Ax);
196
197 Z pfl =order_CutHill_McKee(A->n,Ap,Ai,p);
198 // rebuild perumation
199 if(verb>2|| verbosity>9)
200 std::cout <<" Skyline::fac_symbolic : skyline size " << pfl << std::endl;
201
202
203 }
fac_numeric()204 void fac_numeric(){
205 if(SL) delete SL;
206 SL = new SLMat(A,p,typesolver,verb);
207 if(verbosity>1) std::cout << " size of Skyline mat ="<< SL->size() << " nz :" <<SL->pL[A->n] << " nzz " << A->nnz << " typesolver " << typesolver <<endl;
208 if(verbosity>99)
209 {
210 cout << *SL << endl;
211 cout << *A << endl;
212 }
213 SL->factorize(tol_pivot);
214
215 }
~VirtualSolverSkyLine()216 ~VirtualSolverSkyLine(){
217 if(p) delete [] p;
218 if(v) delete [] v;
219 if(SL) delete SL;
220 }
221
222 };
223
224
225
226 #endif
227