Epetra Package Browser (Single Doxygen Collection) Development
FECrsMatrix_LL/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   long long indexBase = 0;
00120   int ierr = 0;
00121   long long myNodes[4];
00122   double values[9];
00123   values[0] = 2.0;
00124   values[1] = 1.0;
00125   values[2] = 1.0;
00126   values[3] = 1.0;
00127   values[4] = 2.0;
00128   values[5] = 1.0;
00129   values[6] = 1.0;
00130   values[7] = 1.0;
00131   values[8] = 2.0;
00132 
00133   int numMyNodes = 2;
00134 
00135   if (localProc == 0) {
00136     myNodes[0] = 0;
00137     myNodes[1] = 1;
00138   }
00139   else {
00140     myNodes[0] = 2;
00141     myNodes[1] = 3;
00142   }
00143 
00144   Epetra_Map Map((long long)-1, numMyNodes, myNodes, indexBase, map.Comm());
00145 
00146   numMyNodes = 3;
00147 
00148   if (localProc == 0) {
00149     myNodes[0] = 0;
00150     myNodes[1] = 1;
00151     myNodes[2] = 3;
00152   }
00153   else {
00154     myNodes[0] = 1;
00155     myNodes[1] = 2;
00156     myNodes[2] = 3;
00157   }
00158 
00159   int rowLengths = 3;
00160   Epetra_FECrsMatrix A(Copy, Map, rowLengths);
00161 
00162   EPETRA_TEST_ERR( A.InsertGlobalValues(numMyNodes, myNodes,
00163                                         numMyNodes, myNodes, values,
00164                                         Epetra_FECrsMatrix::ROW_MAJOR),ierr);
00165 
00166   EPETRA_TEST_ERR( A.GlobalAssemble(), ierr );
00167   EPETRA_TEST_ERR( A.GlobalAssemble(), ierr );
00168 
00169   //now the test is to check whether the FECrsMatrix data matches the
00170   //epetra_test::matrix_data object...
00171 
00172   bool the_same = matdata.compare_local_data(A);
00173 
00174   if (!the_same) {
00175     return(-1);
00176   }
00177 
00178   return(0);
00179 }
00180 
00181 int Drumm2(const Epetra_Map& map, bool verbose)
00182 {
00183   //Simple 2-element problem (element as in "finite-element") from
00184   //Clif Drumm. Two triangular elements, one per processor, as shown
00185   //here:
00186   //
00187   //   *----*
00188   //  3|\  2|
00189   //   | \  |
00190   //   | 0\1|
00191   //   |   \|
00192   //   *----*
00193   //  0    1
00194   //
00195   //Element 0 on processor 0, element 1 on processor 1.
00196   //Processor 0 will own nodes 0,1,3 and processor 1 will own node 2.
00197   //Each processor will pass a 3x3 element-matrix to Epetra_FECrsMatrix.
00198   //After GlobalAssemble(), the matrix should be as follows:
00199   //
00200   //         row 0: 2  1  0  1
00201   //proc 0   row 1: 1  4  1  2
00202   //         row 2: 0  1  2  1
00203   //----------------------------------
00204   //proc 1   row 3: 1  2  1  4
00205   //
00206 
00207   int numProcs = map.Comm().NumProc();
00208   int localProc = map.Comm().MyPID();
00209 
00210   if (numProcs != 2) return(0);
00211 
00212   long long indexBase = 0;
00213   int ierr = 0;
00214 
00215   double* values = new double[9];
00216   values[0] = 2.0;
00217   values[1] = 1.0;
00218   values[2] = 1.0;
00219   values[3] = 1.0;
00220   values[4] = 2.0;
00221   values[5] = 1.0;
00222   values[6] = 1.0;
00223   values[7] = 1.0;
00224   values[8] = 2.0;
00225 
00226   if (localProc == 0) {
00227     int numMyNodes = 3;
00228     long long* myNodes = new long long[numMyNodes];
00229     myNodes[0] = 0;
00230     myNodes[1] = 1;
00231     myNodes[2] = 3;
00232 
00233     Epetra_Map Map((long long)-1, numMyNodes, myNodes, indexBase, map.Comm());
00234 
00235     int rowLengths = 3;
00236     Epetra_FECrsMatrix A(Copy, Map, rowLengths);
00237 
00238     EPETRA_TEST_ERR( A.InsertGlobalValues(numMyNodes, myNodes,
00239         numMyNodes, myNodes,
00240         values, Epetra_FECrsMatrix::ROW_MAJOR),ierr);
00241 
00242     EPETRA_TEST_ERR( A.GlobalAssemble(), ierr );
00243     EPETRA_TEST_ERR( A.GlobalAssemble(), ierr );
00244 
00245     if (verbose) {
00246     A.Print(cout);
00247     }
00248 
00249     //now let's make sure we can do a matvec with this matrix.
00250     Epetra_Vector x(Map), y(Map);
00251     x.PutScalar(1.0);
00252     EPETRA_TEST_ERR( A.Multiply(false, x, y), ierr);
00253 
00254     if (verbose&&localProc==0) {
00255       cout << "y = A*x, x=1.0's"<<endl;
00256     }
00257 
00258     if (verbose) {
00259     y.Print(cout);
00260     }
00261 
00262     delete [] myNodes;
00263     delete [] values;
00264   }
00265   else {
00266     int numMyNodes = 1;
00267     long long* myNodes = new long long[numMyNodes];
00268     myNodes[0] = 2;
00269 
00270     Epetra_Map Map((long long)-1, numMyNodes, myNodes, indexBase, map.Comm());
00271 
00272     int rowLengths = 3;
00273     Epetra_FECrsMatrix A(Copy, Map, rowLengths);
00274 
00275     delete [] myNodes;
00276     numMyNodes = 3;
00277     myNodes = new long long[numMyNodes];
00278     myNodes[0] = 1;
00279     myNodes[1] = 2;
00280     myNodes[2] = 3;
00281 
00282     EPETRA_TEST_ERR( A.InsertGlobalValues(numMyNodes, myNodes,
00283         numMyNodes, myNodes,
00284         values, Epetra_FECrsMatrix::ROW_MAJOR),ierr);
00285 
00286     EPETRA_TEST_ERR( A.GlobalAssemble(), ierr );
00287     EPETRA_TEST_ERR( A.GlobalAssemble(), ierr );
00288 
00289     if (verbose) {
00290     A.Print(cout);
00291     }
00292 
00293     //now let's make sure we can do a matvec with this matrix.
00294     Epetra_Vector x(Map), y(Map);
00295     x.PutScalar(1.0);
00296     EPETRA_TEST_ERR( A.Multiply(false, x, y), ierr);
00297 
00298     if (verbose) {
00299     y.Print(cout);
00300     }
00301 
00302     delete [] myNodes;
00303     delete [] values;
00304   }
00305 
00306   return(0);
00307 }
00308 
00309 int Drumm3(const Epetra_Map& map, bool verbose)
00310 {
00311   const Epetra_Comm & Comm = map.Comm();
00312 
00313   /* get number of processors and the name of this processor */
00314 
00315   int Numprocs = Comm.NumProc();
00316   int MyPID   = Comm.MyPID();
00317 
00318   if (Numprocs != 2) return(0);
00319 
00320   long long NumGlobalRows = 4;
00321   long long IndexBase = 0;
00322   Epetra_Map Map(NumGlobalRows, IndexBase, Comm);
00323 
00324   // Construct FECrsMatrix
00325 
00326   int NumEntriesPerRow = 3;
00327 
00328   Epetra_FECrsMatrix A(Copy, Map, NumEntriesPerRow);
00329 
00330   double ElementArea = 0.5;
00331 
00332   int NumCols = 3;
00333   long long* Indices = new long long[NumCols];
00334 
00335   if(MyPID==0)  // indices corresponding to element 0 on processor 0
00336   {
00337     Indices[0] = 0;
00338     Indices[1] = 1;
00339     Indices[2] = 3;
00340   }
00341   else if(MyPID==1)  // indices corresponding to element 1 on processor 1
00342   {
00343     Indices[0] = 1;
00344     Indices[1] = 2;
00345     Indices[2] = 3;
00346   }
00347 
00348   double* Values = new double[NumCols*NumCols];
00349 
00350 // removal term
00351   Values[0] = 2*ElementArea/12.;
00352   Values[1] = 1*ElementArea/12.;
00353   Values[2] = 1*ElementArea/12.;
00354   Values[3] = 1*ElementArea/12.;
00355   Values[4] = 2*ElementArea/12.;
00356   Values[5] = 1*ElementArea/12.;
00357   Values[6] = 1*ElementArea/12.;
00358   Values[7] = 1*ElementArea/12.;
00359   Values[8] = 2*ElementArea/12.;
00360 
00361   A.InsertGlobalValues(NumCols, Indices,
00362                         Values,
00363                         Epetra_FECrsMatrix::ROW_MAJOR);
00364 
00365   A.GlobalAssemble();
00366   A.GlobalAssemble();
00367 
00368 //  A.Print(cout);
00369 
00370 // Create vectors for CG algorithm
00371 
00372   Epetra_FEVector* bptr = new Epetra_FEVector(A.RowMap(), 1);
00373   Epetra_FEVector* x0ptr = new Epetra_FEVector(A.RowMap(), 1);
00374 
00375   Epetra_FEVector& b = *bptr;
00376   Epetra_FEVector& x0 = *x0ptr;
00377 
00378   // source terms
00379   NumCols = 2;
00380 
00381   if(MyPID==0)  // indices corresponding to element 0 on processor 0
00382   {
00383     Indices[0] = 0;
00384     Indices[1] = 3;
00385 
00386     Values[0] = 1./2.;
00387     Values[1] = 1./2.;
00388 
00389    }
00390    else
00391    {
00392     Indices[0] = 1;
00393     Indices[1] = 2;
00394 
00395     Values[0] = 0;
00396     Values[1] = 0;
00397    }
00398 
00399   b.SumIntoGlobalValues(NumCols, Indices, Values);
00400 
00401   b.GlobalAssemble();
00402 
00403   if (verbose&&MyPID==0) cout << "b:" << endl;
00404   if (verbose) {
00405   b.Print(cout);
00406   }
00407 
00408   x0 = b;
00409 
00410   if (verbose&&MyPID==0) {
00411   cout << "x:"<<endl;
00412   }
00413 
00414   if (verbose) {
00415   x0.Print(cout);
00416   }
00417 
00418   delete [] Values;
00419   delete [] Indices;
00420 
00421   delete bptr;
00422   delete x0ptr;
00423 
00424   return(0);
00425 }
00426 
00427 int four_quads(const Epetra_Comm& Comm, bool preconstruct_graph, bool verbose)
00428 {
00429   if (verbose) {
00430     cout << "******************* four_quads ***********************"<<endl;
00431   }
00432 
00433   //This function assembles a matrix representing a finite-element mesh
00434   //of four 2-D quad elements. There are 9 nodes in the problem. The
00435   //same problem is assembled no matter how many processors are being used
00436   //(within reason). It may not work if more than 9 processors are used.
00437   //
00438   //  *------*------*
00439   // 6|     7|     8|
00440   //  | E2   | E3   |
00441   //  *------*------*
00442   // 3|     4|     5|
00443   //  | E0   | E1   |
00444   //  *------*------*
00445   // 0      1      2
00446   //
00447   //Nodes are denoted by * with node-numbers below and left of each node.
00448   //E0, E1 and so on are element-numbers.
00449   //
00450   //Each processor will contribute a sub-matrix of size 4x4, filled with 1's,
00451   //for each element. Thus, the coefficient value at position 0,0 should end up
00452   //being 1.0*numProcs, the value at position 4,4 should be 1.0*4*numProcs, etc.
00453   //
00454   //Depending on the number of processors being used, the locations of the
00455   //specific matrix positions (in terms of which processor owns them) will vary.
00456   //
00457 
00458   int numProcs = Comm.NumProc();
00459 
00460   long long numNodes = 9;
00461   int numElems = 4;
00462   int numNodesPerElem = 4;
00463 
00464   long long indexBase = 0;
00465 
00466   //Create a map using epetra-defined linear distribution.
00467   Epetra_Map map(numNodes, indexBase, Comm);
00468 
00469   Epetra_CrsGraph* graph = NULL;
00470 
00471   long long* nodes = new long long[numNodesPerElem];
00472   int i, j, err = 0;
00473 
00474   if (preconstruct_graph) {
00475     graph = new Epetra_CrsGraph(Copy, map, 1);
00476 
00477     //we're going to fill the graph with indices, but remember it will only
00478     //accept indices in rows for which map.MyGID(row) is true.
00479 
00480     for(i=0; i<numElems; ++i) {
00481       switch(i) {
00482       case 0:
00483   nodes[0] = 0; nodes[1] = 1; nodes[2] = 4; nodes[3] = 3;
00484   break;
00485       case 1:
00486   nodes[0] = 1; nodes[1] = 2; nodes[2] = 5; nodes[3] = 4;
00487   break;
00488       case 2:
00489   nodes[0] = 3; nodes[1] = 4; nodes[2] = 7; nodes[3] = 6;
00490   break;
00491       case 3:
00492   nodes[0] = 4; nodes[1] = 5; nodes[2] = 8; nodes[3] = 7;
00493   break;
00494       }
00495 
00496       for(j=0; j<numNodesPerElem; ++j) {
00497   if (map.MyGID(nodes[j])) {
00498     err = graph->InsertGlobalIndices(nodes[j], numNodesPerElem,
00499              nodes);
00500     if (err<0) return(err);
00501   }
00502       }
00503     }
00504 
00505     EPETRA_CHK_ERR( graph->FillComplete() );
00506   }
00507 
00508   Epetra_FECrsMatrix* A = NULL;
00509 
00510   if (preconstruct_graph) {
00511     A = new Epetra_FECrsMatrix(Copy, *graph);
00512   }
00513   else {
00514     A = new Epetra_FECrsMatrix(Copy, map, 1);
00515   }
00516 
00517   EPETRA_CHK_ERR( A->PutScalar(0.0) );
00518 
00519   double* values_1d = new double[numNodesPerElem*numNodesPerElem];
00520   double** values_2d = new double*[numNodesPerElem];
00521 
00522   for(i=0; i<numNodesPerElem*numNodesPerElem; ++i) values_1d[i] = 1.0;
00523 
00524   int offset = 0;
00525   for(i=0; i<numNodesPerElem; ++i) {
00526     values_2d[i] = &(values_1d[offset]);
00527     offset += numNodesPerElem;
00528   }
00529 
00530   int format = Epetra_FECrsMatrix::ROW_MAJOR;
00531   Epetra_LongLongSerialDenseVector epetra_nodes(View, nodes, numNodesPerElem);
00532   Epetra_SerialDenseMatrix epetra_values(View, values_1d, numNodesPerElem,
00533            numNodesPerElem, numNodesPerElem);
00534 
00535   for(i=0; i<numElems; ++i) {
00536     switch(i) {
00537     case 0:
00538       nodes[0] = 0; nodes[1] = 1; nodes[2] = 4; nodes[3] = 3;
00539       if (preconstruct_graph) {
00540   err = A->SumIntoGlobalValues(epetra_nodes,
00541              epetra_values, format);
00542   if (err<0) return(err);
00543       }
00544       else {
00545   err = A->InsertGlobalValues(epetra_nodes,
00546             epetra_values, format);
00547   if (err<0) return(err);
00548       }
00549       break;
00550 
00551     case 1:
00552       nodes[0] = 1; nodes[1] = 2; nodes[2] = 5; nodes[3] = 4;
00553       if (preconstruct_graph) {
00554   err = A->SumIntoGlobalValues(nodes[0], numNodesPerElem, values_2d[0],
00555                                      nodes);
00556   err += A->SumIntoGlobalValues(nodes[1], numNodesPerElem, values_2d[1],
00557                                      nodes);
00558   err += A->SumIntoGlobalValues(nodes[2], numNodesPerElem, values_2d[2],
00559                                      nodes);
00560   err += A->SumIntoGlobalValues(nodes[3], numNodesPerElem, values_2d[3],
00561                                      nodes);
00562   if (err<0) return(err);
00563       }
00564       else {
00565   err = A->InsertGlobalValues(numNodesPerElem, nodes,
00566             values_2d, format);
00567   if (err<0) return(err);
00568       }
00569       break;
00570 
00571     case 2:
00572       nodes[0] = 3; nodes[1] = 4; nodes[2] = 7; nodes[3] = 6;
00573       if (preconstruct_graph) {
00574   err = A->SumIntoGlobalValues(numNodesPerElem, nodes,
00575              numNodesPerElem, nodes,
00576              values_1d, format);
00577   if (err<0) return(err);
00578       }
00579       else {
00580   err = A->InsertGlobalValues(numNodesPerElem, nodes,
00581             numNodesPerElem, nodes,
00582             values_1d, format);
00583   if (err<0) return(err);
00584       }
00585       break;
00586 
00587      case 3:
00588       nodes[0] = 4; nodes[1] = 5; nodes[2] = 8; nodes[3] = 7;
00589       if (preconstruct_graph) {
00590   err = A->SumIntoGlobalValues(numNodesPerElem, nodes,
00591              numNodesPerElem, nodes,
00592              values_2d, format);
00593   if (err<0) return(err);
00594       }
00595       else {
00596   err = A->InsertGlobalValues(numNodesPerElem, nodes,
00597             numNodesPerElem, nodes,
00598             values_2d, format);
00599   if (err<0) return(err);
00600       }
00601       break;
00602     }
00603   }
00604 
00605   err = A->GlobalAssemble();
00606   if (err < 0) {
00607     return(err);
00608   }
00609 
00610   Epetra_Vector x(A->RowMap()), y(A->RowMap());
00611 
00612   x.PutScalar(1.0); y.PutScalar(0.0);
00613 
00614   Epetra_FECrsMatrix Acopy(*A);
00615 
00616   err = Acopy.GlobalAssemble();
00617   if (err < 0) {
00618     return(err);
00619   }
00620 
00621   bool the_same = epetra_test::compare_matrices(*A, Acopy);
00622   if (!the_same) {
00623     return(-1);
00624   }
00625 
00626   Epetra_FECrsMatrix Acopy2(Copy, A->RowMap(), A->ColMap(), 1);
00627 
00628   Acopy2 = Acopy;
00629 
00630   the_same = epetra_test::compare_matrices(*A, Acopy);
00631   if (!the_same) {
00632     return(-1);
00633   }
00634 
00635   int len = 20;
00636   long long* indices = new long long[len];
00637   double* values = new double[len];
00638   int numIndices;
00639 
00640   if (map.MyGID(0)) {
00641     EPETRA_CHK_ERR( A->ExtractGlobalRowCopy(0, len, numIndices,
00642               values, indices) );
00643     if (numIndices != 4) {
00644       return(-1);
00645     }
00646     if (indices[0] != 0) {
00647       return(-2);
00648     }
00649 
00650     if (values[0] != 1.0*numProcs) {
00651       cout << "ERROR: values[0] ("<<values[0]<<") should be "<<numProcs<<endl;
00652       return(-3);
00653     }
00654   }
00655 
00656   if (map.MyGID(4)) {
00657     EPETRA_CHK_ERR( A->ExtractGlobalRowCopy(4, len, numIndices,
00658               values, indices) );
00659 
00660     if (numIndices != 9) {
00661       return(-4);
00662     }
00663     int lcid = A->LCID(4);
00664     if (lcid<0) {
00665       return(-5);
00666     }
00667     if (values[lcid] != 4.0*numProcs) {
00668       cout << "ERROR: values["<<lcid<<"] ("<<values[lcid]<<") should be "
00669      <<4*numProcs<<endl;
00670       return(-6);
00671     }
00672   }
00673 
00674   delete [] values_2d;
00675   delete [] values_1d;
00676   delete [] nodes;
00677   delete [] indices;
00678   delete [] values;
00679 
00680   delete A;
00681   delete graph;
00682 
00683   return(0);
00684 }
00685 
00686 int submatrix_formats(const Epetra_Comm& Comm, bool verbose)
00687 {
00688   (void)verbose;
00689   //
00690   //This function simply verifies that the ROW_MAJOR/COLUMN_MAJOR switch works.
00691   //
00692   int numProcs = Comm.NumProc();
00693   int myPID = Comm.MyPID();
00694 
00695   int numLocalElements = 3;
00696   long long numGlobalElements = numLocalElements*numProcs;
00697   long long indexBase = 0;
00698 
00699   Epetra_Map map(numGlobalElements, numLocalElements, indexBase, Comm);
00700 
00701   Epetra_FECrsMatrix A(Copy, map, numLocalElements);
00702 
00703   Epetra_LongLongSerialDenseVector epetra_indices(numLocalElements);
00704 
00705   long long firstGlobalElement = numLocalElements*myPID;
00706 
00707   int i, j;
00708   for(i=0; i<numLocalElements; ++i) {
00709     epetra_indices[i] = firstGlobalElement+i;
00710   }
00711 
00712   Epetra_SerialDenseMatrix submatrix(numLocalElements, numLocalElements);
00713 
00714   for(i=0; i<numLocalElements; ++i) {
00715     for(j=0; j<numLocalElements; ++j) {
00716       submatrix(i,j) = 1.0*(firstGlobalElement+i);
00717     }
00718   }
00719 
00720   EPETRA_CHK_ERR( A.InsertGlobalValues(epetra_indices, submatrix,
00721                                Epetra_FECrsMatrix::COLUMN_MAJOR) );
00722 
00723   EPETRA_CHK_ERR( A.GlobalAssemble() );
00724 
00725   int len = 20;
00726   int numIndices;
00727   long long* indices = new long long[len];
00728   double* coefs = new double[len];
00729 
00730   for(i=0; i<numLocalElements; ++i) {
00731     long long row = firstGlobalElement+i;
00732 
00733     EPETRA_CHK_ERR( A.ExtractGlobalRowCopy(row, len, numIndices,
00734              coefs, indices) );
00735 
00736     for(j=0; j<numIndices; ++j) {
00737       if (coefs[j] != 1.0*row) {
00738   return(-2);
00739       }
00740     }
00741   }
00742 
00743   //now reset submatrix (transposing the i,j indices)
00744 
00745   for(i=0; i<numLocalElements; ++i) {
00746     for(j=0; j<numLocalElements; ++j) {
00747       submatrix(j,i) = 1.0*(firstGlobalElement+i);
00748     }
00749   }
00750 
00751   //sum these values into the matrix using the ROW_MAJOR switch, which should
00752   //result in doubling what's already there from the above Insert operation.
00753 
00754   EPETRA_CHK_ERR( A.SumIntoGlobalValues(epetra_indices, submatrix,
00755           Epetra_FECrsMatrix::ROW_MAJOR) );
00756 
00757   EPETRA_CHK_ERR( A.GlobalAssemble() );
00758 
00759   for(i=0; i<numLocalElements; ++i) {
00760     long long row = firstGlobalElement+i;
00761 
00762     EPETRA_CHK_ERR( A.ExtractGlobalRowCopy(row, len, numIndices,
00763              coefs, indices) );
00764 
00765     for(j=0; j<numIndices; ++j) {
00766       if (coefs[j] != 2.0*row) {
00767   return(-3);
00768       }
00769     }
00770   }
00771 
00772   delete [] indices;
00773   delete [] coefs;
00774 
00775   return(0);
00776 }
00777 
00778 int rectangular(const Epetra_Comm& Comm, bool verbose)
00779 {
00780   (void)verbose;
00781   int numprocs = Comm.NumProc();
00782   int localproc = Comm.MyPID();
00783   int numMyRows = 2;
00784   long long numGlobalRows = numprocs*numMyRows;
00785   long long* myrows = new long long[numMyRows];
00786 
00787   long long myFirstRow = ((long long)localproc)*numMyRows;
00788   int i;
00789   for(i=0; i<numMyRows; ++i) {
00790     myrows[i] = myFirstRow+i;
00791   }
00792 
00793   Epetra_Map map(numGlobalRows, numMyRows, myrows, 0LL, Comm);
00794 
00795   Epetra_FECrsMatrix A(Copy, map, 30);
00796 
00797   long long numcols = 20;
00798   long long* cols = new long long[(std::size_t) numcols];
00799   for(i=0; i<numcols; ++i) {
00800     cols[i] = i;
00801   }
00802 
00803   double* coefs = new double[(std::size_t) (numGlobalRows*numcols)];
00804   int offset = 0;
00805   for(int j=0; j<numcols; ++j) {
00806     for(i=0; i<numGlobalRows; ++i) {
00807       coefs[offset++] = 1.0*i;
00808     }
00809   }
00810 
00811   long long* globalRows = new long long[(std::size_t) numGlobalRows];
00812   for(i=0; i<numGlobalRows; ++i) globalRows[i] = i;
00813 
00814   EPETRA_CHK_ERR( A.InsertGlobalValues(numGlobalRows, globalRows,
00815                                        numcols, cols, coefs,
00816                                      Epetra_FECrsMatrix::COLUMN_MAJOR));
00817   delete [] coefs;
00818   delete [] globalRows;
00819 
00820   //Since the matrix is rectangular, we need to call GlobalAssemble with
00821   //a domain-map and a range-map. Otherwise, GlobalAssemble has no way of
00822   //knowing what the domain-map and range-map should be.
00823   //We'll use a linear distribution of the columns for a domain-map, and
00824   //our original row-map for the range-map.
00825   int numMyCols = (int) numcols/numprocs;
00826   int rem = (int) (numcols%numprocs);
00827   if (localproc<rem) ++numMyCols;
00828   Epetra_Map domainmap(numcols, numMyCols, 0LL, Comm);
00829 
00830   EPETRA_CHK_ERR( A.GlobalAssemble(domainmap, map) );
00831 
00832   long long numGlobalCols = A.NumGlobalCols64();
00833   long long numGlobalNNZ = A.NumGlobalNonzeros64();
00834 
00835   if (numGlobalCols != numcols ||
00836       numGlobalNNZ != numGlobalRows*numcols) {
00837     return(-1);
00838   }
00839 
00840   delete [] cols;
00841   delete [] myrows;
00842   return(0);
00843 }
00844 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines