C++ Epetra_MultiVector::ConstantStride方法代码示例

Epetra_OskiMultiVector::Epetra_OskiMultiVector(const Epetra_MultiVector& Source) 
  : Epetra_MultiVector(Source),
  Epetra_View_(&Source), 
  Copy_Created_(false) {
    double* A;
    double** Aptr;
    int LDA;
    int* LDAptr;
    LDAptr = new int[1];
    Aptr = new double*[1];
    if(Source.ConstantStride() || (Source.NumVectors() == 1)) {
      if(Source.ExtractView(Aptr, LDAptr))
        std::cerr << "Extract view failed\n";
      else
        Oski_View_ = oski_CreateMultiVecView(*Aptr, Source.MyLength(), Source.NumVectors(), LAYOUT_COLMAJ, *LDAptr);
    }
    else {
      Copy_Created_ = true;
      LDA = Source.MyLength();
      A = new double[LDA*Source.NumVectors()];
      if(Source.ExtractCopy(A, LDA))
        std::cerr << "Extract copy failed\n";
      else
        Oski_View_ = oski_CreateMultiVecView(A, Source.MyLength(), Source.NumVectors(), LAYOUT_COLMAJ, LDA);
    }
    delete [] LDAptr;
    delete [] Aptr;
}

//=========================================================================
// NOTE: This method should be removed and replaced with calls to Epetra_Util_ExtractHbData()
int Epetra_LinearProblemRedistor::ExtractHbData(int & M, int & N, int & nz, int * & ptr,
																								int * & ind, double * & val, int & Nrhs,
																								double * & rhs, int & ldrhs,
																								double * & lhs, int & ldlhs) const {

	Epetra_CrsMatrix * RedistMatrix = dynamic_cast<Epetra_CrsMatrix *>(RedistProblem_->GetMatrix());
	
	if (RedistMatrix==0) EPETRA_CHK_ERR(-1); // This matrix is zero or not an Epetra_CrsMatrix
	if (!RedistMatrix->IndicesAreContiguous()) { // Data must be contiguous for this to work
		EPETRA_CHK_ERR(-2);
	}

	M = RedistMatrix->NumMyRows();
	N = RedistMatrix->NumMyCols();
	nz = RedistMatrix->NumMyNonzeros();
	val = (*RedistMatrix)[0];        // Dangerous, but cheap and effective way to access first element in

	const Epetra_CrsGraph & RedistGraph = RedistMatrix->Graph();
	ind = RedistGraph[0];  // list of values and indices

	Epetra_MultiVector * LHS = RedistProblem_->GetLHS();
	Epetra_MultiVector * RHS = RedistProblem_->GetRHS();
	Nrhs = RHS->NumVectors();
	if (Nrhs>1) {
		if (!RHS->ConstantStride()) {EPETRA_CHK_ERR(-3)}; // Must have strided vectors
		if (!LHS->ConstantStride()) {EPETRA_CHK_ERR(-4)}; // Must have strided vectors
	}
	ldrhs = RHS->Stride();
	rhs = (*RHS)[0]; // Dangerous but effective (again)
	ldlhs = LHS->Stride();
	lhs = (*LHS)[0];

	// Finally build ptr vector

	if (ptr_==0) {
		ptr_ = new int[M+1];
		ptr_[0] = 0;
		for (int i=0; i<M; i++) ptr_[i+1] = ptr_[i] + RedistGraph.NumMyIndices(i);
	}
	ptr = ptr_;
	
  return(0);
}

  int  BuildMatrixTests (Epetra_MultiVector & C,
		       const char TransA, const char TransB, 
		       const double alpha, 
		       Epetra_MultiVector& A, 
		       Epetra_MultiVector& B,
		       const double beta,
		       Epetra_MultiVector& C_GEMM ) {

    // For given values of TransA, TransB, alpha and beta, a (possibly
    // zero) filled Epetra_MultiVector C, and allocated
    // Epetra_MultiVectors A, B and C_GEMM this routine will generate values for 
    // Epetra_MultiVectors A, B and C_GEMM such that, if A, B and (this) are 
    // used with GEMM in this class, the results should match the results 
    // generated by this routine.

    // Test for Strided multivectors (required for GEMM ops)

    if (!A.ConstantStride()   ||
	!B.ConstantStride()      ||
	!C_GEMM.ConstantStride() ||
	!C.ConstantStride()) return(-1); // Error 

    int i, j;
    double fi, fj;  // Used for casting loop variables to floats

    // Get a view of the MultiVectors

    double *Ap      = 0;
    double *Bp      = 0;
    double *Cp      = 0;
    double *C_GEMMp = 0;

    int A_nrows = A.MyLength();
    int A_ncols = A.NumVectors();
    int B_nrows = B.MyLength();
    int B_ncols = B.NumVectors();
    int C_nrows = C.MyLength();
    int C_ncols = C.NumVectors();
    int A_Stride         = 0;
    int B_Stride         = 0;
    int C_Stride         = 0;
    int C_GEMM_Stride    = 0;

    A.ExtractView(&Ap, &A_Stride);
    B.ExtractView(&Bp, &B_Stride);
    C.ExtractView(&Cp, &C_Stride);
    C_GEMM.ExtractView(&C_GEMMp, &C_GEMM_Stride);

      // Define some useful constants


    int opA_ncols = (TransA=='N') ? A.NumVectors() : A.MyLength();
    int opB_nrows = (TransB=='N') ? B.MyLength() : B.NumVectors();
  
    int C_global_inner_dim  = (TransA=='N') ? A.NumVectors() : A.GlobalLength();
  

    bool A_is_local = (!A.DistributedGlobal());
    bool B_is_local = (!B.DistributedGlobal());
    bool C_is_local = (!C.DistributedGlobal());

    int A_IndexBase = A.Map().IndexBase();
    int B_IndexBase = B.Map().IndexBase();
  
    // Build two new maps that we can use for defining global equation indices below
    Epetra_Map * A_Map = new Epetra_Map(-1, A_nrows, A_IndexBase, A.Map().Comm());
    Epetra_Map * B_Map = new Epetra_Map(-1, B_nrows, B_IndexBase, B.Map().Comm());

    int* A_MyGlobalElements = new int[A_nrows];
    A_Map->MyGlobalElements(A_MyGlobalElements);
    int* B_MyGlobalElements = new int[B_nrows];
    B_Map->MyGlobalElements(B_MyGlobalElements);

  // Check for compatible dimensions

    if (C.MyLength()        != C_nrows     ||
	opA_ncols      != opB_nrows   ||
	C.NumVectors()    != C_ncols     ||
	C.MyLength()        != C_GEMM.MyLength()        ||
	C.NumVectors()    != C_GEMM.NumVectors()      ) {
      delete A_Map;
      delete B_Map;
      delete [] A_MyGlobalElements;
      delete [] B_MyGlobalElements;
      return(-2); // Return error
    }

    bool Case1 = ( A_is_local &&  B_is_local &&  C_is_local);  // Case 1 above
    bool Case2 = (!A_is_local && !B_is_local &&  C_is_local && TransA=='T' );// Case 2
    bool Case3 = (!A_is_local &&  B_is_local && !C_is_local && TransA=='N');// Case 3
  
    // Test for meaningful cases

    if (!(Case1 || Case2 || Case3)) {
      delete A_Map;
      delete B_Map;
      delete [] A_MyGlobalElements;
      delete [] B_MyGlobalElements;
      return(-3); // Meaningless case
    }
//.........这里部分代码省略.........