NeoPZ
pzskylmatpar.cpp
Go to the documentation of this file.
1 
6 #include <math.h>
7 #include <stdlib.h>
8 
9 #include <stdio.h>
10 #include <fstream>
11 using namespace std;
12 
13 
14 #ifdef BLAS
15 extern "C"{
16 #include <g2c.h>
17 #include "fblaswr.h"
18 };
19 
20 #endif
21 
22 #include "pzfmatrix.h"
23 #include "pzskylmat.h"
24 #include "pzskylmatpar.h"
25 
26 #include "pzlog.h"
27 
28 #include "pz_pthread.h"
29 
30 #ifdef LOG4CXX
31 static LoggerPtr logger(Logger::getLogger("pz.matrix.tpzskylparmatrix"));
32 #endif
33 
35 pthread_mutex_t skymutex = PTHREAD_MUTEX_INITIALIZER;
37 pthread_cond_t condition = PTHREAD_COND_INITIALIZER;
38 
40 //const int templatedepth = 10;
41 
42 //Constructors
43 template<class TVar>
44 TPZSkylParMatrix<TVar>::TPZSkylParMatrix(const int64_t dim, const TPZVec<int64_t> &skyline,int NumThreads)
45 
46 : TPZRegisterClassId(&TPZSkylParMatrix::ClassId),
47 TPZSkylMatrix<TVar>(dim, skyline),fDec(dim), fCorrectSingular(false)
48 {
49  int i;
50 
51  fSkyline = skyline;
52 
53  for(i=0;i<dim;i++) {
54  fDec[i]=skyline[i]-1;
55 
56  }
57  fNthreads = NumThreads;
58  fEqDec = -1;
59 }
60 
61 template<class TVar>
62 TPZSkylParMatrix<TVar>::TPZSkylParMatrix(const TPZSkylParMatrix<TVar> &copy) : TPZRegisterClassId(&TPZSkylParMatrix::ClassId),
63 TPZSkylMatrix<TVar>(copy), fDec(copy.fDec),fSkyline(copy.fSkyline),fEqDec(copy.fEqDec),
64 fNthreads(copy.fNthreads),fThreadUsed(0),fCorrectSingular(false)
65 {
66 }
67 
68 template<class TVar>
69 TPZSkylParMatrix<TVar>::TPZSkylParMatrix()
70 : TPZRegisterClassId(&TPZSkylParMatrix::ClassId),
71 TPZSkylMatrix<TVar>(),fDec(0),fCorrectSingular(false)
72 {
73  //something
74 }
75 
76 // TPZSkylParMatrix::TPZSkylParMatrix(const TPZSkylParMatrix &cp):TPZSkylMatrix(cp){
77 // fDec = cp.fDec;
78 // fSkyline = cp.fSkyline;
79 // fNthreads = cp.fNthreads;
80 // fEqDec = cp.fEqDec;
81 // }
82 
83 template<class TVar>
84 TPZSkylParMatrix<TVar>::~TPZSkylParMatrix() {
85 }
86 
87 template<class TVar>
88 void TPZSkylParMatrix<TVar>::SetSkyline(const TPZVec<int64_t> &skyline)
89 {
90  fSkyline = skyline;
91  int64_t i,neq=this->Dim();
92  for(i=0;i<neq;i++) fDec[i]=skyline[i]-1;
93  fEqDec=-1;
94  fCorrectSingular = false;
95  TPZSkylMatrix<TVar>::SetSkyline(skyline);
96 }
97 
98 template<class TVar>
99 void TPZSkylParMatrix<TVar>::PrintState() {
100  int i;
101  cout << "Decomposed equations " << fEqDec << ' ' << " Thread 'i' working on equation" << ' ';
102  for(i=0; i<fNthreads; i++) cout << fThreadUsed[i] << ' ';
103  cout << endl;
104  cout << "Columns to work " << ' ';
105  int64_t neq = this->Dim();
106  for(i=0; i<neq; i++) cout << fDec[i] << ' ';
107  cout << endl;
108  cout << "Singular equations ";
109  std::list<int64_t>::iterator it;
110  for (it=fSingular.begin(); it!=fSingular.end(); it++) {
111  cout << *it << ' ';
112  }
113  cout.flush();
114 }
115 
116 template<class TVar>
117 void TPZSkylParMatrix<TVar>::ColumnToWork(int64_t &lcol, int64_t &lprevcol) {
118  int64_t neq = this->Dim();
119  for(lcol=fEqDec+1;lcol<neq;lcol++)
120  {
121  int i;
122  for(i=0;i<fNthreads;i++){
123  if(fThreadUsed[i]==lcol){
124  //lcol=-1;
125  //lprevcol=-1;
126  break;
127  }
128  }
129  if(i < fNthreads && fThreadUsed[i] == lcol) continue;
130  if(lcol == fEqDec+1 && fDec[lcol] == fEqDec){
131  lprevcol=lcol;
132  return;
133  }else if(fDec[lcol]<fEqDec){
134  lprevcol = fDec[lcol] + 1;
135  return;
136  }
137  }
138  if(lcol == neq){
139  lcol = -1;
140  }
141 }
142 
143 template<class TVar>
144 void TPZSkylParMatrix<TVar>::ColumnToWork(int64_t &lcol) {
145  int64_t neq = this->Dim();
146  int64_t lcolentry = lcol;
147 #ifdef PZDEBUG
148  if(lcolentry < 0 || lcolentry >= neq)
149  {
150  LOGPZ_ERROR(logger,"ColumnToWork called with wrong lcol argument")
151  DebugStop();
152  }
153 #endif
154  lcol++;
155  lcol = lcol%neq;
156  while(lcol != lcolentry)
157  {
158  if(fColUsed.find(lcol) != fColUsed.end() || fDec[lcol] == lcol)
159  {
160  lcol++;
161  lcol = lcol%neq;
162  continue;
163  }
164  int decomposeduntil = fDec[lcol]+1;
165  if(fDec[decomposeduntil] == decomposeduntil)
166  {
167 #ifdef LOG4CXX
168  if (logger->isDebugEnabled()) {
169  std::stringstream sout;
170  sout << "Will work column " << lcol;
171  LOGPZ_DEBUG(logger,sout.str())
172  }
173 #endif
174  return;
175  }
176  lcol++;
177  if(lcol == neq) lcol = 0;
178  }
179  if(lcol == lcolentry){
180  lcol = -1;
181  }
182 }
183 
184 
185 template<class TVar>
186 void * TPZSkylParMatrix<TVar>::ParallelCholesky(void *t) {
187 
188  TPZSkylParMatrix<TVar> *loc = (TPZSkylParMatrix<TVar> *) t;
189  int64_t aux_col = -1;//, k;
190  int64_t col, prevcol;
191  prevcol=0;
192  PZ_PTHREAD_MUTEX_LOCK(&skymutex, "TPZSkylParMatrix::ParallelCholesky()");
193  int64_t neq = loc->Dim();
194  while (loc->fEqDec < neq-1) {
195  loc->ColumnToWork(col, prevcol);
196  if(col==-1) {
197  cout.flush();
198  if(neq-loc->fEqDec > loc->fNthreads){
199  PZ_PTHREAD_COND_WAIT(&condition, &skymutex, "TPZSkylParMatrix::ParallelCholesky()");
200  }else
201  {
202  //loc->fNthreads--;
203  //cout << "Terminating thread " << pthread_self() << endl;
204  break;
205  }
206  }
207  else {
208 
209  //Registers the working column number
210  int i;
211  for(i=0;i<loc->fNthreads;i++){
212  if(loc->fThreadUsed[i]==-1) {
213  loc->fThreadUsed[i]=col;
214  aux_col = i;
215  break;
216  }
217  }
218  PZ_PTHREAD_COND_SIGNAL(&condition, "TPZSkylParMatrix::ParallelCholesky()");
219  PZ_PTHREAD_MUTEX_UNLOCK(&skymutex, "TPZSkylParMatrix::ParallelCholesky()");
220  if(loc->fCorrectSingular)
221  {
222  loc->DecomposeColumn(col, prevcol,loc->fSingular);//loc->DecomposeColumn2(col, prevcol);
223  }
224  else {
225  loc->DecomposeColumn(col, prevcol);
226  }
227 
228  PZ_PTHREAD_MUTEX_LOCK(&skymutex, "TPZSkylParMatrix::ParallelCholesky()");
229  loc->fDec[col]=prevcol;
230  loc->fThreadUsed[aux_col]=-1;
231  if (col==prevcol) {
232  if(loc->fEqDec == col-1) {
233  loc->fEqDec = col;
234  if(!(loc->fEqDec%100) && loc->Dim() > 100) {
235  cout << loc->fEqDec << ' ';
236  cout.flush();
237  }
238  if(!(loc->fEqDec%1000)) cout << endl;
239  cout.flush();
240  }
241  else cout << "BIG MISTAKE\n";
242  }
243 
244  }
245  }
246  PZ_PTHREAD_MUTEX_UNLOCK(&skymutex, "TPZSkylParMatrix::ParallelCholesky()");
247  PZ_PTHREAD_COND_BROADCAST(&condition, "TPZSkylParMatrix::ParallelCholesky()");
248  cout << endl;
249  cout.flush();
250  return 0;
251 
252 }
253 
254 template<class TVar>
255 void * TPZSkylParMatrix<TVar>::ParallelLDLt(void *t) {
256 
257  TPZSkylParMatrix<TVar> *loc = (TPZSkylParMatrix<TVar> *) t;
258  int64_t aux_col = -1;//, k;
259  int64_t col, prevcol;
260  prevcol=0;
261  PZ_PTHREAD_MUTEX_LOCK(&skymutex,"TPZSkylParMatrix::ParallelLDLt()");
262  int64_t neq = loc->Dim();
263  while (loc->fEqDec < neq-1) {
264  loc->ColumnToWork(col, prevcol);
265  if(col==-1) {
266  cout.flush();
267  if(neq-loc->fEqDec > loc->fNthreads){
268  PZ_PTHREAD_COND_WAIT(&condition, &skymutex,"TPZSkylParMatrix::ParallelLDLt()");
269  }else
270  {
271  //loc->fNthreads--;
272  //cout << "Terminating thread " << pthread_self() << endl;
273  break;
274  }
275  }
276  else {
277 
278  //Registers the working column number
279  int i;
280  for(i=0;i<loc->fNthreads;i++){
281  if(loc->fThreadUsed[i]==-1) {
282  loc->fThreadUsed[i]=col;
283  aux_col = i;
284  break;
285  }
286  }
287  PZ_PTHREAD_COND_SIGNAL(&condition,"TPZSkylParMatrix::ParallelLDLt()");
288  PZ_PTHREAD_MUTEX_UNLOCK(&skymutex,"TPZSkylParMatrix::ParallelLDLt()");
289  loc->DecomposeColumnLDLt(col, prevcol);//loc->DecomposeColumn2(col, prevcol);
290  PZ_PTHREAD_MUTEX_LOCK(&skymutex,"TPZSkylParMatrix::ParallelLDLt()");
291  loc->fDec[col]=prevcol;
292  loc->fThreadUsed[aux_col]=-1;
293  if (col==prevcol) {
294  if(loc->fEqDec == col-1) {
295  loc->fEqDec = col;
296  if(!(loc->fEqDec%100) && loc->Dim() > 100) {
297  cout << loc->fEqDec << ' ';
298  cout.flush();
299  }
300  if(!(loc->fEqDec%1000)) cout << endl;
301  cout.flush();
302  }
303  else cout << "BIG MISTAKE\n";
304  }
305 
306  }
307  }
308  PZ_PTHREAD_MUTEX_UNLOCK(&skymutex,"TPZSkylParMatrix::ParallelLDLt()");
309  PZ_PTHREAD_COND_BROADCAST(&condition,"TPZSkylParMatrix::ParallelLDLt()");
310  cout << endl;
311  cout.flush();
312  return 0;
313 
314 }
315 
316 template<class TVar>
317 void * TPZSkylParMatrix<TVar>::ParallelLDLt2(void *t) {
318 
319  TPZSkylParMatrix<TVar> *loc = (TPZSkylParMatrix<TVar> *) t;
320  int64_t col = 0, prevcol;
321  prevcol=0;
322  PZ_PTHREAD_MUTEX_LOCK(&skymutex,"TPZSkylParMatrix::ParallelLDLt2()");
323  int64_t neq = loc->Dim();
324  while (loc->fNumDecomposed < neq) {
325  loc->ColumnToWork(col);
326  if(col==-1) {
327  cout.flush();
328  if(neq-loc->fNumDecomposed > loc->fNthreads){
329  PZ_PTHREAD_COND_WAIT(&condition, &skymutex,"TPZSkylParMatrix::ParallelLDLt2()");
330  col = 0;
331  }else
332  {
333  //loc->fNthreads--;
334  break;
335  }
336  }
337  else {
338 
339  //Registers the working column number
340  loc->fColUsed.insert(col);
341  PZ_PTHREAD_COND_SIGNAL(&condition,"TPZSkylParMatrix::ParallelLDLt2()");
342  PZ_PTHREAD_MUTEX_UNLOCK(&skymutex,"TPZSkylParMatrix::ParallelLDLt2()");
343  loc->DecomposeColumnLDLt2(col);//loc->DecomposeColumn2(col, prevcol);
344 
345  PZ_PTHREAD_MUTEX_LOCK(&skymutex,"TPZSkylParMatrix::ParallelLDLt2()");
346  if(loc->fDec[col] == col)
347  {
348  loc->fNumDecomposed++;
349  }
350  loc->fColUsed.erase(col);
351  }
352  }
353  loc->fNthreads--;
354  PZ_PTHREAD_COND_BROADCAST(&condition,"TPZSkylParMatrix::ParallelLDLt2()");
355 #ifdef VC
356  cout << "Terminating thread " << pthread_self().x << " numthreads = " << loc->fNthreads << " numdecomposed " << loc->fNumDecomposed << endl;
357 #else
358  cout << "Terminating thread " << pthread_self() << " numthreads = " << loc->fNthreads << " numdecomposed " << loc->fNumDecomposed << endl;
359 #endif
360  PZ_PTHREAD_MUTEX_UNLOCK(&skymutex,"TPZSkylParMatrix::ParallelLDLt2()");
361  cout << endl;
362  cout.flush();
363  return 0;
364 
365 }
366 
367 template<class TVar>
368 int TPZSkylParMatrix<TVar>::Decompose_Cholesky(std::list<int64_t> &singular)
369 {
370  bool sing = this->fCorrectSingular;
371  this->fCorrectSingular = true;
372  int ans = Decompose_Cholesky();
373  singular = this->fSingular;
374  this->fCorrectSingular = sing;
375  return ans;
376 }
377 
378 template<class TVar>
379 int TPZSkylParMatrix<TVar>::Decompose_Cholesky()
380 {
381 
382  if (this->fDecomposed == ECholesky) return 1;
383 
384  int64_t neq = this->Dim();
385 
386  cout << endl << "TPZSkylParMatrix::Decompose_Cholesky() -> Number of Equations = " << neq << endl;
387 #ifdef BLAS
388  cout << "Using BLAS " << endl;
389 #endif
390  cout.flush();
391 
392  pthread_t *allthreads = new pthread_t[fNthreads-1];
393  int *res = new int[fNthreads];
394 
395  fThreadUsed = new int[fNthreads];
396  int i;
397 
398 
399  //for(i=0; i<neq;i++) fDec[i]=-1;
400  //SetSkyline(&skyline);
401 
402  for(i=0; i<fNthreads;i++) fThreadUsed[i]=-1;
403 
404  fEqDec=-1;
405 
406  for(i=0;i<fNthreads-1;i++) {
407  res[i] = PZ_PTHREAD_CREATE(&allthreads[i], NULL,
408  this->ParallelCholesky, this, __FUNCTION__);
409  }
410 
411  ParallelCholesky(this);
412  for(i=0;i<fNthreads-1;i++) {
413  PZ_PTHREAD_JOIN(allthreads[i], NULL, __FUNCTION__);
414  }
415 
416  delete []allthreads;// fThreadUsed, fDec;
417 #ifdef LOG4CXX
418  if(fSingular.size())
419  {
420  std::stringstream sout;
421  sout << "Singular equations ";
422  std::list<int64_t>::iterator it;
423  for (it=fSingular.begin(); it!=fSingular.end(); it++) {
424  sout << *it << " ";
425  }
426  LOGPZ_WARN(logger,sout.str())
427  }
428 #endif
429  for(i=0;i<this->Dim();i++) {
430  fDec[i]=fSkyline[i]-1;
431  }
432  this->fDecomposed = ECholesky;
433 
434  fEqDec = -1;
435 
436  this->fDecomposed = ECholesky;
437 
438  return 1;
439 }
440 
441 template<class TVar>
442 int TPZSkylParMatrix<TVar>::Decompose_LDLt(std::list<int64_t> &singular)
443 {
444  bool sing = this->fCorrectSingular;
445  this->fCorrectSingular = true;
446  int ans = Decompose_LDLt();
447  singular = this->fSingular;
448  this->fCorrectSingular = sing;
449  return ans;
450 }
451 
452 template<class TVar>
453 int TPZSkylParMatrix<TVar>::Decompose_LDLt()
454 {
455  // return TPZSkylMatrix::Decompose_LDLt();
456  if (this->fDecomposed == ELDLt) return 1;
457 
458  int64_t neq = this->Dim();
459 
460  cout << endl << "TPZSkylParMatrix::Decompose_LDLt() -> Number of Equations = " << neq << endl;
461 #ifdef BLAS
462  cout << "Using BLAS " << endl;
463 #endif
464  cout.flush();
465  int nthreads = fNthreads;
466  // fNthreads=1;
467 
468  pthread_t *allthreads = new pthread_t[fNthreads-1];
469  int *res = new int[fNthreads];
470 
471  fThreadUsed = new int[fNthreads];
472  int i;
473 
474 
475  //for(i=0; i<neq;i++) fDec[i]=-1;
476  //SetSkyline(&skyline);
477 
478  for(i=0; i<fNthreads;i++) fThreadUsed[i]=-1;
479 
480  fNumDecomposed = 0;
481  for(i=0;i<this->Dim();i++) {
482  fDec[i]=fSkyline[i]-1;
483  if(fDec[i] == i-1)
484  {
485  fNumDecomposed++;
486  fDec[i] = i;
487  }
488  }
489  fEqDec=-1;
490 
491  for(i=0;i<nthreads-1;i++) {
492  res[i] = PZ_PTHREAD_CREATE(&allthreads[i], NULL, this->ParallelLDLt2,
493  this, __FUNCTION__);
494  }
495  ParallelLDLt2(this);
496  for(i=0;i<nthreads-1;i++) {
497  PZ_PTHREAD_JOIN(allthreads[i], NULL, __FUNCTION__);
498  }
499 
500  fNthreads = nthreads;
501  delete []allthreads;// fThreadUsed, fDec;
502  delete []fThreadUsed;
503  fThreadUsed = 0;
504 
505  for(i=0;i<this->Dim();i++) {
506  fDec[i]=fSkyline[i];
507  }
508  this->fDecomposed = ELDLt;
509 
510  fEqDec = -1;
511 
512  this->fDecomposed = ELDLt;
513 
514  return 1;
515 }
516 
517 template<>
518 void TPZSkylParMatrix<std::complex<double> >::DecomposeColumnCholesky(int64_t col, int64_t prevcol)
519 {
520  DebugStop();
521 }
522 
523 template<class TVar>
524 void TPZSkylParMatrix<TVar>::DecomposeColumnCholesky(int64_t col, int64_t prevcol){
525 
526  //Cholesky Decomposition
527  TVar *ptrprev; //Pointer to prev column
528  TVar *ptrcol; //Pointer to col column
529  int64_t skprev, skcol; //prev and col Skyline height respectively
530  int64_t minline;
531 
532  skprev = this->SkyHeight(prevcol);
533  skcol = this->SkyHeight(col);
534 
535  ptrprev = this->Diag(prevcol);
536  ptrcol = this->Diag(col);
537 
538  if((prevcol-skprev) > (col-skcol)){
539  minline = prevcol - skprev;
540  }else
541  {
542  minline = col - skcol;
543  }
544  if(minline > prevcol) {
545  cout << "error condition col " << col << " prevcol " << prevcol << "\n";
546  cout.flush();
547  return;
548  }
549  TVar *run1 = ptrprev + (prevcol-minline);
550  TVar *run2 = ptrcol + (col-minline);
551  TVar sum = 0;
552 
553  // int64_t n1=1;
554  // int64_t n2=1;
555 #ifndef BLAS
556 #ifdef USETEMPLATE
557  while(run1-ptrprev > templatedepth) {
558  run1-=templatedepth;
559  run2-=templatedepth;
560  sum += TemplateSum<templatedepth>(run1--,run2--);
561  }
562 #endif
563 
564 
565  while(run1 != ptrprev) sum += (*run1--)*(*run2--);
566 #else
567  //Blas dot product implementation
568  //cout << "Calling BLAS " << endl;
569  int64_t n = run1 - ptrprev;
570  if (n>0){
571  run1-=n-1;
572  run2-=n-1;
573  sum = ddot_(&n, run1, &n1, run2, &n2);
574 
575  /*for(int i=0;i<n;i++) {
576  cout << "run1 "<< run1[i] << " run2 " << run2[i] << "\n";
577  }
578  cout << " sum " << sum << " n " << n <<"\n";
579  run1--;
580  run2--;*/
581  }
582  //cout << "run1-ptrprev " << run1 - ptrprev << "\n";
583 
584 #endif
585  *run2-=sum;
586  if(run1 != run2){
587  *run2 /= *run1;
588  }else{
589  if (this->fCorrectSingular && IsZero(*run2)) {
590  this->fSingular.push_back(col);
591  *run2 = 1.;
592  } else if (IsZero(*run2)) {
593  std::cout << __PRETTY_FUNCTION__ << " Singular pivot at column " << col;
594  DebugStop();
595  }
596 
597  *run2=sqrt(*run2);
598  }
599 
600 }
601 
602 template<class TVar>
603 void TPZSkylParMatrix<TVar>::DecomposeColumnLDLt(int64_t col, int64_t prevcol){
604 
605  //Cholesky Decomposition
606  TVar *ptrprev; //Pointer to prev column
607  TVar *ptrcol; //Pointer to col column
608  int64_t skprev, skcol; //prev and col Skyline height respectively
609  int64_t minline;
610 
611  skprev = this->SkyHeight(prevcol);
612  skcol = this->SkyHeight(col);
613 
614  ptrprev = this->Diag(prevcol);
615  ptrcol = this->Diag(col);
616 
617  if((prevcol-skprev) > (col-skcol)){
618  minline = prevcol - skprev;
619  }else
620  {
621  minline = col - skcol;
622  }
623  if(minline > prevcol) {
624  cout << "error condition col " << col << " prevcol " << prevcol << "\n";
625  cout.flush();
626  return;
627  }
628  int64_t k = minline;
629  TVar *run1 = ptrprev + (prevcol-minline);
630  TVar *run2 = ptrcol + (col-minline);
631  TVar sum = 0;
632 
633  // int64_t n1=1;
634  // int64_t n2=1;
635 
636 
637  while(run1 != ptrprev) sum += (*run1--)*(*run2--)* *this->fElem[k++];
638  *run2-=sum;
639  if(run1 != run2){
640  *run2 /= *run1;
641  }else{
642 #ifdef LOG4CXX
643  if (logger->isDebugEnabled())
644  {
645  std::stringstream sout;
646  sout << "col = " << col << " diagonal " << *run2;
647  LOGPZ_DEBUG(logger,sout.str())
648  }
649 #endif
650  //*run2=sqrt(*run2);
651  }
652 
653 }
654 
655 
656 template<class TVar>
657 void TPZSkylParMatrix<TVar>::DecomposeColumnLDLt2(int64_t col){
658 
659  //Cholesky Decomposition
660  TVar *ptrprev; //Pointer to prev column
661  TVar *ptrcol; //Pointer to col column
662  int64_t skprev, skcol; //prev and col Skyline height respectively
663  int64_t minline;
664  int64_t prevcol;
665  while((fDec[col] < col && fDec[fDec[col]+1] == fDec[col]+1) || fDec[col] == col-1)
666  {
667  prevcol = fDec[col]+1;
668  skprev = this->SkyHeight(prevcol);
669  skcol = this->SkyHeight(col);
670 
671  ptrprev = this->Diag(prevcol);
672  ptrcol = this->Diag(col);
673 
674  if((prevcol-skprev) > (col-skcol)){
675  minline = prevcol - skprev;
676  }else
677  {
678  minline = col - skcol;
679  }
680  if(minline > prevcol) {
681  cout << "error condition col " << col << " prevcol " << prevcol << "\n";
682  cout.flush();
683  return;
684  }
685  int64_t k = minline;
686  TVar *run1 = ptrprev + (prevcol-minline);
687  TVar *run2 = ptrcol + (col-minline);
688  TVar sum = 0;
689 
690  // int64_t n1=1;
691  // int64_t n2=1;
692 
693 
694  while(run1 != ptrprev) sum += (*run1--)*(*run2--)* *this->fElem[k++];
695  *run2-=sum;
696  if(run1 != run2){
697  *run2 /= *run1;
698  }else{
699  if (this->fCorrectSingular && fabs(*run2) < 1.e-10 ) {
700  this->fSingular.push_back(col);
701  *run2 = 1.;
702  }
703 
704 #ifdef LOG4CXX
705  if (logger->isDebugEnabled())
706  {
707  std::stringstream sout;
708  sout << "col = " << col << " diagonal " << *run2;
709  LOGPZ_DEBUG(logger,sout.str())
710  }
711 #endif
712  //*run2=sqrt(*run2);
713  }
714  fDec[col]++;
715  }
716 
717 }
718 
719 template<class TVar>
720 int TPZSkylParMatrix<TVar>::main_nada()
721 {
722 
723  cout << "Calling constructor"<< endl;
724  int64_t tempo1, tempo2, tempo;
725  int64_t dim, nequation, limite, kk, steps;
726  kk=0;
727  cout << "Entre o numero de equacoes\n";
728  cin >> nequation;
729  cout << "Entre com o limite\n";
730  cin >> limite;
731  cout << "Entre com o passo\n";
732  cin >> steps;
733  char filename[20];
734  int tempo_par[50], tempo_sin[50], passos[50];
735  cout << "Entre o nome do Arquivo\n";
736  cin >> filename;
737  ofstream output(filename,ios::app);
738 
739  /*
740  int i;
741  double soma_1=0;
742  double soma_2=0;
743  REAL run1[200];
744  REAL run2[200];
745 
746  const int tempdep = 10;
747 
748 
749  for(i=0;i<tempdep;i++) {
750  run1[i]=i;
751  run2[i]=i;
752  }
753  for(i=0;i<tempdep;i++){
754  soma_1 =soma_1+ run1[i]*run2[i];
755 
756  }
757  soma_2 = soma_2 + TemplateSum<tempdep>(run1,run2);
758 
759 
760 
761  cout << "soma 1 " << soma_1 << endl;
762  cout << "soma 2 " << soma_2 << endl;
763  */
764 
765  for(dim=nequation;dim<=limite;dim=dim+steps)
766  {
767  TPZVec<int64_t> skyline(dim);
768  int64_t i;
769  for(i=0;i<dim;i++) skyline[i] = 0;//random()%(i+1); //0
770  //for(i=0;i<dim;i++) cout << skyline[i] << ' ';
771  cout << endl;
772  TPZSkylParMatrix<TVar> MySky(dim,skyline,2);
773  int64_t j;
774  for(i=0;i<dim;i++) {
775  if(i%100==0) {cout << i << ' '; cout.flush();}
776  for(j=skyline[i];j<=i;j++) {
777  int random = rand();
778  double rnd = (random*10.)/RAND_MAX;
779  MySky(i,j)=rnd;
780  if(i==j) MySky(i,i)=6000.;
781  }
782  }
783  //MySky.SetSkyline(skyline);
784  TPZSkylMatrix<TVar> misael(MySky);
785  // MySky.Print("Gustavo");
786  //misael.Print("Misael");
787 
788  // cout << "Constructed" << endl;
789  //cout << "Calling Decompose_Cholesky routine..."<<endl;
790  time_t t;
791  int nthreads;
792 
793  // cout << "Number of threads to work ? ";
794 
795  // cin >> nthreads;
796  nthreads=2;
797  //cout << endl;
798  if(nthreads < 1) {
799  cout << "VOCE E BURRO MESMO!!!\n";
800  nthreads = 1;
801  }
802  tempo = time(&t);
803  MySky.Decompose_LDLt();
804  tempo1=time(&t)-tempo;
805  MySky.Print("Gustavo",output);
806 
807  tempo = time(&t);
808  misael.Decompose_LDLt();
809  tempo2=time(&t)-tempo;
810 
811 
812  misael.Print("misael",output);
813  passos[kk]=dim;
814  tempo_par[kk]=tempo1;
815  tempo_sin[kk]=tempo2;
816  kk++;
817  }
818  output << "Multi-threaded ";
819  int i;
820  for(i=0;i<kk;i++) output << "{"<< passos[i] << "," << tempo_par[i] << "},";
821  output << endl;
822  output << "Single-Threaded ";
823  for(i=0;i<kk;i++) output << "{"<< passos[i] << "," << tempo_sin[i] << "},";
824  output << endl;
825  output.flush();
826 
827  return 1;
828 }
829 
830 template<class TVar>
831 int TPZSkylParMatrix<TVar>::ClassId() const{
832  return Hash("TPZSkylParMatrix") ^ TPZSkylMatrix<TVar>::ClassId() << 1;
833 }
834 
835 template class TPZSkylParMatrix<float>;
836 template class TPZSkylParMatrix<double>;
837 template class TPZSkylParMatrix<long double>;
838 
839 template class TPZSkylParMatrix<std::complex<float> >;
840 template class TPZSkylParMatrix<std::complex<double> >;
841 template class TPZSkylParMatrix<std::complex<long double> >;
PZ_PTHREAD_MUTEX_UNLOCK
#define PZ_PTHREAD_MUTEX_UNLOCK(mutex, fn)
Definition: pz_pthread.h:79
TPZSkylParMatrix::fColUsed
std::set< int64_t > fColUsed
Definition: pzskylmatpar.h:97
fabs
expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ fabs
Definition: tfadfunc.h:140
TPZSkylParMatrix::ClassId
int ClassId() const override
Define the class id associated with the class.
Definition: pzskylmatpar.cpp:831
PZ_PTHREAD_COND_SIGNAL
#define PZ_PTHREAD_COND_SIGNAL(cond, fn)
Definition: pz_pthread.h:130
stats.filename
filename
Definition: stats.py:82
ELDLt
Definition: pzmatrix.h:52
IsZero
bool IsZero(long double a)
Returns if the value a is close Zero as the allowable tolerance.
Definition: pzreal.h:668
pzlog.h
Contains definitions to LOGPZ_DEBUG, LOGPZ_INFO, LOGPZ_WARN, LOGPZ_ERROR and LOGPZ_FATAL, and the implementation of the inline InitializePZLOG(string) function using log4cxx library or not. It must to be called out of "#ifdef LOG4CXX" scope.
TPZSkylParMatrix::fSkyline
TPZVec< int64_t > fSkyline
Definition: pzskylmatpar.h:92
tempo
TPZTimeTemp tempo
External variable to TPZTimeTemp (to take time)
Definition: TPZTimeTemp.cpp:11
pzskylmatpar.h
Contains TPZSkylParMatrix class which implements a skyline storage format to parallelized process...
TPZSkylParMatrix::Decompose_LDLt
int Decompose_LDLt() override
Decomposes the current matrix using LDLt.
Definition: pzskylmatpar.cpp:453
TPZSkylParMatrix::DecomposeColumnLDLt
void DecomposeColumnLDLt(int64_t lcol, int64_t lprevcol)
Definition: pzskylmatpar.cpp:603
TPZSkylParMatrix::SetSkyline
void SetSkyline(const TPZVec< int64_t > &skyline)
Definition: pzskylmatpar.cpp:88
TPZSkylParMatrix::DecomposeColumnCholesky
void DecomposeColumnCholesky(int64_t lcol, int64_t lprevcol)
Definition: pzskylmatpar.cpp:524
TPZSkylParMatrix::TPZSkylParMatrix
TPZSkylParMatrix()
Default constructor.
Definition: pzskylmatpar.cpp:69
TPZMatrix::fDecomposed
char fDecomposed
Decomposition type used to decompose the current matrix.
Definition: pzmatrix.h:783
skymutex
pthread_mutex_t skymutex
Semaphore.
Definition: pzskylmatpar.cpp:35
TPZSkylMatrix::SkyHeight
int64_t SkyHeight(int64_t col)
return the height of the skyline for a given column
Definition: pzskylmat.h:422
LOGPZ_WARN
#define LOGPZ_WARN(A, B)
Define log for warnings.
Definition: pzlog.h:91
PZ_PTHREAD_COND_WAIT
#define PZ_PTHREAD_COND_WAIT(cond, mutex, fn)
Definition: pz_pthread.h:127
TPZSkylMatrix
Implements a skyline storage format. A Skyline matrix is symmetric so square. Matrix.
Definition: pzskylmat.h:394
TPZSkylMatrix::Diag
TVar * Diag(int64_t col)
This method returns a pointer to the diagonal element of the matrix of the col column.
Definition: pzskylmat.h:577
TPZSkylParMatrix::Decompose_Cholesky
int Decompose_Cholesky() override
Decomposes the current matrix using Cholesky method. The current matrix has to be symmetric...
Definition: pzskylmatpar.cpp:379
PZ_PTHREAD_JOIN
#define PZ_PTHREAD_JOIN(thread, val, fn)
Definition: pz_pthread.h:34
TPZSkylMatrix::SetSkyline
void SetSkyline(const TPZVec< int64_t > &skyline)
modify the skyline of the matrix, throwing away its values skyline indicates the minimum row number w...
Definition: pzskylmat.cpp:1792
nthreads
int nthreads
Definition: numatst.cpp:315
PZ_PTHREAD_MUTEX_LOCK
#define PZ_PTHREAD_MUTEX_LOCK(mutex, fn)
Definition: pz_pthread.h:76
pzfmatrix.h
Contains TPZMatrixclass which implements full matrix (using column major representation).
DebugStop
#define DebugStop()
Returns a message to user put a breakpoint in.
Definition: pzerror.h:20
TPZSkylParMatrix::main_nada
static int main_nada()
Definition: pzskylmatpar.cpp:720
LOGPZ_DEBUG
#define LOGPZ_DEBUG(A, B)
Define log for debug info.
Definition: pzlog.h:87
TPZSkylParMatrix::ParallelLDLt2
static void * ParallelLDLt2(void *t)
Definition: pzskylmatpar.cpp:317
TPZSkylMatrix::Decompose_LDLt
int Decompose_LDLt() override
Decomposes the current matrix using LDLt.
Definition: pzskylmat.cpp:2952
sqrt
expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ sqrt
Definition: tfadfunc.h:120
fblaswr.h
Contains several BLAS rotines for linear algebra.
pz_pthread.h
test.res
string res
Definition: test.py:151
ddot_
doublereal ddot_(integer *N, doublereal *X, integer *incX, doublereal *Y, integer *incY)
pzskylmat.h
Contains TPZSkyline class which implements a skyline storage format.
TPZSkylMatrix::fElem
TPZVec< TVar * > fElem
Storage to keep the first elements to each equation.
Definition: pzskylmat.h:620
LOGPZ_ERROR
#define LOGPZ_ERROR(A, B)
Define log for errors (cout)
Definition: pzlog.h:93
TPZMatrix::Dim
virtual int64_t Dim() const
Returns the dimension of the matrix if the matrix is square.
Definition: pzmatrix.h:892
condition
pthread_cond_t condition
Condition to waiting.
Definition: pzskylmatpar.cpp:37
Hash
int32_t Hash(std::string str)
Definition: TPZHash.cpp:10
TPZSkylParMatrix::~TPZSkylParMatrix
virtual ~TPZSkylParMatrix()
Default destructor.
Definition: pzskylmatpar.cpp:84
TPZSkylParMatrix::Decompose_LDLt
int Decompose_LDLt(std::list< int64_t > &singular) override
Decomposes the current matrix using LDLt. The current matrix has to be symmetric. "L" is lower triangular with 1.0 in its diagonal and "D" is a Diagonal matrix.
Definition: pzskylmatpar.cpp:442
TPZSkylParMatrix::fDec
TPZVec< int > fDec
Definition: pzskylmatpar.h:91
TPZSkylParMatrix
Implements a skyline storage format to parallelized process. Matrix.
Definition: pzskylmatpar.h:35
TPZSkylParMatrix::fSingular
std::list< int64_t > fSingular
Definition: pzskylmatpar.h:100
TPZSkylParMatrix::DecomposeColumnLDLt2
void DecomposeColumnLDLt2(int64_t lcol)
Definition: pzskylmatpar.cpp:657
PZ_PTHREAD_CREATE
#define PZ_PTHREAD_CREATE(thread, attr, routine, args, fn)
Definition: pz_pthread.h:31
TPZSkylParMatrix::ParallelLDLt
static void * ParallelLDLt(void *t)
Definition: pzskylmatpar.cpp:255
TPZSkylMatrix::DecomposeColumn
void DecomposeColumn(int64_t col, int64_t prevcol)
Definition: pzskylmat.cpp:3304
TPZMatrix::Print
virtual void Print(std::ostream &out) const
Definition: pzmatrix.h:253
TPZSkylMatrix::ClassId
int ClassId() const override
Define the class id associated with the class.
Definition: pzskylmat.h:629
ECholesky
Definition: pzmatrix.h:52
TPZSkylParMatrix::ParallelCholesky
static void * ParallelCholesky(void *t)
Definition: pzskylmatpar.cpp:186
TPZSkylParMatrix::ColumnToWork
void ColumnToWork(int64_t &lcol, int64_t &lprevcol)
Determine which column can be decomposed with respect to which column.
Definition: pzskylmatpar.cpp:117
PZ_PTHREAD_COND_BROADCAST
#define PZ_PTHREAD_COND_BROADCAST(cond, fn)
Definition: pz_pthread.h:133
Generated on Sun Aug 16 2020 11:17:05 for NeoPZ by   doxygen 1.8.13