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

Generated on Thu Sep 18 12:37:58 2008 for Epetra Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1