NeoPZ
pzysmp.cpp
Go to the documentation of this file.
1 
6 #include "pzysmp.h"
7 #include "pzfmatrix.h"
8 #include "pzvec.h"
9 
10 #include <memory.h>
11 #include <string>
12 #include <map>
13 #include <pthread.h>
14 
15 #include "tpzverysparsematrix.h"
16 #include "pz_pthread.h"
17 #include "pzstack.h"
18 
19 // #ifdef USING_BLAS
20 // double cblas_ddoti(const int N, const double *X, const int *indx,
21 // const double *Y);
22 // #endif
23 
24 using namespace std;
25 
26 // ****************************************************************************
27 //
28 // Constructors and the destructor
29 //
30 // ****************************************************************************
31 template<class TVar>
32 void TPZFYsmpMatrix<TVar>::InitializeData(){}
33 
34 template<class TVar>
35 void TPZFYsmpMatrix<TVar>::MultiplyDummy(TPZFYsmpMatrix<TVar> & B, TPZFYsmpMatrix<TVar> & Res){
36  int64_t i,j,k;
37  if (B.Rows()!=this->Rows()) return;
38  int64_t rows = this->Rows();
39  REAL aux=0.;
40  for(i=0;i<rows;i++){
41  for(j=0;j<rows;j++){
42  for(k=0;k<rows;k++){
43  // C[i][j] += A[i][k]*B[k][j];
44  aux+=GetVal(i,k)*B.GetVal(i,k);
45  }
46  Res.PutVal(i,j,aux);
47  aux=0.;
48  }
49  }
50 }
51 
52 // ****************************************************************************
53 //
54 // Constructor
55 //
56 // ****************************************************************************
57 
58 template<class TVar>
59 TPZFYsmpMatrix<TVar>::TPZFYsmpMatrix(const TPZVerySparseMatrix<TVar> &cp) : TPZRegisterClassId(&TPZFYsmpMatrix::ClassId),
60 TPZMatrix<TVar>
61 ()
62 {
63  *this = cp;
64 #ifdef USING_MKL
65  fPardisoControl.SetMatrix(this);
66 #endif
67 }
68 
69 template<class TVar>
70 TPZFYsmpMatrix<TVar>::TPZFYsmpMatrix() : TPZRegisterClassId(&TPZFYsmpMatrix::ClassId),
71 TPZMatrix<TVar>(), fIA(1,0),fJA(),fA(),fDiag()
72 {
73 #ifdef USING_MKL
74  fPardisoControl.SetMatrix(this);
75 #endif
76 }
77 
78 template<class TVar>
79 TPZFYsmpMatrix<TVar> &TPZFYsmpMatrix<TVar>::operator=(const TPZFYsmpMatrix<TVar> &cp) {
80  // Deletes everything associated with a TPZFYsmpMatrix
81  TPZMatrix<TVar>::operator=(cp);
82  fIA = cp.fIA;
83  fA = cp.fA;
84  fJA = cp.fJA;
85  fDiag = cp.fDiag;
86 #ifdef USING_MKL
87  fPardisoControl = cp.fPardisoControl;
88  fPardisoControl.SetMatrix(this);
89 #endif
90  return *this;
91 }
92 
93 template<class TVar>
94 TPZFYsmpMatrix<TVar> &TPZFYsmpMatrix<TVar>::operator=(const TPZVerySparseMatrix<TVar> &cp)
95 {
96  // Deletes everything associated with a TPZFYsmpMatrix
97  int64_t nrows = cp.Rows();
98 
99  int64_t count = 0, c = 0, r = 0;
100 
101  count = cp.fExtraSparseData.size();
102  fJA.resize(count);
103  fA.resize(count);
104  fIA.resize(nrows+1);
105  fIA[0] = 0;
106 
107  typename map< pair<int64_t,int64_t>, TVar>::const_iterator it;
108  c = 0;
109  r = 0;
110  for(it=cp.fExtraSparseData.begin(); it!= cp.fExtraSparseData.end(); it++)
111  {
112  int64_t row = it->first.first;
113  if(r != row)
114  {
115  r++;
116  while(r < row)
117  {
118  fIA[r] = c;
119  r++;
120  }
121  fIA[row] = c;
122  r = row;
123  }
124  int64_t col = it->first.second;
125  fJA[c] = col;
126  TVar val = it->second;
127  fA[c] = val;
128  c++;
129  }
130  r++;
131  while(r<=nrows)
132  {
133  fIA[r] = c;
134  r++;
135  }
136  return *this;
137 }
138 
139 template<class TVar>
140 int TPZFYsmpMatrix<TVar>::PutVal(const int64_t row, const int64_t col, const TVar &Value){
141  int64_t k;
142  int flag=0;
143  for(k=fIA[row];k<fIA[row+1];k++){
144  if(fJA[k]==col){
145  flag=1;
146  fA[k]=Value;
147  break;
148  }
149  }
150  if(!flag)
151  {
152  cout << "TPZFYsmpMatrix::PutVal: Non existing position on sparse matrix: line = " << row << " column " << col << endl;
153  DebugStop();
154  return 0;
155  }
156  else
157  {
158  return 1;
159  }
160 }
161 template<class TVar>
162 void TPZFYsmpMatrix<TVar>::AddKel(TPZFMatrix<TVar> & elmat, TPZVec<int64_t> & destinationindex){
163  int64_t i,j,k = 0;
164  TVar value=0.;
165  int64_t ipos,jpos;
166  for(i=0;i<elmat.Rows();i++){
167  for(j=0;j<elmat.Rows();j++){
168  ipos=destinationindex[i];
169  jpos=destinationindex[j];
170  value=elmat.GetVal(i,j);
171  //cout << "j= " << j << endl;
172  if(value != 0.){
173  //cout << "fIA[ipos] " << fIA[ipos] << " fIA[ipos+1] " << fIA[ipos+1] << endl;
174  int flag = 0;
175  k++;
176  if(k >= fIA[ipos] && k < fIA[ipos+1] && fJA[k]==jpos)
177  { // OK -> elements in sequence
178  fA[k]+=value;
179  flag = 1;
180  }else
181  {
182  for(k=fIA[ipos];k<fIA[ipos+1];k++){
183  if(fJA[k]==jpos || fJA[k]==-1){
184  //cout << "fJA[k] " << fJA[k] << " jpos "<< jpos << " " << value << endl;
185  //cout << "k " << k << " "<< jpos << " " << value << endl;
186  flag=1;
187  if(fJA[k]==-1){
188  fJA[k]=jpos;
189  fA[k]=value;
190  // cout << jpos << " " << value << endl;
191  break;
192  }else{
193  fA[k]+=value;
194  break;
195  }
196  }
197  }
198  }
199  if(!flag) cout << "TPZFYsmpMatrix::AddKel: Non existing position on sparse matrix: line =" << ipos << " column =" << jpos << endl; }
200  }
201  }
202 }
203 
204 template<class TVar>
205 void TPZFYsmpMatrix<TVar>::AddKel(TPZFMatrix<TVar> & elmat, TPZVec<int64_t> & sourceindex, TPZVec<int64_t> & destinationindex){
206  int64_t i,j,k = 0;
207  TVar value=0.;
208  int64_t ipos,jpos;
209  for(i=0;i<sourceindex.NElements();i++){
210  for(j=0;j<sourceindex.NElements();j++){
211  ipos=destinationindex[i];
212  jpos=destinationindex[j];
213  value=elmat.GetVal(sourceindex[i],sourceindex[j]);
214  //cout << "j= " << j << endl;
215  if(value != 0.){
216  //cout << "fIA[ipos] " << fIA[ipos] << " fIA[ipos+1] " << fIA[ipos+1] << endl;
217  int flag = 0;
218  k++;
219  if(k >= fIA[ipos] && k < fIA[ipos+1] && fJA[k]==jpos)
220  { // OK -> elements in sequence
221  fA[k]+=value;
222  flag = 1;
223  }else
224  {
225  for(k=fIA[ipos];k<fIA[ipos+1];k++){
226  if(fJA[k]==jpos || fJA[k]==-1){
227  //cout << "fJA[k] " << fJA[k] << " jpos "<< jpos << " " << value << endl;
228  //cout << "k " << k << " "<< jpos << " " << value << endl;
229  flag=1;
230  if(fJA[k]==-1){
231  fJA[k]=jpos;
232  fA[k]=value;
233  // cout << jpos << " " << value << endl;
234  break;
235  }else{
236  fA[k]+=value;
237  break;
238  }
239  }
240  }
241  }
242  if(!flag) cout << "TPZFYsmpMatrix::AddKel: Non existing position on sparse matrix: line =" << ipos << " column =" << jpos << endl; }
243  }
244  }
245 }
246 
247 template<class TVar>
248 void TPZFYsmpMatrix<TVar>::AddKelOld(TPZFMatrix<TVar> & elmat, TPZVec < int > & destinationindex){
249  int64_t i=0;
250  int64_t j=0;
251  int64_t ilocal=0;
252  // int jlocal=0;
253  int64_t nel = destinationindex.NElements();
254  std::multimap<int64_t,int64_t> mapindex;
255  std::multimap<int64_t,int64_t>::iterator hint = mapindex.begin();
256  for(i=0;i<nel;i++){
257  ilocal = destinationindex[i];
258  hint = mapindex.insert(hint,std::make_pair(ilocal,i));
259  // mapindex[ilocal] = i;
260  }
261  for(i=0;i<nel;i++){
262  ilocal = destinationindex[i];
263  int64_t jfirst = fIA[ilocal];
264  int64_t jlast = fIA[ilocal+1];
265  int64_t *Aptr = &fJA[jfirst];
266  int64_t *AptrLast = &fJA[jlast];
267  j=0;
268  std::multimap<int64_t,int64_t>::iterator itelmat = mapindex.begin();
269  while(j<nel && Aptr != AptrLast)
270  {
271  if(*Aptr == (*itelmat).first)
272  {
273  int64_t jel = (*itelmat).second;
274  fA[jfirst] += elmat(i,jel);
275  itelmat++;
276  if(itelmat != mapindex.end() && (*itelmat).second != jel)
277  {
278  Aptr++;
279  jfirst++;
280  }
281  j++;
282  }
283  else if(*Aptr < (*itelmat).first)
284  {
285  Aptr++;
286  jfirst++;
287  }
288  else if(*Aptr > (*itelmat).second)
289  {
290  std::cout << __PRETTY_FUNCTION__ << " inconsistent\n";
291  int64_t *iptr = &fJA[jfirst];
292  while(iptr < AptrLast)
293  {
294  cout << *iptr << " ";
295  iptr++;
296  }
297  cout << endl;
298  std::multimap<int64_t,int64_t>::iterator itelmat2 = mapindex.begin();
299  for(;itelmat2 != mapindex.end(); itelmat2++)
300  {
301  cout << (*itelmat2).first << "/" << (*itelmat2).second << " ";
302  }
303  cout << endl;
304 
305  }
306  }
307  if(j!= nel)
308  {
309  std::cout << __PRETTY_FUNCTION__ << " inconsistent2 j = " << j << " nel " << nel << "\n";
310  int64_t *iptr = &fJA[jfirst];
311  while(iptr < AptrLast)
312  {
313  cout << *iptr << " ";
314  iptr++;
315  }
316  cout << endl;
317  std::multimap<int64_t,int64_t>::iterator itelmat2 = mapindex.begin();
318  for(;itelmat2 != mapindex.end(); itelmat2++)
319  {
320  cout << (*itelmat2).first << "/" << (*itelmat2).second << " ";
321  }
322  cout << endl;
323  }
324  }
325 
326 }
327 
328 template<class TVar>
329 TPZFYsmpMatrix<TVar>::TPZFYsmpMatrix(const int64_t rows,const int64_t cols ) :
330 TPZRegisterClassId(&TPZFYsmpMatrix::ClassId),TPZMatrix<TVar>(rows,cols) {
331  // Constructs an empty TPZFYsmpMatrix
332  // fSolver = -1;
333  fSymmetric = 0;
334  // fMaxIterations = 4;
335  // fSORRelaxation = 1.;
336  fDiag = 0;
337  fA = 0;
338  fIA = 0;
339  fJA = 0;
340 #ifdef USING_MKL
341  fPardisoControl.SetMatrix(this);
342 #endif
343 #ifdef CONSTRUCTOR
344  cerr << "TPZFYsmpMatrix(int rows,int cols)\n";
345 #endif
346 }
347 
348 template<class TVar>
349 TPZFYsmpMatrix<TVar>::~TPZFYsmpMatrix() {
350  // Deletes everything associated with a TPZFYsmpMatrix
351 #ifdef DESTRUCTOR
352  cerr << "~TPZFYsmpMatrix()\n";
353 #endif
354 }
355 
356 // ****************************************************************************
357 //
358 // Find the element of the matrix at (row,col) in the stencil matrix
359 //
360 // ****************************************************************************
361 
362 template<class TVar>
363 const TVar & TPZFYsmpMatrix<TVar>::GetVal(const int64_t row,const int64_t col ) const {
364  // Get the matrix entry at (row,col) without bound checking
365 
366  // Now look through the requested row and see if there is anything
367  // in column col
368  /* int loccol = col+1;
369  for(int ic=fIA[row]-1 ; ic < fIA[row+1]-1; ic++ ) {
370  if ( fJA[ic] == loccol ) return fA[ic];
371  }*/
372  int64_t loccol = col;
373  for(int64_t ic=fIA[row] ; ic < fIA[row+1]; ic++ ) {
374  if ( fJA[ic] == loccol && fJA[ic] != -1 ) return fA[ic];
375  }
376  return this->gZero;
377 }
378 
379 // ****************************************************************************
380 //
381 // Multiply and Multiply-Add
382 //
383 // ****************************************************************************
384 
385 
386 
387 template<class TVar>
388 void TPZFYsmpMatrix<TVar>::MultAddMT(const TPZFMatrix<TVar> &x,const TPZFMatrix<TVar> &y,
389  TPZFMatrix<TVar> &z,
390  const TVar alpha,const TVar beta,const int opt ) {
391  // computes z = beta * y + alpha * opt(this)*x
392  // z and x cannot share storage
393  if(x.Cols() != y.Cols() || x.Cols() != z.Cols() || y.Rows() != z.Rows() )
394  {
395  cout << "\nERROR! in TPZVerySparseMatrix::MultiplyAdd : incompatible dimensions in x, y or z\n";
396  return;
397  }
398 
399  int64_t ir, ic, icol, xcols;
400  xcols = x.Cols();
401  TVar sum;
402  int64_t r = (opt) ? this->Cols() : this->Rows();
403 
404  // Determine how to initialize z
405  for(ic=0; ic<xcols; ic++) {
406  TVar *zp = &(z(0,ic));
407  if(beta != 0) {
408  const TVar *yp = &(y.g(0,0));
409  TVar *zlast = zp+r;
410 
411  if(&z != &y) {
412  memcpy(zp,yp,r*sizeof(TVar));
413  }
414  } else {
415  TVar *zp = &(z(0,0)), *zlast = zp+r;
416  while(zp != zlast) {
417  *zp = 0.;
418  zp ++;
419  }
420  }
421  }
422  // Compute alpha * A * x
423  if(xcols == 1 && opt == 0)
424  {
425  if(this->Cols() != x.Rows() || this->Rows() != y.Rows())
426  {
427  cout << "\nERROR! in TPZFYsmpMatrix::MultiplyAddMT: incompatible dimensions in opt=false\n";
428  return;
429  }
430  for(ir=0; ir<r; ir++) {
431  int64_t icolmin = fIA[ir];
432  int64_t icolmax = fIA[ir+1];
433  const TVar *xptr = &(x.g(0,0));
434  TVar *Aptr = &fA[0];
435  int64_t *JAptr = &fJA[0];
436  for(sum = 0.0, icol=icolmin; icol<icolmax; icol++ ) {
437  sum += Aptr[icol] * xptr[JAptr[icol]];
438  }
439  z(ir,0) += alpha * sum;
440  }
441  }
442  else
443  {
444  for(ic=0; ic<xcols; ic++) {
445  if(opt == 0) {
446 
447  for(ir=0; ir<this->Rows(); ir++) {
448  for(sum = 0.0, icol=fIA[ir]; icol<fIA[ir+1]; icol++ ) {
449  sum += fA[icol] * x.g((fJA[icol]),ic);
450  }
451  z(ir,ic) += alpha * sum;
452  }
453  }
454 
455  // Compute alpha * A^T * x
456  else
457  {
458  if (this->Rows() != x.Rows() || this->Cols() != y.Rows())
459  {
460  cout << "\nERROR! in TPZFYsmpMatrix::MultiplyAddMT: incompatible dimensions in opt=true\n";
461  return;
462  }
463  int64_t jc;
464  int64_t icol;
465  for(ir=0; ir<this->Rows(); ir++) {
466  for(icol=fIA[ir]; icol<fIA[ir+1]; icol++ ) {
467  if(fJA[icol]==-1) break; //Checa a exist�cia de dado ou n�
468  jc = fJA[icol];
469  TVar aval = fA[icol];
470  //cout << "FA["<<icol<<"] = "<<aval<< " * x["<< ir<<"] ="<< x.Get(ir,ic)<< endl;
471  z(jc,ic) += alpha * aval * x.g(ir,ic);
472  }
473  }
474  }
475  }
476  }
477 }
478 
479 // ****************************************************************************
480 //
481 // Multiply and Multiply-Add
482 //
483 // ****************************************************************************
484 
485 template<class TVar>
486 void TPZFYsmpMatrix<TVar>::MultAdd(const TPZFMatrix<TVar> &x,const TPZFMatrix<TVar> &y,
487  TPZFMatrix<TVar> &z,
488  const TVar alpha,const TVar beta,const int opt) const {
489  // computes z = beta * y + alpha * opt(this)*x
490  // z and x cannot share storage
491 
492 #ifdef PZDEBUG
493  if ((!opt && this->Cols() != x.Rows()) || (opt && this->Rows() != x.Rows())) {
494  std::cout << "TPZFMatrix::MultAdd matrix x with incompatible dimensions>" ;
495  return;
496  }
497  if(beta != (double)0. && ((!opt && this->Rows() != y.Rows()) || (opt && this->Cols() != y.Rows()) || y.Cols() != x.Cols())) {
498  std::cout << "TPZFMatrix::MultAdd matrix y with incompatible dimensions>";
499  return;
500  }
501 #endif
502  if(!opt) {
503  if(z.Cols() != x.Cols() || z.Rows() != this->Rows()) {
504  z.Redim(this->Rows(),x.Cols());
505  }
506  } else {
507  if(z.Cols() != x.Cols() || z.Rows() != this->Cols()) {
508  z.Redim(this->Cols(),x.Cols());
509  }
510  }
511  if(this->Cols() == 0) {
512  z.Zero();
513  return;
514  }
515 
516  int64_t ic, xcols;
517  xcols = x.Cols();
518  int64_t r = (opt) ? this->Cols() : this->Rows();
519 
520  // Determine how to initialize z
521  for(ic=0; ic<xcols; ic++) {
522  TVar *zp = &(z(0,ic));
523  if(beta != 0) {
524  const TVar *yp = &(y.g(0,0));
525  TVar *zlast = zp+r;
526  if(beta != 1.) {
527  while(zp < zlast) {
528  *zp = beta * (*yp);
529  zp ++;
530  yp ++;
531  }
532  }
533  else if(&z != &y) {
534  memcpy(zp,yp,r*sizeof(REAL));
535  }
536  } else {
537  TVar *zp = &(z(0,ic)), *zlast = zp+r;
538  while(zp != zlast) {
539  *zp = 0.;
540  zp ++;
541  }
542  }
543  }
544  /*
545  TPZFYsmpMatrix *target;
546  int fFirsteq;
547  int fLasteq;
548  TPZFMatrix<>*fX;
549  TPZFMatrix<>*fZ;
550  REAL fAlpha;
551  int fOpt;
552  */
553  int numthreads = 2;
554  if(opt) numthreads = 1;
555  TPZVec<pthread_t> allthreads(numthreads);
556  TPZVec<TPZMThread> alldata(numthreads);
557  TPZVec<int> res(numthreads);
558  int i;
559  int eqperthread = r/numthreads;
560  int firsteq = 0;
561  for(i=0;i<numthreads;i++)
562  {
563  alldata[i].target = this;
564  alldata[i].fFirsteq = firsteq;
565  alldata[i].fLasteq = firsteq+eqperthread;
566  firsteq += eqperthread;
567  if(i==numthreads-1) alldata[i].fLasteq = this->Rows();
568  alldata[i].fX = &x;
569  alldata[i].fZ = &z;
570  alldata[i].fAlpha = alpha;
571  alldata[i].fOpt = opt;
572  res[i] = PZ_PTHREAD_CREATE(&allthreads[i], NULL,
573  ExecuteMT, &alldata[i], __FUNCTION__);
574  }
575  for(i=0;i<numthreads;i++) {
576  PZ_PTHREAD_JOIN(allthreads[i], NULL, __FUNCTION__);
577  }
578 
579 }
580 
581 // ****************************************************************************
582 //
583 // Print the matrix
584 //
585 // ****************************************************************************
586 
587 template<class TVar>
588 void TPZFYsmpMatrix<TVar>::Print(const char *title, ostream &out ,const MatrixOutputFormat form) const {
589  // Print the matrix along with a identification title
590  if(form != EInputFormat) {
591  out << "\nTFYsmpMatrix Print: " << title << '\n'
592  << "\tRows = " << this->Rows() << '\n'
593  << "\tColumns = " << this->Cols() << '\n';
594  int64_t i;
595  out << "\tIA\tJA\tA\n"
596  << "\t--\t--\t-\n";
597  for(i=0; i<this->Rows(); i++) {
598  out << i << '\t'
599  << fIA[i] << '\t'
600  << fJA[i] << '\t'
601  << fA[i] << '\n';
602  }
603  for(i=this->Rows()+1; i<fIA[this->Rows()]; i++) {
604  out << i << "\t\t"
605  << fJA[i] << '\t'
606  << fA[i] << '\n';
607  }
608  }
609 }
610 
611 
612 // ****************************************************************************
613 //
614 // Various solvers
615 //
616 // ****************************************************************************
617 
618 template<class TVar>
619 void TPZFYsmpMatrix<TVar>::ComputeDiagonal() {
620  if(fDiag.size()) return;
621  int64_t rows = this->Rows();
622  fDiag.resize(rows);
623  for(int64_t ir=0; ir<rows; ir++) {
624  fDiag[ir] = GetVal(ir,ir);
625  }
626 }
627 
628 template<class TVar>
629 void TPZFYsmpMatrix<TVar>::SolveSOR( int64_t &numiterations, const TPZFMatrix<TVar> &rhs, TPZFMatrix<TVar> &x,
630  TPZFMatrix<TVar> *residual, TPZFMatrix<TVar> &/*scratch*/,
631  const REAL overrelax, REAL &tol,
632  const int FromCurrent,const int direction ) {
633 
634  // if(!fDiag) ComputeDiagonal();
635  int64_t irStart = 0,irLast = this->Rows(),irInc = 1;
636  if(direction < 0) {
637  irStart = this->Rows()-1;
638  irLast = -1;
639  irInc = -1;
640  }
641  if(!FromCurrent) x.Zero();
642  TVar eqres = 2.*tol;
643  int64_t iteration;
644  for(iteration=0; iteration<numiterations && eqres >= tol; iteration++) {
645  eqres = 0.;
646  int64_t ir=irStart;
647  while(ir != irLast) {
648  TVar xnewval=rhs.g(ir,0);
649  for(int64_t ic=fIA[ir]; ic<fIA[ir+1]; ic++) {
650  xnewval -= fA[ic] * x(fJA[ic],0);
651  }
652  eqres += xnewval*xnewval;
653  x(ir,0) += overrelax*(xnewval/fDiag[ir]);
654  ir += irInc;
655  }
656  eqres = sqrt(eqres);
657  }
658  tol = eqres;
659  numiterations = iteration;
660  if(residual) this->Residual(x,rhs,*residual);
661 }
662 
663 template<class TVar>
664 int TPZFYsmpMatrix<TVar>::Zero()
665 {
666  fA.Fill(TVar(0.));
667  fDiag.Fill(TVar(0.));
668  return 1;
669 }
670 
671 
682 template<class TVar>
683 void TPZFYsmpMatrix<TVar>::SolveJacobi(int64_t & numiterations, const TPZFMatrix<TVar> & F, TPZFMatrix<TVar> & result, TPZFMatrix<TVar> * residual, TPZFMatrix<TVar> & scratch, REAL & tol, const int FromCurrent)
684 {
685  if(!fDiag.size()) {
686  cout << "TPZSYsmpMatrix::Jacobi cannot be called without diagonal\n";
687  numiterations = 0;
688  if(residual) {
689  this->Residual(result,F,*residual);
690  tol = sqrt(Norm(*residual));
691  }
692  return;
693  }
694  int64_t c = F.Cols();
695  int64_t r = this->Rows();
696  int64_t it=0;
697  if(FromCurrent) {
698  this->Residual(result,F,scratch);
699  for(int64_t ic=0; ic<c; ic++) {
700  for(int64_t i=0; i<r; i++) {
701  result(i,ic) += scratch(i,ic)/(fDiag)[i];
702  }
703  }
704  } else
705  {
706  for(int64_t ic=0; ic<c; ic++) {
707  for(int64_t i=0; i<r; i++) {
708  result(i,ic) = F.GetVal(i,ic)/(fDiag)[i];
709  }
710  }
711  }
712  if(it<numiterations)
713  {
714  this->Residual(result,F,scratch);
715  TVar res = Norm(scratch);
716  for(int64_t it=1; it<numiterations && res > tol; it++) {
717  for(int64_t ic=0; ic<c; ic++) {
718  for(int64_t i=0; i<r; i++) {
719  result(i,ic) += (scratch)(i,ic)/(fDiag)[i];
720  }
721  }
722  this->Residual(result,F,scratch);
723  res = Norm(scratch);
724  }
725  }
726  if(residual) *residual = scratch;
727 }
728 
729 template<class TVar>
730 void *TPZFYsmpMatrix<TVar>::ExecuteMT(void *entrydata)
731 {
732  //DebugStop();
733  TPZMThread *data = (TPZMThread *) entrydata;
734  const TPZFYsmpMatrix *mat = data->target;
735  TVar sum;
736  int64_t xcols = data->fX->Cols();
737  int64_t ic,ir,icol;
738  // Compute alpha * A * x
739  if(xcols == 1 && data->fOpt == 0)
740  {
741  for(ir=data->fFirsteq; ir<data->fLasteq; ir++) {
742  int64_t icolmin = mat->fIA[ir];
743  int64_t icolmax = mat->fIA[ir+1];
744  const TVar *xptr = &(data->fX->g(0,0));
745  TVar *Aptr = &mat->fA[0];
746  int64_t *JAptr = &mat->fJA[0];
747  for(sum = 0.0, icol=icolmin; icol<icolmax; icol++ ) {
748  sum += Aptr[icol] * xptr[JAptr[icol]];
749  }
750  data->fZ->operator()(ir,0) += data->fAlpha * sum;
751  }
752  }
753  else
754  {
755  for(ic=0; ic<xcols; ic++) {
756  if(data->fOpt == 0) {
757 
758  for(ir=data->fFirsteq; ir<data->fLasteq; ir++) {
759  for(sum = 0.0, icol=mat->fIA[ir]; icol<mat->fIA[ir+1]; icol++ ) {
760  sum += mat->fA[icol] * data->fX->g((mat->fJA[icol]),ic);
761  }
762  data->fZ->operator()(ir,ic) += data->fAlpha * sum;
763  }
764  }
765 
766  // Compute alpha * A^T * x
767  else
768  {
769  int64_t jc;
770  int64_t icol;
771  for(ir=data->fFirsteq; ir<data->fLasteq; ir++)
772  {
773  for(icol=mat->fIA[ir]; icol<mat->fIA[ir+1]; icol++ )
774  {
775  if(mat->fJA[icol]==-1) break; //Checa a exist�cia de dado ou n�
776  jc = mat->fJA[icol];
777  data->fZ->operator()(jc,ic) += data->fAlpha * mat->fA[icol] * data->fX->g(ir,ic);
778  }
779  }
780 
781  }
782  }
783  }
784  return 0;
785 }
786 static int FindCol(int64_t *colf, int64_t *coll, int64_t col)
787 {
788  if(col == *colf) return 0;
789  int64_t *begin = colf;
790  int64_t *end = coll;
791  while (begin != end)
792  {
793  int64_t dist = (end-begin)/2;
794  int64_t *mid = begin+dist;
795  if(*mid == col) return (mid-colf);
796  else if(*mid > col) end=mid;
797  else begin = mid;
798  }
799  return -1;
800 }
801 
802 template<class TVar>
803 int TPZFYsmpMatrix<TVar>::GetSub(const int64_t sRow,const int64_t sCol,const int64_t rowSize,
804  const int64_t colSize, TPZFMatrix<TVar> & A ) const {
805  int64_t ir;
806  for(ir=0; ir<rowSize; ir++)
807  {
808  int64_t row = sRow+ir;
809  int64_t colfirst = fIA[row];
810  int64_t collast = fIA[row+1];
811  int64_t iacol = FindCol(&fJA[0]+colfirst,&fJA[0]+collast-1,sCol);
812  int64_t ic;
813  for(ic=0; ic<colSize; ic++) A(ir,ic) = fA[iacol+colfirst];
814  }
815  return 0;
816 }
817 
818 
819 template<class TVar>
820 void TPZFYsmpMatrix<TVar>::GetSub(const TPZVec<int64_t> &indices,TPZFMatrix<TVar> &block) const
821 {
822  std::map<int64_t,int64_t> indord;
823  int64_t i,size = indices.NElements();
824  for(i=0; i<size; i++)
825  {
826  indord[indices[i]] = i;
827  }
828  std::map<int64_t,int64_t>::iterator itset,jtset;
829  for(itset = indord.begin(); itset != indord.end(); itset++)
830  {
831  int64_t *jfirst = &fJA[0]+fIA[(*itset).first];
832  int64_t *jlast = &fJA[0]+fIA[(*itset).first+1]-1;
833  // int row = (*itset).first;
834  for(jtset = indord.begin(); jtset != indord.end(); jtset++)
835  {
836  int64_t col = FindCol(jfirst,jlast,(*jtset).first);
837  int64_t dist = jfirst+col-&fJA[0];
838  block((*itset).second,(*jtset).second) = fA[dist];
839  jfirst += col+1;
840  }
841  }
842 }
843 
844 /*
845  * Perform row update of the sparse matrix
846  */
847 template<class TVar>
848 void TPZFYsmpMatrix<TVar>::RowLUUpdate(int64_t sourcerow, int64_t destrow)
849 {
850  int64_t *sourcefirst = &fJA[0]+fIA[sourcerow];
851  int64_t *sourcelast = &fJA[0]+fIA[sourcerow+1]-1;
852  int64_t sourcecol = FindCol(sourcefirst,sourcelast,sourcerow);
853  if(sourcecol < 0)
854  {
855  cout << __PRETTY_FUNCTION__ << " at line " << __LINE__ << " source not found\n";
856  return;
857  }
858  int64_t sourcedist = sourcefirst+sourcecol-&fJA[0];
859  int64_t *destfirst = &fJA[0]+fIA[destrow];
860  int64_t *destlast = &fJA[0]+fIA[destrow+1]-1;
861  int64_t destcol = FindCol(destfirst,destlast,destrow);
862  int64_t destdist = destfirst+destcol-&fJA[0];
863  if(destcol < 0)
864  {
865  cout << __PRETTY_FUNCTION__ << " at line " << __LINE__ << " destrow not found\n";
866  return;
867  }
868  if(fA[sourcedist] < 1.e-15)
869  {
870  cout << __PRETTY_FUNCTION__ << " at line " << __LINE__ << " small pivot " << fA[sourcedist] << "\n";
871  return;
872  }
873  TVar mult = fA[destdist]/fA[sourcedist];
874  if(mult == 0.) return;
875  destdist++;
876  sourcedist++;
877  while(destdist < fIA[destrow+1] && sourcedist < fIA[sourcerow+1])
878  {
879  if(fJA[destdist] == fJA[sourcedist])
880  {
881  fA[destdist] -= fA[sourcedist]*mult;
882  destdist++;
883  sourcedist++;
884  }
885  else if(fJA[destdist] < fJA[sourcedist])
886  {
887  destdist++;
888  }
889  else
890  {
891  sourcedist++;
892  }
893  }
894 
895 }
898 template<class TVar>
899 void TPZFYsmpMatrix<TVar>::AutoFill(int64_t nrow, int64_t ncol, int symmetric)
900 {
901  if (symmetric && nrow != ncol) {
902  DebugStop();
903  }
904  TPZFMatrix<TVar> orig;
905  orig.AutoFill(nrow,ncol,symmetric);
906 
907  TPZVec<int64_t> IA(nrow+1);
908  TPZStack<int64_t> JA;
909  TPZStack<TVar> A;
910  IA[0] = 0;
911  TPZVec<std::set<int64_t> > eqs(nrow);
912  for (int64_t row=0; row<nrow; row++) {
913  if(nrow == ncol) eqs[row].insert(row);
914  for (int64_t col = 0; col<ncol; col++) {
915  REAL test = rand()*1./RAND_MAX;
916  if (test > 0.5) {
917  eqs[row].insert(col);
918  if (symmetric) {
919  eqs[col].insert(row);
920  }
921  }
922  }
923  }
924  int64_t pos=0;
925  for (int64_t row=0; row< nrow; row++) {
926  for (std::set<int64_t>::iterator col = eqs[row].begin(); col != eqs[row].end(); col++) {
927  JA.Push(*col);
928  A.Push(orig(row,*col));
929  }
930  IA[row+1] = JA.size();
931  }
932  TPZMatrix<TVar>::Resize(nrow,ncol);
933  SetData(IA, JA, A);
934 }
935 
936 
940 template<class TVar>
941 int TPZFYsmpMatrix<TVar>::Decompose_LU(std::list<int64_t> &singular)
942 {
943  return Decompose_LU();
944 }
945 template<class TVar>
946 int TPZFYsmpMatrix<TVar>::Decompose_LU()
947 {
948 
949 #ifdef USING_MKL
950  if(this->IsDecomposed() == ELU) return 1;
951  if (this->IsDecomposed() != ENoDecompose) {
952  DebugStop();
953  }
954  fPardisoControl.SetMatrixType(TPZPardisoControl<TVar>::ENonSymmetric,TPZPardisoControl<TVar>::EIndefinite);
955  fPardisoControl.Decompose();
956  this->SetIsDecomposed(ELU);
957  return 1;
958 #endif
959 
960  int64_t row;
961  int64_t neq = this->Rows();
962  for(row=1; row<neq; row++)
963  {
964  // int firstcol = fIA[row];
965  int64_t lastcol = fIA[row+1];
966  int64_t colind = 0;
967  if(fJA[lastcol-1] < row) continue;
968  while(fJA[colind] < row)
969  {
970  RowLUUpdate(fJA[colind],row);
971  colind++;
972  }
973  }
974  this->fDecomposed=1;
975  return 1;
976 }
977 
978 template<class TVar>
979 int TPZFYsmpMatrix<TVar>::Substitution( TPZFMatrix<TVar> *B ) const
980 {
981 
982 #ifdef USING_MKL
983  TPZFMatrix<TVar> x(*B);
984  fPardisoControl.Solve(*B,x);
985  *B = x;
986  return 1;
987 #endif
988 
989  int64_t row;
990  int64_t bcol = B->Cols();
991  int64_t col;
992  int64_t neq = this->Rows();
993 
994  // forward substitution
995  for(row=0; row<neq; row++)
996  {
997  int64_t firstrow = fIA[row];
998  int64_t lastrow = fIA[row+1];
999  if(fJA[firstrow] > row || fJA[lastrow-1] < row)
1000  {
1001  cout << __PRETTY_FUNCTION__ << " " << __LINE__ << " inconsistent column information for row " << row << endl;
1002  continue;
1003  }
1004  int64_t rowcounter = firstrow;
1005  while(fJA[rowcounter] < row)
1006  {
1007  for(col=0; col<bcol; col++)
1008  {
1009  (*B)(row,col) -= fA[rowcounter]*(*B)(fJA[rowcounter],col);
1010  }
1011  }
1012  for(col=0; col<bcol; col++)
1013  {
1014  (*B)(row,col) /= fA[rowcounter];
1015  }
1016  }
1017  // backward substitution
1018  for(row = neq-1; row >= 0; row--)
1019  {
1020  int64_t firstrow = fIA[row];
1021  int64_t lastrow = fIA[row+1];
1022  int64_t col = FindCol(&fJA[0]+firstrow,&fJA[0]+lastrow-1,row);
1023  if(col < 0)
1024  {
1025  cout << __PRETTY_FUNCTION__ << " " << __LINE__ << " inconsistent column information for row " << row << endl;
1026  continue;
1027  }
1028  int64_t coldist = firstrow+col+1;
1029  while(coldist < lastrow)
1030  {
1031  for(col=0; col<bcol; col++)
1032  {
1033  (*B)(row,col) -= fA[coldist]*(*B)(fJA[coldist],col);
1034  }
1035  }
1036  }
1037  return 1;
1038 }
1039 
1040 
1041 template<class TVar>
1042 int TPZFYsmpMatrix<TVar>::ClassId() const{
1043  return Hash("TPZFYsmpMatrix") ^ TPZMatrix<TVar>::ClassId() << 1;
1044 }
1045 template class TPZFYsmpMatrix<long double>;
1046 template class TPZFYsmpMatrix<double>;
1047 template class TPZFYsmpMatrix<float>;
TPZFYsmpMatrix::AutoFill
void AutoFill(int64_t nrow, int64_t ncol, int symmetric)
Fill matrix storage with randomic values.
Definition: pzysmp.cpp:899
TPZFYsmpMatrix::fA
TPZVec< TVar > fA
Definition: pzysmp.h:208
TPZFYsmpMatrix::MultiplyDummy
void MultiplyDummy(TPZFYsmpMatrix< TVar > &B, TPZFYsmpMatrix< TVar > &Res)
Definition: pzysmp.cpp:35
TPZFYsmpMatrix::Substitution
virtual int Substitution(TPZFMatrix< TVar > *B) const override
Computes Forward and Backward substitution for a "LU" decomposed matrix.
Definition: pzysmp.cpp:979
TPZFYsmpMatrix::TPZMThread
An auxiliary structure to hold the data of the subset of equations used to multiply in a multi-thre...
Definition: pzysmp.h:57
TPZMatrix::Residual
virtual void Residual(const TPZFMatrix< TVar > &x, const TPZFMatrix< TVar > &rhs, TPZFMatrix< TVar > &res)
Computes res = rhs - this * x.
Definition: pzmatrix.h:827
ENoDecompose
Definition: pzmatrix.h:52
TPZFYsmpMatrix::Print
virtual void Print(const char *title, std::ostream &out=std::cout, const MatrixOutputFormat form=EFormatted) const override
Print the matrix along with a identification title.
Definition: pzysmp.cpp:588
TPZFYsmpMatrix::AddKel
virtual void AddKel(TPZFMatrix< TVar > &elmat, TPZVec< int64_t > &destinationindex) override
Add a contribution of a stiffness matrix.
Definition: pzysmp.cpp:162
TPZVec::resize
virtual void resize(const int64_t newsize)
Definition: pzvec.h:213
TPZVerySparseMatrix::fExtraSparseData
std::map< std::pair< int64_t, int64_t >, TVar > fExtraSparseData
Save elements different from zero, of Sparse matrix.
Definition: tpzverysparsematrix.h:122
TPZFYsmpMatrix::AddKelOld
virtual void AddKelOld(TPZFMatrix< TVar > &elmat, TPZVec< int > &destinationindex)
Add a contribution of a stiffness matrix putting it on destination indexes position.
Definition: pzysmp.cpp:248
EInputFormat
Definition: pzmatrix.h:55
pzvec.h
Templated vector implementation.
test
Definition: test.py:1
TPZFYsmpMatrix::PutVal
int PutVal(const int64_t row, const int64_t col, const TVar &Value) override
Put values without bounds checking This method is faster than "Put" if DEBUG is defined.
Definition: pzysmp.cpp:140
MatrixOutputFormat
MatrixOutputFormat
Defines output format.
Definition: pzmatrix.h:55
TPZMatrix::Resize
virtual int Resize(const int64_t newRows, const int64_t newCols)
Redimensions a matriz keeping the previous values.
Definition: pzmatrix.h:278
TPZFYsmpMatrix::RowLUUpdate
void RowLUUpdate(int64_t sourcerow, int64_t destrow)
Definition: pzysmp.cpp:848
TPZMatrix::SetIsDecomposed
void SetIsDecomposed(int val)
Sets current matrix to decomposed state.
Definition: pzmatrix.h:410
TPZFYsmpMatrix::~TPZFYsmpMatrix
virtual ~TPZFYsmpMatrix()
Definition: pzysmp.cpp:349
TPZMatrix::fDecomposed
char fDecomposed
Decomposition type used to decompose the current matrix.
Definition: pzmatrix.h:783
TPZFYsmpMatrix
Implements a non symmetric sparse matrix (Yale Sparse Matrix Storage). Matrix.
Definition: pzysmp.h:45
TPZTransform::ClassId
int ClassId() const override
Define the class id associated with the class.
Definition: pztrnsform.h:86
val
REAL val(STATE &number)
Returns value of the variable.
Definition: pzartdiff.h:23
TPZFYsmpMatrix::Zero
virtual int Zero() override
Zeroes the matrix.
Definition: pzysmp.cpp:664
TPZFYsmpMatrix::SetData
virtual void SetData(int64_t *IA, int64_t *JA, TVar *A)
Pass the data to the class.
Definition: pzysmp.h:232
TPZMatrix::gZero
static TVar gZero
Initializing value to static variable.
Definition: pzmatrix.h:786
TPZFYsmpMatrix::MultAdd
virtual void MultAdd(const TPZFMatrix< TVar > &x, const TPZFMatrix< TVar > &y, TPZFMatrix< TVar > &z, const TVar alpha=1., const TVar beta=0., const int opt=0) const override
It computes z = beta * y + alpha * opt(this)*x but z and x can not overlap in memory.
Definition: pzysmp.cpp:486
PZ_PTHREAD_JOIN
#define PZ_PTHREAD_JOIN(thread, val, fn)
Definition: pz_pthread.h:34
TPZFMatrix::Zero
int Zero() override
Makes Zero all the elements.
Definition: pzfmatrix.h:651
TPZVec::size
int64_t size() const
Returns the number of elements of the vector.
Definition: pzvec.h:196
TPZFYsmpMatrix::InitializeData
void InitializeData()
Implements a initialization method for the sparse structure. It sets the initial value for the fIA an...
Definition: pzysmp.cpp:32
pzgeom::tol
static const double tol
Definition: pzgeoprism.cpp:23
TPZStack::Push
void Push(const T object)
Pushes a copy of the object on the stack.
Definition: pzstack.h:80
TPZFYsmpMatrix::Decompose_LU
virtual int Decompose_LU() override
Definition: pzysmp.cpp:946
TPZMatrix::IsDecomposed
int IsDecomposed() const
Checks if current matrix is already decomposed.
Definition: pzmatrix.h:405
Norm
TVar Norm(const TPZFMatrix< TVar > &A)
Returns the norm of the matrix A.
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
TPZVerySparseMatrix
Implements a matrix whose nonzero elements are stored in binary tree. Matrix.
Definition: pzfmatrix.h:33
substruct_tst14.test.numthreads
list numthreads
Definition: test.py:210
ELU
Definition: pzmatrix.h:52
TPZFYsmpMatrix::ExecuteMT
static void * ExecuteMT(void *entrydata)
Definition: pzysmp.cpp:730
TPZMatrix::Rows
int64_t Rows() const
Returns number of rows.
Definition: pzmatrix.h:803
pzysmp.h
Contains the TPZFYsmpMatrix class which implements a non symmetric sparse matrix. ...
TPZFYsmpMatrix::fDiag
TPZVec< TVar > fDiag
Definition: pzysmp.h:210
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
pz_pthread.h
test.res
string res
Definition: test.py:151
TPZFYsmpMatrix::fSymmetric
int fSymmetric
Definition: pzysmp.h:212
TPZVec::begin
T * begin() const
Casting operator. Returns The fStore pointer.
Definition: pzvec.h:450
TPZFMatrix
Full matrix class. Matrix.
Definition: pzfmatrix.h:32
TPZFMatrix::Redim
int Redim(const int64_t newRows, const int64_t newCols) override
Redimension a matrix and ZERO your elements.
Definition: pzfmatrix.h:616
TPZFYsmpMatrix::MultAddMT
virtual void MultAddMT(const TPZFMatrix< TVar > &x, const TPZFMatrix< TVar > &y, TPZFMatrix< TVar > &z, const TVar alpha=1., const TVar beta=0., const int opt=0)
Definition: pzysmp.cpp:388
pzstack.h
A simple stack.
Hash
int32_t Hash(std::string str)
Definition: TPZHash.cpp:10
TPZFYsmpMatrix::GetSub
virtual int GetSub(const int64_t sRow, const int64_t sCol, const int64_t rowSize, const int64_t colSize, TPZFMatrix< TVar > &A) const override
Gets submatrix storing it on Target.
Definition: pzysmp.cpp:803
dist
REAL dist(TPZVec< T1 > &vec1, TPZVec< T1 > &vec2)
Definition: pzvec_extras.h:124
TPZFYsmpMatrix::GetVal
virtual const TVar & GetVal(const int64_t row, const int64_t col) const override
Get the matrix entry at (row,col) without bound checking.
Definition: pzysmp.cpp:363
TPZFYsmpMatrix::fJA
TPZVec< int64_t > fJA
Definition: pzysmp.h:207
TPZMatrix::AutoFill
void AutoFill(int64_t nrow, int64_t ncol, int symmetric)
Fill matrix storage with randomic values.
Definition: pzmatrix.cpp:1960
PZ_PTHREAD_CREATE
#define PZ_PTHREAD_CREATE(thread, attr, routine, args, fn)
Definition: pz_pthread.h:31
TPZVec::Fill
void Fill(const T &copy, const int64_t from=0, const int64_t numelem=-1)
Will fill the elements of the vector with a copy object.
Definition: pzvec.h:460
tpzverysparsematrix.h
Contains TPZVerySparseMatrix class which implements a matrix whose nonzero elements are stored in bin...
TPZMatrix::Cols
int64_t Cols() const
Returns number of cols.
Definition: pzmatrix.h:809
FindCol
static int FindCol(int64_t *colf, int64_t *coll, int64_t col)
Definition: pzysmp.cpp:786
TPZVec::NElements
int64_t NElements() const
Returns the number of elements of the vector.
Definition: pzvec.h:190
TPZFYsmpMatrix::fIA
TPZVec< int64_t > fIA
Definition: pzysmp.h:206
TPZFYsmpMatrix::TPZMThread::fX
const TPZFMatrix< TVar > * fX
Definition: pzysmp.h:61
TPZFYsmpMatrix::SolveJacobi
virtual void SolveJacobi(int64_t &numiterations, const TPZFMatrix< TVar > &F, TPZFMatrix< TVar > &result, TPZFMatrix< TVar > *residual, TPZFMatrix< TVar > &scratch, REAL &tol, const int FromCurrent=0) override
Solves the linear system using Jacobi method. .
Definition: pzysmp.cpp:683
TPZFYsmpMatrix::ClassId
int ClassId() const override
Define the class id associated with the class.
Definition: pzysmp.cpp:1042
TPZFMatrix::GetVal
const TVar & GetVal(const int64_t row, const int64_t col) const override
Get values without bounds checking This method is faster than "Get" if DEBUG is defined.
Definition: pzfmatrix.h:566
TPZFYsmpMatrix::SolveSOR
void SolveSOR(int64_t &numiterations, const TPZFMatrix< TVar > &rhs, TPZFMatrix< TVar > &x, TPZFMatrix< TVar > *residual, TPZFMatrix< TVar > &scratch, const REAL overrelax, REAL &tol, const int FromCurrent=0, const int direction=1) override
Solves the linear system using Successive Over Relaxation method (Gauss Seidel). ...
Definition: pzysmp.cpp:629
TPZFYsmpMatrix::TPZMThread::target
const TPZFYsmpMatrix< TVar > * target
Definition: pzysmp.h:58
TPZFYsmpMatrix::operator=
TPZFYsmpMatrix & operator=(const TPZFYsmpMatrix< TVar > &copy)
Definition: pzysmp.cpp:79
TPZVec::end
T * end() const
Returns a pointer to the last+1 element.
Definition: pzvec.h:455
TPZFYsmpMatrix::TPZMThread::fZ
TPZFMatrix< TVar > * fZ
Definition: pzysmp.h:62
TPZFMatrix::g
TVar & g(const int64_t row, const int64_t col) const
Definition: pzfmatrix.h:594
TPZMatrix::ClassId
int ClassId() const override
Define the class id associated with the class.
Definition: pzmatrix.h:957
TPZMatrix
Root matrix class (abstract). Matrix.
Definition: pzmatrix.h:60
TPZFYsmpMatrix::ComputeDiagonal
void ComputeDiagonal()
Definition: pzysmp.cpp:619
Generated on Sun Aug 16 2020 11:17:05 for NeoPZ by   doxygen 1.8.13