Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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"
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 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
181int 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
309int 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
427int 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
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
686int 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
778int 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
#define EPETRA_CHK_ERR(a)
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 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 ...
long long NumGlobalNonzeros64() const
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.
long long NumGlobalCols64() const
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_LongLongSerialDenseVector: 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.