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

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