Epetra Package Browser (Single Doxygen Collection)  Development
FEVector/ExecuteTestProblems.cpp
Go to the documentation of this file.
1 //@HEADER
2 // ************************************************************************
3 //
4 // Epetra: Linear Algebra Services Package
5 // Copyright 2011 Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ************************************************************************
40 //@HEADER
41 
42 
43 #include "Epetra_BLAS.h"
44 #include "ExecuteTestProblems.h"
45 #include "Epetra_Comm.h"
46 #include "Epetra_Vector.h"
49 #include <vector>
50 
51 int MultiVectorTests(const Epetra_BlockMap & Map, int NumVectors, bool verbose)
52 {
53  (void)NumVectors;
54  const Epetra_Comm & Comm = Map.Comm();
55  int ierr = 0;
56 
57  /* get number of processors and the name of this processor */
58 
59  // int NumProc = Comm.getNumProc();
60  int MyPID = Comm.MyPID();
61 
62  // Construct FEVector
63 
64  if (verbose&&MyPID==0) cout << "constructing Epetra_FEVector" << endl;
65 
66  Epetra_FEVector A(Map, 1);
67 
68  //For an extreme test, we'll have each processor sum-in a 1.0 for All
69  //global ids.
70 
71  int minGID = Map.MinAllGID();
72  int numGlobalIDs = Map.NumGlobalElements();
73 
74  //For now we're going to have just one point associated with
75  //each GID (element).
76 
77  int* ptIndices = new int[numGlobalIDs];
78  double* ptCoefs = new double[numGlobalIDs];
79 
80  Epetra_IntSerialDenseVector epetra_indices(View, ptIndices, numGlobalIDs);
81  Epetra_SerialDenseVector epetra_coefs(View, ptCoefs, numGlobalIDs);
82 
83  {for(int i=0; i<numGlobalIDs; ++i) {
84  ptIndices[i] = minGID+i;
85  ptCoefs[i] = 1.0;
86  }}
87 
88  if (verbose&&MyPID==0) {
89  cout << "calling A.SumIntoGlobalValues with " << numGlobalIDs << " values"<<endl;
90  }
91  EPETRA_TEST_ERR( A.SumIntoGlobalValues(numGlobalIDs, ptIndices, ptCoefs), ierr);
92 
93  if (verbose&&MyPID==0) {
94  cout << "calling A.SumIntoGlobalValues with " << numGlobalIDs << " values"<<endl;
95  }
96  EPETRA_TEST_ERR( A.SumIntoGlobalValues(epetra_indices, epetra_coefs), ierr);
97 
98  if (verbose&&MyPID==0) {
99  cout << "calling A.GlobalAssemble()" << endl;
100  }
101 
102  EPETRA_TEST_ERR( A.GlobalAssemble(), ierr );
103 
104  if (verbose&&MyPID==0) {
105  cout << "after globalAssemble"<<endl;
106  }
107  if (verbose) {
108  A.Print(cout);
109  }
110 
111  //now do a quick test of the copy constructor
113 
114  double nrm2a, nrm2b;
115  A.Norm2(&nrm2a);
116  B.Norm2(&nrm2b);
117 
118  if (nrm2a != nrm2b) {
119  cerr << "copy-constructor test failed, norm of copy doesn't equal"
120  << " norm of original."<<endl;
121  return(-1);
122  }
123 
124  delete [] ptIndices;
125  delete [] ptCoefs;
126 
127  return(ierr);
128 }
129 
130 int fevec0(Epetra_Comm& Comm, bool verbose)
131 {
132  int ierr = 0;
133  int NumGlobalRows = 4;
134  int indexBase = 0;
135  Epetra_Map Map(NumGlobalRows, indexBase, Comm);
136 
137  int Numprocs = Comm.NumProc();
138  int MyPID = Comm.MyPID();
139 
140  if (Numprocs != 2) return(0);
141 
142 
143  int NumCols = 3;
144  int* Indices = new int[NumCols];
145 
146  double* Values = new double[NumCols];
147 
148 // Create vectors
149 
150  Epetra_FEVector b(Map, 1);
151  Epetra_FEVector x0(Map, 1);
152 
153 // source terms
154  NumCols = 2;
155 
156  if(MyPID==0) // indices corresponding to element 0 on processor 0
157  {
158  Indices[0] = 0;
159  Indices[1] = 3;
160 
161  Values[0] = 1./2.;
162  Values[1] = 1./2.;
163 
164  }
165  else
166  {
167  Indices[0] = 1;
168  Indices[1] = 2;
169 
170  Values[0] = 0;
171  Values[1] = 0;
172  }
173 
174  EPETRA_TEST_ERR( b.SumIntoGlobalValues(NumCols, Indices, Values),
175  ierr);
176 
177  EPETRA_TEST_ERR( b.GlobalAssemble(), ierr);
178 
179  if (verbose&&MyPID==0) {
180  cout << "b:"<<endl;
181  }
182 
183  if (verbose) {
184  b.Print(cout);
185  }
186 
187  x0 = b;
188 
189  if (verbose&&MyPID==0) {
190  cout << "x:"<<endl;
191  }
192 
193  if (verbose) {
194  x0.Print(cout);
195  }
196 
197  delete [] Values;
198  delete [] Indices;
199 
200  return(0);
201 }
202 
203 int fevec1(Epetra_Comm& Comm, bool verbose)
204 {
205  int Numprocs = Comm.NumProc();
206 
207  if (Numprocs != 2) return(0);
208  int MyPID = Comm.MyPID();
209 
210  int ierr = 0;
211  int NumGlobalRows = 6;
212  const int NumVectors = 4;
213  int indexBase = 0;
214  Epetra_Map Map(NumGlobalRows, indexBase, Comm);
215 
216  const int Num = 4;
217  int Indices[Num];
218 
219  double Values[Num];
220 
221 // Create vectors
222 
223  Epetra_FEVector b(Map, NumVectors);
224  Epetra_FEVector x0(Map, NumVectors);
225 
226 // source terms
227 
228  if(MyPID==0) // indices corresponding to element 0 on processor 0
229  {
230  Indices[0] = 0;
231  Indices[1] = 1;
232  Indices[2] = 4;
233  Indices[3] = 5;
234 
235  Values[0] = 1./2.;
236  Values[1] = 1./2.;
237  Values[2] = 1./2.;
238  Values[3] = 1./2.;
239 
240  }
241  else
242  {
243  Indices[0] = 1;
244  Indices[1] = 2;
245  Indices[2] = 3;
246  Indices[3] = 4;
247 
248  Values[0] = 0;
249  Values[1] = 0;
250  Values[2] = 0;
251  Values[3] = 0;
252  }
253 
254  for(int i=0; i<NumVectors; ++i) {
255  EPETRA_TEST_ERR( b.SumIntoGlobalValues(Num, Indices, Values, i),
256  ierr);
257  }
258 
259  EPETRA_TEST_ERR( b.GlobalAssemble(), ierr);
260 
261  double nrm2[NumVectors];
262 
263  b.Norm2(nrm2);
264 
265  for(int i=1; i<NumVectors; ++i) {
266  if (fabs(nrm2[i]-nrm2[0]) > 1.e-12) {
267  EPETRA_TEST_ERR(-1, ierr);
268  return(-1);
269  }
270  }
271 
272 
273  //now sum-in again, to make sure the previous call to GlobalAssemble
274  //didn't do something nasty to internal non-local data structures.
275  //(This is a specific case that has bitten me. Hence this test...)
276  for(int i=0; i<NumVectors; ++i) {
277  EPETRA_TEST_ERR( b.SumIntoGlobalValues(Num, Indices, Values, i),
278  ierr);
279  }
280 
281  //and now GlobalAssemble again...
282  EPETRA_TEST_ERR( b.GlobalAssemble(), ierr);
283 
284 
285  if (verbose&&MyPID==0) {
286  cout << "b:"<<endl;
287  }
288 
289  if (verbose) {
290  b.Print(cout);
291  }
292 
293  x0 = b;
294 
295  if (verbose&&MyPID==0) {
296  cout << "x:"<<endl;
297  }
298 
299  if (verbose) {
300  x0.Print(cout);
301  }
302 
303  return(0);
304 }
305 
306 int fevec2(Epetra_Comm& Comm, bool verbose)
307 {
308  int ierr = 0;
309  int NumGlobalElems = 4;
310  int elemSize = 3;
311  int indexBase = 0;
312  Epetra_BlockMap Map(NumGlobalElems, elemSize, indexBase, Comm);
313 
314  int Numprocs = Comm.NumProc();
315  int MyPID = Comm.MyPID();
316 
317  if (Numprocs != 2) return(0);
318 
319  int NumCols = 3;
320  int* Indices = new int[NumCols];
321  int* numValuesPerID = new int[NumCols];
322  for(int i=0; i<NumCols; ++i) {
323  numValuesPerID[i] = elemSize;
324  }
325 
326  double* Values = new double[NumCols*elemSize];
327 
328 // Create vectors
329 
330  Epetra_FEVector b(Map, 1);
331  Epetra_FEVector x0(Map, 1);
332 
333 // source terms
334  NumCols = 2;
335 
336  if(MyPID==0) // indices corresponding to element 0 on processor 0
337  {
338  Indices[0] = 0;
339  Indices[1] = 3;
340 
341  Values[0] = 1./2.;
342  Values[1] = 1./2.;
343  Values[2] = 1./2.;
344  Values[3] = 1./2.;
345  Values[4] = 1./2.;
346  Values[5] = 1./2.;
347 
348  }
349  else
350  {
351  Indices[0] = 1;
352  Indices[1] = 2;
353 
354  Values[0] = 0;
355  Values[1] = 0;
356  Values[2] = 0;
357  Values[3] = 0;
358  Values[4] = 0;
359  Values[5] = 0;
360  }
361 
362  EPETRA_TEST_ERR( b.SumIntoGlobalValues(NumCols, Indices,
363  numValuesPerID, Values),
364  ierr);
365 
366  EPETRA_TEST_ERR( b.GlobalAssemble(), ierr);
367 
368  if (verbose&&MyPID==0) {
369  cout << "b:"<<endl;
370  }
371 
372  if (verbose) {
373  b.Print(cout);
374  }
375 
376  x0 = b;
377 
378  if (verbose&&MyPID==0) {
379  cout << "x:"<<endl;
380  }
381 
382  if (verbose) {
383  x0.Print(cout);
384  }
385 
386  delete [] Values;
387  delete [] Indices;
388  delete [] numValuesPerID;
389 
390  return(0);
391 }
392 
393 int fevec3(Epetra_Comm& Comm, bool verbose)
394 {
395  int ierr = 0;
396  int NumGlobalElems = 4;
397  int elemSize = 40;
398  int indexBase = 0;
399  Epetra_BlockMap Map(NumGlobalElems, elemSize, indexBase, Comm);
400 
401  int Numprocs = Comm.NumProc();
402  int MyPID = Comm.MyPID();
403 
404  if (Numprocs != 2) return(0);
405 
406  int NumCols = 3;
407  int* Indices = new int[NumCols];
408  int* numValuesPerID = new int[NumCols];
409  for(int i=0; i<NumCols; ++i) {
410  numValuesPerID[i] = elemSize;
411  }
412 
413  double* Values = new double[NumCols*elemSize];
414 
415 // Create vectors
416 
417  Epetra_FEVector b(Map, 1);
418  Epetra_FEVector x0(Map, 1);
419 
420 // source terms
421  NumCols = 2;
422 
423  if(MyPID==0) // indices corresponding to element 0 on processor 0
424  {
425  Indices[0] = 0;
426  Indices[1] = 3;
427 
428  for(int ii=0; ii<NumCols*elemSize; ++ii) {
429  Values[ii] = 1./2.;
430  }
431 
432  }
433  else
434  {
435  Indices[0] = 1;
436  Indices[1] = 2;
437 
438  for(int ii=0; ii<NumCols*elemSize; ++ii) {
439  Values[ii] = 0.;
440  }
441 
442  }
443 
444  EPETRA_TEST_ERR( b.SumIntoGlobalValues(NumCols, Indices,
445  numValuesPerID, Values),
446  ierr);
447 
448  EPETRA_TEST_ERR( b.GlobalAssemble(), ierr);
449 
450  if (verbose&&MyPID==0) {
451  cout << "b:"<<endl;
452  }
453 
454  if (verbose) {
455  b.Print(cout);
456  }
457 
458  x0 = b;
459 
460  if (verbose&&MyPID==0) {
461  cout << "x:"<<endl;
462  }
463 
464  if (verbose) {
465  x0.Print(cout);
466  }
467 
468  delete [] Values;
469  delete [] Indices;
470  delete [] numValuesPerID;
471 
472  return(0);
473 }
474 
475 int fevec4(Epetra_Comm& Comm, bool verbose)
476 {
477  int NumElements = 4;
478  Epetra_Map Map(NumElements, 0, Comm);
479  Epetra_FEVector x1(Map);
480  const double value = 1.;
481  x1.PutScalar (value);
482  // replace one element by itself. processor 0
483  // does not own this element
484  const int GID = 3;
485  x1.ReplaceGlobalValues(1, &GID, &value);
486  x1.GlobalAssemble (Insert);
487 
488  if (Map.MyGID(3)) {
489  //insist that the value for GID==3 is 1:
490  if (std::abs(x1.Values()[Map.LID(3)] - 1) > 1.e-9) return -1;
491  }
492 
493  std::cout << x1;
494 
495  Comm.Barrier();
496 
497  // re-apply GlobalAssemble. Nothing should
498  // happen
499  x1.GlobalAssemble (Insert);
500  std::cout << x1;
501  if (Map.MyGID(3)) {
502  //insist that the value for GID==3 is 1:
503  if (std::abs(x1.Values()[Map.LID(3)] - 1) > 1.e-9) return -1;
504  }
505 
506  return 0;
507 }
508 
509 int fevec5(Epetra_Comm& Comm, bool verbose)
510 {
511  int NumElements = 4;
512  Epetra_Map Map(NumElements, 0, Comm);
513  Epetra_FEVector x1(Map);
514  x1.PutScalar (0);
515 
516  // let all processors set global entry 0 to 1
517  const int GID = 0;
518  const double value = 1;
519  x1.ReplaceGlobalValues(1, &GID, &value);
520  x1.GlobalAssemble (Insert);
521  if (Comm.MyPID()==0)
522  std::cout << "Entry " << GID << " after construct & set: "
523  << x1[0][0] << std::endl;
524 
525  // copy vector
526  Epetra_FEVector x2 (x1);
527 
528  x2.PutScalar(0);
529 
530  // re-apply 1 to the vector, but only on the
531  // owning processor. should be enough to set
532  // the value (as non-local data in x1 should
533  // have been eliminated after calling
534  // GlobalAssemble).
535  if (Comm.MyPID()==0)
536  x2.ReplaceGlobalValues(1, &GID, &value);
537  x2.GlobalAssemble (Insert);
538 
539  if (Comm.MyPID()==0)
540  std::cout << "Entry " << GID << " after copy & set: "
541  << x2[0][0] << std::endl;
542 
543  return 0;
544 }
545 
546 int fevec6(Epetra_Comm& Comm, bool verbose)
547 {
548  int NumElements = 4;
549  Epetra_Map Map(NumElements, 0, Comm);
550  Epetra_FEVector x1(Map);
551  x1.PutScalar (0);
552 
553  // let all processors set global entry 0 to 1
554  const int GID = 0;
555  const double value = 1;
556  x1.ReplaceGlobalValues(1, &GID, &value);
557  x1.GlobalAssemble (Insert);
558  if (Comm.MyPID()==0)
559  std::cout << "Entry " << GID << " after construct & set: "
560  << x1[0][0] << std::endl;
561 
562  x1.PutScalar(0);
563 
564  // re-apply 1 to the vector, but only on the
565  // owning processor. should be enough to set
566  // the value (as non-local data in x1 should
567  // have been eliminated after calling
568  // GlobalAssemble).
569  if (Comm.MyPID()==0)
570  x1.ReplaceGlobalValues(1, &GID, &value);
571  x1.GlobalAssemble (Insert);
572 
573  if (Comm.MyPID()==0) {
574  std::cout << "Entry " << GID << " after PutScalar & set: "
575  << x1[0][0] << std::endl;
576  if (x1[0][0] != value) return -1;
577  }
578 
579  return 0;
580 }
581 
582 int fevec7(Epetra_Comm& Comm, bool verbose)
583 {
584  const int NumVectors = 4;
585  const int NumElements = 4;
586  Epetra_Map Map(NumElements, 0, Comm);
587  std::vector<double> mydata(NumElements*NumVectors, 1.0);
588  Epetra_FEVector x1(View, Map, &mydata[0], NumElements, NumVectors);
589 
590  x1.PutScalar (0);
591 
592  // let all processors set global entry 0 to 1
593  const int GID = 0;
594  const double value = 1;
595  x1.ReplaceGlobalValues(1, &GID, &value);
596  x1.GlobalAssemble (Insert);
597 
598  if (Comm.MyPID()==0 && x1[0][0] != value) return -1;
599  return 0;
600 }
601 
int Norm2(double *Result) const
Compute 2-norm of each vector in multi-vector.
double * Values() const
Get pointer to MultiVector values.
Epetra_Map: A class for partitioning vectors and matrices.
Definition: Epetra_Map.h:119
Epetra_IntSerialDenseVector: A class for constructing and using dense vectors.
int fevec1(Epetra_Comm &Comm, bool verbose)
int MinAllGID() const
Returns the minimum global ID across the entire map.
virtual void Print(std::ostream &os) const
Print method.
#define EPETRA_TEST_ERR(a, b)
int fevec7(Epetra_Comm &Comm, bool verbose)
int PutScalar(double ScalarConstant)
Initialize all values in a multi-vector with constant value.
int fevec6(Epetra_Comm &Comm, bool verbose)
virtual void Barrier() const =0
Epetra_Comm Barrier function.
virtual int MyPID() const =0
Return my process ID.
Epetra_SerialDenseVector: A class for constructing and using dense vectors.
bool MyGID(int GID_in) const
Returns true if the GID passed in belongs to the calling processor in this map, otherwise returns fal...
Epetra_Comm: The Epetra Communication Abstract Base Class.
Definition: Epetra_Comm.h:73
int fevec5(Epetra_Comm &Comm, bool verbose)
Epetra_BlockMap: A class for partitioning block element vectors and matrices.
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
int fevec2(Epetra_Comm &Comm, bool verbose)
int LID(int GID) const
Returns local ID of global ID, return -1 if not found on this processor.
int MultiVectorTests(const Epetra_BlockMap &Map, int NumVectors, bool verbose)
int ReplaceGlobalValues(int numIDs, const int *GIDs, const double *values, int vectorIndex=0)
Copy values into the vector overwriting any values that already exist for the specified indices...
virtual int NumProc() const =0
Returns total number of processes.
int GlobalAssemble(Epetra_CombineMode mode=Add, bool reuse_map_and_exporter=false)
Gather any overlapping/shared data into the non-overlapping partitioning defined by the Map that was ...
Epetra Finite-Element Vector.
int NumGlobalElements() const
Number of elements across all processors.
int fevec0(Epetra_Comm &Comm, bool verbose)
int fevec3(Epetra_Comm &Comm, bool verbose)
int fevec4(Epetra_Comm &Comm, bool verbose)
int SumIntoGlobalValues(int numIDs, const int *GIDs, const double *values, int vectorIndex=0)
Accumulate values into the vector, adding them to any values that already exist for the specified ind...