2016-02-14

疎行列クラスの基礎 Suitesparse の cholmodを使ってみる

Ubuntu Suitesparse cholmod

前の記事の疎行列用クラスのラッパークラス(spinxform)のヘッダー定義と関数コードを眺めて、疎行列クラスがどのようになっているかを見てみる
疎行列クラスの基本
- nxm行列の要素を全部並べると大きすぎる。疎であることを利用して、「番地」と「値」の組だけにしたい
- 「番地」は整数ペア
- 「値」は複素数を取らせることにすると実数ペア
- これをc++、STLらしくイテレータを使って実装している
ヘッダー定義を見よう
- その中でstd::pair, std::mapを使っているところを抜粋する
  - 行列の(列、行)整数ペアをEntryIndexとし
  - 行列のセルの複素数値をdoubleのペアとしてEntryValueとしている
  - その上で、(行、列) -> (double,double)というmapを作りEntryMapとしている
  - このままだとイテレータが使えないので、EntryMpaをイテレーションできるようにiterator, const_iteratorの２つを作っている。片方は変更可能に、もう片方は固定
  - その２つのイテレータが使えるように、先頭と終端とが取れるように関数も定義してある

         typedef std::pair<int,int> EntryIndex; // NOTE: column THEN row! (makes it easier to build compressed format)
         typedef std::pair<double,double> EntryValue;
         typedef std::map<EntryIndex,EntryValue> EntryMap;
         typedef EntryMap::iterator       iterator;
         typedef EntryMap::const_iterator const_iterator;

               iterator begin( void );
         const_iterator begin( void ) const;
               iterator   end( void );
         const_iterator   end( void ) const;

- その実用を、疎行列の転置処理関数で見てみる
  - 転置後には、転置前と同じ数の番地ペアが必要なので、それ用にEntryMap transposedを作り、そこの番地(整数ペア)と値(実数ペア)を納めていく
  - 転置前の疎行列のすべての要素をイテレーションするためにconst_iteratorを使っている
  - dataは、Sparseクラスに持たせているcholmod_sparseクラスなので、その番地の始点・終点取り出しをしていることもわかる
  - std::mapの要素取り出しは <,>と二重入れ子になっているので、first.firtsと言えば、最初のintが、first.secondと言えば二番目のintが、secondと言えば、というペアが取り出される
  - 転置処理なので、int,intの前後をひっくり返して登録している

   void Sparse :: transpose( void )
   {
      EntryMap transposed;

      for( const_iterator e = data.begin(); e != data.end(); e++ )
      {
         int i = e->first.first;
         int j = e->first.second;
         transposed[ EntryIndex( j, i ) ] = e->second;
      }

      data = transposed;
      swap( m, n );
   }

- 以下がSparseクラスの定義全体

#include <cholmod.h>
#include <map>
   class Sparse
   {
      public:
         Sparse( Common& common, int m = 0, int n = 0, int xtype = CHOLMOD_REAL );
         // initialize an mxn matrix of doubles
         // xtype is either CHOLMOD_REAL or CHOLMOD_ZOMPLEX

         ~Sparse( void );

         void resize( int m, int n );
         // clears and resizes to mxn matrix
         
         cholmod_sparse* operator*( void );
         // dereference operator gets pointer to underlying cholmod_sparse data structure

         int size( int dim ) const;
         // returns the size of the dimension specified by scalar dim

         int length( void ) const;
         // returns the size of the largest dimension

         void zero( double rVal = 0., double iVal = 0. );
         // sets all nonzero elements to rVal+iVal*i
         
         void transpose( void );
         // replaces this matrix with its transpose

         void horzcat( const Sparse& A, const Sparse& B );
         // replaces the current matrix with [ A, B ]

         void vertcat( const Sparse& A, const Sparse& B );
         // replaces the current matrix with [ A; B ]

         double& operator()( int row, int col );
         double  operator()( int row, int col ) const;
         double& r( int row, int col );
         double  r( int row, int col ) const;
         // access real part of element (row,col)
         // note: uses 0-based indexing

         double& i( int row, int col );
         double  i( int row, int col ) const;
         // access imaginary part of element (row,col)
         // note: uses 0-based indexing
         
         typedef std::pair<int,int> EntryIndex; // NOTE: column THEN row! (makes it easier to build compressed format)
         typedef std::pair<double,double> EntryValue;
         typedef std::map<EntryIndex,EntryValue> EntryMap;
         typedef EntryMap::iterator       iterator;
         typedef EntryMap::const_iterator const_iterator;

               iterator begin( void );
         const_iterator begin( void ) const;
               iterator   end( void );
         const_iterator   end( void ) const;

      protected:
         Common& common;
         int m, n;
         int xtype, stype;
         cholmod_sparse* A;
         EntryMap data;

         double& retrieveEntry( int row, int col, int c );
         double  retrieveEntry( int row, int col, int c ) const;
   };