Epetra Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
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 
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.
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 Multiply(bool TransA, const Epetra_Vector &x, Epetra_Vector &y) const
Returns the result of a Epetra_CrsMatrix multiplied by a Epetra_Vector x in y.
#define EPETRA_TEST_ERR(a, b)
int PutScalar(double ScalarConstant)
Initialize all values in a multi-vector with constant value.
virtual void Print(std::ostream &os) const
Print method.
#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...
const Epetra_Map & ColMap() const
Returns the Epetra_Map object that describes the set of column-indices that appear in each processor&#39;...
int Drumm1(const Epetra_Map &map, bool verbose)
virtual int MyPID() const =0
Return my process ID.
int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
override base-class Epetra_CrsMatrix::InsertGlobalValues method
int PutScalar(double ScalarConstant)
Initialize all values in the matrix with constant value.
const Epetra_Map & RowMap() const
Returns the Epetra_Map object associated with the rows of this matrix.
Epetra Finite-Element CrsMatrix.
Epetra_Comm: The Epetra Communication Abstract Base Class.
Definition: Epetra_Comm.h:73
int ExtractGlobalRowCopy(int GlobalRow, int Length, int &NumEntries, double *Values, int *Indices) const
Returns a copy of the specified global row in user-provided arrays.
int LCID(int GCID_in) const
Returns the local column index for given global column index, returns -1 if no local column for this ...
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 ...
virtual void Print(std::ostream &os) const
Print method.
int NumGlobalNonzeros() const
Returns the number of nonzero entries in the global matrix.
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 NumGlobalCols() const
Returns the number of global matrix columns.
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.
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
matrix_data is a very simple data source to be used for filling test matrices.
int SumIntoGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
override base-class Epetra_CrsMatrix::SumIntoGlobalValues method
virtual int NumProc() const =0
Returns total number of processes.
Epetra Finite-Element Vector.
int Drumm3(const Epetra_Map &map, bool verbose)
int submatrix_formats(const Epetra_Comm &Comm, bool verbose)
int rectangular(const Epetra_Comm &Comm, bool verbose)
int Drumm2(const Epetra_Map &map, bool verbose)
Epetra_CrsGraph: A class for constructing and using sparse compressed row graphs. ...
int four_quads(const Epetra_Comm &Comm, bool preconstruct_graph, 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...