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

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