NeoPZ
tpzeqnarray.cpp
Go to the documentation of this file.
1 
6 #include <iostream>
7 #include <math.h>
8 #include <stdlib.h>
9 #include <stdio.h>
10 #include <fstream>
11 
12 #include "tpzeqnarray.h"
13 #include "pzreal.h"
14 
15 using namespace std;
16 
17 template<class TVar>
18 void TPZEqnArray<TVar>::SetNonSymmetric(){
19  fSymmetric=EIsNonSymmetric;
20 }
21 template<class TVar>
22 int TPZEqnArray<TVar>::IsSymmetric(){
23  return fSymmetric;
24 }
25 template<class TVar>
26 void TPZEqnArray<TVar>::SetSymmetric(){
27  fSymmetric = EIsSymmetric;
28 }
29 template<class TVar>
30 TPZEqnArray<TVar>::~TPZEqnArray(){
31 
32 }
33 template<class TVar>
34 TPZEqnArray<TVar>::TPZEqnArray() : fEqStart(), fEqNumber(), fEqValues(), fIndex() {
35  fEqStart.Push(0);
36  fNumEq=0;
37  fLastTerm=0;
38  fSymmetric=EIsUndefined;
39 
40 }
41 template<class TVar>
42 void TPZEqnArray<TVar>::Print(const char *name, std::ostream& out)
43 {
44  if(name) out << name << endl;
45  int i,j;
46  for(i=0;i<fNumEq;i++){
47  int aux_limit;
48 
49  if(i==fNumEq-1){
50  aux_limit=fEqValues.NElements();
51  }else aux_limit=fEqStart[i+1];
52  for(j=fEqStart[i]; j<aux_limit ; j++) {
53  out << "col = " << fIndex[j] << " ";
54  }
55  out << endl;
56  for(j=fEqStart[i]; j<aux_limit ; j++) {
57  out << fEqValues[j] << " ";
58  }
59  out << endl;
60  out << endl;
61  }
62 }
63 
64 template<class TVar>
65 void TPZEqnArray<TVar>::Write(char * outputfile) { }
66 template<class TVar>
67 void TPZEqnArray<TVar>::Read(char * inputfile) { }
68 
69 template<class TVar>
70 void TPZEqnArray<TVar>::EndEquation() {
71  fEqStart.Push(fLastTerm);
72 }
73 
74 template<class TVar>
75 void TPZEqnArray<TVar>::BeginEquation(int eq) {
76  fEqNumber.Push(eq);
77  fNumEq++;
78 }
79 
80 template<class TVar>
81 void TPZEqnArray<TVar>::Reset(){
82  fEqStart.Resize(0);
83  fEqStart.Push(0);
84  fEqValues.Resize(0);
85  fIndex.Resize(0);
86  fNumEq=0;
87  fLastTerm=0;
88  fSymmetric=EIsUndefined;
89 
90 }
91 
92 template<class TVar>
93 void TPZEqnArray<TVar>::EqnBackward(TPZFMatrix<TVar> & U, DecomposeType dec) {
94  int n;
95  if(IsSymmetric()==EIsSymmetric){
96  for(n=fNumEq-1; n>=0; n--) {
97  int index = fEqStart[n];
98  int last_term = fEqStart[n + 1];
99  TVar acum = 0.;
100  int i;
101  for(i = index + 1; i < last_term; i++) acum -= U(fIndex[i],0) * fEqValues[i];
102 
104  U(fIndex[index],0) += acum;
105 
107  if(dec == ECholesky){
108  U(fIndex[index],0) /= fEqValues[index];
109  }
110  }
111  }else if(IsSymmetric()==EIsNonSymmetric){
112  for(n=fNumEq-2; n>=0; n-=2) {
113  int index = fEqStart[n];
114  int last_term = fEqStart[n + 1];
115  TVar acum = 0.;
116  int i;
117  for(i = index + 1; i < last_term; i++) acum -= U(fIndex[i],0) * fEqValues[i];
118 
120  U(fIndex[index],0) += acum;
121 
123  if(dec == ECholesky || dec == ELU){
124  U(fIndex[index],0) /= fEqValues[index];
125  }
126  }
127  }
128 }
129 
130 template<class TVar>
131 void TPZEqnArray<TVar>::EqnForward(TPZFMatrix<TVar> & F, DecomposeType dec) {
132  int j;
133  if(IsSymmetric()==EIsSymmetric){
134 
135  for(j=0; j<fNumEq; j++) {
136  int index = fEqStart[j];
137  if(fEqValues[index] == (TVar)0.) {
138  cout << "Diagonal Value = 0 >> Aborting on Index = " << index << " Equation = " << j << endl;
139  DebugStop();
140  }
141  int last_term;
142  if(j==fNumEq-1){
143  last_term=fEqValues.NElements();
144  }else last_term = fEqStart[j+1];
145 
147  if(dec==ECholesky || dec==ELU) {
148  F(fIndex[index],0) /= fEqValues[index];
149  }
151  TVar udiag = F(fIndex[index],0);
152 
153  int i;
154  //+1 ou +2
155  for(i = index + 1; i < last_term; i++) F(fIndex[i],0) -= udiag * fEqValues[i];
157  if(dec == ELDLt) F(fIndex[index],0) /= fEqValues[index];
158  }
159  } else if(IsSymmetric()==EIsNonSymmetric) {
160  for(j=1; j<fNumEq; j+=2) {
161  int index = fEqStart[j];
162 
163  if(fEqValues[index] == (TVar)0.){
164  cout << "Diagonal Value = 0 >> Aborting on Index = " << index << " Equation = " << j << endl;
165  DebugStop();
166  }
167 
168  int last_term;
169  if(j==fNumEq-1){
170  last_term=fEqValues.NElements();
171  }else last_term = fEqStart[j+1];
172  //if(fEqStart.NElements()
173 
174 
176  if(dec==ECholesky || dec==ELU) {
177  F(fIndex[index],0) /= fEqValues[index];
178  }
180  TVar udiag = F(fIndex[index],0);
181 
182 
183  int i;
184  for(i = index + 1; i < last_term; i++) F(fIndex[i],0) -= udiag * fEqValues[i];
185  }
186  }
187 }
188 
189 template<class TVar>
190 void TPZEqnArray<TVar>::Write(FILE * outputfile){
192  fwrite(&fNumEq,sizeof(int),1,outputfile);
193  //cout << ftell(outputfile) << endl;
195  fwrite(&fLastTerm,sizeof(int),1,outputfile);
196  //cout << ftell(outputfile) << endl;
198  int aux = fEqStart.NElements();
199  if(aux == 0)
200  {
201  std::cout << __PRETTY_FUNCTION__ << __LINE__ << std::endl;
202  }
203  fwrite(&aux,sizeof(int),1,outputfile);
204  fwrite(&fEqStart[0],sizeof(int), aux ,outputfile);
205 
207  aux = fEqNumber.NElements();
208  if(aux == 0)
209  {
210  std::cout << __PRETTY_FUNCTION__ << __LINE__ << std::endl;
211  }
212  fwrite(&aux,sizeof(int),1,outputfile);
213  fwrite(&fEqNumber[0],sizeof(int), aux ,outputfile);
214 
216  aux = fIndex.NElements();
217  fwrite(&aux,sizeof(int),1,outputfile);
218  fwrite(&fIndex[0],sizeof(int), aux ,outputfile);
219 
221  aux = fEqValues.NElements();
222  fwrite(&aux,sizeof(int),1,outputfile);
223  fwrite(&fEqValues[0],sizeof(REAL), aux ,outputfile);
224 }
225 
226 template<class TVar>
227 void TPZEqnArray<TVar>::Read(FILE * inputfile) {
228  int64_t sizereturn;
229  sizereturn = 0;
231  sizereturn = fread(&fNumEq,sizeof(int),1,inputfile);
232 #ifdef PZDEBUG
233  if (sizereturn != 1) DebugStop();
234 #endif
235 
236  sizereturn = fread(&fLastTerm,sizeof(int),1,inputfile);
237 #ifdef PZDEBUG
238  if (sizereturn != 1) DebugStop();
239 #endif
240  int aux=0;
242  sizereturn = fread(&aux,sizeof(int),1,inputfile);
243 #ifdef PZDEBUG
244  if (sizereturn != 1) DebugStop();
245 #endif
246  fEqStart.Resize(aux);
247  sizereturn = fread(&fEqStart[0],sizeof(int), aux ,inputfile);
248 #ifdef PZDEBUG
249  if (sizereturn != aux) DebugStop();
250 #endif
251 
252  sizereturn = fread(&aux,sizeof(int),1,inputfile);
253 #ifdef PZDEBUG
254  if (sizereturn != 1) DebugStop();
255 #endif
256  fEqNumber.Resize(aux);
257  sizereturn = fread(&fEqNumber[0],sizeof(int), aux ,inputfile);
258 #ifdef PZDEBUG
259  if (sizereturn != aux) DebugStop();
260 #endif
261 
262  sizereturn = fread(&aux,sizeof(int),1,inputfile);
263 #ifdef PZDEBUG
264  if (sizereturn != 1) DebugStop();
265 #endif
266  fIndex.Resize(aux);
267  sizereturn = fread(&fIndex[0],sizeof(int), aux ,inputfile);
268 #ifdef PZDEBUG
269  if (sizereturn != aux) DebugStop();
270 #endif
271 
272  sizereturn = fread(&aux,sizeof(int),1,inputfile);
273 #ifdef PZDEBUG
274  if (sizereturn != 1) DebugStop();
275 #endif
276  fEqValues.Resize(aux);
277  sizereturn = fread(&fEqValues[0],sizeof(REAL), aux ,inputfile);
278 #ifdef PZDEBUG
279  if (sizereturn != aux) DebugStop();
280 #endif
281 
282  this->fSymmetric = EIsSymmetric;
283  if(fNumEq && fNumEq%2==0 && fEqNumber[0]==fEqNumber[1]) fSymmetric = EIsNonSymmetric;
284 }
285 
286 template<class TVar>
287 void TPZEqnArray<TVar>::main()
288 {
289  char filename[20];
290 
291  cout << "Entre o nome do Arquivo\n";
292  cin >> filename;
293  ofstream output(filename,ios::app);
294 
295  TPZFMatrix<TVar> MatrixA(10,10);
296  int i, j;
297  for(i=0;i<10;i++) {
298  for(j=i;j<10;j++) {
299  int random = rand();
300  double rnd = (random*10.)/RAND_MAX;
301  MatrixA(i,j)=rnd;
302  MatrixA(j,i)=MatrixA(i,j);
303  if(i==j) MatrixA(i,j)=6000.;
304  }
305  }
306 
307  MatrixA.Print("Teste 1",std::cout);
308  MatrixA.Decompose_Cholesky();
309  MatrixA.Print("Decomposta",std::cout);
310 
311  TPZEqnArray<TVar> Test;
312  Test.Reset();
313  for(i=0;i<10;i++) {
314  Test.BeginEquation(i);
315  for(j=i;j<10;j++) {
316  Test.AddTerm(j, MatrixA(i,j));
317  }
318  Test.EndEquation();
319  }
320 
321  TPZFMatrix<TVar> rhs(10,1), rhs2(10,1);
322  //Inicializar rhs:
323  rhs.Zero();
324  rhs(0,0) = 1.;
325  rhs2=rhs;
326 
327  Test.Print("Result", cout);
328 
329  MatrixA.Subst_Forward(&rhs);
330  DecomposeType decType = ECholesky;
331  Test.EqnForward(rhs2, decType);
332 
333  rhs.Print("FMatrix Decomposition");
334  rhs2.Print("FrontalMatrix Decomposition");
335 
336 }
337 
338 
339 template class TPZEqnArray<float>;
340 template class TPZEqnArray<double>;
341 template class TPZEqnArray<long double>;
342 
343 template class TPZEqnArray<std::complex<float> >;
344 template class TPZEqnArray<std::complex<double> >;
345 template class TPZEqnArray<std::complex<long double> >;
TPZEqnArray::fEqNumber
TPZStack< int, 100 > fEqNumber
Definition: tpzeqnarray.h:126
TPZEqnArray::EqnForward
void EqnForward(TPZFMatrix< TVar > &F, DecomposeType dec)
Forward substitution on equations stored in EqnArray.
Definition: tpzeqnarray.cpp:131
stats.filename
filename
Definition: stats.py:82
TPZEqnArray::fEqStart
TPZStack< int, 100 > fEqStart
Equation start point index.
Definition: tpzeqnarray.h:125
ELDLt
Definition: pzmatrix.h:52
TPZEqnArray::fLastTerm
int fLastTerm
Indicates the last used position in fEqValues.
Definition: tpzeqnarray.h:135
TPZEqnArray
It is an equation array, generally in its decomposed form. Frontal.
Definition: tpzeqnarray.h:36
TPZEqnArray::main
static void main()
Static main function for testing.
Definition: tpzeqnarray.cpp:287
TPZEqnArray::Write
void Write(char *outputfile)
Writes on disk.
Definition: tpzeqnarray.cpp:65
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
TPZEqnArray::BeginEquation
void BeginEquation(int eq)
It starts an equation storage.
Definition: tpzeqnarray.cpp:75
TPZManVector::Resize
virtual void Resize(const int64_t newsize, const T &object)
Resizes the vector object.
Definition: pzmanvector.h:426
tpzeqnarray.h
Contains the TPZEqnArray class which implements an equation array.
TPZEqnArray::Print
void Print(const char *name, std::ostream &out)
It prints all terms stored in TPZEqnArray.
Definition: tpzeqnarray.cpp:42
TPZEqnArray::Read
void Read(char *inputfile)
Reads from disk.
Definition: tpzeqnarray.cpp:67
TPZStack::Push
void Push(const T object)
Pushes a copy of the object on the stack.
Definition: pzstack.h:80
TPZEqnArray::fNumEq
int fNumEq
Number of equations.
Definition: tpzeqnarray.h:122
DebugStop
#define DebugStop()
Returns a message to user put a breakpoint in.
Definition: pzerror.h:20
TPZEqnArray::EqnBackward
void EqnBackward(TPZFMatrix< TVar > &U, DecomposeType dec)
Backward substitution on equations stored in EqnArray.
Definition: tpzeqnarray.cpp:93
ELU
Definition: pzmatrix.h:52
TPZEqnArray::Reset
void Reset()
Resets data structure.
Definition: tpzeqnarray.cpp:81
TPZEqnArray::SetNonSymmetric
void SetNonSymmetric()
Sets EqnArray to a non symmetric form.
Definition: tpzeqnarray.cpp:18
TPZEqnArray::SetSymmetric
void SetSymmetric()
Sets fSymmetric to the symmetric value.
Definition: tpzeqnarray.cpp:26
TPZEqnArray::AddTerm
void AddTerm(int col, TVar val)
Add a term to the current equation.
Definition: tpzeqnarray.h:83
TPZFMatrix
Full matrix class. Matrix.
Definition: pzfmatrix.h:32
TPZEqnArray::fSymmetric
MSymmetric fSymmetric
Indicates the symetry or not of the equationarray.
Definition: tpzeqnarray.h:118
TPZEqnArray::EndEquation
void EndEquation()
Ends the current equation.
Definition: tpzeqnarray.cpp:70
TPZEqnArray::fEqValues
TPZStack< TVar, 1000 > fEqValues
Equations coefficients values.
Definition: tpzeqnarray.h:129
TPZEqnArray::TPZEqnArray
TPZEqnArray()
Simple constructor.
Definition: tpzeqnarray.cpp:34
TPZEqnArray::~TPZEqnArray
~TPZEqnArray()
Simple desctructor.
Definition: tpzeqnarray.cpp:30
TPZMatrix::Print
virtual void Print(std::ostream &out) const
Definition: pzmatrix.h:253
TPZFMatrix::Decompose_Cholesky
virtual int Decompose_Cholesky() override
Cholesky Decomposition Optmized. for walks in the direction of the vector that composes the matrix...
Definition: pzfmatrix.cpp:1574
TPZVec::NElements
int64_t NElements() const
Returns the number of elements of the vector.
Definition: pzvec.h:190
pzreal.h
Contains the declaration of TPZFlopCounter class and TPZCounter struct.
ECholesky
Definition: pzmatrix.h:52
TPZEqnArray::fIndex
TPZStack< int, 1000 > fIndex
Line/Column number associated to each fEqValues values.
Definition: tpzeqnarray.h:132
TPZEqnArray::IsSymmetric
int IsSymmetric()
Gets the symetry situation of EqnArray.
Definition: tpzeqnarray.cpp:22
DecomposeType
DecomposeType
Defines decomposition type for any matrix classes.
Definition: pzmatrix.h:52
Generated on Sun Aug 16 2020 11:17:01 for NeoPZ by   doxygen 1.8.13