test/CrsGraph/cxx_main.cpp

Go to the documentation of this file.

/*
//@HEADER
// ************************************************************************
//
//               Epetra: Linear Algebra Services Package
//                 Copyright 2011 Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
//
// ************************************************************************
//@HEADER
*/

#include "Epetra_ConfigDefs.h"
#include "Epetra_CrsGraph.h"
#include "Epetra_Map.h"
#ifdef EPETRA_MPI
#include "Epetra_MpiComm.h"
#include <mpi.h>
#else
#include "Epetra_SerialComm.h"
#endif
#include "../epetra_test_err.h"
#include "Epetra_Version.h"

// Prototype
int check(Epetra_CrsGraph& A, int NumMyRows1, int NumGlobalRows1, int NumMyNonzeros1,
          int NumGlobalNonzeros1, int* MyGlobalElements, bool verbose);

int checkSharedOwnership(Epetra_Comm& Comm, bool verbose);
int checkCopyAndAssignment(Epetra_Comm& Comm, bool verbose);

int main(int argc, char *argv[]) {
  int ierr = 0;
  int i;
  int forierr = 0;
  int* Indices;
  bool debug = true;

#ifdef EPETRA_MPI

  // Initialize MPI

  MPI_Init(&argc,&argv);
  int rank; // My process ID

  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
  Epetra_MpiComm Comm( MPI_COMM_WORLD );

#else

  int rank = 0;
  Epetra_SerialComm Comm;

#endif

  bool verbose = false;

  // Check if we should print results to standard out
  if(argc > 1) {
    if(argv[1][0]=='-' && argv[1][1]=='v') {
      verbose = true;
    }
  }

  //char tmp;
  //if (rank==0) cout << "Press any key to continue..."<< endl;
  //if (rank==0) cin >> tmp;
  //Comm.Barrier();

  Comm.SetTracebackMode(0); // This should shut down any error traceback reporting
  int MyPID = Comm.MyPID();
  int NumProc = Comm.NumProc();

  if(verbose && MyPID==0)
    cout << Epetra_Version() << endl << endl;

  if(verbose) cout << "Processor "<<MyPID<<" of "<< NumProc << " is alive." << endl;

  bool verbose1 = verbose;

  // Redefine verbose to only print on PE 0
  if(verbose && rank != 0)
    verbose = false;

  int NumMyEquations = 5;
  int NumGlobalEquations = NumMyEquations*NumProc+EPETRA_MIN(NumProc,3);
  if(MyPID < 3)
    NumMyEquations++;

  // Construct a Map that puts approximately the same Number of equations on each processor

  Epetra_Map& Map = *new Epetra_Map(NumGlobalEquations, NumMyEquations, 0, Comm);

  // Get update list and number of local equations from newly created Map
  int* MyGlobalElements = new int[Map.NumMyElements()];
  Map.MyGlobalElements(MyGlobalElements);

  // Create an integer vector NumNz that is used to build the Petra Matrix.
  // NumNz[i] is the Number of OFF-DIAGONAL term for the ith global equation on this processor

  int* NumNz = new int[NumMyEquations];

  // We are building a tridiagonal matrix where each row has (-1 2 -1)
  // So we need 2 off-diagonal terms (except for the first and last equation)

  for(i=0; i<NumMyEquations; i++)
    if(MyGlobalElements[i]==0 || MyGlobalElements[i] == NumGlobalEquations-1)
      NumNz[i] = 1;
    else
      NumNz[i] = 2;

  // Create a Epetra_CrsGraph

  Epetra_CrsGraph& A = *new Epetra_CrsGraph(Copy, Map, NumNz);
  EPETRA_TEST_ERR(A.IndicesAreGlobal(),ierr);
  EPETRA_TEST_ERR(A.IndicesAreLocal(),ierr);

  // Add  rows one-at-a-time
  // Need some vectors to help
  // Off diagonal Values will always be -1

  Indices = new int[2];
  int NumEntries;

  forierr = 0;
  for(i = 0; i < NumMyEquations; i++) {
    if(MyGlobalElements[i] == 0) {
      Indices[0] = 1;
      NumEntries = 1;
    }
    else if(MyGlobalElements[i] == NumGlobalEquations-1) {
      Indices[0] = NumGlobalEquations-2;
      NumEntries = 1;
    }
    else {
      Indices[0] = MyGlobalElements[i]-1;
      Indices[1] = MyGlobalElements[i]+1;
      NumEntries = 2;
    }
    forierr += !(A.InsertGlobalIndices(MyGlobalElements[i], NumEntries, Indices)==0);
    forierr += !(A.InsertGlobalIndices(MyGlobalElements[i], 1, MyGlobalElements+i)>0); // Put in the diagonal entry (should cause realloc)
  }
  EPETRA_TEST_ERR(forierr,ierr);
  //A.PrintGraphData(cout);
  delete[] Indices;

  int gRID = A.GRID(0);
  int numIndices = A.NumGlobalIndices(gRID);
  std::vector<int> indices_vec(numIndices);
  A.ExtractGlobalRowCopy(gRID, numIndices, numIndices, &indices_vec[0]);
  A.RemoveGlobalIndices(gRID);
  EPETRA_TEST_ERR(!(A.NumGlobalIndices(gRID)==0),ierr);
  if (ierr != 0) cout << "tests FAILED" << std::endl;

  A.InsertGlobalIndices(gRID, numIndices, &indices_vec[0]);
  EPETRA_TEST_ERR(!(A.NumGlobalIndices(gRID)==numIndices),ierr);
  if (ierr != 0) cout << "tests FAILED" << std::endl;

  // Finish up
  EPETRA_TEST_ERR(!(A.IndicesAreGlobal()),ierr);
  EPETRA_TEST_ERR(!(A.FillComplete()==0),ierr);
  EPETRA_TEST_ERR(!(A.IndicesAreLocal()),ierr);
  EPETRA_TEST_ERR(A.StorageOptimized(),ierr);

  A.OptimizeStorage();

  EPETRA_TEST_ERR(!(A.StorageOptimized()),ierr);
  EPETRA_TEST_ERR(A.UpperTriangular(),ierr);
  EPETRA_TEST_ERR(A.LowerTriangular(),ierr);

  if (ierr != 0) cout << "tests FAILED" << std::endl;

  EPETRA_TEST_ERR(A.NumMyDiagonals()-NumMyEquations,ierr);
  if (ierr != 0) cout << "tests FAILED" << std::endl;
  EPETRA_TEST_ERR(A.NumMyBlockDiagonals()-NumMyEquations,ierr);
  if (ierr != 0) cout << "tests FAILED" << std::endl;

  EPETRA_TEST_ERR(A.NumGlobalDiagonals()-NumGlobalEquations,ierr);
  if (ierr != 0) cout << "tests FAILED" << std::endl;
  EPETRA_TEST_ERR(A.NumGlobalBlockDiagonals()-NumGlobalEquations,ierr);
  if (ierr != 0) cout << "tests FAILED" << std::endl;

  if(verbose) cout << "\n*****Testing variable entry constructor\n" << endl;

  int NumMyNonzeros = 3 * NumMyEquations;
  if(A.LRID(0) >= 0)
    NumMyNonzeros--; // If I own first global row, then there is one less nonzero
  if(A.LRID(NumGlobalEquations-1) >= 0)
    NumMyNonzeros--; // If I own last global row, then there is one less nonzero

  EPETRA_TEST_ERR(check(A, NumMyEquations, NumGlobalEquations, NumMyNonzeros, 3*NumGlobalEquations-2,
                        MyGlobalElements, verbose),ierr);
  forierr = 0;
  for(i = 0; i < NumMyEquations; i++)
    forierr += !(A.NumGlobalIndices(MyGlobalElements[i])==NumNz[i]+1);
  EPETRA_TEST_ERR(forierr,ierr);
  for(i = 0; i < NumMyEquations; i++)
    forierr += !(A.NumMyIndices(i)==NumNz[i]+1);
  EPETRA_TEST_ERR(forierr,ierr);

  if(verbose) cout << "NumIndices function check OK" << endl;

  delete &A;

  if(debug) Comm.Barrier();

  if(verbose) cout << "\n*****Testing constant entry constructor\n" << endl;

  Epetra_CrsGraph& AA = *new Epetra_CrsGraph(Copy, Map, 5);

  if(debug) Comm.Barrier();

  for(i = 0; i < NumMyEquations; i++)
    AA.InsertGlobalIndices(MyGlobalElements[i], 1, MyGlobalElements+i);

  // Note:  All processors will call the following Insert routines, but only the processor
  //        that owns it will actually do anything

  int One = 1;
  if(AA.MyGlobalRow(0)) {
    EPETRA_TEST_ERR(!(AA.InsertGlobalIndices(0, 0, &One)==0),ierr);
  }
  else
    EPETRA_TEST_ERR(!(AA.InsertGlobalIndices(0, 1, &One)==-2),ierr);
  EPETRA_TEST_ERR(!(AA.FillComplete()==0),ierr);
  EPETRA_TEST_ERR(AA.StorageOptimized(),ierr);
  EPETRA_TEST_ERR(!(AA.UpperTriangular()),ierr);
  EPETRA_TEST_ERR(!(AA.LowerTriangular()),ierr);

  if(debug) Comm.Barrier();
  EPETRA_TEST_ERR(check(AA, NumMyEquations, NumGlobalEquations, NumMyEquations, NumGlobalEquations,
                        MyGlobalElements, verbose),ierr);

  if(debug) Comm.Barrier();

  forierr = 0;
  for(i = 0; i < NumMyEquations; i++)
    forierr += !(AA.NumGlobalIndices(MyGlobalElements[i])==1);
  EPETRA_TEST_ERR(forierr,ierr);

  if(verbose) cout << "NumIndices function check OK" << endl;

  if(debug) Comm.Barrier();

  if(verbose) cout << "\n*****Testing copy constructor\n" << endl;

  Epetra_CrsGraph& B = *new Epetra_CrsGraph(AA);
  delete &AA;

  EPETRA_TEST_ERR(check(B, NumMyEquations, NumGlobalEquations, NumMyEquations, NumGlobalEquations,
                        MyGlobalElements, verbose),ierr);

  forierr = 0;
  for(i = 0; i < NumMyEquations; i++)
    forierr += !(B.NumGlobalIndices(MyGlobalElements[i])==1);
  EPETRA_TEST_ERR(forierr,ierr);

  if(verbose) cout << "NumIndices function check OK" << endl;

  if(debug) Comm.Barrier();

  if(verbose) cout << "\n*****Testing post construction modifications\n" << endl;

  EPETRA_TEST_ERR(!(B.InsertGlobalIndices(0, 1, &One)==-2),ierr);
  delete &B;

  // Release all objects
  delete[] MyGlobalElements;
  delete[] NumNz;
  delete &Map;
      

  if (verbose1) {
    // Test ostream << operator (if verbose1)
    // Construct a Map that puts 2 equations on each PE

    int NumMyElements1 = 4;
    int NumMyEquations1 = NumMyElements1;
    int NumGlobalEquations1 = NumMyEquations1*NumProc;

    Epetra_Map& Map1 = *new Epetra_Map(NumGlobalEquations1, NumMyElements1, 1, Comm);

    // Get update list and number of local equations from newly created Map
    int* MyGlobalElements1 = new int[Map1.NumMyElements()];
    Map1.MyGlobalElements(MyGlobalElements1);

    // Create an integer vector NumNz that is used to build the Petra Matrix.
    // NumNz[i] is the Number of OFF-DIAGONAL term for the ith global equation on this processor

    int* NumNz1 = new int[NumMyEquations1];

    // We are building a tridiagonal matrix where each row has (-1 2 -1)
    // So we need 2 off-diagonal terms (except for the first and last equation)

    for(i = 0; i < NumMyEquations1; i++)
      if(MyGlobalElements1[i]==1 || MyGlobalElements1[i] == NumGlobalEquations1)
        NumNz1[i] = 1;
      else
        NumNz1[i] = 2;

    // Create a Epetra_Graph using 1-based arithmetic

    Epetra_CrsGraph& A1 = *new Epetra_CrsGraph(Copy, Map1, NumNz1);

    // Add  rows one-at-a-time
    // Need some vectors to help
    // Off diagonal Values will always be -1


    int* Indices1 = new int[2];
    int NumEntries1;

    forierr = 0;
    for(i = 0; i < NumMyEquations1; i++) {
      if(MyGlobalElements1[i]==1) {
        Indices1[0] = 2;
        NumEntries1 = 1;
      }
      else if(MyGlobalElements1[i] == NumGlobalEquations1) {
      Indices1[0] = NumGlobalEquations1-1;
      NumEntries1 = 1;
      }
      else {
        Indices1[0] = MyGlobalElements1[i]-1;
        Indices1[1] = MyGlobalElements1[i]+1;
        NumEntries1 = 2;
      }
      forierr += !(A1.InsertGlobalIndices(MyGlobalElements1[i], NumEntries1, Indices1)==0);
      forierr += !(A1.InsertGlobalIndices(MyGlobalElements1[i], 1, MyGlobalElements1+i)>0); // Put in the diagonal entry
    }
    EPETRA_TEST_ERR(forierr,ierr);
    
    // Finish up
    EPETRA_TEST_ERR(!(A1.FillComplete()==0),ierr);

    if(verbose) cout << "Print out tridiagonal matrix, each part on each processor. Index base is one.\n" << endl;
    cout << A1 << endl;

  // Release all objects
  delete[] NumNz1;
  delete[] Indices1;
  delete[] MyGlobalElements1;

  delete &A1;
  delete &Map1;
  }

  // Test copy constructor, op=, and reference-counting
  int tempierr = 0;
  if(verbose) cout << "\n*****Checking cpy ctr, op=, and reference counting." << endl;
  tempierr = checkCopyAndAssignment(Comm, verbose);
  EPETRA_TEST_ERR(tempierr, ierr);
  if(verbose && (tempierr == 0)) cout << "Checked OK." << endl;

  // Test shared-ownership code (not implemented yet)
  tempierr = 0;
  if(verbose) cout << "\n*****Checking shared-ownership tests." << endl;
  tempierr = checkSharedOwnership(Comm, verbose);
  EPETRA_TEST_ERR(tempierr, ierr);
  if(verbose && (tempierr == 0)) cout << "Checked OK." << endl;

      
#ifdef EPETRA_MPI
  MPI_Finalize() ;
#endif
/* end main */
  return(ierr);
}

//==============================================================================
int checkSharedOwnership(Epetra_Comm& Comm, bool verbose) {
  // check to make sure each function returns 1 when it should
  // check to make sure each function doesn't return 1 when it shouldn't
  int ierr = 0;

  // initialize Map
  const int NumMyElements = 10;
  const int IndexBase = 0;
  Epetra_Map Map1(-1, NumMyElements, IndexBase, Comm);
  // initialize Graphs
  const int NumIndicesPerRow = 5;
  Epetra_CrsGraph * SoleOwner = new Epetra_CrsGraph(Copy, Map1, Map1, NumIndicesPerRow);
  Epetra_CrsGraph SharedOrig(Copy, Map1, Map1, NumIndicesPerRow);
  Epetra_CrsGraph SharedOwner(SharedOrig);
  // arrays used by Insert & Remove
  Epetra_IntSerialDenseVector array1(2);
  array1[0] = NumIndicesPerRow / 2;
  array1[1] = array1[0] + 1;
  Epetra_IntSerialDenseVector array2(NumIndicesPerRow);
  for(int i = 0; i < NumIndicesPerRow; i++)
    array2[i] = i;
  // output variables (declaring them here lets us comment out indiv. tests)
  int soleOutput, sharedOutput;

  // InsertMyIndices
  if(verbose) cout << "InsertMyIndices..." << endl;
  soleOutput = SoleOwner->InsertMyIndices(0, 2, array1.Values());
  sharedOutput = SharedOwner.InsertMyIndices(0, 2, array1.Values());
  EPETRA_TEST_ERR(!(soleOutput == 0), ierr);
  EPETRA_TEST_ERR(!(sharedOutput == 1), ierr);
  if(verbose && ierr > 0) cout << "soleOutput = " << soleOutput << " sharedOutput = " << sharedOutput << endl;

  // RemoveMyIndices (#0)
  if(verbose) cout << "RemoveMyIndices(#0)..." << endl;
  soleOutput = SoleOwner->RemoveMyIndices(0);
  EPETRA_TEST_ERR(!(soleOutput == 0), ierr);

  EPETRA_TEST_ERR(!(SoleOwner->NumMyIndices(0)==0),ierr);
  if (ierr != 0) cout << "tests FAILED" << std::endl;

  soleOutput = SoleOwner->InsertMyIndices(0, 2, array1.Values());
  EPETRA_TEST_ERR(!(soleOutput == 0), ierr);

  // SortIndices
  //if(verbose) cout << "SortIndices..." << endl;
  //soleOutput = SoleOwner.SortIndices();
  //sharedOutput = SharedOwner.SortIndices();
  //EPETRA_TEST_ERR(!(soleOutput == 0), ierr);
  //EPETRA_TEST_ERR(!(sharedOutput == 1), ierr);
  //if(verbose && ierr > 0) cout << "soleOutput = " << soleOutput << " sharedOutput = " << sharedOutput << endl;

  // RemoveRedundantIndices
  //if(verbose) cout << "RemoveRedundantIndices..." << endl;
  //SoleOwner.InsertGlobalIndices(0, 1, array1.Values());
  //SharedOwner.InsertGlobalIndices(0, 1, array1.Values());
  //soleOutput = SoleOwner.RemoveRedundantIndices();
  //sharedOutput = SharedOwner.RemoveRedundantIndices();
  //EPETRA_TEST_ERR(!(soleOutput == 0), ierr);
  //EPETRA_TEST_ERR(!(sharedOutput == 1), ierr);
  //if(verbose && ierr > 0) cout << "soleOutput = " << soleOutput << " sharedOutput = " << sharedOutput << endl;

  // FillComplete (#1)
  if(verbose) cout << "FillComplete..." << endl;
  soleOutput = SoleOwner->FillComplete();
  sharedOutput = SharedOwner.FillComplete();
  EPETRA_TEST_ERR(!(soleOutput == 0), ierr);
  EPETRA_TEST_ERR(!(sharedOutput == 1), ierr);
  if(verbose && ierr > 0) cout << "soleOutput = " << soleOutput << " sharedOutput = " << sharedOutput << endl;

  // OptimizeStorage
  if(verbose) cout << "OptimizeStorage..." << endl;
  soleOutput = SoleOwner->OptimizeStorage();
  sharedOutput = SharedOwner.OptimizeStorage();
  EPETRA_TEST_ERR(!(soleOutput == 0), ierr);
  EPETRA_TEST_ERR(!(sharedOutput == 0), ierr);
  if(verbose && ierr > 0) cout << "soleOutput = " << soleOutput << " sharedOutput = " << sharedOutput << endl;

  // RemoveMyIndices (#1)
  if(verbose) cout << "RemoveMyIndices..." << endl;
  soleOutput = SoleOwner->RemoveMyIndices(0, 1, &array1[1]);
  sharedOutput = SharedOwner.RemoveMyIndices(0, 1, &array1[1]);
  EPETRA_TEST_ERR(!(soleOutput == -1), ierr);
  EPETRA_TEST_ERR(!(sharedOutput == -1), ierr);
  if(verbose && ierr > 0) cout << "soleOutput = " << soleOutput << " sharedOutput = " << sharedOutput << endl;

  // RemoveMyIndices (#2)
  if(verbose) cout << "RemoveMyIndices(#2)..." << endl;
  soleOutput = SoleOwner->RemoveMyIndices(0);
  sharedOutput = SharedOwner.RemoveMyIndices(0);
  EPETRA_TEST_ERR(!(soleOutput == -1), ierr);
  EPETRA_TEST_ERR(!(sharedOutput == -1), ierr);
  if(verbose && ierr > 0) cout << "soleOutput = " << soleOutput << " sharedOutput = " << sharedOutput << endl;

  // FillComplete (#2)
  if(verbose) cout << "FillComplete(#2)..." << endl;
  soleOutput = SoleOwner->FillComplete(SoleOwner->DomainMap(), SoleOwner->RangeMap());
  sharedOutput = SharedOwner.FillComplete(SharedOwner.DomainMap(), SharedOwner.RangeMap());
  EPETRA_TEST_ERR(!(soleOutput == 0), ierr);
  EPETRA_TEST_ERR(!(sharedOutput == 1), ierr);
  if(verbose && ierr > 0) cout << "soleOutput = " << soleOutput << " sharedOutput = " << sharedOutput << endl;

  {
    // make new Graphs so that we can insert Global instead of Local
    // inside of new scope so that we can use same names
    Epetra_CrsGraph SoleOwnerG(Copy, Map1, NumIndicesPerRow);
    Epetra_CrsGraph SharedOrigG(Copy, Map1, NumIndicesPerRow);
    Epetra_CrsGraph SharedOwnerG(SharedOrig);
    
    int GlobalRow = SoleOwnerG.GRID(0);

    // InsertGlobalIndices
    if(verbose) cout << "InsertGlobalIndices..." << endl;
    soleOutput = SoleOwnerG.InsertGlobalIndices(GlobalRow, 2, array2.Values());
    sharedOutput = SharedOwnerG.InsertGlobalIndices(GlobalRow, 2, array2.Values());
    EPETRA_TEST_ERR(!(soleOutput == 0), ierr);
    EPETRA_TEST_ERR(!(sharedOutput == 1), ierr);
    if(verbose && ierr > 0) cout << "soleOutput = " << soleOutput << " sharedOutput = " << sharedOutput << endl;
    
    // RemoveGlobalIndices (#1)
    if(verbose) cout << "RemoveGlobalIndices..." << endl;
    soleOutput = SoleOwnerG.RemoveGlobalIndices(GlobalRow, 1, &array2[1]);
    sharedOutput = SharedOwnerG.RemoveGlobalIndices(GlobalRow, 1, &array2[1]);
    EPETRA_TEST_ERR(!(soleOutput == 0), ierr);
    EPETRA_TEST_ERR(!(sharedOutput == 1), ierr);
    if(verbose && ierr > 0) cout << "soleOutput = " << soleOutput << " sharedOutput = " << sharedOutput << endl;
    
    // RemoveGlobalIndices (#2)
    if(verbose) cout << "RemoveGlobalIndices(#2)..." << endl;
    soleOutput = SoleOwnerG.RemoveGlobalIndices(GlobalRow);
    sharedOutput = SharedOwnerG.RemoveGlobalIndices(GlobalRow);
    EPETRA_TEST_ERR(!(soleOutput == 0), ierr);
    EPETRA_TEST_ERR(!(sharedOutput == 1), ierr);
    if(verbose && ierr > 0) cout << "soleOutput = " << soleOutput << " sharedOutput = " << sharedOutput << endl;
    }


  // *PROT* InsertIndices
  // *PROT* MakeIndicesLocal
  delete SoleOwner;
  return(ierr);

}

//==============================================================================
int checkCopyAndAssignment(Epetra_Comm& Comm, bool verbose) {
  int ierr = 0;
  // check to make sure that copy ctr and op= are working correctly
  // (making level 1 deep copies)

  // create initial Map and Graph
  const int NumIndicesPerRow = 10;
  const int NumGlobalElements = 50;
  const int IndexBase = 0;
  Epetra_Map Map1(NumGlobalElements, IndexBase, Comm);
  Epetra_CrsGraph Graph1(Copy, Map1, NumIndicesPerRow);
  int g1count = Graph1.ReferenceCount();
  const Epetra_CrsGraphData* g1addr = Graph1.DataPtr();
  EPETRA_TEST_ERR(!(g1count == 1), ierr);
  if(verbose) cout << "Graph1 created (def ctr). data addr = " << g1addr << " ref. count = " << g1count << endl;

  //test copy constructor
  {
    Epetra_CrsGraph Graph2(Graph1);
    int g2count = Graph2.ReferenceCount();
    const Epetra_CrsGraphData* g2addr = Graph2.DataPtr();
    EPETRA_TEST_ERR(!(g2count == g1count+1), ierr);
    EPETRA_TEST_ERR(!(g2addr == g1addr), ierr);
    if(verbose) cout << "Graph2 created (cpy ctr). data addr = " << g2addr << " ref. count = " << g2count << endl;
  }
  int g1newcount = Graph1.ReferenceCount();
  const Epetra_CrsGraphData* g1newaddr = Graph1.DataPtr();
  EPETRA_TEST_ERR(!(g1newcount == g1count), ierr);
  EPETRA_TEST_ERR(!(g1newaddr = g1addr), ierr);
  if(verbose) cout << "Graph2 destroyed. Graph1 data addr = " << g1newaddr << " ref. count = " << g1newcount << endl;

  //test op=
  Epetra_CrsGraph Graph3(Copy, Map1, NumIndicesPerRow+1);
  int g3count = Graph3.ReferenceCount();
  const Epetra_CrsGraphData* g3addr = Graph3.DataPtr();
  EPETRA_TEST_ERR(!(g3addr != g1addr), ierr);
  if(verbose) cout << "Graph3 created (op= before). data addr = " << g3addr << " ref. count = " << g3count << endl;
  Graph3 = Graph1;
  g3count = Graph3.ReferenceCount();
  g3addr = Graph3.DataPtr();
  EPETRA_TEST_ERR(!(g3count == g1count+1), ierr);
  EPETRA_TEST_ERR(!(g3addr == g1addr), ierr);
  if(verbose) cout << "Graph3 set equal to Graph1 (op= after). data addr = " << g3addr << " ref. count = " << g3count << endl;

  return(ierr);
}

//==============================================================================
int check(Epetra_CrsGraph& A, int NumMyRows1, int NumGlobalRows1, int NumMyNonzeros1,
    int NumGlobalNonzeros1, int* MyGlobalElements, bool verbose)
{

  int ierr = 0;
  int i;
  int j;
  int forierr = 0;
  int NumGlobalIndices;
  int NumMyIndices;
  int* MyViewIndices;
  int MaxNumIndices = A.MaxNumIndices();
  int* MyCopyIndices = new int[MaxNumIndices];
  int* GlobalCopyIndices = new int[MaxNumIndices];

  // Test query functions

  int NumMyRows = A.NumMyRows();
  if(verbose) cout << "Number of local Rows = " << NumMyRows << endl;

  EPETRA_TEST_ERR(!(NumMyRows==NumMyRows1),ierr);

  int NumMyNonzeros = A.NumMyNonzeros();
  if(verbose) cout << "Number of local Nonzero entries = " << NumMyNonzeros << endl;

  EPETRA_TEST_ERR(!(NumMyNonzeros==NumMyNonzeros1),ierr);

  int NumGlobalRows = A.NumGlobalRows();
  if(verbose) cout << "Number of global Rows = " << NumGlobalRows << endl;

  EPETRA_TEST_ERR(!(NumGlobalRows==NumGlobalRows1),ierr);

  int NumGlobalNonzeros = A.NumGlobalNonzeros();
  if(verbose) cout << "Number of global Nonzero entries = " << NumGlobalNonzeros << endl;

  EPETRA_TEST_ERR(!(NumGlobalNonzeros==NumGlobalNonzeros1),ierr);

  // GlobalRowView should be illegal (since we have local indices)

  EPETRA_TEST_ERR(!(A.ExtractGlobalRowView(A.RowMap().MaxMyGID(), NumGlobalIndices, GlobalCopyIndices)==-2),ierr);

  // Other binary tests

  EPETRA_TEST_ERR(A.NoDiagonal(),ierr);
  EPETRA_TEST_ERR(!(A.Filled()),ierr);
  EPETRA_TEST_ERR(!(A.MyGRID(A.RowMap().MaxMyGID())),ierr);
  EPETRA_TEST_ERR(!(A.MyGRID(A.RowMap().MinMyGID())),ierr);
  EPETRA_TEST_ERR(A.MyGRID(1+A.RowMap().MaxMyGID()),ierr);
  EPETRA_TEST_ERR(A.MyGRID(-1+A.RowMap().MinMyGID()),ierr);
  EPETRA_TEST_ERR(!(A.MyLRID(0)),ierr);
  EPETRA_TEST_ERR(!(A.MyLRID(NumMyRows-1)),ierr);
  EPETRA_TEST_ERR(A.MyLRID(-1),ierr);
  EPETRA_TEST_ERR(A.MyLRID(NumMyRows),ierr);

  forierr = 0;
  for(i = 0; i < NumMyRows; i++) {
    int Row = A.GRID(i);
    A.ExtractGlobalRowCopy(Row, MaxNumIndices, NumGlobalIndices, GlobalCopyIndices);
    A.ExtractMyRowView(i, NumMyIndices, MyViewIndices);
    forierr += !(NumGlobalIndices==NumMyIndices);
    for(j = 1; j < NumMyIndices; j++) EPETRA_TEST_ERR(!(MyViewIndices[j-1]<MyViewIndices[j]),ierr);
    for(j = 0; j < NumGlobalIndices; j++) {
      forierr += !(GlobalCopyIndices[j]==A.GCID(MyViewIndices[j]));
      forierr += !(A.LCID(GlobalCopyIndices[j])==MyViewIndices[j]);
    }
  }
  EPETRA_TEST_ERR(forierr,ierr);
  forierr = 0;
  for(i = 0; i < NumMyRows; i++) {
    int Row = A.GRID(i);
    A.ExtractGlobalRowCopy(Row, MaxNumIndices, NumGlobalIndices, GlobalCopyIndices);
    A.ExtractMyRowCopy(i, MaxNumIndices, NumMyIndices, MyCopyIndices);
    forierr += !(NumGlobalIndices==NumMyIndices);
    for(j = 1; j < NumMyIndices; j++)
      EPETRA_TEST_ERR(!(MyCopyIndices[j-1]<MyCopyIndices[j]),ierr);
    for(j = 0; j < NumGlobalIndices; j++) {
      forierr += !(GlobalCopyIndices[j]==A.GCID(MyCopyIndices[j]));
      forierr += !(A.LCID(GlobalCopyIndices[j])==MyCopyIndices[j]);
    }

  }
  EPETRA_TEST_ERR(forierr,ierr);

  delete[] MyCopyIndices;
  delete[] GlobalCopyIndices;

  if(verbose) cout << "Rows sorted check OK" << endl;

  return(ierr);
}