C++ Epetra_SerialDenseMatrix类代码示例

// ============================================================================
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);
}

int Epetra_FECrsMatrix::ReplaceGlobalValues(const Epetra_LongLongSerialDenseVector& indices,
              const Epetra_SerialDenseMatrix& values,
              int format)
{
  if (indices.Length() != values.M() || indices.Length() != values.N()) {
    return(-1);
  }

  return( ReplaceGlobalValues(indices.Length(), indices.Values(),
            values.A(), format) );
}

  void ISVDUDV::shrink(const int down, std::vector<double> &S, Epetra_SerialDenseMatrix &U, Epetra_SerialDenseMatrix &V) 
  {
    //
    // put RU1 into an Epetra MultiVector
    Epetra_LocalMap LocalMap(curRank_, 0, U_->Map().Comm());
    Epetra_MultiVector Uh1(Epetra_DataAccess::View, LocalMap, U.A(), U.LDA(), curRank_-down);
    Epetra_MultiVector Vh1(Epetra_DataAccess::View, LocalMap, V.A(), V.LDA(), curRank_-down);

    //
    // update bases
    Teuchos::RCP<Epetra_MultiVector> newwU, fullU, newU, newwV, fullV, newV;
    fullU = Teuchos::rcp( new Epetra_MultiVector(Epetra_DataAccess::View,*U_,0,curRank_) );
    newwU = Teuchos::rcp( new Epetra_MultiVector(Epetra_DataAccess::View,*workU_,0,curRank_-down) );
    // multiply by Uh1
    int info = newwU->Multiply('N','N',1.0,*fullU,Uh1,0.0);
    TEUCHOS_TEST_FOR_EXCEPTION(info != 0,std::logic_error,"ISVDUDV::shrink(): Error calling EMV::Multiply(U).");
    fullU = Teuchos::null;
    newU = Teuchos::rcp( new Epetra_MultiVector(Epetra_DataAccess::View,*U_,0,curRank_-down) );
    *newU = *newwU;
    newU = Teuchos::null;
    newwU = Teuchos::null;

    // multiply by Vh1
    // get multivector views of V(1:numProc,1:curRank) and workV(1:numProc,1:curRank-down)
    double *V_A, *workV_A;
    int V_LDA, workV_LDA;
    info = V_->ExtractView(&V_A,&V_LDA);
    TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
        "RBGen::ISVDUDV::shrink(): Error calling Epetra_MultiVector::ExtractView() on V_.");
    info = workV_->ExtractView(&workV_A,&workV_LDA);
    TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
        "RBGen::ISVDUDV::shrink(): Error calling Epetra_MultiVector::ExtractView() on workV_.");
    Epetra_LocalMap lclmap(numProc_,0,A_->Comm());
    fullV = Teuchos::rcp( new Epetra_MultiVector(Epetra_DataAccess::View,lclmap,    V_A,    V_LDA,curRank_     ) );
    newwV = Teuchos::rcp( new Epetra_MultiVector(Epetra_DataAccess::View,lclmap,workV_A,workV_LDA,curRank_-down) );
    // multiply workV = fullV * Vh1
    info = newwV->Multiply('N','N',1.0,*fullV,Vh1,0.0);
    TEUCHOS_TEST_FOR_EXCEPTION(info != 0,std::logic_error,"ISVDUDV::shrink(): Error calling EMV::Multiply(V).");
    fullV = Teuchos::null;
    // now set newV = workV
    newV = Teuchos::rcp( new Epetra_MultiVector(Epetra_DataAccess::View,lclmap, V_A, V_LDA, curRank_-down) );
    *newV = *newwV;
    newV = Teuchos::null;
    newwV = Teuchos::null;

    // save new singular values
    for (int i=0; i<curRank_-down; i++) {
      sigma_[i] = S[i];
    }

    curRank_ = curRank_-down;
  }

//=============================================================================
int Epetra_SerialDenseSVD::SetMatrix(Epetra_SerialDenseMatrix & A_in) {
  ResetMatrix();
  Matrix_ = &A_in;
//  Factor_ = &A_in;
  M_ = A_in.M();
  N_ = A_in.N();
  Min_MN_ = EPETRA_MIN(M_,N_);
  LDA_ = A_in.LDA();
//  LDAF_ = LDA_;
  A_ = A_in.A();
//  AF_ = A_in.A();
  return(0);
}

//=========================================================================
// checks if two matrices are independent or not
bool seperateData(Epetra_SerialDenseMatrix& a, Epetra_SerialDenseMatrix& b) {
	bool seperate;

	int r = EPETRA_MIN(a.M(),b.M()) / 2; // ensures (r,c) is valid
	int c = EPETRA_MIN(a.N(),b.N()) / 2; // in both matrices

	double orig_a = a(r,c);
	double new_value = a(r,c) + 1;
	if(b(r,c) == new_value) // there's a chance b could be independent, but
		new_value++;          // already have new_value in (r,c).
	
	a(r,c) = new_value;
	if(b(r,c) == new_value)
		seperate = false;
	else
		seperate = true;

	a(r,c) = orig_a; // undo change we made to a

	return(seperate);
}

//=========================================================================
// checks the signatures of two matrices
bool identicalSignatures(Epetra_SerialDenseMatrix& a, Epetra_SerialDenseMatrix& b, bool testLDA) {

	if((a.M()  != b.M()  )|| // check properties first
		 (a.N()  != b.N()  )||
		 (a.CV() != b.CV() ))
		return(false);

	if(testLDA == true)      // if we are coming from op= c->c #2 (have enough space)
		if(a.LDA() != b.LDA()) // then we don't check LDA (but we do check it in the test function)
			return(false);

	if(a.CV() == View) { // if we're still here, we need to check the data
		if(a.A() != b.A()) // for a view, this just means checking the pointers
			return(false);   // for a copy, this means checking each element
	}
	else { // CV == Copy
		const int m = a.M();
		const int n = a.N();
		for(int i = 0; i < m; i++)
			for(int j = 0; j < n; j++) {
				if(a(i,j) != b(i,j))
					return(false);
			}
	}

	return(true); // if we're still here, signatures are identical
}

int main(int argc, char *argv[])
{
  int ierr = 0, i, j, k;
  bool debug = false;

#ifdef EPETRA_MPI
  MPI_Init(&argc,&argv);
  Epetra_MpiComm Comm( MPI_COMM_WORLD );
#else
  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;

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

  int rank = Comm.MyPID();
  //  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
  if (verbose) cout << Comm <<endl;

  //  bool verbose1 = verbose;

  // Redefine verbose to only print on PE 0
  if (verbose && rank!=0) verbose = false;
	
  int N = 20;
  int NRHS = 4;
  double * A = new double[N*N];
  double * A1 = new double[N*N];
  double * X = new double[(N+1)*NRHS];
  double * X1 = new double[(N+1)*NRHS];
  int LDX = N+1;
  int LDX1 = N+1;
  double * B = new double[N*NRHS];
  double * B1 = new double[N*NRHS];
  int LDB = N;
  int LDB1 = N;

  int LDA = N;
  int LDA1 = LDA;
  double OneNorm1;
  bool Transpose = false;

  Epetra_SerialDenseSolver solver;
  Epetra_SerialDenseMatrix * Matrix;
  for (int kk=0; kk<2; kk++) {
    for (i=1; i<=N; i++) {
      GenerateHilbert(A, LDA, i);
      OneNorm1 = 0.0;
      for (j=1; j<=i; j++) OneNorm1 += 1.0/((double) j); // 1-Norm = 1 + 1/2 + ...+1/n

      if (kk==0) {
	Matrix = new Epetra_SerialDenseMatrix(View, A, LDA, i, i);
	LDA1 = LDA;
      }
      else {
	Matrix = new Epetra_SerialDenseMatrix(Copy, A, LDA, i, i);
	LDA1 = i;
      }

      GenerateHilbert(A1, LDA1, i);
	
      if (kk==1) {
	solver.FactorWithEquilibration(true);
	solver.SolveWithTranspose(true);
	Transpose = true;
	solver.SolveToRefinedSolution(true);
      }

      for (k=0; k<NRHS; k++)
	for (j=0; j<i; j++) {
	  B[j+k*LDB] = 1.0/((double) (k+3)*(j+3));
	  B1[j+k*LDB1] = B[j+k*LDB1];
	}
      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;
      }
    }
  }

//.........这里部分代码省略.........

void BilinearFormUtility::computeStiffnessMatrix(FieldContainer<double> &stiffness, 
                                                 FieldContainer<double> &innerProductMatrix,
                                                 FieldContainer<double> &optimalTestWeights) {
  // stiffness has dimensions (numCells, numTrialDofs, numTrialDofs)
  // innerProductMatrix has dim. (numCells, numTestDofs, numTestDofs)
  // optimalTestWeights has dim. (numCells, numTrialDofs, numTestDofs)
  // all this does is computes stiffness = weights^T * innerProductMatrix * weights
  int numCells = stiffness.dimension(0);
  int numTrialDofs = stiffness.dimension(1);
  int numTestDofs = innerProductMatrix.dimension(1);
  
  // check that all the dimensions are compatible:
  TEUCHOS_TEST_FOR_EXCEPTION( ( optimalTestWeights.dimension(0) != numCells ),
                     std::invalid_argument,
                     "stiffness.dimension(0) and optimalTestWeights.dimension(0) (numCells) do not match.");
  TEUCHOS_TEST_FOR_EXCEPTION( ( optimalTestWeights.dimension(1) != numTrialDofs ),
                     std::invalid_argument,
                     "numTrialDofs and optimalTestWeights.dimension(1) do not match.");
  TEUCHOS_TEST_FOR_EXCEPTION( ( optimalTestWeights.dimension(2) != numTestDofs ),
                     std::invalid_argument,
                     "numTestDofs and optimalTestWeights.dimension(2) do not match.");
  TEUCHOS_TEST_FOR_EXCEPTION( ( innerProductMatrix.dimension(2) != innerProductMatrix.dimension(1) ),
                     std::invalid_argument,
                     "innerProductMatrix.dimension(1) and innerProductMatrix.dimension(2) do not match.");
  
  TEUCHOS_TEST_FOR_EXCEPTION( ( stiffness.dimension(1) != stiffness.dimension(2) ),
                     std::invalid_argument,
                     "stiffness.dimension(1) and stiffness.dimension(2) do not match.");
  
  stiffness.initialize(0);
  
  for (int cellIndex=0; cellIndex < numCells; cellIndex++) {
    Epetra_SerialDenseMatrix weightsT(Copy,
                                     &optimalTestWeights(cellIndex,0,0),
                                     optimalTestWeights.dimension(2), // stride
                                     optimalTestWeights.dimension(2),optimalTestWeights.dimension(1));
    
    Epetra_SerialDenseMatrix ipMatrixT(Copy,
                                      &innerProductMatrix(cellIndex,0,0),
                                      innerProductMatrix.dimension(2), // stride
                                      innerProductMatrix.dimension(2),innerProductMatrix.dimension(1));
    
    Epetra_SerialDenseMatrix   stiffT (View,
                                      &stiffness(cellIndex,0,0),
                                      stiffness.dimension(2), // stride
                                      stiffness.dimension(2),stiffness.dimension(1));
    
    Epetra_SerialDenseMatrix intermediate( numTrialDofs, numTestDofs );
    
    // account for the fact that SDM is column-major and FC is row-major: 
    //   (weightsT) * (ipMatrixT)^T * (weightsT)^T
    int success = intermediate.Multiply('T','T',1.0,weightsT,ipMatrixT,0.0);
    
    if (success != 0) {
      cout << "computeStiffnessMatrix: intermediate.Multiply() failed with error code " << success << endl;
    }
    
    success = stiffT.Multiply('N','N',1.0,intermediate,weightsT,0.0);
    // stiffT is technically the transpose of stiffness, but the construction A^T * B * A is symmetric even in general...
    
    if (success != 0) {
      cout << "computeStiffnessMatrix: stiffT.Multiply() failed with error code " << success << endl;
    }
  }
  
  if ( ! checkForZeroRowsAndColumns("stiffness",stiffness) ) {
    //cout << "stiffness: " << stiffness;
  }
  
  bool enforceNumericalSymmetry = false;
  if (enforceNumericalSymmetry) {
    for (unsigned int c=0; c < numCells; c++)
      for (unsigned int i=0; i < numTrialDofs; i++)
        for (unsigned int j=i+1; j < numTrialDofs; j++)
        {
          stiffness(c,i,j) = (stiffness(c,i,j) + stiffness(c,j,i)) / 2.0;
          stiffness(c,j,i) = stiffness(c,i,j);
        }
  }
}

//=============================================================================
int Epetra_SerialDenseSVD::SetVectors(Epetra_SerialDenseMatrix & X_in, Epetra_SerialDenseMatrix & B_in)
{
  if (B_in.M()!=X_in.M() || B_in.N() != X_in.N()) EPETRA_CHK_ERR(-1);
  if (B_in.A()==0) EPETRA_CHK_ERR(-2);
  if (B_in.LDA()<1) EPETRA_CHK_ERR(-3);
  if (X_in.A()==0) EPETRA_CHK_ERR(-4);
  if (X_in.LDA()<1) EPETRA_CHK_ERR(-5);

  ResetVectors();
  LHS_ = &X_in;
  RHS_ = &B_in;
  NRHS_ = B_in.N();

  B_ = B_in.A();
  LDB_ = B_in.LDA();
  X_ = X_in.A();
  LDX_ = X_in.LDA();
  return(0);
}

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

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

int NSIDE = 1;
int NPIX ;
int NSTOKES = 3;

NPIX = 12. * NSIDE * NSIDE +1; //total pixel size, each pixel is an element which contains 3 floats which are IQU
Epetra_BlockMap PixMap(NPIX,NSTOKES,0,Comm);

int * PixMyGlobalElements = PixMap.MyGlobalElements();

cout <<  PixMap << endl;
Epetra_FEVbrMatrix invM(Copy, PixMap, 1);

int BlockIndices[1];
BlockIndices[0] = 2;
int RowDim, NumBlockEntries;
int err;
Epetra_SerialDenseMatrix Mpp(NSTOKES, NSTOKES);
Mpp[0][0] = 1.;
cout << Mpp << endl;

Epetra_SerialDenseMatrix * Zero;

for( int i=0 ; i<PixMap.NumMyElements(); ++i ) { //loop on local pixel
    BlockIndices[0] = PixMyGlobalElements[i];
    Zero = new Epetra_SerialDenseMatrix(NSTOKES, NSTOKES);
    invM.BeginInsertGlobalValues(BlockIndices[0], 1, BlockIndices);
    err = invM.SubmitBlockEntry(Zero->A(), Zero->LDA(), NSTOKES, NSTOKES);
            if (err != 0) {
                cout << "PID:" << Comm.MyPID() << "Error in inserting init zero values in M, error code:" << err << endl;
                }
    err = invM.EndSubmitEntries();
    }

BlockIndices[0] = 2;
cout << invM << endl;

int NumHits = 2*Comm.MyPID() + 5;
for( int i=0 ; i<NumHits; ++i ) { //loop on local pointing

    invM.BeginSumIntoGlobalValues(BlockIndices[0], 1, BlockIndices);

    err = invM.SubmitBlockEntry(Mpp.A(), Mpp.LDA(), NSTOKES, NSTOKES); //FIXME check order
            if (err != 0) {
                cout << "PID:" << Comm.MyPID() << "Error in inserting values in M, error code:" << err << endl;
                }

    err = invM.EndSubmitEntries();
            if (err != 0) {
                cout << "PID:" << Comm.MyPID() << " LocalRow[i]:" << i << " Error in ending submit entries in M, error code:" << err << endl;
                }

}
invM.GlobalAssemble();

cout << invM << endl;

if (Comm.MyPID() == 0) {

Epetra_SerialDenseMatrix * blockM;
int * BlockIndicesBlock;
invM.BeginExtractMyBlockRowView(2, RowDim, NumBlockEntries, BlockIndicesBlock);
invM.ExtractEntryView(blockM);

    cout << *blockM << endl;
    cout << "*blockM[0][0]" << endl;
    cout << *blockM[0][0] << endl;
    cout << "*blockM[0][4]" << endl;
    cout << *blockM[0][4] << endl;
    cout << "*blockM[1][1]" << endl;
    cout << *blockM[1][1] << endl;
}

#ifdef HAVE_MPI
  MPI_Finalize();
#endif

return(0);

};

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

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

  bool verbose = false;

  // Check if we should print results to standard out
  verbose = true;

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

  int rank = Comm.MyPID();

  if (verbose) cout << Comm <<endl;

  // Redefine verbose to only print on PE 0
  if (verbose && rank!=0) verbose = false;
	
  int N = 5;
  int NRHS = 5;
  double * X = new double[NRHS];
  double * ed_X = new double[NRHS];

  double * B = new double[NRHS];
  double * ed_B = new double[NRHS];

  Ifpack_SerialTriDiSolver solver;
  Ifpack_SerialTriDiMatrix * Matrix;

  Epetra_SerialDenseSolver ed_solver;
  Epetra_SerialDenseMatrix * ed_Matrix;

  bool Transpose = false;
  bool Refine = false;
  solver.SolveWithTranspose(Transpose);
  solver.SolveToRefinedSolution(Refine);

  ed_solver.SolveWithTranspose(Transpose);
  ed_solver.SolveToRefinedSolution(Refine);

  Matrix = new Ifpack_SerialTriDiMatrix(5,true);
  ed_Matrix = new Epetra_SerialDenseMatrix(5,5);

  for(int i=0;i<N;++i) {
    B[i] = ed_B[i] =2;
    Matrix->D()[i]=2.0;
    if(i<(N-1)) {
      Matrix->DL()[i]=-1.0;
      if(i!=2) Matrix->DU()[i]=-1.0;
    }
  }

  Matrix->Print(std::cout);

  double * ed_a = ed_Matrix->A();
  for(int i=0;i<N;++i)
    for(int j=0;j<N;++j) {
      if(i==j) ed_a[j*N+i] = 2.0;
      else if(abs(i-j) == 1)   ed_a[j*N+i] = -1.0;
      else  ed_a[j*N + i] = 0;
      if(i==2&&j==3) ed_a[j*N+i] = 0.0;
    }


  Epetra_SerialDenseVector LHS(Copy, X, N);
  Epetra_SerialDenseVector RHS(Copy, B, N);

  Epetra_SerialDenseVector ed_LHS(Copy, ed_X, N);
  Epetra_SerialDenseVector ed_RHS(Copy, ed_B, N);

  solver.SetMatrix(*Matrix);
  solver.SetVectors(LHS, RHS);
  
  ed_solver.SetMatrix(*ed_Matrix);
  ed_solver.SetVectors(ed_LHS, ed_RHS);

  solver.Solve();  
  ed_solver.Solve();

  std::cout << " LHS vals are: "<<std::endl;
  bool TestPassed=true;
  for(int i=0;i<N;++i) { 
    std::cout << "["<<i<<"] "<< LHS(i)<<"  "<<ed_LHS(i)<<" delta "<<LHS(i)-ed_LHS(i)<<std::endl;
    if( fabs( (LHS(i)- ed_LHS(i))/(LHS(i)+ ed_LHS(i)) ) > 1.0e-12 ) {
       TestPassed = false;
       std::cout << " not equal for "<<i<<" delta "<< LHS(i)- ed_LHS(i)<<std::endl;
    }
  }

  Ifpack_SerialTriDiMatrix * tdfac = solver.FactoredMatrix();
  Epetra_SerialDenseMatrix * sdfac = ed_solver.FactoredMatrix();

  std::cout << " Tri Di Factored "<<std::endl;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	}
      }
    }
  }

  int * Indices = new int [MaxNnzRow];
  double * Values  = new double [MaxNnzRow];
  for( i=0 ; i<MaxNnzRow ; ++i ) Values[i] = 0.0;

  // now use Struct to fill build the matrix structure
  for( int Row=0 ; Row<NumMyElements ; ++Row ) {
    int Length = 0;
    for( int j=0 ; j<MaxNnzRow ; ++j ) {
      if( Struct(Row,j) == -1 ) break;
      Indices[Length] = Struct(Row,j);
      Length++;
    }
    GlobalRow = MyGlobalElements[Row];    
    A.InsertGlobalValues(GlobalRow, Length, Values, Indices);
  }

  // replace global numbering with local one in T
  for( int Element=0 ; Element<FE_NumMyElements ; ++Element ) {
    for( int i=0 ; i<3 ; ++i ) {
      int global = T(Element,i);
      int local = find(MyGlobalElements,NumMyTotalElements,
			global);
      if( global == -1 ) {
	cerr << "ERROR\n";
	return( EXIT_FAILURE );
      }
      T(Element,i) = local;
    }
  }

  // - - - - - - - - - - - - - - //
  // M A T R I X   F I L L - I N //
  // - - - - - - - - - - - - - - //

  // room for the local matrix
  Epetra_SerialDenseMatrix Ke;
  Ke.Shape(3,3);

  // now fill the matrix
  for(  int Element=0 ; Element<FE_NumMyElements ; ++Element ) {
    // variables used inside
    int GlobalRow;
    int MyRow;
    int GlobalCol;
    double x_triangle[3];
    double y_triangle[3];
    // get the spatial coordinate of each local node
    for( int i=0 ; i<3 ; ++i ) {
      MyRow = T(Element,i);
      y_triangle[i] = CoordX[MyRow]; 
      x_triangle[i] = CoordY[MyRow]; 
    }
    // compute the local matrix for Element

    compute_loc_matrix( x_triangle, y_triangle,Ke ); 

    // insert it in the global one
    // cycle over each row
    for( int i=0 ; i<3 ; ++i ) {
      // get the global and local number of this row
      MyRow = T(Element,i);
      if( MyRow < NumMyElements ) {
	for( int j=0 ; j<3 ; ++j ) {
	  // get global column number
	  GlobalRow = MyGlobalElements[MyRow];
	  GlobalCol = MyGlobalElements[T(Element,j)];
	  A.SumIntoGlobalValues(GlobalRow,1,&(Ke(i,j)),&GlobalCol);
	}
      }
    }
  }

  A.FillComplete();

  // - - - - - - - - - - - - - //
  // R H S  &  S O L U T I O N //
  // - - - - - - - - - - - - - //

  Epetra_Vector x(Map), b(Map);
  x.Random(); b.PutScalar(0.0);

  // Solution can be obtained using Aztecoo

  // free memory before leaving
  delete MyGlobalElements;
  delete Indices;
  delete Values;

#ifdef HAVE_MPI
  MPI_Finalize();
#endif

  return( EXIT_SUCCESS );

} /* main */