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