NeoPZ
pzmatrix.cpp
Go to the documentation of this file.
1 
6 #include <stdio.h>
7 #include <stdlib.h>
8 #include <math.h>
9 #include <memory.h>
10 #include <sstream>
11 #include <set>
12 
13 //#include "pzfmatrix.h"
14 #include "pzmatrix.h"
15 #include "pzsolve.h"
16 #include "pzvec.h"
17 #include "pzextractval.h"
18 
19 #ifdef _AUTODIFF
20 #include "fadType.h"
21 #include "fad.h"
22 #endif
23 
24 #include <sstream>
25 #include <exception>
26 #include "pzlog.h"
27 #include <complex>
28 
29 #ifdef LOG4CXX
30 static LoggerPtr logger(Logger::getLogger("pz.matrix.tpzmatrix"));
31 static LoggerPtr loggerCheck(Logger::getLogger("pz.checkconsistency"));
32 #endif
33 
34 #ifdef PZDEBUG
35 #define DEBUG2
36 #endif
37 
38 using namespace std;
39 //template <class TVar>
40 //TVar TPZMatrix<TVar>::gZero = TVar(0);
41 
42 template <class TVar>
43 TPZMatrix<TVar>::~TPZMatrix()
44 {
45  fDecomposed = 0;
46  fDefPositive = 0;
47  fRow = 0;
48  fCol = 0;
49 }
50 
55 template<class TVar>
56 const TVar &TPZMatrix<TVar>::GetVal(const int64_t /*row*/, const int64_t /*col*/ ) const
57 {
58  return gZero;
59 }
60 
61 
62 
63 template<class TVar>
64 void TPZMatrix<TVar>::Add(const TPZMatrix<TVar>&A,TPZMatrix<TVar>&B) const {
65  if ((Rows() != A.Rows()) || (Cols() != A.Cols()) ) {
66  Error( "Add(TPZMatrix<>&, TPZMatrix) <different dimensions>" );
67  }
68 
69  B.Redim( A.Rows(), A.Cols() );
70 
71  for ( int64_t r = 0; r < Rows(); r++ )
72  for ( int64_t c = 0; c < Cols(); c++ ) B.PutVal( r, c, GetVal(r,c)+A.GetVal(r,c) );
73 }
74 template<class TVar>
75 void TPZMatrix<TVar>::Substract(const TPZMatrix<TVar> &A,TPZMatrix<TVar> &result) const {
76 
77  if ((Rows() != A.Rows()) || (Cols() != A.Cols()) ) {
78  Error( "Add(TPZMatrix<>&, TPZMatrix<>&) <different dimensions>" );
79  }
80 
81  result.Resize( Rows(), Cols() );
82  for ( int64_t r = 0; r < Rows(); r++ ) {
83  for ( int64_t c = 0; c < Cols(); c++ ) {
84  result.PutVal( r, c, GetVal(r,c)-A.GetVal(r,c) );
85  }
86  }
87 }
88 
89 template<class TVar>
90 void TPZMatrix<TVar>::Simetrize() {
91 
92  if ( Rows() != Cols() ) {
93  Error( "Simetrize only work for square matrices" );
94  }
95 
96  int64_t row,col;
97  int64_t fDim1 = Rows();
98  for(row=0; row<fDim1; row++) {
99  for(col=row+1; col<fDim1; col++) {
100  this->s(col,row) = this->s(row,col);
101  }
102  }
103 
104 }
105 
107 template<class TVar>
108 TPZFMatrix<TVar> operator+(const TPZMatrix<TVar> &A, const TPZMatrix<TVar> &B ) {
109  TPZFMatrix<TVar> temp;
110  temp.Redim( A.Rows(), A.Cols() );
111  A.Add(B,temp);
112  return temp;
113 }
114 
115 
117 template<class TVar>
118 TPZFMatrix<TVar> operator-(const TPZMatrix<TVar> &A, const TPZMatrix<TVar> &B ) {
119  TPZFMatrix<TVar> temp;
120  TPZFMatrix<TVar> res;
121  res.Redim( A.Rows(), A.Cols() );
122  A.Substract(B,res);
123  return temp;
124 }
125 
127 template<class TVar>
128 TPZFMatrix<TVar> operator*( TPZMatrix<TVar> &A, const TPZFMatrix<TVar> &B ) {
129  TPZFMatrix<TVar> res;
130  res.Redim( A.Rows(), B.Cols() );
131  A.Multiply(B,res);
132  return res;
133 }
134 template<class TVar>
135 void TPZMatrix<TVar>::PrepareZ(const TPZFMatrix<TVar> &y, TPZFMatrix<TVar> &z,const TVar beta,const int opt) const
136 {
137  int64_t numeq = (opt) ? Cols() : Rows();
138  int64_t xcols = y.Cols();
139  int64_t ic;
140  if(!z.Rows()) return;
141  for (ic = 0; ic < xcols; ic++)
142  {
143  TVar *zp = &z(0,ic), *zlast = zp+numeq;
144  if(beta != (TVar)0)
145  {
146  const TVar *yp = &y.g(0,ic);
147  if(&z != &y)
148  {
149  memcpy((void *)(zp),(void *)(yp),numeq*sizeof(TVar));
150  }
151  for (int64_t i=0; i<numeq; i++) {
152  z(i,ic) *= beta;
153  }
154  } else
155  {
156  while(zp != zlast)
157  {
158  *zp = 0.;
159  zp ++;
160  }
161  }
162  }
163 }
164 
165 template<class TVar>
166 void TPZMatrix<TVar>::MultAdd(const TPZFMatrix<TVar> &x,const TPZFMatrix<TVar> &y, TPZFMatrix<TVar> &z, const TVar alpha,const TVar beta,const int opt) const {
167  if ((!opt && Cols() != x.Rows()) || Rows() != x.Rows())
168  Error( "Operator* <matrixs with incompatible dimensions>" );
169  if(x.Cols() != y.Cols() || x.Cols() != z.Cols() || x.Rows() != y.Rows() || x.Rows() != z.Rows()) {
170  Error ("TPZFMatrix::MultiplyAdd incompatible dimensions\n");
171  }
172  int64_t rows = Rows();
173  int64_t cols = Cols();
174  int64_t xcols = x.Cols();
175  int64_t ic, c, r;
176  PrepareZ(y,z,beta,opt);
177  TVar val = 0.;
178  for (ic = 0; ic < xcols; ic++) {
179  if(!opt) {
180  for ( c = 0; c<cols; c++) {
181  for ( r = 0; r < rows; r++ ) {
182  val = z(r,ic) + alpha * GetVal(r,c) * x.GetVal(c,ic);
183  z.PutVal(r,ic,val);
184  }
185  }
186  } else {
187  for (r = 0; r<rows; r++) {
188  val = 0.;
189  for(c = 0; c<cols; c++) {
190  val += GetVal(c,r)* x.GetVal(c,ic);
191  }
192  z.PutVal(r,ic,alpha*val);
193  }
194  }
195  }
196 }
197 
198 template<class TVar>
199 void TPZMatrix<TVar>::Identity() {
200 
201  if ( Cols() != Rows() ) {
202  Error( "identity (TPZMatrix<>*) <TPZMatrix<>must be square>" );
203  }
204  for ( int64_t row = 0; row < Rows(); row++) {
205  for ( int64_t col = 0; col < Cols(); col++ ) {
206  (row == col)? PutVal(row,col,1.):PutVal(row,col,0.);
207  }
208  }
209 }
210 
211 /*************/
212 /*** Input ***/
213 #ifdef _AUTODIFF
214 template<>
215 void TPZMatrix<TFad<6,REAL> >::Input(std::istream& in )
216 {
217  DebugStop();
218 }
219 template<>
220 void TPZMatrix<Fad<REAL> >::Input(std::istream& in )
221 {
222  DebugStop();
223 }
224 #endif
225 
226 template<class TVar>
227 void TPZMatrix<TVar>::Input(std::istream& in )
228 {
229 
230  int64_t newRow, newCol;
231  in >> newRow;
232  in >> newCol;
233  Redim( newRow, newCol );
234  int64_t i,j;
235  TVar elem;
236  for(i=0;i<Rows();i++)
237  for(j=0;j<Cols();j++)
238  {
239  in >> elem;
240  Put( i,j, elem );
241  }
242 }
243 
245 template<class TVar>
246 std::istream & operator>>(std::istream& in,TPZMatrix<TVar> &A)
247 {
248  A.Input(in);
249  return in;
250 }
251 
252 
253 /*************/
254 /*** Print ***/
255 #ifdef _AUTODIFF
256 template<>
257 void TPZMatrix<TFad<6,REAL> >::Print(const char *name, std::ostream& out,const MatrixOutputFormat form) const
258 {
259  out << name << std::endl;
260  for (int64_t i=0; i<fRow; i++) {
261  for (int64_t j=0; j<fCol; j++) {
262  out << "i = " << i << " j = " << j << " val " << Get(i, j) << std::endl;
263  }
264  }
265 }
266 #endif
267 
268 
269 template<class TVar>
270 void TPZMatrix<TVar>::Print(const char *name, std::ostream& out,const MatrixOutputFormat form) const {
271 
272  if(form == EFormatted) {
273  out << "Writing matrix '";
274  if(name) out << name;
275  out << "' (" << Rows() << " x " << Cols() << "):\n";
276 
277  for ( int64_t row = 0; row < Rows(); row++) {
278  out << "\t";
279  for ( int64_t col = 0; col < Cols(); col++ ) {
280  out << Get( row, col) << " ";
281  }
282  out << "\n";
283  }
284  out << "\n";
285  } else if (form == EInputFormat) {
286  out << Rows() << " " << Cols() << endl;
287  for ( int64_t row = 0; row < Rows(); row++) {
288  for ( int64_t col = 0; col < Cols(); col++ ) {
289  TVar val = Get (row, col);
290  if(val != (TVar)0.) out << row << ' ' << col << ' ' << val << std::endl;
291  }
292  }
293  out << "-1 -1 0.\n";
294  } else if( form == EMathematicaInput)
295  {
296  char number[256];
297  out << name << "\n{ ";
298  for ( int64_t row = 0; row < Rows(); row++) {
299  out << "\n{ ";
300  for ( int64_t col = 0; col < Cols(); col++ ) {
301  TVar val = Get (row, col);
302  #ifdef STATE_COMPLEX
303  sprintf(number, "%16.16Lf",(long double) fabs(val));
304  #else
305  sprintf(number, "%16.14Lf",(long double) TPZExtractVal::val(val) );
306  #endif
307  out << number;
308  if(col < Cols()-1)
309  out << ", ";
310  if((col+1) % 6 == 0)out << std::endl;
311  }
312  out << " }";
313  if(row < Rows()-1)
314  out << ",";
315  }
316 
317  out << " };\n";
318 
319  }else if( form == EMatlabNonZeros)
320  {
321  out << name;
322  for ( int64_t row = 0; row < Rows(); row++) {
323  out << "\n|";
324  for ( int64_t col = 0; col < Cols(); col++ )
325  if(IsZero(Get (row, col)) ){
326  out << ".";
327  }else{
328  out << "#";
329  }
330  out << "|";
331  }
332  out << "\n";
333  }
334  else if( form == EMatrixMarket)
335  {
336  bool sym = IsSimetric();
337  int64_t numzero = 0;
338  int64_t nrow = Rows();
339  for ( int64_t row = 0; row < Rows(); row++) {
340  int64_t colmax = nrow;
341  if (sym) colmax = row+1;
342  for (int64_t col = 0; col < colmax; col++ )
343  {
344  TVar val = GetVal(row, col);
345  if (val != (TVar)0.) {
346  numzero++;
347  }
348  }
349  }
350  out << "%%"<< name << std::endl;
351  out << Rows() << ' ' << Cols() << ' ' << numzero << std::endl;
352  for ( int64_t row = 0; row < Rows(); row++) {
353  int64_t colmax = nrow;
354  if (sym) colmax = row+1;
355  for (int64_t col = 0; col < colmax; col++ )
356  {
357  TVar val = GetVal(row, col);
358  if (val != (TVar)0.) {
359  out << row+1 << ' ' << col+1 << ' ' << val << std::endl;
360  }
361  }
362  }
363  }
364 }
365 
366 template<>
367 void TPZMatrix<std::complex<float> >::Print(const char *name, std::ostream& out,const MatrixOutputFormat form) const {
368  if(form == EFormatted) {
369  out << "Writing matrix '";
370  if(name) out << name;
371  out << "' (" << Rows() << " x " << Cols() << "):\n";
372 
373  for ( int64_t row = 0; row < Rows(); row++) {
374  out << "\t";
375  for ( int64_t col = 0; col < Cols(); col++ ) {
376  out << Get( row, col).real() << " " << Get(row, col).imag() << " ";
377  }
378  out << "\n";
379  }
380  out << "\n";
381  } else if (form == EInputFormat) {
382  out << Rows() << " " << Cols() << endl;
383  for ( int64_t row = 0; row < Rows(); row++) {
384  for ( int64_t col = 0; col < Cols(); col++ ) {
385  std::complex<double> val = Get (row, col);
386  if(val != 0.) out << row << ' ' << col << ' ' << val.real() << " " << val.imag() << std::endl;
387  }
388  }
389  out << "-1 -1 0.\n";
390  } else if( form == EMathematicaInput)
391  {
392  char number[128];
393  out << name << "\n{ ";
394  for ( int64_t row = 0; row < Rows(); row++) {
395  out << "\n{ ";
396  for ( int64_t col = 0; col < Cols(); col++ ) {
397  std::complex<double> val = Get (row, col);
398 
399  sprintf(number, "%16.16lf + I %16.16lf", val.real(), val.imag());
400 
401 
402  out << number;
403  if(col < Cols()-1)
404  out << ", ";
405  if((col+1) % 6 == 0)out << std::endl;
406  }
407  out << " }";
408  if(row < Rows()-1)
409  out << ",";
410  }
411 
412  out << " };\n";
413 
414  }else if( form == EMatlabNonZeros)
415  {
416  out << name;
417  for ( int64_t row = 0; row < Rows(); row++) {
418  out << "\n|";
419  for ( int64_t col = 0; col < Cols(); col++ )
420  if(IsZero(Get (row, col)) ){
421  out << ".";
422  }else{
423  out << "#";
424  }
425  out << "|";
426  }
427  out << "\n";
428  }
429 }
430 
431 template<>
432 void TPZMatrix<std::complex<double> >::Print(const char *name, std::ostream& out, const MatrixOutputFormat form) const {
433  //std::cout << __PRETTY_FUNCTION__ << " please implement me!\n";
434  //DebugStop();
435  if(form == EFormatted) {
436  out << "Writing matrix '";
437  if(name) out << name;
438  out << "' (" << Rows() << " x " << Cols() << "):\n";
439 
440  for ( int64_t row = 0; row < Rows(); row++) {
441  out << "\t";
442  for ( int64_t col = 0; col < Cols(); col++ ) {
443  out << Get( row, col).real() << " " << Get(row, col).imag() << " ";
444  }
445  out << "\n";
446  }
447  out << "\n";
448  } else if (form == EInputFormat) {
449  out << Rows() << " " << Cols() << endl;
450  for ( int64_t row = 0; row < Rows(); row++) {
451  for ( int64_t col = 0; col < Cols(); col++ ) {
452  std::complex<double> val = Get (row, col);
453  if(val != 0.) out << row << ' ' << col << ' ' << val.real() << " " << val.imag() << std::endl;
454  }
455  }
456  out << "-1 -1 0.\n";
457  } else if( form == EMathematicaInput)
458  {
459  char number[128];
460  out << name << "\n{ ";
461  for ( int64_t row = 0; row < Rows(); row++) {
462  out << "\n{ ";
463  for ( int64_t col = 0; col < Cols(); col++ ) {
464  std::complex<double> val = Get (row, col);
465 
466  sprintf(number, "%16.16lf + I %16.16lf", val.real(), val.imag());
467 
468 
469  out << number;
470  if(col < Cols()-1)
471  out << ", ";
472  if((col+1) % 6 == 0)out << std::endl;
473  }
474  out << " }";
475  if(row < Rows()-1)
476  out << ",";
477  }
478 
479  out << " };\n";
480 
481  }else if( form == EMatlabNonZeros)
482  {
483  out << name;
484  for ( int64_t row = 0; row < Rows(); row++) {
485  out << "\n|";
486  for ( int64_t col = 0; col < Cols(); col++ )
487  if(IsZero(Get (row, col)) ){
488  out << ".";
489  }else{
490  out << "#";
491  }
492  out << "|";
493  }
494  out << "\n";
495  }
496 
497 }
498 
499 template<>
500 void TPZMatrix<std::complex<long double> >::Print(const char *name, std::ostream& out,const MatrixOutputFormat form) const {
501  std::cout << __PRETTY_FUNCTION__ << " please implement me!\n";
502  DebugStop();
503 }
504 
505 template<class TVar>
506 void TPZMatrix<TVar>::AddKel(TPZFMatrix<TVar> &elmat, TPZVec<int64_t> &destinationindex) {
507 
508  int64_t nelem = elmat.Rows();
509  int64_t icoef,jcoef,ieq,jeq;
510  if(IsSimetric()) {
511  for(icoef=0; icoef<nelem; icoef++) {
512  ieq = destinationindex[icoef];
513  for(jcoef=icoef; jcoef<nelem; jcoef++) {
514  jeq = destinationindex[jcoef];
515  TVar prevval = GetVal(ieq,jeq);
516  prevval += elmat(icoef,jcoef);
517  PutVal(ieq,jeq,prevval);
518  }
519  }
520  } else {
521  for(icoef=0; icoef<nelem; icoef++) {
522  ieq = destinationindex[icoef];
523  for(jcoef=0; jcoef<nelem; jcoef++) {
524  jeq = destinationindex[jcoef];
525  TVar prevval = GetVal(ieq,jeq);
526  prevval += elmat(icoef,jcoef);
527  PutVal(ieq,jeq,prevval);
528  }
529  }
530  }
531 }
532 
533 template<class TVar>
534 void TPZMatrix<TVar>::AddKel(TPZFMatrix<TVar> &elmat, TPZVec<int64_t> &source, TPZVec<int64_t> &destinationindex) {
535 
536  int64_t nelem = source.NElements();
537  int64_t icoef,jcoef,ieq,jeq,ieqs,jeqs;
538  TVar prevval;
539  if(IsSimetric()) {
540  for(icoef=0; icoef<nelem; icoef++) {
541  ieq = destinationindex[icoef];
542  ieqs = source[icoef];
543  for(jcoef=icoef; jcoef<nelem; jcoef++) {
544  jeq = destinationindex[jcoef];
545  jeqs = source[jcoef];
546  prevval = GetVal(ieq,jeq);
547  prevval += elmat(ieqs,jeqs);
548  PutVal(ieq,jeq,prevval);
549  }
550  }
551  } else {
552  for(icoef=0; icoef<nelem; icoef++) {
553  ieq = destinationindex[icoef];
554  ieqs = source[icoef];
555  for(jcoef=0; jcoef<nelem; jcoef++) {
556  jeq = destinationindex[jcoef];
557  jeqs = source[jcoef];
558  prevval = GetVal(ieq,jeq);
559  prevval += elmat(ieqs,jeqs);
560  PutVal(ieq,jeq,prevval);
561  }
562  }
563  }
564 }
565 
566 
567 /***************/
568 /*** Put Sub ***/
569 //
570 // Put submatriz A como sub matriz a partir do ponto
571 // (sRow, sCol).
572 //
573 template<class TVar>
574 int TPZMatrix<TVar>::PutSub(const int64_t sRow,const int64_t sCol,const TPZFMatrix<TVar> &A ) {
575  if(sRow >= this->Rows()){
576  LOGPZ_ERROR(logger,"incompatible dimensions")
577  DebugStop();
578  }
579 
580  int64_t minRow = MIN( A.Rows(), Rows() - sRow );
581  int64_t minCol = MIN( A.Cols(), Cols() - sCol );
582 
583  int64_t row = sRow;
584  for ( int64_t r = 0; r < minRow; r++, row++ ) {
585  int64_t col = sCol;
586  for ( int64_t c = 0; c < minCol; c++, col++ ) {
587  PutVal( row, col, A.GetVal( r, c ) );
588  }
589  }
590  return( 1 );
591 }
592 
593 
594 
595 /***************/
596 /*** Get Sub ***/
597 //
598 // Le a sub-matriz de dimensoes (colSize, rowSize) para a matriz A.
599 // O inicio da sub-matriz e' dado por (sRow, sCol).
600 //
601 template<class TVar>
602 int TPZMatrix<TVar>::GetSub(const int64_t sRow,const int64_t sCol,const int64_t rowSize,
603  const int64_t colSize, TPZFMatrix<TVar> & A ) const {
604  if ( ((sRow + rowSize) > Rows()) || ((sCol + colSize) > Cols()) ) {
605  return( Error( "GetSub <t.he sub-matrix is too big>" ) );
606  }
607  A.Resize( rowSize, colSize );
608  int64_t row = sRow;
609  for ( int64_t r = 0; r < rowSize; r++, row++ ) {
610  int64_t col = sCol;
611  for ( int64_t c = 0; c < colSize; c++, col++ ) {
612  A.PutVal( r, c, GetVal( row, col ) );
613  }
614  }
615  return( 1 );
616 }
617 
618 
619 /***************/
620 /*** Add Sub ***/
621 //
622 // Adds a matriz A to a sub-matriz which starts at (row, col).
623 //
624 template<class TVar>
625 int TPZMatrix<TVar>::AddSub(const int64_t sRow,const int64_t sCol,const TPZFMatrix<TVar> &A ) {
626 
627  int64_t minRow = MIN( A.Rows(), Rows() - sRow );
628  int64_t minCol = MIN( A.Cols(), Cols() - sCol );
629  // REAL v;
630 
631  int64_t row = sRow;
632  for ( int64_t r = 0; r < minRow; r++, row++ ) {
633  int64_t col = sCol;
634  for ( int64_t c = 0; c < minCol; c++, col++ ) {
635  PutVal( row, col, GetVal( row, col ) + A.GetVal( r, c ) );
636  }
637  }
638  return( 1 );
639 }
640 
641 /***************/
642 /*** InsertSub ***/
643 //
644 // Inserts a submatriz with the current matrix without changing the dimensions
645 //
646 template<class TVar>
647 int TPZMatrix<TVar>::InsertSub(const int64_t sRow,const int64_t sCol,const int64_t rowSize,
648  const int64_t colSize,const int64_t pRow,const int64_t pCol, TPZMatrix<TVar> *pA ) const {
649 
650 
651  if ( ((pRow + rowSize) > pA->Rows()) || ((pCol + colSize) > pA->Cols())) {
652  return( Error( "InsertSub <the sub-matrix is too big that target>" ) );
653  }
654 
655  int64_t NewRowSize = rowSize+pRow;
656  int64_t NewColSize = colSize+pCol;
657 
658 
659  int64_t row = sRow;
660  for ( int64_t r = pRow; r < NewRowSize; r++, row++ ) {
661  int64_t col = sCol;
662  for ( int64_t c = pCol ; c < NewColSize; c++, col++ ) {
663  pA->PutVal( r, c, GetVal( row, col ) );
664  }
665  }
666  return( 1 );
667 }
668 
669 
670 /***************/
671 /*** AddSub ***/
672 //
673 // Adds the submatrix to *pA at the given place
674 //
675 template<class TVar>
676 int TPZMatrix<TVar>::AddSub(const int64_t sRow, const int64_t sCol, const int64_t rowSize,
677  const int64_t colSize,const int64_t pRow,const int64_t pCol, TPZMatrix<TVar> *pA ) const {
678  if ( ((pRow + rowSize) > pA->Rows()) || ((pCol + colSize) > pA->Cols())) {
679  Error( "AddSub <the sub-matrix is too big that target>" );
680  }
681  int64_t NewRowSize = rowSize+pRow;
682  int64_t NewColSize = colSize+pCol;
683 
684  int64_t row = sRow;
685  for ( int64_t r = pRow; r < NewRowSize; r++, row++ ) {
686  int64_t col = sCol;
687  for ( int64_t c = pCol ; c < NewColSize; c++, col++ ) {
688  pA->PutVal( r, c, GetVal( row, col )+pA->GetVal(r,c));
689  }
690  }
691  return( 1 );
692 }
693 
694 /*****************/
695 /*** Transpose ***/
696 template<class TVar>
697 void TPZMatrix<TVar>::Transpose(TPZMatrix<TVar> *T) const {
698  T->Resize( Cols(), Rows() );
699 
700  for ( int64_t r = 0; r < Rows(); r++ ) {
701  for ( int64_t c = 0; c < Cols(); c++ ) {
702  T->PutVal( c, r, GetVal( r, c ) );
703  }
704  }
705 }
706 
707 /*************/
708 /*** Solve ***/
709 template<class TVar>
710 int TPZMatrix<TVar>::SolveDirect( TPZFMatrix<TVar> &B , DecomposeType dt, std::list<int64_t> &singular) {
711 
712  switch ( dt ) {
713  case ELU:
714  return( Solve_LU( &B ,singular) );
715  case ECholesky:
716  return( Solve_Cholesky( &B , singular) );
717  case ELDLt:
718  return( Solve_LDLt( &B, singular ) );
719  default:
720  Error( "Solve < Unknow decomposition type >" );
721  break;
722  }
723  return ( 0 );
724 }
725 template<class TVar>
726 int TPZMatrix<TVar>::SolveDirect( TPZFMatrix<TVar> &B , DecomposeType dt) {
727 
728  switch ( dt ) {
729  case ELU:
730  return( Solve_LU( &B) );
731  case ECholesky:
732  return( Solve_Cholesky( &B ) );
733  case ELDLt:
734  return( Solve_LDLt( &B ) );
735  default:
736  Error( "Solve < Unknow decomposition type >" );
737  break;
738  }
739  return ( 0 );
740 }
741 
742 template<class TVar>
743 void TPZMatrix<TVar>::SolveJacobi(int64_t &numiterations,const TPZFMatrix<TVar> &F, TPZFMatrix<TVar> &result,
744  TPZFMatrix<TVar> *residual, TPZFMatrix<TVar> &scratch, REAL &tol,const int FromCurrent) {
745 
746 
747  if(FromCurrent) {
748  Residual(result,F,scratch);
749  } else {
750  scratch = F;
751  result.Zero();
752  }
753  REAL res;
754  res = TPZExtractVal::val(Norm(scratch));
755  int64_t r = Dim();
756  int64_t c = F.Cols();
757  for(int64_t it=0; it<numiterations && (fabs(res)) > tol; it++) {
758  for(int64_t ic=0; ic<c; ic++) {
759  for(int64_t i=0; i<r; i++) {
760  result(i,ic) += (scratch)(i,ic)/GetVal(i,i);
761  }
762  }
763  Residual(result,F,scratch);
764  res = TPZExtractVal::val(Norm(scratch));
765  }
766  if(residual) *residual = scratch;
767 }
768 
769 #ifdef _AUTODIFF
770 template<>
771 void TPZMatrix<TFad<6,REAL> >::SolveSOR(int64_t & numiterations, const TPZFMatrix<TFad<6,REAL> > &F,
772  TPZFMatrix<TFad<6,REAL> > &result, TPZFMatrix<TFad<6,REAL> > *residual, TPZFMatrix<TFad<6,REAL> > &/*scratch*/, const REAL overrelax,
773  REAL &tol,const int FromCurrent,const int direction) {
774  DebugStop();
775 }
776 template<>
777 void TPZMatrix<Fad<REAL> >::SolveSOR(int64_t & numiterations, const TPZFMatrix<Fad<REAL> > &F,
778  TPZFMatrix<Fad<REAL> > &result, TPZFMatrix<Fad<REAL> > *residual, TPZFMatrix<Fad<REAL> > &/*scratch*/, const REAL overrelax,
779  REAL &tol,const int FromCurrent,const int direction) {
780  DebugStop();
781 }
782 #endif
783 
784 template<class TVar>
785 void TPZMatrix<TVar>::SolveSOR(int64_t & numiterations, const TPZFMatrix<TVar> &F,
786  TPZFMatrix<TVar> &result, TPZFMatrix<TVar> *residual, TPZFMatrix<TVar> &/*scratch*/, const REAL overrelax,
787  REAL &tol,const int FromCurrent,const int direction) {
788 
789  if(residual == &F) {
790  cout << "TPZMatrix::SolveSOR called with residual and F equal, no solution\n";
791  return;
792  }
793  TVar res = (TVar)2*(TVar)tol+(TVar)1.;
794  if(residual) res = Norm(*residual);
795  if(!FromCurrent) {
796  result.Zero();
797  }
798  int64_t r = Dim();
799  int64_t c = F.Cols();
800  int64_t i,ifirst = 0, ilast = r, iinc = 1;
801  int64_t it;
802  if(direction == -1) {
803  ifirst = r-1; //misael
804  ilast = -1;
805  iinc = -1;
806  }
807  TVar eqres;
808  for(it=0; it<numiterations && (REAL)(fabs(res)) > fabs(tol); it++) {
809  res = 0.;
810  for(int64_t ic=0; ic<c; ic++) {
811  for(i=ifirst; i!=ilast; i+= iinc) {
812  eqres = F.GetVal(i,ic);
813  for(int64_t j=0; j<r; j++) {
814  eqres -= GetVal(i,j)*result(j,ic);
815  }
816  res += eqres*eqres;
817  result(i,ic) += (TVar)overrelax*eqres/GetVal(i,i);
818  }
819  }
820  res = sqrt(res);
821  }
822  if(residual) Residual(result,F,*residual);
823  numiterations = it;
824  tol = fabs(res);
825 }
826 
827 template <class TVar>
828 void TPZMatrix<TVar>::SolveSSOR(int64_t &numiterations, const TPZFMatrix<TVar> &F,
829  TPZFMatrix<TVar> &result, TPZFMatrix<TVar> *residual, TPZFMatrix<TVar> &scratch, const REAL overrelax,
830  REAL &tol,const int FromCurrent) {
831  REAL res = tol*2.+1.;
832  int64_t i, one = 1;
833  int fromcurrent = FromCurrent;
834  for(i=0; i<numiterations && fabs(res) > fabs(tol); i++) {
835  one = 1;
836  res = tol;
837  SolveSOR(one,F,result,residual,scratch,overrelax,res,fromcurrent,1);
838  one = 1;
839  res = tol;
840  fromcurrent = 1;
841  SolveSOR(one,F,result,residual,scratch,overrelax,res,fromcurrent,-1);
842  // cout << "SSOR iter = " << i << " res = " << res << endl;
843  }
844  numiterations = i;
845  tol = res;
846 }
847 
848 #include "cg.h"
849 template <class TVar>
850 void TPZMatrix<TVar>::SolveCG(int64_t &numiterations, TPZSolver<TVar> &preconditioner,
851  const TPZFMatrix<TVar> &F, TPZFMatrix<TVar> &result,
852  TPZFMatrix<TVar> *residual, REAL &tol, const int FromCurrent) {
853  CG(*this, result, F, preconditioner, residual, numiterations, tol, FromCurrent);
854 }
855 
856 #ifdef _AUTODIFF
857 template <>
858 void TPZMatrix< TFad<6,REAL> >::SolveCG(int64_t &numiterations, TPZSolver< TFad<6,REAL> > &preconditioner,
859  const TPZFMatrix< TFad<6,REAL> > &F, TPZFMatrix< TFad<6,REAL> > &result,
860  TPZFMatrix< TFad<6,REAL> > *residual, REAL &tol, const int FromCurrent) {
861  DebugStop(); // Does not work with complex numbers. To be implemented in the future.
862 }
863 template <>
864 void TPZMatrix< Fad<REAL> >::SolveCG(int64_t &numiterations, TPZSolver< Fad<REAL> > &preconditioner,
865  const TPZFMatrix< Fad<REAL> > &F, TPZFMatrix< Fad<REAL> > &result,
866  TPZFMatrix< Fad<REAL> > *residual, REAL &tol, const int FromCurrent) {
867  DebugStop(); // Does not work with complex numbers. To be implemented in the future.
868 }
869 
870 #endif
871 
872 template <>
873 void TPZMatrix< std::complex<float> >::SolveCG(int64_t &numiterations, TPZSolver< std::complex<float> > &preconditioner,
874  const TPZFMatrix< std::complex<float> > &F, TPZFMatrix< std::complex<float> > &result,
875  TPZFMatrix< std::complex<float> > *residual, REAL &tol, const int FromCurrent) {
876  DebugStop(); // Does not work with complex numbers. To be implemented in the future.
877 }
878 
879 template <>
880 void TPZMatrix< std::complex<double> >::SolveCG(int64_t &numiterations, TPZSolver< std::complex<double> > &preconditioner,
881  const TPZFMatrix< std::complex<double> > &F, TPZFMatrix< std::complex<double> > &result,
882  TPZFMatrix< std::complex<double> > *residual, REAL &tol, const int FromCurrent) {
883  DebugStop(); // Does not work with complex numbers. To be implemented in the future.
884 }
885 
886 template <>
887 void TPZMatrix< std::complex<long double> >::SolveCG(int64_t &numiterations, TPZSolver< std::complex<long double> > &preconditioner,
888  const TPZFMatrix< std::complex<long double> > &F, TPZFMatrix< std::complex<long double> > &result,
889  TPZFMatrix< std::complex<long double> > *residual, REAL &tol, const int FromCurrent) {
890  DebugStop(); // Does not work with complex numbers. To be implemented in the future.
891 }
892 
893 #include "gmres.h"
894 #include "pzstepsolver.h"
895 template <class TVar>
896 void TPZMatrix<TVar>::SolveGMRES(int64_t &numiterations, TPZSolver<TVar> &preconditioner,
897  TPZFMatrix<TVar> &H, int &numvectors,
898  const TPZFMatrix<TVar> &F, TPZFMatrix<TVar> &result,
899  TPZFMatrix<TVar> *residual, REAL &tol,const int FromCurrent)
900 {
901  if(F.Cols() > 1)
902  {
903  int64_t locnumiter = numiterations;
904  TVar loctol = tol;
905  int64_t nrow = F.Rows();
906  int64_t ncol = F.Cols();
907  int64_t col;
908  //preconditioner.Solve(F, result);
909  for (col=0; col<ncol; col++) {
910 // std::cout << "Column " << col << std::endl;
911  numiterations = locnumiter;
912  tol = TPZExtractVal::val(loctol);
913  TPZFMatrix<TVar> FCol(nrow,1);
914  memcpy((void *)(&FCol(0,0)), (void *)(&F.GetVal(0,col)), nrow*sizeof(REAL));
915  TPZFMatrix<TVar> resultCol(nrow,1,&result(0,col),nrow);
916  GMRES(*this,resultCol,FCol,preconditioner,H,numvectors,numiterations,tol,residual,FromCurrent);
917  }
918  }
919  else
920  {
921  GMRES(*this,result,F,preconditioner,H,numvectors,numiterations,tol,residual,FromCurrent);
922  }
923 }
924 
925 #ifdef _AUTODIFF
926 template <>
927 void TPZMatrix< TFad<6,REAL> >::SolveGMRES(int64_t &numiterations, TPZSolver< TFad<6,REAL> > &preconditioner,
928  TPZFMatrix< TFad<6,REAL> > &H, int &numvectors,
929  const TPZFMatrix< TFad<6,REAL> > &F, TPZFMatrix< TFad<6,REAL> > &result,
930  TPZFMatrix< TFad<6,REAL> > *residual, REAL &tol,const int FromCurrent)
931 {
932  DebugStop(); // Does not work with complex numbers. To be implemented in the future.
933 }
934 template <>
935 void TPZMatrix< Fad<REAL> >::SolveGMRES(int64_t &numiterations, TPZSolver< Fad<REAL> > &preconditioner,
936  TPZFMatrix< Fad<REAL> > &H, int &numvectors,
937  const TPZFMatrix< Fad<REAL> > &F, TPZFMatrix< Fad<REAL> > &result,
938  TPZFMatrix< Fad<REAL> > *residual, REAL &tol,const int FromCurrent)
939 {
940  DebugStop(); // Does not work with complex numbers. To be implemented in the future.
941 }
942 #endif
943 
944 template <>
945 void TPZMatrix< std::complex<float> >::SolveGMRES(int64_t &numiterations, TPZSolver< std::complex<float> > &preconditioner,
946  TPZFMatrix< std::complex<float> > &H, int &numvectors,
947  const TPZFMatrix< std::complex<float> > &F, TPZFMatrix< std::complex<float> > &result,
948  TPZFMatrix< std::complex<float> > *residual, REAL &tol,const int FromCurrent)
949 {
950  DebugStop(); // Does not work with complex numbers. To be implemented in the future.
951 }
952 
953 template <>
954 void TPZMatrix< std::complex<double> >::SolveGMRES(int64_t &numiterations, TPZSolver< std::complex<double> > &preconditioner,
955  TPZFMatrix< std::complex<double> > &H, int &numvectors,
956  const TPZFMatrix< std::complex<double> > &F, TPZFMatrix< std::complex<double> > &result,
957  TPZFMatrix< std::complex<double> > *residual, REAL &tol,const int FromCurrent)
958 {
959  DebugStop(); // Does not work with complex numbers. To be implemented in the future.
960 }
961 
962 template <>
963 void TPZMatrix< std::complex<long double> >::SolveGMRES(int64_t &numiterations, TPZSolver< std::complex<long double> > &preconditioner,
964  TPZFMatrix< std::complex<long double> > &H, int &numvectors,
965  const TPZFMatrix< std::complex<long double> > &F, TPZFMatrix< std::complex<long double> > &result,
966  TPZFMatrix< std::complex<long double> > *residual, REAL &tol,const int FromCurrent)
967 {
968  DebugStop(); // Does not work with complex numbers. To be implemented in the future.
969 }
970 
971 
972 #include "bicg.h"
973 template<class TVar>
974 void TPZMatrix<TVar>::SolveBICG(int64_t &numiterations, TPZSolver<TVar> &preconditioner,
975  const TPZFMatrix<TVar> &F,
976  TPZFMatrix<TVar> &result,
977  REAL &tol) {
978  BiCG (*this, result,F,preconditioner,numiterations,tol);
979 }
980 
981 #ifdef _AUTODIFF
982 template<>
983 void TPZMatrix< TFad<6,REAL> >::SolveBICG(int64_t &numiterations, TPZSolver< TFad<6,REAL> > &preconditioner,
984  const TPZFMatrix< TFad<6,REAL> > &F,
985  TPZFMatrix< TFad<6,REAL> > &result,
986  REAL &tol) {
987  DebugStop(); // Does not work with complex numbers. To be implemented in the future.
988 }
989 template<>
990 void TPZMatrix< Fad<REAL> >::SolveBICG(int64_t &numiterations, TPZSolver< Fad<REAL> > &preconditioner,
991  const TPZFMatrix< Fad<REAL> > &F,
992  TPZFMatrix< Fad<REAL> > &result,
993  REAL &tol) {
994  DebugStop(); // Does not work with complex numbers. To be implemented in the future.
995 }
996 #endif
997 
998 template<>
999 void TPZMatrix< std::complex<float> >::SolveBICG(int64_t &numiterations, TPZSolver< std::complex<float> > &preconditioner,
1000  const TPZFMatrix< std::complex<float> > &F,
1001  TPZFMatrix< std::complex<float> > &result,
1002  REAL &tol) {
1003  DebugStop(); // Does not work with complex numbers. To be implemented in the future.
1004 }
1005 
1006 template<>
1007 void TPZMatrix< std::complex<double> >::SolveBICG(int64_t &numiterations, TPZSolver< std::complex<double> > &preconditioner,
1008  const TPZFMatrix< std::complex<double> > &F,
1009  TPZFMatrix< std::complex<double> > &result,
1010  REAL &tol) {
1011  DebugStop(); // Does not work with complex numbers. To be implemented in the future.
1012 }
1013 
1014 template<>
1015 void TPZMatrix< std::complex<long double> >::SolveBICG(int64_t &numiterations, TPZSolver< std::complex<long double> > &preconditioner,
1016  const TPZFMatrix< std::complex<long double> > &F,
1017  TPZFMatrix< std::complex<long double> > &result,
1018  REAL &tol) {
1019  DebugStop(); // Does not work with complex numbers. To be implemented in the future.
1020 }
1021 
1022 #include "bicgstab.h"
1023 template <class TVar>
1024 void TPZMatrix<TVar>::SolveBICGStab(int64_t &numiterations, TPZSolver<TVar> &preconditioner,
1025  const TPZFMatrix<TVar> &F, TPZFMatrix<TVar> &result,
1026  TPZFMatrix<TVar> *residual, REAL &tol,const int FromCurrent) {
1027  BiCGSTAB(*this,result,F,preconditioner,numiterations,tol,residual,FromCurrent);
1028 }
1029 
1030 #ifdef _AUTODIFF
1031 template <>
1032 void TPZMatrix< TFad<6,REAL> >::SolveBICGStab(int64_t &numiterations, TPZSolver< TFad<6,REAL> > &preconditioner,
1033  const TPZFMatrix< TFad<6,REAL> > &F, TPZFMatrix< TFad<6,REAL> > &result,
1034  TPZFMatrix< TFad<6,REAL> > *residual, REAL &tol,const int FromCurrent) {
1035  DebugStop(); // Does not work with complex numbers. To be implemented in the future.
1036 }
1037 template <>
1038 void TPZMatrix< Fad<REAL> >::SolveBICGStab(int64_t &numiterations, TPZSolver< Fad<REAL> > &preconditioner,
1039  const TPZFMatrix< Fad<REAL> > &F, TPZFMatrix< Fad<REAL> > &result,
1040  TPZFMatrix< Fad<REAL> > *residual, REAL &tol,const int FromCurrent) {
1041  DebugStop(); // Does not work with complex numbers. To be implemented in the future.
1042 }
1043 #endif
1044 
1045 template <>
1046 void TPZMatrix< std::complex<float> >::SolveBICGStab(int64_t &numiterations, TPZSolver< std::complex<float> > &preconditioner,
1047  const TPZFMatrix< std::complex<float> > &F, TPZFMatrix< std::complex<float> > &result,
1048  TPZFMatrix< std::complex<float> > *residual, REAL &tol,const int FromCurrent) {
1049  DebugStop(); // Does not work with complex numbers. To be implemented in the future.
1050 }
1051 
1052 template <>
1053 void TPZMatrix< std::complex<double> >::SolveBICGStab(int64_t &numiterations, TPZSolver< std::complex<double> > &preconditioner,
1054  const TPZFMatrix< std::complex<double> > &F, TPZFMatrix< std::complex<double> > &result,
1055  TPZFMatrix< std::complex<double> > *residual, REAL &tol,const int FromCurrent) {
1056  DebugStop(); // Does not work with complex numbers. To be implemented in the future.
1057 }
1058 
1059 template <>
1060 void TPZMatrix< std::complex<long double> >::SolveBICGStab(int64_t &numiterations, TPZSolver< std::complex<long double> > &preconditioner,
1061  const TPZFMatrix< std::complex<long double> > &F, TPZFMatrix< std::complex<long double> > &result,
1062  TPZFMatrix< std::complex<long double> > *residual, REAL &tol,const int FromCurrent) {
1063  DebugStop(); // Does not work with complex numbers. To be implemented in the future.
1064 }
1065 
1066 
1067 
1068 #include "ir.h"
1069 template <class TVar>
1070 void TPZMatrix<TVar>::SolveIR(int64_t &numiterations, TPZSolver<TVar> &preconditioner,
1071  const TPZFMatrix<TVar> &F, TPZFMatrix<TVar> &result,
1072  TPZFMatrix<TVar> *residual, REAL &tol,
1073  const int FromCurrent) {
1074  IR(*this,result,F,preconditioner,residual,numiterations,tol,FromCurrent);
1075 }
1076 
1077 #ifdef _AUTODIFF
1078 template <>
1079 void TPZMatrix< TFad<6,REAL> >::SolveIR(int64_t &numiterations, TPZSolver< TFad<6,REAL> > &preconditioner,
1080  const TPZFMatrix< TFad<6,REAL> > &F, TPZFMatrix< TFad<6,REAL> > &result,
1081  TPZFMatrix< TFad<6,REAL> > *residual, REAL &tol,
1082  const int FromCurrent) {
1083  DebugStop(); // Does not work with complex numbers. To be implemented in the future.
1084 }
1085 template <>
1086 void TPZMatrix< Fad<REAL> >::SolveIR(int64_t &numiterations, TPZSolver< Fad<REAL> > &preconditioner,
1087  const TPZFMatrix< Fad<REAL> > &F, TPZFMatrix< Fad<REAL> > &result,
1088  TPZFMatrix< Fad<REAL> > *residual, REAL &tol,
1089  const int FromCurrent) {
1090  DebugStop(); // Does not work with complex numbers. To be implemented in the future.
1091 }
1092 #endif
1093 
1094 template <>
1095 void TPZMatrix< std::complex<float> >::SolveIR(int64_t &numiterations, TPZSolver< std::complex<float> > &preconditioner,
1096  const TPZFMatrix< std::complex<float> > &F, TPZFMatrix< std::complex<float> > &result,
1097  TPZFMatrix< std::complex<float> > *residual, REAL &tol,
1098  const int FromCurrent) {
1099  DebugStop(); // Does not work with complex numbers. To be implemented in the future.
1100 }
1101 
1102 template <>
1103 void TPZMatrix< std::complex<double> >::SolveIR(int64_t &numiterations, TPZSolver< std::complex<double> > &preconditioner,
1104  const TPZFMatrix< std::complex<double> > &F, TPZFMatrix< std::complex<double> > &result,
1105  TPZFMatrix< std::complex<double> > *residual, REAL &tol,
1106  const int FromCurrent) {
1107  DebugStop(); // Does not work with complex numbers. To be implemented in the future.
1108 }
1109 
1110 template <>
1111 void TPZMatrix< std::complex<long double> >::SolveIR(int64_t &numiterations, TPZSolver< std::complex<long double> > &preconditioner,
1112  const TPZFMatrix< std::complex<long double> > &F, TPZFMatrix< std::complex<long double> > &result,
1113  TPZFMatrix< std::complex<long double> > *residual, REAL &tol,
1114  const int FromCurrent) {
1115  DebugStop(); // Does not work with complex numbers. To be implemented in the future.
1116 }
1117 
1118 
1119 /*****************/
1120 /*** Decompose_LU ***/
1121 template <class TVar>
1122 int TPZMatrix<TVar>::Decompose_LU(std::list<int64_t> &singular) {
1123  return Decompose_LU();
1124 }
1125 template <class TVar>
1126 int TPZMatrix<TVar>::Decompose_LU() {
1127 
1128  if ( fDecomposed && fDecomposed != ELU) Error( "Decompose_LU <Matrix already Decomposed with other scheme>" );
1129  if (fDecomposed) return 1;
1130 
1131  TVar nn, pivot;
1132  int64_t min = ( Cols() < (Rows()) ) ? Cols() : Rows();
1133 
1134  for ( int64_t k = 0; k < min ; k++ ) {
1135  if (IsZero( pivot = GetVal(k, k))) Error( "Decompose_LU <matrix is singular>" );
1136  for ( int64_t i = k+1; i < Rows(); i++ ) {
1137  nn = GetVal( i, k ) / pivot;
1138  PutVal( i, k, nn );
1139  for ( int64_t j = k+1; j < Cols(); j++ ) PutVal(i,j,GetVal(i,j)-nn*GetVal(k,j));
1140  }
1141  }
1142  fDecomposed=ELU;
1143  return 1;
1144 }
1145 
1146 
1147 
1148 /****************/
1149 /*** Substitution ***/
1150 template <class TVar>
1151 int TPZMatrix<TVar>::Substitution( TPZFMatrix<TVar> *B ) const{
1152 
1153  int64_t rowb = B->Rows();
1154  int64_t colb = B->Cols();
1155  if ( rowb != Rows() )
1156  Error( "SubstitutionLU <incompatible dimensions>" );
1157  int64_t i;
1158  for ( i = 0; i < rowb; i++ ) {
1159  for ( int64_t col = 0; col < colb; col++ ) {
1160  for ( int64_t j = 0; j < i; j++ ) {
1161  B->PutVal( i, col, B->GetVal(i, col)-GetVal(i, j) * B->GetVal(j, col) );
1162  }
1163  }
1164  }
1165  for (int64_t col=0; col<colb; col++) {
1166  for ( i = rowb-1; i >= 0; i-- ) {
1167  for ( int64_t j = i+1; j < rowb ; j++ ) {
1168  B->PutVal( i, col, B->GetVal(i, col) -
1169  GetVal(i, j) * B->GetVal(j, col) );
1170  }
1171  if ( IsZero( GetVal(i, i) ) ) {
1172  Error( "BackSub( SubstitutionLU ) <Matrix is singular" );
1173  }
1174  B->PutVal( i, col, B->GetVal( i, col) / GetVal(i, i) );
1175  }
1176  }
1177  return( 1 );
1178 }
1179 template <class TVar>
1180 int TPZMatrix<TVar>::Decompose_LDLt(std::list<int64_t> &singular) {
1181  return Decompose_LDLt();
1182 }
1183 template <class TVar>
1184 int TPZMatrix<TVar>::Decompose_LDLt() {
1185 
1186  if ( fDecomposed && fDecomposed != ELDLt) {
1187  Error( "Decompose_LDLt <Matrix already Decomposed with other scheme> " );
1188  } else if(fDecomposed ) {
1189  return ELDLt;
1190  }
1191  if ( Rows()!=Cols() ) Error( "Decompose_LDLt <Matrix must be square>" );
1192 
1193  int64_t j,k,l,dim=Rows();
1194 
1195  for ( j = 0; j < dim; j++ ) {
1196  for ( k=0; k<j; k++) {
1197  PutVal( j,j,GetVal(j,j) - GetVal(k,k)*GetVal(k,j)*GetVal(k,j) );
1198  }
1199  for ( k=0; k<j; k++) {
1200  for( l=j+1; l<dim;l++) {
1201  PutVal(l,j, GetVal(l,j)-GetVal(k,k)*GetVal(j,k)*GetVal(l,k) );
1202  PutVal(j,l,GetVal(l,j) );
1203  }
1204  }
1205  TVar tmp = GetVal(j,j);
1206  if ( IsZero(tmp) ) Error( "Decompose_LDLt <Zero on diagonal>" );
1207  for( l=j+1; l<dim;l++) {
1208  PutVal(l,j, GetVal(l,j)/GetVal(j,j) ) ;
1209  PutVal(j,l, GetVal(l,j) );
1210  }
1211  }
1212  fDecomposed = ELDLt;
1213  fDefPositive = 0;
1214  return( 1 );
1215 }
1216 template <class TVar>
1217 int TPZMatrix<TVar>::Decompose_Cholesky(std::list<int64_t> &singular) {
1218  if ( fDecomposed && fDecomposed != ECholesky) Error( "Decompose_Cholesky <Matrix already Decomposed>" );
1219  if ( fDecomposed ) return ECholesky;
1220  if ( Rows()!=Cols() ) Error( "Decompose_Cholesky <Matrix must be square>" );
1221  //return 0;
1222 
1223  int64_t dim=Dim();
1224  for (int64_t i=0 ; i<dim; i++) {
1225  for(int64_t k=0; k<i; k++) { //elementos da diagonal
1226  PutVal( i,i,GetVal(i,i)-GetVal(i,k)*GetVal(i,k) );
1227  }
1228  if((fabs(GetVal(i,i))) <= fabs((TVar)1.e-12))
1229  {
1230  singular.push_back(i);
1231  PutVal(i,i,1.);
1232  }
1233  TVar tmp = sqrt(GetVal(i,i));
1234  PutVal( i,i,tmp );
1235  for (int64_t j=i+1;j<dim; j++) { //elementos fora da diagonal
1236  for(int64_t k=0; k<i; k++) {
1237  PutVal( i,j,GetVal(i,j)-GetVal(i,k)*GetVal(k,j) );
1238  }
1239  TVar tmp2 = GetVal(i,i);
1240  if ( IsZero(tmp2) ) {
1241  Error( "Decompose_Cholesky <Zero on diagonal>" );
1242  }
1243  PutVal(i,j,GetVal(i,j)/GetVal(i,i) );
1244  PutVal(j,i,GetVal(i,j));
1245 
1246  }
1247  }
1248  fDecomposed = ECholesky;
1249  return ECholesky;
1250 }
1251 template <class TVar>
1252 int TPZMatrix<TVar>::Decompose_Cholesky() {
1253  if ( fDecomposed && fDecomposed != ECholesky) Error( "Decompose_Cholesky <Matrix already Decomposed>" );
1254  if ( fDecomposed ) return ECholesky;
1255  if ( Rows()!=Cols() ) Error( "Decompose_Cholesky <Matrix must be square>" );
1256  //return 0;
1257 
1258  int64_t dim=Dim();
1259  for (int64_t i=0 ; i<dim; i++) {
1260  for(int64_t k=0; k<i; k++) { //elementos da diagonal
1261  PutVal( i,i,GetVal(i,i)-GetVal(i,k)*GetVal(i,k) );
1262  }
1263  TVar tmp = sqrt(GetVal(i,i));
1264  PutVal( i,i,tmp );
1265  for (int64_t j=i+1;j<dim; j++) { //elementos fora da diagonal
1266  for(int64_t k=0; k<i; k++) {
1267  PutVal( i,j,GetVal(i,j)-GetVal(i,k)*GetVal(k,j) );
1268  }
1269  TVar tmp2 = GetVal(i,i);
1270  if ( IsZero(tmp2) ) {
1271  Error( "Decompose_Cholesky <Zero on diagonal>" );
1272  }
1273  PutVal(i,j,GetVal(i,j)/GetVal(i,i) );
1274  PutVal(j,i,GetVal(i,j));
1275 
1276  }
1277  }
1278  fDecomposed = ECholesky;
1279  return ECholesky;
1280 
1281 }
1282 
1283 template <class TVar>
1284 int TPZMatrix<TVar>::Subst_Forward( TPZFMatrix<TVar> *B ) const {
1285  if ( (B->Rows() != Dim()) || !fDecomposed || fDecomposed != ECholesky)
1286  return( 0 );
1287  for ( int64_t r = 0; r < Dim(); r++ ) {
1288  TVar pivot = GetVal( r, r );
1289  for ( int64_t c = 0; c < B->Cols(); c++ ) {
1290  // Faz sum = SOMA( A[r,i] * B[i,c] ); i = 0, ..., r-1.
1291  //
1292  TVar sum = 0.0;
1293  for ( int64_t i = 0; i < r; i++ ) sum += GetVal(r, i) * B->GetVal(i, c);
1294 
1295  // Faz B[r,c] = (B[r,c] - sum) / A[r,r].
1296  //
1297  B->PutVal( r, c, (B->GetVal(r, c) - sum) / pivot );
1298  }
1299  }
1300  return( 1 );
1301 }
1302 
1303 /**********************/
1304 /*** Subst Backward ***/
1305 //
1306 // Faz Ax = b, onde A(NxN) e' triangular superior.
1307 //
1308 template <class TVar>
1309 int TPZMatrix<TVar>::Subst_Backward( TPZFMatrix<TVar> *B ) const {
1310 
1311  if ( (B->Rows() != Dim()) || !fDecomposed || fDecomposed != ECholesky) return( 0 );
1312  for ( int64_t r = Dim()-1; r >= 0; r-- ) {
1313  TVar pivot = GetVal( r, r );
1314  for ( int64_t c = 0; c < B->Cols(); c++ ) {
1315  // Faz sum = SOMA( A[r,i] * B[i,c] ); i = N, ..., r+1.
1316  //
1317  TVar sum = 0.0;
1318  for ( int64_t i = Dim()-1; i > r; i-- ) sum += GetVal(r, i) * B->GetVal(i, c);
1319  // Faz B[r,c] = (B[r,c] - sum) / A[r,r].
1320  //
1321  B->PutVal( r, c, (B->GetVal(r, c) - sum) / pivot );
1322  }
1323  }
1324  return( 1 );
1325 }
1326 
1327 /***********************/
1328 /*** Subst L Forward ***/
1329 //
1330 // Faz Ax = b, onde A e' triangular inferior e A(i,i) = 1.
1331 //
1332 template <class TVar>
1333 int TPZMatrix<TVar>::Subst_LForward( TPZFMatrix<TVar> *B ) const {
1334  if ( (B->Rows() != Dim()) || !fDecomposed || fDecomposed != ELDLt) {
1335  Error("TPZMatrix::Subst_LForward incompatible dimensions\n");
1336  }
1337  for ( int64_t r = 0; r < Dim(); r++ ) {
1338  for ( int64_t c = 0; c < B->Cols(); c++ ) {
1339  // Faz sum = SOMA( A[r,i] * B[i,c] ); i = 0, ..., r-1.
1340  //
1341  TVar sum = 0.0;
1342  for ( int64_t i = 0; i < r; i++ ) sum += GetVal(r, i) * B->GetVal(i, c);
1343 
1344  // Faz B[r,c] = (B[r,c] - sum) / A[r,r].
1345  //
1346  B->PutVal( r, c, B->GetVal(r, c) - sum );
1347  }
1348  }
1349  return( 1 );
1350 }
1351 
1352 /************************/
1353 /*** Subst L Backward ***/
1354 //
1355 // Faz Ax = b, onde A e' triangular superior e A(i,i) = 1.
1356 //
1357 template <class TVar>
1358 int TPZMatrix<TVar>::Subst_LBackward( TPZFMatrix<TVar> *B ) const {
1359  if ( (B->Rows() != Dim()) || !fDecomposed || fDecomposed != ELDLt){
1360  Error("TPZMatrix::Subst_LBackward incompatible dimensions \n");
1361  }
1362 
1363  for ( int64_t r = Dim()-1; r >= 0; r-- ) {
1364  for ( int64_t c = 0; c < B->Cols(); c++ ) {
1365  // Faz sum = SOMA( A[r,i] * B[i,c] ); i = N, ..., r+1.
1366  //
1367  TVar sum = 0.0;
1368  for ( int64_t i = Dim()-1; i > r; i-- ) sum += GetVal(r, i) * B->GetVal(i, c);
1369 
1370  // Faz B[r,c] = B[r,c] - sum.
1371  //
1372  B->PutVal( r, c, B->GetVal(r, c) - sum );
1373  }
1374  }
1375  return( 1 );
1376 }
1377 
1378 /******************/
1379 /*** Subst Diag ***/
1380 //
1381 // Faz Ax = b, sendo que A e' assumida ser uma matriz diagonal.
1382 //
1383 template <class TVar>
1384 int TPZMatrix<TVar>::Subst_Diag( TPZFMatrix<TVar> *B ) const {
1385  if ( (B->Rows() != Dim())) {
1386  Error("TPZMatrix::Subst_Diag incompatible dimensions\n");
1387  }
1388  for ( int64_t r = 0; r < Dim(); r++ ) {
1389  TVar pivot = GetVal( r, r );
1390  for ( int64_t c = 0; c < B->Cols(); c++ ) {
1391  B->PutVal( r, c, B->GetVal( r, c ) / pivot );
1392  }
1393  }
1394  return( 1 );
1395 }
1396 
1397 /************************** Private **************************/
1398 
1399 /*************/
1400 /*** Error ***/
1401 template <class TVar>
1402 int TPZMatrix<TVar>::Error(const char *msg ,const char *msg2) {
1403  ostringstream out;
1404  out << "TPZMatrix::" << msg;
1405  if(msg2) out << msg2;
1406  out << ".\n";
1407  LOGPZ_ERROR (logger, out.str().c_str());
1408  DebugStop();
1409  std::bad_exception myex;
1410  throw myex;
1411 }
1412 template <class TVar>
1413 void TPZMatrix<TVar>::Read( TPZStream &buf, void *context ){ //ok
1414  buf.Read(&fRow);
1415  buf.Read(&fCol);
1416  buf.Read(&fDecomposed);
1417  buf.Read(&fDefPositive);
1418 }
1419 template <class TVar>
1420 void TPZMatrix<TVar>::Write( TPZStream &buf, int withclassid ) const { //ok
1421  buf.Write(&fRow);
1422  buf.Write(&fCol);
1423  buf.Write(&fDecomposed);
1424  buf.Write(&fDefPositive);
1425 }
1426 
1428 
1432 template <class TVar>
1433 bool TPZMatrix<TVar>::Compare(TPZSavable *copy, bool override)
1434 {
1435  TPZMatrix<TVar> *copmat = dynamic_cast<TPZMatrix<TVar> *> (copy);
1436  if(!copmat) return false;
1437  bool result = true;
1438  if(fRow != copmat->fRow || fCol != copmat->fCol || fDecomposed != copmat->fDecomposed) result = false;
1439  if(!result)
1440  {
1441  std::stringstream sout;
1442  sout << __PRETTY_FUNCTION__ << " did not compare ";
1443  LOGPZ_ERROR(loggerCheck,sout.str())
1444  }
1445  if(override && !result)
1446  {
1447  this->operator=(*copmat);
1448  }
1449  return result;
1450 }
1451 
1453 
1457 template <class TVar>
1458 bool TPZMatrix<TVar>::Compare(TPZSavable *copy, bool override) const
1459 {
1460  TPZMatrix<TVar> *copmat = dynamic_cast<TPZMatrix<TVar> *> (copy);
1461  if(!copmat) return false;
1462  bool result = true;
1463  if(fRow != copmat->fRow || fCol != copmat->fCol || fDecomposed != copmat->fDecomposed) result = false;
1464  if(!result)
1465  {
1466  std::stringstream sout;
1467  sout << __PRETTY_FUNCTION__ << " did not compare ";
1468  LOGPZ_ERROR(loggerCheck,sout.str())
1469  }
1470  if(override && !result)
1471  {
1472  DebugStop();
1473  }
1474  return result;
1475 }
1476 
1477 template <class TVar>
1478 void TPZMatrix<TVar>::GetSub(const TPZVec<int64_t> &indices,TPZFMatrix<TVar> &block) const
1479 {
1480  int64_t nel = indices.NElements();
1481  int64_t iel,jel;
1482  for(iel=0; iel<nel; iel++)
1483  {
1484  for(jel=0; jel<nel; jel++)
1485  {
1486  block(iel,jel) = GetVal(indices[iel],indices[jel]);
1487  }
1488  }
1489 
1490 }
1491 
1492 #ifdef _AUTODIFF
1493 template<>
1494 int TPZMatrix<TFad<6,REAL> >::VerifySymmetry(REAL tol) const{
1495  DebugStop();
1496  return -1;
1497 }
1498 #endif
1499 
1500 template <class TVar>
1501 int TPZMatrix<TVar>::VerifySymmetry(REAL tol) const{
1502  int64_t nrows = this->Rows();
1503  int64_t ncols = this->Cols();
1504  if (nrows != ncols) return 0;
1505 
1506  for( int64_t i = 0; i < nrows; i++){
1507  for(int64_t j = 0; j <= i; j++){
1508  TVar exp = this->Get(i,j) - this->Get(j,i);
1509  if ( (REAL)(fabs( exp )) > tol ) {
1510  #ifdef STATE_COMPLEX
1511  cout << "Elemento: " << i << ", " << j << " -> " << fabs( exp ) << "/" <<
1512  this->Get(i,j) << endl;
1513  #else
1514  cout << "Elemento: " << i << ", " << j << " -> " << exp << "/" <<
1515  this->Get(i,j) << endl;
1516  #endif
1517  return 0;
1518  }
1519  }
1520  }
1521  return 1;
1522 }
1523 
1524 #ifdef _AUTODIFF
1525 template<>
1526 bool TPZMatrix<TFad<6,REAL> >::CompareValues(TPZMatrix<TFad<6,REAL> > &M, TFad<6,REAL> tol){
1527  DebugStop();
1528  return false;
1529 }
1530 #endif
1531 template <class TVar>
1532 bool TPZMatrix<TVar>::CompareValues(TPZMatrix<TVar> &M, TVar tol){
1533 
1534  int64_t nrows = this->Rows();
1535  int64_t ncols = this->Cols();
1536  if ( (M.Rows() != nrows) || (M.Cols() != ncols) ) return false;
1537 
1538  REAL diff;
1539  for( int64_t i = 0; i < nrows; i++)
1540  for( int64_t j = 0; j < ncols; j++){
1541  TVar exps = this->Get(i,j) - M.Get(i,j);
1542  diff = fabs( exps );
1543  if (diff > fabs(tol)) return false;
1544  }
1545 
1546  return true;
1547 }
1548 
1549 #ifdef _AUTODIFF
1550 template<>
1551 TFad<6,REAL> TPZMatrix<TFad<6,REAL> >::ReturnNearestValue(TFad<6,REAL> val, TPZVec<TFad<6,REAL> >& Vec, TFad<6,REAL> tol)
1552 {
1553  DebugStop();
1554  TFad<6,REAL> res;
1555  return res;
1556 }
1557 #endif
1558 
1559 template <class TVar>
1560 TVar TPZMatrix<TVar>::ReturnNearestValue(TVar val, TPZVec<TVar>& Vec, TVar tol)
1561 {
1562  TVar exps = val - (TVar)Vec[0];
1563  TVar diff0 = fabs(exps) >= fabs(tol) ? (val - (TVar)Vec[0]) : (TVar)1.E10;
1564  TVar diff1, res = (TVar)Vec[0];
1565  for(int64_t i = 1; i < Vec.NElements(); i++)
1566  {
1567  diff1 = val - (TVar)Vec[i];
1568  diff0 = ( fabs(diff1) >= fabs(tol) && fabs(diff1) < fabs(diff0) ) ? res = Vec[i],diff1 : diff0;
1569  }
1570  return res;
1571 }
1572 
1573 #ifdef _AUTODIFF
1574 template<>
1575 bool TPZMatrix<TFad<6,REAL> >::SolveEigensystemJacobi(int64_t &numiterations, REAL & tol, TPZVec<TFad<6,REAL> > & Eigenvalues, TPZFMatrix<TFad<6,REAL> > & Eigenvectors) const{
1576  DebugStop();
1577  return false;
1578 }
1579 
1580 #endif
1581 
1582 template <class TVar>
1583 bool TPZMatrix<TVar>::SolveEigensystemJacobi(int64_t &numiterations, REAL & tol, TPZVec<TVar> & Eigenvalues, TPZFMatrix<TVar> & Eigenvectors) const{
1584 
1585  int64_t NumIt = numiterations;
1586  REAL tolerance = tol;
1587 
1588 #ifdef PZDEBUG2
1589  if (this->Rows() != this->Cols())
1590  {
1591  PZError << __PRETTY_FUNCTION__ <<
1592  " - Jacobi method of computing eigensystem requires a symmetric square matrix. this->Rows = " << this->Rows() << " - this->Cols() = " << this->Cols() << endl;
1593  return false;
1594  }
1595 
1596  if (this->VerifySymmetry(1.e-8) == false)
1597  {
1598  PZError << __PRETTY_FUNCTION__ <<
1599  " - Jacobi method of computing eigensystem requires a symmetric square matrix. This matrix is not symmetric." << endl;
1600  return false;
1601  }
1602 #endif
1603 
1604  const int64_t size = this->Rows();
1605 
1607  TPZFNMatrix<9,TVar> Matrix(size,size); //fast constructor in case of this is a stress or strain tensor.
1608  for(int64_t i = 0; i < size; i++) for(int64_t j = 0; j < size; j++) Matrix(i,j) = this->Get(i,j);
1609 
1611  bool result = Matrix.SolveEigenvaluesJacobi(numiterations, tol, &Eigenvalues);
1612  if (result == false) return false;
1613 
1615  TPZFNMatrix<3, TVar> VecIni(size,1,0.), VecIni_cp(size,1,0.);
1616 
1617  Eigenvectors.Resize(size, size);
1618  Eigenvectors.Zero();
1619  for(int64_t eigen = 0; eigen < size; eigen++)
1620  {
1621  for(int64_t i = 0; i < size; i++) VecIni.PutVal(i,0,rand());
1622 
1623  TPZFNMatrix<9,TVar> Matrix(*this);
1624 
1625  TVar answ = ReturnNearestValue(Eigenvalues[eigen], Eigenvalues,((TVar)1.E-5));
1626  TVar exp = answ - Eigenvalues[eigen];
1627  if((REAL)(fabs(exp)) > 1.E-5)
1628  {
1629  for(int64_t i = 0; i < size; i++) Matrix.PutVal(i,i, this->GetVal(i,i) - (TVar)(Eigenvalues[eigen] - (TVar)(0.01 * fabs(exp))) );
1630  }
1631  else
1632  {
1633  for(int64_t i = 0; i < size; i++) Matrix.PutVal(i,i, this->GetVal(i,i) - (Eigenvalues[eigen] - TVar(0.01)) );
1634  }
1635 
1637  REAL norm1 = 0.;
1638  for(int64_t i = 0; i < size; i++) norm1 += fabs(VecIni(i,0)) * fabs(VecIni(i,0)); norm1 = sqrt(norm1);
1639  for(int64_t i = 0; i < size; i++) VecIni(i,0) = VecIni(i,0)/(TVar)norm1;
1640 
1641  int64_t count = 0;
1642  double dif = 10., difTemp = 0.;
1643 
1644  while(dif > tolerance && count <= NumIt)
1645  {
1646  for(int64_t i = 0; i < size; i++) VecIni_cp(i,0) = VecIni(i,0);
1647 
1649  Matrix.Solve_LU(&VecIni);
1650 
1652  REAL norm2 = 0.;
1653  for(int64_t i = 0; i < size; i++) norm2 += fabs(VecIni(i,0)) * fabs(VecIni(i,0));
1654  norm2 = sqrt(norm2);
1655 
1656  difTemp = 0.;
1657  for(int64_t i = 0; i < size; i++)
1658  {
1659  VecIni(i,0) = VecIni(i,0)/(TVar)norm2;
1660  TVar exp = (VecIni_cp(i,0) - VecIni(i,0)) * (VecIni_cp(i,0) - VecIni(i,0));
1661  difTemp += fabs(exp);
1662  }
1663  dif = sqrt(difTemp);
1664  count++;
1665  }
1666 
1668  for(int64_t i = 0; i < size; i++)
1669  {
1670  TVar val = VecIni(i,0);
1671  if((REAL)(fabs(val)) < 1.E-5) val = 0.;
1672  Eigenvectors(eigen,i) = val;
1673  }
1674 
1675 #ifdef PZDEBUG2
1676  double norm = 0.;
1677  for(int64_t i = 0; i < size; i++) norm += fabs(VecIni(i,0)) * fabs(VecIni(i,0));
1678  if (fabs(norm - 1.) > 1.e-10)
1679  {
1680  PZError << __PRETTY_FUNCTION__ << endl;
1681  }
1682  if(count > NumIt-1)// O metodo nao convergiu !!!
1683  {
1684  PZError << __PRETTY_FUNCTION__ << endl;
1685 #ifdef LOG4CXX
1686  {
1687  std::stringstream sout;
1688  Print("Matrix for SolveEigensystemJacobi did not converge",sout);
1689  LOGPZ_DEBUG(logger,sout.str());
1690  }
1691 #endif
1692  }
1693 #endif
1694  }
1695 
1696  return true;
1697 
1698 }//method
1699 
1700 template <>
1701 bool TPZMatrix< std::complex< float > >::SolveEigenvaluesJacobi(int64_t &numiterations, REAL & tol, TPZVec< std::complex< float > > * Sort){
1702  DebugStop(); // Does not work with complex numbers. To be implemented in the future.
1703  return false;
1704 }
1705 
1706 template <>
1707 bool TPZMatrix< std::complex< double > >::SolveEigenvaluesJacobi(int64_t &numiterations, REAL & tol, TPZVec< std::complex< double > > * Sort){
1708  DebugStop(); // Does not work with complex numbers. To be implemented in the future.
1709  return false;
1710 }
1711 
1712 template <>
1713 bool TPZMatrix< std::complex< long double > >::SolveEigenvaluesJacobi(int64_t &numiterations, REAL & tol, TPZVec< std::complex< long double > > * Sort){
1714  DebugStop(); // Does not work with complex numbers. To be implemented in the future.
1715  return false;
1716 }
1717 
1718 #ifdef _AUTODIFF
1719 template<>
1720 bool TPZMatrix<TFad<6,REAL> >::SolveEigenvaluesJacobi(int64_t &numiterations, REAL & tol, TPZVec<TFad<6,REAL> > * Sort){
1721  DebugStop();
1722  return false;
1723 }
1724 #endif
1725 
1726 template <class TVar>
1727 bool TPZMatrix<TVar>::SolveEigenvaluesJacobi(int64_t &numiterations, REAL & tol, TPZVec<TVar> * Sort){
1728 
1729 #ifdef PZDEBUG2
1730  if (this->Rows() != this->Cols()){
1731  PZError << __PRETTY_FUNCTION__ << " - Jacobi method of computing eigenvalues requires a symmetric square matrix. this->Rows = " << this->Rows() << " - this->Cols() = " << this->Cols() << endl;
1732  return false;
1733  }
1734 
1735  if (this->VerifySymmetry(1.e-8) == false){
1736  PZError << __PRETTY_FUNCTION__ << " - Jacobi method of computing eigenvalues requires a symmetric square matrix. This matrix is not symmetric." << endl;
1737  return false;
1738  }
1739 #endif
1740 
1741  int64_t iter = 0;
1742  TVar res = (TVar)2. * (TVar)tol;
1743  const int64_t size = this->Rows();
1744  int64_t i, j;
1745  int64_t p = -1, q = -1;
1746  TVar maxval, theta, cost, sint, aux, aux2, Spp, Sqq, Spq;
1747  while (iter < numiterations){
1749  maxval = 0.;
1750  for(i = 0; i < size; i++){
1751  for(j = 0; j < i; j++) {
1752  if( fabs(this->operator ( )(i,j) ) > fabs(maxval) ) {
1753  p = i;
1754  q = j;
1755  maxval = fabs( this->operator ( )(i,j) );
1756  }//if
1757  }//for j
1758  }//for i
1759 
1761  res = maxval;
1762  if ((REAL)(fabs(res)) < tol) break;
1763 
1765  theta = 0.5 * atan((TVar)2. * this->operator ( )(p,q) / (this->operator ( )(q,q) - this->operator ( )(p,p) ) );
1766  cost = cos(theta);
1767  sint = sin(theta);
1768 
1770  for(i = 0; i < size; i++){
1771  if (i != p && i != q){
1772 
1773  aux = this->operator ( )(i,p) * (TVar)cost - this->operator ( )(i,q) * (TVar)sint;
1774  aux2 = this->operator ( )(i,p) * (TVar)sint + this->operator ( )(i,q) * (TVar)cost;
1775 
1776  this->operator ( )(i,p) = aux;
1777  this->operator ( )(p,i) = aux;
1778 
1779  this->operator ( )(i,q) = aux2;
1780  this->operator ( )(q,i) = aux2;
1781 
1782  }//if
1783  }//for i
1784 
1785  Spp = this->operator ( )(p,p) * cost * cost -(TVar)2. * this->operator ( )(p,q) * sint * cost + this->operator ( )(q,q) * sint * sint;
1786  Sqq = this->operator ( )(p,p) * sint * sint +(TVar)2. * this->operator ( )(p,q) * sint * cost + this->operator ( )(q,q) * cost * cost;
1787  Spq = ( this->operator ( )(p,p) - this->operator ( )(q,q) ) * cost * sint + this->operator ( )(p,q)*( cost*cost - sint*sint );
1788 
1789  this->operator ( )(p,p) = Spp;
1790  this->operator ( )(q,q) = Sqq;
1791  this->operator ( )(p,q) = Spq;
1792  this->operator ( )(q,p) = Spq;
1793 
1794  iter++;
1795  }//while
1796 
1798  if (Sort){
1799  multiset< REAL > myset;
1800  for(i = 0; i < size; i++)
1801  { TVar exps = this->operator ( )(i,i);
1802  myset.insert( (TPZExtractVal::val(exps)) );
1803  }
1804 
1805 #ifdef PZDEBUG2
1806  if ((int64_t)myset.size() != size) PZError << __PRETTY_FUNCTION__ << " - ERROR!" << endl;
1807 #endif
1808 
1809  Sort->Resize(size);
1810  multiset< REAL >::iterator w, e = myset.end();
1811  for(i = size - 1, w = myset.begin(); w != e; w++, i--){
1812  Sort->operator [ ](i) = *w;
1813  }//for
1814  }//if (Sort)
1815 
1816 
1817  if ((REAL)(fabs(res)) < tol){
1818  tol = fabs(res);
1819  numiterations = iter;
1820  return true;
1821  }
1822 
1823  tol = fabs(res);
1824  numiterations = iter;
1825  return false;
1826 
1827 }//method
1828 
1829 #ifdef _AUTODIFF
1830 template<>
1831 TFad<6,REAL> TPZMatrix<TFad<6,REAL> >::MatrixNorm(int p, int64_t numiter, REAL tol) const{
1832  DebugStop();
1833  TFad<6,REAL> res;
1834  return res;
1835 }
1836 #endif
1837 
1838 template <class TVar>
1839 TVar TPZMatrix<TVar>::MatrixNorm(int p, int64_t numiter, REAL tol) const{
1840  const int64_t n = this->Rows();
1841  if (!n) return 0.;
1842  if (n != this->Cols()){
1843  PZError << __PRETTY_FUNCTION__
1844  << " - matrix must be square - Rows() = "
1845  << this->Rows() << " - Cols() = "
1846  << this->Cols() << std::endl;
1847  }
1848  switch(p){
1849  case 0:{
1850  TVar max = 0.;
1851  int64_t i, j;
1852  for(i = 0; i < n; i++){
1853  REAL sum = 0.;
1854  for(j = 0; j < n; j++) sum += fabs( this->Get(i,j) );
1855  if (sum > fabs(max)) max = sum;
1856  }
1857  return max;
1858  }
1859  case 1:{
1860  TVar max = 0.;
1861  int64_t i, j;
1862  for(i = 0; i < n; i++){
1863  REAL sum = 0.;
1864  for(j = 0; j < n; j++) sum += fabs( this->Get(j,i) );
1865  if (sum > fabs(max)) max = sum;
1866  }
1867  return max;
1868  }
1869  case 2:{
1870  TPZFMatrix<TVar> transp(n,n);
1871  int64_t i, j, k;
1872  //TRANSPOSE
1873  for(i = 0; i < n; i++) for(j = 0; j < n; j++) transp(i,j) = this->Get(j,i);
1874  //MULTIPLY transp = transp.this
1875  TPZVec<TVar> ROW(n);
1876  for(i = 0; i < n; i++){
1877  for(j = 0; j < n; j++) ROW[j] = transp(i,j);
1878  for(j = 0; j < n; j++){
1879  TVar sum = 0.;
1880  for(k = 0; k < n; k++){
1881  sum += ROW[k] * this->Get(k,j);
1882  }//for k
1883  transp(i,j) = sum;
1884  }//for j
1885  }//for i
1886 
1887  TPZVec<TVar> EigenVal;
1888  bool result = transp.SolveEigenvaluesJacobi(numiter, tol, &EigenVal);
1889  if (result == false) PZError << __PRETTY_FUNCTION__
1890  << " - it was not possible to find Eigenvalues. Iterations = "
1891  << numiter << " - error found = " << tol << std::endl;
1892  return sqrt(EigenVal[0]);
1893  }//case 2
1894  default:{
1895  PZError << __PRETTY_FUNCTION__ << " p = " << p << " is not a correct option" << std::endl;
1896  }
1897  }//switch
1898  return 0.;
1899 }//method
1900 
1901 template <class TVar>
1902 TVar TPZMatrix<TVar>::ConditionNumber(int p, int64_t numiter, REAL tol){
1903  int64_t localnumiter = numiter;
1904  REAL localtol = tol;
1905  TPZFMatrix<TVar> Inv;
1906  TVar thisnorm = this->MatrixNorm(p, localnumiter, localtol);
1907  if (!this->Inverse(Inv,ENoDecompose)){
1908  PZError << __PRETTY_FUNCTION__ << " - it was not possible to compute the inverse matrix." << std:: endl;
1909  return 0.;
1910  }
1911  TVar invnorm = Inv.MatrixNorm(p, numiter, tol);
1912  return thisnorm * invnorm;
1913 }
1914 
1915 template<class TVar>
1916 void TPZMatrix<TVar>::Multiply(const TPZFMatrix<TVar> &A, TPZFMatrix<TVar>&B, int opt) const {
1917  if ((opt==0 && Cols() != A.Rows()) || (opt ==1 && Rows() != A.Rows()))
1918  Error( "Multiply (TPZMatrix<>&,TPZMatrix<TVar> &) <incompatible dimensions>" );
1919  if(!opt && (B.Rows() != Rows() || B.Cols() != A.Cols())) {
1920  B.Redim(Rows(),A.Cols());
1921  }
1922  else if (opt && (B.Rows() != Cols() || B.Cols() != A.Cols())) {
1923  B.Redim(Cols(),A.Cols());
1924  }
1925  MultAdd( A, B, B, 1.0, 0.0, opt);
1926 }
1927 
1928 template<class TVar>
1929 int TPZMatrix<TVar>::Inverse(TPZFMatrix<TVar>&Inv, DecomposeType dec){
1930  const int64_t n = this->Rows();
1931  if (!n) return 0;
1932  if (n != this->Cols()){
1933  PZError << __PRETTY_FUNCTION__
1934  << " - matrix must be square - Rows() = "
1935  << this->Rows() << " - Cols() = "
1936  << this->Cols() << std::endl;
1937  return 0;
1938  }
1939 
1940  Inv.Redim(n,n);
1941  Inv.Zero();
1942  int64_t i;
1943  for(i = 0; i < n; i++) Inv(i,i) = 1.;
1944 
1945  if(dec != ENoDecompose)
1946  {
1947  this->SolveDirect(Inv,dec);
1948  }
1949  else
1950  {
1951  const int issimetric = this->IsSimetric();
1952  if (issimetric) return this->SolveDirect(Inv, ELDLt);
1953  if (!issimetric) return this->SolveDirect(Inv, ELU);
1954  }
1955  return 0;
1956 }//method
1957 
1959 template <class TVar>
1960 void TPZMatrix<TVar>::AutoFill(int64_t nrow, int64_t ncol, int symmetric) {
1961 
1962  Resize(nrow,ncol);
1963  int64_t i, j;
1964  TVar val, sum;
1966  for(i=0;i<Rows();i++) {
1967  sum = 0.0;
1968  j=0;
1969  if (symmetric) {
1970  for (j=0; j<i; j++) {
1971  PutVal(i, j, GetVal(j,i));
1972  sum += fabs(GetVal(i, j));
1973  }
1974  }
1975  for(;j<Cols();j++) {
1976  val = ((TVar)rand())/((TVar)RAND_MAX);
1977  if(!PutVal(i,j,val))
1978  Error("AutoFill (TPZMatrix) failed.");
1979  if(i!=j) sum += fabs(val);
1980  }
1981  if (Rows() == Cols()) {
1983  if(fabs(sum) > fabs(GetVal(i,i))) // Deve satisfazer: |Aii| > SUM( |Aij| ) sobre j != i
1984  PutVal(i,i,sum+(TVar)1.);
1985  // To sure diagonal is not zero.
1986  if(IsZero(sum) && IsZero(GetVal(i,i)))
1987  PutVal(i,i,1.);
1988  }
1989  }
1990 }
1991 
1992 template<class TVar>
1993 int TPZMatrix<TVar>::Solve_LDLt( TPZFMatrix<TVar>* B, std::list<int64_t> &singular ) {
1994 
1995  int result = Decompose_LDLt(singular);
1996  if (result == 0) {
1997  return result;
1998  }
1999 #ifdef LOG4CXX
2000  if (logger->isDebugEnabled())
2001  {
2002  std::stringstream sout;
2003  B->Print("On input " , sout);
2004  LOGPZ_DEBUG(logger, sout.str())
2005  }
2006 #endif
2007  Subst_LForward( B );
2008 #ifdef LOG4CXX
2009  if (logger->isDebugEnabled())
2010  {
2011  std::stringstream sout;
2012  B->Print("Only forward " , sout);
2013  LOGPZ_DEBUG(logger, sout.str())
2014  }
2015 #endif
2016  Subst_Diag( B );
2017 #ifdef LOG4CXX
2018  if (logger->isDebugEnabled())
2019  {
2020  std::stringstream sout;
2021  B->Print("After forward and diagonal " , sout);
2022  LOGPZ_DEBUG(logger, sout.str())
2023  }
2024 #endif
2025  result = Subst_LBackward( B );
2026 #ifdef LOG4CXX
2027  if (logger->isDebugEnabled())
2028  {
2029  std::stringstream sout;
2030  B->Print("Final result " , sout);
2031  LOGPZ_DEBUG(logger, sout.str())
2032  }
2033 #endif
2034  return result;
2035 }
2036 
2037 #include <complex>
2038 template class TPZMatrix< std::complex<float> >;
2039 template class TPZMatrix< std::complex<double> >;
2040 template class TPZMatrix< std::complex<long double> >;
2041 
2042 template class TPZMatrix<long double>;
2043 template class TPZMatrix<double>;
2044 template class TPZMatrix<int64_t>;
2045 template class TPZMatrix<int>;
2046 template class TPZMatrix<float>;
2047 template class TPZMatrix<TPZFlopCounter>;
2048 
2049 #ifdef _AUTODIFF
2050 template class TPZMatrix<TFad<6,REAL> >;
2051 template class TPZMatrix<Fad<float> >;
2052 template class TPZMatrix<Fad<double> >;
2053 template class TPZMatrix<Fad<long double> >;
2054 #endif
2055 
2057 template<class TVar>
2058 std::ostream &operator<<(std::ostream& out,const TPZMatrix<TVar> &A) {
2059  A.Print("operator << ",out);
2060  return out;
2061 }
2062 
2063 template std::ostream &operator<<(std::ostream &out, const TPZMatrix<int> &A);
2064 template std::ostream &operator<<(std::ostream &out, const TPZMatrix<float> &A);
2065 template std::ostream &operator<<(std::ostream &out, const TPZMatrix<double> &A);
2066 template std::ostream &operator<<(std::ostream &out, const TPZMatrix<long double> &A);
2067 template std::ostream &operator<<(std::ostream &out, const TPZMatrix<std::complex<float> > &A);
2068 template std::ostream &operator<<(std::ostream &out, const TPZMatrix<std::complex<double> > &A);
TPZMatrix::SolveSOR
virtual void SolveSOR(int64_t &numiterations, const TPZFMatrix< TVar > &F, TPZFMatrix< TVar > &result, TPZFMatrix< TVar > *residual, TPZFMatrix< TVar > &scratch, const REAL overrelax, REAL &tol, const int FromCurrent=0, const int direction=1)
Solves the linear system using Successive Over Relaxation method (Gauss Seidel). ...
Definition: pzmatrix.cpp:785
TPZMatrix::AddSub
virtual int AddSub(const int64_t sRow, const int64_t sCol, const TPZFMatrix< TVar > &Source)
It adds Source matrix on current matrix from position (sRow, sCol)
Definition: pzmatrix.cpp:625
TPZMatrix::SolveEigenvaluesJacobi
virtual bool SolveEigenvaluesJacobi(int64_t &numiterations, REAL &tol, TPZVec< TVar > *Sort=0)
Transforms this matrix in a diagonal matrix, where the diagonal values are its eigenvalues. This method is efficient only for small matrices.
Definition: pzmatrix.cpp:1727
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
GMRES
int GMRES(Operator &A, Vector &x, const Vector &b, Preconditioner &M, Matrix &H, int &m, int64_t &max_iter, Real &tol, Vector *residual, const int FromCurrent)
GMRES solves the unsymmetric linear system Ax = b using the Generalized Minimum Residual method...
Definition: gmres.h:71
TPZMatrix::Transpose
virtual void Transpose(TPZMatrix< TVar > *const T) const
It makes *T the transpose of current matrix.
Definition: pzmatrix.cpp:697
TPZMatrix::PutSub
virtual int PutSub(const int64_t sRow, const int64_t sCol, const TPZFMatrix< TVar > &Source)
It puts submatrix Source on actual matrix structure.
Definition: pzmatrix.cpp:574
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
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
TPZMatrix::Solve_LU
int Solve_LU(TPZFMatrix< TVar > *B, std::list< int64_t > &singular)
Solves the linear system using LU method .
Definition: pzmatrix.h:900
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.
TPZMatrix::CompareValues
bool CompareValues(TPZMatrix< TVar > &M, TVar tol)
Compare values of this to B, with a precision tolerance tol.
Definition: pzmatrix.cpp:1532
TPZMatrix::SolveSSOR
virtual void SolveSSOR(int64_t &numiterations, const TPZFMatrix< TVar > &F, TPZFMatrix< TVar > &result, TPZFMatrix< TVar > *residual, TPZFMatrix< TVar > &scratch, const REAL overrelax, REAL &tol, const int FromCurrent=0)
Solves the linear system using Symmetric Successive Over Relaxation method (Gauss Seidel)...
Definition: pzmatrix.cpp:828
ENoDecompose
Definition: pzmatrix.h:52
TPZMatrix::Write
void Write(TPZStream &buf, int withclassid) const override
Packs the object structure in a stream of bytes.
Definition: pzmatrix.cpp:1420
tolerance
const double tolerance
Tolerance value (Is zero)
Definition: tpzellipse3d.cpp:29
TPZMatrix::ReturnNearestValue
static TVar ReturnNearestValue(TVar val, TPZVec< TVar > &Vec, TVar tol)
Retorna o valor mais proximo a "val" (exceto valores no intervalo -tol <= val <= +tol) contido no vet...
Definition: pzmatrix.cpp:1560
EFormatted
Definition: pzmatrix.h:55
TPZMatrix::Error
static int Error(const char *msg, const char *msg2=0)
Returns error messages.
Definition: pzmatrix.cpp:1402
Sort
TPZVec< T > & Sort(TPZVec< T > &v)
Sorting the elements into v.
Definition: pzvec_extras.h:142
EInputFormat
Definition: pzmatrix.h:55
TPZSolver
Defines a abstract class of solvers which will be used by matrix classes. Solver. ...
Definition: pzmatrix.h:32
operator*
TPZFMatrix< TVar > operator*(TPZMatrix< TVar > &A, const TPZFMatrix< TVar > &B)
Implements product of matrices: .
Definition: pzmatrix.cpp:128
pzvec.h
Templated vector implementation.
TPZMatrix::SolveEigensystemJacobi
virtual bool SolveEigensystemJacobi(int64_t &numiterations, REAL &tol, TPZVec< TVar > &Eigenvalues, TPZFMatrix< TVar > &Eigenvectors) const
Compute Eigenvalues and Eigenvectors of this matrix. This method is efficient only for small matrice...
Definition: pzmatrix.cpp:1583
MatrixOutputFormat
MatrixOutputFormat
Defines output format.
Definition: pzmatrix.h:55
TPZMatrix::Identity
virtual void Identity()
Converts the matrix in an identity matrix.
Definition: pzmatrix.cpp:199
TPZMatrix::Resize
virtual int Resize(const int64_t newRows, const int64_t newCols)
Redimensions a matriz keeping the previous values.
Definition: pzmatrix.h:278
TPZMatrix::fRow
int64_t fRow
Number of rows in matrix.
Definition: pzmatrix.h:779
TPZMatrix::Get
virtual const TVar & Get(const int64_t row, const int64_t col) const
Get value with bound checking.
Definition: pzmatrix.h:855
TPZMatrix::PrepareZ
void PrepareZ(const TPZFMatrix< TVar > &y, TPZFMatrix< TVar > &z, const TVar beta, const int opt) const
Is an auxiliar method used by MultiplyAdd.
Definition: pzmatrix.cpp:135
TPZMatrix::Inverse
int Inverse(TPZFMatrix< TVar > &Inv, DecomposeType dec)
It makes Inv =[this]. IMPORTANT OBSERVATION –> The original matrix (calling object) no is more equal...
Definition: pzmatrix.cpp:1929
bicg.h
Contains the implementation of the BiCG function which solves the unsymmetric linear system using Pre...
TPZMatrix::fDecomposed
char fDecomposed
Decomposition type used to decompose the current matrix.
Definition: pzmatrix.h:783
TPZMatrix::Subst_LBackward
virtual int Subst_LBackward(TPZFMatrix< TVar > *b) const
Computes B = Y, where A*Y = B, A is upper triangular with A(i,i)=1.
Definition: pzmatrix.cpp:1358
Fad
Definition: fad.h:54
BiCGSTAB
int BiCGSTAB(const Matrix &A, Vector &x, const Vector &b, Preconditioner &M, int64_t &max_iter, Real &tol, Vector *residual, const int FromCurrent)
BiCGSTAB solves the unsymmetric linear system using the Preconditioned BiConjugate Gradient Stabiliz...
Definition: bicgstab.h:27
val
REAL val(STATE &number)
Returns value of the variable.
Definition: pzartdiff.h:23
TPZMatrix::Subst_Forward
virtual int Subst_Forward(TPZFMatrix< TVar > *b) const
Computes B = Y, where A*Y = B, A is lower triangular.
Definition: pzmatrix.cpp:1284
tmp
AutoPointerMutexArrayInit tmp
Definition: tpzautopointer.cpp:42
TPZMatrix::~TPZMatrix
virtual ~TPZMatrix()
Simple destructor.
Definition: pzmatrix.cpp:43
TPZMatrix::fDefPositive
char fDefPositive
Definite Posistiveness of current matrix.
Definition: pzmatrix.h:785
TPZMatrix::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
It computes z = beta * y + alpha * opt(this)*x but z and x can not overlap in memory.
Definition: pzmatrix.cpp:166
gmres.h
Contains the implementation of the CMRES function which solves the unsymmetric linear system using th...
operator+
TPZFMatrix< TVar > operator+(const TPZMatrix< TVar > &A, const TPZMatrix< TVar > &B)
Implements sum of matrices: .
Definition: pzmatrix.cpp:108
TPZMatRed::MultAdd
void MultAdd(const TPZFMatrix< TVar > &x, const TPZFMatrix< TVar > &y, TPZFMatrix< TVar > &z, const TVar alpha, const TVar beta, const int opt) const override
It computes z = beta * y + alpha * opt(this)*x but z and x can not overlap in memory.
Definition: pzmatred.cpp:521
TPZMatrix::gZero
static TVar gZero
Initializing value to static variable.
Definition: pzmatrix.h:786
bicgstab.h
Contains the implementation of the BiCGSTAB function which solves the unsymmetric linear system using...
TPZMatrix::GetSub
virtual int GetSub(const int64_t sRow, const int64_t sCol, const int64_t rowSize, const int64_t colSize, TPZFMatrix< TVar > &Target) const
Gets submatrix storing it on Target.
Definition: pzmatrix.cpp:602
TPZMatrix::Decompose_LDLt
virtual int Decompose_LDLt()
Decomposes the current matrix using LDLt.
Definition: pzmatrix.cpp:1184
TPZMatrix::GetVal
virtual const TVar & GetVal(const int64_t, const int64_t) const
Get values without bounds checking This method is faster than "Get" if DEBUG is defined.
Definition: pzmatrix.cpp:56
BiCG
int BiCG(Matrix &A, Vector &x, const Vector &b, Preconditioner &M, int64_t &max_iter, Real &tol)
Definition: bicg.h:27
TPZMatRed::s
virtual TVar & s(const int64_t row, const int64_t col) override
The operators check on the bounds if the DEBUG variable is defined.
Definition: pzmatred.cpp:113
TPZMatRed::Print
void Print(const char *name=NULL, std::ostream &out=std::cout, const MatrixOutputFormat=EFormatted) const override
Prints the object data structure.
Definition: pzmatred.cpp:454
TPZFMatrix::Zero
int Zero() override
Makes Zero all the elements.
Definition: pzfmatrix.h:651
TPZVec::Resize
virtual void Resize(const int64_t newsize, const T &object)
Resizes the vector object reallocating the necessary storage, copying the existing objects to the new...
Definition: pzvec.h:373
TPZMatrix::Solve_LDLt
int Solve_LDLt(TPZFMatrix< TVar > *B, std::list< int64_t > &singular)
Solves the linear system using LDLt method .
Definition: pzmatrix.cpp:1993
sin
sin
Definition: fadfunc.h:63
pzgeom::tol
static const double tol
Definition: pzgeoprism.cpp:23
operator>>
std::istream & operator>>(std::istream &in, TPZMatrix< TVar > &A)
Overload >> operator to input data of the matrix.
Definition: pzmatrix.cpp:246
TPZMatrix::SolveGMRES
virtual void SolveGMRES(int64_t &numiterations, TPZSolver< TVar > &preconditioner, TPZFMatrix< TVar > &H, int &numvectors, const TPZFMatrix< TVar > &F, TPZFMatrix< TVar > &result, TPZFMatrix< TVar > *residual, REAL &tol, const int FromCurrent)
Solves the linear system using Generalized Minimal Residual (GMRES) method. .
Definition: pzmatrix.cpp:896
TPZStream::Write
virtual void Write(const bool val)
Definition: TPZStream.cpp:8
TPZMatrix::SolveBICGStab
virtual void SolveBICGStab(int64_t &numiterations, TPZSolver< TVar > &preconditioner, const TPZFMatrix< TVar > &F, TPZFMatrix< TVar > &result, TPZFMatrix< TVar > *residual, REAL &tol, const int FromCurrent=0)
Solves the linear system using Bi-Conjugate Gradient stabilized method. .
Definition: pzmatrix.cpp:1024
Norm
TVar Norm(const TPZFMatrix< TVar > &A)
Returns the norm of the matrix A.
DebugStop
#define DebugStop()
Returns a message to user put a breakpoint in.
Definition: pzerror.h:20
LOGPZ_DEBUG
#define LOGPZ_DEBUG(A, B)
Define log for debug info.
Definition: pzlog.h:87
TPZMatrix::Add
virtual void Add(const TPZMatrix< TVar > &A, TPZMatrix< TVar > &res) const
It adds itself to TPZMatrix<TVar>A putting the result in res.
Definition: pzmatrix.cpp:64
TFad
Definition: tfad.h:64
TPZMatrix::SolveIR
virtual void SolveIR(int64_t &numiterations, TPZSolver< TVar > &preconditioner, const TPZFMatrix< TVar > &F, TPZFMatrix< TVar > &result, TPZFMatrix< TVar > *residual, REAL &tol, const int FromCurrent=0)
Solves the linear system using IR method. .
Definition: pzmatrix.cpp:1070
ELU
Definition: pzmatrix.h:52
TPZMatrix::Compare
virtual bool Compare(TPZSavable *copy, bool override=false) override
Compare the object for identity with the object pointed to, eventually copy the object.
Definition: pzmatrix.cpp:1433
TPZMatrix::SolveCG
virtual void SolveCG(int64_t &numiterations, TPZSolver< TVar > &preconditioner, const TPZFMatrix< TVar > &F, TPZFMatrix< TVar > &result, TPZFMatrix< TVar > *residual, REAL &tol, const int FromCurrent=0)
Solves the linear system using Conjugate Gradient method. .
Definition: pzmatrix.cpp:850
TPZMatrix::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)
Solves the linear system using Jacobi method. .
Definition: pzmatrix.cpp:743
TPZMatrix::Rows
int64_t Rows() const
Returns number of rows.
Definition: pzmatrix.h:803
TPZMatrix::Subst_Diag
virtual int Subst_Diag(TPZFMatrix< TVar > *b) const
Computes B = Y, where A*Y = B, A is diagonal matrix.
Definition: pzmatrix.cpp:1384
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
TPZMatrix::Put
virtual int Put(const int64_t row, const int64_t col, const TVar &value)
Put values with bounds checking if DEBUG variable is defined.
Definition: pzmatrix.h:836
TPZMatrix::AddKel
virtual void AddKel(TPZFMatrix< TVar > &elmat, TPZVec< int64_t > &destinationindex)
Add a contribution of a stiffness matrix.
Definition: pzmatrix.cpp:506
TPZMatrix::InsertSub
virtual int InsertSub(const int64_t sRow, const int64_t sCol, const int64_t rowSize, const int64_t colSize, const int64_t pRow, const int64_t pCol, TPZMatrix< TVar > *Target) const
Inserts a submatrix from current object on matrix *Target with no redimentioning.
Definition: pzmatrix.cpp:647
test.res
string res
Definition: test.py:151
TPZMatrix::Read
void Read(TPZStream &buf, void *context) override
Unpacks the object structure from a stream of bytes.
Definition: pzmatrix.cpp:1413
MIN
#define MIN(a, b)
Gets minime value between a and b.
Definition: pzreal.h:85
TPZMatrix::Input
virtual void Input(std::istream &in=std::cin)
Input operation.
Definition: pzmatrix.cpp:227
LOGPZ_ERROR
#define LOGPZ_ERROR(A, B)
Define log for errors (cout)
Definition: pzlog.h:93
TPZMatrix::Redim
virtual int Redim(const int64_t newRows, const int64_t newCols)
Redimensions the matrix reinitializing it with zero.
Definition: pzmatrix.h:289
TPZExtractVal::val
static REAL val(const int number)
Definition: pzextractval.h:21
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
TPZMatRed::Redim
int Redim(const int64_t dim, const int64_t dim00) override
Redim: Set the dimension of the complete matrix and reduced matrix.
Definition: pzmatred.cpp:481
TPZMatrix::SolveBICG
virtual void SolveBICG(int64_t &numiterations, TPZSolver< TVar > &preconditioner, const TPZFMatrix< TVar > &F, TPZFMatrix< TVar > &result, REAL &tol)
Solves the linear system using Bi-Conjugate Gradient method. .
Definition: pzmatrix.cpp:974
TPZMatrix::Subst_LForward
virtual int Subst_LForward(TPZFMatrix< TVar > *b) const
Computes B = Y, where A*Y = B, A is lower triangular with A(i,i)=1.
Definition: pzmatrix.cpp:1333
TPZMatrix::Dim
virtual int64_t Dim() const
Returns the dimension of the matrix if the matrix is square.
Definition: pzmatrix.h:892
pzmatrix.h
Contains TPZMatrix<TVar>class, root matrix class.
TPZMatrix::fCol
int64_t fCol
Number of cols in matrix.
Definition: pzmatrix.h:781
TPZMatRed::PutVal
virtual int PutVal(const int64_t row, const int64_t col, const TVar &value) override
Put and Get values without bounds checking these methods are faster than "Put" e "Get" if DEBUG is de...
Definition: pzmatred.cpp:88
TPZMatrix::Decompose_Cholesky
virtual int Decompose_Cholesky()
Decomposes the current matrix using Cholesky method. The current matrix has to be symmetric...
Definition: pzmatrix.cpp:1252
TPZMatrix::PutVal
virtual int PutVal(const int64_t, const int64_t, const TVar &val)
Put values without bounds checking This method is faster than "Put" if DEBUG is defined.
Definition: pzmatrix.h:154
CG
int CG(Matrix &A, Vector &x, const Vector &b, Preconditioner &M, Vector *residual, int64_t &max_iter, Real &tol, const int FromCurrent)
CG solves the symmetric positive definite linear system using the Conjugate Gradient method...
Definition: cg.h:31
TPZMatrix::Substract
virtual void Substract(const TPZMatrix< TVar > &A, TPZMatrix< TVar > &result) const
It substracts A from storing the result in result.
Definition: pzmatrix.cpp:75
IR
int IR(Matrix &A, Vector &x, const Vector &b, Preconditioner &M, Vector *residual, int64_t &max_iter, Real &tol, const int FromCurrent)
IR solves the unsymmetric linear system Ax = b using Iterative Refinement (preconditioned Richardson ...
Definition: ir.h:24
TPZMatrix::VerifySymmetry
virtual int VerifySymmetry(REAL tol=1.e-13) const
Checks if current matrix value is symmetric.
Definition: pzmatrix.cpp:1501
atan
expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ atan
Definition: tfadfunc.h:85
TPZMatrix::AutoFill
void AutoFill(int64_t nrow, int64_t ncol, int symmetric)
Fill matrix storage with randomic values.
Definition: pzmatrix.cpp:1960
EMatlabNonZeros
Definition: pzmatrix.h:55
EMatrixMarket
Definition: pzmatrix.h:55
TPZMatrix::Multiply
virtual void Multiply(const TPZFMatrix< TVar > &A, TPZFMatrix< TVar > &res, int opt=0) const
It mutiplies itself by TPZMatrix<TVar>A putting the result in res.
Definition: pzmatrix.cpp:1916
pzsolve.h
Contains TPZSolver class which defines a abstract class of solvers which will be used by matrix class...
TPZMatRed::IsSimetric
virtual int IsSimetric() const override
returns 1 or 0 depending on whether the fK00 matrix is zero or not
Definition: pzmatred.cpp:61
TPZMatrix::Decompose_LU
virtual int Decompose_LU()
Definition: pzmatrix.cpp:1126
EMathematicaInput
Definition: pzmatrix.h:55
pzstepsolver.h
Contains TPZStepSolver class which defines step solvers class.
TPZMatRed::fDim1
int64_t fDim1
Definition: pzmatred.h:253
TPZMatrix::Cols
int64_t Cols() const
Returns number of cols.
Definition: pzmatrix.h:809
exp
expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ exp
Definition: tfadfunc.h:125
TPZMatrix::Print
virtual void Print(std::ostream &out) const
Definition: pzmatrix.h:253
TPZFMatrix::Resize
int Resize(const int64_t newRows, const int64_t wCols) override
Redimension a matrix, but maintain your elements.
Definition: pzfmatrix.cpp:1016
TPZStream
Defines the interface for saving and reading data. Persistency.
Definition: TPZStream.h:50
TPZVec::NElements
int64_t NElements() const
Returns the number of elements of the vector.
Definition: pzvec.h:190
TPZMatrix::Substitution
virtual int Substitution(TPZFMatrix< TVar > *B) const
Computes Forward and Backward substitution for a "LU" decomposed matrix.
Definition: pzmatrix.cpp:1151
TPZMatrix::Solve_Cholesky
virtual int Solve_Cholesky(TPZFMatrix< TVar > *B)
Solves the linear system using Cholesky method .
Definition: pzmatrix.h:919
TPZMatrix::MatrixNorm
TVar MatrixNorm(int p, int64_t numiter=2000000, REAL tol=1.e-10) const
Computes the matrix norm of this.
Definition: pzmatrix.cpp:1839
TPZMatRed::GetVal
virtual 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: pzmatred.cpp:101
TPZMatrix::ConditionNumber
TVar ConditionNumber(int p, int64_t numiter=2000000, REAL tol=1.e-10)
Computes the matrix condition number of this.
Definition: pzmatrix.cpp:1902
ECholesky
Definition: pzmatrix.h:52
TPZSavable
This class defines the interface to save and restore objects from TPZStream objects. Persistency.
Definition: TPZSavable.h:67
cos
TPZFlopCounter cos(const TPZFlopCounter &orig)
Returns the cosine in radians and increments the counter of the Cosine.
Definition: pzreal.h:514
TPZFMatrix::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: pzfmatrix.h:548
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
TPZMatrix::Simetrize
virtual void Simetrize()
Simetrizes copies upper plan to the lower plan, making its data simetric.
Definition: pzmatrix.cpp:90
TPZFMatrix::g
TVar & g(const int64_t row, const int64_t col) const
Definition: pzfmatrix.h:594
cg.h
Contains the implementation of the CG function which solves the symmetric positive definite linear sy...
TPZMatrix::SolveDirect
virtual int SolveDirect(TPZFMatrix< TVar > &F, const DecomposeType dt, std::list< int64_t > &singular)
Solves the linear system using Direct methods.
Definition: pzmatrix.cpp:710
ir.h
Contains the implementation of the IR function which solves the unsymmetric linear system using the I...
TPZFNMatrix
Non abstract class which implements full matrices with preallocated storage with (N+1) entries...
Definition: pzfmatrix.h:716
PZError
#define PZError
Defines the output device to error messages and the DebugStop() function.
Definition: pzerror.h:15
DecomposeType
DecomposeType
Defines decomposition type for any matrix classes.
Definition: pzmatrix.h:52
TPZStream::Read
virtual void Read(bool &val)
Definition: TPZStream.cpp:91
TPZMatrix::Subst_Backward
virtual int Subst_Backward(TPZFMatrix< TVar > *b) const
Computes B = Y, where A*Y = B, A is upper triangular.
Definition: pzmatrix.cpp:1309
TPZMatrix
Root matrix class (abstract). Matrix.
Definition: pzmatrix.h:60
operator-
TPZFMatrix< TVar > operator-(const TPZMatrix< TVar > &A, const TPZMatrix< TVar > &B)
Implements difference of matrices: .
Definition: pzmatrix.cpp:118
Generated on Sun Aug 16 2020 11:17:05 for NeoPZ by   doxygen 1.8.13