ラッパーを借用して手軽にcholmod UbuntuでSuitesparse Suitesparse の cholmodを使ってみる

• cholmod_solve, cholmode_l_solveは結構大変
  • どう大変か、というと
  • こちらからもわかるように

cholmod_dense * 	cholmod_solve (int sys, cholmod_factor *L, cholmod_dense *B, cholmod_common *Common)

• • solveする式のパターンを決めるint sys (たとえばAx =B を解くならsys=0)、AをFactorizationした結果のL、それにBを与え、最期に、cholmod的に扱うために与える *Commonを引数として、cholmodeでの密行列(密ベクトル)が返る
  • factorizationして、それをfactorizatio結果に特化したクラス(cholmod_factor)に入れなおして…
  • と、大変
• この大変さをspinxformライブラリでは、ラッパーで簡単に扱えるようにしている
  • もちろん、cholmod, suitesparseのすべてのクラスをラッピングするのは大変なので、ごく一部だけだが
    • ラッピングされているのは
      • 密行列(DENSE):cholmod_denseをラッピング
      • 疎行列(SPARSE):cholmod_sparse
      • Factorizationの結果(FACTOR):cholmod_factor
      • 行列(UPPER):cholmod_upper。これはSPARSEを継承している。基本は行列の形があって、値が全部ゼロの完全に疎な行列。それに、非ゼロセルだけ値を持たせやすい、実(疎)行列
      • cholmodが持つ、cholmod扱いするためのオブジェクト(COMMON):cholmod_common
  • このラッパーの特徴は
    • ラッピングしているcholmodのクラスを *hogeと呼び出せる(*が、そのような機能を持つオペレータ)
    • 行列のサイズ情報を取り出したり変更したりしやすい

#include "cholmod.h"
#include "CMWrapperY.h"
#include <iostream>

using namespace std;



int main()
{
	int n = 3; // n x n 行列を作る
	cm::Common cc; // CMWrapperのnamespaceはcm。そこでcholmod的ハンドリングのためにCOMMON対応オブジェクトを作る
	cm::Upper A(cc,n,n); // n x n 行列を作成し、値を付置する
	cm::Sparse S(cc,3,4);
	cout << "S size:" << S.size(1) << "\n" << endl;
	A(0,0) = 2.0;
	A(0,1) = -1.0;
	A(0,2) = 0.0;
	A(1,0) = -1.0;
	A(1,1) = 3.0;
	A(1,2) = -1.0;
	A(2,0) = 0.0;
	A(2,1) = -1.0;
	A(2,2) = 2.0;
	//A.transpose();
	for(int i=0;i<n;i++){
		for(int j=0;j<n;j++){
			cout << A(i,j) << "\n" << endl;
		}
	}
	cm::Factor B(cc); // FACTORを作るのもccを介する
	B.build(A); // FactorはUpperを引数にbuild()関数を使うようにラッパーが作られているのでこのようにする

	cm::Dense x(cc,n,1); // Ax = bを解くべく、xを１列の密行列として作る
	cm::Dense b(cc,n,1); // bを１列の密行列として作る、そして付置する
	b(0,0) = 1.0;
	b(1,0) = 2.0;
	b(2,0) = 3.0;
	x = cholmod_l_solve(CHOLMOD_A, *B, *b, cc);
	for(int i=0;i<n;i++){
		cout << "i:" << i << " " << x(i,0) << "\n" << endl;
		cout << "i:" << i << " " << x(i,1) << "\n" << endl;
	}
	return 0;
}

• • これを可能にするCMWrapperY.h,CMWrapperY.cppが

// ============================================================================
// SpinXForm -- CMWrapper.h
// Keenan Crane
// August 16, 2011
//
// CMWrapper provides an object-oriented interface for building C++
// applications that interface with CHOLMOD or other SuiteSparse packages using
// the CHOLMOD matrix format.  It includes wrapper classes for CHOLMOD dense
// and sparse matrices as well as cholmod_common objects.  Wrapper objects are
// mainly for convenience and may not always be your most efficient option.
// CHOLMOD has been tested with SuiteSparse 3.4.0 with METIS 4.0.1 and Intel
// TBB 2.2 on Max OS X 10.6.2.
// 
// In order to link against CHOLMOD you must:
// 
//    1. (optional:) download the Intel Threading Building Blocks library from
//       http://www.threadingbuildingblocks.org/
// 
//    2. download METIS from http://glaros.dtc.umn.edu/gkhome/metis/metis/overview
// 
//    3. download SuiteSparse from http://www.cise.ufl.edu/research/sparse/SuiteSparse/
// 
//    4. copy (and build, if necessary) the following Intel TBB libraries into
//       your library path:
//          libtbb.dylib
//          libtbb_debug.dylib
//          libtbbmalloc.dylib
//          libtbbmalloc_debug.dylib
// 
//    5. build and copy the following METIS libraries into your library path:
//          libmetis.a
// 
//    6. build and copy the following SuiteSparse libraries into your library path:
//          libamd.a
//          libcamd.a
//          libccolamd.a
//          libcholmod.a
//          libcolamd.a
//       Note: may be simpler to just type "make" in each of the appropriate
//       SuiteSparse/*/Lib subdirectories -- you do not have to build all of
//       SuiteSparse
//

#ifndef SPINYFORM_CMWRAPPER_H
#define SPINYFORM_CMWRAPPER_H

#include <cholmod.h>
#include <map>

// Object-oriented wrapper for CHOLMOD sparse matrix format.

namespace cm
{
   class Common
   {
      public:
         Common( void );
         ~Common( void );

         operator cholmod_common*( void );
         // allows cm::Common to be treated as a cholmod_common*

      protected:
         cholmod_common common;
   };

   class Dense
   {
      public:
         Dense( 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
         
         Dense( const Dense& A );
         // copy constructor
         
         ~Dense( void );
         // destructor

         cholmod_dense* operator*( void );
         // dereference operator gets pointer to underlying cholmod_dense data structure

         const Dense& operator=( const Dense& A );
         // copies A

         const Dense& operator=( cholmod_dense* A );
         // gets pointer to A; will deallocate A upon destruction

         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 elements to rVal+iVal*i

         double norm( void );
         // returns the infinity norm

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

         void vertcat( const Dense& A, const Dense& 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
         
         double& operator()( int index );
         double  operator()( int index ) const;
         // access real part of element ind of a vector
         // note: uses 0-based indexing
         
         double& r( int index );
         double  r( int index ) const;
         // access real part of element ind of a vector
         // note: uses 0-based indexing
         
         double& i( int index );
         double  i( int index ) const;
         // access imaginary part of element ind of a vector
         // note: uses 0-based indexing

      protected:
         void initializeFromCopy( void );

         Common& common;
         int m, n;
         int xtype;
         cholmod_dense* data;
         double* rData;
         double* iData;
   };

   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;
   };

   class Upper : public Sparse
   {
      public:
         Upper( 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
   };

   class Factor
   {
      public:
         Factor( Common& common );
         ~Factor( void );

         void build( Upper& A );
         // factorizes positive-definite matrix A using CHOLMOD

         cholmod_factor* operator*( void );
         // dereference operator gets pointer to underlying cholmod_factor data structure

      protected:
         Common& common;
         cholmod_factor *L;
   };
}

#endif

// ============================================================================
// SpinXForm -- CMWrapper.cpp
// Keenan Crane
// August 16, 2011
//

#include "CMWrapperY.h"
#include <algorithm>
#include <cassert>
#include <cmath>

using namespace std;

namespace cm
{
   Common :: Common( void )
   {
      cholmod_l_start( &common );
   }

   Common :: ~Common( void )
   {
      cholmod_l_finish( &common );
   }

   Common :: operator cholmod_common*( void )
   {
      return &common;
   }

   Sparse :: Sparse( Common& _common, int _m, int _n, int _xtype )
   : common( _common ),
     m( _m ),
     n( _n ),
     xtype( _xtype ),
     stype( 0 ),
     A( NULL )
   {}

   Upper :: Upper( Common& _common, int _m, int _n, int _xtype )
   : Sparse( _common, _m, _n, _xtype )
   {
      stype = 1;
   }

   Sparse :: ~Sparse( void )
   {
      if( A )
      {
         cholmod_l_free_sparse( &A, common );
      }
   }

   void Sparse :: resize( int _m, int _n )
   {
      m = _m;
      n = _n;

      data.clear();
   }

   int Sparse :: size( int dim ) const
   {
      if( dim == 1 ) return m;
      if( dim == 2 ) return n;
      return 0;
   }

   int Sparse :: length( void ) const
   {
      return max( m, n );
   }

   void Sparse :: zero( double rVal, double iVal )
   {
      EntryValue val( rVal, iVal );

      for( EntryMap::iterator i = data.begin(); i != data.end(); i++ )
      {
         i->second = val;
      }
   }

   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 );
   }

   void Sparse :: horzcat( const Sparse& A, const Sparse& B )
   {
      assert( A.m == B.m );
      assert( A.xtype == B.xtype );

      data.clear();
      m = A.m;
      n = A.n + B.n;
      xtype = A.xtype;

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

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

   void Sparse :: vertcat( const Sparse& A, const Sparse& B )
   {
      assert( A.n == B.n );
      assert( A.xtype == B.xtype );

      data.clear();
      m = A.m + B.m;
      n = A.n;
      xtype = A.xtype;

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

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

   cholmod_sparse* Sparse :: operator*( void )
   {
      if( A )
      {
         cholmod_l_free_sparse( &A, common );
         A = NULL;
      }

      int nzmax = data.size();
      int sorted = true;
      int packed = true;
      A = cholmod_l_allocate_sparse( m, n, nzmax, sorted, packed, stype, xtype, common );

      // build compressed matrix (note that EntryMap stores entries in column-major order)
      double* pr = (double*) A->x;
      double* pi = (double*) A->z;
      UF_long* ir = (UF_long*) A->i;
      UF_long* jc = (UF_long*) A->p;
      int i = 0;
      int j = -1;
      for( EntryMap::const_iterator e = data.begin(); e != data.end(); e++ )
      {
         int c = e->first.first;
         if( c != j )
         {
            for( int k = j+1; k <= c; k++ )
            {
               jc[k] = i;
            }
            j = c;
         }

         ir[i] = e->first.second;
         pr[i] = e->second.first;
         if( xtype == CHOLMOD_ZOMPLEX )
         {
            pi[i] = e->second.second;
         }
         i++;
      }
      for( int k = j+1; k <= n; k++ )
      {
         jc[k] = i;
      }

      return A;
   }

   double& Sparse :: operator()( int row, int col )       { return retrieveEntry( row, col, CHOLMOD_REAL    ); }
   double  Sparse :: operator()( int row, int col ) const { return retrieveEntry( row, col, CHOLMOD_REAL    ); }
   double& Sparse ::          r( int row, int col )       { return retrieveEntry( row, col, CHOLMOD_REAL    ); }
   double  Sparse ::          r( int row, int col ) const { return retrieveEntry( row, col, CHOLMOD_REAL    ); }
   double& Sparse ::          i( int row, int col )       { return retrieveEntry( row, col, CHOLMOD_ZOMPLEX ); }
   double  Sparse ::          i( int row, int col ) const { return retrieveEntry( row, col, CHOLMOD_ZOMPLEX ); }

   Sparse::iterator Sparse :: begin( void )
   {
      return data.begin();
   }

   Sparse::const_iterator Sparse :: begin( void ) const
   {
      return data.begin();
   }

   Sparse::iterator Sparse :: end( void )
   {
      return data.end();
   }

   Sparse::const_iterator Sparse :: end( void ) const
   {
      return data.end();
   }

   double& Sparse :: retrieveEntry( int row, int col, int c )
   {
      EntryIndex index( col, row );

      EntryMap::const_iterator entry = data.find( index );
      if( entry == data.end())
      {
         data[ index ] = EntryValue( 0., 0. );
      }

      if( c == CHOLMOD_REAL ) return data[ index ].first;
      else             return data[ index ].second;
   }

   double Sparse :: retrieveEntry( int row, int col, int c ) const
   {
      EntryIndex index( col, row );

      EntryMap::const_iterator entry = data.find( index );

      if( entry == data.end())
      {
         return 0;
      }

      if( c == CHOLMOD_REAL ) return entry->second.first;
      else                    return entry->second.second;
   }

   Dense :: Dense( Common& _common, int _m, int _n, int _xtype )
   : common( _common ),
     m( _m ),
     n( _n ),
     xtype( _xtype ),
     data( NULL ),
     rData( NULL ),
     iData( NULL )
   {
      int d = m; // leading dimension
      data = cholmod_l_allocate_dense( m, n, d, xtype, common );
      rData = (double*) data->x;
      if( xtype == CHOLMOD_ZOMPLEX )
      {
         iData = (double*) data->z;
      }
   }

   Dense :: Dense( const Dense& A )
   : common( A.common ),
     data( NULL ),
     rData( NULL ),
     iData( NULL )
   {
      *this = A;
   }

   Dense :: ~Dense( void )
   {
      if( data ) cholmod_l_free_dense( &data, common );
   }

   cholmod_dense* Dense :: operator*( void )
   {
      return data;
   }

   void Dense :: initializeFromCopy( void )
   {
      assert( data->dtype == CHOLMOD_DOUBLE );

      xtype = data->xtype;
      m = data->nrow;
      n = data->ncol;
      rData = (double*) data->x;
      if( data->xtype == CHOLMOD_ZOMPLEX )
      {
         iData = (double*) data->z;
      }
   }

   const Dense& Dense :: operator=( const Dense& A )
   {
      if( data ) { cholmod_l_free_dense( &data, common ); data = NULL; }

      data = cholmod_l_copy_dense( A.data, common );

      initializeFromCopy();

      return *this;
   }

   const Dense& Dense :: operator=( cholmod_dense* A )
   {
      if( data ) { cholmod_l_free_dense( &data, common ); data = NULL; }

      data = A;

      initializeFromCopy();

      return *this;
   }

   int Dense :: size( int dim ) const
   {
      if( dim == 1 ) return m;
      if( dim == 2 ) return n;
      return 0;
   }

   int Dense :: length( void ) const
   {
      return max( m, n );
   }

   void Dense :: zero( double rVal, double iVal )
   {
      if( rData )
      for( int i = 0; i < m*n; i++ )
      {
         rData[i] = rVal;
      }

      if( iData )
      for( int i = 0; i < m*n; i++ )
      {
         iData[i] = iVal;
      }
   }

   double Dense :: norm( void )
   // returns the infinity norm
   {
      return cholmod_l_norm_dense( **this, 0, common );
   }

   void Dense :: horzcat( const Dense& A, const Dense& B )
   {
      assert( A.m == B.m );
      assert( A.xtype == B.xtype );

      m = A.m;
      n = A.n + B.n;
      xtype = A.xtype;

      if( data ) { cholmod_l_free_dense( &data, common ); data = NULL; }
      data = cholmod_l_allocate_dense( m, n, m, xtype, common );
      rData = (double*) data->x;
      if( xtype == CHOLMOD_ZOMPLEX )
      {
         iData = (double*) data->z;
      }

      for( int i = 0; i < A.m; i++ )
      for( int j = 0; j < A.n; j++ )
      {
         rData[i+m*j] = A.rData[i+A.m*j];
         if( xtype == CHOLMOD_ZOMPLEX )
         {
            iData[i+m*j] = A.iData[i+A.m*j];
         }
      }

      for( int i = 0; i < B.m; i++ )
      for( int j = 0; j < B.n; j++ )
      {
         rData[i+m*(j+A.n)] = B.rData[i+B.m*j];
         if( xtype == CHOLMOD_ZOMPLEX )
         {
            iData[i+m*(j+A.n)] = B.iData[i+B.m*j];
         }
      }
   }

   void Dense :: vertcat( const Dense& A, const Dense& B )
   {
      assert( A.n == B.n );
      assert( A.xtype == B.xtype );

      m = A.m + B.m;
      n = A.n;
      xtype = A.xtype;

      if( data ) { cholmod_l_free_dense( &data, common ); data = NULL; }
      data = cholmod_l_allocate_dense( m, n, m, xtype, common );
      rData = (double*) data->x;
      if( xtype == CHOLMOD_ZOMPLEX )
      {
         iData = (double*) data->z;
      }

      for( int i = 0; i < A.m; i++ )
      for( int j = 0; j < A.n; j++ )
      {
         rData[i+m*j] = A.rData[i+A.m*j];
         if( xtype == CHOLMOD_ZOMPLEX )
         {
            iData[i+m*j] = A.iData[i+A.m*j];
         }
      }

      for( int i = 0; i < B.m; i++ )
      for( int j = 0; j < B.n; j++ )
      {
         rData[(i+A.m)+m*j] = B.rData[i+B.m*j];
         if( xtype == CHOLMOD_ZOMPLEX )
         {
            iData[(i+A.m)+m*j] = B.iData[i+B.m*j];
         }
      }
   }

   double& Dense :: operator()( int row, int col )       { return rData[ row + m*col ]; }
   double  Dense :: operator()( int row, int col ) const { return rData[ row + m*col ]; }
   double& Dense ::          r( int row, int col )       { return rData[ row + m*col ]; }
   double  Dense ::          r( int row, int col ) const { return rData[ row + m*col ]; }
   double& Dense ::          i( int row, int col )       { return iData[ row + m*col ]; }
   double  Dense ::          i( int row, int col ) const { return iData[ row + m*col ]; }
   double& Dense :: operator()( int index )              { return rData[ index ]; }
   double  Dense :: operator()( int index ) const        { return rData[ index ]; }
   double& Dense ::          r( int index )              { return rData[ index ]; }
   double  Dense ::          r( int index ) const        { return rData[ index ]; }
   double& Dense ::          i( int index )              { return iData[ index ]; }
   double  Dense ::          i( int index ) const        { return iData[ index ]; }

   Factor :: Factor( Common& _common )
   : common( _common ),
     L( NULL )
   {}

   Factor :: ~Factor( void )
   {
      if( L )
      {
         cholmod_l_free_factor( &L, common );
      }
   }

   void Factor :: build( Upper& A )
   {
      if( L )
      {
         cholmod_l_free_factor( &L, common );
         L = NULL;
      }

      L = cholmod_l_analyze( *A, common );
      cholmod_l_factorize( *A, L, common );
   }

   cholmod_factor* Factor :: operator*( void )
   {
      return L;
   }
}

• Rで検算

A <- matrix(c(2,-1,0,-1,3,-1,0,-1,2),byrow=TRUE,3,3)
A. <- t(A)
b <- c(1,2,3)
library(Matrix)
x <- Matrix::solve(A,b,sys="A")
x
A %*% x

> A <- matrix(c(2,-1,0,-1,3,-1,0,-1,2),byrow=TRUE,3,3)
> A. <- t(A)
> b <- c(1,2,3)
> library(Matrix)
> x <- Matrix::solve(A,b,sys="A")
> x
[1] 1.5 2.0 2.5
> A %*% x
     [,1]
[1,]    1
[2,]    2
[3,]    3