takapack.h (疎行列用 LU分解, CG法の素朴(遅い)な実装)


戻る topへ

要旨

  修士で大学院に入って以来、疎な連立方程式はumfpackが解いてくれて、本当に感謝している.
  企業とのコラボレーションのため,GPLのumfpackと一時お別れしなくてはならなくなった.
  
  疎な連立方程式を解くtakapack.h/takapack.cppを作ったので公開しておく.
  + 与える行列は, umfpackと同様の comressed row form (umfpackではcompressed column formだった)
  + CG法による  連立方程式 Ax = b の 求解をサポート
  + LU分解による 連立方程式 Ax = b の 求解をサポート
  
  NYSL(Version 0.9982)ライセンスなので、好きに使ってください(大学のレポートとかの参考になれば良いなと思っています.).
  


ソースコード

   takapack.c takapack.cpp


使い方

  1) 解きたい連立方程式 Ax = b の行列Aをcompressed row formで作る。
    配列bもつくる。

/*----------------------------
    2  3  0  0  0        8
    3  0  4  0  6       45
A =   0 -1 -3  2  0    b =-3
    0  0  1  0  0        3
    0  4  2  0  1       19
----------------------------*/

    //(umfpackは complessed columnだけど、実装の都合上 Compressed Rowを採用 ())
    // umfpackと併用するときは umfpack_solve の第一引数に UMFPACK_At を食わせて転置すればOK
    //compressed row form of A 
    const int N = 5;
    int    Ap [ ] = { 0,        2,          5,            8,   9,      12} ;
    int    Ai [ ] = { 0,  1,    0, 2, 4,   1,  2, 3,    2,   1, 2, 4   } ;
    double Ax [ ] = { 2., 3.,   3, 4, 6,  -1, -3, 2,    1,   4, 2, 1   } ;

    double b  [ ] = { 8, 45, -3, 3, 19} ;

  2) CG法で解くなら CG法の関数を呼ぶ
    double x1 [5] = {0,0,0,0,0};
    takapack_CG_sparse_solve(N, Ap, Ai, Ax, b, x2, 0.00000001);
    fprintf( stderr, "%f %f %f %f %f\n", x2[0], x2[1], x2[2], x2[3], x2[4]);}
  3)LU分解で解くなら i)LU分解関数 ii)解く関数 iii)解放する関数 を順に呼ぶ
    double x1 [5] = {0,0,0,0,0};
    
    int    *LUp, *LUi, *LU_rowflip;
    double *LUx;
    
    takapack_LU_factorization( N, Ap, Ai, Ax, LUp, LUi, LUx, LU_rowflip);     
    takapack_LU_solve        ( N, LUp, LUi, LUx, LU_rowflip, b, x1 ); 
    fprintf( stderr, "%f %f %f %f %f\n", x1[0], x1[1], x1[2], x1[3], x1[4]);//適当にprintfしてください...
    takapack_LU_free( LUp, LUi, LUx, LU_rowflip);