Epetra Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
FEVbrMatrix/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_Map.h"
48 #include "Epetra_FEVbrMatrix.h"
49 #include <Epetra_FEVector.h>
50 
51 #include "../epetra_test_err.h"
52 
53 
54 int quad1(const Epetra_Map& map, bool verbose)
55 {
56  const Epetra_Comm & Comm = map.Comm();
57  int numProcs = Comm.NumProc();
58  int localProc = Comm.MyPID();
59 
60  Comm.Barrier();
61  if (verbose && localProc == 0) {
62  cout << "====================== quad1 =============================="<<endl;
63  }
64 
65  //Set up a simple finite-element mesh containing 2-D quad elements, 1 per proc.
66  //
67  // *-----*-----*-----*
68  // 0| 2| 4| 6|
69  // | 0 | 1 | np-1|
70  // | | | |
71  // *-----*-----*-----*
72  // 1 3 5 7
73  //
74  //In the above drawing, 'np' means num-procs. node-numbers are to the
75  //lower-left of each node (*).
76  //
77  //Each processor owns element 'localProc', and each processor owns
78  //nodes 'localProc*2' and 'localProc*2+1' except for the last processor,
79  //which also owns the last two nodes.
80  //
81  //There will be 3 degrees-of-freedom per node, so each element-matrix is
82  //of size 12x12. (4 nodes per element, X 3 dof per node)
83  //
84 
85  int myFirstNode = localProc*2;
86  int myLastNode = localProc*2+1;
87  if (localProc == numProcs-1) {
88  myLastNode += 2;
89  }
90 
91  int numMyNodes = myLastNode - myFirstNode + 1;
92  int* myNodes = new int[numMyNodes];
93  int i, j;
94  int ierr = 0;
95  for(i=0; i<numMyNodes; ++i) {
96  myNodes[i] = myFirstNode + i;
97  }
98 
99  int dofPerNode = 3; //degrees-of-freedom per node
100  int indexBase = 0;
101 
102  Epetra_BlockMap blkMap(-1, numMyNodes, myNodes, dofPerNode,
103  indexBase, Comm);
104 
105  int rowLengths = 3; //each element-matrix will have 4 block-columns.
106  //the rows of the assembled matrix will be longer than
107  //this, but we don't need to worry about that because the
108  //VbrMatrix will add memory as needed. For a real
109  //application where efficiency is a concern, better
110  //performance would be obtained by giving more accurate
111  //row-lengths here.
112 
113  Epetra_FEVbrMatrix A(Copy, blkMap, rowLengths);
114 
115  int nodesPerElem = 4;
116  int* elemNodes = new int[nodesPerElem];
117  for(i=0; i<nodesPerElem; ++i) elemNodes[i] = myFirstNode+i;
118 
119  int elemMatrixDim = nodesPerElem*dofPerNode;
120  int len = elemMatrixDim*elemMatrixDim;
121  double* elemMatrix = new double[len];
122 
123  //In an actual finite-element problem, we would calculate and fill
124  //meaningful element stiffness matrices. But for this simple matrix assembly
125  //test, we're just going to fill our element matrix with 1.0's. This will
126  //make it easy to see whether the matrix is correct after it's assembled.
127 
128  for(i=0; i<len; ++i) elemMatrix[i] = 1.0;
129 
130  //For filling in the matrix block-entries, we would ordinarily have to
131  //carefully copy, or set pointers to, appropriate sections of the
132  //element-matrix. But for this simple case we know that the element-matrix
133  //is all 1's, so we'll just set our block-entry pointer to point to the
134  //beginning of the element-matrix and leave it at that.
135  //Note that the matrix class will refer to dofPerNode X dofPerNode (==9)
136  //positions in the memory pointed to by 'blockEntry'.
137 
138  double* blockEntry = elemMatrix;
139 
140  //The element-matrix is a 4x4 (nodesPerElem X nodesPerElem) matrix of
141  //3x3 block-entries. We'll now load our element-matrix into the global
142  //matrix by looping over it and loading block-entries individually.
143 
144  for(i=0; i<nodesPerElem; ++i) {
145  int blkrow = myFirstNode+i;
146  EPETRA_TEST_ERR( A.BeginInsertGlobalValues(blkrow, nodesPerElem, elemNodes),
147  ierr);
148 
149  for(j=0; j<nodesPerElem; ++j) {
150  for(int ii=0; ii<dofPerNode*dofPerNode; ii++) {
151  blockEntry[ii] = blkrow+elemNodes[j];
152  }
153  EPETRA_TEST_ERR( A.SubmitBlockEntry( blockEntry, dofPerNode,
154  dofPerNode, dofPerNode), ierr);
155  }
156 
157  int err = A.EndSubmitEntries();
158  if (err < 0) {
159  cout << "quad1: error in A.EndSubmitEntries: "<<err<<endl;
160  return(-1);
161  }
162  }
163 
164  EPETRA_TEST_ERR( A.GlobalAssemble(), ierr);
165 
166  if (verbose && localProc==0) {
167  cout << "after globalAssemble"<<endl;
168  }
169  if (verbose) {
170  A.Print(cout);
171  }
172 
173  int numMyRows = A.NumMyRows();
174  int correct_numMyRows = dofPerNode*numMyNodes;
175 
176  if (numMyRows != correct_numMyRows) {
177  cout << "proc " << localProc << ", numMyRows("<<numMyRows<<") doesn't match"
178  << " correct_numMyRows("<<correct_numMyRows<<")."<<endl;
179  return(-1);
180  }
181 
182  int numMyNonzeros = A.NumMyNonzeros();
183  int correct_numMyNonzeros = nodesPerElem*nodesPerElem*dofPerNode*dofPerNode;
184 
185  if (numProcs > 1) {
186  if (localProc == numProcs-1) {
187  correct_numMyNonzeros += dofPerNode*dofPerNode*4;
188  }
189  else if (localProc > 0) {
190  correct_numMyNonzeros -= dofPerNode*dofPerNode*4;
191  }
192  else {
193  //localProc==0 && numProcs > 1
194  correct_numMyNonzeros -= dofPerNode*dofPerNode*8;
195  }
196  }
197 
198  if (numMyNonzeros != correct_numMyNonzeros) {
199  cout << "proc " << localProc << ", numMyNonzeros(" << numMyNonzeros
200  <<") != correct_numMyNonzeros("<<correct_numMyNonzeros<<")"<<endl;
201  return(-1);
202  }
203 
204 
205  delete [] elemMatrix;
206  delete [] myNodes;
207  delete [] elemNodes;
208 
209  Comm.Barrier();
210 
211  return(0);
212 }
213 
214 int quad2(const Epetra_Map& map, bool verbose)
215 {
216  const Epetra_Comm & Comm = map.Comm();
217  int numProcs = Comm.NumProc();
218  int localProc = Comm.MyPID();
219  Comm.Barrier();
220  if (verbose && localProc == 0) {
221  cout << "====================== quad2 =============================="<<endl;
222  }
223 
224  //Set up a simple finite-element mesh containing 2-D quad elements, 2 per proc.
225  //(This test is similar to quad1() above, except for having 2 elements-per-proc
226  // rather than 1.)
227  //
228  // *-----*-----*-----*-------*
229  // 0| 2| 4| 6| 8|
230  // | 0 | 1 | 2 | 2*np-1|
231  // | | | | |
232  // *-----*-----*-----*-------*
233  // 1 3 5 7 9
234  //
235  //In the above drawing, 'np' means num-procs. node-numbers are to the
236  //lower-left of each node (*).
237  //
238  //Each processor owns element 'localProc' and 'localProc+1', and each processor
239  //owns nodes 'localProc*4' through 'localProc*4+3' except for the last
240  //processor, which also owns the last two nodes in the mesh.
241  //
242  //There will be 3 degrees-of-freedom per node, so each element-matrix is
243  //of size 12x12.
244  //
245 
246  int myFirstNode = localProc*4;
247  int myLastNode = localProc*4+3;
248  if (localProc == numProcs-1) {
249  myLastNode += 2;
250  }
251 
252  int numMyElems = 2;
253 
254  int numMyNodes = myLastNode - myFirstNode + 1;
255  int* myNodes = new int[numMyNodes];
256  int i, j;
257  int ierr = 0;
258  for(i=0; i<numMyNodes; ++i) {
259  myNodes[i] = myFirstNode + i;
260  }
261 
262  int dofPerNode = 3; //degrees-of-freedom per node
263  int indexBase = 0;
264 
265  Epetra_BlockMap blkMap(-1, numMyNodes, myNodes, dofPerNode,
266  indexBase, Comm);
267 
268  int rowLengths = 4; //each element-matrix will have 4 block-columns.
269  //the rows of the assembled matrix will be longer than
270  //this, but we don't need to worry about that because the
271  //VbrMatrix will add memory as needed. For a real
272  //application where efficiency is a concern, better
273  //performance would be obtained by giving a more accurate
274  //row-length here.
275 
276  Epetra_FEVbrMatrix A(Copy, blkMap, rowLengths);
277 
278  int nodesPerElem = 4;
279  int* elemNodes = new int[nodesPerElem];
280 
281  int elemMatrixDim = nodesPerElem*dofPerNode;
282  int len = elemMatrixDim*elemMatrixDim;
283  double* elemMatrix = new double[len];
284 
285  //In an actual finite-element problem, we would calculate and fill
286  //meaningful element stiffness matrices. But for this simple matrix assembly
287  //test, we're just going to fill our element matrix with 1.0's. This will
288  //make it easy to see whether the matrix is correct after it's assembled.
289 
290  for(i=0; i<len; ++i) elemMatrix[i] = 1.0;
291 
292  //For filling in the matrix block-entries, we would ordinarily have to
293  //carefully copy, or set pointers to, appropriate sections of the
294  //element-matrix. But for this simple case we know that the element-matrix
295  //is all 1's, so we'll just set our block-entry pointer to point to the
296  //beginning of the element-matrix and leave it at that.
297  //Note that the matrix class will refer to dofPerNode X dofPerNode (==9)
298  //positions in the memory pointed to by 'blockEntry'.
299 
300  double* blockEntry = elemMatrix;
301 
302  //Each element-matrix is a 4x4 (nodesPerElem X nodesPerElem) matrix of
303  //3x3 block-entries. We'll now load our element-matrices into the global
304  //matrix by looping over them and loading block-entries individually.
305 
306  int firstNode = myFirstNode;
307  for(int el=0; el<numMyElems; ++el) {
308  for(i=0; i<nodesPerElem; ++i) {
309 
310  for(int n=0; n<nodesPerElem; ++n) elemNodes[n] = firstNode+n;
311 
312  int blkrow = firstNode+i;
313  EPETRA_TEST_ERR( A.BeginInsertGlobalValues(blkrow, nodesPerElem, elemNodes),
314  ierr);
315 
316  for(j=0; j<nodesPerElem; ++j) {
317  EPETRA_TEST_ERR( A.SubmitBlockEntry( blockEntry, dofPerNode,
318  dofPerNode, dofPerNode), ierr);
319  }
320 
321  int this_err = A.EndSubmitEntries();
322  if (this_err < 0) {
323  cerr << "error in quad2, A.EndSubmitEntries(): " << this_err << endl;
324  return(this_err);
325  }
326  }
327 
328  firstNode += 2;
329  }
330 
331  EPETRA_TEST_ERR( A.GlobalAssemble(), ierr);
332 
333  if (verbose && localProc==0) {
334  cout << "after globalAssemble"<<endl;
335  }
336 
337  if (verbose) {
338  A.Print(cout);
339  }
340 
341  //now let's make sure that we can perform a matvec...
342  Epetra_FEVector x(blkMap, 1), y(blkMap, 1);
343  x.PutScalar(1.0);
344 
345  EPETRA_TEST_ERR( A.Multiply(false, x, y), ierr);
346 
347  if (verbose && localProc==0) {
348  cout << "quad2, y:"<<endl;
349  }
350  if (verbose) {
351  y.Print(cout);
352  }
353 
354  delete [] elemMatrix;
355  delete [] myNodes;
356  delete [] elemNodes;
357 
358  return(0);
359 }
360 
361 int MultiVectorTests(const Epetra_Map & Map, int NumVectors, bool verbose)
362 {
363  const Epetra_Comm & Comm = Map.Comm();
364  int ierr = 0, i, j;
365 
366  /* get number of processors and the name of this processor */
367 
368  int MyPID = Comm.MyPID();
369 
370  // Construct FEVbrMatrix
371 
372  if (verbose && MyPID==0) cout << "constructing Epetra_FEVbrMatrix" << endl;
373 
374  //
375  //we'll set up a tri-diagonal matrix.
376  //
377 
378  int numGlobalRows = Map.NumGlobalElements();
379  int minLocalRow = Map.MinMyGID();
380  int rowLengths = 3;
381 
382  Epetra_FEVbrMatrix A(Copy, Map, rowLengths);
383 
384  if (verbose && MyPID==0) {
385  cout << "calling A.InsertGlobalValues with 1-D data array"<<endl;
386  }
387 
388  int numCols = 3;
389  int* ptIndices = new int[numCols];
390  for(int k=0; k<numCols; ++k) {
391  ptIndices[k] = minLocalRow+k;
392  }
393 
394  double* values_1d = new double[numCols*numCols];
395  for(j=0; j<numCols*numCols; ++j) {
396  values_1d[j] = 3.0;
397  }
398 
399  //For an extreme test, we'll have all processors sum into all rows.
400 
401  int minGID = Map.MinAllGID();
402 
403  //For now we're going to assume that there's just one point associated with
404  //each GID (element).
405 
406  double* ptCoefs = new double[3];
407 
408  {for(i=0; i<numGlobalRows; ++i) {
409  if (i>0 && i<numGlobalRows-1) {
410  ptIndices[0] = minGID+i-1;
411  ptIndices[1] = minGID+i;
412  ptIndices[2] = minGID+i+1;
413  ptCoefs[0] = -1.0;
414  ptCoefs[1] = 2.0;
415  ptCoefs[2] = -1.0;
416  numCols = 3;
417  }
418  else if (i == 0) {
419  ptIndices[0] = minGID+i;
420  ptIndices[1] = minGID+i+1;
421  ptIndices[2] = minGID+i+2;
422  ptCoefs[0] = 2.0;
423  ptCoefs[1] = -1.0;
424  ptCoefs[2] = -1.0;
425  numCols = 3;
426  }
427  else {
428  ptIndices[0] = minGID+i-2;
429  ptIndices[1] = minGID+i-1;
430  ptIndices[2] = minGID+i;
431  ptCoefs[0] = -1.0;
432  ptCoefs[1] = -1.0;
433  ptCoefs[2] = 2.0;
434  numCols = 3;
435  }
436 
437  int row = minGID+i;
438 
439  EPETRA_TEST_ERR( A.BeginInsertGlobalValues(row, rowLengths, ptIndices), ierr);
440 
441  for(j=0; j<rowLengths; ++j) {
442  EPETRA_TEST_ERR( A.SubmitBlockEntry(&(ptCoefs[j]), 1, 1, 1), ierr);
443  }
444 
445  EPETRA_TEST_ERR( A.EndSubmitEntries(), ierr);
446 
447  }}
448 
449  if (verbose&&MyPID==0) {
450  cout << "calling A.GlobalAssemble()" << endl;
451  }
452 
453  EPETRA_TEST_ERR( A.GlobalAssemble(), ierr );
454 
455  if (verbose&&MyPID==0) {
456  cout << "after globalAssemble"<<endl;
457  }
458  if (verbose) {
459  A.Print(cout);
460  }
461 
462  delete [] values_1d;
463  delete [] ptIndices;
464  delete [] ptCoefs;
465 
466  return(ierr);
467 }
468 
469 int four_quads(const Epetra_Comm& Comm, bool preconstruct_graph, bool verbose)
470 {
471  if (verbose) {
472  cout << "******************* four_quads ***********************"<<endl;
473  }
474 
475  //This function assembles a matrix representing a finite-element mesh
476  //of four 2-D quad elements. There are 9 nodes in the problem. The
477  //same problem is assembled no matter how many processors are being used
478  //(within reason). It may not work if more than 9 processors are used.
479  //
480  // *------*------*
481  // 6| 7| 8|
482  // | E2 | E3 |
483  // *------*------*
484  // 3| 4| 5|
485  // | E0 | E1 |
486  // *------*------*
487  // 0 1 2
488  //
489  //Nodes are denoted by * with node-numbers below and left of each node.
490  //E0, E1 and so on are element-numbers.
491  //
492  //Each processor will contribute a sub-matrix of size 4x4, filled with 1's,
493  //for each element. Thus, the coefficient value at position 0,0 should end up
494  //being 1.0*numProcs, the value at position 4,4 should be 1.0*4*numProcs, etc.
495  //
496  //Depending on the number of processors being used, the locations of the
497  //specific matrix positions (in terms of which processor owns them) will vary.
498  //
499 
500  int numProcs = Comm.NumProc();
501 
502  int numNodes = 9;
503  int numElems = 4;
504  int numNodesPerElem = 4;
505 
506  int blockSize = 1;
507  int indexBase = 0;
508 
509  //Create a map using epetra-defined linear distribution.
510  Epetra_BlockMap map(numNodes, blockSize, indexBase, Comm);
511 
512  Epetra_CrsGraph* graph = NULL;
513 
514  int* nodes = new int[numNodesPerElem];
515  int i, j, k, err = 0;
516 
517  if (preconstruct_graph) {
518  graph = new Epetra_CrsGraph(Copy, map, 1);
519 
520  //we're going to fill the graph with indices, but remember it will only
521  //accept indices in rows for which map.MyGID(row) is true.
522 
523  for(i=0; i<numElems; ++i) {
524  switch(i) {
525  case 0:
526  nodes[0] = 0; nodes[1] = 1; nodes[2] = 4; nodes[3] = 3;
527  break;
528  case 1:
529  nodes[0] = 1; nodes[1] = 2; nodes[2] = 5; nodes[3] = 4;
530  break;
531  case 2:
532  nodes[0] = 3; nodes[1] = 4; nodes[2] = 7; nodes[3] = 6;
533  break;
534  case 3:
535  nodes[0] = 4; nodes[1] = 5; nodes[2] = 8; nodes[3] = 7;
536  break;
537  }
538 
539  for(j=0; j<numNodesPerElem; ++j) {
540  if (map.MyGID(nodes[j])) {
541  err = graph->InsertGlobalIndices(nodes[j], numNodesPerElem,
542  nodes);
543  if (err<0) return(err);
544  }
545  }
546  }
547 
548  EPETRA_CHK_ERR( graph->FillComplete() );
549  }
550 
551  Epetra_FEVbrMatrix* A = NULL;
552 
553  if (preconstruct_graph) {
554  A = new Epetra_FEVbrMatrix(Copy, *graph);
555  }
556  else {
557  A = new Epetra_FEVbrMatrix(Copy, map, 1);
558  }
559 
560  //EPETRA_CHK_ERR( A->PutScalar(0.0) );
561 
562  double* values_1d = new double[numNodesPerElem*numNodesPerElem];
563  double** values_2d = new double*[numNodesPerElem];
564 
565  for(i=0; i<numNodesPerElem*numNodesPerElem; ++i) values_1d[i] = 1.0;
566 
567  int offset = 0;
568  for(i=0; i<numNodesPerElem; ++i) {
569  values_2d[i] = &(values_1d[offset]);
570  offset += numNodesPerElem;
571  }
572 
573  for(i=0; i<numElems; ++i) {
574  switch(i) {
575  case 0:
576  nodes[0] = 0; nodes[1] = 1; nodes[2] = 4; nodes[3] = 3;
577  break;
578 
579  case 1:
580  nodes[0] = 1; nodes[1] = 2; nodes[2] = 5; nodes[3] = 4;
581  break;
582 
583  case 2:
584  nodes[0] = 3; nodes[1] = 4; nodes[2] = 7; nodes[3] = 6;
585  break;
586 
587  case 3:
588  nodes[0] = 4; nodes[1] = 5; nodes[2] = 8; nodes[3] = 7;
589  break;
590  }
591 
592  for(j=0; j<numNodesPerElem; ++j) {
593  if (preconstruct_graph) {
594  err = A->BeginSumIntoGlobalValues(nodes[j], numNodesPerElem, nodes);
595  if (err<0) return(err);
596  }
597  else {
598  err = A->BeginInsertGlobalValues(nodes[j], numNodesPerElem, nodes);
599  if (err<0) return(err);
600  }
601 
602  for(k=0; k<numNodesPerElem; ++k) {
603  err = A->SubmitBlockEntry(values_1d, blockSize, blockSize, blockSize);
604  if (err<0) return(err);
605  }
606 
607  err = A->EndSubmitEntries();
608  if (err<0) return(err);
609  }
610  }
611 
613 
614  Epetra_FEVbrMatrix* Acopy = new Epetra_FEVbrMatrix(*A);
615 
616  if (verbose) {
617  cout << "A:"<<*A << endl;
618  cout << "Acopy:"<<*Acopy<<endl;
619  }
620 
621  Epetra_Vector x(A->RowMap()), y(A->RowMap());
622 
623  x.PutScalar(1.0); y.PutScalar(0.0);
624 
625  Epetra_Vector x2(Acopy->RowMap()), y2(Acopy->RowMap());
626 
627  x2.PutScalar(1.0); y2.PutScalar(0.0);
628 
629  A->Multiply(false, x, y);
630 
631  Acopy->Multiply(false, x2, y2);
632 
633  double ynorm2, y2norm2;
634 
635  y.Norm2(&ynorm2);
636  y2.Norm2(&y2norm2);
637  if (ynorm2 != y2norm2) {
638  cerr << "norm2(A*ones) != norm2(*Acopy*ones)"<<endl;
639  return(-99);
640  }
641 
642  Epetra_FEVbrMatrix* Acopy2 =
643  new Epetra_FEVbrMatrix(Copy, A->RowMap(), A->ColMap(), 1);
644 
645  *Acopy2 = *Acopy;
646 
647  Epetra_Vector x3(Acopy->RowMap()), y3(Acopy->RowMap());
648 
649  x3.PutScalar(1.0); y3.PutScalar(0.0);
650 
651  Acopy2->Multiply(false, x3, y3);
652 
653  double y3norm2;
654  y3.Norm2(&y3norm2);
655 
656  if (y3norm2 != y2norm2) {
657  cerr << "norm2(Acopy*ones) != norm2(Acopy2*ones)"<<endl;
658  return(-999);
659  }
660 
661  int len = 20;
662  int* indices = new int[len];
663  double* values = new double[len];
664  int numIndices;
665 
666  if (map.MyGID(0)) {
667  int lid = map.LID(0);
668  EPETRA_CHK_ERR( A->ExtractMyRowCopy(lid, len, numIndices,
669  values, indices) );
670  if (numIndices != 4) {
671  return(-1);
672  }
673  if (indices[0] != lid) {
674  return(-2);
675  }
676 
677  if (values[0] != 1.0*numProcs) {
678  cout << "ERROR: values[0] ("<<values[0]<<") should be "<<numProcs<<endl;
679  return(-3);
680  }
681  }
682 
683  if (map.MyGID(4)) {
684  int lid = map.LID(4);
685  EPETRA_CHK_ERR( A->ExtractMyRowCopy(lid, len, numIndices,
686  values, indices) );
687 
688  if (numIndices != 9) {
689  return(-4);
690  }
691  int lcid = A->LCID(4);
692 // if (indices[lcid] != 4) {
693 // cout << "ERROR: indices[4] ("<<indices[4]<<") should be "
694 // <<A->LCID(4)<<endl;
695 // return(-5);
696 // }
697  if (values[lcid] != 4.0*numProcs) {
698  cout << "ERROR: values["<<lcid<<"] ("<<values[lcid]<<") should be "
699  <<4*numProcs<<endl;
700  return(-6);
701  }
702  }
703 
704  delete [] values_2d;
705  delete [] values_1d;
706  delete [] nodes;
707  delete [] indices;
708  delete [] values;
709 
710  delete A;
711  delete Acopy2;
712  delete Acopy;
713  delete graph;
714 
715  return(0);
716 }
int NumGlobalElements() const
Number of elements across all processors.
Epetra_Map: A class for partitioning vectors and matrices.
Definition: Epetra_Map.h:127
int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const
Returns a copy of the specified local row in user-provided arrays.
int Multiply(bool TransA, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Epetra_VbrMatrix multiplied by a Epetra_MultiVector X in Y.
int FillComplete()
Tranform to local index space. Perform other operations to allow optimal matrix operations.
int quad1(const Epetra_Map &map, bool verbose)
int GlobalAssemble(bool callFillComplete=true)
#define EPETRA_TEST_ERR(a, b)
int PutScalar(double ScalarConstant)
Initialize all values in a multi-vector with constant value.
const Epetra_BlockMap & ColMap() const
Returns the ColMap as an Epetra_BlockMap (the Epetra_Map base class) needed for implementing Epetra_R...
virtual void Print(std::ostream &os) const
Print method.
#define EPETRA_CHK_ERR(a)
int NumMyRows() const
Returns the number of matrix rows owned by the calling processor.
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.
virtual void Barrier() const =0
Epetra_Comm Barrier function.
virtual int MyPID() const =0
Return my process ID.
int quad2(const Epetra_Map &map, bool verbose)
int MultiVectorTests(const Epetra_Map &Map, int NumVectors, bool verbose)
Epetra_Comm: The Epetra Communication Abstract Base Class.
Definition: Epetra_Comm.h:81
int NumMyNonzeros() const
Returns the number of nonzero entriesowned by the calling processor .
int BeginSumIntoGlobalValues(int BlockRow, int NumBlockEntries, int *BlockIndices)
Initiate summing into current values with this list of entries for a given global row of the matrix...
Epetra_BlockMap: A class for partitioning block element vectors and matrices.
const Epetra_BlockMap & RowMap() const
Returns the RowMap object as an Epetra_BlockMap (the Epetra_Map base class) needed for implementing E...
int EndSubmitEntries()
Completes processing of all data passed in for the current block row.
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 InsertGlobalIndices(int GlobalRow, int NumIndices, int *Indices)
Enter a list of elements in a specified global row of the graph.
int MinMyGID() const
Returns the minimum global ID owned by this processor.
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
int LCID(int GCID_in) const
Returns the local column index for given global column index, returns -1 if no local column for this ...
Epetra Finite-Element VbrMatrix.
virtual int NumProc() const =0
Returns total number of processes.
int SubmitBlockEntry(double *Values, int LDA, int NumRows, int NumCols)
Submit a block entry to the indicated block row and column specified in the Begin routine...
Epetra Finite-Element Vector.
virtual void Print(std::ostream &os) const
Print method.
int n
Epetra_CrsGraph: A class for constructing and using sparse compressed row graphs. ...
int BeginInsertGlobalValues(int BlockRow, int NumBlockEntries, int *BlockIndices)
Initiate insertion of a list of elements in a given global row of the matrix, values are inserted via...
int four_quads(const Epetra_Comm &Comm, bool preconstruct_graph, bool verbose)