C++ matrix_type类代码示例

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

  /** 
   * 13. 
   * @brief Copy the contents of a matrix.
   * @param[out] A A is overwritten with B.
   * @param[in] B The matrix to be copied.
   */
  static void copy (matrix_type &A, const matrix_type &B) { 

    /** Result matrices should always be of the right size */
    A.Resize (B.Height(), B.Width());

    /** This one is pretty simple, but the order is different */
    elem::Copy (B, A);
  }

  /** 
   * 15. 
   * @brief Extract the diagonal elements of a matrix.
   * @param[in] A The (square) matrix whose diagonal is to be extracted.
   * @param[out] B A diagonal matrix containing entries from A.
   */
  static void diag(const matrix_type& A,
                   matrix_type& B) {
    /* create a zeros matrix of the right dimension and assign to B */
    const int n = A.Width();
    B = zeros(n, n);

    /* Set the diagonals of B */
    for (int i=0; i<n; ++i) B.Set(i, i, A.Get(i,i));
  }

  /** 
   * 6. 
   * @brief Compute the matrix transpose.
   * @param[in] A The matrix to be transposed.
   * @param[out] B B is overwritten with A^{T}
   */
  static void transpose (const matrix_type& A, 
                         matrix_type& B) {

    /** Result matrices should always have sufficient space */
    B.Resize(A.Width(), A.Height());

    /** Compute transpose */
    elem::Transpose (A, B);
  }

  /** 
   * 5. 
   * @brief Multiply one matrix with another.
   * @param[in] A The first matrix
   * @param[in] B The second matrix
   * @param[out] C C is overwritten with (A*B)
   */
  static void multiply (const matrix_type& A,
                        const matrix_type& B,
                        matrix_type& C) {

    /** Result matrices should always have sufficient space */
    C.Resize(A.Height(), B.Width());

    /** We have to do a Gemm */
    elem::Gemm (elem::NORMAL, elem::NORMAL, 1.0, A, B, 0.0, C);
  }

	/** @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_);
		}
	}

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

  /** 
   * 7. 
   * @brief Compute element-wise negation of a matrix.
   * @param[in] A The matrix to be negated.
   * @param[out] B B is overwritten with -1.0*A
   */
  static void negation (const matrix_type& A,
                        matrix_type& B) {

    /** Result matrices should always have sufficient space */
    B.Resize(A.Height(), A.Width());

    /** Copy over the matrix */
    elem::Copy (A, B);

    /** Multiply by -1.0 */
    elem::Scal(-1.0, B);
  }

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

    /**
     * Apply columnwise the sketching transform that is described by the
     * the transform with output sketch_of_A.
     */
    void apply (const matrix_type& A,
                output_matrix_type& sketch_of_A,
                columnwise_tag dimension) const {

        const value_type *a = A.LockedBuffer();
        El::Int lda = A.LDim();
        value_type *sa = sketch_of_A.Buffer();
        El::Int ldsa = sketch_of_A.LDim();

        for (El::Int j = 0; j < A.Width(); j++)
            for (El::Int i = 0; i < data_type::_S; i++)
                sa[j * ldsa + i] = a[j * lda + data_type::_samples[i]];
    }

  /** 
   * 4. 
   * @brief Subtract one matrix from another.
   * @param[in] A The first matrix
   * @param[in] B The second matrix
   * @param[out] C C is overwritten with (A-B)
   */
  static void minus (const matrix_type& A,
                     const matrix_type& B, 
                     matrix_type& C) {

    /** Result matrices should always have sufficient space */
    C.Resize(A.Height(), A.Width());

    /** first copy the matrix over */
    elem::Copy (A, C);

    /** now, subtract the other matrix in */
    elem::Axpy (-1.0, B, C);
  }

    bool is_positive_definite( const matrix<T,N,A>& m )
    {
        typedef matrix<T,N,A> matrix_type;
        typedef typename matrix_type::range_type range_type;

        if ( m.row() != m.col() ) return false;
        
        for ( std::size_t i = 1; i != m.row(); ++i )
        {
            const matrix_type a{ m, range_type{0,i}, range_type{0,i} };
            if ( a.det() <= T(0) ) return false;
        }

       return true; 
    }

        // --------------------------------------------------------------------
        static matrix_type _compute_rooted_fa(const matrix_type & fa)
        {
            const size_t K = fa.get_height();
            const size_t J = fa.get_width();

            assert(K > 1);

            matrix_type rooted_fa (K - 1, J);

            for (size_t k = 0; k + 1 < K; k++)
                for (size_t j = 0; j < J; j++)
                    rooted_fa(k, j) = fa(k + 1, j) - fa(0, j);

            return rooted_fa;
        }

    /**
     * Apply the sketching transform on A and write to sketch_of_A.
     * Implementation for columnwise.
     */
    void apply_impl_vdist(const matrix_type& A,
        output_matrix_type& sketch_of_A,
        skylark::sketch::columnwise_tag tag) const {

        // Just a local operation on the Matrix
        _local.apply(A.LockedMatrix(), sketch_of_A.Matrix(), tag);
    }

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