Epetra Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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
112  Epetra_FEVector B(A);
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 NumGlobalElements() const
Number of elements across all processors.
Epetra_Map: A class for partitioning vectors and matrices.
Definition: Epetra_Map.h:127
Epetra_IntSerialDenseVector: A class for constructing and using dense vectors.
int fevec1(Epetra_Comm &Comm, bool verbose)
#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.
virtual void Print(std::ostream &os) const
Print method.
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.
int MultiVectorTests(const Epetra_Map &Map, int NumVectors, bool verbose)
Epetra_Comm: The Epetra Communication Abstract Base Class.
Definition: Epetra_Comm.h:81
int fevec5(Epetra_Comm &Comm, bool verbose)
Epetra_BlockMap: A class for partitioning block element vectors and matrices.
int MinAllGID() const
Returns the minimum global ID across the entire map.
int LID(int GID) const
Returns local ID of global ID, return -1 if not found on this processor.
bool MyGID(int GID_in) const
Returns true if the GID passed in belongs to the calling processor in this map, otherwise returns fal...
int fevec2(Epetra_Comm &Comm, bool verbose)
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
double * Values() const
Get pointer to MultiVector values.
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 Norm2(double *Result) const
Compute 2-norm of each vector in multi-vector.
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...