C++ MatrixT::size1方法代码示例

 inline DCMatrix inv(const MatrixT &mat, double regularizationCoeff = 0.0) {
   BOOST_ASSERT(mat.size1() == mat.size2());
   unsigned int n = mat.size1();
   DCMatrix inv = mat; // copy data, as it will be modified below
   if (regularizationCoeff != 0.0)
     inv += regularizationCoeff * ublas::identity_matrix<double>(n);
   std::vector<int> ipiv(n); // pivot vector, is first filled by trf, then used by tri to inverse matrix
   lapack::getrf(inv,ipiv); // inv and ipiv will both be modified
   lapack::getri(inv,ipiv); // afterwards, "inv" is the inverse
   return inv;
 }

 inline DCMatrix invSym(const MatrixT &mat, double regularizationCoeff = 0.0) {
   BOOST_ASSERT(mat.size1() == mat.size2());
   unsigned int n = mat.size1();
   DCMatrix inv = mat; // copy data, as it will be modified below
   if (regularizationCoeff != 0.0)
     inv += regularizationCoeff * ublas::identity_matrix<double>(n);
   std::vector<int> ipiv(n); // pivot vector, is first filled by trf, then used by tri to inverse matrix
   // TODO (9): use "po..." (po=positive definite matrix) instead if "sy..." (symmetric indefinite matrix) to make it faster
   lapack::sytrf('U',inv,ipiv); // inv and ipiv will both be modified
   lapack::sytri('U',inv,ipiv); // afterwards, "inv" is the real inverse, but only the upper elements are valid!!!
   ublas::symmetric_adaptor<DCMatrix, ublas::upper> iSym(inv);
   return iSym; // copies upper matrix to lower
 }

 ilut_precond(MatrixT const & mat, ilut_tag const & tag) : tag_(tag), LU_(mat.size1(), mat.size2())
 {
   //initialize preconditioner:
   //std::cout << "Start CPU precond" << std::endl;
   init(mat);
   //std::cout << "End CPU precond" << std::endl;
 }

 ichol0_precond(MatrixT const & mat, ichol0_tag const & tag) : tag_(tag), LLT(mat.size1(), mat.size2(), viennacl::context(viennacl::MAIN_MEMORY))
 {
     //initialize preconditioner:
     //std::cout << "Start CPU precond" << std::endl;
     init(mat);
     //std::cout << "End CPU precond" << std::endl;
 }

    typename viennacl::result_of::cpu_value_type<typename MatrixT::value_type>::type
    eig(MatrixT const& A, power_iter_tag const & tag)
    {
      typedef typename viennacl::result_of::vector_for_matrix<MatrixT>::type    VectorT;

      VectorT eigenvec(A.size1());
      return eig(A, tag, eigenvec);
    }

void run_solver_vector(MatrixT const & matrix, VectorT2 const & vector2,VectorT & result, SolverTag)
{
  std::cout << "------- Solver tag: " <<SolverTag::name()<<" ----------" << std::endl;
  result = viennacl::linalg::solve(matrix, vector2, SolverTag());

  Timer timer;
  viennacl::backend::finish();

  timer.start();
  for (int runs=0; runs<BENCHMARK_RUNS; ++runs)
  {
    result = viennacl::linalg::solve(matrix, vector2, SolverTag());
  }
  double exec_time = timer.get();
  viennacl::backend::finish();
  std::cout << "GPU: ";printOps(double(matrix.size1() * matrix.size1()),(static_cast<double>(exec_time) / static_cast<double>(BENCHMARK_RUNS)));
  std::cout << "GPU: "<< double(matrix.size1() * matrix.size1() * sizeof(NumericT)) / (static_cast<double>(exec_time) / static_cast<double>(BENCHMARK_RUNS)) / 1e9 << " GB/sec" << std::endl;
  std::cout << "Execution time: " << exec_time/BENCHMARK_RUNS << std::endl;
  std::cout << "------- Finished: " << SolverTag::name() << " ----------" << std::endl;
}

typename viennacl::result_of::cpu_value_type<typename MatrixT::value_type>::type
eig(MatrixT const& matrix, power_iter_tag const & tag)
{

    typedef typename viennacl::result_of::value_type<MatrixT>::type           ScalarType;
    typedef typename viennacl::result_of::cpu_value_type<ScalarType>::type    CPU_ScalarType;
    typedef typename viennacl::result_of::vector_for_matrix<MatrixT>::type    VectorT;

    CPU_ScalarType eigenvalue;
    vcl_size_t matrix_size = matrix.size1();
    VectorT r(matrix_size);
    VectorT r2(matrix_size);
    std::vector<CPU_ScalarType> s(matrix_size);

    for (vcl_size_t i=0; i<s.size(); ++i)
        s[i] = (i % 3) * CPU_ScalarType(0.1234) - CPU_ScalarType(0.5);   //'random' starting vector

    detail::copy_vec_to_vec(s,r);

    //std::cout << s << std::endl;

    double epsilon = tag.factor();
    CPU_ScalarType norm = norm_2(r);
    CPU_ScalarType norm_prev = 0;
    long numiter = 0;

    for (vcl_size_t i=0; i<tag.max_iterations(); ++i)
    {
        if (std::fabs(norm - norm_prev) / std::fabs(norm) < epsilon)
            break;

        r /= norm;
        r2 = viennacl::linalg::prod(matrix, r);  //using helper vector r2 for the computation of r <- A * r in order to avoid the repeated creation of temporaries
        r = r2;
        norm_prev = norm;
        norm = norm_2(r);
        numiter++;
    }

    eigenvalue = norm;
    return eigenvalue;
}

    typename viennacl::result_of::cpu_value_type<typename MatrixT::value_type>::type
    eig(MatrixT const& matrix, power_iter_tag const & tag)
    {

      typedef typename viennacl::result_of::value_type<MatrixT>::type           ScalarType;
      typedef typename viennacl::result_of::cpu_value_type<ScalarType>::type    CPU_ScalarType;
      typedef typename viennacl::result_of::vector_for_matrix<MatrixT>::type    VectorT;
    
      CPU_ScalarType eigenvalue;
      long matrix_size = matrix.size1();
      VectorT r(matrix_size);
      std::vector<CPU_ScalarType> s(matrix_size);
      
      for(std::size_t i=0; i<s.size(); ++i)
        s[i] = (i % 3) * CPU_ScalarType(0.1234) - CPU_ScalarType(0.5);   //'random' starting vector

      detail::copy_vec_to_vec(s,r);

      //std::cout << s << std::endl;
      
      double epsilon = tag.factor();
      CPU_ScalarType norm = norm_2(r);
      CPU_ScalarType norm_prev = 0;
      long numiter = 0;

      for (std::size_t i=0; i<tag.max_iterations(); ++i)
      {
        if (std::abs<CPU_ScalarType>(norm - norm_prev) / std::abs<CPU_ScalarType>(norm) < epsilon)
          break; 
           
        r /= norm;
        r = viennacl::linalg::prod(matrix, r);
        norm_prev = norm;
        norm = norm_2(r);
        numiter++;
      }

      eigenvalue = norm;
      return eigenvalue;
    }

    typename viennacl::result_of::cpu_value_type<typename MatrixT::value_type>::type
    eig(MatrixT const& A, power_iter_tag const & tag, VectorT & eigenvec)
    {

      typedef typename viennacl::result_of::value_type<MatrixT>::type           ScalarType;
      typedef typename viennacl::result_of::cpu_value_type<ScalarType>::type    CPU_ScalarType;

      vcl_size_t matrix_size = A.size1();
      VectorT r(eigenvec);
      std::vector<CPU_ScalarType> s(matrix_size);

      for (vcl_size_t i=0; i<s.size(); ++i)
        s[i] = CPU_ScalarType(i % 3) * CPU_ScalarType(0.1234) - CPU_ScalarType(0.5);   //'random' starting vector

      detail::copy_vec_to_vec(s, eigenvec);

      double epsilon = tag.factor();
      CPU_ScalarType norm = norm_2(eigenvec);
      CPU_ScalarType norm_prev = 0;
      long numiter = 0;

      for (vcl_size_t i=0; i<tag.max_iterations(); ++i)
      {
        if (std::fabs(norm - norm_prev) / std::fabs(norm) < epsilon)
          break;

        eigenvec /= norm;
        r = viennacl::linalg::prod(A, eigenvec);  //using helper vector r for the computation of x <- A * x in order to avoid the repeated creation of temporaries
        eigenvec = r;
        norm_prev = norm;
        norm = norm_2(eigenvec);
        numiter++;
      }

      return norm;
    }

      void operator()(SystemType pde_system,
                      DomainType & domain,
                      MatrixT    & system_matrix,
                      VectorT    & load_vector
                     ) const
      {
        typedef typename viennagrid::result_of::cell_tag<DomainType>::type CellTag;

        typedef typename viennagrid::result_of::point<DomainType>::type                                   PointType;
        typedef typename viennagrid::result_of::element<DomainType, CellTag>::type                        CellType;

        typedef typename viennagrid::result_of::element_range<DomainType, CellTag>::type                  CellContainer;
        typedef typename viennagrid::result_of::iterator<CellContainer>::type                             CellIterator;

        typedef typename SystemType::equation_type                  EquationType;

     #ifdef VIENNAFEM_DEBUG
        std::cout << "Strong form: " << pde_system.pde(0) << std::endl;
     #endif
        log_strong_form(pde_system);
        EquationType weak_form_general = viennafem::make_weak_form(pde_system.pde(0));
     #ifdef VIENNAFEM_DEBUG
        std::cout << "* pde_solver::operator(): Using weak form general: " << weak_form_general << std::endl;
     #endif
        std::vector<EquationType> temp(1); temp[0] = weak_form_general;
        log_weak_form(temp, pde_system);
        EquationType weak_form = viennamath::apply_coordinate_system(viennamath::cartesian< PointType::dim >(), weak_form_general);
        //EquationType weak_form = viennamath::apply_coordinate_system(viennamath::cartesian<Config::coordinate_system_tag::dim>(), weak_form_general);
        temp[0] = weak_form;
        log_coordinated_weak_form(temp, pde_system);

     #ifdef VIENNAFEM_DEBUG
        std::cout << "* pde_solver::operator(): Using weak form " << weak_form << std::endl;
        std::cout << "* pde_solver::operator(): Write dt_dx coefficients" << std::endl;
     #endif

        typedef typename reference_cell_for_basis<CellTag, viennafem::lagrange_tag<1> >::type    ReferenceCell;

        //
        // Create accessors for performance in the subsequent dt_dx_handler step
        //

        //viennafem::dtdx_assigner<DomainType, StorageType, ReferenceCell>::apply(domain, storage);

        viennafem::dt_dx_handler<DomainType, StorageType, ReferenceCell>  dt_dx_handler(domain, storage);

        //fill with cell quantities
        CellContainer cells = viennagrid::elements<CellType>(domain);
        for (CellIterator cell_iter = cells.begin();
            cell_iter != cells.end();
            ++cell_iter)
        {
          //cell_iter->print_short();
          //viennadata::access<example_key, double>()(*cell_iter) = i;
          //viennafem::dt_dx_handler<ReferenceCell>::apply(storage, *cell_iter);
          dt_dx_handler(*cell_iter);
        }

     #ifdef VIENNAFEM_DEBUG
        std::cout << "* pde_solver::operator(): Create Mapping:" << std::endl;
     #endif
        std::size_t map_index = create_mapping(storage, pde_system, domain);

     #ifdef VIENNAFEM_DEBUG
        std::cout << "* pde_solver::operator(): Assigned degrees of freedom in domain so far: " << map_index << std::endl;
     #endif
        // resize global system matrix and load vector if needed:
        // TODO: This can be a performance bottleneck for large numbers of segments! (lots of resize operations...)
        if (map_index > system_matrix.size1())
        {
          MatrixT temp = system_matrix;
          ////std::cout << "Resizing system matrix..." << std::endl;
          system_matrix.resize(map_index, map_index, false);
          system_matrix.clear();
          system_matrix.resize(map_index, map_index, false);
          for (typename MatrixT::iterator1 row_it = temp.begin1();
               row_it != temp.end1();
               ++row_it)
          {
            for (typename MatrixT::iterator2 col_it = row_it.begin();
                 col_it != row_it.end();
                 ++col_it)
                 system_matrix(col_it.index1(), col_it.index2()) = *col_it;
          }
        }
        if (map_index > load_vector.size())
        {
          VectorT temp = load_vector;
       #ifdef VIENNAFEM_DEBUG
          std::cout << "Resizing load vector..." << std::endl;
       #endif
          load_vector.resize(map_index, false);
          load_vector.clear();
          load_vector.resize(map_index, false);
          for (std::size_t i=0; i<temp.size(); ++i)
            load_vector(i) = temp(i);
        }

     #ifdef VIENNAFEM_DEBUG
        std::cout << "* pde_solver::operator(): Transform to reference element" << std::endl;
//.........这里部分代码省略.........