C++ MatrixT类代码示例

    void pretty_print_matrix(std::ostream &ostr, MatrixT const &mat)
    {
        IndexType rows, cols;
        mat.get_shape(rows, cols);
        typename MatrixT::ScalarType zero(mat.get_zero());

        for (IndexType row = 0; row < rows; ++row)
        {
            ostr << ((row == 0) ? "[[" : " [");
            if (cols > 0)
            {
                auto val = mat.get_value_at(row, 0);
                if (val == zero)
                    ostr << " ";
                else
                    ostr << val;
            }

            for (IndexType col = 1; col < cols; ++col)
            {
                auto val = mat.get_value_at(row, col);
                if (val == zero)
                    ostr << ",  ";
                else
                    ostr << ", " << val;
            }
            ostr << ((row == rows - 1) ? "]]\n" : "]\n");
        }
    }

void
test_transpose_readonly(MatrixT view)
{
  typedef typename MatrixT::value_type T;

  // Check that view is initialized
  check_matrix(view, 0);

  length_type const size1 = view.size(1);

  typename MatrixT::const_transpose_type trans = view.transpose();

  test_assert(trans.size(0) == view.size(1));
  test_assert(trans.size(1) == view.size(0));

  for (index_type idx0=0; idx0<trans.size(0); ++idx0)
    for (index_type idx1=0; idx1<trans.size(1); ++idx1)
    {
      T expected = T(idx1 * size1 + idx0 + 0);
      test_assert(equal(trans.get(idx0, idx1), expected));
      test_assert(equal(trans.get(idx0,  idx1),
		   view. get(idx1, idx0)));
      }

  // Check that view is unchanged
  check_matrix(view, 0);
}

void LDLDecomposition<T>::getPseudoInverse(MatrixT& Ainv) const
{
  Ainv.resize(LDL.n,LDL.n);
  VectorT temp(LDL.n,Zero),y,x;
  for(int i=0;i<LDL.n;i++) {
    temp(i)=One;
    LBackSub(temp,y);
    for(int j=0;j<y.n;j++) {
      if(!FuzzyZero(LDL(j,j),zeroTolerance))
	y(j) = y(j)/LDL(j,j);
      else
	y(j) = 0.0;
    }
    LTBackSub(y,x);
    //fill in a column
    for(int j=0;j<LDL.n;j++)
      Ainv(j,i)=x(j);
    temp(i)=Zero;
  }

  T tol = Ainv.maxAbsElement()*Epsilon;
  for(int i=0;i<LDL.n;i++)
    for(int j=0;j<i;j++) {
      if(!FuzzyEquals(Ainv(i,j),Ainv(j,i),tol))
	LOG4CXX_INFO(KrisLibrary::logger(),Ainv);
      Assert(FuzzyEquals(Ainv(i,j),Ainv(j,i),tol));
      Ainv(i,j)=Ainv(j,i) = 0.5*(Ainv(i,j)+Ainv(j,i));
    }
}

      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

bool QRDecomposition<T>::set(const MatrixT& A)
{
  QR.copy(A);
  tau.resize(Min(A.m,A.n));

  for (int i=0;i<Min(A.m,A.n);i++) {
    /* Compute the Householder transformation to reduce the j-th
       column of the matrix to a multiple of the j-th unit vector */
    
    VectorT c_full,c;
    QR.getColRef(i,c_full);
    c.setRef(c_full,i);
    T tau_i = HouseholderTransform (c);
    tau(i)=tau_i;
    
    /* Apply the transformation to the remaining columns and
       update the norms */
    
    if (i+1 < A.n) {
      MatrixT m;
      m.setRef(QR,i,i+1);
      HouseholderPreMultiply (tau_i, c, m);
    }
  }
  return true;
}

 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;
 }

void DiagonalMatrixTemplate<T>::postMultiplyInverse(const MatrixT& a,MatrixT& x) const
{
  Assert(this->n == a.n);
  x.resize(a.m,this->n);
  MyT xrow,arow;
  for(int i=0;i<a.m;i++) {
    x.getRowRef(i,xrow);
    a.getRowRef(i,arow);
    xrow.componentDiv(arow,*this);
  }
}

void NRQRDecomposition<T>::getQ(MatrixT& Q) const
{
  int n=c.n;
  Q.resize(n,n);
  VectorT Qi;
  Q.set(0);
  for(int i=0;i<n;i++) {
    Q.getRowRef(i,Qi);
    Qi(i)=1;
    QBackSub(Qi,Qi);
  }
}

void backwardSolve(MatrixT const & R, VectorT const & y, VectorT & x)
{
  for (long i2 = static_cast<long>(R.size2())-1; i2 >= 0; i2--)
  {
    vcl_size_t i = static_cast<vcl_size_t>(i2);
    x(i) = y(i);
    for (vcl_size_t j = static_cast<vcl_size_t>(i)+1; j < R.size2(); ++j)
      x(i) -= R(i,j)*x(j);

    x(i) /= R(i,i);
  }
}

void CholeskyDecomposition<T>::getInverse(MatrixT& Ainv) const
{
  Ainv.resize(L.n,L.n);
  VectorT temp(L.n,Zero),y,x;
  for(int i=0;i<L.n;i++) {
    Ainv.getColRef(i,x);  //x &= col i of A
    temp(i)=One;
    LBackSub(temp,y);
    LTBackSub(y,x);
    temp(i)=Zero;
  }
}

void
fill_matrix(MatrixT view, int offset=0)
{
  typedef typename MatrixT::value_type T;

  length_type const size1 = view.size(1);

  // Initialize view

  for (index_type idx0=0; idx0<view.size(0); ++idx0)
    for (index_type idx1=0; idx1<view.size(1); ++idx1)
      view.put(idx0, idx1, T(idx0 * size1 + idx1 + offset));
}

 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
 }

void convolute_3d_out_of_place(MatrixT& _image, MatrixT& _kernel) {

  if (_image.size() != _kernel.size()) {
    std::cerr << "received image and kernel of mismatching size!\n";
    return;
  }

  unsigned M, N, K;
  M = _image.shape()[0];
  N = _image.shape()[1];
  K = _image.shape()[2];

  unsigned fft_size = M * N * (K / 2 + 1);

  // setup fourier space arrays
  fftwf_complex* image_fourier = static_cast<fftwf_complex*>(
      fftwf_malloc(sizeof(fftwf_complex) * fft_size));
  fftwf_complex* kernel_fourier = static_cast<fftwf_complex*>(
      fftwf_malloc(sizeof(fftwf_complex) * fft_size));
  float scale = 1.0 / (M * N * K);

  // define+run forward plans
  fftwf_plan image_fwd_plan = fftwf_plan_dft_r2c_3d(
      M, N, K, _image.data(), image_fourier, FFTW_ESTIMATE);
  fftwf_execute(image_fwd_plan);

  fftwf_plan kernel_fwd_plan = fftwf_plan_dft_r2c_3d(
      M, N, K, _kernel.data(), kernel_fourier, FFTW_ESTIMATE);
  fftwf_execute(kernel_fwd_plan);

  // multiply
  for (unsigned index = 0; index < fft_size; ++index) {
    float real = image_fourier[index][0] * kernel_fourier[index][0] -
                 image_fourier[index][1] * kernel_fourier[index][1];
    float imag = image_fourier[index][0] * kernel_fourier[index][1] +
                 image_fourier[index][1] * kernel_fourier[index][0];
    image_fourier[index][0] = real;
    image_fourier[index][1] = imag;
  }

  fftwf_destroy_plan(kernel_fwd_plan);
  fftwf_destroy_plan(image_fwd_plan);

  fftwf_plan image_rev_plan = fftwf_plan_dft_c2r_3d(
      M, N, K, image_fourier, _image.data(), FFTW_ESTIMATE);
  fftwf_execute(image_rev_plan);

  for (unsigned index = 0; index < _image.num_elements(); ++index) {
    _image.data()[index] *= scale;
  }

  fftwf_destroy_plan(image_rev_plan);
  fftwf_free(image_fourier);
  fftwf_free(kernel_fourier);
}

 void MatrixT<ValueType>::transform_point(ValueType out[4], const MatrixT& m, const ValueType in[4])
 {
    out[0] = m.get(0, 0) * in[0] + m.get(0, 1) * in[1] + m.get(0, 2) * in[2] + m.get(0, 3) * in[3];
    out[1] = m.get(1, 0) * in[0] + m.get(1, 1) * in[1] + m.get(1, 2) * in[2] + m.get(1, 3) * in[3];
    out[2] = m.get(2, 0) * in[0] + m.get(2, 1) * in[1] + m.get(2, 2) * in[2] + m.get(2, 3) * in[3];
    out[3] = m.get(3, 0) * in[0] + m.get(3, 1) * in[1] + m.get(3, 2) * in[2] + m.get(3, 3) * in[3];
 }

void
check_matrix(MatrixT view, int offset=0)
{
  typedef typename MatrixT::value_type T;

  length_type const size1 = view.size(1);

  for (index_type idx0=0; idx0<view.size(0); ++idx0)
    for (index_type idx1=0; idx1<view.size(1); ++idx1)
    {
      test_assert(equal(view.get(idx0, idx1),
		   T(idx0 * size1 + idx1 + offset)));
    }
}