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

    vector_type solve(const matrix_type& A, const vector_type& y)
    {
        typedef typename matrix_type::size_type size_type;
        typedef typename matrix_type::value_type value_type;

        namespace ublas = boost::numeric::ublas;

        matrix_type Q(A.size1(), A.size2()), R(A.size1(), A.size2());

        qr (A, Q, R);

        vector_type b = prod(trans(Q), y);

        vector_type result;
        if (R.size1() > R.size2())
        {
            size_type min = (R.size1() < R.size2() ? R.size1() : R.size2());

            result = ublas::solve(subrange(R, 0, min, 0, min),
                                  subrange(b, 0, min),
                                  ublas::upper_tag());
        }
        else
        {
            result = ublas::solve(R, b, ublas::upper_tag());
        }
        return result;
    }

// Vectorizes (stacks columns of) an input matrix.
state_type Vec(const matrix_type &A)
{
    state_type vectorized(A.size1()*A.size2());
    
    for (int n = 0; n < A.size1(); n++) {
        for (int m = 0; m < A.size2(); m++) {
            vectorized(n*A.size1()+m) = A(m,n);
        }
    }
    
    return vectorized;
}

// Tensor product. Probably a better way of doing this.
matrix_type Kronecker(const matrix_type &A, const matrix_type &B)
{
    matrix_type C(A.size1()*B.size1(), A.size2()*B.size2());
    
    for (int i=0; i < A.size1(); i++) {
        for (int j=0; j < A.size2(); j++) {
            for (int k=0; k < B.size1(); k++) {
                for (int l=0; l < B.size2(); l++) {
                    C(i*B.size1()+k, j*B.size2()+l) = A(i,j)*B(k,l);
                }
            }
        }
    }
    
    return C;
}

void qr(const matrix_type& A, matrix_type& Q, matrix_type& R)
{
    using namespace boost::numeric::ublas;

    typedef typename matrix_type::size_type size_type;
    typedef typename matrix_type::value_type value_type;

    // TODO resize Q and R to match the needed size.
    int m=A.size1();
    int n=A.size2();

    identity_matrix<value_type> ident(m);
    if (Q.size1() != ident.size1() || Q.size2() != ident.size2())
        Q = matrix_type(m, m);
    Q.assign(ident);

    R.clear();
    R = A;

    for (size_type k=0; k< R.size1() && k<R.size2(); k++)
    {
        slice s1(k, 1, m - k);
        slice s2(k, 0, m - k);
        unit_vector<value_type> e1(m - k, 0);

        // x = A(k:m, k);
        matrix_vector_slice<matrix_type> x(R, s1, s2);
        matrix_type F(x.size(), x.size());
        
        Reflector(x, F);

        matrix_type temp = subrange(R, k, m, k, n);
        //F = prod(F, temp);
        subrange(R, k, m, k, n) = prod(F, temp);

        // <<---------------------------------------------->>
        // forming Q
        identity_matrix<value_type> iqk(A.size1());
        matrix_type Q_k(iqk);
        
        subrange(Q_k, Q_k.size1() - F.size1(), Q_k.size1(),
                 Q_k.size2() - F.size2(), Q_k.size2()) = F;

        Q = prod(Q, Q_k);
    }
}

	/** @brief Constructs a tensor with a matrix
	 *
	 * \note Initially the tensor will be two-dimensional.
	 *
	 *  @param v matrix to be copied.
	 */
	BOOST_UBLAS_INLINE
	tensor (const matrix_type &v)
		: tensor_expression_type<self_type>()
		, extents_ ()
		, strides_ ()
		, data_    (v.data())
	{
		if(!data_.empty()){
			extents_ = extents_type{v.size1(),v.size2()};
			strides_ = strides_type(extents_);
		}
	}

	static void calc(matrix_type& a, vector_type& v, const char& jobz) {
		namespace impl = boost::numeric::bindings::lapack::detail;
		//impl::syev(jobz,a,v);

		char      uplo  = 'L';
		integer_t n     = a.size2();
		integer_t lda   = a.size1();
		integer_t lwork =-1;
		integer_t info  = 0;
		value_type dlwork;
		impl::syev(jobz, uplo, n, mtraits::data(a), lda, vtraits::data(v), &dlwork, lwork, info);
	}

void test_getrf_getrs(matrix_type& lu, matrix_type& x)
{
  typedef typename matrix_type::value_type value_type ;

  numerics::matrix< value_type > a( lu ) ; // tmp to verify result
  numerics::matrix< value_type > b( x ) ;  // tmp to verify result
  std::vector< int > ipiv( x.size1() ) ;

  boost::numeric::bindings::lapack::getrf( lu, ipiv ) ;
  matrix_type ia( lu );
  boost::numeric::bindings::lapack::getrs( 'N', lu, ipiv, x ) ;
  boost::numeric::bindings::lapack::getri( ia, ipiv ) ;

  std::cout << prod(a,x) - b << std::endl ; 
  std::cout << prod(a,ia) << std::endl ; 
}

box<dimension>::box(matrix_type const& edges)
{
    if (edges.size1() != dimension || edges.size2() != dimension) {
        throw std::invalid_argument("edge vectors have invalid dimensionality");
    }
    edges_ = zero_matrix_type(dimension, dimension);
    for (unsigned int i = 0; i < dimension; ++i) {
        edges_(i, i) = edges(i, i);
    }
    if (norm_inf(edges_ - edges) != 0) {
        throw std::invalid_argument("non-cuboid geomtries are not implemented");
    }
    for (unsigned int i = 0; i < dimension; ++i) {
        length_[i] = edges_(i, i);
    }
    length_half_ = 0.5 * length_;

    LOG("edge lengths of simulation domain: " << length_);
}

// Right multiplying operator to superoperator.
matrix_type RightOp(const matrix_type &A)
{
    return Kronecker(trans(A), identity_matrix< complex<double> >(A.size1()));
}

// Left multiplying operator to superoperator.
matrix_type LeftOp(const matrix_type &A)
{
    return Kronecker(identity_matrix< complex<double> >(A.size1()),A);
}

 static std::ptrdiff_t num_rows (matrix_type& cm) { return cm.size1(); }

 static int num_rows (matrix_type& cm) { return cm.size1(); }

 static int size1 (matrix_type& cm) { return cm.size1(); } 

 static index_pointer index1_storage (matrix_type& cm) {
   assert (cm.filled1() == layout_type::size1 (cm.size1(), cm.size2()) + 1);
   return vector_traits<idx_array_t>::storage (cm.index1_data()); 
 }

 static int size1 (matrix_type& hm) { return hm.size1(); }