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

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