NeoPZ
pzstencil.cpp
Go to the documentation of this file.
1 
6 #include "pzstencil.h"
7 #include "pzfmatrix.h"
8 
9 #include <memory.h>
10 #include <string>
11 using namespace std;
12 
13 // ****************************************************************************
14 //
15 // Constructors and the destructor
16 //
17 // ****************************************************************************
18 template<class TVar>
19 TPZStencilMatrix<TVar>::TPZStencilMatrix( int rows, int cols ) {
20  // Constructs an empty TPZStencilMatrix
21  fRows = rows;
22  fCols = cols;
23  fMystencils = 0;
24  fNumberOfStencilPointers = 0;
25  fStencilNumbers = 0;
26  fSolver = -1;
27  fMaxIterations = 4;
28  fSORRelaxation = 1.;
29  fDiag = 0;
30 #ifdef CONSTRUCTOR
31  cerr << "TPZStencilMatrix(int rows,int cols)\n";
32 #endif
33 }
34 
35 template<class TVar>
36 TPZStencilMatrix<TVar>::~TPZStencilMatrix() {
37  // Deletes everything associated with a TPZStencilMatrix
38  for(int i=0; i<fNumberOfStencilPointers; i++) delete fMystencils[i];
39  delete fMystencils;
40  fMystencils = 0;
41  fNumberOfStencilPointers = 0;
42  delete fStencilNumbers;
43  fStencilNumbers = 0;
44  delete fDiag;
45  fDiag = 0;
46  fSolver = -1;
47 #ifdef DESTRUCTOR
48  cerr << "~TPZStencilMatrix()\n";
49 #endif
50 }
51 
52 // ****************************************************************************
53 //
54 // Find the element of the matrix at (row,col) in the stencil matrix
55 //
56 // ****************************************************************************
57 
58 template<class TVar>
59 const TVar & TPZStencilMatrix<TVar>::GetVal(const int row,const int col ) const {
60  // Get the matrix entry at (row,col) without bound checking
61  int stenindex;
62  MPStencil *st;
63 
64  // Walk past each row of the matrix trying to determine where the
65  // first stencil in row begins
66  static TPZStencilMatrix const *whatsthis = 0;
67  static int ix=0;
68  static int ir = 0;
69  if( whatsthis != this || row < ir ) {
70  whatsthis = this;
71  for( ir=ix=0; ir<row; ir++) {
72  stenindex = fStencilNumbers[ir];
73  st = fMystencils[stenindex];
74  ix += st->fInc;
75  }
76  }
77  else {
78  for( ; ir<row; ir++) {
79  stenindex = fStencilNumbers[ir];
80  st = fMystencils[stenindex];
81  ix += st->fInc;
82  }
83  }
84 
85  // Now look through the requested row and see if there is anything
86  // in column col
87  st = fMystencils[fStencilNumbers[row]];
88  int nitems = st->fNumberOfItems;
89  int *ia = st->fIA;
90  TVar *a = st->fA;
91  int it=0;
92 
93  while(++it<nitems) {
94  int numintegers = *ia++;
95  int *ialast = ia+numintegers;
96  while(ia < ialast) {
97  if ( ix+(*ia++) == col ) return *a;
98  it++;
99  }
100  a+=numintegers+1;
101  }
102  return this->gZero;
103 }
104 
105 // ****************************************************************************
106 //
107 // Multiply and Multiply-Add
108 //
109 // ****************************************************************************
110 
111 template<class TVar>
112 void TPZStencilMatrix<TVar>::MultAdd(const TPZFMatrix<TVar> &x,const TPZFMatrix<TVar> &y,
113  TPZFMatrix<TVar> &z,
114  const TVar alpha,const TVar beta,const int opt) const {
115  // computes z = beta * y + alpha * opt(this)*x
116  // z and x cannot share storage
117  int ix=0;
118  int r = (opt) ? Cols() : Rows();
119  if(beta != 0) {
120  TVar *zp = &(z(0,0)), *zlast = zp+r;
121  const TVar *yp = &(y.GetVal(0,0));
122 
123  if(&z != &y) {
124  memcpy(zp,yp,r*sizeof(REAL));
125  }
126  } else {
127  TVar *zp = &(z(0,0)), *zlast = zp+r;
128  while(zp != zlast) {
129  *zp = 0.;
130  zp ++;
131  }
132  }
133  if(opt == 0) {
134 
135  for(int ir=0; ir<fRows; ir++) {
136  int stenindex = fStencilNumbers[ir];
137  MPStencil *st = fMystencils[stenindex];
138  int nitems = st->fNumberOfItems;
139  int *ia = st->fIA;
140  TVar *a = st->fA;
141  int it=0;
142  while(++it<nitems) {
143  int numintegers = *ia++;
144  int *ialast = ia+numintegers;
145  TVar val;
146  val = 0.;
147  const TVar *xp = &(x.GetVal(ix,0));
148  while(ia < ialast) {
149  val += *(xp+((*ia++)));
150  it++;
151  }
152  z(ir,0) += alpha*val*(*a);
153  a+=numintegers+1;
154  }
155  ix += st->fInc;
156  }
157  } else {
158 
159  int ir;
160  for(ir=0; ir<fRows; ir++) {
161  MPStencil *st = fMystencils[fStencilNumbers[ir]];
162  int nitems = st->fNumberOfItems;
163  int *ia = st->fIA;
164  TVar *a = st->fA;
165  int it=0;
166  while(++it<nitems) {
167  int numintegers = *ia++;
168  int *ialast = ia + numintegers;
169  TVar xval = alpha*x.GetVal(ir,0)*(*a);
170  while(ia<ialast) {
171  z((ix+(*ia++)),0) += xval;
172  it++;
173  }
174  a += numintegers+1;
175  }
176  ix += st->fInc;
177  }
178  }
179 }
180 
181 // ****************************************************************************
182 //
183 // Print the matrix
184 //
185 // ****************************************************************************
186 
187 template<class TVar>
188 void TPZStencilMatrix<TVar>::Print(const char *title, ostream &out,const MatrixOutputFormat form ) const {
189  // Print the matrix along with a identification title
190  if(form != EInputFormat) {
191  out << "\nTStencilMatrix Print: " << title << '\n'
192  << "\tRows = " << fRows << '\n'
193  << "\tColumns = " << fCols << '\n';
194  out << '\t';
195  for(int ir=0; ir<fNumberOfStencilPointers; ir++) {
196  MPStencil *st = fMystencils[ir];
197  if (st != 0) {
198  int nitems = st->fNumberOfItems;
199  int *ia = st->fIA;
200  TVar *a = st->fA;
201  out << "\nStencil " << ir
202  << ", increment " << st->fInc << ", "
203  << nitems << " items\n";
204  out << "\tIA\tA\n";
205  out << "\t--\t-\n";
206  int it=0;
207  while(++it<nitems) {
208  int numintegers = *ia;
209  out << '\t' << *ia++ << '\t' << *a << '\n';
210  int *ialast = ia+numintegers;
211  while(ia < ialast) {
212  out << '\t' << *ia++ << '\n';
213  it++;
214  }
215  a+=numintegers+1;
216  }
217  }
218  }
219 
220  if(fStencilNumbers) {
221  out << "\nIndex\tStencil number\n";
222  out << "-----\t--------------\n";
223  for(int i=0; i<fRows; i++) {
224  out << i << '\t' << fStencilNumbers[i] << '\n';
225  }
226  }
227  }
228 }
229 
230 
231 // ****************************************************************************
232 //
233 // Various solvers
234 //
235 // ****************************************************************************
236 
237 template<class TVar>
238 void TPZStencilMatrix<TVar>::ComputeDiagonal() {
239 
240  if(fDiag) return;
241  fDiag = new REAL[Rows()];
242  for(int ir=0; ir<fRows; ir++) {
243  MPStencil *st = fMystencils[fStencilNumbers[ir]];
244  int nitems = st->fNumberOfItems;
245  int *ia = st->fIA;
246  TVar *a = st->fA;
247  if(st->fInc != 1) {
248  cout << "TPZStencilMatrix::ComputeDiagonal "
249  "Computing the diagonal of a nonsquare matrix?";
250  }
251  int it=0;
252  while(++it<nitems) {
253  int numintegers = *ia++;
254  int *ialast = ia+numintegers;
255  while(ia < ialast) {
256  if(*ia++ == 0) {
257  fDiag[ir] = *a;
258  it=nitems;
259  break;
260  }
261  it++;
262  }
263  a+=numintegers+1;
264  }
265  }
266 }
267 
268 template<class TVar>
269 void TPZStencilMatrix<TVar>::SolveSOR( int &numiterations,const TPZFMatrix<TVar> &rhs, TPZFMatrix<TVar> &x,
270  TPZFMatrix<TVar> *residual, TPZFMatrix<TVar>&/*scratch*/,
271  const TVar overrelax, TVar &tol,
272  const int FromCurrent, const int direction ) {
273 
274  if(!fDiag) ComputeDiagonal();
275  int irStart = 0,irLast = fRows,irInc = 1;
276  if(direction < 0) {
277  irStart = fRows-1;
278  irLast = -1;
279  irInc = -1;
280  }
281  if(!FromCurrent) x.Zero();
282  TVar eqres = 2.*tol;
283  int iteration;
284  for(iteration=0; iteration<numiterations && eqres >= tol; iteration++) {
285  eqres = 0.;
286  int ir=irStart;
287  while(ir != irLast) {
288  int stenindex = fStencilNumbers[ir];
289  MPStencil *st = fMystencils[stenindex];
290  int nitems = st->fNumberOfItems;
291  int *ia = st->fIA;
292  TVar *a = st->fA;
293  int it=0;
294  TVar xnewval=rhs.GetVal(ir,0);
295  while(++it<nitems) {
296  int numintegers = *ia++;
297  int *ialast = ia+numintegers;
298  TVar val =0.;
299  TVar *xp = &x(ir,0);
300  while(ia < ialast) {
301  val -= *(xp+(*ia++));
302  it++;
303  }
304  xnewval += *a * val;
305  a+=numintegers+1;
306  }
307  eqres += xnewval*xnewval;
308  x(ir,0) += overrelax*(xnewval/fDiag[ir]);
309  ir += irInc;
310  }
311  eqres = sqrt(eqres);
312  }
313  tol = eqres;
314  numiterations = iteration;
315  if(residual) Residual(x,rhs,*residual);
316 }
317 
318 
319 // ****************************************************************************
320 //
321 // Construct pieces of the stencil matrices
322 //
323 // ****************************************************************************
324 
325 template<class TVar>
326 void TPZStencilMatrix<TVar>::SetStencil( int stencilnumber, int inc,
327  int *IA, TVar *A ) {
328  // initiates Stencil number "stencilnumber" with the data
329 
330  if(stencilnumber < 0) return;
331  if(stencilnumber >= fNumberOfStencilPointers)
332  IncreaseStencilPointers(stencilnumber);
333  if(fMystencils[stencilnumber]) delete fMystencils[stencilnumber];
334  fMystencils[stencilnumber] = new MPStencil(inc,IA,A);
335 }
336 
337 template<class TVar>
338 void TPZStencilMatrix<TVar>::IncreaseStencilPointers( int numsten ) {
339  int newnum = fNumberOfStencilPointers+10;
340  if(newnum < numsten) newnum = numsten+1;
341  MPStencil **newptr = new MPStencil*[newnum];
342  memcpy(newptr,fMystencils,sizeof(MPStencil *)*fNumberOfStencilPointers);
343  int i=fNumberOfStencilPointers;
344  while(i< newnum) newptr[i++] = 0;
345  delete fMystencils;
346  fMystencils = newptr;
347  fNumberOfStencilPointers = newnum;
348 }
349 
350 template<class TVar>
351 void TPZStencilMatrix<TVar>::SetNodeStencils( int *stencilnumber ) {
352  // Associates the given stencil number with each row
353  if(fStencilNumbers) delete fStencilNumbers;
354  fStencilNumbers = new int[fRows];
355  memcpy(fStencilNumbers,stencilnumber,sizeof(int)*fRows);
356 }
357 
358 
359 // ****************************************************************************
360 //
361 // Subclass MPStencil: Constructors and the destructor
362 //
363 // ****************************************************************************
364 
365 template<class TVar>
366 TPZStencilMatrix<TVar>::MPStencil::MPStencil( int inc, int *IA, TVar *A ) {
367  fInc = inc;
368  fNumberOfItems=0;
369  int *IAPtr = IA;
370  while(*IAPtr) {
371  fNumberOfItems += *IAPtr +1;
372  IAPtr += *IAPtr +1;
373  }
374  fNumberOfItems++;
375  fIA = new int[fNumberOfItems];
376  fA = new TVar[fNumberOfItems];
377  memcpy(fIA,IA,sizeof(int)*fNumberOfItems);
378  memcpy(fA,A,sizeof(TVar)*fNumberOfItems);
379 }
380 template<class TVar>
381 TPZStencilMatrix<TVar>::MPStencil::~MPStencil() {
382  fNumberOfItems = 0;
383  delete fIA; fIA = 0;
384  delete fA; fA = 0;
385 }
TPZStencilMatrix::IncreaseStencilPointers
void IncreaseStencilPointers(int stencilnumber)
Definition: pzstencil.cpp:338
TPZStencilMatrix
Implements a sparse matrix defined by a stencil. Matrix.
Definition: pzstencil.h:19
EInputFormat
Definition: pzmatrix.h:55
TPZStencilMatrix::MPStencil::MPStencil
MPStencil(int inc, int *IA, TVar *A)
Definition: pzstencil.cpp:366
MatrixOutputFormat
MatrixOutputFormat
Defines output format.
Definition: pzmatrix.h:55
pzstencil.h
Contains TPZStencilMatrix class which implements a sparse matrix defined by a stencil. Purpose: Defines operations on sparse matrices stored by stencils. Solvers: SOR and SSOR.
val
REAL val(STATE &number)
Returns value of the variable.
Definition: pzartdiff.h:23
TPZStencilMatrix::SolveSOR
void SolveSOR(int &numiterations, const TPZFMatrix< TVar > &rhs, TPZFMatrix< TVar > &x, TPZFMatrix< TVar > *residual, TPZFMatrix< TVar > &scratch, const TVar overrelax, TVar &tol, const int FromCurrent=0, const int direction=1)
Definition: pzstencil.cpp:269
TPZFMatrix::Zero
int Zero() override
Makes Zero all the elements.
Definition: pzfmatrix.h:651
pzgeom::tol
static const double tol
Definition: pzgeoprism.cpp:23
pzfmatrix.h
Contains TPZMatrixclass which implements full matrix (using column major representation).
TPZStencilMatrix::TPZStencilMatrix
TPZStencilMatrix(int rows, int cols)
sets up the StencilMatrix based on the stencil
Definition: pzstencil.cpp:19
TPZStencilMatrix::GetVal
virtual const TVar & GetVal(const int row, const int col) const
Get the matrix entry at (row,col) without bound checking.
Definition: pzstencil.cpp:59
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
TPZStencilMatrix::ComputeDiagonal
void ComputeDiagonal()
Definition: pzstencil.cpp:238
TPZStencilMatrix::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
computes
Definition: pzstencil.cpp:112
TPZFMatrix
Full matrix class. Matrix.
Definition: pzfmatrix.h:32
TPZStencilMatrix::SetStencil
void SetStencil(int stencilnumber, int inc, int *IA, TVar *A)
initiates Stencil number "stencilnumber" with the data
Definition: pzstencil.cpp:326
TPZStencilMatrix::SetNodeStencils
void SetNodeStencils(int *stencilnumber)
associates the given stencil number with each row
Definition: pzstencil.cpp:351
TPZStencilMatrix::Print
virtual void Print(const char *title, std::ostream &out=std::cout, const MatrixOutputFormat form=EFormatted) const
Print the matrix along with a identification title.
Definition: pzstencil.cpp:188
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
TPZStencilMatrix::~TPZStencilMatrix
virtual ~TPZStencilMatrix()
Definition: pzstencil.cpp:36
Generated on Sun Aug 16 2020 11:17:05 for NeoPZ by   doxygen 1.8.13