TAUCS を使う

2018年10月16日

はじめに

マトリックスソルバー TAUCS を使う。

バージョン

  • TAUCS Version 2.2
  • Ubuntu 16.04 LTS

TAUCS のコンパイル

ここ からパッケージを入手する。いくつかあるが、"Version 2.2 of the code, no external libraries, tgz format" でよい。

そのまま展開するとファイルが散らかるので、ディレクトリを作成してから展開する。

$ mkdir taucs
$ mv taucs.tgz taucs
$ cd taucs
$ tar xvzf taucs.tgz

基本的な操作は、configure, make のみ。

$ ./configure
cc -o configurator/configurator configurator/taucs_config.c
configurator/taucs_config.c: In function ‘emit_configfile’:
configurator/taucs_config.c:244:5: warning: implicit declaration of function ‘mkdir’ [-Wimplicit-function-declaration]
     mkdir(configdir);
     ^

警告が出ているが、気にしなくてよさそう。

$ make

cilk のエラーが出るが無視してよい。

コンパイルした環境では、いつくかエラーが出たので、以下の処置をとった。

パッケージのインストール。

$ sudo apt-get install libatlas-base-dev
$ sudo apt-get install libmetis-dev

config/linux.mk を編集。FC を修正。

FC        = gfortran #g77

CFLAGS がたくさんあるので、それらの最後に以下を追加 (それより上のものを上書き)。

CFLAGS    = -O3 -Wall

LIBF77 を修正。

LIBF77 = -lgfortran #-lg2c

これでコンパイルはできるが、テストプログラムでリンクして実行してみると、METIS のエラーらしきものが出る。

$ ./a.out
taucs_linsolve: preparing to factor
taucs_linsolve: ordering (llt=1, lu=0, ordering=metis)
 Runtime parameters:
   Objective type: Unknown!
   Coarsening type: Unknown!
   Initial partitioning type: Unknown!
   Refinement type: METIS_RTYPE_FM
   Perform a 2-hop matching: Yes
   Number of balancing constraints: 1
   Number of refinement iterations: 0
   Random number seed: 0
   Number of separators: 1094145072
   Compress graph prior to ordering: No
   Detect & order connected components separately: Yes
   Prunning factor for high degree vertices: 0.000000
   Allowed maximum load imbalance: 33.765

METIS のバージョンが新しいとダメらしい (Ver. 5 だとダメ?)。

METIS のコンパイル

TAUCS Version 2.2 は 2003 年にリリースされているので、その近辺にリリースされたバージョンの METIS を用いる。ここでは METIS 4.0.3 を用いた。こちらからダウンロードする。

$ tar xvzf metis-4.0.3.tar.gz
$ cd metis-4.0.3
$ make

TAUCS の再コンパイル

config/linux.mk を修正。コンパイルした METIS を使うようにする。

#LIBMETIS  = -L external/lib/linux -lmetis
LIBMETIS  = -L../metis-4.0.3 -lmetis

コンパイル。

$ make clean
$ make

インストール

インストールスクリプトはなさそうなので、以下のようなものを作った。

install

#!/bin/sh
INSTALL_DIR=../build
INCLUDE=$INSTALL_DIR/include
LIB=$INSTALL_DIR/lib

mkdir -p $INCLUDE
cp build/linux/*.h $INCLUDE
cp src/taucs.h $INCLUDE
cp src/taucs_private.h $INCLUDE

mkdir -p $LIB
cp lib/linux/*.a $LIB
cp ../metis-4.0.3/libmetis.a $LIB

コンパイルした METIS の静的ライブラリもコピーしている。

テスト

テストプログラムをコンパイル、実行してみる。作業ディレクトリを test とし、テストプログラムを test.c とする。

Makefile は、雑だが次のようにする。

test/Makefile

all:
        gcc test.c -I../build/include -L../build/lib -ltaucs -lf77blas -lcblas -latlas -llapack -lmetis-edf -lf2c -lgfortran -lm

コンパイル、実行。

$ cd test
$ make
$ ./a.out
taucs_linsolve: preparing to factor
taucs_linsolve: ordering (llt=1, lu=0, ordering=metis)
taucs_linsolve: ordering time 0.00e+00 seconds (0.00e+00 seconds CPU time)
taucs_linsolve: starting LLT factorization
taucs_linsolve: pre-factorization permuting of A
taucs_linsolve: starting IC LLT factorization
taucs_linsolve: starting IC LLT MF factorization
		Elimination tree depth is 3
		Symbolic Analysis of LL^T: 8.00e+00 nonzeros, 2.20e+01 flops, 2.85e+02 bytes in L
		Relaxed  Analysis of LL^T: 1.00e+01 nonzeros, 3.40e+01 flops, 3.33e+02 bytes in L
		Symbolic Analysis            =      0.000 seconds (0.000 cpu)
		Supernodal Multifrontal LL^T =      0.000 seconds (0.000 cpu)
taucs_linsolve: preparing to solve
taucs_linsolve: pre-solve permuting of A
x =
      1.000e+00
      3.000e+00
      4.000e+00
      2.000e+00

問題なし。

テストプログラム

test.c

/*
 * TAUCS Test
 *
 * [ 2 -1  0  0] [ ]   [-1]
 * [-1  3 -1  0] [x] = [ 4]
 * [ 0 -1  3 -1] [ ]   [ 7]
 * [ 0  0 -1  2] [ ]   [ 0]
 */
#include <stdio.h>

#define TAUCS_CORE_DOUBLE
#include <taucs.h>


int main(void)
{
	taucs_ccs_matrix *a;
	taucs_double b[] = {-1., 4., 7., 0.};
	taucs_double x[4];
	char *options[] = {"taucs.factor.LLT=true", "taucs.factor.mf=true", NULL};
	//char *options[] = {"taucs.factor.LLT=true", "taucs.solve.cg=true", NULL};
	int i;

	taucs_logfile("stdout");

	a = taucs_ccs_create(4, 4, 7, TAUCS_DOUBLE|TAUCS_SYMMETRIC|TAUCS_LOWER);

	a->colptr[0] = 0;
	a->colptr[1] = 2;
	a->colptr[2] = 4;
	a->colptr[3] = 6;
	a->colptr[4] = 7;

	a->rowind[0] = 0;
	a->rowind[1] = 1;
	a->rowind[2] = 1;
	a->rowind[3] = 2;
	a->rowind[4] = 2;
	a->rowind[5] = 3;
	a->rowind[6] = 3;

	a->taucs_values[0] =  2.;
	a->taucs_values[1] = -1.;
	a->taucs_values[2] =  3.;
	a->taucs_values[3] = -1.;
	a->taucs_values[4] =  3.;
	a->taucs_values[5] = -1.;
	a->taucs_values[6] =  2.;

	//taucs_ccs_write_ijv(a, "amat.txt");

	if(taucs_linsolve(a, NULL, 1, x, b, options, NULL) != TAUCS_SUCCESS){
		printf("solution failed\n");
		return 0;
	}

	printf("x =\n");
	for(i = 0; i < 4; i++) printf("%15.3e\n", x[i]);

	taucs_ccs_free(a);

	return 0;
}

CCS

TAUCS は、マトリックスの圧縮格納形式を扱うことができる。対称行列を CCS (Compressed Column Storage) で配列に格納する場合には、次のようになる。

a =

  [ 2 -1  0  0]
  [-1  3 -1  0]
  [ 0 -1  3 -1]
  [ 0  0 -1  2]

※インデックスは 0 から開始

value = [2, -1, 3, -1, 3, -1, 2]  非 0 要素
column = [0, 2, 4, 6, 7]          対角要素のインデックス。最後の値は非0要素の数
row = [0, 1, 1, 2, 2, 3, 3]       非 0 要素の行のインデックス

要素 a(i, j) へのアクセスは、次のようになる。

k: row(k) = i and column(j) <= k < column(j+1)
a(i, j): value(k)

例) a(2, 1) の場合

column(1) = 2, column(2) = 4
2 <= k < 4 and row(k) = 2 -> k = 3
a(2, 1): value(3) = -1

例) a(3, 3) の場合

column(3) = 6, colume(4) = 7
6 <= k < 7 and row(k) = 3 -> k = 6
a(3, 3): value(6) = 2