Epetra Package Browser (Single Doxygen Collection) Development
FECrsMatrix/ExecuteTestProblems.cpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ************************************************************************
00003 //
00004 //               Epetra: Linear Algebra Services Package
00005 //                 Copyright 2011 Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ************************************************************************
00040 //@HEADER
00041 
00042 
00043 #include "Epetra_BLAS.h"
00044 #include "ExecuteTestProblems.h"
00045 #include "Epetra_Comm.h"
00046 #include "Epetra_Vector.h"
00047 #include "Epetra_FEVector.h"
00048 #include "Epetra_IntSerialDenseVector.h"
00049 #include "Epetra_SerialDenseMatrix.h"
00050 #include <../src/Epetra_matrix_data.h>
00051 #include <../src/Epetra_test_functions.h>
00052 
00053 
00054 int Drumm1(const Epetra_Map& map, bool verbose)
00055 {
00056   (void)verbose;
00057   //Simple 2-element problem (element as in "finite-element") from
00058   //Clif Drumm. Two triangular elements, one per processor, as shown
00059   //here:
00060   //
00061   //   *----*
00062   //  3|\  2|
00063   //   | \  |
00064   //   | 0\1|
00065   //   |   \|
00066   //   *----*
00067   //  0    1
00068   //
00069   //Element 0 on processor 0, element 1 on processor 1.
00070   //Processor 0 will own nodes 0,1 and processor 1 will own nodes 2,3.
00071   //Each processor will pass a 3x3 element-matrix to Epetra_FECrsMatrix.
00072   //After GlobalAssemble(), the matrix should be as follows:
00073   //
00074   //         row 0: 2  1  0  1
00075   //proc 0   row 1: 1  4  1  2
00076   //----------------------------------
00077   //         row 2: 0  1  2  1
00078   //proc 1   row 3: 1  2  1  4
00079   //
00080 
00081   int numProcs = map.Comm().NumProc();
00082   int localProc = map.Comm().MyPID();
00083 
00084   if (numProcs != 2) return(0);
00085 
00086   //so first we'll set up a epetra_test::matrix_data object with
00087   //contents that match the above-described matrix. (but the
00088   //matrix_data object will have all 4 rows on each processor)
00089 
00090   int i;
00091   int rowlengths[4];
00092   rowlengths[0] = 3;
00093   rowlengths[1] = 4;
00094   rowlengths[2] = 3;
00095   rowlengths[3] = 4;
00096 
00097   epetra_test::matrix_data matdata(4, rowlengths);
00098   for(i=0; i<4; ++i) {
00099     for(int j=0; j<matdata.rowlengths()[i]; ++j) {
00100       matdata.colindices()[i][j] = j;
00101     }
00102   }
00103 
00104   matdata.colindices()[0][2] = 3;
00105 
00106   matdata.colindices()[2][0] = 1;
00107   matdata.colindices()[2][1] = 2;
00108   matdata.colindices()[2][2] = 3;
00109 
00110   double** coefs = matdata.coefs();
00111   coefs[0][0] = 2.0; coefs[0][1] = 1.0;                    coefs[0][2] = 1.0;
00112   coefs[1][0] = 1.0; coefs[1][1] = 4.0; coefs[1][2] = 1.0; coefs[1][3] = 2.0;
00113                      coefs[2][0] = 1.0; coefs[2][1] = 2.0; coefs[2][2] = 1.0;
00114   coefs[3][0] = 1.0; coefs[3][1] = 2.0; coefs[3][2] = 1.0; coefs[3][3] = 4.0;
00115 
00116   //now we'll load a Epetra_FECrsMatrix with data that matches the
00117   //above-described finite-element problem.
00118 
00119   int indexBase = 0, ierr = 0;
00120   int myNodes[4];
00121   double values[9];
00122   values[0] = 2.0;
00123   values[1] = 1.0;
00124   values[2] = 1.0;
00125   values[3] = 1.0;
00126   values[4] = 2.0;
00127   values[5] = 1.0;
00128   values[6] = 1.0;
00129   values[7] = 1.0;
00130   values[8] = 2.0;
00131 
00132   int numMyNodes = 2;
00133 
00134   if (localProc == 0) {
00135     myNodes[0] = 0;
00136     myNodes[1] = 1;
00137   }
00138   else {
00139     myNodes[0] = 2;
00140     myNodes[1] = 3;
00141   }
00142 
00143   Epetra_Map Map(-1, numMyNodes, myNodes, indexBase, map.Comm());
00144 
00145   numMyNodes = 3;
00146 
00147   if (localProc == 0) {
00148     myNodes[0] = 0;
00149     myNodes[1] = 1;
00150     myNodes[2] = 3;
00151   }
00152   else {
00153     myNodes[0] = 1;
00154     myNodes[1] = 2;
00155     myNodes[2] = 3;
00156   }
00157 
00158   int rowLengths = 3;
00159   Epetra_FECrsMatrix A(Copy, Map, rowLengths);
00160 
00161   EPETRA_TEST_ERR( A.InsertGlobalValues(numMyNodes, myNodes,
00162                                         numMyNodes, myNodes, values,
00163                                         Epetra_FECrsMatrix::ROW_MAJOR),ierr);
00164 
00165   EPETRA_TEST_ERR( A.GlobalAssemble(), ierr );
00166   EPETRA_TEST_ERR( A.GlobalAssemble(), ierr );
00167 
00168   //now the test is to check whether the FECrsMatrix data matches the
00169   //epetra_test::matrix_data object...
00170 
00171   bool the_same = matdata.compare_local_data(A);
00172 
00173   if (!the_same) {
00174     return(-1);
00175   }
00176 
00177   return(0);
00178 }
00179 
00180 int Drumm2(const Epetra_Map& map, bool verbose)
00181 {
00182   //Simple 2-element problem (element as in "finite-element") from
00183   //Clif Drumm. Two triangular elements, one per processor, as shown
00184   //here:
00185   //
00186   //   *----*
00187   //  3|\  2|
00188   //   | \  |
00189   //   | 0\1|
00190   //   |   \|
00191   //   *----*
00192   //  0    1
00193   //
00194   //Element 0 on processor 0, element 1 on processor 1.
00195   //Processor 0 will own nodes 0,1,3 and processor 1 will own node 2.
00196   //Each processor will pass a 3x3 element-matrix to Epetra_FECrsMatrix.
00197   //After GlobalAssemble(), the matrix should be as follows:
00198   //
00199   //         row 0: 2  1  0  1
00200   //proc 0   row 1: 1  4  1  2
00201   //         row 2: 0  1  2  1
00202   //----------------------------------
00203   //proc 1   row 3: 1  2  1  4
00204   //
00205 
00206   int numProcs = map.Comm().NumProc();
00207   int localProc = map.Comm().MyPID();
00208 
00209   if (numProcs != 2) return(0);
00210 
00211   int indexBase = 0, ierr = 0;
00212 
00213   double* values = new double[9];
00214   values[0] = 2.0;
00215   values[1] = 1.0;
00216   values[2] = 1.0;
00217   values[3] = 1.0;
00218   values[4] = 2.0;
00219   values[5] = 1.0;
00220   values[6] = 1.0;
00221   values[7] = 1.0;
00222   values[8] = 2.0;
00223 
00224   if (localProc == 0) {
00225     int numMyNodes = 3;
00226     int* myNodes = new int[numMyNodes];
00227     myNodes[0] = 0;
00228     myNodes[1] = 1;
00229     myNodes[2] = 3;
00230 
00231     Epetra_Map Map(-1, numMyNodes, myNodes, indexBase, map.Comm());
00232 
00233     int rowLengths = 3;
00234     Epetra_FECrsMatrix A(Copy, Map, rowLengths);
00235 
00236     EPETRA_TEST_ERR( A.InsertGlobalValues(numMyNodes, myNodes,
00237         numMyNodes, myNodes,
00238         values, Epetra_FECrsMatrix::ROW_MAJOR),ierr);
00239 
00240     EPETRA_TEST_ERR( A.GlobalAssemble(), ierr );
00241     EPETRA_TEST_ERR( A.GlobalAssemble(), ierr );
00242 
00243     if (verbose) {
00244     A.Print(cout);
00245     }
00246 
00247     //now let's make sure we can do a matvec with this matrix.
00248     Epetra_Vector x(Map), y(Map);
00249     x.PutScalar(1.0);
00250     EPETRA_TEST_ERR( A.Multiply(false, x, y), ierr);
00251 
00252     if (verbose&&localProc==0) {
00253       cout << "y = A*x, x=1.0's"<<endl;
00254     }
00255 
00256     if (verbose) {
00257     y.Print(cout);
00258     }
00259 
00260     delete [] myNodes;
00261     delete [] values;
00262   }
00263   else {
00264     int numMyNodes = 1;
00265     int* myNodes = new int[numMyNodes];
00266     myNodes[0] = 2;
00267 
00268     Epetra_Map Map(-1, numMyNodes, myNodes, indexBase, map.Comm());
00269 
00270     int rowLengths = 3;
00271     Epetra_FECrsMatrix A(Copy, Map, rowLengths);
00272 
00273     delete [] myNodes;
00274     numMyNodes = 3;
00275     myNodes = new int[numMyNodes];
00276     myNodes[0] = 1;
00277     myNodes[1] = 2;
00278     myNodes[2] = 3;
00279 
00280     EPETRA_TEST_ERR( A.InsertGlobalValues(numMyNodes, myNodes,
00281         numMyNodes, myNodes,
00282         values, Epetra_FECrsMatrix::ROW_MAJOR),ierr);
00283 
00284     EPETRA_TEST_ERR( A.GlobalAssemble(), ierr );
00285     EPETRA_TEST_ERR( A.GlobalAssemble(), ierr );
00286 
00287     if (verbose) {
00288     A.Print(cout);
00289     }
00290 
00291     //now let's make sure we can do a matvec with this matrix.
00292     Epetra_Vector x(Map), y(Map);
00293     x.PutScalar(1.0);
00294     EPETRA_TEST_ERR( A.Multiply(false, x, y), ierr);
00295 
00296     if (verbose) {
00297     y.Print(cout);
00298     }
00299 
00300     delete [] myNodes;
00301     delete [] values;
00302   }
00303 
00304   return(0);
00305 }
00306 
00307 int Drumm3(const Epetra_Map& map, bool verbose)
00308 {
00309   const Epetra_Comm & Comm = map.Comm();
00310 
00311   /* get number of processors and the name of this processor */
00312 
00313   int Numprocs = Comm.NumProc();
00314   int MyPID   = Comm.MyPID();
00315 
00316   if (Numprocs != 2) return(0);
00317 
00318   int NumGlobalRows = 4;
00319   int IndexBase = 0;
00320   Epetra_Map Map(NumGlobalRows, IndexBase, Comm);
00321 
00322   // Construct FECrsMatrix
00323 
00324   int NumEntriesPerRow = 3;
00325 
00326   Epetra_FECrsMatrix A(Copy, Map, NumEntriesPerRow);
00327 
00328   double ElementArea = 0.5;
00329 
00330   int NumCols = 3;
00331   int* Indices = new int[NumCols];
00332 
00333   if(MyPID==0)  // indices corresponding to element 0 on processor 0
00334   {
00335     Indices[0] = 0;
00336     Indices[1] = 1;
00337     Indices[2] = 3;
00338   }
00339   else if(MyPID==1)  // indices corresponding to element 1 on processor 1
00340   {
00341     Indices[0] = 1;
00342     Indices[1] = 2;
00343     Indices[2] = 3;
00344   }
00345 
00346   double* Values = new double[NumCols*NumCols];
00347 
00348 // removal term
00349   Values[0] = 2*ElementArea/12.;
00350   Values[1] = 1*ElementArea/12.;
00351   Values[2] = 1*ElementArea/12.;
00352   Values[3] = 1*ElementArea/12.;
00353   Values[4] = 2*ElementArea/12.;
00354   Values[5] = 1*ElementArea/12.;
00355   Values[6] = 1*ElementArea/12.;
00356   Values[7] = 1*ElementArea/12.;
00357   Values[8] = 2*ElementArea/12.;
00358 
00359   A.InsertGlobalValues(NumCols, Indices,
00360                         Values,
00361                         Epetra_FECrsMatrix::ROW_MAJOR);
00362 
00363   A.GlobalAssemble();
00364   A.GlobalAssemble();
00365 
00366 //  A.Print(cout);
00367 
00368 // Create vectors for CG algorithm
00369 
00370   Epetra_FEVector* bptr = new Epetra_FEVector(A.RowMap(), 1);
00371   Epetra_FEVector* x0ptr = new Epetra_FEVector(A.RowMap(), 1);
00372 
00373   Epetra_FEVector& b = *bptr;
00374   Epetra_FEVector& x0 = *x0ptr;
00375 
00376   // source terms
00377   NumCols = 2;
00378 
00379   if(MyPID==0)  // indices corresponding to element 0 on processor 0
00380   {
00381     Indices[0] = 0;
00382     Indices[1] = 3;
00383 
00384     Values[0] = 1./2.;
00385     Values[1] = 1./2.;
00386 
00387    }
00388    else
00389    {
00390     Indices[0] = 1;
00391     Indices[1] = 2;
00392 
00393     Values[0] = 0;
00394     Values[1] = 0;
00395    }
00396 
00397   b.SumIntoGlobalValues(NumCols, Indices, Values);
00398 
00399   b.GlobalAssemble();
00400 
00401   if (verbose&&MyPID==0) cout << "b:" << endl;
00402   if (verbose) {
00403   b.Print(cout);
00404   }
00405 
00406   x0 = b;
00407 
00408   if (verbose&&MyPID==0) {
00409   cout << "x:"<<endl;
00410   }
00411 
00412   if (verbose) {
00413   x0.Print(cout);
00414   }
00415 
00416   delete [] Values;
00417   delete [] Indices;
00418 
00419   delete bptr;
00420   delete x0ptr;
00421 
00422   return(0);
00423 }
00424 
00425 int four_quads(const Epetra_Comm& Comm, bool preconstruct_graph, bool verbose)
00426 {
00427   if (verbose) {
00428     cout << "******************* four_quads ***********************"<<endl;
00429   }
00430 
00431   //This function assembles a matrix representing a finite-element mesh
00432   //of four 2-D quad elements. There are 9 nodes in the problem. The
00433   //same problem is assembled no matter how many processors are being used
00434   //(within reason). It may not work if more than 9 processors are used.
00435   //
00436   //  *------*------*
00437   // 6|     7|     8|
00438   //  | E2   | E3   |
00439   //  *------*------*
00440   // 3|     4|     5|
00441   //  | E0   | E1   |
00442   //  *------*------*
00443   // 0      1      2
00444   //
00445   //Nodes are denoted by * with node-numbers below and left of each node.
00446   //E0, E1 and so on are element-numbers.
00447   //
00448   //Each processor will contribute a sub-matrix of size 4x4, filled with 1's,
00449   //for each element. Thus, the coefficient value at position 0,0 should end up
00450   //being 1.0*numProcs, the value at position 4,4 should be 1.0*4*numProcs, etc.
00451   //
00452   //Depending on the number of processors being used, the locations of the
00453   //specific matrix positions (in terms of which processor owns them) will vary.
00454   //
00455 
00456   int numProcs = Comm.NumProc();
00457 
00458   int numNodes = 9;
00459   int numElems = 4;
00460   int numNodesPerElem = 4;
00461 
00462   int indexBase = 0;
00463 
00464   //Create a map using epetra-defined linear distribution.
00465   Epetra_Map map(numNodes, indexBase, Comm);
00466 
00467   Epetra_CrsGraph* graph = NULL;
00468 
00469   int* nodes = new int[numNodesPerElem];
00470   int i, j, err = 0;
00471 
00472   if (preconstruct_graph) {
00473     graph = new Epetra_CrsGraph(Copy, map, 1);
00474 
00475     //we're going to fill the graph with indices, but remember it will only
00476     //accept indices in rows for which map.MyGID(row) is true.
00477 
00478     for(i=0; i<numElems; ++i) {
00479       switch(i) {
00480       case 0:
00481   nodes[0] = 0; nodes[1] = 1; nodes[2] = 4; nodes[3] = 3;
00482   break;
00483       case 1:
00484   nodes[0] = 1; nodes[1] = 2; nodes[2] = 5; nodes[3] = 4;
00485   break;
00486       case 2:
00487   nodes[0] = 3; nodes[1] = 4; nodes[2] = 7; nodes[3] = 6;
00488   break;
00489       case 3:
00490   nodes[0] = 4; nodes[1] = 5; nodes[2] = 8; nodes[3] = 7;
00491   break;
00492       }
00493 
00494       for(j=0; j<numNodesPerElem; ++j) {
00495   if (map.MyGID(nodes[j])) {
00496     err = graph->InsertGlobalIndices(nodes[j], numNodesPerElem,
00497              nodes);
00498     if (err<0) return(err);
00499   }
00500       }
00501     }
00502 
00503     EPETRA_CHK_ERR( graph->FillComplete() );
00504   }
00505 
00506   Epetra_FECrsMatrix* A = NULL;
00507 
00508   if (preconstruct_graph) {
00509     A = new Epetra_FECrsMatrix(Copy, *graph);
00510   }
00511   else {
00512     A = new Epetra_FECrsMatrix(Copy, map, 1);
00513   }
00514 
00515   EPETRA_CHK_ERR( A->PutScalar(0.0) );
00516 
00517   double* values_1d = new double[numNodesPerElem*numNodesPerElem];
00518   double** values_2d = new double*[numNodesPerElem];
00519 
00520   for(i=0; i<numNodesPerElem*numNodesPerElem; ++i) values_1d[i] = 1.0;
00521 
00522   int offset = 0;
00523   for(i=0; i<numNodesPerElem; ++i) {
00524     values_2d[i] = &(values_1d[offset]);
00525     offset += numNodesPerElem;
00526   }
00527 
00528   int format = Epetra_FECrsMatrix::ROW_MAJOR;
00529   Epetra_IntSerialDenseVector epetra_nodes(View, nodes, numNodesPerElem);
00530   Epetra_SerialDenseMatrix epetra_values(View, values_1d, numNodesPerElem,
00531            numNodesPerElem, numNodesPerElem);
00532 
00533   for(i=0; i<numElems; ++i) {
00534     switch(i) {
00535     case 0:
00536       nodes[0] = 0; nodes[1] = 1; nodes[2] = 4; nodes[3] = 3;
00537       if (preconstruct_graph) {
00538   err = A->SumIntoGlobalValues(epetra_nodes,
00539              epetra_values, format);
00540   if (err<0) return(err);
00541       }
00542       else {
00543   err = A->InsertGlobalValues(epetra_nodes,
00544             epetra_values, format);
00545   if (err<0) return(err);
00546       }
00547       break;
00548 
00549     case 1:
00550       nodes[0] = 1; nodes[1] = 2; nodes[2] = 5; nodes[3] = 4;
00551       if (preconstruct_graph) {
00552   err = A->SumIntoGlobalValues(nodes[0], numNodesPerElem, values_2d[0],
00553                                      nodes);
00554   err += A->SumIntoGlobalValues(nodes[1], numNodesPerElem, values_2d[1],
00555                                      nodes);
00556   err += A->SumIntoGlobalValues(nodes[2], numNodesPerElem, values_2d[2],
00557                                      nodes);
00558   err += A->SumIntoGlobalValues(nodes[3], numNodesPerElem, values_2d[3],
00559                                      nodes);
00560   if (err<0) return(err);
00561       }
00562       else {
00563   err = A->InsertGlobalValues(numNodesPerElem, nodes,
00564             values_2d, format);
00565   if (err<0) return(err);
00566       }
00567       break;
00568 
00569     case 2:
00570       nodes[0] = 3; nodes[1] = 4; nodes[2] = 7; nodes[3] = 6;
00571       if (preconstruct_graph) {
00572   err = A->SumIntoGlobalValues(numNodesPerElem, nodes,
00573              numNodesPerElem, nodes,
00574              values_1d, format);
00575   if (err<0) return(err);
00576       }
00577       else {
00578   err = A->InsertGlobalValues(numNodesPerElem, nodes,
00579             numNodesPerElem, nodes,
00580             values_1d, format);
00581   if (err<0) return(err);
00582       }
00583       break;
00584 
00585      case 3:
00586       nodes[0] = 4; nodes[1] = 5; nodes[2] = 8; nodes[3] = 7;
00587       if (preconstruct_graph) {
00588   err = A->SumIntoGlobalValues(numNodesPerElem, nodes,
00589              numNodesPerElem, nodes,
00590              values_2d, format);
00591   if (err<0) return(err);
00592       }
00593       else {
00594   err = A->InsertGlobalValues(numNodesPerElem, nodes,
00595             numNodesPerElem, nodes,
00596             values_2d, format);
00597   if (err<0) return(err);
00598       }
00599       break;
00600     }
00601   }
00602 
00603   err = A->GlobalAssemble();
00604   if (err < 0) {
00605     return(err);
00606   }
00607 
00608   Epetra_Vector x(A->RowMap()), y(A->RowMap());
00609 
00610   x.PutScalar(1.0); y.PutScalar(0.0);
00611 
00612   Epetra_FECrsMatrix Acopy(*A);
00613 
00614   err = Acopy.GlobalAssemble();
00615   if (err < 0) {
00616     return(err);
00617   }
00618 
00619   bool the_same = epetra_test::compare_matrices(*A, Acopy);
00620   if (!the_same) {
00621     return(-1);
00622   }
00623 
00624   Epetra_FECrsMatrix Acopy2(Copy, A->RowMap(), A->ColMap(), 1);
00625 
00626   Acopy2 = Acopy;
00627 
00628   the_same = epetra_test::compare_matrices(*A, Acopy);
00629   if (!the_same) {
00630     return(-1);
00631   }
00632 
00633   int len = 20;
00634   int* indices = new int[len];
00635   double* values = new double[len];
00636   int numIndices;
00637 
00638   if (map.MyGID(0)) {
00639     EPETRA_CHK_ERR( A->ExtractGlobalRowCopy(0, len, numIndices,
00640               values, indices) );
00641     if (numIndices != 4) {
00642       return(-1);
00643     }
00644     if (indices[0] != 0) {
00645       return(-2);
00646     }
00647 
00648     if (values[0] != 1.0*numProcs) {
00649       cout << "ERROR: values[0] ("<<values[0]<<") should be "<<numProcs<<endl;
00650       return(-3);
00651     }
00652   }
00653 
00654   if (map.MyGID(4)) {
00655     EPETRA_CHK_ERR( A->ExtractGlobalRowCopy(4, len, numIndices,
00656               values, indices) );
00657 
00658     if (numIndices != 9) {
00659       return(-4);
00660     }
00661     int lcid = A->LCID(4);
00662     if (lcid<0) {
00663       return(-5);
00664     }
00665     if (values[lcid] != 4.0*numProcs) {
00666       cout << "ERROR: values["<<lcid<<"] ("<<values[lcid]<<") should be "
00667      <<4*numProcs<<endl;
00668       return(-6);
00669     }
00670   }
00671 
00672   delete [] values_2d;
00673   delete [] values_1d;
00674   delete [] nodes;
00675   delete [] indices;
00676   delete [] values;
00677 
00678   delete A;
00679   delete graph;
00680 
00681   return(0);
00682 }
00683 
00684 int submatrix_formats(const Epetra_Comm& Comm, bool verbose)
00685 {
00686   (void)verbose;
00687   //
00688   //This function simply verifies that the ROW_MAJOR/COLUMN_MAJOR switch works.
00689   //
00690   int numProcs = Comm.NumProc();
00691   int myPID = Comm.MyPID();
00692 
00693   int numLocalElements = 3;
00694   int numGlobalElements = numLocalElements*numProcs;
00695   int indexBase = 0;
00696 
00697   Epetra_Map map(numGlobalElements, numLocalElements, indexBase, Comm);
00698 
00699   Epetra_FECrsMatrix A(Copy, map, numLocalElements);
00700 
00701   Epetra_IntSerialDenseVector epetra_indices(numLocalElements);
00702 
00703   int firstGlobalElement = numLocalElements*myPID;
00704 
00705   int i, j;
00706   for(i=0; i<numLocalElements; ++i) {
00707     epetra_indices[i] = firstGlobalElement+i;
00708   }
00709 
00710   Epetra_SerialDenseMatrix submatrix(numLocalElements, numLocalElements);
00711 
00712   for(i=0; i<numLocalElements; ++i) {
00713     for(j=0; j<numLocalElements; ++j) {
00714       submatrix(i,j) = 1.0*(firstGlobalElement+i);
00715     }
00716   }
00717 
00718   EPETRA_CHK_ERR( A.InsertGlobalValues(epetra_indices, submatrix,
00719                                Epetra_FECrsMatrix::COLUMN_MAJOR) );
00720 
00721   EPETRA_CHK_ERR( A.GlobalAssemble() );
00722 
00723   int len = 20;
00724   int numIndices;
00725   int* indices = new int[len];
00726   double* coefs = new double[len];
00727 
00728   for(i=0; i<numLocalElements; ++i) {
00729     int row = firstGlobalElement+i;
00730 
00731     EPETRA_CHK_ERR( A.ExtractGlobalRowCopy(row, len, numIndices,
00732              coefs, indices) );
00733 
00734     for(j=0; j<numIndices; ++j) {
00735       if (coefs[j] != 1.0*row) {
00736   return(-2);
00737       }
00738     }
00739   }
00740 
00741   //now reset submatrix (transposing the i,j indices)
00742 
00743   for(i=0; i<numLocalElements; ++i) {
00744     for(j=0; j<numLocalElements; ++j) {
00745       submatrix(j,i) = 1.0*(firstGlobalElement+i);
00746     }
00747   }
00748 
00749   //sum these values into the matrix using the ROW_MAJOR switch, which should
00750   //result in doubling what's already there from the above Insert operation.
00751 
00752   EPETRA_CHK_ERR( A.SumIntoGlobalValues(epetra_indices, submatrix,
00753           Epetra_FECrsMatrix::ROW_MAJOR) );
00754 
00755   EPETRA_CHK_ERR( A.GlobalAssemble() );
00756 
00757   for(i=0; i<numLocalElements; ++i) {
00758     int row = firstGlobalElement+i;
00759 
00760     EPETRA_CHK_ERR( A.ExtractGlobalRowCopy(row, len, numIndices,
00761              coefs, indices) );
00762 
00763     for(j=0; j<numIndices; ++j) {
00764       if (coefs[j] != 2.0*row) {
00765   return(-3);
00766       }
00767     }
00768   }
00769 
00770   delete [] indices;
00771   delete [] coefs;
00772 
00773   return(0);
00774 }
00775 
00776 int rectangular(const Epetra_Comm& Comm, bool verbose)
00777 {
00778   (void)verbose;
00779   int numprocs = Comm.NumProc();
00780   int localproc = Comm.MyPID();
00781   int numMyRows = 2;
00782   int numGlobalRows = numprocs*numMyRows;
00783   int* myrows = new int[numMyRows];
00784 
00785   int myFirstRow = localproc*numMyRows;
00786   int i;
00787   for(i=0; i<numMyRows; ++i) {
00788     myrows[i] = myFirstRow+i;
00789   }
00790 
00791   Epetra_Map map(numGlobalRows, numMyRows, myrows, 0, Comm);
00792 
00793   Epetra_FECrsMatrix A(Copy, map, 30);
00794 
00795   int numcols = 20;
00796   int* cols = new int[numcols];
00797   for(i=0; i<numcols; ++i) {
00798     cols[i] = i;
00799   }
00800 
00801   double* coefs = new double[numGlobalRows*numcols];
00802   int offset = 0;
00803   for(int j=0; j<numcols; ++j) {
00804     for(i=0; i<numGlobalRows; ++i) {
00805       coefs[offset++] = 1.0*i;
00806     }
00807   }
00808 
00809   int* globalRows = new int[numGlobalRows];
00810   for(i=0; i<numGlobalRows; ++i) globalRows[i] = i;
00811 
00812   EPETRA_CHK_ERR( A.InsertGlobalValues(numGlobalRows, globalRows,
00813                                        numcols, cols, coefs,
00814                                      Epetra_FECrsMatrix::COLUMN_MAJOR));
00815   delete [] coefs;
00816   delete [] globalRows;
00817 
00818   //Since the matrix is rectangular, we need to call GlobalAssemble with
00819   //a domain-map and a range-map. Otherwise, GlobalAssemble has no way of
00820   //knowing what the domain-map and range-map should be.
00821   //We'll use a linear distribution of the columns for a domain-map, and
00822   //our original row-map for the range-map.
00823   int numMyCols = numcols/numprocs;
00824   int rem = numcols%numprocs;
00825   if (localproc<rem) ++numMyCols;
00826   Epetra_Map domainmap(numcols, numMyCols, 0, Comm);
00827 
00828   EPETRA_CHK_ERR( A.GlobalAssemble(domainmap, map) );
00829 
00830   int numGlobalCols = A.NumGlobalCols();
00831   int numGlobalNNZ = A.NumGlobalNonzeros();
00832 
00833   if (numGlobalCols != numcols ||
00834       numGlobalNNZ != numGlobalRows*numcols) {
00835     return(-1);
00836   }
00837 
00838   delete [] cols;
00839   delete [] myrows;
00840   return(0);
00841 }
00842 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines