C++ Epetra_SerialDenseMatrix::Shape方法代码示例

// ============================================================================
void
Solve(const Epetra_RowMatrix* Matrix, const Epetra_MultiVector* LHS,
      const Epetra_MultiVector* RHS)
{
  if (Matrix->Comm().NumProc() != 1)
    throw(Exception(__FILE__, __LINE__,
                    "Solve() works only in serial"));
  if (LHS->NumVectors() != RHS->NumVectors())
    throw(Exception(__FILE__, __LINE__,
                    "number of vectors in multivectors not consistent"));

  if(Matrix->NumGlobalRows64() > std::numeric_limits<int>::max())
    throw(Exception(__FILE__, __LINE__,
                    "Matrix->NumGlobalRows64() > std::numeric_limits<int>::max()"));

  int n = static_cast<int>(Matrix->NumGlobalRows64());
  int NumVectors = LHS->NumVectors();

  Epetra_SerialDenseMatrix DenseMatrix;
  DenseMatrix.Shape(n, n);

  for (int i = 0 ; i < n ; ++i)
    for (int j = 0 ; j < n ; ++j)
      DenseMatrix(i,j) = 0.0;

  // allocate storage to extract matrix rows.
  int Length = Matrix->MaxNumEntries();
  vector<double> Values(Length);
  vector<int>    Indices(Length);

  for (int j = 0 ; j < Matrix->NumMyRows() ; ++j)
  {
    int NumEntries;
    int ierr = Matrix->ExtractMyRowCopy(j, Length, NumEntries,
                                        &Values[0], &Indices[0]);

    for (int k = 0 ; k < NumEntries ; ++k)
      DenseMatrix(j,Indices[k]) = Values[k];
  }

  Epetra_SerialDenseMatrix DenseX(n, NumVectors);
  Epetra_SerialDenseMatrix DenseB(n, NumVectors);

  for (int i = 0 ; i < n ; ++i)
    for (int j = 0 ; j < NumVectors ; ++j)
      DenseB(i,j) = (*RHS)[j][i];

  Epetra_SerialDenseSolver DenseSolver;

  DenseSolver.SetMatrix(DenseMatrix);
  DenseSolver.SetVectors(DenseX,DenseB);

  DenseSolver.Factor();
  DenseSolver.Solve();

  for (int i = 0 ; i < n ; ++i)
    for (int j = 0 ; j < NumVectors ; ++j)
       (*LHS)[j][i] = DenseX(i,j);
}

//.........这里部分代码省略.........
      Epetra_SerialDenseMatrix Epetra_B(View, B, LDB, i, NRHS);
      Epetra_SerialDenseMatrix Epetra_X(View, X, LDX, i, NRHS);

      solver.SetMatrix(*Matrix);
      solver.SetVectors(Epetra_X, Epetra_B);

      ierr = check(solver, A1, LDA1,  i, NRHS, OneNorm1, B1, LDB1,  X1, LDX1, Transpose, verbose);
      assert (ierr>-1);
      delete Matrix;
      if (ierr!=0) {
	if (verbose) cout << "Factorization failed due to bad conditioning.  This is normal if RCOND is small."
			  << endl;
	break;
      }
    }
  }

  delete [] A;
  delete [] A1;
  delete [] X;
  delete [] X1;
  delete [] B;
  delete [] B1;

  /////////////////////////////////////////////////////////////////////
  // Now test norms and scaling functions
  /////////////////////////////////////////////////////////////////////

  Epetra_SerialDenseMatrix D;
  double ScalarA = 2.0;

  int DM = 10;
  int DN = 8;
  D.Shape(DM, DN);
  for (j=0; j<DN; j++)
    for (i=0; i<DM; i++) D[j][i] = (double) (1+i+j*DM) ;

  //cout << D << endl;

  double NormInfD_ref = (double)(DM*(DN*(DN+1))/2);
  double NormOneD_ref = (double)((DM*DN*(DM*DN+1))/2 - (DM*(DN-1)*(DM*(DN-1)+1))/2 );

  double NormInfD = D.NormInf();
  double NormOneD = D.NormOne();

  if (verbose) {
    cout << " *** Before scaling *** " << endl
	 << " Computed one-norm of test matrix = " << NormOneD << endl
	 << " Expected one-norm                = " << NormOneD_ref << endl
	 << " Computed inf-norm of test matrix = " << NormInfD << endl
	 << " Expected inf-norm                = " << NormInfD_ref << endl;
  }
  D.Scale(ScalarA); // Scale entire D matrix by this value
  NormInfD = D.NormInf();
  NormOneD = D.NormOne();
  if (verbose) {
    cout << " *** After scaling *** " << endl
	 << " Computed one-norm of test matrix = " << NormOneD << endl
	 << " Expected one-norm                = " << NormOneD_ref*ScalarA << endl
	 << " Computed inf-norm of test matrix = " << NormInfD << endl
	 << " Expected inf-norm                = " << NormInfD_ref*ScalarA << endl;
  }


  /////////////////////////////////////////////////////////////////////
  // Now test that A.Multiply(false, x, y) produces the same result

// main driver
int main(int argc, char *argv[])
{

#ifdef HAVE_MPI
  MPI_Init(&argc, &argv);
  Epetra_MpiComm Comm(MPI_COMM_WORLD);
#else
  Epetra_SerialComm Comm;
#endif

  if (Comm.NumProc() != 2) {
#ifdef HAVE_MPI
    MPI_Finalize();
#endif
    return(0);
  }

  int NumMyElements = 0;         // NODES assigned to this processor
  int NumMyExternalElements = 0; // nodes used by this proc, but not hosted
  int NumMyTotalElements = 0;
  int FE_NumMyElements = 0;      // TRIANGLES assigned to this processor
  int * MyGlobalElements = 0;    // nodes assigned to this processor
  Epetra_IntSerialDenseMatrix T; // store the grid connectivity

  int MyPID=Comm.MyPID();

  cout << MyPID << endl;

  switch( MyPID ) {

  case 0:
    NumMyElements = 3;
    NumMyExternalElements = 2;
    NumMyTotalElements = NumMyElements + NumMyExternalElements;
    FE_NumMyElements = 3;

    MyGlobalElements = new int[NumMyTotalElements];
    MyGlobalElements[0] = 0;
    MyGlobalElements[1] = 4;
    MyGlobalElements[2] = 3;
    MyGlobalElements[3] = 1;
    MyGlobalElements[4] = 5;

    break;
  case 1:
    NumMyElements = 3;
    NumMyExternalElements = 2;
    NumMyTotalElements = NumMyElements + NumMyExternalElements;
    FE_NumMyElements = 3;

    MyGlobalElements = new int[NumMyTotalElements];
    MyGlobalElements[0] = 1;
    MyGlobalElements[1] = 2;
    MyGlobalElements[2] = 5;
    MyGlobalElements[3] = 0;
    MyGlobalElements[4] = 4;
    break;

  }

  // build Map corresponding to update
  Epetra_Map Map(-1,NumMyElements,MyGlobalElements,0,Comm);

  // vector containing coordinates BEFORE exchanging external nodes
  Epetra_Vector CoordX_noExt(Map);
  Epetra_Vector CoordY_noExt(Map);

  switch( MyPID ) {

  case 0:
    T.Shape(3,FE_NumMyElements);

    // fill x-coordinates
    CoordX_noExt[0] = 0.0; 
    CoordX_noExt[1] = 1.0; 
    CoordX_noExt[2] = 0.0;
    // fill y-coordinates
    CoordY_noExt[0] = 0.0; 
    CoordY_noExt[1] = 1.0; 
    CoordY_noExt[2] = 1.0;
    // fill connectivity
    T(0,0) = 0; T(0,1) = 4; T(0,2) = 3;
    T(1,0) = 0; T(1,1) = 1; T(1,2) = 4;
    T(2,0) = 4; T(2,1) = 1; T(2,2) = 5;
    break;
    
  case 1:

    T.Shape(3,FE_NumMyElements);

    // fill x-coordinates
    CoordX_noExt[0] = 1.0; 
    CoordX_noExt[1] = 2.0; 
    CoordX_noExt[2] = 2.0;
    // fill y-coordinates
    CoordY_noExt[0] = 0.0; 
    CoordY_noExt[1] = 0.0; 
    CoordY_noExt[2] = 1.0;
    // fill connectivity
//.........这里部分代码省略.........