cholesky.hpp

来自「矩阵运算源码最新版本」· HPP 代码 · 共 580 行 · 第 1/2 页

// Software License for MTL
// 
// Copyright (c) 2007 The Trustees of Indiana University. All rights reserved.
// Authors: Peter Gottschling and Andrew Lumsdaine
// 
// This file is part of the Matrix Template Library
// 
// See also license.mtl.txt in the distribution.

#ifndef MTL_CHOLESKY_INCLUDE
#define MTL_CHOLESKY_INCLUDE

#include <boost/numeric/mtl/recursion/matrix_recursator.hpp>
#include <boost/numeric/mtl/utility/glas_tag.hpp>
#include <boost/numeric/mtl/operation/dmat_dmat_mult.hpp>
#include <boost/numeric/mtl/operation/assign_mode.hpp>
#include <boost/numeric/mtl/matrix/transposed_view.hpp>
#include <boost/numeric/mtl/recursion/base_case_cast.hpp>

namespace mtl {

namespace with_bracket {

    // ============================================================================
    // Generic Cholesky factorization and operands for Cholesky on with submatrices
    // ============================================================================
	
    template < typename Matrix > 
    void cholesky_base (Matrix & matrix)
    {
	for (int k = 0; k < matrix.num_cols(); k++) {
	    matrix[k][k] = sqrt (matrix[k][k]);
	    
		for (int i = k + 1; i < matrix.num_rows(); i++) {
		    matrix[i][k] /= matrix[k][k];
		    typename Matrix::value_type d = matrix[i][k];

		    for (int j = k + 1; j <= i; j++)
			matrix[i][j] -= d * matrix[j][k];
		}
	}
    }
	
    
    template < typename MatrixSW, typename MatrixNW > 
    void tri_solve_base(MatrixSW & SW, const MatrixNW & NW)
    {
	for (int k = 0; k < NW.num_rows (); k++) {
	    
	    for (int i = 0; i < SW.num_rows (); i++) {
		SW[i][k] /= NW[k][k];
		typename MatrixSW::value_type d = SW[i][k];
		
		for (int j = k + 1; j < SW.num_cols (); j++)
		    SW[i][j] -= d * NW[j][k];
	    }
	}
    }
    

    // Lower(SE) -= SW * SW^T
    template < typename MatrixSE, typename MatrixSW > 
    void tri_schur_base(MatrixSE & SE, const MatrixSW & SW)
    {
	for (int k = 0; k < SW.num_cols (); k++)
	    
	    for (int i = 0; i < SE.num_rows (); i++) {
		    typename MatrixSW::value_type d = SW[i][k];
		    for (int j = 0; j <= i; j++)
			SE[i][j] -= d * SW[j][k];
	    }
    }


    template < typename MatrixNE, typename MatrixNW, typename MatrixSW >
    void schur_update_base(MatrixNE & NE, const MatrixNW & NW, const MatrixSW & SW)
    {
	for (int k = 0; k < NW.num_cols (); k++) 
	    for (int i = 0; i < NE.num_rows (); i++) {
		typename MatrixNW::value_type d = NW[i][k];
		for (int j = 0; j < NE.num_cols (); j++)
		    NE[i][j] -= d * SW[j][k];
	    }
    }


    // ======================
    // Corresponding functors
    // ======================
    
    struct cholesky_base_t
    {
	template < typename Matrix > 
	void operator() (Matrix & matrix)
	{
	    cholesky_base(matrix);
	}
    };
    
    struct tri_solve_base_t
    {
	template < typename MatrixSW, typename MatrixNW > 
	void operator() (MatrixSW & SW, const MatrixNW & NW)
	{
	    tri_solve_base(SW, NW);
	}
    };

    struct tri_schur_base_t
    {
	template < typename MatrixSE, typename MatrixSW > 
	void operator() (MatrixSE & SE, const MatrixSW & SW)
	{
	    tri_schur_base(SE, SW);
	}
    };

    struct schur_update_base_t
    {
	template < typename MatrixNE, typename MatrixNW, typename MatrixSW >
	void operator() (MatrixNE & NE, const MatrixNW & NW, const MatrixSW & SW)
	{
	    schur_update_base(NE, NW, SW);
	}
    };

} // namespace with_bracket


namespace with_iterator {

    // ============================================================================
    // Generic Cholesky factorization and operands for Cholesky on with submatrices
    // ============================================================================
    

    template < typename Matrix > 
    void cholesky_base (Matrix& matrix)
    {
	using namespace glas::tag; using traits::range_generator;
	typedef tag::iter::all    all_it;

        typedef typename Matrix::value_type                         value_type;
        typedef typename range_generator<col, Matrix>::type       cur_type;             
        typedef typename range_generator<all_it, cur_type>::type    iter_type;            

	typedef typename range_generator<row, Matrix>::type       rcur_type;
	typedef typename range_generator<all_it, rcur_type>::type   riter_type;   
	
	typename Matrix::size_type k= 0;
	for (cur_type kb= begin<col>(matrix), kend= end<col>(matrix); kb != kend; ++kb, ++k) {

	    iter_type ib= begin<all_it>(kb), iend= end<all_it>(kb); 
	    ib+= k; // points now to matrix[k][k]

	    value_type root= sqrt (*ib);
	    *ib= root;

	    ++ib; // points now to matrix[k+1][k]
	    rcur_type rb= begin<row>(matrix); rb+= k+1; // to row k+1
	    for (int i= k + 1; ib != iend; ++ib, ++rb, ++i) {
		*ib = *ib / root;
		typename Matrix::value_type d = *ib;
		riter_type it1= begin<all_it>(rb);    it1+= k+1;      // matrix[i][k+1]
		riter_type it1end= begin<all_it>(rb); it1end+= i+1;   // matrix[i][i+1]
		iter_type it2= begin<all_it>(kb);     it2+= k+1;      // matrix[k+1][k]
		for (; it1 != it1end; ++it1, ++it2)
		    *it1 = *it1 - d * *it2;
	    }
	}
    }

    
    template < typename MatrixSW, typename MatrixNW > 
    void tri_solve_base(MatrixSW & SW, const MatrixNW & NW)
    {
	using namespace glas::tag; using traits::range_generator;
	typedef tag::iter::all        all_it;
	typedef tag::const_iter::all  all_cit;

        typedef typename range_generator<col, MatrixNW>::type       ccur_type;             
        typedef typename range_generator<all_cit, ccur_type>::type    citer_type;            

	typedef typename range_generator<row, MatrixSW>::type       rcur_type;
	typedef typename range_generator<all_it, rcur_type>::type     riter_type;   

	for (int k = 0; k < NW.num_rows (); k++) 
	    for (int i = 0; i < SW.num_rows (); i++) {

		typename MatrixSW::value_type d = SW[i][k] /= NW[k][k];

		rcur_type sw_i= begin<row>(SW);     sw_i+= i;  // row i
		riter_type it1= begin<all_it>(sw_i);  it1+= k+1; // SW[i][k+1]
		riter_type it1end= end<all_it>(sw_i);    
	
		ccur_type nw_k= begin<col>(NW);     nw_k+= k;  // column k
		citer_type it2= begin<all_cit>(nw_k); it2+= k+1; // NW[k+1][k]

		for(; it1 != it1end; ++it1, ++it2)
		    *it1 = *it1 - d * *it2;
	    }
    }
    

    // Lower(SE) -= SW * SW^T
    template < typename MatrixSE, typename MatrixSW > 
    void tri_schur_base(MatrixSE & SE, const MatrixSW & SW)
    {
	using namespace glas::tag; using traits::range_generator;
	typedef tag::iter::all        all_it;
	typedef tag::const_iter::all  all_cit;

        typedef typename range_generator<col, MatrixSW>::type       ccur_type;             
        typedef typename range_generator<all_cit, ccur_type>::type    citer_type;            

	typedef typename range_generator<row, MatrixSE>::type       rcur_type;
	typedef typename range_generator<all_it, rcur_type>::type     riter_type;   

	for (int k = 0; k < SW.num_cols (); k++)
	    for (int i = 0; i < SE.num_rows (); i++) {
		typename MatrixSW::value_type d = SW[i][k];

		rcur_type se_i= begin<row>(SE);       se_i+= i;      // row i
		riter_type it1= begin<all_it>(se_i);                   // SE[i][0]
		riter_type it1end= begin<all_it>(se_i); it1end+= i+1;  // SE[i][i+i]

		ccur_type sw_k= begin<col>(SW);     sw_k+= k;        // column k
		citer_type it2= begin<all_cit>(sw_k);                  // SW[0][k]

		for(; it1 != it1end; ++it1, ++it2)
		    *it1 = *it1 - d * *it2;
	    }
    }


    template < typename MatrixNE, typename MatrixNW, typename MatrixSW >
    void schur_update_base(MatrixNE & NE, const MatrixNW & NW, const MatrixSW & SW)
    {
	using namespace glas::tag; using traits::range_generator;
	typedef tag::iter::all        all_it;
	typedef tag::const_iter::all  all_cit;

        typedef typename range_generator<col, MatrixSW>::type       ccur_type;             
        typedef typename range_generator<all_cit, ccur_type>::type    citer_type;            

	typedef typename range_generator<row, MatrixNE>::type       rcur_type;
	typedef typename range_generator<all_it, rcur_type>::type     riter_type;   

	for (int k = 0; k < NW.num_cols (); k++) 
	    for (int i = 0; i < NE.num_rows (); i++) {
		typename MatrixNW::value_type d = NW[i][k];

		rcur_type ne_i= begin<row>(NE);       ne_i+= i;      // row i
		riter_type it1= begin<all_it>(ne_i);                   // NE[i][0]
		riter_type it1end= end<all_it>(ne_i);                  // NE[i][num_col]

		ccur_type sw_k= begin<col>(SW);     sw_k+= k;        // column k
		citer_type it2= begin<all_cit>(sw_k);                  // SW[0][k]

		for (int j = 0; j < NE.num_cols (); j++)
		    NE[i][j] -= d * SW[j][k];
	    }
    }


    // ======================
    // Corresponding functors
    // ======================
    
    struct cholesky_base_t
    {
	template < typename Matrix > 
	void operator() (Matrix & matrix)
	{
	    cholesky_base(matrix);
	}
    };
    
    struct tri_solve_base_t
    {
	template < typename MatrixSW, typename MatrixNW > 
	void operator() (MatrixSW & SW, const MatrixNW & NW)
	{
	    tri_solve_base(SW, NW);
	}
    };

    struct tri_schur_base_t
    {
	template < typename MatrixSE, typename MatrixSW >