matrix.h

强大的矩阵模版类

  // [ 1 4 7 2 5 8 3 6 9 ]
  // </codeblock>
  //
  // <h4>OneD Storage Type Selectors</h4>
  // This specifies a normal dense vector to be used as the OneD
  // part of matrix with array storage.
  //
  //  
  //!tparam: MemLoc - Specify whether the memory used is "owned" by the matrix or if it was provided to the matrix from some external source (with a pointer to some data) - internal
  //!component: type
  //!category: containers, selectors
  //!example: swap_rows.cc, general_matvec_mult.cc
  //!definition: matrix.h
  template <int MemLoc=internal>
  struct dense {
    typedef int size_type;
    enum { id = DENSE, oned_id, ext=MemLoc, issparse=0, index };
  };


  //: Packed storage type selectors
  // This storage type is equivalent to the BLAS/LAPACK packed
  // storage format.
  //
  // The packed storage format is similar to the banded format,
  // except that the storage for each row/column of the band
  // is variable so there is no wasted space. This is better
  // for efficiently storing triangular matrices.
  //
  // <codeblock>
  // [  1   2   3   4   5  ]
  // [  0   6   7   8   9  ]
  // [  0   0  10  11  12  ]
  // [  0   0   0  13  14  ]
  // [  0   0   0   0  15  ]
  //
  // [  1   2   3   4   5  ]
  // [  6   7   8   9  ]
  // [  10  11  12  ]
  // [  13  14  ]
  // [  15  ]
  //
  // mapped to linear memory with row-major order:
  //
  // [  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  ]
  // </codeblock>
  //
  //!tparam: MemLoc - Specify whether the memory used is "owned" by the matrix or if it was provided to the matrix from some external source (with a pointer to some data) - internal
  //!component: type
  //!category: containers, selectors
  //!example: tri_pack_vect.cc
  template <int MemLoc=internal>
  struct packed { 
    typedef int size_type;
    enum { id = PACKED, oned_id, ext=MemLoc, issparse=0, index }; 
  };

  //: Banded view storage type selectors
  // This storage type is used for creating matrices that
  // are "views" into existing full dense matrices.
  // For instance, one could create a triangular view
  // of a full matrix.
  //!tparam: MemLoc - Specify whether the memory used is "owned" by the matrix or if it was provided to the matrix from some external source (with a pointer to some data) - internal
  //!component: type
  //!category: containers, selectors
  //!definition: matrix.h
  template <int MemLoc=internal>
  struct banded_view {
    typedef int size_type;
    enum{id = BAND_VIEW,oned_id, ext=MemLoc, issparse=0, index }; 
  };

  //: Compressed sparse storage type selectors (for both OneD and TwoD)
  //
  // <h4> TwoD Storage Type Selectors</h4>
  // This storage type is the traditional compressed row
  // or compressed column format.
  //
  // The storage consists of three arrays, one array for
  // all of the elements, one array consisting of the
  // row or column index (row for column-major and column
  // for row-major matrices), and one array consisting
  // of pointers to the start of each row/column.
  //
  // The following is an example sparse matrix in compressed for
  // format, with the stored indices specified as
  // <tt>index_from_one</tt>. Note that the MTL interface is still
  // indexed from zero whether or not the underlying stored indices
  // are from one.
  //
  // <codeblock>
  // [  1      2   3      ]
  // [     4           5  ]
  // [  6      7   8      ]
  // [  9         10      ]
  // [        11      12  ]
  //
  // row pointer array
  // [  1  4  6  9  11 13 ]
  // 
  // element value array
  // [  1  2  3  4  5  6  7  8  9 10 11 12 ]
  //
  // element column index array
  // [  1  3  4  2  5  1  3  4  1  4  3  5 ]
  // </codeblock>
  //
  // Of course, the user of the MTL sparse matrix does not
  // need to concern his or herself with the implementation
  // details of this matrix storage format. The interface
  // to an MTL compressed row matrix is the same as that
  // of any MTL matrix, as described in Matrix.
  //
  // <h4> OneD Storage Type Selectors</h4>
  // This is a OneD type used to construct array matrices.
  // The compressed OneD format uses two arrays, one to hold
  // the elements of the vector, and the other to hold the
  // indices that coorespond to each element (either their
  // row or column number).
  //
  //
  //!tparam: SizeType - The type used in the index and pointer array - int
  //!tparam: MemLoc - Specify whether the memory used is "owned" by the matrix or if it was provided to the matrix from some external source (with a pointer to some data) - internal
  //!tparam: IndexStyle - Specify whether the underlying index array stores indices starting from one (fortan style) or from zero (c-style)  - index_from_zero
  //!component: type
  //!category: containers, selectors
  //!example: sparse_matrix.cc
  //!definition: matrix.h
  template <class SizeType = int, int MemLoc=internal, 
            int IndexStyle = index_from_zero>
  struct compressed {
    typedef SizeType size_type;
    enum { id = COMPRESSED, oned_id, ext=MemLoc, 
           issparse=1, index=IndexStyle }; 
  };

  //: Array storage type selectors
  //
  // This storage type gives an "array of pointers" style
  // implementation of a matrix. Each row or column of the matrix
  // is allocated separately. The type of vector used for the rows
  // or columns is very flexible, and one can choose from any of the
  // OneD storage types, which include <tt>dense</tt>,
  // <tt>compressed</tt>, <tt>sparse_pair</tt>, <tt>tree</tt>,
  // and <tt>linked_list</tt>.
  //
  // <codeblock>
  // matrix < double,
  //          rectangle<>,
  //          array< dense<> >,
  //          row_major >::type
  // [ ] -> [  1  0  0  4  0 ]
  // [ ] -> [  0  7  8  0  0 ]
  // [ ] -> [ 11  0 13 14  0 ]
  // [ ] -> [ 16  0 18  0 20 ]
  // [ ] -> [  0 22  0 24  0 ]
  //
  // matrix < double, 
  //          rectangle<>,
  //          array< sparse_pair<> >,
  //          row_major >::type
  // [ ] -> [ (1,0) (4,3) ]
  // [ ] -> [ (7,1) (8,2) ]
  // [ ] -> [ (11,0) (13,2) (14,3) ]
  // [ ] -> [ (16,0) (18,2) (20,4) ]
  // [ ] -> [ (22,1) (24,3) ]
  // </codeblock>
  //
  // <p>One advantage of this type of storage is that rows can be
  // swapped in constant time. For instance, one could swap the
  // row 3 and 4 of a matrix in the following way.
  //
  // <codeblock>
  // Matrix::OneD tmp = A[3];
  // A[3] = A[4];
  // A[4] = tmp;
  // </codeblock>
  //
  // The rows are individually reference counted so that the user
  // does not have to worry about deallocating the rows.
  //
  //!tparam: OneD - The storage type used for each row/column of the matrix - dense
  //!tparam: MemLoc - Specify whether the memory used is "owned" by the matrix or if it was provided to the matrix from some external source (with a pointer to some data) - internal
  //!component: type
  //!category: containers, selectors
  //!models: TwoDStorage
  //!definition: matrix.h
  template <class OneD = dense<>, int MemLoc=internal>
  struct array {
    typedef typename OneD::size_type size_type;
    enum { id=ARRAY, oned_id=OneD::id, ext=MemLoc,
           issparse = OneD::issparse, index=index_from_zero }; 
  };

  //: Envelope storage type selectors
  //
  // The storage scheme is for sparse symmetric matrices, where most
  // of the non-zero elements fall near the main diagonal.  The
  // storage format is useful in certain factorizations since the
  // fill-ins fall in areas already allocated.  This scheme is
  // different than most sparse matrices since the row containers are
  // actually dense, similar to a banded matrix.
  //
  // <codeblock>
  // [  1              ]
  // [  2  3           ]
  // [  4     5        ]
  // [        6  7     ]
  // [     8     9  10 ]
  //
  // 
  //          [ 0  2  5  7 11 ]    Diagonals pointer array
  //    _______/__/   |  |__\___________
  //   V     V        V     V           V
  // [ 1  2  3  4  0  5  6  7  8  0  9 10 ] Element values array
  //
  // </codeblock>
  //!tparam: MemLoc - Specify whether the memory used is "owned" by the matrix or if it was provided to the matrix from some external source (with a pointer to some data) - internal
  //!component: type
  //!category: containers, selectors
  //!definition: matrix.h
  template <int MemLoc=internal>
  struct envelope { 
    typedef int size_type;
    enum { id = ENVELOPE, oned_id, ext=MemLoc, issparse = 0, index };
  };

  /* OneD Storage*/


  //: Index-Value Pair Sparse Vector implemented with mtl::dense1D
  // This is a OneD type for constructing array matrices.
  // The implementation is a mtl::dense1D consisting of index-value
  // pairs.
  //!component: type
  //!category: containers, selectors
  //!definition: matrix.h
  struct sparse_pair {
    typedef int size_type;
    enum { id = SPARSE_PAIR, issparse = 1, index=index_from_zero }; 
  };

  //: Index-Value Pair Sparse Vector implemented with std::set
  // This is a OneD type for constructing array matrices.
  // The implementation is a std::set consisting of index-value
  // pairs.
  // 
  //!component: type
  //!category: containers, selectors
  //!definition: matrix.h
  struct tree {
    typedef int size_type;
    enum { id = TREE, issparse = 1, index=index_from_zero}; 
  };

  //: Index-Value Pair Sparse Vector implemented with std::list
  // This is a OneD type for constructing array matrices.
  // The implementation is a std::list consisting of index-value
  // pairs.
  // 
  //!component: type
  //!category: containers, selectors
  //!definition: matrix.h
  struct linked_list { 
    typedef int size_type;
    enum { id = LINKED_LIST, issparse = 1, index=index_from_zero }; 
  };

  //: Matrix Type Generators Error
  // If you see this in a compiler error message it means
  // you made a mistake in selecting arguments for your matrix type.
  //!noindex:
  struct generators_error { };

  //: generators for oned
  //!noindex:
  template <class T, class Storage>
  struct generate_oned {
    enum { Storage_oned_id = Storage::oned_id, 
           Storage_index = Storage::index };
    typedef typename IF< EQUAL< Storage_oned_id,DENSE>::RET,
                         dense1D<T>,
                     typename IF< EQUAL< Storage_oned_id,COMPRESSED>::RET,
                         compressed1D<T, typename Storage::size_type, 
                                         Storage_index>,
                     typename IF< EQUAL< Storage_oned_id,SPARSE_PAIR>::RET,
                         sparse1D< mtl::dense1D< entry1<T> > >,
                     typename IF< EQUAL< Storage_oned_id,TREE>::RET,
                         sparse1D< std::set< entry1<T> > >,
                     typename IF< EQUAL< Storage_oned_id,LINKED_LIST>::RET,
                         sparse1D< std::list< entry1<T> > >,
                         generators_error
                       >::RET
                       >::RET
                       >::RET
                       >::RET
                       >::RET RET;
  };

  //: blah
  //!noindex:
  template <class T, class Storage, int M, int N>
  struct generate_internal {
    enum { Storage_id = Storage::id, Storage_index = Storage::index };
    typedef typename IF< EQUAL< Storage_id, DENSE>::RET,
                         gen_dense2D<T, gen_rect_offset<M,N>,M,N>,
                     typename IF< EQUAL< Storage_id, COMPRESSED >::RET,
                         gen_compressed2D<T, typename Storage::size_type,
                                          Storage_index,M,N>,
                     typename IF< EQUAL< Storage_id, ENVELOPE >::RET,
                         gen_envelope2D<T,M,N>,
                     typename IF< EQUAL< Storage_id, PACKED >::RET,
                         gen_dense2D<T, gen_packed_offset<M,N>,M,N>,
                     typename IF< EQUAL< Storage_id, BAND >::RET,
                         gen_dense2D<T, gen_banded_offset<M,N>,M,N>,
                     typename IF< EQUAL< Storage_id, BAND_VIEW >::RET,
                         gen_dense2D<T, gen_banded_view_offset<M,N>,M,N>,
                     typename IF< EQUAL< Storage_id, ARRAY >::RET,
                         gen_array2D< typename generate_oned<T,Storage>::RET >,
                         generators_error
                     >::RET
                     >::RET
                     >::RET
                     >::RET
                     >::RET
                     >::RET
                     >::RET RET;
  };

  //: blah
  //!noindex:
  template <class T, class Storage, int M, int N>
  struct generate_external {
    enum { Storage_id = Storage::id, Storage_index = Storage::index };
    typedef typename IF< EQUAL< Storage_id, DENSE>::RET,
                         gen_external2D<T, gen_rect_offset<M,N>, M,N>,
                     typename IF< EQUAL< Storage_id, COMPRESSED >::RET,
                         gen_ext_comp2D<T, typename Storage::size_type,
                                        Storage_index,M,N>,
                     typename IF< EQUAL< Storage_id, ENVELOPE >::RET,
                         gen_envelope2D<T,M,N>,
                     typename IF< EQUAL< Storage_id, PACKED >::RET,
                         gen_external2D<T, gen_packed_offset<M,N>,M,N>,
                     typename IF< EQUAL< Storage_id, BAND >::RET,
                         gen_external2D<T, gen_banded_offset<M,N>,M,N>,
                     typename IF< EQUAL< Storage_id, BAND_VIEW >::RET,
                         gen_external2D<T, gen_banded_view_offset<M,N>,M,N>,
                     typename IF< EQUAL< Storage_id, ARRAY >::RET,
                         gen_array2D< typename generate_oned<T,Storage>::RET >,
                         generators_error
                     >::RET
                     >::RET
                     >::RET
                     >::RET
                     >::RET
                     >::RET
                     >::RET RET;
  };

  //: blah
  //!noindex:
  template <class T, class Storage, int M, int N>
  struct generate_storage {
    enum { Storage_ext = Storage::ext };
    typedef typename IF< Storage_ext, 
                         typename generate_external<T,Storage,M,N>::RET,
                         typename generate_internal<T, Storage,M,N>::RET 
                       >::RET RET;
  };

  //: blah