Epetra Package Browser (Single Doxygen Collection) Development
FEVbrMatrix/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_Map.h"
00048 #include "Epetra_FEVbrMatrix.h"
00049 #include <Epetra_FEVector.h>
00050 
00051 #include "../epetra_test_err.h"
00052 
00053 
00054 int quad1(const Epetra_Map& map, bool verbose)
00055 {
00056   const Epetra_Comm & Comm = map.Comm();
00057   int numProcs = Comm.NumProc();
00058   int localProc = Comm.MyPID();
00059 
00060   Comm.Barrier();
00061   if (verbose && localProc == 0) {
00062     cout << "====================== quad1 =============================="<<endl;
00063   }
00064 
00065   //Set up a simple finite-element mesh containing 2-D quad elements, 1 per proc.
00066   //
00067   //   *-----*-----*-----*
00068   //  0|    2|    4|    6|
00069   //   | 0   | 1   | np-1|
00070   //   |     |     |     |
00071   //   *-----*-----*-----*
00072   //  1     3     5     7
00073   //
00074   //In the above drawing, 'np' means num-procs. node-numbers are to the
00075   //lower-left of each node (*).
00076   //
00077   //Each processor owns element 'localProc', and each processor owns
00078   //nodes 'localProc*2' and 'localProc*2+1' except for the last processor,
00079   //which also owns the last two nodes.
00080   //
00081   //There will be 3 degrees-of-freedom per node, so each element-matrix is
00082   //of size 12x12. (4 nodes per element, X 3 dof per node)
00083   //
00084 
00085   int myFirstNode = localProc*2;
00086   int myLastNode = localProc*2+1;
00087   if (localProc == numProcs-1) {
00088     myLastNode += 2;
00089   }
00090 
00091   int numMyNodes = myLastNode - myFirstNode + 1;
00092   int* myNodes = new int[numMyNodes];
00093   int i, j, ierr;
00094   for(i=0; i<numMyNodes; ++i) {
00095     myNodes[i] = myFirstNode + i;
00096   }
00097 
00098   int dofPerNode = 3; //degrees-of-freedom per node
00099   int indexBase = 0;
00100 
00101   Epetra_BlockMap blkMap(-1, numMyNodes, myNodes, dofPerNode,
00102        indexBase, Comm);
00103 
00104   int rowLengths = 3; //each element-matrix will have 4 block-columns.
00105                       //the rows of the assembled matrix will be longer than
00106                       //this, but we don't need to worry about that because the
00107                       //VbrMatrix will add memory as needed. For a real
00108                       //application where efficiency is a concern, better
00109                       //performance would be obtained by giving more accurate
00110                       //row-lengths here.
00111 
00112   Epetra_FEVbrMatrix A(Copy, blkMap, rowLengths);
00113 
00114   int nodesPerElem = 4;
00115   int* elemNodes = new int[nodesPerElem];
00116   for(i=0; i<nodesPerElem; ++i) elemNodes[i] = myFirstNode+i;
00117 
00118   int elemMatrixDim = nodesPerElem*dofPerNode;
00119   int len = elemMatrixDim*elemMatrixDim;
00120   double* elemMatrix = new double[len];
00121 
00122   //In an actual finite-element problem, we would calculate and fill 
00123   //meaningful element stiffness matrices. But for this simple matrix assembly
00124   //test, we're just going to fill our element matrix with 1.0's. This will
00125   //make it easy to see whether the matrix is correct after it's assembled.
00126 
00127   for(i=0; i<len; ++i) elemMatrix[i] = 1.0;
00128 
00129   //For filling in the matrix block-entries, we would ordinarily have to
00130   //carefully copy, or set pointers to, appropriate sections of the
00131   //element-matrix. But for this simple case we know that the element-matrix
00132   //is all 1's, so we'll just set our block-entry pointer to point to the
00133   //beginning of the element-matrix and leave it at that.
00134   //Note that the matrix class will refer to dofPerNode X dofPerNode (==9)
00135   //positions in the memory pointed to by 'blockEntry'.
00136 
00137   double* blockEntry = elemMatrix;
00138 
00139   //The element-matrix is a 4x4 (nodesPerElem X nodesPerElem) matrix of
00140   //3x3 block-entries. We'll now load our element-matrix into the global
00141   //matrix by looping over it and loading block-entries individually.
00142 
00143   for(i=0; i<nodesPerElem; ++i) {
00144     int blkrow = myFirstNode+i;
00145     EPETRA_TEST_ERR( A.BeginInsertGlobalValues(blkrow, nodesPerElem, elemNodes),
00146         ierr);
00147 
00148     for(j=0; j<nodesPerElem; ++j) {
00149       for(int ii=0; ii<dofPerNode*dofPerNode; ii++) {
00150         blockEntry[ii] = blkrow+elemNodes[j];
00151       }
00152       EPETRA_TEST_ERR( A.SubmitBlockEntry( blockEntry, dofPerNode,
00153             dofPerNode, dofPerNode), ierr);
00154     }
00155 
00156     int err = A.EndSubmitEntries();
00157     if (err < 0) {
00158       cout << "quad1: error in A.EndSubmitEntries: "<<err<<endl;
00159       return(-1);
00160     }
00161   }
00162 
00163   EPETRA_TEST_ERR( A.GlobalAssemble(), ierr);
00164 
00165   if (verbose && localProc==0) {
00166     cout << "after globalAssemble"<<endl;
00167   }
00168   if (verbose) {
00169     A.Print(cout);
00170   }
00171 
00172   int numMyRows = A.NumMyRows();
00173   int correct_numMyRows = dofPerNode*numMyNodes;
00174 
00175   if (numMyRows != correct_numMyRows) {
00176     cout << "proc " << localProc << ", numMyRows("<<numMyRows<<") doesn't match"
00177    << " correct_numMyRows("<<correct_numMyRows<<")."<<endl;
00178     return(-1);
00179   }
00180 
00181   int numMyNonzeros = A.NumMyNonzeros();
00182   int correct_numMyNonzeros = nodesPerElem*nodesPerElem*dofPerNode*dofPerNode;
00183 
00184   if (numProcs > 1) {
00185     if (localProc == numProcs-1) {
00186       correct_numMyNonzeros += dofPerNode*dofPerNode*4;
00187     }
00188     else if (localProc > 0) {
00189       correct_numMyNonzeros -= dofPerNode*dofPerNode*4;
00190     }
00191     else {
00192       //localProc==0 && numProcs > 1
00193       correct_numMyNonzeros -= dofPerNode*dofPerNode*8;
00194     }
00195   }
00196 
00197   if (numMyNonzeros != correct_numMyNonzeros) {
00198     cout << "proc " << localProc << ", numMyNonzeros(" << numMyNonzeros
00199    <<") != correct_numMyNonzeros("<<correct_numMyNonzeros<<")"<<endl;
00200     return(-1);
00201   }
00202 
00203 
00204   delete [] elemMatrix;
00205   delete [] myNodes;
00206   delete [] elemNodes;
00207 
00208   Comm.Barrier();
00209 
00210   return(0);
00211 }
00212 
00213 int quad2(const Epetra_Map& map, bool verbose)
00214 {
00215   const Epetra_Comm & Comm = map.Comm();
00216   int numProcs = Comm.NumProc();
00217   int localProc = Comm.MyPID();
00218   Comm.Barrier();
00219   if (verbose && localProc == 0) {
00220     cout << "====================== quad2 =============================="<<endl;
00221   }
00222 
00223   //Set up a simple finite-element mesh containing 2-D quad elements, 2 per proc.
00224   //(This test is similar to quad1() above, except for having 2 elements-per-proc
00225   // rather than 1.)
00226   //
00227   //   *-----*-----*-----*-------*
00228   //  0|    2|    4|    6|      8|
00229   //   | 0   | 1   | 2   | 2*np-1|
00230   //   |     |     |     |       |
00231   //   *-----*-----*-----*-------*
00232   //  1     3     5     7       9
00233   //
00234   //In the above drawing, 'np' means num-procs. node-numbers are to the
00235   //lower-left of each node (*).
00236   //
00237   //Each processor owns element 'localProc' and 'localProc+1', and each processor
00238   //owns nodes 'localProc*4' through 'localProc*4+3' except for the last
00239   //processor, which also owns the last two nodes in the mesh.
00240   //
00241   //There will be 3 degrees-of-freedom per node, so each element-matrix is
00242   //of size 12x12.
00243   //
00244 
00245   int myFirstNode = localProc*4;
00246   int myLastNode = localProc*4+3;
00247   if (localProc == numProcs-1) {
00248     myLastNode += 2;
00249   }
00250 
00251   int numMyElems = 2;
00252 
00253   int numMyNodes = myLastNode - myFirstNode + 1;
00254   int* myNodes = new int[numMyNodes];
00255   int i, j, ierr;
00256   for(i=0; i<numMyNodes; ++i) {
00257     myNodes[i] = myFirstNode + i;
00258   }
00259 
00260   int dofPerNode = 3; //degrees-of-freedom per node
00261   int indexBase = 0;
00262 
00263   Epetra_BlockMap blkMap(-1, numMyNodes, myNodes, dofPerNode,
00264        indexBase, Comm);
00265 
00266   int rowLengths = 4; //each element-matrix will have 4 block-columns.
00267                       //the rows of the assembled matrix will be longer than
00268                       //this, but we don't need to worry about that because the
00269                       //VbrMatrix will add memory as needed. For a real
00270                       //application where efficiency is a concern, better
00271                       //performance would be obtained by giving a more accurate
00272                       //row-length here.
00273 
00274   Epetra_FEVbrMatrix A(Copy, blkMap, rowLengths);
00275 
00276   int nodesPerElem = 4;
00277   int* elemNodes = new int[nodesPerElem];
00278 
00279   int elemMatrixDim = nodesPerElem*dofPerNode;
00280   int len = elemMatrixDim*elemMatrixDim;
00281   double* elemMatrix = new double[len];
00282 
00283   //In an actual finite-element problem, we would calculate and fill 
00284   //meaningful element stiffness matrices. But for this simple matrix assembly
00285   //test, we're just going to fill our element matrix with 1.0's. This will
00286   //make it easy to see whether the matrix is correct after it's assembled.
00287 
00288   for(i=0; i<len; ++i) elemMatrix[i] = 1.0;
00289 
00290   //For filling in the matrix block-entries, we would ordinarily have to
00291   //carefully copy, or set pointers to, appropriate sections of the
00292   //element-matrix. But for this simple case we know that the element-matrix
00293   //is all 1's, so we'll just set our block-entry pointer to point to the
00294   //beginning of the element-matrix and leave it at that.
00295   //Note that the matrix class will refer to dofPerNode X dofPerNode (==9)
00296   //positions in the memory pointed to by 'blockEntry'.
00297 
00298   double* blockEntry = elemMatrix;
00299 
00300   //Each element-matrix is a 4x4 (nodesPerElem X nodesPerElem) matrix of
00301   //3x3 block-entries. We'll now load our element-matrices into the global
00302   //matrix by looping over them and loading block-entries individually.
00303 
00304   int firstNode = myFirstNode;
00305   for(int el=0; el<numMyElems; ++el) {
00306     for(i=0; i<nodesPerElem; ++i) {
00307 
00308       for(int n=0; n<nodesPerElem; ++n) elemNodes[n] = firstNode+n;
00309 
00310       int blkrow = firstNode+i;
00311       EPETRA_TEST_ERR( A.BeginInsertGlobalValues(blkrow, nodesPerElem, elemNodes),
00312            ierr);
00313 
00314       for(j=0; j<nodesPerElem; ++j) {
00315   EPETRA_TEST_ERR( A.SubmitBlockEntry( blockEntry, dofPerNode,
00316                dofPerNode, dofPerNode), ierr);
00317       }
00318 
00319       int this_err = A.EndSubmitEntries();
00320       if (this_err < 0) {
00321   cerr << "error in quad2, A.EndSubmitEntries(): " << this_err << endl;
00322   return(this_err);
00323       }
00324     }
00325 
00326     firstNode += 2;
00327   }
00328 
00329   EPETRA_TEST_ERR( A.GlobalAssemble(), ierr);
00330 
00331   if (verbose && localProc==0) {
00332     cout << "after globalAssemble"<<endl;
00333   }
00334 
00335   if (verbose) {
00336     A.Print(cout);
00337   }
00338 
00339   //now let's make sure that we can perform a matvec...
00340   Epetra_FEVector x(blkMap, 1), y(blkMap, 1);
00341   x.PutScalar(1.0);
00342 
00343   EPETRA_TEST_ERR( A.Multiply(false, x, y), ierr);
00344 
00345   if (verbose && localProc==0) {
00346     cout << "quad2, y:"<<endl;
00347   }
00348   if (verbose) {
00349     y.Print(cout);
00350   }
00351 
00352   delete [] elemMatrix;
00353   delete [] myNodes;
00354   delete [] elemNodes;
00355 
00356   return(0);
00357 }
00358 
00359 int MultiVectorTests(const Epetra_Map & Map, int NumVectors, bool verbose)
00360 {
00361   const Epetra_Comm & Comm = Map.Comm();
00362   int ierr = 0, i, j;
00363   
00364   /* get number of processors and the name of this processor */
00365   
00366   int MyPID   = Comm.MyPID();
00367   
00368   // Construct FEVbrMatrix
00369   
00370   if (verbose && MyPID==0) cout << "constructing Epetra_FEVbrMatrix" << endl;
00371 
00372   //
00373   //we'll set up a tri-diagonal matrix.
00374   //
00375 
00376   int numGlobalRows = Map.NumGlobalElements();
00377   int minLocalRow = Map.MinMyGID();
00378   int rowLengths = 3;
00379 
00380   Epetra_FEVbrMatrix A(Copy, Map, rowLengths);
00381  
00382   if (verbose && MyPID==0) {
00383     cout << "calling A.InsertGlobalValues with 1-D data array"<<endl;
00384   }
00385 
00386   int numCols = 3;
00387   int* ptIndices = new int[numCols];
00388   for(int k=0; k<numCols; ++k) {
00389     ptIndices[k] = minLocalRow+k;
00390   }
00391 
00392   double* values_1d = new double[numCols*numCols];
00393   for(j=0; j<numCols*numCols; ++j) {
00394     values_1d[j] = 3.0;
00395   }
00396 
00397   //For an extreme test, we'll have all processors sum into all rows.
00398 
00399   int minGID = Map.MinAllGID();
00400 
00401   //For now we're going to assume that there's just one point associated with
00402   //each GID (element).
00403 
00404   double* ptCoefs = new double[3];
00405 
00406   {for(i=0; i<numGlobalRows; ++i) {
00407     if (i>0 && i<numGlobalRows-1) {
00408       ptIndices[0] = minGID+i-1;
00409       ptIndices[1] = minGID+i;
00410       ptIndices[2] = minGID+i+1;
00411       ptCoefs[0] = -1.0;
00412       ptCoefs[1] = 2.0;
00413       ptCoefs[2] = -1.0;
00414       numCols = 3;
00415     }
00416     else if (i == 0) {
00417       ptIndices[0] = minGID+i;
00418       ptIndices[1] = minGID+i+1;
00419       ptIndices[2] = minGID+i+2;
00420       ptCoefs[0] = 2.0;
00421       ptCoefs[1] = -1.0;
00422       ptCoefs[2] = -1.0;
00423       numCols = 3;
00424     }
00425     else {
00426       ptIndices[0] = minGID+i-2;
00427       ptIndices[1] = minGID+i-1;
00428       ptIndices[2] = minGID+i;
00429       ptCoefs[0] = -1.0;
00430       ptCoefs[1] = -1.0;
00431       ptCoefs[2] = 2.0;
00432       numCols = 3;
00433     }
00434 
00435     int row = minGID+i;
00436 
00437     EPETRA_TEST_ERR( A.BeginInsertGlobalValues(row, rowLengths, ptIndices), ierr);
00438 
00439     for(j=0; j<rowLengths; ++j) {
00440       EPETRA_TEST_ERR( A.SubmitBlockEntry(&(ptCoefs[j]), 1, 1, 1), ierr);
00441     }
00442 
00443     EPETRA_TEST_ERR( A.EndSubmitEntries(), ierr);
00444 
00445   }}
00446 
00447   if (verbose&&MyPID==0) {
00448     cout << "calling A.GlobalAssemble()" << endl;
00449   }
00450 
00451   EPETRA_TEST_ERR( A.GlobalAssemble(), ierr );
00452 
00453   if (verbose&&MyPID==0) {
00454   cout << "after globalAssemble"<<endl;
00455   }
00456   if (verbose) {
00457   A.Print(cout);
00458   }
00459 
00460   delete [] values_1d;
00461   delete [] ptIndices;
00462   delete [] ptCoefs;
00463 
00464   return(ierr);
00465 }
00466 
00467 int four_quads(const Epetra_Comm& Comm, bool preconstruct_graph, bool verbose)
00468 {   
00469   if (verbose) {
00470     cout << "******************* four_quads ***********************"<<endl;
00471   }
00472 
00473   //This function assembles a matrix representing a finite-element mesh
00474   //of four 2-D quad elements. There are 9 nodes in the problem. The
00475   //same problem is assembled no matter how many processors are being used
00476   //(within reason). It may not work if more than 9 processors are used.
00477   //
00478   //  *------*------*
00479   // 6|     7|     8|
00480   //  | E2   | E3   |
00481   //  *------*------*
00482   // 3|     4|     5|
00483   //  | E0   | E1   |
00484   //  *------*------*
00485   // 0      1      2
00486   //
00487   //Nodes are denoted by * with node-numbers below and left of each node.
00488   //E0, E1 and so on are element-numbers.
00489   //
00490   //Each processor will contribute a sub-matrix of size 4x4, filled with 1's,
00491   //for each element. Thus, the coefficient value at position 0,0 should end up
00492   //being 1.0*numProcs, the value at position 4,4 should be 1.0*4*numProcs, etc.
00493   //
00494   //Depending on the number of processors being used, the locations of the
00495   //specific matrix positions (in terms of which processor owns them) will vary.
00496   //
00497   
00498   int numProcs = Comm.NumProc();
00499   
00500   int numNodes = 9;
00501   int numElems = 4;
00502   int numNodesPerElem = 4;
00503 
00504   int blockSize = 1;
00505   int indexBase = 0;
00506 
00507   //Create a map using epetra-defined linear distribution.
00508   Epetra_BlockMap map(numNodes, blockSize, indexBase, Comm);
00509 
00510   Epetra_CrsGraph* graph = NULL;
00511 
00512   int* nodes = new int[numNodesPerElem];
00513   int i, j, k, err = 0;
00514 
00515   if (preconstruct_graph) {
00516     graph = new Epetra_CrsGraph(Copy, map, 1);
00517 
00518     //we're going to fill the graph with indices, but remember it will only
00519     //accept indices in rows for which map.MyGID(row) is true.
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         break;
00526       case 1:
00527         nodes[0] = 1; nodes[1] = 2; nodes[2] = 5; nodes[3] = 4;
00528         break;
00529       case 2:
00530         nodes[0] = 3; nodes[1] = 4; nodes[2] = 7; nodes[3] = 6;
00531         break;
00532       case 3:
00533         nodes[0] = 4; nodes[1] = 5; nodes[2] = 8; nodes[3] = 7;
00534         break;
00535       }
00536 
00537       for(j=0; j<numNodesPerElem; ++j) {
00538         if (map.MyGID(nodes[j])) {
00539           err = graph->InsertGlobalIndices(nodes[j], numNodesPerElem,
00540                                            nodes);
00541           if (err<0) return(err);
00542         }
00543       }
00544     }
00545 
00546     EPETRA_CHK_ERR( graph->FillComplete() );
00547   }
00548 
00549   Epetra_FEVbrMatrix* A = NULL;
00550 
00551   if (preconstruct_graph) {
00552     A = new Epetra_FEVbrMatrix(Copy, *graph);
00553   }
00554   else {
00555     A = new Epetra_FEVbrMatrix(Copy, map, 1);
00556   }
00557 
00558   //EPETRA_CHK_ERR( A->PutScalar(0.0) );
00559 
00560   double* values_1d = new double[numNodesPerElem*numNodesPerElem];
00561   double** values_2d = new double*[numNodesPerElem];
00562 
00563   for(i=0; i<numNodesPerElem*numNodesPerElem; ++i) values_1d[i] = 1.0;
00564 
00565   int offset = 0;
00566   for(i=0; i<numNodesPerElem; ++i) {
00567     values_2d[i] = &(values_1d[offset]);
00568     offset += numNodesPerElem;
00569   }
00570 
00571   for(i=0; i<numElems; ++i) {
00572     switch(i) {
00573     case 0:
00574       nodes[0] = 0; nodes[1] = 1; nodes[2] = 4; nodes[3] = 3;
00575       break;
00576 
00577     case 1:
00578       nodes[0] = 1; nodes[1] = 2; nodes[2] = 5; nodes[3] = 4;
00579       break;
00580 
00581     case 2:
00582       nodes[0] = 3; nodes[1] = 4; nodes[2] = 7; nodes[3] = 6;
00583       break;
00584 
00585      case 3:
00586       nodes[0] = 4; nodes[1] = 5; nodes[2] = 8; nodes[3] = 7;
00587       break;
00588     }
00589 
00590     for(j=0; j<numNodesPerElem; ++j) {
00591       if (preconstruct_graph) {
00592   err = A->BeginSumIntoGlobalValues(nodes[j], numNodesPerElem, nodes);
00593   if (err<0) return(err);
00594       }
00595       else {
00596   err = A->BeginInsertGlobalValues(nodes[j], numNodesPerElem, nodes);
00597   if (err<0) return(err);
00598       }
00599     
00600       for(k=0; k<numNodesPerElem; ++k) {
00601   err = A->SubmitBlockEntry(values_1d, blockSize, blockSize, blockSize);
00602   if (err<0) return(err);
00603       }
00604 
00605       err = A->EndSubmitEntries();
00606       if (err<0) return(err);
00607     }
00608   }
00609 
00610   EPETRA_CHK_ERR( A->GlobalAssemble() );
00611 
00612   Epetra_FEVbrMatrix* Acopy = new Epetra_FEVbrMatrix(*A);
00613 
00614   if (verbose) {
00615     cout << "A:"<<*A << endl;
00616     cout << "Acopy:"<<*Acopy<<endl;
00617   }
00618 
00619   Epetra_Vector x(A->RowMap()), y(A->RowMap());
00620 
00621   x.PutScalar(1.0); y.PutScalar(0.0);
00622 
00623   Epetra_Vector x2(Acopy->RowMap()), y2(Acopy->RowMap());
00624 
00625   x2.PutScalar(1.0); y2.PutScalar(0.0);
00626 
00627   A->Multiply(false, x, y);
00628 
00629   Acopy->Multiply(false, x2, y2);
00630 
00631   double ynorm2, y2norm2;
00632 
00633   y.Norm2(&ynorm2);
00634   y2.Norm2(&y2norm2);
00635   if (ynorm2 != y2norm2) {
00636     cerr << "norm2(A*ones) != norm2(*Acopy*ones)"<<endl;
00637     return(-99);
00638   }
00639 
00640   Epetra_FEVbrMatrix* Acopy2 =
00641     new Epetra_FEVbrMatrix(Copy, A->RowMap(), A->ColMap(), 1);
00642 
00643   *Acopy2 = *Acopy;
00644 
00645   Epetra_Vector x3(Acopy->RowMap()), y3(Acopy->RowMap());
00646 
00647   x3.PutScalar(1.0); y3.PutScalar(0.0);
00648 
00649   Acopy2->Multiply(false, x3, y3);
00650 
00651   double y3norm2;
00652   y3.Norm2(&y3norm2);
00653 
00654   if (y3norm2 != y2norm2) {
00655     cerr << "norm2(Acopy*ones) != norm2(Acopy2*ones)"<<endl;
00656     return(-999);
00657   }
00658 
00659   int len = 20;
00660   int* indices = new int[len];
00661   double* values = new double[len];
00662   int numIndices;
00663 
00664   if (map.MyGID(0)) {
00665     int lid = map.LID(0);
00666     EPETRA_CHK_ERR( A->ExtractMyRowCopy(lid, len, numIndices,
00667           values, indices) );
00668     if (numIndices != 4) {
00669       return(-1);
00670     }
00671     if (indices[0] != lid) {
00672       return(-2);
00673     }
00674 
00675     if (values[0] != 1.0*numProcs) {
00676       cout << "ERROR: values[0] ("<<values[0]<<") should be "<<numProcs<<endl;
00677       return(-3);
00678     }
00679   }
00680 
00681   if (map.MyGID(4)) {
00682     int lid = map.LID(4);
00683     EPETRA_CHK_ERR( A->ExtractMyRowCopy(lid, len, numIndices,
00684           values, indices) );
00685 
00686     if (numIndices != 9) {
00687       return(-4);
00688     }
00689     int lcid = A->LCID(4);
00690 //     if (indices[lcid] != 4) {
00691 //       cout << "ERROR: indices[4] ("<<indices[4]<<") should be "
00692 //     <<A->LCID(4)<<endl;
00693 //       return(-5);
00694 //     }
00695     if (values[lcid] != 4.0*numProcs) {
00696       cout << "ERROR: values["<<lcid<<"] ("<<values[lcid]<<") should be "
00697      <<4*numProcs<<endl;
00698       return(-6);
00699     }
00700   }
00701 
00702   delete [] values_2d;
00703   delete [] values_1d;
00704   delete [] nodes;
00705   delete [] indices;
00706   delete [] values;
00707 
00708   delete A;
00709   delete Acopy2;
00710   delete Acopy;
00711   delete graph;
00712 
00713   return(0);
00714 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines