C++ Epetra_RowMatrix::NumGlobalCols方法代码示例

void show_matrix(const char *txt, const Epetra_RowMatrix &matrix, const Epetra_Comm &comm)
{
  int me = comm.MyPID();
  if (comm.NumProc() > 10){
    if (me == 0){
      std::cout << txt << std::endl;
      std::cout << "Printed matrix format only works for 10 or fewer processes" << std::endl;
    }
    return;
  }

  int numRows = matrix.NumGlobalRows();
  int numCols = matrix.NumGlobalCols();

  if ((numRows > 200) || (numCols > 500)){
    if (me == 0){
      std::cerr << txt << std::endl;
      std::cerr << "show_matrix: problem is too large to display" << std::endl;
    }
    return;
  }

  int *myA = new int [numRows * numCols];

  make_my_A(matrix, myA, comm);

  printMatrix(txt, myA, NULL, NULL, numRows, numCols, comm);

  delete [] myA;
}

static int make_my_A(const Epetra_RowMatrix &matrix, int *myA, const Epetra_Comm &comm)
{
  int me = comm.MyPID();

  const Epetra_Map &rowmap = matrix.RowMatrixRowMap();
  const Epetra_Map &colmap = matrix.RowMatrixColMap();

  int myRows = matrix.NumMyRows();
  int numRows = matrix.NumGlobalRows();
  int numCols = matrix.NumGlobalCols();
  int base = rowmap.IndexBase();
  int maxRow = matrix.MaxNumEntries();

  memset(myA, 0, sizeof(int) * numRows * numCols);

  int *myIndices = new int [maxRow];
  double *tmp = new double [maxRow];

  int rowLen = 0;

  for (int i=0; i< myRows; i++){

    int rc = matrix.ExtractMyRowCopy(i, maxRow, rowLen, tmp, myIndices);

    if (rc){
      if (me == 0){
        std::cout << "Error in make_my_A" << std::endl;
      }
       return 1;
    }

    int *row = myA + (numCols * (rowmap.GID(i) - base));

    for (int j=0; j < rowLen; j++){

      int colGID = colmap.GID(myIndices[j]);
      
      row[colGID - base] = me + 1;
    }
  }

  if (maxRow){
    delete [] myIndices;
    delete [] tmp;
  }
  return 0;
}

int compute_hypergraph_metrics(const Epetra_RowMatrix &matrix,
            Isorropia::Epetra::CostDescriber &costs,
            double &myGoalWeight,
            double &balance, double &cutn, double &cutl)  // output
{
  const Epetra_BlockMap &rmap = 
    static_cast<const Epetra_BlockMap &>(matrix.RowMatrixRowMap());

  const Epetra_BlockMap &cmap = 
    static_cast<const Epetra_BlockMap &>(matrix.RowMatrixColMap());

  return compute_hypergraph_metrics(rmap, cmap,
                                     matrix.NumGlobalCols(),
                                     costs,
                                     myGoalWeight,
                                     balance, cutn, cutl);

}

//=============================================================================
int Amesos_Dscpack::Solve()
{
  if (IsNumericFactorizationOK_ == false) 
    AMESOS_CHK_ERR(NumericFactorization());

  ResetTimer(0);
  ResetTimer(1);
  
  Epetra_RowMatrix *RowMatrixA = Problem_->GetMatrix();
  if (RowMatrixA == 0)
    AMESOS_CHK_ERR(-1);

  // MS // some checks on matrix size
  if (RowMatrixA->NumGlobalRows() != RowMatrixA->NumGlobalCols())
    AMESOS_CHK_ERR(-1);

  //  Convert vector b to a vector in the form that DSCPACK needs it
  //
  Epetra_MultiVector* vecX = Problem_->GetLHS(); 
  Epetra_MultiVector* vecB = Problem_->GetRHS(); 

  if ((vecX == 0) || (vecB == 0))
    AMESOS_CHK_ERR(-1); // something wrong with input

  int NumVectors = vecX->NumVectors(); 
  if (NumVectors != vecB->NumVectors())
    AMESOS_CHK_ERR(-2);

  double *dscmapXvalues ;
  int dscmapXlda ;
  Epetra_MultiVector dscmapX(DscRowMap(),NumVectors) ; 
  int ierr;
  AMESOS_CHK_ERR(dscmapX.ExtractView(&dscmapXvalues,&dscmapXlda));
  assert (dscmapXlda == NumLocalCols); 

  double *dscmapBvalues ;
  int dscmapBlda ;
  Epetra_MultiVector dscmapB(DscRowMap(), NumVectors ) ; 
  ierr = dscmapB.ExtractView( &dscmapBvalues, &dscmapBlda );
  AMESOS_CHK_ERR(ierr);
  assert( dscmapBlda == NumLocalCols ) ; 

  AMESOS_CHK_ERR(dscmapB.Import(*vecB, Importer(), Insert));

  VecRedistTime_ = AddTime("Total vector redistribution time", VecRedistTime_, 0);
  ResetTimer(0);
  
  // MS // now solve the problem
  
  std::vector<double> ValuesInNewOrder( NumLocalCols ) ;

  OverheadTime_ = AddTime("Total Amesos overhead time", OverheadTime_, 1);

  if ( MyDscRank >= 0 ) {
    for ( int j =0 ; j < NumVectors; j++ ) { 
      for ( int i = 0; i < NumLocalCols; i++ ) { 
	ValuesInNewOrder[i] = dscmapBvalues[DscColMap().LID( LocalStructOldNum[i] ) +j*dscmapBlda ] ;
      }
      AMESOS_CHK_ERR( DSC_InputRhsLocalVec ( PrivateDscpackData_->MyDSCObject_, &ValuesInNewOrder[0], NumLocalCols ) ) ;
      AMESOS_CHK_ERR( DSC_Solve ( PrivateDscpackData_->MyDSCObject_ ) ) ; 
      AMESOS_CHK_ERR( DSC_GetLocalSolution ( PrivateDscpackData_->MyDSCObject_, &ValuesInNewOrder[0], NumLocalCols ) ) ; 
      for ( int i = 0; i < NumLocalCols; i++ ) { 
	dscmapXvalues[DscColMap().LID( LocalStructOldNum[i] ) +j*dscmapXlda ] = ValuesInNewOrder[i];
      }
    }
  }

  SolveTime_ = AddTime("Total solve time", SolveTime_, 0);
  ResetTimer(0);
  ResetTimer(1);

  vecX->Export( dscmapX, Importer(), Insert ) ;

  VecRedistTime_ = AddTime("Total vector redistribution time", VecRedistTime_, 0);

  if (ComputeTrueResidual_)
    ComputeTrueResidual(*(GetProblem()->GetMatrix()), *vecX, *vecB, 
                        false, "Amesos_Dscpack");

  if (ComputeVectorNorms_)
    ComputeVectorNorms(*vecX, *vecB, "Amesos_Dscpack");
  OverheadTime_ = AddTime("Total Amesos overhead time", OverheadTime_, 1);
  
  NumSolve_++;

  return(0) ; 
}

//=============================================================================
int Amesos_Dscpack::PerformNumericFactorization()
{
  ResetTimer(0);
  ResetTimer(1);

  Epetra_RowMatrix* RowMatrixA = Problem_->GetMatrix();
  if (RowMatrixA == 0)
    AMESOS_CHK_ERR(-1);

  const Epetra_Map& OriginalMap = RowMatrixA->RowMatrixRowMap() ; 

  int numrows = RowMatrixA->NumGlobalRows();
  assert( numrows == RowMatrixA->NumGlobalCols() );
  
  //
  //  Call Dscpack to perform Numeric Factorization
  //  
  std::vector<double> MyANonZ;
#if 0
    if ( IsNumericFactorizationOK_ ) { 
      DSC_ReFactorInitialize(PrivateDscpackData_->MyDSCObject);
    }
#endif

  DscRowMap_ = Teuchos::rcp(new Epetra_Map(numrows, NumLocalCols, 
                                           LocalStructOldNum, 0, Comm()));

  if (DscRowMap_.get() == 0)
    AMESOS_CHK_ERR(-1);
  
  Importer_ = rcp(new Epetra_Import(DscRowMap(), OriginalMap));
    
  //
  //  Import from the CrsMatrix
  //
  Epetra_CrsMatrix DscMat(Copy, DscRowMap(), 0);
  AMESOS_CHK_ERR(DscMat.Import(*RowMatrixA, Importer(), Insert));
  AMESOS_CHK_ERR(DscMat.FillComplete()); 

  DscColMap_ = Teuchos::rcp(new Epetra_Map(DscMat.RowMatrixColMap()));

  assert( MyDscRank >= 0 || NumLocalNonz == 0 ) ;
  assert( MyDscRank >= 0 || NumLocalCols == 0 ) ;
  assert( MyDscRank >= 0 || NumGlobalCols == 0  ) ; 
  MyANonZ.resize( NumLocalNonz ) ; 
  int NonZIndex = 0 ;

  int max_num_entries = DscMat.MaxNumEntries() ; 
  std::vector<int> col_indices( max_num_entries ) ; 
  std::vector<double> mat_values( max_num_entries ) ; 
  assert( NumLocalCols == DscRowMap().NumMyElements() ) ;
  std::vector<int> my_global_elements( NumLocalCols ) ; 
  AMESOS_CHK_ERR(DscRowMap().MyGlobalElements( &my_global_elements[0] ) ) ;
  
  std::vector<int> GlobalStructOldColNum( NumGlobalCols ) ; 
  
  typedef std::pair<int, double> Data; 
  std::vector<Data> sort_array(max_num_entries); 
  std::vector<int>  sort_indices(max_num_entries);
  
  for ( int i = 0; i < NumLocalCols ; i++ ) { 
    assert( my_global_elements[i] == LocalStructOldNum[i] ) ; 
    int num_entries_this_row; 
#ifdef USE_LOCAL
    AMESOS_CHK_ERR( DscMat.ExtractMyRowCopy( i, max_num_entries, num_entries_this_row, 
					     &mat_values[0], &col_indices[0] ) ) ; 
#else
    AMESOS_CHK_ERR( DscMat.ExtractGlobalRowCopy( DscMat.GRID(i), max_num_entries, num_entries_this_row, 
						 &mat_values[0], &col_indices[0] ) ) ; 
#endif
    int OldRowNumber =  LocalStructOldNum[i] ;
    if (GlobalStructOwner[ OldRowNumber ] == -1)
      AMESOS_CHK_ERR(-1);
    
    int NewRowNumber = GlobalStructNewColNum[ my_global_elements[ i ] ] ; 
    
    //
    //  Sort the column elements 
    //
    
    for ( int j = 0; j < num_entries_this_row; j++ ) { 
#ifdef USE_LOCAL
      sort_array[j].first = GlobalStructNewColNum[ DscMat.GCID( col_indices[j])] ; 
      sort_indices[j] =  GlobalStructNewColNum[ DscMat.GCID( col_indices[j])] ; 
#else
      sort_array[j].first = GlobalStructNewColNum[ col_indices[j] ]; 
#endif
      sort_array[j].second = mat_values[j] ; 
    }
    sort(&sort_array[0], &sort_array[num_entries_this_row]);
    
    for ( int j = 0; j < num_entries_this_row; j++ ) { 
      int NewColNumber = sort_array[j].first ; 
      if ( NewRowNumber <= NewColNumber ) MyANonZ[ NonZIndex++ ] = sort_array[j].second ; 

#ifndef USE_LOCAL
      assert( NonZIndex <= NumLocalNonz );  // This assert can fail on non-symmetric matrices
#endif
    }
//.........这里部分代码省略.........

int check(Epetra_RowMatrix& A, Epetra_RowMatrix & B, bool verbose)  {

    int ierr = 0;
    EPETRA_TEST_ERR(!A.Comm().NumProc()==B.Comm().NumProc(),ierr);
    EPETRA_TEST_ERR(!A.Comm().MyPID()==B.Comm().MyPID(),ierr);
    EPETRA_TEST_ERR(!A.Filled()==B.Filled(),ierr);
    EPETRA_TEST_ERR(!A.HasNormInf()==B.HasNormInf(),ierr);
    EPETRA_TEST_ERR(!A.LowerTriangular()==B.LowerTriangular(),ierr);
    EPETRA_TEST_ERR(!A.Map().SameAs(B.Map()),ierr);
    EPETRA_TEST_ERR(!A.MaxNumEntries()==B.MaxNumEntries(),ierr);
    EPETRA_TEST_ERR(!A.NumGlobalCols()==B.NumGlobalCols(),ierr);
    EPETRA_TEST_ERR(!A.NumGlobalDiagonals()==B.NumGlobalDiagonals(),ierr);
    EPETRA_TEST_ERR(!A.NumGlobalNonzeros()==B.NumGlobalNonzeros(),ierr);
    EPETRA_TEST_ERR(!A.NumGlobalRows()==B.NumGlobalRows(),ierr);
    EPETRA_TEST_ERR(!A.NumMyCols()==B.NumMyCols(),ierr);
    EPETRA_TEST_ERR(!A.NumMyDiagonals()==B.NumMyDiagonals(),ierr);
    EPETRA_TEST_ERR(!A.NumMyNonzeros()==B.NumMyNonzeros(),ierr);
    for (int i=0; i<A.NumMyRows(); i++) {
        int nA, nB;
        A.NumMyRowEntries(i,nA);
        B.NumMyRowEntries(i,nB);
        EPETRA_TEST_ERR(!nA==nB,ierr);
    }
    EPETRA_TEST_ERR(!A.NumMyRows()==B.NumMyRows(),ierr);
    EPETRA_TEST_ERR(!A.OperatorDomainMap().SameAs(B.OperatorDomainMap()),ierr);
    EPETRA_TEST_ERR(!A.OperatorRangeMap().SameAs(B.OperatorRangeMap()),ierr);
    EPETRA_TEST_ERR(!A.RowMatrixColMap().SameAs(B.RowMatrixColMap()),ierr);
    EPETRA_TEST_ERR(!A.RowMatrixRowMap().SameAs(B.RowMatrixRowMap()),ierr);
    EPETRA_TEST_ERR(!A.UpperTriangular()==B.UpperTriangular(),ierr);
    EPETRA_TEST_ERR(!A.UseTranspose()==B.UseTranspose(),ierr);

    int NumVectors = 5;
    {   // No transpose case
        Epetra_MultiVector X(A.OperatorDomainMap(), NumVectors);
        Epetra_MultiVector YA1(A.OperatorRangeMap(), NumVectors);
        Epetra_MultiVector YA2(YA1);
        Epetra_MultiVector YB1(YA1);
        Epetra_MultiVector YB2(YA1);
        X.Random();

        bool transA = false;
        A.SetUseTranspose(transA);
        B.SetUseTranspose(transA);
        A.Apply(X,YA1);
        A.Multiply(transA, X, YA2);
        EPETRA_TEST_ERR(checkMultiVectors(YA1,YA2,"A Multiply and A Apply", verbose),ierr);
        B.Apply(X,YB1);
        EPETRA_TEST_ERR(checkMultiVectors(YA1,YB1,"A Multiply and B Multiply", verbose),ierr);
        B.Multiply(transA, X, YB2);
        EPETRA_TEST_ERR(checkMultiVectors(YA1,YB2,"A Multiply and B Apply", verbose), ierr);

    }
    {   // transpose case
        Epetra_MultiVector X(A.OperatorRangeMap(), NumVectors);
        Epetra_MultiVector YA1(A.OperatorDomainMap(), NumVectors);
        Epetra_MultiVector YA2(YA1);
        Epetra_MultiVector YB1(YA1);
        Epetra_MultiVector YB2(YA1);
        X.Random();

        bool transA = true;
        A.SetUseTranspose(transA);
        B.SetUseTranspose(transA);
        A.Apply(X,YA1);
        A.Multiply(transA, X, YA2);
        EPETRA_TEST_ERR(checkMultiVectors(YA1,YA2, "A Multiply and A Apply (transpose)", verbose),ierr);
        B.Apply(X,YB1);
        EPETRA_TEST_ERR(checkMultiVectors(YA1,YB1, "A Multiply and B Multiply (transpose)", verbose),ierr);
        B.Multiply(transA, X,YB2);
        EPETRA_TEST_ERR(checkMultiVectors(YA1,YB2, "A Multiply and B Apply (transpose)", verbose),ierr);

    }

    Epetra_Vector diagA(A.RowMatrixRowMap());
    EPETRA_TEST_ERR(A.ExtractDiagonalCopy(diagA),ierr);
    Epetra_Vector diagB(B.RowMatrixRowMap());
    EPETRA_TEST_ERR(B.ExtractDiagonalCopy(diagB),ierr);
    EPETRA_TEST_ERR(checkMultiVectors(diagA,diagB, "ExtractDiagonalCopy", verbose),ierr);

    Epetra_Vector rowA(A.RowMatrixRowMap());
    EPETRA_TEST_ERR(A.InvRowSums(rowA),ierr);
    Epetra_Vector rowB(B.RowMatrixRowMap());
    EPETRA_TEST_ERR(B.InvRowSums(rowB),ierr)
    EPETRA_TEST_ERR(checkMultiVectors(rowA,rowB, "InvRowSums", verbose),ierr);

    Epetra_Vector colA(A.RowMatrixColMap());
    EPETRA_TEST_ERR(A.InvColSums(colA),ierr);
    Epetra_Vector colB(B.RowMatrixColMap());
    EPETRA_TEST_ERR(B.InvColSums(colB),ierr);
    EPETRA_TEST_ERR(checkMultiVectors(colA,colB, "InvColSums", verbose),ierr);

    EPETRA_TEST_ERR(checkValues(A.NormInf(), B.NormInf(), "NormInf before scaling", verbose), ierr);
    EPETRA_TEST_ERR(checkValues(A.NormOne(), B.NormOne(), "NormOne before scaling", verbose),ierr);

    EPETRA_TEST_ERR(A.RightScale(colA),ierr);
    EPETRA_TEST_ERR(B.RightScale(colB),ierr);


    EPETRA_TEST_ERR(A.LeftScale(rowA),ierr);
    EPETRA_TEST_ERR(B.LeftScale(rowB),ierr);
//.........这里部分代码省略.........

int Amesos_Scalapack::RedistributeA( ) {

  if( debug_ == 1 ) std::cout << "Entering `RedistributeA()'" << std::endl;

  Time_->ResetStartTime();
  
  Epetra_RowMatrix *RowMatrixA = dynamic_cast<Epetra_RowMatrix *>(Problem_->GetOperator());
  EPETRA_CHK_ERR( RowMatrixA == 0 ) ; 

  const Epetra_Map &OriginalMap = RowMatrixA->RowMatrixRowMap() ; 
  int NumberOfProcesses = Comm().NumProc() ; 

  //
  //  Compute a uniform distribution as ScaLAPACK would want it
  //    MyFirstElement - The first element which this processor would have
  //    NumExpectedElemetns - The number of elements which this processor would have
  //

  int NumRows_ = RowMatrixA->NumGlobalRows() ; 
  int NumColumns_  = RowMatrixA->NumGlobalCols() ; 
  if ( MaxProcesses_ > 0 ) {
    NumberOfProcesses = EPETRA_MIN( NumberOfProcesses, MaxProcesses_ ) ; 
  }
  else {
    int ProcessNumHeuristic = (1+NumRows_/200)*(1+NumRows_/200);
    NumberOfProcesses = EPETRA_MIN( NumberOfProcesses,  ProcessNumHeuristic );
  }
  
  if ( debug_ == 1) std::cout  << "iam_ = " << iam_  << "Amesos_Scalapack.cpp:171" << std::endl;
  //
  // Create the ScaLAPACK data distribution.
  // The TwoD data distribution is created in a completely different
  // manner and is not transposed (whereas the SaLAPACK 1D data
  // distribution was transposed) 
  //
  if ( debug_ == 1) std::cout  << "iam_ = " << iam_  << "Amesos_Scalapack.cpp:163" << std::endl;
  Comm().Barrier(); 
  if ( TwoD_distribution_ ) { 
    if ( debug_ == 1) std::cout  << "iam_ = " << iam_  << "Amesos_Scalapack.cpp:166" << std::endl;
    Comm().Barrier(); 
    npcol_ = EPETRA_MIN( NumberOfProcesses, 
			 EPETRA_MAX ( 2, (int) sqrt( NumberOfProcesses * 0.5 ) ) ) ; 
    nprow_ = NumberOfProcesses / npcol_ ;

    //
    //  Create the map for FatA - our first intermediate matrix
    //
    int NumMyElements = RowMatrixA->RowMatrixRowMap().NumMyElements() ;
    std::vector<int> MyGlobalElements( NumMyElements );
    RowMatrixA->RowMatrixRowMap().MyGlobalElements( &MyGlobalElements[0] ) ;

    int NumMyColumns = RowMatrixA->RowMatrixColMap().NumMyElements() ;
    std::vector<int> MyGlobalColumns( NumMyColumns );
    RowMatrixA->RowMatrixColMap().MyGlobalElements( &MyGlobalColumns[0] ) ;

    if ( debug_ == 1) std::cout  << "iam_ = " << iam_  << "Amesos_Scalapack.cpp:194" << std::endl;

    std::vector<int> MyFatElements( NumMyElements * npcol_ ); 
    
    for( int LocalRow=0; LocalRow<NumMyElements; LocalRow++ ) {
      for (int i = 0 ; i < npcol_; i++ ){
	MyFatElements[LocalRow*npcol_+i] = MyGlobalElements[LocalRow]*npcol_+i;
      } 
    }
    
    Epetra_Map FatInMap( npcol_*NumRows_, NumMyElements*npcol_, 
			 &MyFatElements[0], 0, Comm() ); 
    
    //
    //  Create FatIn, our first intermediate matrix
    //
    Epetra_CrsMatrix FatIn( Copy, FatInMap, 0 );
    
    
    std::vector<std::vector<int> > FatColumnIndices(npcol_,std::vector<int>(1));
    std::vector<std::vector<double> > FatMatrixValues(npcol_,std::vector<double>(1));
    std::vector<int> FatRowPtrs(npcol_);  // A FatRowPtrs[i] = the number 
    // of entries in local row LocalRow*npcol_ + i 
    
    if ( debug_ == 1) std::cout  << "iam_ = " << iam_  << "Amesos_Scalapack.cpp:219" << std::endl;
    //
    mypcol_ = iam_%npcol_;
    myprow_ = (iam_/npcol_)%nprow_;
    if ( iam_ >= nprow_ * npcol_ ) {
      myprow_ = nprow_; 
      mypcol_ = npcol_; 
    }
    //  Each row is split into npcol_ rows, with each of the 
    //  new rows containing only those elements belonging to 
    //  its process column (in the ScaLAPACK 2D process grid)
    //
    int MaxNumIndices = RowMatrixA->MaxNumEntries();
    int NumIndices;
    std::vector<int> ColumnIndices(MaxNumIndices);
    std::vector<double> MatrixValues(MaxNumIndices); 
    
    if ( debug_ == 1) std::cout  << "iam_ = " << iam_  << "Amesos_Scalapack.cpp:232 NumMyElements = " 
			    << NumMyElements 
			    << std::endl;
    
//.........这里部分代码省略.........