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

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

      //std::cout << fe_matrix_bak << std::endl;

      for (int x_order=0;x_order<basis_order;x_order++) {
        for (int y_order=0;y_order<basis_order-x_order;y_order++) {
          for (int z_order=0;z_order<basis_order-x_order-y_order;z_order++) {
	    for (int comp=0;comp<cellDim;comp++) {
	      fe_matrix.initialize();
	      // copy mass matrix 
	      for (int i=0;i<numFields;i++) {
		for (int j=0;j<numFields;j++) {
		  fe_matrix(0,i,j) = fe_matrix_bak(0,i,j);
		}
	      }
	      
	      // clear old vector data
	      rhs_and_soln_vec.initialize();
	      
	      // now get rhs vector
	      
	      cub_points_cell.resize(1,numCubPointsCell,cellDim);
	     
	      rhs_at_cub_points_cell.initialize();
	      rhsFunc(rhs_at_cub_points_cell,
		      cub_points_cell,
		      comp, 
		      x_order,
		      y_order,
		      z_order);
            
	      cub_points_cell.resize(numCubPointsCell,cellDim);

	      cub_weights_cell.resize(numCubPointsCell);

	      FunctionSpaceTools::integrate<double>(rhs_and_soln_vec,
						    rhs_at_cub_points_cell,
						    w_value_of_basis_at_cub_points_cell,
						    COMP_BLAS);
	      
	      // solve linear system

// 	      solver.GESV(numFields, 1, &fe_matrix[0], numFields, &ipiv(0), &rhs_and_soln_vec[0], 
// 			  numFields, &info);
	      solver.POTRF('L',numFields,&fe_matrix[0],numFields,&info);
	      solver.POTRS('L',numFields,1,&fe_matrix[0],numFields,&rhs_and_soln_vec[0],numFields,&info);
	      
	      interp_points_ref.resize(1,numInterpPoints,cellDim);
	      // get exact solution for comparison
	      FieldContainer<double> exact_solution(1,numInterpPoints,cellDim);
	      exact_solution.initialize();
	      u_exact( exact_solution , interp_points_ref , comp , x_order, y_order, z_order);
	      interp_points_ref.resize(numInterpPoints,cellDim);

	      // compute interpolant
	      // first evaluate basis at interpolation points
	      value_of_basis_at_interp_points.resize(1,numFields,numInterpPoints,cellDim);
	      FunctionSpaceTools::evaluate<double>( interpolant , 
						    rhs_and_soln_vec ,
						    value_of_basis_at_interp_points );
	      value_of_basis_at_interp_points.resize(numFields,numInterpPoints,cellDim);
	      
	      RealSpaceTools<double>::subtract(interpolant,exact_solution);
	      
	      double nrm= RealSpaceTools<double>::vectorNorm(&interpolant[0],interpolant.dimension(1), NORM_TWO);
	      
	      *outStream << "\nNorm-2 error between scalar components of exact solution of order ("
			 << x_order << ", " << y_order << ", " << z_order
			 << ") in component " << comp 
			 << " and finite element interpolant of order " << basis_order << ": "
			 << nrm << "\n";
	      
	      if (nrm > zero) {
		*outStream << "\n\nPatch test failed for solution polynomial order ("
			   << x_order << ", " << y_order << ", " << z_order << ") and basis order (scalar, vector)  ("
			   << basis_order << ", " << basis_order+1 << ")\n\n";
		errorFlag++;
	      }
            }
          }
        }
      }
    }
    
  }
  
  catch (std::logic_error err) {
    *outStream << err.what() << "\n\n";
    errorFlag = -1000;
  };
  
  if (errorFlag != 0)
    std::cout << "End Result: TEST FAILED\n";
  else
    std::cout << "End Result: TEST PASSED\n";
  
  // reset format state of std::cout
  std::cout.copyfmt(oldFormatState);
  
  return errorFlag;
}

  void BlockCGIter<ScalarType,MV,OP>::iterate()
  {
    //
    // Allocate/initialize data structures
    //
    if (initialized_ == false) {
      initialize();
    }
    // Allocate data needed for LAPACK work.
    int info = 0;
    char UPLO = 'U';
    Teuchos::LAPACK<int,ScalarType> lapack;

    // Allocate memory for scalars.
    Teuchos::SerialDenseMatrix<int,ScalarType> alpha( blockSize_, blockSize_ );
    Teuchos::SerialDenseMatrix<int,ScalarType> beta( blockSize_, blockSize_ );
    Teuchos::SerialDenseMatrix<int,ScalarType> rHz( blockSize_, blockSize_ ), 
      rHz_old( blockSize_, blockSize_ ), pAp( blockSize_, blockSize_ );

    // Create convenience variables for zero and one.
    const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
    
    // Get the current solution std::vector.
    Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();

    // Check that the current solution std::vector has blockSize_ columns. 
    TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != blockSize_, CGIterateFailure,
                        "Belos::BlockCGIter::iterate(): current linear system does not have the right number of vectors!" );
    int rank = ortho_->normalize( *P_, Teuchos::null );
    TEUCHOS_TEST_FOR_EXCEPTION(rank != blockSize_,CGIterationOrthoFailure,
                         "Belos::BlockCGIter::iterate(): Failed to compute initial block of orthonormal direction vectors.");


    ////////////////////////////////////////////////////////////////
    // Iterate until the status test tells us to stop.
    //
    while (stest_->checkStatus(this) != Passed) {
        
      // Increment the iteration
      iter_++;
    
      // Multiply the current direction std::vector by A and store in Ap_
      lp_->applyOp( *P_, *AP_ );
      
      // Compute alpha := <P_,R_> / <P_,AP_>
      // 1) Compute P^T * A * P = pAp and P^T * R 
      // 2) Compute the Cholesky Factorization of pAp
      // 3) Back and forward solves to compute alpha
      //
      MVT::MvTransMv( one, *P_, *R_, alpha );
      MVT::MvTransMv( one, *P_, *AP_, pAp );      
     
      // Compute Cholesky factorization of pAp
      lapack.POTRF(UPLO, blockSize_, pAp.values(), blockSize_, &info);
      TEUCHOS_TEST_FOR_EXCEPTION(info != 0,CGIterationLAPACKFailure,
                         "Belos::BlockCGIter::iterate(): Failed to compute Cholesky factorization using LAPACK routine POTRF.");

      // Compute alpha by performing a back and forward solve with the Cholesky factorization in pAp.
      lapack.POTRS(UPLO, blockSize_, blockSize_, pAp.values(), blockSize_, 
		   alpha.values(), blockSize_, &info);
      TEUCHOS_TEST_FOR_EXCEPTION(info != 0,CGIterationLAPACKFailure,
                         "Belos::BlockCGIter::iterate(): Failed to compute alpha using Cholesky factorization (POTRS).");
      
      //
      // Update the solution std::vector X := X + alpha * P_
      //
      MVT::MvTimesMatAddMv( one, *P_, alpha, one, *cur_soln_vec );
      lp_->updateSolution();
      //
      // Compute the new residual R_ := R_ - alpha * AP_
      //
      MVT::MvTimesMatAddMv( -one, *AP_, alpha, one, *R_ );
      //
      // Compute the new preconditioned residual, Z_.
      if ( lp_->getLeftPrec() != Teuchos::null ) {
        lp_->applyLeftPrec( *R_, *Z_ );
        if ( lp_->getRightPrec() != Teuchos::null ) {
          Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, blockSize_ );
          lp_->applyRightPrec( *Z_, *tmp );
          Z_ = tmp;
        }
      }
      else if ( lp_->getRightPrec() != Teuchos::null ) {
        lp_->applyRightPrec( *R_, *Z_ );
      } 
      else {
        Z_ = R_;
      }
      //
      // Compute beta := <AP_,Z_> / <P_,AP_> 
      // 1) Compute AP_^T * Z_ 
      // 2) Compute the Cholesky Factorization of pAp (already have)
      // 3) Back and forward solves to compute beta

      // Compute <AP_,Z>
      MVT::MvTransMv( -one, *AP_, *Z_, beta );
      //
      lapack.POTRS(UPLO, blockSize_, blockSize_, pAp.values(), blockSize_, 
		   beta.values(), blockSize_, &info);
      TEUCHOS_TEST_FOR_EXCEPTION(info != 0,CGIterationLAPACKFailure,
//.........这里部分代码省略.........