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

      void operator()(LinPdeSysT const & pde_system,
                      SegmentT   const & segment,
                      StorageType      & storage,
                      MatrixT          & system_matrix,
                      VectorT          & load_vector)
      {
        typedef viennamath::equation                          equ_type;
        typedef viennamath::expr                              expr_type;
        typedef typename expr_type::interface_type            interface_type;
        typedef typename expr_type::numeric_type              numeric_type;

        typedef typename viennagrid::result_of::cell_tag<SegmentT>::type CellTag;

        std::size_t map_index = viennafvm::create_mapping(pde_system, segment, storage);

        system_matrix.clear();
        system_matrix.resize(map_index, map_index, false);
        load_vector.clear();
        load_vector.resize(map_index);


        for (std::size_t pde_index = 0; pde_index < pde_system.size(); ++pde_index)
        {
#ifdef VIENNAFVM_DEBUG
          std::cout << std::endl;
          std::cout << "//" << std::endl;
          std::cout << "//   Equation " << pde_index << std::endl;
          std::cout << "//" << std::endl;
#endif
          assemble(pde_system, pde_index,
                   segment, storage,
                   system_matrix, load_vector);

        } // for pde_index
      } // functor

void make_rotation_matrix(MatrixT & mat,
                          vcl_size_t new_size,
                          vcl_size_t off_diagonal_distance = 4)
{
  mat.resize(new_size, new_size, false);
  mat.clear();

  double val = 1.0 / std::sqrt(2.0);

  for (vcl_size_t i=0; i<new_size; ++i)
    mat(i,i) = val;

  for (vcl_size_t i=off_diagonal_distance; i<new_size; ++i)
  {
    mat(i-off_diagonal_distance, i)                       = val;
    mat(i,                       i-off_diagonal_distance) = -val;
  }

}

      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;
//.........这里部分代码省略.........

void RowEchelon<T>::getNullspace(MatrixT& N) const
{
  if(R.isEmpty()) {
    N.clear();
    return;
  }
  Assert((int)firstEntry.size() == R.m+1);

  int nullspace_dims=getNull();
  N.resize(R.n,nullspace_dims);

  //first get nullspace vectors from 0 to firstEntry[0]
  int i,j,numVecs=0;
  int m=R.m,n=R.n;
  for(j=0;j<firstEntry[0];j++) {
    N.setCol(numVecs,0); N(j,numVecs)=1;
    numVecs++;
  }
  for(i=0;i<m;i++) {
    //cancel out the i'th entry
    for(j=firstEntry[i]+1;j<firstEntry[i+1];j++) {
      if(numVecs >= N.n) {
	LOG4CXX_INFO(KrisLibrary::logger(),"Num nullspace vectors "<<numVecs);
	LOG4CXX_INFO(KrisLibrary::logger(),"Found more nullspace vectors than found dims, row "<<i);
	LOG4CXX_INFO(KrisLibrary::logger(),MatrixPrinter(R));
      }
      Assert(numVecs < N.n);
      VectorT xn; N.getColRef(numVecs,xn);
      xn.setZero();
      xn[firstEntry[i]] = R(i,j);
      xn[j] = -R(i,firstEntry[i]);

      //cancel out all the entries prior to i 
      int isave=i;
      i--;
      for(;i>=0;i--) {
        VectorT ri; R.getRowRef(i,ri);
	//calculate alpha
	int ji=firstEntry[i];
	Assert(ji != n);
	int ji2=firstEntry[i+1]; //(i+1==m?n:firstEntry[i+1]);
	T alpha;
	if(ji2 == n) alpha = Zero;
        else {
          VectorT rji2; rji2.setRef(ri,ji2,1,R.n-ji2);
          VectorT xji2; xji2.setRef(xn,ji2,1,R.n-ji2);
          alpha = xji2.dot(rji2);
        }
	xn[ji] = -alpha/ri[ji];
      }
      i=isave;
      numVecs++;
    }
  }
  if(numVecs != nullspace_dims) {
    LOG4CXX_ERROR(KrisLibrary::logger(),"Error in counting rank in row-eschelon decomposition");
    LOG4CXX_INFO(KrisLibrary::logger(),"Num nullspace vectors "<<numVecs);
    LOG4CXX_INFO(KrisLibrary::logger(),MatrixPrinter(R));
  }
  Assert(numVecs == nullspace_dims);

  /*
  VectorT temp;
  for(int i=0;i<numVecs;i++) {
    VectorT xi; N.getColRef(i,xi);
    xi.print();
    R.mul(xi,temp);
    if(temp.maxAbsElement() > 1e-4) {
      LOG4CXX_INFO(KrisLibrary::logger(),"Nullspace vector "<<i<<" not in null space!");
      xi.print();
      LOG4CXX_INFO(KrisLibrary::logger(),"Result = "); temp.print();
      KrisLibrary::loggerWait();
    }
  }
  */

  VectorT* N0 = new VectorT[nullspace_dims];
  for(int i=0;i<nullspace_dims;i++) N.getColRef(i,N0[i]);
  int num=OrthonormalBasis(N0,N0,nullspace_dims);
  Assert(num == nullspace_dims);
  delete [] N0;
}