Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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
54int 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
180int 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
307int 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
425int 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
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
684int 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
776int 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
#define EPETRA_CHK_ERR(a)
@ View
@ Copy
int Drumm3(const Epetra_Map &map, bool verbose)
int four_quads(const Epetra_Comm &Comm, bool preconstruct_graph, bool verbose)
int Drumm1(const Epetra_Map &map, bool verbose)
int rectangular(const Epetra_Comm &Comm, bool verbose)
int submatrix_formats(const Epetra_Comm &Comm, bool verbose)
int Drumm2(const Epetra_Map &map, bool verbose)
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
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
virtual int NumProc() const =0
Returns total number of processes.
virtual int MyPID() const =0
Return my process ID.
Epetra_CrsGraph: A class for constructing and using sparse compressed row graphs.
int InsertGlobalIndices(int GlobalRow, int NumIndices, int *Indices)
Enter a list of elements in a specified global row of the graph.
int FillComplete()
Tranform to local index space. Perform other operations to allow optimal matrix operations.
const Epetra_Map & RowMap() const
Returns the Epetra_Map object associated with the rows of this matrix.
int NumGlobalCols() const
Returns the number of global matrix columns.
int PutScalar(double ScalarConstant)
Initialize all values in the matrix with constant value.
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 NumGlobalNonzeros() const
Returns the number of nonzero entries in the global matrix.
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.
const Epetra_Map & ColMap() const
Returns the Epetra_Map object that describes the set of column-indices that appear in each processor'...
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.
Epetra Finite-Element CrsMatrix.
int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
override base-class Epetra_CrsMatrix::InsertGlobalValues method
int SumIntoGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
override base-class Epetra_CrsMatrix::SumIntoGlobalValues method
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 ...
void Print(std::ostream &os) const
Print method.
Epetra Finite-Element Vector.
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...
int GlobalAssemble(Epetra_CombineMode mode=Add, bool reuse_map_and_exporter=false)
Gather any overlapping/shared data into the non-overlapping partitioning defined by the Map that was ...
Epetra_IntSerialDenseVector: A class for constructing and using dense vectors.
Epetra_Map: A class for partitioning vectors and matrices.
Definition: Epetra_Map.h:119
virtual void Print(std::ostream &os) const
Print method.
int PutScalar(double ScalarConstant)
Initialize all values in a multi-vector with constant value.
Epetra_SerialDenseMatrix: A class for constructing and using real double precision general dense matr...
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.
matrix_data is a very simple data source to be used for filling test matrices.
bool compare_local_data(const Epetra_CrsMatrix &A)
The portion of this matrix_data object's data that corresponds to the locally-owned rows of A,...
#define EPETRA_TEST_ERR(a, b)
bool compare_matrices(const Epetra_CrsMatrix &A, const Epetra_CrsMatrix &B)
Check whether the two CrsMatrix arguments have the same size, structure and coefs.