NeoPZ
pzsfulmat.cpp
Go to the documentation of this file.
1 
6 #include <math.h>
7 #include <stdlib.h>
8 
9 #include "pzfmatrix.h"
10 #include "pzsfulmat.h"
11 //#include "pzerror.h"
12 
13 #include <sstream>
14 #include "pzlog.h"
15 #ifdef LOG4CXX
16 static LoggerPtr logger(Logger::getLogger("pz.matrix.tpzsfmatrix"));
17 #endif
18 
19 using namespace std;
20 
21 /*******************/
22 /*** Constructor ***/
23 template<class TVar>
24 TPZSFMatrix<TVar> ::TPZSFMatrix(const int64_t dim )
25 : TPZRegisterClassId(&TPZSFMatrix::ClassId),
26 TPZMatrix<TVar>( dim, dim )
27 {
28  int64_t size = Size();
29  fElem = new TVar[ size ] ;
30 
31  if ( fElem == NULL )
32  TPZMatrix<TVar> ::Error(__PRETTY_FUNCTION__, "Constructor <memory allocation error>." );
33 
34  // Zera a Matriz.
35  Zero();
36 }
37 
38 /*********************************/
39 /*** Constructor( TPZSFMatrix& ) ***/
40 template<class TVar>
41 TPZSFMatrix<TVar> ::TPZSFMatrix (const TPZSFMatrix<TVar> & A)
42 : TPZRegisterClassId(&TPZSFMatrix::ClassId),
43 TPZMatrix<TVar> ( A.Dim(), A.Dim() )
44 {
45  int64_t size = Size();
46  fElem = new TVar[size] ;
47  if ( fElem == NULL )
48  TPZMatrix<TVar> ::Error(__PRETTY_FUNCTION__, "Constructor <memory allocation error>." );
49 
50  // Copia a matriz
51  TVar *src = A.fElem;
52  TVar *dst = fElem;
53  TVar *end = &fElem[Size()];
54  while ( dst < end )
55  *dst++ = *src++;
56 }
57 
58 
59 
60 /*******************************/
61 /*** Constructor( TPZMatrix<TVar> & ) ***/
62 
63 template<class TVar>
64 TPZSFMatrix<TVar> ::TPZSFMatrix(const TPZMatrix<TVar> &A )
65 : TPZRegisterClassId(&TPZSFMatrix::ClassId),
66 TPZMatrix<TVar> ( A )
67 {
68  int64_t size = Size();
69  fElem = new TVar[size] ;
70 
71  if ( fElem == NULL )
72  TPZMatrix<TVar> ::Error(__PRETTY_FUNCTION__, "Constructor <memory allocation error>." );
73 
74  // Copia a matriz
75  TVar *dst = fElem;
76  for ( int64_t r = 0; r < this->Dim(); r++ )
77  for ( int64_t c = 0; c <= r; c++ )
78  *dst++ = A.GetVal( r, c );
79 }
80 
81 
82 
83 /******************/
84 /*** Destructor ***/
85 
86 template<class TVar>
87 TPZSFMatrix<TVar> ::~TPZSFMatrix ()
88 {
89  if ( fElem != NULL )
90  delete []fElem;
91 }
92 
93 
94 
95 /******** Operacoes com matrizes FULL simetricas ********/
96 
97 /******************/
98 /*** Operator = ***/
99 
100 template<class TVar>
101 TPZSFMatrix<TVar> &
102 TPZSFMatrix<TVar> ::operator=(const TPZSFMatrix<TVar> &A )
103 {
104  if ( this->Dim() != A.Dim() )
105  {
106  if ( fElem != NULL )
107  {
108  delete( fElem );
109  }
110  int64_t size = A.Size();
111  fElem = new TVar[ size ] ;
112  if ( fElem == NULL )
113  TPZMatrix<TVar> ::Error(__PRETTY_FUNCTION__, "Operator= <memory allocation error>." );
114  }
115 
116  this->fRow = this->fCol = A.Dim();
117 
118  // Copia a matriz
119  TVar *src = A.fElem;
120  TVar *dst = fElem;
121  TVar *end = &fElem[ Size() ];
122  while ( dst < end )
123  *dst++ = *src++;
124 
125  this->fDecomposed = A.fDecomposed;
126 
127  return *this;
128 }
129 
130 
131 
132 /*******************************/
133 /*** Operator+( TPZSFMatrix & ) ***/
134 
135 template<class TVar>
136 TPZSFMatrix<TVar>
137 TPZSFMatrix<TVar> ::operator+(const TPZSFMatrix<TVar> &A ) const
138 {
139  if ( A.Dim() != this->Dim() )
140  TPZMatrix<TVar> ::Error(__PRETTY_FUNCTION__, "Operator+ <matrixs with different dimensions>" );
141 
142  TPZSFMatrix<TVar> res( this->Dim() );
143  TVar *pm = fElem;
144  TVar *pa = A.fElem;
145  TVar *pr = res.fElem;
146  TVar *end = &(res.fElem[Size()]);
147 
148  while ( pr < end )
149  *pr++ = (*pm++) + (*pa++);
150 
151  return( res );
152 }
153 
154 
155 
156 /*******************************/
157 /*** Operator-( TPZSFMatrix & ) ***/
158 template<class TVar>
159 TPZSFMatrix<TVar>
160 TPZSFMatrix<TVar> ::operator-(const TPZSFMatrix<TVar> &A ) const
161 {
162  if ( A.Dim() != this->Dim() )
163  TPZMatrix<TVar> ::Error(__PRETTY_FUNCTION__, "Operator- <matrixs with different dimensions>" );
164 
165  TPZSFMatrix<TVar> res( this->Dim() );
166  TVar *pm = fElem;
167  TVar *pa = A.fElem;
168  TVar *pr = res.fElem;
169  TVar *end = &(res.fElem[Size()]);
170 
171  while ( pr < end )
172  *pr++ = (*pm++) - (*pa++);
173 
174  return( res );
175 }
176 
177 
178 
179 /********************************/
180 /*** Operator+=( TPZSFMatrix & ) ***/
181 
182 template<class TVar>
183 TPZSFMatrix<TVar> &
184 TPZSFMatrix<TVar> ::operator+=(const TPZSFMatrix<TVar> &A )
185 {
186  if ( A.Dim() != this->Dim() )
187  TPZMatrix<TVar> ::Error(__PRETTY_FUNCTION__, "Operator+= <matrixs with different dimensions>" );
188 
189  TVar *src = A.fElem;
190  TVar *dst = fElem;
191  TVar *end = &fElem[ Size() ];
192  while ( dst < end )
193  *dst++ += *src++;
194 
195  this->fDecomposed = 0;
196  return( *this );
197 }
198 
199 
200 
201 /*******************************/
202 /*** Operator-=( TPZSFMatrix & ) ***/
203 
204 template<class TVar>
205 TPZSFMatrix<TVar> &
206 TPZSFMatrix<TVar> ::operator-=(const TPZSFMatrix<TVar> &A )
207 {
208  if ( A.Dim() != this->Dim() )
209  TPZMatrix<TVar> ::Error(__PRETTY_FUNCTION__,"Operator-= <matrixs with different dimensions>" );
210 
211  TVar *src = A.fElem;
212  TVar *dst = fElem;
213  TVar *end = &fElem[ Size() ];
214  while ( dst < end )
215  *dst++ -= *src++;
216 
217  this->fDecomposed = 0;
218  return( *this );
219 }
220 
221 
222 
223 /******** Operacoes com matrizes GENERICAS ********/
224 
225 /******************/
226 /*** Operator = ***/
227 
228 template<class TVar>
229 TPZSFMatrix<TVar> &
230 TPZSFMatrix<TVar> ::operator=(const TPZMatrix<TVar> &A )
231 {
232  int64_t newDim = Min( A.Rows(), A.Cols() );
233  int64_t size = newDim * (newDim + 1) / 2;
234 
235  if ( newDim != this->Dim() )
236  {
237  if ( fElem != NULL )
238  delete( fElem );
239  fElem = new TVar[ size ] ;
240  }
241 
242  // Copia a matriz.
243  TVar *dst = fElem;
244  for ( int64_t c = 0; c < newDim; c++ )
245  for ( int64_t r = 0; r <= c; r++ )
246  *dst++ = A.Get( r, c );
247 
248  this->fRow = this->fCol = newDim;
249 
250  // Ajusta Status de decomposicao.
251 
252  this->SetIsDecomposed(A.IsDecomposed());
253  this->fDefPositive = (char)A.IsDefPositive();
254 
255  return( *this );
256 }
257 
258 
259 
260 /******************/
261 /*** Operator + ***/
262 
263 template<class TVar>
264 TPZSFMatrix<TVar>
265 TPZSFMatrix<TVar> ::operator+(const TPZMatrix<TVar> &A ) const
266 {
267  if ( this->Dim() != A.Dim() )
268  TPZMatrix<TVar> ::Error(__PRETTY_FUNCTION__,"Operator+ <matrixs with incoompatible dimensions" );
269 
270  TPZSFMatrix<TVar> res( this->Dim() );
271  TVar *pm = fElem;
272  TVar *pr = res.fElem;
273  for ( int64_t c = 0; c < this->Dim(); c++ )
274  for ( int64_t r = 0; r <= c; r++ )
275  *pr++ = (*pm++) + A.Get( r, c );
276 
277  return( *this );
278 }
279 
280 
281 
282 /******************/
283 /*** Operator - ***/
284 
285 template<class TVar>
286 TPZSFMatrix<TVar>
287 TPZSFMatrix<TVar> ::operator-( const TPZMatrix<TVar> &A ) const
288 {
289  if ( this->Dim() != A.Dim() )
290  TPZMatrix<TVar> ::Error(__PRETTY_FUNCTION__,"Operator- <matrixs with incoompatible dimensions" );
291 
292  TPZSFMatrix<TVar> res( this->Dim() );
293  TVar *pm = fElem;
294  TVar *pr = res.fElem;
295  for ( int64_t c = 0; c < this->Dim(); c++ )
296  for ( int64_t r = 0; r <= c; r++ )
297  *pr++ = (*pm++) - A.Get( r, c );
298 
299  return( *this );
300 }
301 
302 
303 
304 /*******************/
305 /*** Operator += ***/
306 template<class TVar>
307 TPZSFMatrix<TVar> &
308 TPZSFMatrix<TVar> ::operator+=(const TPZMatrix<TVar> &A )
309 {
310  if (this->Dim() != A.Dim() )
311  TPZMatrix<TVar> ::Error(__PRETTY_FUNCTION__,"Operator+= (TPZMatrix<>&) <different dimensions>" );
312 
313  TVar *pm = fElem;
314  for ( int64_t c = 0; c < this->Dim(); c++ )
315  for ( int64_t r = 0; r <= c; r++ )
316  *pm++ += A.Get( r, c );
317 
318  this->fDecomposed = 0;
319  return( *this );
320 }
321 
322 
323 
324 /*******************/
325 /*** Operator -= ***/
326 
327 template<class TVar>
328 TPZSFMatrix<TVar> &
329 TPZSFMatrix<TVar> ::operator-=(const TPZMatrix<TVar> &A )
330 {
331  if ( this->Dim() != A.Dim() )
332  TPZMatrix<TVar> ::Error(__PRETTY_FUNCTION__,"Operator-= (TPZMatrix<>&) <different dimensions>" );
333 
334  TVar *pm = fElem;
335  for ( int64_t c = 0; c < this->Dim(); c++ )
336  for ( int64_t r = 0; r <= c; r++ )
337  *pm++ -= A.Get( r, c );
338 
339  this->fDecomposed = 0;
340  return( *this );
341 }
342 
343 
344 
345 /******** Operacoes com valores NUMERICOS ********/
346 
347 /******************/
348 /*** Operator = ***/
349 
350 template<class TVar>
351 TPZSFMatrix<TVar> &
352 TPZSFMatrix<TVar> ::operator=(const TVar value )
353 {
354  TVar *dst = fElem;
355  TVar *end = &fElem[ Size() ];
356  while ( dst < end )
357  *dst++ = value;
358  this->fDecomposed = 0;
359  this->fDefPositive = 0;
360  return( *this );
361 }
362 
363 
364 
365 /**************************/
366 /*** Operator+( value ) ***/
367 
368 template<class TVar>
369 TPZSFMatrix<TVar>
370 TPZSFMatrix<TVar> ::operator+(const TVar value ) const
371 {
372  TPZSFMatrix<TVar> res( this->Dim() );
373 
374  TVar *dst = res.fElem;
375  TVar *src = fElem;
376  TVar *end = &fElem[ Size() ];
377  while ( src < end )
378  *dst++ = (*src++) + value;
379 
380  return( res );
381 }
382 
383 
384 
385 /**************************/
386 /*** Operator*( value ) ***/
387 
388 template<class TVar>
389 TPZSFMatrix<TVar>
390 TPZSFMatrix<TVar> ::operator*(const TVar value ) const
391 {
392  TPZSFMatrix<TVar> res( this->Dim() );
393 
394  TVar *dst = res.fElem;
395  TVar *src = fElem;
396  TVar *end = &fElem[ Size() ];
397  while ( src < end )
398  *dst++ = (*src++) * value;
399 
400  return( res );
401 }
402 
403 
404 
405 /***************************/
406 /*** Operator+=( value ) ***/
407 
408 template<class TVar>
409 TPZSFMatrix<TVar> &
410 TPZSFMatrix<TVar> ::operator+=( TVar value )
411 {
412  TVar *dst = fElem;
413  TVar *end = &fElem[ Size() ];
414  while ( dst < end )
415  *dst++ += value;
416  this->fDecomposed = 0;
417  return( *this );
418 }
419 
420 
421 
422 /***************************/
423 /*** Operator*=( value ) ***/
424 
425 template<class TVar>
426 TPZSFMatrix<TVar> &
427 TPZSFMatrix<TVar> ::operator*=( TVar value )
428 {
429  TVar *dst = fElem;
430  TVar *end = &fElem[ Size() ];
431  while ( dst < end )
432  *dst++ *= value;
433  this->fDecomposed = 0;
434  return( *this );
435 }
436 
437 
438 
439 /**************/
440 /*** Resize ***/
441 template<class TVar>
442 int
443 TPZSFMatrix<TVar> ::Resize( int64_t newDim , int64_t )
444 {
445  if ( newDim == this->Dim() )
446  return( 1 );
447 
448  int64_t newSize = newDim * (newDim + 1) / 2;
449  int64_t oldSize = Size();
450  TVar *newElem = new TVar[newSize] ;
451  if ( newElem == NULL )
452  return TPZMatrix<TVar> ::Error(__PRETTY_FUNCTION__, "Resize <memory allocation error>." );
453 
454  int64_t minSize = Min( newSize, oldSize );
455  TVar *src = fElem;
456  TVar *dst = newElem;
457  TVar *end = &fElem[ minSize ];
458 
459  // Copia os elementos antigos para a nova matriz.
460  while ( src < end )
461  *dst++ = *src++;
462 
463  // Preenche os elementos que sobrarem (se sobrarem) com ZEROS.
464  end = &newElem[ newSize ];
465  while ( dst < end )
466  *dst++ = 0.0;
467 
468  if ( fElem != NULL )
469  delete( fElem );
470  fElem = newElem;
471  this->fRow = this->fCol = newDim;
472  this->fDecomposed = 0;
473  return( 1 );
474 }
475 
476 
477 
478 /*************/
479 /*** Redim ***/
480 
481 template<class TVar>
482 int
483 TPZSFMatrix<TVar> ::Redim( int64_t newDim , int64_t)
484 {
485  // Se for preciso, desaloca a matriz antiga e aloca uma
486  // nova com o novo tamanho.
487  if ( newDim != this->Dim() )
488  {
489  this->fRow = this->fCol = newDim;
490  if ( fElem != NULL )
491  delete( fElem );
492  fElem = new TVar[Size()] ;
493  }
494 
495  // Zera a matriz.
496  TVar *dst = fElem;
497  TVar *end = &fElem[ Size() ];
498  while ( dst < end )
499  *dst++ = 0.;
500  this->fDecomposed = 0;
501  this->fDefPositive = 0;
502 
503 
504  return( 1 );
505 }
506 
507 
508 template<class TVar>
509 int
510 TPZSFMatrix<TVar> ::Zero()
511 {
512  TVar *dst = fElem;
513  TVar *end = &fElem[ Size() ];
514  while ( dst < end )
515  *dst++ = 0.;
516  this->fDecomposed = 0;
517  this->fDefPositive = 0;
518  return( 1 );
519 }
520 /******** Resolucao de Sistemas ********/
521 
522 /**************************/
523 /*** Decompose Cholesky ***/
524 template<class TVar>
525 int
526 TPZSFMatrix<TVar> ::Decompose_Cholesky(std::list<int64_t> &singular)
527 {
528  return Decompose_Cholesky();
529 }
530 
531 template<class TVar>
532 int
533 TPZSFMatrix<TVar> ::Decompose_Cholesky()
534 {
535  if ( this->fDecomposed ) TPZMatrix<TVar> ::Error(__PRETTY_FUNCTION__, "Decompose_Cholesky <Matrix already Decomposed>" );
536  TVar *ptr_k = fElem;
537  for ( int64_t k = 0; k < this->Dim(); k++, ptr_k += k )
538  {
539  // Faz sum = SOMA( A(k,p) * A(k,p) ), p = 1, ..., k-1.
540  //
541  TVar sum = 0.0;
542  TVar *pk = ptr_k;
543  TVar *pkk = ptr_k + k;
544  for ( ; pk < pkk; pk++ )
545  sum += (*pk) * (*pk);
546 
547  // Faz A(k,k) = sqrt( A(k,k) - sum ).
548  //
549  if ( (*pk -= sum) < 1.e-10 )
550  return( 0 );
551  *pk = sqrt( *pk );
552 
553  // Loop para i = k+1 ... Dim().
554  //
555  TVar *ptr_i = ptr_k;
556  TVar *pi;
557  for ( int64_t i = k+1; i <this->Dim(); i++ )
558  {
559  // Faz sum = SOMA( A(i,p) * A(k,p) ), p = 1, ..., k-1.
560  //
561  sum = 0.0;
562  ptr_i += i;
563  pk = ptr_k;
564  pi = ptr_i;
565  for ( int64_t p = 0; p < k; p++ )
566  sum += (*pk++) * (*pi++);
567 
568  // Faz A(i,k) = (A(i,k) - sum) / A(k,k)
569  //
570  *pi = (*pi - sum) / *pk;
571  }
572  }
573 
574  this->fDecomposed = ECholesky;
575  this->fDefPositive = 1;
576  return( 1 );
577 }
578 
579 
580 
581 /**********************/
582 /*** Decompose LDLt ***/
583 template<class TVar>
584 int
585 TPZSFMatrix<TVar> ::Decompose_LDLt(std::list<int64_t> &singular)
586 {
587  return Decompose_LDLt();
588 }
589 
590 template<class TVar>
591 int
592 TPZSFMatrix<TVar> ::Decompose_LDLt()
593 {
594  if ( this->fDecomposed && this->fDecomposed != ELDLt)
595  TPZMatrix<TVar> ::Error(__PRETTY_FUNCTION__,"Decompose_LDLt <Matrix already Decomposed with a different scheme>" );
596  else if ( this->fDecomposed ) return 0;
597 
598  int64_t j,k,l;
599  // TVar sum;
600 
601 
602  for ( j = 0; j < this->Dim(); j++ )
603  {
604  // sum=0.;
605  for ( k=0; k<j; k++)
606  // sum=sum-GetVal(k,k)*GetVal(k,j)*GetVal(k,j);
607  PutVal( j,j,GetVal(j,j)-GetVal(k,k)*GetVal(k,j)*GetVal(k,j) );
608  // PutVal(j,j,GetVal(j,j)+sum);
609  for ( k=0; k<j; k++)
610  {
611 
612  for( l=j+1; l<this->Dim();l++)
613  PutVal(l,j, GetVal(l,j)-GetVal(k,k)*GetVal(j,k)*GetVal(l,k) );
614  }
615 
616  if ( IsZero(GetVal(j,j)) )TPZMatrix<TVar> ::Error(__PRETTY_FUNCTION__,"Decompose_LDLt <Zero on diagonal>" );
617 
618  for( l=j+1; l<this->Dim();l++) PutVal( l,j,GetVal(l,j)/GetVal(j,j) ) ;
619 
620  }
621  this->fDecomposed = ELDLt;
622  this->fDefPositive = 0;
623 
624 
625  return( 1 );
626 }
627 
628 
629 
630 /*********************/
631 /*** Subst Forward ***/
632 
633 template<class TVar>
634 int
635 TPZSFMatrix<TVar> ::Subst_Forward( TPZFMatrix<TVar> *B ) const
636 {
637  if ( (B->Rows() != this->Dim()) || !this->fDecomposed )
638  return( 0 );
639 
640  if ( B->IsSimetric() )
641  TPZMatrix<TVar> ::Error(__PRETTY_FUNCTION__, "Subst_Forward <the matrix result can not be simetric>" );
642 
643  TVar *ptr_k = fElem;
644  for ( int64_t k = 0; k < this->Dim(); k++ )
645  {
646  for ( int64_t j = 0; j < B->Cols(); j++ )
647  {
648  // Faz sum = SOMA( A[k,i] * B[i,j] ), para i = 1,.., k-1.
649  //
650  TVar *pk = ptr_k;
651  TVar sum = 0.0;
652  for ( int64_t i = 0; i < k; i++ )
653  sum += (*pk++) * B->GetVal( i, j );
654 
655  // Faz B[k,j] = (B[k,j] - sum) / A[k,k].
656  //
657  B->PutVal( k, j, (B->GetVal( k, j ) - sum) / *pk );
658  }
659  ptr_k += k + 1;
660  }
661 
662  return( 1 );
663 }
664 
665 
666 
667 /**********************/
668 /*** Subst Backward ***/
669 template<class TVar>
670 int
671 TPZSFMatrix<TVar> ::Subst_Backward( TPZFMatrix<TVar> *B ) const
672 {
673  if ( (B->Rows() != this->Dim()) || !this->fDecomposed )
674  return( 0 );
675 
676  if ( B->IsSimetric() )
677  TPZMatrix<TVar> ::Error(__PRETTY_FUNCTION__, "Subst_Backward <the matrix result can not be simetric>" );
678 
679  TVar *ptr_k = &fElem[ Size()-1 ];
680  for ( int64_t k = this->Dim()-1; k >= 0; k--, ptr_k-- )
681  for ( int64_t j = 0; j < B->Cols(); j++ )
682  {
683  // Faz sum = SOMA( A[k,i] * B[i,j] ); i = N, ..., k+1.
684  //
685  TVar sum = 0.0;
686  TVar *pk = ptr_k;
687  for ( int64_t i = this->Dim()-1; i > k; i-- )
688  {
689  sum += *pk * B->GetVal( i, j );
690  pk -= i;
691  }
692 
693  // Faz B[k,j] = (B[k,j] - sum) / A[k,k].
694  //
695  B->PutVal( k, j, (B->GetVal(k, j) - sum) / *pk );
696  }
697 
698  return( 1 );
699 }
700 
701 
702 
703 /***********************/
704 /*** Subst L Forward ***/
705 template<class TVar>
706 int
707 TPZSFMatrix<TVar> ::Subst_LForward( TPZFMatrix<TVar> *B ) const
708 {
709  if ( (B->Rows() != this->Dim()) || !this->fDecomposed )
710  return( 0 );
711 
712  if ( B->IsSimetric() )
713  TPZMatrix<TVar> ::Error(__PRETTY_FUNCTION__, "Subst_LForward <the matrix result can not be simetric>" );
714 
715  TVar *ptr_k = fElem;
716  for ( int64_t k = 0; k < this->Dim(); k++, ptr_k += k )
717  for ( int64_t j = 0; j < B->Cols(); j++ )
718  {
719  // Faz sum = SOMA( A[k,i] * B[i,j] ), para i = 1, ..., k-1.
720  //
721  TVar *pk = ptr_k;
722  TVar sum = 0.0;
723  for ( int64_t i = 0; i < k; i++ )
724  sum += (*pk++) * B->GetVal( i, j );
725 
726  // Faz b[k] = (b[k] - sum).
727  //
728  B->PutVal( k, j, B->GetVal( k, j ) - sum );
729  }
730 
731  return( 1 );
732 }
733 
734 
735 
736 /************************/
737 /*** Subst L Backward ***/
738 
739 template<class TVar>
740 int
741 TPZSFMatrix<TVar> ::Subst_LBackward( TPZFMatrix<TVar> *B ) const
742 {
743  if ( (B->Rows() != this->Dim()) || !this->fDecomposed )
744  return( 0 );
745 
746  if ( B->IsSimetric() )
747  TPZMatrix<TVar> ::Error(__PRETTY_FUNCTION__, "Subst_LBackward <the matrix result can not be simetric>" );
748 
749  TVar *ptr_k = &fElem[ Size()-1 ];
750  for ( int64_t k = this->Dim()-1; k >= 0; k--, ptr_k-- )
751  for ( int64_t j = 0; j < B->Cols(); j++ )
752  {
753  // Faz sum = SOMA( A[k,i] * B[i,j] ); i = N, ..., k+1.
754  //
755  TVar sum = 0.0;
756  TVar *pk = ptr_k;
757  for ( int64_t i = this->Dim()-1; i > k; i-- )
758  {
759  sum += *pk * B->GetVal( i, j );
760  pk -= i;
761  }
762 
763  // Faz B[k,j] = B[k,j] - sum.
764  //
765  B->PutVal( k, j, B->GetVal(k, j) - sum );
766  }
767 
768  return( 1 );
769 }
770 
771 
772 
773 /******************/
774 /*** Subst Diag ***/
775 
776 template<class TVar>
777 int
778 TPZSFMatrix<TVar> ::Subst_Diag( TPZFMatrix<TVar> *B ) const
779 {
780  if ( (B->Rows() != this->Dim()) || !this->fDecomposed )
781  return( 0 );
782 
783  TVar *pDiag = fElem;
784  for ( int64_t k = 0; k < this->Dim(); k++ )
785  {
786  for ( int64_t j = 0; j < B->Cols(); j++ )
787  B->PutVal( k, j, B->GetVal( k, j ) / *pDiag );
788  pDiag += (k + 2);
789  }
790 
791  return( 1 );
792 }
793 
794 
795 
796 
797 /************************** Private **************************/
798 
799 /*************/
800 /*** Error ***/
801 
802 /*int
803  TPZSFMatrix::Error(const char *msg1,const char *msg2 )
804  {
805  ostringstream out;
806  out << "TPZSFMatrix::" << msg1 << msg2 << ".\n";
807  //pzerror.Show();
808  LOGPZ_ERROR (logger, out.str().c_str());
809  DebugStop();
810  return 0;
811  }*/
812 
813 
814 
815 /*************/
816 /*** Clear ***/
817 
818 template<class TVar>
819 int
820 TPZSFMatrix<TVar> ::Clear()
821 {
822  if ( fElem != NULL )
823  delete( fElem );
824 
825  fElem = NULL;
826  this->fRow = this->fCol = 0;
827  return( 1 );
828 }
829 
830 template class TPZSFMatrix<float>;
831 template class TPZSFMatrix<double>;
832 template class TPZSFMatrix<long double>;
833 
834 
835 #ifdef OOPARLIB
836 
837 template<class TVar>
838 int TPZSFMatrix<TVar> ::Unpack( TReceiveStorage *buf ){
839  TSaveable::Unpack(buf);
840  buf->UpkDouble(fElem);
841  return 1;
842 }
843 
844 
845 template<class TVar>
846 TSaveable *TPZSFMatrix<TVar>::CreateInstance(TReceiveStorage *buf) {
847  TPZSFMatrix<TVar> *m = new TPZSFMatrix<TVar>();
848  m->Unpack(buf);
849  return m;
850 }
851 
852 template<class TVar>
853 int TPZSFMatrix<TVar>::Pack( TSendStorage *buf ){
854  TSaveable::Pack(buf);
855  buf->PkDouble(fElem);
856  return 1;
857 }
858 
859 template<class TVar>
860 int TPZSFMatrix<TVar>::DerivedFrom(int64_t Classid){
861  return TSaveable::DerivedFrom(Classid);
862 }
863 template<class TVar>
864 int TPZSFMatrix<TVar>::DerivedFrom(char *classname){
865 
866  if(!strcmp(ClassName(),classname)) return 1;
867  return TSaveable::DerivedFrom(classname);
868 }
869 
870 #endif
871 
872 template<class TVar>
873 int TPZSFMatrix<TVar>::ClassId() const{
874  return Hash("TPZSFMatrix") ^ TPZMatrix<TVar>::ClassId() << 1;
875 }
TPZSFMatrix::Size
int64_t Size() const
Definition: pzsfulmat.h:142
TPZSFMatrix::Zero
int Zero() override
Resets all elements.
Definition: pzsfulmat.cpp:510
TPZSFMatrix::Subst_Forward
virtual int Subst_Forward(TPZFMatrix< TVar > *B) const override
Computes B = Y, where A*Y = B, A is lower triangular.
Definition: pzsfulmat.cpp:635
ELDLt
Definition: pzmatrix.h:52
IsZero
bool IsZero(long double a)
Returns if the value a is close Zero as the allowable tolerance.
Definition: pzreal.h:668
pzlog.h
Contains definitions to LOGPZ_DEBUG, LOGPZ_INFO, LOGPZ_WARN, LOGPZ_ERROR and LOGPZ_FATAL, and the implementation of the inline InitializePZLOG(string) function using log4cxx library or not. It must to be called out of "#ifdef LOG4CXX" scope.
TPZSFMatrix::operator*
TPZSFMatrix operator*(const TVar val) const
Definition: pzsfulmat.cpp:390
TPZSFMatrix::Clear
int Clear() override
It clears data structure.
Definition: pzsfulmat.cpp:820
TPZMatrix::Error
static int Error(const char *msg, const char *msg2=0)
Returns error messages.
Definition: pzmatrix.cpp:1402
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::SetIsDecomposed
void SetIsDecomposed(int val)
Sets current matrix to decomposed state.
Definition: pzmatrix.h:410
TPZSFMatrix::Decompose_LDLt
virtual int Decompose_LDLt() override
Decomposes the current matrix using LDLt.
Definition: pzsfulmat.cpp:592
TPZMatrix::fDecomposed
char fDecomposed
Decomposition type used to decompose the current matrix.
Definition: pzmatrix.h:783
TPZMatrix::fDefPositive
char fDefPositive
Definite Posistiveness of current matrix.
Definition: pzmatrix.h:785
TPZSFMatrix::Subst_Backward
virtual int Subst_Backward(TPZFMatrix< TVar > *B) const override
Computes B = Y, where A*Y = B, A is upper triangular.
Definition: pzsfulmat.cpp:671
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
TPZSFMatrix::operator-=
TPZSFMatrix & operator-=(const TPZSFMatrix< TVar > &A)
Definition: pzsfulmat.cpp:206
TPZSFMatrix::operator-
TPZSFMatrix operator-() const
Definition: pzsfulmat.h:98
TPZSFMatrix::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: pzsfulmat.h:153
TPZMatrix::IsSimetric
virtual int IsSimetric() const
Checks if the current matrix is symmetric.
Definition: pzmatrix.h:394
TPZMatrix::IsDecomposed
int IsDecomposed() const
Checks if current matrix is already decomposed.
Definition: pzmatrix.h:405
pzfmatrix.h
Contains TPZMatrixclass which implements full matrix (using column major representation).
TPZSFMatrix::Subst_Diag
virtual int Subst_Diag(TPZFMatrix< TVar > *B) const override
Computes B = Y, where A*Y = B, A is diagonal matrix.
Definition: pzsfulmat.cpp:778
TPZSFMatrix::operator*=
TPZSFMatrix & operator*=(const TVar val)
Definition: pzsfulmat.cpp:427
TPZSFMatrix::operator+=
TPZSFMatrix & operator+=(const TPZSFMatrix< TVar > &A)
Definition: pzsfulmat.cpp:184
TPZMatrix::Rows
int64_t Rows() const
Returns number of rows.
Definition: pzmatrix.h:803
TPZSavable::CreateInstance
static TPZSavable * CreateInstance(const int &classId)
Definition: TPZSavable.cpp:120
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
TPZSFMatrix::Subst_LBackward
virtual int Subst_LBackward(TPZFMatrix< TVar > *B) const override
Computes B = Y, where A*Y = B, A is upper triangular with A(i,i)=1.
Definition: pzsfulmat.cpp:741
test.res
string res
Definition: test.py:151
Min
const T & Min(const T &a, const T &b)
Returns the minimum value between a and b.
Definition: pzreal.h:729
TPZFMatrix
Full matrix class. Matrix.
Definition: pzfmatrix.h:32
TPZMatrix::Dim
virtual int64_t Dim() const
Returns the dimension of the matrix if the matrix is square.
Definition: pzmatrix.h:892
Hash
int32_t Hash(std::string str)
Definition: TPZHash.cpp:10
TPZSFMatrix::Decompose_Cholesky
virtual int Decompose_Cholesky() override
Decomposes the current matrix using Cholesky method. The current matrix has to be symmetric...
Definition: pzsfulmat.cpp:533
TPZSFMatrix::operator=
TPZSFMatrix & operator=(const TPZSFMatrix< TVar > &A)
Definition: pzsfulmat.cpp:102
TPZMatrix::fCol
int64_t fCol
Number of cols in matrix.
Definition: pzmatrix.h:781
TPZSFMatrix::Resize
int Resize(const int64_t newDim, const int64_t) override
Resize the array but keeps its entirety.
Definition: pzsfulmat.cpp:443
TPZSFMatrix::fElem
TVar * fElem
Definition: pzsfulmat.h:146
TPZSFMatrix::Subst_LForward
virtual int Subst_LForward(TPZFMatrix< TVar > *B) const override
Computes B = Y, where A*Y = B, A is lower triangular with A(i,i)=1.
Definition: pzsfulmat.cpp:707
TPZSFMatrix::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: pzsfulmat.h:168
TPZMatrix::Cols
int64_t Cols() const
Returns number of cols.
Definition: pzmatrix.h:809
TPZSFMatrix::operator+
TPZSFMatrix operator+(const TPZSFMatrix< TVar > &A) const
Definition: pzsfulmat.cpp:137
TPZSFMatrix
Implements a symmetric full matrix. Matrix.
Definition: pzsfulmat.h:24
TPZMatrix::IsDefPositive
virtual int IsDefPositive() const
Checks if current matrix is definite positive.
Definition: pzmatrix.h:403
m
clarg::argString m("-m", "input matrix file name (text format)", "matrix.txt")
ECholesky
Definition: pzmatrix.h:52
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
TPZSFMatrix::Redim
int Redim(const int64_t newRows, const int64_t) override
Resize the array and resets ist entirety.
Definition: pzsfulmat.cpp:483
pzsfulmat.h
Contains TPZSFMatrix class which implements a symmetric full matrix.
TPZMatrix::ClassId
int ClassId() const override
Define the class id associated with the class.
Definition: pzmatrix.h:957
TPZSFMatrix::ClassId
int ClassId() const override
Define the class id associated with the class.
Definition: pzsfulmat.cpp:873
TPZMatrix
Root matrix class (abstract). Matrix.
Definition: pzmatrix.h:60
Generated on Sun Aug 16 2020 11:17:05 for NeoPZ by   doxygen 1.8.13