C++ LAPACK::GETRS方法代码示例

int main(int argc, char* argv[])
{
#ifdef HAVE_MPI
  MPI_Init(&argc,&argv);
#endif

  // Creating an instance of the LAPACK class for double-precision routines looks like:
  Teuchos::LAPACK<int, double> lapack;

  // This instance provides the access to all the LAPACK routines.
  Teuchos::SerialDenseMatrix<int, double> My_Matrix(4,4);
  Teuchos::SerialDenseVector<int, double> My_Vector(4);
  My_Matrix.random();
  My_Vector.random();

  // Perform an LU factorization of this matrix. 
  int ipiv[4], info;
  char TRANS = 'N';
  lapack.GETRF( 4, 4, My_Matrix.values(), My_Matrix.stride(), ipiv, &info ); 
  
  // Solve the linear system.
  lapack.GETRS( TRANS, 4, 1, My_Matrix.values(), My_Matrix.stride(), 
		ipiv, My_Vector.values(), My_Vector.stride(), &info );  

  // Print out the solution.
  cout << My_Vector << endl;

#ifdef HAVE_MPI
  MPI_Finalize();
#endif
  return 0;
}

int MxDimMatrix<T, DIM>::solve(MxDimVector<T, DIM> & x, const MxDimVector<T, DIM> & b) const {
  x = b;

  MxDimMatrix<T, DIM> copy(*this);

  Teuchos::LAPACK<int, T> lapack;
  MxDimVector<int, DIM> ipiv;
  //int ipiv[DIM];
  int info;

  lapack.GETRF(DIM, DIM, &copy(0, 0), DIM, &ipiv[0], &info);
  if (info != 0)
    std::cout << "MxDimMatrix::solve(...): error in lapack routine getrf. Return code: " << info << ".\n";

  lapack.GETRS('T', DIM, 1, &copy(0, 0), DIM, &ipiv[0], &x[0], DIM, &info);
  if (info != 0)
    std::cout << "MxDimMatrix::solve(...): error in lapack routine getrs. Return code: " << info << ".\n";

  return info;
}

void updateGuess(Teuchos::SerialDenseVector<int, std::complex<double> >& myCurrentGuess,
		Teuchos::SerialDenseVector<int, std::complex<double> >& myTargetsCalculated,
		Teuchos::SerialDenseMatrix<int, std::complex<double> >& myJacobian, 
		Teuchos::LAPACK<int, std::complex<double> >& myLAPACK
		 )
{
	//v = J(inverse) * (-F(x))
	//new guess = v + old guess
	myTargetsCalculated *= -1.0;

	//Perform an LU factorization of this matrix. 
	int ipiv[NUMDIMENSIONS], info;
	char TRANS = 'N';
	myLAPACK.GETRF( NUMDIMENSIONS, NUMDIMENSIONS, myJacobian.values(), myJacobian.stride(), ipiv, &info ); 

	// Solve the linear system.
	myLAPACK.GETRS( TRANS, NUMDIMENSIONS, 1, myJacobian.values(), myJacobian.stride(),
		       	ipiv, myTargetsCalculated.values(), myTargetsCalculated.stride(), &info );  

	//We have overwritten myTargetsCalculated with guess update values
	//myBLAS.AXPY(NUMDIMENSIONS, 1.0, myGuessAdjustment.values(), 1, myCurrentGuess.values(), 1);
	myCurrentGuess += myTargetsCalculated;
}

NOX::Abstract::Group::ReturnType 
LOCA::BorderedSolver::LowerTriangularBlockElimination::
solve(Teuchos::ParameterList& params,
      const LOCA::BorderedSolver::AbstractOperator& op,
      const LOCA::MultiContinuation::ConstraintInterface& B,
      const NOX::Abstract::MultiVector::DenseMatrix& C,
      const NOX::Abstract::MultiVector* F,
      const NOX::Abstract::MultiVector::DenseMatrix* G,
      NOX::Abstract::MultiVector& X,
      NOX::Abstract::MultiVector::DenseMatrix& Y) const
{
  string callingFunction = 
    "LOCA::BorderedSolver::LowerTriangularBlockElimination::solve()";
  NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;
  NOX::Abstract::Group::ReturnType status;

  // Determine if X or Y is zero
  bool isZeroF = (F == NULL);
  bool isZeroG = (G == NULL);
  bool isZeroB = B.isDXZero();
  bool isZeroX = isZeroF;
  bool isZeroY = isZeroG && (isZeroB  || isZeroX);

  // First compute X
  if (isZeroX)
    X.init(0.0);
  else {
    // Solve X = J^-1 F, note F must be nonzero
    status = op.applyInverse(params, *F, X);
    finalStatus = 
      globalData->locaErrorCheck->combineAndCheckReturnTypes(status, 
							     finalStatus,
							     callingFunction);
  }

  // Now compute Y
  if (isZeroY)
    Y.putScalar(0.0);
  else {
    // Compute G - B^T*X and store in Y
    if (isZeroG) 
      B.multiplyDX(-1.0, X, Y);
    else {
      Y.assign(*G);
      if (!isZeroB && !isZeroX) {
	NOX::Abstract::MultiVector::DenseMatrix T(Y.numRows(),Y.numCols());
	B.multiplyDX(1.0, X, T);
	Y -= T;
      }
    }

    // Overwrite Y with Y = C^-1 * (G - B^T*X)
    NOX::Abstract::MultiVector::DenseMatrix M(C);
    int *ipiv = new int[M.numRows()];
    Teuchos::LAPACK<int,double> L;
    int info;
    L.GETRF(M.numRows(), M.numCols(), M.values(), M.stride(), ipiv, &info);
    if (info != 0) {
      status = NOX::Abstract::Group::Failed;
      finalStatus = 
	globalData->locaErrorCheck->combineAndCheckReturnTypes(
							      status, 
							      finalStatus,
							      callingFunction);
    }
    L.GETRS('N', M.numRows(), Y.numCols(), M.values(), M.stride(), ipiv, 
	    Y.values(), Y.stride(), &info);
    delete [] ipiv;
    if (info != 0) {
      status = NOX::Abstract::Group::Failed;
      finalStatus = 
	globalData->locaErrorCheck->combineAndCheckReturnTypes(
							     status, 
							     finalStatus,
							     callingFunction);
    }
  }

  return finalStatus;
}

void DenseContainer<MatrixType, LocalScalarType>::
applyImpl (const local_mv_type& X,
           local_mv_type& Y,
           Teuchos::ETransp mode,
           LocalScalarType alpha,
           LocalScalarType beta) const
{
    using Teuchos::ArrayRCP;
    using Teuchos::RCP;
    using Teuchos::rcp;
    using Teuchos::rcpFromRef;

    TEUCHOS_TEST_FOR_EXCEPTION(
        X.getLocalLength () != Y.getLocalLength (),
        std::logic_error, "Ifpack2::DenseContainer::applyImpl: X and Y have "
        "incompatible dimensions (" << X.getLocalLength () << " resp. "
        << Y.getLocalLength () << ").  Please report this bug to "
        "the Ifpack2 developers.");
    TEUCHOS_TEST_FOR_EXCEPTION(
        localMap_->getNodeNumElements () != X.getLocalLength (),
        std::logic_error, "Ifpack2::DenseContainer::applyImpl: The inverse "
        "operator and X have incompatible dimensions (" <<
        localMap_->getNodeNumElements () << " resp. "
        << X.getLocalLength () << ").  Please report this bug to "
        "the Ifpack2 developers.");
    TEUCHOS_TEST_FOR_EXCEPTION(
        localMap_->getNodeNumElements () != Y.getLocalLength (),
        std::logic_error, "Ifpack2::DenseContainer::applyImpl: The inverse "
        "operator and Y have incompatible dimensions (" <<
        localMap_->getNodeNumElements () << " resp. "
        << Y.getLocalLength () << ").  Please report this bug to "
        "the Ifpack2 developers.");
    TEUCHOS_TEST_FOR_EXCEPTION(
        X.getLocalLength () != static_cast<size_t> (diagBlock_.numRows ()),
        std::logic_error, "Ifpack2::DenseContainer::applyImpl: The input "
        "multivector X has incompatible dimensions from those of the "
        "inverse operator (" << X.getLocalLength () << " vs. "
        << (mode == Teuchos::NO_TRANS ? diagBlock_.numCols () : diagBlock_.numRows ())
        << ").  Please report this bug to the Ifpack2 developers.");
    TEUCHOS_TEST_FOR_EXCEPTION(
        X.getLocalLength () != static_cast<size_t> (diagBlock_.numRows ()),
        std::logic_error, "Ifpack2::DenseContainer::applyImpl: The output "
        "multivector Y has incompatible dimensions from those of the "
        "inverse operator (" << Y.getLocalLength () << " vs. "
        << (mode == Teuchos::NO_TRANS ? diagBlock_.numRows () : diagBlock_.numCols ())
        << ").  Please report this bug to the Ifpack2 developers.");

    typedef Teuchos::ScalarTraits<local_scalar_type> STS;
    const int numVecs = static_cast<int> (X.getNumVectors ());
    if (alpha == STS::zero ()) { // don't need to solve the linear system
        if (beta == STS::zero ()) {
            // Use BLAS AXPY semantics for beta == 0: overwrite, clobbering
            // any Inf or NaN values in Y (rather than multiplying them by
            // zero, resulting in NaN values).
            Y.putScalar (STS::zero ());
        }
        else { // beta != 0
            Y.scale (beta);
        }
    }
    else { // alpha != 0; must solve the linear system
        Teuchos::LAPACK<int, local_scalar_type> lapack;
        // If beta is nonzero or Y is not constant stride, we have to use
        // a temporary output multivector.  It gets a (deep) copy of X,
        // since GETRS overwrites its (multi)vector input with its output.
        RCP<local_mv_type> Y_tmp;
        if (beta == STS::zero () ) {
            Tpetra::deep_copy (Y, X);
            Y_tmp = rcpFromRef (Y);
        }
        else {
            Y_tmp = rcp (new local_mv_type (X, Teuchos::Copy));
        }
        const int Y_stride = static_cast<int> (Y_tmp->getStride ());
        ArrayRCP<local_scalar_type> Y_view = Y_tmp->get1dViewNonConst ();
        local_scalar_type* const Y_ptr = Y_view.getRawPtr ();

        int INFO = 0;
        const char trans =
            (mode == Teuchos::CONJ_TRANS ? 'C' : (mode == Teuchos::TRANS ? 'T' : 'N'));
        lapack.GETRS (trans, diagBlock_.numRows (), numVecs,
                      diagBlock_.values (), diagBlock_.stride (),
                      ipiv_.getRawPtr (), Y_ptr, Y_stride, &INFO);
        TEUCHOS_TEST_FOR_EXCEPTION(
            INFO != 0, std::runtime_error, "Ifpack2::DenseContainer::applyImpl: "
            "LAPACK's _GETRS (solve using LU factorization with partial pivoting) "
            "failed with INFO = " << INFO << " != 0.");

        if (beta != STS::zero ()) {
            Y.update (alpha, *Y_tmp, beta);
        }
        else if (! Y.isConstantStride ()) {
            Tpetra::deep_copy (Y, *Y_tmp);
        }
    }
}

  void RCGIter<ScalarType,MV,OP>::iterate()
  {
    TEUCHOS_TEST_FOR_EXCEPTION( initialized_ == false, RCGIterFailure,
                        "Belos::RCGIter::iterate(): RCGIter class not initialized." );
    
    // We'll need LAPACK
    Teuchos::LAPACK<int,ScalarType> lapack;

    // Create convenience variables for zero and one.
    ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
    ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();

    // Allocate memory for scalars
    std::vector<int> index(1);
    Teuchos::SerialDenseMatrix<int,ScalarType> pAp(1,1), rTz(1,1);

    // Get the current solution std::vector.
    Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
 
    // Check that the current solution std::vector only has one column.
    TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1, RCGIterFailure,
                        "Belos::RCGIter::iterate(): current linear system has more than one std::vector!" );
    
    // Compute the current search dimension. 
    int searchDim = numBlocks_+1;

    // index of iteration within current cycle
    int i_ = 0;

    ////////////////////////////////////////////////////////////////
    // iterate until the status test tells us to stop.
    //
    // also break if our basis is full
    //
    Teuchos::RCP<const MV> p_ = Teuchos::null;
    Teuchos::RCP<MV> pnext_ = Teuchos::null;
    while (stest_->checkStatus(this) != Passed && curDim_+1 <= searchDim) {

      // Ap = A*p;
      index.resize( 1 );
      index[0] = i_;
      p_  = MVT::CloneView( *P_,  index );
      lp_->applyOp( *p_, *Ap_ );

      // d = p'*Ap;
      MVT::MvTransMv( one, *p_, *Ap_, pAp );
      (*D_)(i_,0) = pAp(0,0);

      // alpha = rTz_old / pAp
      (*Alpha_)(i_,0) = (*rTz_old_)(0,0) / pAp(0,0);

      // Check that alpha is a positive number
      TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(pAp(0,0)) <= zero, RCGIterFailure, "Belos::RCGIter::iterate(): non-positive value for p^H*A*p encountered!" );

      // x = x + (alpha * p);
      MVT::MvAddMv( one, *cur_soln_vec, (*Alpha_)(i_,0), *p_, *cur_soln_vec );
      lp_->updateSolution();

      // r = r - (alpha * Ap);
      MVT::MvAddMv( one, *r_, -(*Alpha_)(i_,0), *Ap_, *r_ );

      std::vector<MagnitudeType> norm(1);
      MVT::MvNorm( *r_, norm );
//printf("i = %i\tnorm(r) = %e\n",i_,norm[0]);

      // z = M\r
      if ( lp_->getLeftPrec() != Teuchos::null ) {
        lp_->applyLeftPrec( *r_, *z_ );
      } 
      else if ( lp_->getRightPrec() != Teuchos::null ) {
        lp_->applyRightPrec( *r_, *z_ );
      }
      else {
        z_ = r_;
      }

      // rTz_new = r'*z;
      MVT::MvTransMv( one, *r_, *z_, rTz );

      // beta = rTz_new/rTz_old;
      (*Beta_)(i_,0) = rTz(0,0) / (*rTz_old_)(0,0);

      // rTz_old = rTz_new;
      (*rTz_old_)(0,0) = rTz(0,0);

      // get pointer for next p
      index.resize( 1 );
      index[0] = i_+1;
      pnext_ = MVT::CloneViewNonConst( *P_,  index );

      if (existU_) {
        // mu = UTAU \ (AU'*z);
        Teuchos::SerialDenseMatrix<int,ScalarType> mu( Teuchos::View, *Delta_, recycleBlocks_, 1, 0, i_ );
        MVT::MvTransMv( one, *AU_, *z_, mu );
        char TRANS = 'N';
        int info;
        lapack.GETRS( TRANS, recycleBlocks_, 1, LUUTAU_->values(), LUUTAU_->stride(), &(*ipiv_)[0], mu.values(), mu.stride(), &info );
        TEUCHOS_TEST_FOR_EXCEPTION(info != 0, RCGIterLAPACKFailure,
                           "Belos::RCGIter::solve(): LAPACK GETRS failed to compute a solution.");
        // p = -(U*mu) + (beta*p) + z (in two steps)
//.........这里部分代码省略.........