Epetra Package Browser (Single Doxygen Collection)  Development
FECrsMatrix/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"
47 #include "Epetra_FEVector.h"
50 #include <../src/Epetra_matrix_data.h>
51 #include <../src/Epetra_test_functions.h>
52 
53 
54 int Drumm1(const Epetra_Map& map, bool verbose)
55 {
56  (void)verbose;
57  //Simple 2-element problem (element as in "finite-element") from
58  //Clif Drumm. Two triangular elements, one per processor, as shown
59  //here:
60  //
61  // *----*
62  // 3|\ 2|
63  // | \ |
64  // | 0\1|
65  // | \|
66  // *----*
67  // 0 1
68  //
69  //Element 0 on processor 0, element 1 on processor 1.
70  //Processor 0 will own nodes 0,1 and processor 1 will own nodes 2,3.
71  //Each processor will pass a 3x3 element-matrix to Epetra_FECrsMatrix.
72  //After GlobalAssemble(), the matrix should be as follows:
73  //
74  // row 0: 2 1 0 1
75  //proc 0 row 1: 1 4 1 2
76  //----------------------------------
77  // row 2: 0 1 2 1
78  //proc 1 row 3: 1 2 1 4
79  //
80 
81  int numProcs = map.Comm().NumProc();
82  int localProc = map.Comm().MyPID();
83 
84  if (numProcs != 2) return(0);
85 
86  //so first we'll set up a epetra_test::matrix_data object with
87  //contents that match the above-described matrix. (but the
88  //matrix_data object will have all 4 rows on each processor)
89 
90  int i;
91  int rowlengths[4];
92  rowlengths[0] = 3;
93  rowlengths[1] = 4;
94  rowlengths[2] = 3;
95  rowlengths[3] = 4;
96 
97  epetra_test::matrix_data matdata(4, rowlengths);
98  for(i=0; i<4; ++i) {
99  for(int j=0; j<matdata.rowlengths()[i]; ++j) {
100  matdata.colindices()[i][j] = j;
101  }
102  }
103 
104  matdata.colindices()[0][2] = 3;
105 
106  matdata.colindices()[2][0] = 1;
107  matdata.colindices()[2][1] = 2;
108  matdata.colindices()[2][2] = 3;
109 
110  double** coefs = matdata.coefs();
111  coefs[0][0] = 2.0; coefs[0][1] = 1.0; coefs[0][2] = 1.0;
112  coefs[1][0] = 1.0; coefs[1][1] = 4.0; coefs[1][2] = 1.0; coefs[1][3] = 2.0;
113  coefs[2][0] = 1.0; coefs[2][1] = 2.0; coefs[2][2] = 1.0;
114  coefs[3][0] = 1.0; coefs[3][1] = 2.0; coefs[3][2] = 1.0; coefs[3][3] = 4.0;
115 
116  //now we'll load a Epetra_FECrsMatrix with data that matches the
117  //above-described finite-element problem.
118 
119  int indexBase = 0, ierr = 0;
120  int myNodes[4];
121  double values[9];
122  values[0] = 2.0;
123  values[1] = 1.0;
124  values[2] = 1.0;
125  values[3] = 1.0;
126  values[4] = 2.0;
127  values[5] = 1.0;
128  values[6] = 1.0;
129  values[7] = 1.0;
130  values[8] = 2.0;
131 
132  int numMyNodes = 2;
133 
134  if (localProc == 0) {
135  myNodes[0] = 0;
136  myNodes[1] = 1;
137  }
138  else {
139  myNodes[0] = 2;
140  myNodes[1] = 3;
141  }
142 
143  Epetra_Map Map(-1, numMyNodes, myNodes, indexBase, map.Comm());
144 
145  numMyNodes = 3;
146 
147  if (localProc == 0) {
148  myNodes[0] = 0;
149  myNodes[1] = 1;
150  myNodes[2] = 3;
151  }
152  else {
153  myNodes[0] = 1;
154  myNodes[1] = 2;
155  myNodes[2] = 3;
156  }
157 
158  int rowLengths = 3;
159  Epetra_FECrsMatrix A(Copy, Map, rowLengths);
160 
161  EPETRA_TEST_ERR( A.InsertGlobalValues(numMyNodes, myNodes,
162  numMyNodes, myNodes, values,
164 
165  EPETRA_TEST_ERR( A.GlobalAssemble(), ierr );
166  EPETRA_TEST_ERR( A.GlobalAssemble(), ierr );
167 
168  //now the test is to check whether the FECrsMatrix data matches the
169  //epetra_test::matrix_data object...
170 
171  bool the_same = matdata.compare_local_data(A);
172 
173  if (!the_same) {
174  return(-1);
175  }
176 
177  return(0);
178 }
179 
180 int Drumm2(const Epetra_Map& map, bool verbose)
181 {
182  //Simple 2-element problem (element as in "finite-element") from
183  //Clif Drumm. Two triangular elements, one per processor, as shown
184  //here:
185  //
186  // *----*
187  // 3|\ 2|
188  // | \ |
189  // | 0\1|
190  // | \|
191  // *----*
192  // 0 1
193  //
194  //Element 0 on processor 0, element 1 on processor 1.
195  //Processor 0 will own nodes 0,1,3 and processor 1 will own node 2.
196  //Each processor will pass a 3x3 element-matrix to Epetra_FECrsMatrix.
197  //After GlobalAssemble(), the matrix should be as follows:
198  //
199  // row 0: 2 1 0 1
200  //proc 0 row 1: 1 4 1 2
201  // row 2: 0 1 2 1
202  //----------------------------------
203  //proc 1 row 3: 1 2 1 4
204  //
205 
206  int numProcs = map.Comm().NumProc();
207  int localProc = map.Comm().MyPID();
208 
209  if (numProcs != 2) return(0);
210 
211  int indexBase = 0, ierr = 0;
212 
213  double* values = new double[9];
214  values[0] = 2.0;
215  values[1] = 1.0;
216  values[2] = 1.0;
217  values[3] = 1.0;
218  values[4] = 2.0;
219  values[5] = 1.0;
220  values[6] = 1.0;
221  values[7] = 1.0;
222  values[8] = 2.0;
223 
224  if (localProc == 0) {
225  int numMyNodes = 3;
226  int* myNodes = new int[numMyNodes];
227  myNodes[0] = 0;
228  myNodes[1] = 1;
229  myNodes[2] = 3;
230 
231  Epetra_Map Map(-1, numMyNodes, myNodes, indexBase, map.Comm());
232 
233  int rowLengths = 3;
234  Epetra_FECrsMatrix A(Copy, Map, rowLengths);
235 
236  EPETRA_TEST_ERR( A.InsertGlobalValues(numMyNodes, myNodes,
237  numMyNodes, myNodes,
238  values, Epetra_FECrsMatrix::ROW_MAJOR),ierr);
239 
240  EPETRA_TEST_ERR( A.GlobalAssemble(), ierr );
241  EPETRA_TEST_ERR( A.GlobalAssemble(), ierr );
242 
243  if (verbose) {
244  A.Print(cout);
245  }
246 
247  //now let's make sure we can do a matvec with this matrix.
248  Epetra_Vector x(Map), y(Map);
249  x.PutScalar(1.0);
250  EPETRA_TEST_ERR( A.Multiply(false, x, y), ierr);
251 
252  if (verbose&&localProc==0) {
253  cout << "y = A*x, x=1.0's"<<endl;
254  }
255 
256  if (verbose) {
257  y.Print(cout);
258  }
259 
260  delete [] myNodes;
261  delete [] values;
262  }
263  else {
264  int numMyNodes = 1;
265  int* myNodes = new int[numMyNodes];
266  myNodes[0] = 2;
267 
268  Epetra_Map Map(-1, numMyNodes, myNodes, indexBase, map.Comm());
269 
270  int rowLengths = 3;
271  Epetra_FECrsMatrix A(Copy, Map, rowLengths);
272 
273  delete [] myNodes;
274  numMyNodes = 3;
275  myNodes = new int[numMyNodes];
276  myNodes[0] = 1;
277  myNodes[1] = 2;
278  myNodes[2] = 3;
279 
280  EPETRA_TEST_ERR( A.InsertGlobalValues(numMyNodes, myNodes,
281  numMyNodes, myNodes,
282  values, Epetra_FECrsMatrix::ROW_MAJOR),ierr);
283 
284  EPETRA_TEST_ERR( A.GlobalAssemble(), ierr );
285  EPETRA_TEST_ERR( A.GlobalAssemble(), ierr );
286 
287  if (verbose) {
288  A.Print(cout);
289  }
290 
291  //now let's make sure we can do a matvec with this matrix.
292  Epetra_Vector x(Map), y(Map);
293  x.PutScalar(1.0);
294  EPETRA_TEST_ERR( A.Multiply(false, x, y), ierr);
295 
296  if (verbose) {
297  y.Print(cout);
298  }
299 
300  delete [] myNodes;
301  delete [] values;
302  }
303 
304  return(0);
305 }
306 
307 int Drumm3(const Epetra_Map& map, bool verbose)
308 {
309  const Epetra_Comm & Comm = map.Comm();
310 
311  /* get number of processors and the name of this processor */
312 
313  int Numprocs = Comm.NumProc();
314  int MyPID = Comm.MyPID();
315 
316  if (Numprocs != 2) return(0);
317 
318  int NumGlobalRows = 4;
319  int IndexBase = 0;
320  Epetra_Map Map(NumGlobalRows, IndexBase, Comm);
321 
322  // Construct FECrsMatrix
323 
324  int NumEntriesPerRow = 3;
325 
326  Epetra_FECrsMatrix A(Copy, Map, NumEntriesPerRow);
327 
328  double ElementArea = 0.5;
329 
330  int NumCols = 3;
331  int* Indices = new int[NumCols];
332 
333  if(MyPID==0) // indices corresponding to element 0 on processor 0
334  {
335  Indices[0] = 0;
336  Indices[1] = 1;
337  Indices[2] = 3;
338  }
339  else if(MyPID==1) // indices corresponding to element 1 on processor 1
340  {
341  Indices[0] = 1;
342  Indices[1] = 2;
343  Indices[2] = 3;
344  }
345 
346  double* Values = new double[NumCols*NumCols];
347 
348 // removal term
349  Values[0] = 2*ElementArea/12.;
350  Values[1] = 1*ElementArea/12.;
351  Values[2] = 1*ElementArea/12.;
352  Values[3] = 1*ElementArea/12.;
353  Values[4] = 2*ElementArea/12.;
354  Values[5] = 1*ElementArea/12.;
355  Values[6] = 1*ElementArea/12.;
356  Values[7] = 1*ElementArea/12.;
357  Values[8] = 2*ElementArea/12.;
358 
359  A.InsertGlobalValues(NumCols, Indices,
360  Values,
362 
363  A.GlobalAssemble();
364  A.GlobalAssemble();
365 
366 // A.Print(cout);
367 
368 // Create vectors for CG algorithm
369 
370  Epetra_FEVector* bptr = new Epetra_FEVector(A.RowMap(), 1);
371  Epetra_FEVector* x0ptr = new Epetra_FEVector(A.RowMap(), 1);
372 
373  Epetra_FEVector& b = *bptr;
374  Epetra_FEVector& x0 = *x0ptr;
375 
376  // source terms
377  NumCols = 2;
378 
379  if(MyPID==0) // indices corresponding to element 0 on processor 0
380  {
381  Indices[0] = 0;
382  Indices[1] = 3;
383 
384  Values[0] = 1./2.;
385  Values[1] = 1./2.;
386 
387  }
388  else
389  {
390  Indices[0] = 1;
391  Indices[1] = 2;
392 
393  Values[0] = 0;
394  Values[1] = 0;
395  }
396 
397  b.SumIntoGlobalValues(NumCols, Indices, Values);
398 
399  b.GlobalAssemble();
400 
401  if (verbose&&MyPID==0) cout << "b:" << endl;
402  if (verbose) {
403  b.Print(cout);
404  }
405 
406  x0 = b;
407 
408  if (verbose&&MyPID==0) {
409  cout << "x:"<<endl;
410  }
411 
412  if (verbose) {
413  x0.Print(cout);
414  }
415 
416  delete [] Values;
417  delete [] Indices;
418 
419  delete bptr;
420  delete x0ptr;
421 
422  return(0);
423 }
424 
425 int four_quads(const Epetra_Comm& Comm, bool preconstruct_graph, bool verbose)
426 {
427  if (verbose) {
428  cout << "******************* four_quads ***********************"<<endl;
429  }
430 
431  //This function assembles a matrix representing a finite-element mesh
432  //of four 2-D quad elements. There are 9 nodes in the problem. The
433  //same problem is assembled no matter how many processors are being used
434  //(within reason). It may not work if more than 9 processors are used.
435  //
436  // *------*------*
437  // 6| 7| 8|
438  // | E2 | E3 |
439  // *------*------*
440  // 3| 4| 5|
441  // | E0 | E1 |
442  // *------*------*
443  // 0 1 2
444  //
445  //Nodes are denoted by * with node-numbers below and left of each node.
446  //E0, E1 and so on are element-numbers.
447  //
448  //Each processor will contribute a sub-matrix of size 4x4, filled with 1's,
449  //for each element. Thus, the coefficient value at position 0,0 should end up
450  //being 1.0*numProcs, the value at position 4,4 should be 1.0*4*numProcs, etc.
451  //
452  //Depending on the number of processors being used, the locations of the
453  //specific matrix positions (in terms of which processor owns them) will vary.
454  //
455 
456  int numProcs = Comm.NumProc();
457 
458  int numNodes = 9;
459  int numElems = 4;
460  int numNodesPerElem = 4;
461 
462  int indexBase = 0;
463 
464  //Create a map using epetra-defined linear distribution.
465  Epetra_Map map(numNodes, indexBase, Comm);
466 
467  Epetra_CrsGraph* graph = NULL;
468 
469  int* nodes = new int[numNodesPerElem];
470  int i, j, err = 0;
471 
472  if (preconstruct_graph) {
473  graph = new Epetra_CrsGraph(Copy, map, 1);
474 
475  //we're going to fill the graph with indices, but remember it will only
476  //accept indices in rows for which map.MyGID(row) is true.
477 
478  for(i=0; i<numElems; ++i) {
479  switch(i) {
480  case 0:
481  nodes[0] = 0; nodes[1] = 1; nodes[2] = 4; nodes[3] = 3;
482  break;
483  case 1:
484  nodes[0] = 1; nodes[1] = 2; nodes[2] = 5; nodes[3] = 4;
485  break;
486  case 2:
487  nodes[0] = 3; nodes[1] = 4; nodes[2] = 7; nodes[3] = 6;
488  break;
489  case 3:
490  nodes[0] = 4; nodes[1] = 5; nodes[2] = 8; nodes[3] = 7;
491  break;
492  }
493 
494  for(j=0; j<numNodesPerElem; ++j) {
495  if (map.MyGID(nodes[j])) {
496  err = graph->InsertGlobalIndices(nodes[j], numNodesPerElem,
497  nodes);
498  if (err<0) return(err);
499  }
500  }
501  }
502 
503  EPETRA_CHK_ERR( graph->FillComplete() );
504  }
505 
506  Epetra_FECrsMatrix* A = NULL;
507 
508  if (preconstruct_graph) {
509  A = new Epetra_FECrsMatrix(Copy, *graph);
510  }
511  else {
512  A = new Epetra_FECrsMatrix(Copy, map, 1);
513  }
514 
515  EPETRA_CHK_ERR( A->PutScalar(0.0) );
516 
517  double* values_1d = new double[numNodesPerElem*numNodesPerElem];
518  double** values_2d = new double*[numNodesPerElem];
519 
520  for(i=0; i<numNodesPerElem*numNodesPerElem; ++i) values_1d[i] = 1.0;
521 
522  int offset = 0;
523  for(i=0; i<numNodesPerElem; ++i) {
524  values_2d[i] = &(values_1d[offset]);
525  offset += numNodesPerElem;
526  }
527 
528  int format = Epetra_FECrsMatrix::ROW_MAJOR;
529  Epetra_IntSerialDenseVector epetra_nodes(View, nodes, numNodesPerElem);
530  Epetra_SerialDenseMatrix epetra_values(View, values_1d, numNodesPerElem,
531  numNodesPerElem, numNodesPerElem);
532 
533  for(i=0; i<numElems; ++i) {
534  switch(i) {
535  case 0:
536  nodes[0] = 0; nodes[1] = 1; nodes[2] = 4; nodes[3] = 3;
537  if (preconstruct_graph) {
538  err = A->SumIntoGlobalValues(epetra_nodes,
539  epetra_values, format);
540  if (err<0) return(err);
541  }
542  else {
543  err = A->InsertGlobalValues(epetra_nodes,
544  epetra_values, format);
545  if (err<0) return(err);
546  }
547  break;
548 
549  case 1:
550  nodes[0] = 1; nodes[1] = 2; nodes[2] = 5; nodes[3] = 4;
551  if (preconstruct_graph) {
552  err = A->SumIntoGlobalValues(nodes[0], numNodesPerElem, values_2d[0],
553  nodes);
554  err += A->SumIntoGlobalValues(nodes[1], numNodesPerElem, values_2d[1],
555  nodes);
556  err += A->SumIntoGlobalValues(nodes[2], numNodesPerElem, values_2d[2],
557  nodes);
558  err += A->SumIntoGlobalValues(nodes[3], numNodesPerElem, values_2d[3],
559  nodes);
560  if (err<0) return(err);
561  }
562  else {
563  err = A->InsertGlobalValues(numNodesPerElem, nodes,
564  values_2d, format);
565  if (err<0) return(err);
566  }
567  break;
568 
569  case 2:
570  nodes[0] = 3; nodes[1] = 4; nodes[2] = 7; nodes[3] = 6;
571  if (preconstruct_graph) {
572  err = A->SumIntoGlobalValues(numNodesPerElem, nodes,
573  numNodesPerElem, nodes,
574  values_1d, format);
575  if (err<0) return(err);
576  }
577  else {
578  err = A->InsertGlobalValues(numNodesPerElem, nodes,
579  numNodesPerElem, nodes,
580  values_1d, format);
581  if (err<0) return(err);
582  }
583  break;
584 
585  case 3:
586  nodes[0] = 4; nodes[1] = 5; nodes[2] = 8; nodes[3] = 7;
587  if (preconstruct_graph) {
588  err = A->SumIntoGlobalValues(numNodesPerElem, nodes,
589  numNodesPerElem, nodes,
590  values_2d, format);
591  if (err<0) return(err);
592  }
593  else {
594  err = A->InsertGlobalValues(numNodesPerElem, nodes,
595  numNodesPerElem, nodes,
596  values_2d, format);
597  if (err<0) return(err);
598  }
599  break;
600  }
601  }
602 
603  err = A->GlobalAssemble();
604  if (err < 0) {
605  return(err);
606  }
607 
608  Epetra_Vector x(A->RowMap()), y(A->RowMap());
609 
610  x.PutScalar(1.0); y.PutScalar(0.0);
611 
612  Epetra_FECrsMatrix Acopy(*A);
613 
614  err = Acopy.GlobalAssemble();
615  if (err < 0) {
616  return(err);
617  }
618 
619  bool the_same = epetra_test::compare_matrices(*A, Acopy);
620  if (!the_same) {
621  return(-1);
622  }
623 
624  Epetra_FECrsMatrix Acopy2(Copy, A->RowMap(), A->ColMap(), 1);
625 
626  Acopy2 = Acopy;
627 
628  the_same = epetra_test::compare_matrices(*A, Acopy);
629  if (!the_same) {
630  return(-1);
631  }
632 
633  int len = 20;
634  int* indices = new int[len];
635  double* values = new double[len];
636  int numIndices;
637 
638  if (map.MyGID(0)) {
639  EPETRA_CHK_ERR( A->ExtractGlobalRowCopy(0, len, numIndices,
640  values, indices) );
641  if (numIndices != 4) {
642  return(-1);
643  }
644  if (indices[0] != 0) {
645  return(-2);
646  }
647 
648  if (values[0] != 1.0*numProcs) {
649  cout << "ERROR: values[0] ("<<values[0]<<") should be "<<numProcs<<endl;
650  return(-3);
651  }
652  }
653 
654  if (map.MyGID(4)) {
655  EPETRA_CHK_ERR( A->ExtractGlobalRowCopy(4, len, numIndices,
656  values, indices) );
657 
658  if (numIndices != 9) {
659  return(-4);
660  }
661  int lcid = A->LCID(4);
662  if (lcid<0) {
663  return(-5);
664  }
665  if (values[lcid] != 4.0*numProcs) {
666  cout << "ERROR: values["<<lcid<<"] ("<<values[lcid]<<") should be "
667  <<4*numProcs<<endl;
668  return(-6);
669  }
670  }
671 
672  delete [] values_2d;
673  delete [] values_1d;
674  delete [] nodes;
675  delete [] indices;
676  delete [] values;
677 
678  delete A;
679  delete graph;
680 
681  return(0);
682 }
683 
684 int submatrix_formats(const Epetra_Comm& Comm, bool verbose)
685 {
686  (void)verbose;
687  //
688  //This function simply verifies that the ROW_MAJOR/COLUMN_MAJOR switch works.
689  //
690  int numProcs = Comm.NumProc();
691  int myPID = Comm.MyPID();
692 
693  int numLocalElements = 3;
694  int numGlobalElements = numLocalElements*numProcs;
695  int indexBase = 0;
696 
697  Epetra_Map map(numGlobalElements, numLocalElements, indexBase, Comm);
698 
699  Epetra_FECrsMatrix A(Copy, map, numLocalElements);
700 
701  Epetra_IntSerialDenseVector epetra_indices(numLocalElements);
702 
703  int firstGlobalElement = numLocalElements*myPID;
704 
705  int i, j;
706  for(i=0; i<numLocalElements; ++i) {
707  epetra_indices[i] = firstGlobalElement+i;
708  }
709 
710  Epetra_SerialDenseMatrix submatrix(numLocalElements, numLocalElements);
711 
712  for(i=0; i<numLocalElements; ++i) {
713  for(j=0; j<numLocalElements; ++j) {
714  submatrix(i,j) = 1.0*(firstGlobalElement+i);
715  }
716  }
717 
718  EPETRA_CHK_ERR( A.InsertGlobalValues(epetra_indices, submatrix,
720 
721  EPETRA_CHK_ERR( A.GlobalAssemble() );
722 
723  int len = 20;
724  int numIndices;
725  int* indices = new int[len];
726  double* coefs = new double[len];
727 
728  for(i=0; i<numLocalElements; ++i) {
729  int row = firstGlobalElement+i;
730 
731  EPETRA_CHK_ERR( A.ExtractGlobalRowCopy(row, len, numIndices,
732  coefs, indices) );
733 
734  for(j=0; j<numIndices; ++j) {
735  if (coefs[j] != 1.0*row) {
736  return(-2);
737  }
738  }
739  }
740 
741  //now reset submatrix (transposing the i,j indices)
742 
743  for(i=0; i<numLocalElements; ++i) {
744  for(j=0; j<numLocalElements; ++j) {
745  submatrix(j,i) = 1.0*(firstGlobalElement+i);
746  }
747  }
748 
749  //sum these values into the matrix using the ROW_MAJOR switch, which should
750  //result in doubling what's already there from the above Insert operation.
751 
752  EPETRA_CHK_ERR( A.SumIntoGlobalValues(epetra_indices, submatrix,
754 
755  EPETRA_CHK_ERR( A.GlobalAssemble() );
756 
757  for(i=0; i<numLocalElements; ++i) {
758  int row = firstGlobalElement+i;
759 
760  EPETRA_CHK_ERR( A.ExtractGlobalRowCopy(row, len, numIndices,
761  coefs, indices) );
762 
763  for(j=0; j<numIndices; ++j) {
764  if (coefs[j] != 2.0*row) {
765  return(-3);
766  }
767  }
768  }
769 
770  delete [] indices;
771  delete [] coefs;
772 
773  return(0);
774 }
775 
776 int rectangular(const Epetra_Comm& Comm, bool verbose)
777 {
778  (void)verbose;
779  int numprocs = Comm.NumProc();
780  int localproc = Comm.MyPID();
781  int numMyRows = 2;
782  int numGlobalRows = numprocs*numMyRows;
783  int* myrows = new int[numMyRows];
784 
785  int myFirstRow = localproc*numMyRows;
786  int i;
787  for(i=0; i<numMyRows; ++i) {
788  myrows[i] = myFirstRow+i;
789  }
790 
791  Epetra_Map map(numGlobalRows, numMyRows, myrows, 0, Comm);
792 
793  Epetra_FECrsMatrix A(Copy, map, 30);
794 
795  int numcols = 20;
796  int* cols = new int[numcols];
797  for(i=0; i<numcols; ++i) {
798  cols[i] = i;
799  }
800 
801  double* coefs = new double[numGlobalRows*numcols];
802  int offset = 0;
803  for(int j=0; j<numcols; ++j) {
804  for(i=0; i<numGlobalRows; ++i) {
805  coefs[offset++] = 1.0*i;
806  }
807  }
808 
809  int* globalRows = new int[numGlobalRows];
810  for(i=0; i<numGlobalRows; ++i) globalRows[i] = i;
811 
812  EPETRA_CHK_ERR( A.InsertGlobalValues(numGlobalRows, globalRows,
813  numcols, cols, coefs,
815  delete [] coefs;
816  delete [] globalRows;
817 
818  //Since the matrix is rectangular, we need to call GlobalAssemble with
819  //a domain-map and a range-map. Otherwise, GlobalAssemble has no way of
820  //knowing what the domain-map and range-map should be.
821  //We'll use a linear distribution of the columns for a domain-map, and
822  //our original row-map for the range-map.
823  int numMyCols = numcols/numprocs;
824  int rem = numcols%numprocs;
825  if (localproc<rem) ++numMyCols;
826  Epetra_Map domainmap(numcols, numMyCols, 0, Comm);
827 
828  EPETRA_CHK_ERR( A.GlobalAssemble(domainmap, map) );
829 
830  int numGlobalCols = A.NumGlobalCols();
831  int numGlobalNNZ = A.NumGlobalNonzeros();
832 
833  if (numGlobalCols != numcols ||
834  numGlobalNNZ != numGlobalRows*numcols) {
835  return(-1);
836  }
837 
838  delete [] cols;
839  delete [] myrows;
840  return(0);
841 }
842 
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 FillComplete()
Tranform to local index space. Perform other operations to allow optimal matrix operations.
int rectangular(const Epetra_Comm &Comm, bool verbose)
bool compare_local_data(const Epetra_CrsMatrix &A)
The portion of this matrix_data object&#39;s data that corresponds to the locally-owned rows of A...
int Drumm1(const Epetra_Map &map, bool verbose)
virtual void Print(std::ostream &os) const
Print method.
#define EPETRA_TEST_ERR(a, b)
int Drumm2(const Epetra_Map &map, bool verbose)
int PutScalar(double ScalarConstant)
Initialize all values in a multi-vector with constant value.
int Drumm3(const Epetra_Map &map, bool verbose)
#define EPETRA_CHK_ERR(a)
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.
Epetra_SerialDenseMatrix: A class for constructing and using real double precision general dense matr...
virtual int MyPID() const =0
Return my process ID.
Epetra Finite-Element CrsMatrix.
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 GlobalAssemble(bool callFillComplete=true, Epetra_CombineMode combineMode=Add, bool save_off_and_reuse_map_exporter=false)
Gather any overlapping/shared data into the non-overlapping partitioning defined by the Map that was ...
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
bool compare_matrices(const Epetra_CrsMatrix &A, const Epetra_CrsMatrix &B)
Check whether the two CrsMatrix arguments have the same size, structure and coefs.
int InsertGlobalIndices(int GlobalRow, int NumIndices, int *Indices)
Enter a list of elements in a specified global row of the graph.
matrix_data is a very simple data source to be used for filling test matrices.
virtual int NumProc() const =0
Returns total number of processes.
Epetra Finite-Element Vector.
int four_quads(const Epetra_Comm &Comm, bool preconstruct_graph, bool verbose)
int submatrix_formats(const Epetra_Comm &Comm, bool verbose)
Epetra_CrsGraph: A class for constructing and using sparse compressed row graphs. ...
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...