matrix.h

强大的矩阵模版类

  //!noindex:
  template <class T, class Orientation, class Storage, int M, int N>
  struct generate_rect {
    enum { Orientation_id = Orientation::id };
    typedef typename generate_storage<T,Storage,M,N>::RET storage_t;
    typedef typename storage_t::type::size_type size_type;
    typedef typename IF< EQUAL<Orientation_id, ROW_MAJOR>::RET,
                           row_matrix< storage_t,
                              gen_rect_indexer<row_orien,M,N, size_type> >,
                     typename IF< EQUAL<Orientation_id, COL_MAJOR>::RET,
                           column_matrix< storage_t,
                               gen_rect_indexer<column_orien,M,N,size_type> >,
                        generators_error
                       >::RET
                       >::RET RET;
  };

  //: blah
  //!noindex:
  template <class T, class Orientation, class Storage, int M, int N>
  struct generate_banded {
    enum { Orientation_id = Orientation::id };
    typedef typename generate_storage<T,Storage,M,N>::RET storage_t;
    typedef typename storage_t::type::size_type size_type;
    typedef typename IF< EQUAL<Orientation_id, ROW_MAJOR>::RET,
                        row_matrix< storage_t,
                             gen_banded_indexer<row_orien,M,N,size_type> >,
                   typename IF< EQUAL<Orientation_id, COL_MAJOR>::RET,
                          column_matrix< storage_t,
                           gen_banded_indexer<column_orien,N,M,size_type> >,
                           generators_error
                       >::RET
                       >::RET RET;
  };

  //: blah
  //!noindex:
  template <class T, class Storage, int M, int N>
  struct generate_diagonal {
    typedef typename generate_storage<T,Storage,M,N>::RET storage_type;
    typedef typename storage_type::type::size_type size_type;
    typedef gen_diagonal_indexer<row_orien, M, N, size_type> indexer_gen;
    typedef diagonal_matrix<storage_type, indexer_gen> RET;
  };

  //: blah
  //!noindex:
  template <int Uplo>
  struct generate_uplo {
    typedef typename IF< EQUAL< Uplo, upper >::RET, upper__,
          typename IF< EQUAL< Uplo, unit_upper >::RET, unit_upper__,
          typename IF< EQUAL< Uplo, lower>::RET, lower__, 
          typename IF< EQUAL< Uplo, unit_lower>::RET, unit_lower__,
#if defined __SUNPRO_CC
          typename IF< EQUAL< Uplo, 100>::RET, dynamic_uplo__,
#endif
          dynamic_uplo__
#if defined __SUNPRO_CC
      >::RET
#endif
      >::RET
      >::RET
      >::RET
      >::RET RET;
  };

  //: blah
  //!noindex:
  template <class T, class Shape, class Orien, class Storage, int M, int N>
  struct generate_triangle {
    enum { Orien_id = Orien::id, Shape_uplo = Shape::uplo, 
           Storage_issparse = Storage::issparse };
    typedef typename generate_storage<T,Storage,M,N>::RET storage_type;
    typedef typename storage_type::type::size_type size_type;
    typedef typename
     IF< Storage_issparse, 
       typename IF< EQUAL< Orien_id, ROW_MAJOR>::RET,
         triangle_matrix<row_matrix< storage_type,
            gen_banded_indexer<row_orien,M,N,size_type> >, 
            typename generate_uplo<Shape_uplo>::RET >,
         triangle_matrix<column_matrix< storage_type,
            gen_banded_indexer<column_orien,M,N,size_type> >,
            typename generate_uplo<Shape_uplo>::RET >
                     >::RET
     ,
       typename IF< EQUAL< Orien_id, ROW_MAJOR>::RET,
         triangle_matrix<row_matrix< storage_type,
            gen_banded_indexer<row_orien,M,N,size_type> >,
            typename generate_uplo<Shape_uplo>::RET >,
         triangle_matrix<column_matrix< storage_type,
            gen_banded_indexer<column_orien,M,N,size_type> >, 
            typename generate_uplo<Shape_uplo>::RET >
                     >::RET
     >::RET RET;
  };


  /* if storage is sparse, use rect indexer, otherwise use
   banded indexer */

  //: blah
  //!noindex:
  template <class T, class Shape, class Orien, class Storage, int M, int N>
  struct generate_symmetric {
    enum { Shape_uplo = Shape::uplo, Orien_id = Orien::id, 
           Storage_issparse = Storage::issparse };
    typedef typename generate_storage<T,Storage,M,N>::RET storage_type;
    typedef typename storage_type::type::size_type size_type;
    typedef typename
     IF< Storage_issparse, 
       typename IF< EQUAL< Orien_id, ROW_MAJOR>::RET,
       symmetric_matrix<row_matrix< storage_type,
           gen_rect_indexer<row_orien,M,N,size_type> >, 
           typename generate_uplo<Shape_uplo>::RET >,
       symmetric_matrix<column_matrix< storage_type,
           gen_rect_indexer<column_orien,N,M,size_type> >, 
           typename generate_uplo<Shape_uplo>::RET >
                     >::RET
     ,
       typename IF< EQUAL< Orien_id, ROW_MAJOR>::RET,
       symmetric_matrix<row_matrix< storage_type,
           gen_banded_indexer<row_orien,M,N,size_type> >, 
           typename generate_uplo<Shape_uplo>::RET >,
       symmetric_matrix<column_matrix< storage_type,
           gen_banded_indexer<column_orien,N,M,size_type> >, 
           typename generate_uplo<Shape_uplo>::RET >
                     >::RET 
     >::RET RET;
  };

  //: Matrix type generators class.
  //  
  // Matrices that occur in real engineering and scientific applications
  // often have special structure, especially in terms of how many zeros
  // are in the matrix, and where the non-zeros are located in the matrix.
  // This means that space and time saving can be acheived by using various
  // types of compressed storage. There are a multitude of matrix storage
  // formats in use today, and the MTL tries to support many of the more
  // common storage formats. The following discussion will describe how the
  // user of MTL can select the type of matrix he or she wishes to use.
  //
  // To create a MTL matrix, one first needs to construct the
  // appropriate matrix type. This is done using the <tt>matrix</tt>
  // type generation class, which is easier to think of as a function. It
  // takes as input the characteristics of the matrix type that you want
  // and then returns the appropriate MTL matrix. The <tt>matrix</tt> type
  // generators ``function'' has defaults defined, so in order to create a
  // normal rectangular matrix type, one merely does the following:
  //
  // <codeblock>
  // typedef matrix< double >::type MyMatrix;
  // MyMatrix A(M, N);
  // </codeblock>
  //
  // The matrix type generators can take up to four arguments, the
  // element type, the matrix shape, the storage
  // format, and the orientation. The following is the
  // ``prototype'' for the <tt>matrix</tt> type generators.
  //
  // <codeblock>
  // matrix< EltType, Shape, Storage, Orientation >::type 
  // </codeblock>
  //
  // This type of "generative" interface technique was developed by by
  // <a href="http://nero.prakinf.tu-ilmenau.de:80/~czarn/">Krzysztof
  // Czarnecki</a> and <a
  // href="http://home.t-online.de/home/Ulrich.Eisenecker/">Ulrich
  // Eisenecker</a> in their work on the <a href="http://nero.prakinf.tu-ilmenau.de:80/~czarn/gmcl/">Generative Matrix Computation Library</a>.
  //  <p>
  //  *Storage can be made external by specifying such in the storage
  //   parameter. eg. dense&lt;external&gt;, packed&lt;external&gt;.
  //
  //!tparam: EltType - Valid choices for this argument include <tt>double</tt>, <tt>complex<float></tt>, and <tt>bool</tt>.  In essence, any builtin or user defined type can be used for the <tt>EltType</tt>, however, if one uses the matrix with a particular algorithm, the <tt>EltType</tt> must support the operations required by the algorithm. For MTL algorithms these typically include the usual numerical operators such as addition and multiplication. The <tt>std::complex</tt> class is a good example of what is required in a numerical type. The documentation for each algorithm will include the requirements on the element type.
  //  
  //!tparam: Shape - This argument specifies the general positioning of the non zero elements in the matrix, but does not specify the actual storage format.  In addition it specifies certain properties such as symmetry.  The choices for this argument include <tt>rectangle</tt>, <tt>banded</tt>, <tt>diagonal</tt>, <tt>triangle</tt>, and <tt>symmetric</tt>. Hermitian is not yet implemented.
  //  
  //!tparam: Storage - The argument specifies the storage scheme used to lay out the matrix elements (and sometimes the element indices) in memory. The storage formats include <tt>dense</tt> , <tt>banded</tt>, <tt>packed</tt> , <tt>banded_view</tt>, <tt>compressed</tt>, <tt>envelope</tt>, and <tt>array</tt>.
  //
  //!tparam: Orientation - The storage order for an MTL matrix can either be <tt>row_major</tt> or <tt>column_major</tt>.
  //
  //!category: containers, generators
  //!component: type
  //!definition: matrix.h
  //
  template < class T, class Shape = rectangle<>, class Storage = dense<>, class Orientation = row_major >
  struct matrix {
    enum { Shape_id = Shape::id, Shape_M = Shape::M, Shape_N = Shape::N };
    //: The generated type
    typedef typename IF< EQUAL< Shape_id, RECT>::RET,
                        typename generate_rect<T,Orientation,Storage,
                                               Shape_M, Shape_N>::RET,
                    typename IF< EQUAL< Shape_id, DIAG>::RET,
                        typename generate_diagonal<T,Storage,
                                               Shape_M, Shape_N>::RET,
                    typename IF< EQUAL< Shape_id, BAND>::RET,
                        typename generate_banded<T,Orientation,Storage,
                                               Shape_M, Shape_N>::RET,
                    typename IF< EQUAL< Shape_id, TRI>::RET,
                        typename generate_triangle<T,Shape,
                                               Orientation,Storage,
                                               Shape_M, Shape_N>::RET,
                    typename IF< EQUAL< Shape_id, SYMM>::RET,
                        typename generate_symmetric<T,Shape,
                                              Orientation,Storage,
                                               Shape_M, Shape_N>::RET,
                        generators_error
                        >::RET
                        >::RET
                        >::RET
                        >::RET
                        >::RET type;
  };



#ifndef MTL_DISABLE_BLOCKING

  //: Block View Matrix Type Constructor
  //
  // <codeblock>
  // block_view<Matrix> bA = blocked<>(A, 16, 16);
  // or
  // block_view<Matrix,BM,BN> bA = blocked<BM,BN>(A);
  // </codeblock>
  //
  // Note: currently not supported for egcs (internal compiler error).
  //
  //!category: containers, generators
  //!component: type
  //!example: blocked_matrix.cc
  //!definition: matrix.h
  //!tparam: Matrix - The type of the Matrix to be blocked, must be dense
  //!tparam: BM - The blocking factor for the rows (M dimension) - 0 for dynamic size
  //!tparam: BN - The blocking factor for the columns (N dimension) - 0 for dynamic size
  template <class Matrix, int BM = 0, int BN = 0>
  struct block_view {
    //: The generated type
    typedef typename Matrix:: MTL_TEMPLATE blocked_view<BM,BN>::type type;
  };

  template <int BM, int BN>
  struct blk { enum { M = BM, N = BN }; };

  //: Blocked Matrix Generator
  //
  // <codeblock>
  // block_view<Matrix> bA = blocked<>(A, 16, 16);
  // </codeblock>
  //
  // Note: currently not supported for egcs (internal compiler error).
  //
  //!category: containers
  //!component: function
  //!example: blocked_matrix.cc
  //!definition: matrix.h
  //!tparam: Matrix - The type of the Matrix to be blocked, must be dense
  template <class Matrix>
  typename block_view<Matrix, 0, 0>::type 
  blocked(const Matrix& A, int bm, int bn) {
    typedef dimension<typename Matrix::size_type, 0,0> bdt;
    return typename block_view<Matrix, 0, 0>::type(A, bdt(bm, bn));
  }

  //: Blocked Matrix Generator
  //
  // This version of the blocked matrix generator is for statically
  // sized blocks.
  //
  // <codeblock>
  // block_view<Matrix,BM,BN> bA = blocked<BM,BN>(A);
  // </codeblock>
  //
  // Note: currently not supported for egcs (internal compiler error).
  //
  //!category: containers
  //!component: function
  //!example: blocked_matrix.cc
  //!definition: matrix.h
  //!tparam: Matrix - The type of the Matrix to be blocked, must be dense
  //!tparam: BM - The blocking factor for the rows (M dimension)
  //!tparam: BN - The blocking factor for the columns (N dimension)
  template <class Matrix, int BM, int BN>
  typename block_view<Matrix, BM, BN>::type
  blocked(const Matrix& A, blk<BM,BN>) {
    typedef dimension<typename Matrix::size_type, BM, BN> bdt;
    return typename block_view<Matrix, BM, BN>::type(A, bdt(BM, BN));
  }
#endif

  //: Band View Matrix Type Constructor
  // A helper class for creating a banded_view into an existing dense matrix.
  //!category: containers, generators
  //!component: type
  //!example: banded_view_test.cc
  //!definition: matrix.h
  //!tparam: Matrix - The type of the Matrix to be viewed, must be dense
  template <class Matrix>
  struct band_view {
    //: The generated type
    typedef typename Matrix::banded_view_type type;
  };

  //: Triangle View Matrix Type Constructor
  // A helper class for creating a triangle view into an existing dense
  //   or sparse matrix. For sparse matrices, the matrix must already
  //   have elements in the appropriate triangular portion of the matrix.
  //   This just provides the proper triangular matrix interface.
  //!category: containers, generators
  //!component: type
  //!example: banded_view_test.cc
  //!definition: matrix.h
  //!tparam: Matrix - The type of the Matrix to be viewed, must be dense
  //!tparam: Uplo - Whether to view the upper or lower triangle of the matrix
  template <class Matrix, int UL>
  struct triangle_view {
    typedef typename Matrix::sparsity Sparsity;
    enum { Sparsity_id = Sparsity::id };
    //: The generated type
    typedef typename generate_uplo<UL>::RET up_or_lower;
    typedef typename IF< EQUAL< Sparsity_id, DENSE>::RET,
			 triangle_matrix<typename Matrix::banded_view_type, up_or_lower>,
                         triangle_matrix< Matrix, up_or_lower >
                       >::RET type;
  };

  //: Triangle View Creation Helper Fuctor
  // Example: tri_view<upper>()(A)
  //!category: containers, generators
  //!component: type
  //!example: banded_view_test.cc
  //!definition: matrix.h
  //!tparam: Uplo - Whether to view the upper or lower triangle of the matrix
  template <int Uplo>
  struct tri_view {
    template <class Matrix>
    inline typename triangle_view<Matrix, Uplo>::type operator()(Matrix x) const {
      typedef typename triangle_view<Matrix, Uplo>::type TriView;
      return TriView(x);
    }
  };

  

  //: Symmetric View Matrix Type Constructor
  // A helper class for creating a symmetric view into an existing dense
  //   or sparse matrix. For sparse matrices, the matrix must already
  //   have elements in the appropriate lower/upper portion of the matrix.
  //   This just provides the proper symmetric matrix interface.
  //!category: containers, generators
  //!component: type
  //!definition: matrix.h
  //!tparam: Matrix - The type of the Matrix to be viewed, must be dense
  //!tparam: Uplo - Whether to view the upper or lower triangle of the matrix
  template <class Matrix, int Uplo>
  struct symmetric_view {
    typedef typename Matrix::sparsity Sparsity;
          enum { Sparsity_id = Sparsity::id };
    //: The generated type
    typedef typename IF< EQUAL< Sparsity_id, DENSE>::RET,
      symmetric_matrix< typename Matrix::banded_view_type, 
                        typename generate_uplo<Uplo>::RET>,
                symmetric_matrix< Matrix, typename generate_uplo<Uplo>::RET>
              >::RET type;
  };


} /* namespace mtl */

#endif /* _MTL_MATRIX_H_ */