matrix.h

强大的矩阵模版类

// -*- c++ -*-
//
// Copyright 1997, 1998, 1999 University of Notre Dame.
// Authors: Andrew Lumsdaine, Jeremy G. Siek, Lie-Quan Lee
//
// This file is part of the Matrix Template Library
//
// You should have received a copy of the License Agreement for the
// Matrix Template Library along with the software;  see the
// file LICENSE.  If not, contact Office of Research, University of Notre
// Dame, Notre Dame, IN  46556.
//
// Permission to modify the code and to distribute modified code is
// granted, provided the text of this NOTICE is retained, a notice that
// the code was modified is included with the above COPYRIGHT NOTICE and
// with the COPYRIGHT NOTICE in the LICENSE file, and that the LICENSE
// file is distributed with the modified code.
//
// LICENSOR MAKES NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED.
// By way of example, but not limitation, Licensor MAKES NO
// REPRESENTATIONS OR WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY
// PARTICULAR PURPOSE OR THAT THE USE OF THE LICENSED SOFTWARE COMPONENTS
// OR DOCUMENTATION WILL NOT INFRINGE ANY PATENTS, COPYRIGHTS, TRADEMARKS
// OR OTHER RIGHTS.
//
//===========================================================================


#ifndef _MTL_MATRIX_H_
#define _MTL_MATRIX_H_


/* namespace polution from <sys/sysmacros.h> */
#undef major
#undef minor

#include <set>
#include <list>
#include <vector>

#include "mtl/meta_equal.h"
#include "mtl/meta_if.h"

#include "mtl/matrix_traits.h"
#include "mtl/compressed1D.h"

#include "mtl/matrix_implementation.h"
#include "mtl/rect_indexer.h"
#include "mtl/banded_indexer.h"
#include "mtl/diagonal_indexer.h"

#include "mtl/dense2D.h"
#include "mtl/array2D.h"
#include "mtl/compressed2D.h"
#include "mtl/envelope2D.h"

#include "mtl/linalg_vec.h"
#include "mtl/sparse1D.h"
#include "mtl/entry.h"
#include "mtl/uplo.h"

namespace mtl {

  /* Shapes */

  //: rectangle shape type selector
  //
  // A MTL rectangular matrix is one in which elements could appear in
  // any position in the matrix, i.e., there can be any element
  // <tt>A(i,j)</tt> where <tt>0 <= i <= M</tt> and <tt>0 <= i <=
  // M</tt>.  Both dense and sparse matrices can fit into this
  // category.  The compatible storage types are <tt>dense</tt>,
  // <tt>compressed</tt>, and <tt>array</tt>.  If the matrix size is known
  // at compile time, the matrix size can be built in to the shape
  // (for dense external matrices). This allows some algorithms to
  // optimize for small matrices.  Here are a few examples of creating
  // rectangular matrix types:
  //  
  // <codeblock>
  // typedef matrix < double,
  //                  rectangle<>, 
  //                  dense<>,
  //                  column_major >::type FortranMatrix;
  //  
  // typedef matrix < double,
  //                  rectangle<>, 
  //                  compressed<>,
  //                  row_major >::type CSR;
  //  
  // typedef matrix < double,
  //                  rectangle<4,4>, 
  //                  dense<>,
  //                  row_major >::type Quaternion;
  //  
  // typedef matrix < double,
  //                  rectangle<>, 
  //                  array< compressed<> >,
  //                  row_major >::type SparseArrayMat;
  // </codeblock>
  //
  //!tparam: MM - The number of rows of the matrix, if the matrix has static size (known at compile time) - not static
  //!tparam: NN - The number of columns of the matrix, if the matrix has static size (known at compile time) - not static
  //!component: type
  //!category: containers,selectors
  //!definition: matrix.h
  template <int MM = 0, int NN = 0> 
  class rectangle {
  public:
    enum { M = MM, N = NN, id = RECT, uplo };
  };

  //: banded shape type selectors, also banded storage type selectors
  //
  // <h4>Shape Type Selectors</h4>
  // A banded shape matrix is one in which non-zero matrix elements only
  // appear within a ``band'' of a matrix. The bandwidth of a matrix is
  // described by the number of diagonals in the band that are below the
  // main diagonal, refered to as the <i>sub</i> diagonals, and the number
  // of diagonals in the band above the main diagonal, refered to as the
  // <i>super</i> diagonals. The following is an example of a matrix with a
  // bandwidth of (1,2).
  //
  // <codeblock>
  // [  1   2   3   0   0 ]
  // [  4   5   6   7   0 ]
  // [  0   8   9  10  11 ]
  // [  0   0  12  13  14 ]
  // [  0   0   0  15  16 ]
  // </codeblock>
  //
  // There are many storage types that can be used to efficiently
  // represent banded matrices. The MTL storage types that can be
  // used are <tt>banded</tt>, <tt>packed</tt>, <tt>banded_view</tt>,
  // and <tt>array</tt>. Here are some examples of creating
  // banded matrix types:
  //
  // <codeblock>
  // typedef matrix < double,
  //                  banded<>, 
  //                  packed<>,
  //                  column_major >::type BLAS_Packed;
  //
  // typedef matrix < double,
  //                  banded<>, 
  //                  banded<>,
  //                  column_major >::type BLAS_Banded;
  // </codeblock>
  //
  // <h4>Storage Type Selectors</h4>
  // <tt>banded</tt> is also the type selectors for the banded storage format.
  // This storage format is equivalent to the banded storage used in the
  // BLAS and LAPACK. Similar to the <tt>dense</tt> storage format, a
  // single contiguous chunk of memory is allocated. The banded storage
  // format maps the bands of the matrix to a twod-array of dimension (sub
  // + super + 1) by min(M, N + sub). In MTL the 2D array can be row or
  // column major (for the BLAS it is always column major). The twod-array
  // is then in turn mapped the the linear memory space of the single chunk
  // of memory. The following is an example banded matrix with the mapping
  // to the row-major and column-major 2D arrays. The x's represent
  // memory locations that are not used.
  //
  // <codeblock>
  // [  1   2   3   0   0   0  ]
  // [  4   5   6   7   0   0  ]
  // [  0   8   9  10  11   0  ]
  // [  0   0  12  13  14  15  ]
  // [  0   0   0  16  17  18  ]
  // [  0   0   0   0  19  20  ]
  //
  // row-major
  // [  1   2   3   x  ]
  // [  4   5   6   7  ]
  // [  8   9  10  11  ]
  // [ 12  13  14  15  ]
  // [  x  16  17  18  ]
  // [  x   x  19  20  ]
  //
  // column-major
  // [  x   x   3   7  11  15  ]
  // [  x   2   6  10  14  18  ]
  // [  1   5   9  13  17  20  ]
  // [  4   8  12  16  19   x  ]
  // </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 banded {
    typedef int size_type;
    enum { id = BAND, oned_id, uplo, ext=MemLoc, M=0, N=0, 
           issparse=0, index }; 
  };

  //: diagonal shape type selectors
  //
  // The diagonal matrix shape is similar to the banded matrix in that
  // there is a bandwidth that describes the area of the matrix in which
  // non-zero matrix elements can reside. The difference between the banded
  // matrix shape lies in how the MTL iterators traverse the matrix, which
  // is explained in DiagonalMatrix. The MTL storage types that
  // can be used are <tt>banded</tt>, <tt>packed</tt>,
  // <tt>banded_view</tt>, and <tt>array</tt>. To get the traditional
  // tridiagonal matrix format, one just has to specify the bandwith to be
  // (1,1) and use the <tt> array &lt; dense &lt; &gt; &gt;</tt> storage format.
  //
  //!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 diagonal { 
    enum { uplo, id = DIAG, ext=MemLoc, M=0, N=0 }; 
  };

  //: triangle shape type selectors
  //
  // The triangular shape is a special case of the banded shape. There
  // are four kinds of triangular matrices in MTL, based on the
  // <tt>Uplo</tt> argument:
  //
  // <table border=1>
  // <tr> <td> Uplo type </td> <td> Sub </td> <td> Super </td> </tr>
  // <tr> <td> upper </td> <td> 0 </td> <td> N - 1 </td> </tr>
  // <tr> <td> unit_upper </td> <td> -1 </td> <td> N - 1 </td> </tr>
  // <tr> <td> lower </td> <td> M - 1 </td> <td> 0 </td> </tr>
  // <tr> <td> unit_lower </td> <td> M - 1 </td> <td> -1 </td> </tr>
  // </table>
  //
  // The following is an example of a <tt>triangle&lt;upper&gt;</tt> shaped
  // matrix:
  //
  // <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  ]
  // </codeblock>
  //
  // The next example is of a <tt>triangle&lt;unit_lower&gt;</tt>
  // matrix.  The main diagonal is not stored, since it consists of
  // all ones.  The MTL algorithms recognize when a matrix is ``unit''
  // and perform a slightly different operation to take this into
  // account.  The ones will not show up in an iteration of the
  // matrix, and access to the A(i,i) element of a unit lower/upper
  // matrix is an error.
  //
  // <codeblock>
  // [  1   0   0   0   0  ]
  // [  1   1   0   0   0  ]
  // [  2   3   1   0   0  ]
  // [  4   5   6   1   0  ]
  // [  7   8   9  10   1  ]
  // </codeblock>
  //
  //
  // Here are a couple examples of creating some triangular matrix types:
  //
  // <codeblock>
  // typedef matrix < double,
  //                  triangle<upper>, 
  //                  banded<>,
  //                  column_major >::type UpperTriangle;
  //
  // typedef matrix < double,
  //                 triangle<unit_lower>, 
  //                 packed<>,
  //                 row_major >::type UnitLowerTriangle;
  // </codeblock>
  //
  //!tparam: Uplo - The type of triangular matrix. Either upper, lower, unit_upper, or unit_lower.
  //!component: type
  //!category: containers, selectors
  //!definition: matrix.h
  template <int Uplo = dynamic_uplo> 
  struct triangle { 
    enum { id = TRI, uplo = Uplo, M=0, N=0 }; 
  };



  //: symmetric shape type selectors
  //
  // Symmetric matrices are similar to banded matrices in that there
  // is only access to a particular band of the matrix. The difference
  // is that in an MTL symmetric matrix, <tt>A(i,j)</tt> and
  // <tt>A(j,i)</tt> refer to the same element. The following is an
  // example of a symmetric matrix:
  //
  // <codeblock>
  // the full symmetric matrix
  // [  1   2   3   4   5  ]
  // [  2   6   7   8   9  ]
  // [  3   7  10  11  12  ]
  // [  4   8  11  13  14  ]
  // [  5   9  12  14  15  ]
  //
  // the symmetric matrix in packed storage
  // [  1  ]
  // [  2   6  ]
  // [  3   7  10  ]
  // [  4   8  11  13  ]
  // [  5   9  12  14  15  ]
  // </codeblock>
  //
  // Similar to the triangle shape, the user must provide an
  // <tt>Uplo</tt> argument which specifies which part of the matrix is
  // actually stored.  The valid choices are <tt>upper</tt> and
  // <tt>lower</tt> for symmetric matrices.
  //
  // <codeblock>
  // typedef matrix < double,
  //                  symmetric<lower>, 
  //                  packed<>,
  //                  row_major >::type SymmMatrix;
  // </codeblock>
  //
  //!component: type
  //!category: containers, selectors
  //!tparam: Uplo - The portion of the matrix that is stored. Either upper or lower.
  //!example: symm_packed_vec_prod.cc, symm_banded_vec_prod.cc
  template <int Uplo = dynamic_uplo>
  struct symmetric { 
    enum { id = SYMM, uplo = Uplo, M=0, N=0 };
  };

  //: hermitian shape type selectors
  // Hermitian matrices are not yet implemented.
  //!component: type
  //!category: containers, selectors
  //!tparam: Uplo - The portion of the matrix that is stored. Either upper or lower.
  //!definition: matrix.h
  template <int Uplo = dynamic_uplo> 
  struct hermitian {
    enum { id = HERM, uplo = Uplo, M=0, N=0 };
  };

  /* Ordering */
  //: Row major ordering type selectors
  //!category: containers, selectors
  struct row_major {
    enum { id = ROW_MAJOR }; 
  };

  //: Column major ordering type selectors
  //!category: containers, selectors
  struct column_major { 
    enum { id = COL_MAJOR };
  };

  /* Storage */

  //: Dense storage type selectors (for both TwoD and OneD storage)
  //
  // <h4>TwoD Storage Type Selectors</h4>
  //This is the most common way of storing matrices, and consists of
  //one contiguous piece of memory that is divided up into rows or
  //columns of equal length.  The following example shows how a matrix
  //can be mapped to linear memory in either a row-major or
  //column-major fashion.
  //  
  // <codeblock>
  // [ 1 2 3 ]
  // [ 4 5 6 ]
  // [ 7 8 9 ]
  //  
  // row major:
  // [ 1 2 3 4 5 6 7 8 9 ]
  //  
  // column major: