FECrsGraph/ExecuteTestProblems.cpp

Go to the documentation of this file.
00001 //@HEADER
00002 // ************************************************************************
00003 // 
00004 //               Epetra: Linear Algebra Services Package 
00005 //                 Copyright (2001) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ************************************************************************
00027 //@HEADER
00028 
00029 #include "Epetra_BLAS.h"
00030 #include "ExecuteTestProblems.h"
00031 #include "Epetra_Comm.h"
00032 #include "Epetra_Map.h"
00033 #include "Epetra_Vector.h"
00034 #include "Epetra_FECrsGraph.h"
00035 #include "Epetra_FECrsMatrix.h"
00036 #include "Epetra_IntSerialDenseVector.h"
00037 #include "Epetra_SerialDenseMatrix.h"
00038 
00039 
00040 int Drumm1(const Epetra_Map& map, bool verbose)
00041 {
00042   //Simple 2-element problem (element as in "finite-element") from
00043   //Clif Drumm. Two triangular elements, one per processor, as shown
00044   //here:
00045   //
00046   //   *----*
00047   //  3|\  2|
00048   //   | \  |
00049   //   | 0\1|
00050   //   |   \|
00051   //   *----*
00052   //  0    1
00053   //
00054   //Element 0 on processor 0, element 1 on processor 1.
00055   //Processor 0 will own nodes 0,1 and processor 1 will own nodes 2,3.
00056   //Each processor will pass a 3x3 element-connectivity-matrix to
00057   //Epetra_FECrsGraph.
00058   //After GlobalAssemble(), the graph should be as follows:
00059   //
00060   //         row 0: 2  1  0  1
00061   //proc 0   row 1: 1  4  1  2
00062   //----------------------------------
00063   //         row 2: 0  1  2  1
00064   //proc 1   row 3: 1  2  1  4
00065   //
00066 
00067   int numProcs = map.Comm().NumProc();
00068   int localProc = map.Comm().MyPID();
00069 
00070   if (numProcs != 2) return(0);
00071 
00072   int indexBase = 0, ierr = 0;
00073 
00074   int numMyNodes = 2;
00075   int* myNodes = new int[numMyNodes];
00076 
00077   if (localProc == 0) {
00078     myNodes[0] = 0;
00079     myNodes[1] = 1;
00080   }
00081   else {
00082     myNodes[0] = 2;
00083     myNodes[1] = 3;
00084   }
00085 
00086   Epetra_Map Map(-1, numMyNodes, myNodes, indexBase, map.Comm());
00087 
00088   delete [] myNodes;
00089   numMyNodes = 3;
00090   myNodes = new int[numMyNodes];
00091 
00092   if (localProc == 0) {
00093     myNodes[0] = 0;
00094     myNodes[1] = 1;
00095     myNodes[2] = 3;
00096   }
00097   else {
00098     myNodes[0] = 1;
00099     myNodes[1] = 2;
00100     myNodes[2] = 3;
00101   }
00102 
00103   int rowLengths = 3;
00104   Epetra_FECrsGraph A(Copy, Map, rowLengths);
00105 
00106   EPETRA_TEST_ERR( A.InsertGlobalIndices(numMyNodes, myNodes,
00107            numMyNodes, myNodes),ierr);
00108 
00109   EPETRA_TEST_ERR( A.GlobalAssemble(), ierr );
00110 
00111   if (verbose) {
00112     A.Print(std::cout);
00113   }
00114 
00115   delete [] myNodes;
00116 
00117   return(0);
00118 }
00119 
00120 int Drumm2(const Epetra_Map& map, bool verbose)
00121 {
00122   //Simple 2-element problem (element as in "finite-element") from
00123   //Clif Drumm. Two triangular elements, one per processor, as shown
00124   //here:
00125   //
00126   //   *----*
00127   //  3|\  2|
00128   //   | \  |
00129   //   | 0\1|
00130   //   |   \|
00131   //   *----*
00132   //  0    1
00133   //
00134   //Element 0 on processor 0, element 1 on processor 1.
00135   //Processor 0 will own nodes 0,1,3 and processor 1 will own node 2.
00136   //Each processor will pass a 3x3 element-connectivity-matrix to
00137   //Epetra_FECrsGraph.
00138   //After GlobalAssemble(), the graph should be as follows:
00139   //
00140   //         row 0: 2  1  0  1
00141   //proc 0   row 1: 1  4  1  2
00142   //         row 2: 0  1  2  1
00143   //----------------------------------
00144   //proc 1   row 3: 1  2  1  4
00145   //
00146 
00147   int numProcs = map.Comm().NumProc();
00148   int localProc = map.Comm().MyPID();
00149 
00150   if (numProcs != 2) return(0);
00151 
00152   int indexBase = 0, ierr = 0;
00153   int numMyNodes = 3;
00154   int* myNodes = new int[numMyNodes];
00155 
00156   if (localProc == 0) {
00157     myNodes[0] = 0;
00158     myNodes[1] = 1;
00159     myNodes[2] = 3;
00160   }
00161   else {
00162     numMyNodes = 1;
00163     myNodes[0] = 2;
00164   }
00165 
00166   Epetra_Map Map(-1, numMyNodes, myNodes, indexBase, map.Comm());
00167 
00168   int rowLengths = 3;
00169   Epetra_FECrsGraph A(Copy, Map, rowLengths);
00170 
00171   if (localProc != 0) {
00172     numMyNodes = 3;
00173     myNodes[0] = 1;
00174     myNodes[1] = 2;
00175     myNodes[2] = 3;
00176   }
00177 
00178   EPETRA_TEST_ERR( A.InsertGlobalIndices(numMyNodes, myNodes,
00179            numMyNodes, myNodes),ierr);
00180 
00181   EPETRA_TEST_ERR( A.GlobalAssemble(), ierr );
00182 
00183   if (verbose) {
00184     A.Print(std::cout);
00185   }
00186 
00187   delete [] myNodes;
00188 
00189   return(0);
00190 }
00191 
00192 int four_quads(const Epetra_Comm& Comm, bool preconstruct_graph, bool verbose)
00193 {
00194   if (verbose) {
00195     std::cout << "******************* four_quads ***********************"<<std::endl;
00196   }
00197 
00198   //This function assembles a matrix representing a finite-element
00199   //mesh of four 2-D quad elements. There are 9 nodes in the problem. The
00200   //same problem is assembled no matter how many processors are being used
00201   //(within reason). It may not work if more than 9 processors are used.
00202   //
00203   //  *------*------*
00204   // 6|     7|     8|
00205   //  | E2   | E3   |
00206   //  *------*------*
00207   // 3|     4|     5|
00208   //  | E0   | E1   |
00209   //  *------*------*
00210   // 0      1      2
00211   //
00212   //Nodes are denoted by * with node-numbers below and left of each node.
00213   //E0, E1 and so on are element-numbers.
00214   //
00215   //Each processor will contribute a sub-matrix of size 4x4, filled with 1's,
00216   //for each element. Thus, the coefficient value at position 0,0 should end up
00217   //being 1.0*numProcs, the value at position 4,4 should be 1.0*4*numProcs, etc.
00218   //
00219   //Depending on the number of processors being used, the locations of the
00220   //specific matrix positions (in terms of which processor owns them) will vary.
00221   //
00222 
00223   int numProcs = Comm.NumProc();
00224 
00225   int numNodes = 9;
00226   int numElems = 4;
00227   int numNodesPerElem = 4;
00228 
00229   int indexBase = 0;
00230 
00231   //Create a map using epetra-defined linear distribution.
00232   Epetra_Map map(numNodes, indexBase, Comm);
00233 
00234   Epetra_FECrsGraph* graph = NULL;
00235 
00236   int* nodes = new int[numNodesPerElem];
00237   int i, err = 0;
00238 
00239   if (preconstruct_graph) {
00240     graph = new Epetra_FECrsGraph(Copy, map, 1);
00241 
00242     //we're going to fill the graph with indices, by passing our
00243     //connectivity lists.
00244     //FECrsGraph should accept indices in all rows, regardless of
00245     //whether map.MyGID(row) is true.
00246 
00247     for(i=0; i<numElems; ++i) {
00248       switch(i) {
00249       case 0:
00250   nodes[0] = 0; nodes[1] = 1; nodes[2] = 4; nodes[3] = 3;
00251   break;
00252       case 1:
00253   nodes[0] = 1; nodes[1] = 2; nodes[2] = 5; nodes[3] = 4;
00254   break;
00255       case 2:
00256   nodes[0] = 3; nodes[1] = 4; nodes[2] = 7; nodes[3] = 6;
00257   break;
00258       case 3:
00259   nodes[0] = 4; nodes[1] = 5; nodes[2] = 8; nodes[3] = 7;
00260   break;
00261       }
00262 
00263       err = graph->InsertGlobalIndices(numNodesPerElem, nodes,
00264                                        numNodesPerElem, nodes);
00265       if (err < 0) {
00266         std::cerr << "ERROR, FECrsGraph error in InsertGlobalIndices, err="
00267           << err << std::endl;
00268         return(-1);
00269       }
00270     }
00271 
00272     EPETRA_CHK_ERR( graph->GlobalAssemble() );
00273   }
00274 
00275   Epetra_FECrsMatrix* A = NULL;
00276 
00277   if (preconstruct_graph) {
00278     A = new Epetra_FECrsMatrix(Copy, *graph);
00279   }
00280   else {
00281     A = new Epetra_FECrsMatrix(Copy, map, 1);
00282   }
00283 
00284   EPETRA_CHK_ERR( A->PutScalar(0.0) );
00285 
00286   double* values_1d = new double[numNodesPerElem*numNodesPerElem];
00287   double** values_2d = new double*[numNodesPerElem];
00288 
00289   for(i=0; i<numNodesPerElem*numNodesPerElem; ++i) values_1d[i] = 1.0;
00290 
00291   int offset = 0;
00292   for(i=0; i<numNodesPerElem; ++i) {
00293     values_2d[i] = &(values_1d[offset]);
00294     offset += numNodesPerElem;
00295   }
00296 
00297   int format = Epetra_FECrsMatrix::ROW_MAJOR;
00298   Epetra_IntSerialDenseVector epetra_nodes(View, nodes, numNodesPerElem);
00299   Epetra_SerialDenseMatrix epetra_values(View, values_1d, numNodesPerElem,
00300            numNodesPerElem, numNodesPerElem);
00301 
00302   for(i=0; i<numElems; ++i) {
00303     switch(i) {
00304     case 0:
00305       nodes[0] = 0; nodes[1] = 1; nodes[2] = 4; nodes[3] = 3;
00306       if (preconstruct_graph) {
00307   err = A->SumIntoGlobalValues(epetra_nodes,
00308              epetra_values, format);
00309   if (err<0) return(err);
00310       }
00311       else {
00312   err = A->InsertGlobalValues(epetra_nodes,
00313             epetra_values, format);
00314   if (err<0) return(err);
00315       }
00316       break;
00317 
00318     case 1:
00319       nodes[0] = 1; nodes[1] = 2; nodes[2] = 5; nodes[3] = 4;
00320       if (preconstruct_graph) {
00321   err = A->SumIntoGlobalValues(numNodesPerElem, nodes,
00322              values_2d, format);
00323   if (err<0) return(err);
00324       }
00325       else {
00326   err = A->InsertGlobalValues(numNodesPerElem, nodes,
00327             values_2d, format);
00328   if (err<0) return(err);
00329       }
00330       break;
00331 
00332     case 2:
00333       nodes[0] = 3; nodes[1] = 4; nodes[2] = 7; nodes[3] = 6;
00334       if (preconstruct_graph) {
00335   err = A->SumIntoGlobalValues(numNodesPerElem, nodes,
00336              numNodesPerElem, nodes,
00337              values_1d, format);
00338   if (err<0) return(err);
00339       }
00340       else {
00341   err = A->InsertGlobalValues(numNodesPerElem, nodes,
00342             numNodesPerElem, nodes,
00343             values_1d, format);
00344   if (err<0) return(err);
00345       }
00346       break;
00347 
00348      case 3:
00349       nodes[0] = 4; nodes[1] = 5; nodes[2] = 8; nodes[3] = 7;
00350       if (preconstruct_graph) {
00351   err = A->SumIntoGlobalValues(numNodesPerElem, nodes,
00352              numNodesPerElem, nodes,
00353              values_2d, format);
00354   if (err<0) return(err);
00355       }
00356       else {
00357   err = A->InsertGlobalValues(numNodesPerElem, nodes,
00358             numNodesPerElem, nodes,
00359             values_2d, format);
00360   if (err<0) return(err);
00361       }
00362       break;
00363     }
00364   }
00365 
00366   err = A->GlobalAssemble();
00367   if (err < 0) {
00368     return(err);
00369   }
00370 
00371   Epetra_Vector x(A->RowMap()), y(A->RowMap());
00372 
00373   x.PutScalar(1.0); y.PutScalar(0.0);
00374 
00375   Epetra_FECrsMatrix Acopy(*A);
00376 
00377   Epetra_Vector x2(Acopy.RowMap()), y2(Acopy.RowMap());
00378 
00379   x2.PutScalar(1.0); y2.PutScalar(0.0);
00380 
00381   A->Multiply(false, x, y);
00382 
00383   Acopy.Multiply(false, x2, y2);
00384 
00385   double ynorm2, y2norm2;
00386 
00387   y.Norm2(&ynorm2);
00388   y2.Norm2(&y2norm2);
00389   if (ynorm2 != y2norm2) {
00390     std::cerr << "norm2(A*ones) != norm2(Acopy*ones)"<<std::endl;
00391     return(-99);
00392   }
00393 
00394   err = Acopy.GlobalAssemble();
00395   if (err < 0) {
00396     return(err);
00397   }
00398 
00399   if (verbose) {
00400     std::cout << "A:"<<std::endl<<*A << std::endl;
00401     std::cout << "Acopy:"<<std::endl<<Acopy << std::endl;
00402   }
00403 
00404   Epetra_FECrsMatrix Acopy2(Copy, A->RowMap(), A->ColMap(), 1);
00405 
00406   Acopy2 = Acopy;
00407 
00408   Epetra_Vector x3(Acopy.RowMap()), y3(Acopy.RowMap());
00409 
00410   x3.PutScalar(1.0); y3.PutScalar(0.0);
00411 
00412   Acopy2.Multiply(false, x3, y3);
00413 
00414   double y3norm2;
00415   y3.Norm2(&y3norm2);
00416 
00417   if (y3norm2 != y2norm2) {
00418     std::cerr << "norm2(Acopy*ones) != norm2(Acopy2*ones)"<<std::endl;
00419     return(-999);
00420   }
00421 
00422   int len = 20;
00423   int* indices = new int[len];
00424   double* values = new double[len];
00425   int numIndices;
00426 
00427   if (map.MyGID(0)) {
00428     EPETRA_CHK_ERR( A->ExtractGlobalRowCopy(0, len, numIndices,
00429               values, indices) );
00430     if (numIndices != 4) {
00431       return(-1);
00432     }
00433     if (indices[0] != 0) {
00434       return(-2);
00435     }
00436 
00437     if (values[0] != 1.0*numProcs) {
00438       std::cout << "ERROR: values[0] ("<<values[0]<<") should be "<<numProcs<<std::endl;
00439       return(-3);
00440     }
00441   }
00442 
00443   if (map.MyGID(4)) {
00444     EPETRA_CHK_ERR( A->ExtractGlobalRowCopy(4, len, numIndices,
00445               values, indices) );
00446 
00447     if (numIndices != 9) {
00448       return(-4);
00449     }
00450     int lcid = A->LCID(4);
00451     if (lcid<0) {
00452       return(-5);
00453     }
00454     if (values[lcid] != 4.0*numProcs) {
00455       std::cout << "ERROR: values["<<lcid<<"] ("<<values[lcid]<<") should be "
00456      <<4*numProcs<<std::endl;
00457       return(-6);
00458     }
00459   }
00460 
00461 // now let's do the checks for Acopy...
00462 
00463   if (map.MyGID(0)) {
00464     EPETRA_CHK_ERR( Acopy.ExtractGlobalRowCopy(0, len, numIndices,
00465                                             values, indices) );
00466     if (numIndices != 4) {
00467       return(-1);
00468     }
00469     if (indices[0] != 0) {
00470       return(-2);
00471     }
00472 
00473     if (values[0] != 1.0*numProcs) {
00474       std::cout << "ERROR: Acopy.values[0] ("<<values[0]<<") should be "<<numProcs<<std::endl;
00475       return(-3);
00476     }
00477   }
00478 
00479   if (map.MyGID(4)) {
00480     EPETRA_CHK_ERR( Acopy.ExtractGlobalRowCopy(4, len, numIndices,
00481                                             values, indices) );
00482 
00483     if (numIndices != 9) {
00484       return(-4);
00485     }
00486     int lcid = A->LCID(4);
00487     if (lcid<0) {
00488       return(-5);
00489     }
00490     if (values[lcid] != 4.0*numProcs) {
00491       std::cout << "ERROR: Acopy.values["<<lcid<<"] ("<<values[lcid]<<") should be "
00492            <<4*numProcs<<std::endl;
00493       return(-6);
00494     }
00495   }
00496 
00497 // now let's do the checks for Acopy2...
00498 
00499   if (map.MyGID(0)) {
00500     EPETRA_CHK_ERR( Acopy2.ExtractGlobalRowCopy(0, len, numIndices,
00501                                             values, indices) );
00502     if (numIndices != 4) {
00503       return(-1);
00504     }
00505     if (indices[0] != 0) {
00506       return(-2);
00507     }
00508 
00509     if (values[0] != 1.0*numProcs) {
00510       std::cout << "ERROR: Acopy2.values[0] ("<<values[0]<<") should be "<<numProcs<<std::endl;
00511       return(-3);
00512     }
00513   }
00514 
00515   if (map.MyGID(4)) {
00516     EPETRA_CHK_ERR( Acopy2.ExtractGlobalRowCopy(4, len, numIndices,
00517                                             values, indices) );
00518 
00519     if (numIndices != 9) {
00520       return(-4);
00521     }
00522     int lcid = A->LCID(4);
00523     if (lcid<0) {
00524       return(-5);
00525     }
00526     if (values[lcid] != 4.0*numProcs) {
00527       std::cout << "ERROR: Acopy2.values["<<lcid<<"] ("<<values[lcid]<<") should be "
00528            <<4*numProcs<<std::endl;
00529       return(-6);
00530     }
00531   }
00532 
00533   delete [] values_2d;
00534   delete [] values_1d;
00535   delete [] nodes;
00536   delete [] indices;
00537   delete [] values;
00538 
00539   delete A;
00540   delete graph;
00541 
00542   return(0);
00543 }
00544 
00545 int Young1(const Epetra_Comm& Comm, bool verbose)
00546 {
00547   //This is a test case submitted by Joe Young with bug 2421. It runs
00548   //only on 2 processors.
00549   if (Comm.NumProc() != 2) {
00550     return(0);
00551   }
00552 
00553   // Give rows 0-2 to proc 0 and 3-5 to proc 1
00554   int             RowIndices[3];
00555   if (Comm.MyPID() == 0) {
00556     RowIndices[0] = 0;
00557     RowIndices[1] = 1;
00558     RowIndices[2] = 2;
00559   } else {
00560     RowIndices[0] = 3;
00561     RowIndices[1] = 4;
00562     RowIndices[2] = 5;
00563   }
00564   Epetra_Map      RangeMap(-1, 3, RowIndices, 0, Comm);
00565   Epetra_Map & RowMap = RangeMap;
00566 
00567   // Define a second map that gives col 0 to proc 0 and col 1 to proc 1 
00568   int             ColIndices[1];
00569   if (Comm.MyPID() == 0) {
00570     ColIndices[0] = 0;
00571   }
00572   else {
00573     ColIndices[0] = 1;
00574   }
00575   Epetra_Map      DomainMap(-1, 1, ColIndices, 0, Comm);
00576 
00577   // Construct a graph where both processors only insert into local
00578   // elements
00579   Epetra_FECrsGraph BrokenGraph(Copy, RowMap, 2);
00580   for (int i = 0; i < RangeMap.NumMyElements(); i++) {
00581     int             ig = RowIndices[i];
00582     int             jgs[2] = { 0, 1 };
00583     BrokenGraph.InsertGlobalIndices(1, &ig, 2, jgs);
00584   }
00585   BrokenGraph.GlobalAssemble(DomainMap, RangeMap);
00586 
00587   // Check the size of the matrix that would be created from the graph 
00588   int numCols1 = BrokenGraph.NumGlobalCols();
00589   if (verbose) {
00590     std::cout << "Number of global rows in the graph where only "
00591         "local elements were inserted: " << BrokenGraph.NumGlobalRows()
00592       << std::endl;
00593     std::cout << "Number of global cols in the graph where only "
00594         "local elements were inserted: " << BrokenGraph.NumGlobalCols()
00595       << std::endl;
00596   }
00597   // Construct a graph where both processors insert into global elements
00598   Epetra_FECrsGraph Graph(Copy, RowMap, 2);
00599   for (int i = 0; i < 6; i++) {
00600     int             ig = i;
00601     int             jgs[2] = { 0, 1 };
00602     Graph.InsertGlobalIndices(1, &ig, 2, jgs);
00603   }
00604   Graph.GlobalAssemble(DomainMap, RangeMap);
00605 
00606   // Check the size of the matrix that would be created from the graph 
00607   int numCols2 = Graph.NumGlobalCols();
00608   if (verbose) {
00609     std::cout << "Number of global rows in the graph where "
00610         "global elements were inserted: " << Graph.NumGlobalRows()
00611        << std::endl;
00612     std::cout << "Number of global cols in the graph where "
00613         "global elements were inserted: " << Graph.NumGlobalCols()
00614       << std::endl;
00615   }
00616 
00617   if (numCols1 != numCols2) return(-1);
00618   return(0);
00619 }
00620 
00621 int rectangular(const Epetra_Comm& Comm, bool verbose)
00622 {
00623   int mypid = Comm.MyPID();
00624   int numlocalrows = 3;
00625   Epetra_Map rowmap(-1, numlocalrows, 0, Comm);
00626 
00627   int numglobalrows = numlocalrows*Comm.NumProc();
00628 
00629   int numcols = 2*numglobalrows;
00630 
00631   Epetra_FECrsGraph fegraph(Copy, rowmap, numcols);
00632 
00633   int* cols = new int[numcols];
00634   for(int j=0; j<numcols; ++j) cols[j] = j;
00635 
00636   Epetra_Map domainmap(-1, numcols, 0, Comm);
00637 
00638   int firstlocalrow = numlocalrows*mypid;
00639   int lastlocalrow = numlocalrows*(mypid+1)-1;
00640 
00641   for(int i=0; i<numglobalrows; ++i) {
00642     //if i is a local row, then skip it. We want each processor to only
00643     //load rows that belong on other processors.
00644     if (i >= firstlocalrow && i <= lastlocalrow) continue;
00645 
00646     EPETRA_CHK_ERR( fegraph.InsertGlobalIndices(1, &i, numcols, &(cols[0])) );
00647   }
00648 
00649   EPETRA_CHK_ERR( fegraph.GlobalAssemble(domainmap, rowmap) );
00650 
00651   if (verbose) {
00652     std::cout << "********************** fegraph **********************" << std::endl;
00653     std::cout << fegraph << std::endl;
00654   }
00655 
00656   delete [] cols;
00657 
00658   return(0);
00659 }
00660 

Generated on Wed May 12 21:41:05 2010 for Epetra Package Browser (Single Doxygen Collection) by  doxygen 1.4.7