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

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

  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::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 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);
}

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;
      }
    }
  }

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

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);

};