EigenからMAGMAに行列データを受け渡してGPUで連立一次方程式を解く

# はじめに
<!-- 発端や概要を記載 -->
GPUの性能が急速に伸びている中，数値計算界隈でもGPUでソルバーを解く機運が高まってきています．

しかし，cuda関連のプログラミングはレガシーな書き方が必要なため，自作はなかなか大変です．

そこで今回は，MAGMAと呼ばれるGPUによる線形代数ライブラリをインストールして連立一次方程式を解きたいと思います．
インストールだけだと，[公式](https://icl.cs.utk.edu/projectsfiles/magma/doxygen/installing.html)を見たほうが良いということになってしまうため，ここではEigen(CPUによる線形代数ライブラリ)からMAGMAへのデータ受け渡しについても書きたいと思います．

本記事で行う操作のイメージとしては，

Eigenで行列操作
↓
MAGMAに受け渡してAx=bを解く．

という感じです．

- makefileの編集ができること
- cudaがインストール済みであること

を前提に書いていきます．
私は，こちらの方の[サイト](https://rgb.sh/blog/magma)を参考にCUDAを入れました．

~~~
nvcc -V
~~~

と打ってバージョンが出れば問題ないと思います．

# 目次
<!-- タイトルとアンカー名を編集 -->
1. [MAGMAとは](#Chapter1)
1. [マシン情報](#Chapter2)
1. [インストール手順](#Chapter3)
1. [Eigen+MAGMAコード例](#Chapter4)
1. [所感](#Chapter5)
1. [参考文献](#reference)

<!-- 各チャプター -->
<a id="#Chapter1"></a>
# MAGMAとは
[参考](https://qiita.com/Hiroki11x/items/c751869e934d26a6b5af)


<a id="#Chapter2"></a>
# マシン情報
ディストリビューション: Ubuntu
バージョン: 18.04.5 LTS

CPU
Intel(R) Core(TM) i9-10900X CPU @ 3.70GHz
物理コア数：1
総スレッド数：２0

GPU
NVIDIA GeForce RTX 2080 Ti

CUDA version: 11.4
Driver version: 470.57.02

メモリ：　32GB

<a id="#Chapter3"></a>
# インストール手順

### openblasのインストール
MAGMAは，intel-MKL,openblas...等幅広いLAPACKに対応しています．自分は最初intelMKL(OneAPI)のライブラリを使おうとしましたが，OneAPIで入れたmklは対応しておらず（不可能ではないと思いますが），代わりにopenblasを用いることにしました．インストールは簡単で，

~~~
sudo apt install libopenblas-dev gfortran
~~~

だけです．

~10月7日追記~
ローカルに落としたい方は[こちら](https://qiita.com/JeJeNeNo/items/6c48401845c2e392a91a)が参考になります．

### MAGMAのダウンロード
2021年10月4日地点での最新版を落とします．

~~~
wget http://icl.utk.edu/projectsfiles/magma/downloads/magma-2.6.1.tar.gz
tar -xzf magma-2.6.1.tar.gz
~~~

### ビルド用のPATHをつなぐ
人によって違うと思いますが参考に．．．
~~~
export CUDADIR=/usr/local/cuda
export OPENBLASDIR=/usr/lib/x86_64-linux-gnu/openblas
~~~

### make
make.inc-examplesディレクトリの中に様々なパターンのmake.incファイルがあります．
ここでは，"make.inc.openblas"を使います．

~~~
cp make.inc-examples/make.inc.openblas make.inc
~~~

このmake.incを開くと設定ファイルが出てくるので，修正していきます．
まず，openblasのパスは以下の通り．

~~~
OPENBLASDIR ?= /usr/lib/x86_64-linux-gnu/openblas
~~~

"set our GPU targets"では，マシンによって決まるGPUのアーキテクチャを選びます．
2080Tiは"Turing"なので，

~~~
 GPU_TARGET = Turing
~~~

あとはmakeします．
自分は上記の操作で通りましたが，環境によってはmake.incをうまく編集する必要があります．

~~~
make
sudo mkdir /usr/local/magma
sudo chown $(whoami) /usr/local/magma
make install
sudo chown -R root /usr/local/magma
~~~

~10月7日追記~
ローカルに落としたい人は

~~~
make
make install prefix=<your path>/magma
~~~


magma/example/内でビルドがうまくいっているかテストすることができます．

<a id="#Chapter4"></a>
# Eigen+MAGMAコード例
以下のような"magmacg.cu"ファイル及び"Makefile"を作りました．
Eigen::SparseMatrixのouterIndexPtr(), innerIndexPtr(),valuePtr()メソッドで，ダイレクトに行列情報を送ることができるのが素晴らしいです．

~~~c++:magmacg.cu
#include <Eigen/Core>
#include <Eigen/Sparse>
#include <iostream>
#include <magma/include/magma_v2.h>
#include <magma/sparse/include/magmasparse.h>

using namespace Eigen;
using namespace std;

// solve Ax=b
int main()
{
  int m = 10;
  int n = 1;

  // diagonal matrix A
  SparseMatrix<double> AEigen(m, m);
  for(int i = 0; i < m; i++)
    AEigen.coeffRef(i, i) = 0.05;

  // x,b
  VectorXd rhs(m), sol(m);
  rhs.setOnes();
  sol.setZero();

  // Initialize MAGMA and create some LA structures.
  magma_init();
  magma_dopts opts;
  magma_queue_t queue;
  magma_queue_create(0, &queue);

  memset(&opts, 0, sizeof(magma_dopts));

  magma_d_matrix A = {Magma_CSR}, dA = {Magma_CSR};
  magma_d_matrix b = {Magma_CSR}, db = {Magma_CSR};
  magma_d_matrix x = {Magma_CSR}, dx = {Magma_CSR};

  // Pass the system to MAGMA.
  magma_dcsrset(m, m, AEigen.outerIndexPtr(), AEigen.innerIndexPtr(),
                AEigen.valuePtr(), &A, queue);
  magma_dvset(m, 1, rhs.data(), &b, queue);

  // Choose a solver, preconditioner, etc. - see documentation for options.
  opts.solver_par.solver = Magma_PCGMERGE;
  opts.solver_par.maxiter = 5000;
  opts.solver_par.rtol = 1e-10;
  opts.precond_par.solver = Magma_JACOBI;
  opts.precond_par.levels = 0;
  opts.precond_par.trisolver = Magma_CUSOLVE;

  // Initialize the solver.
  magma_dsolverinfo_init(&opts.solver_par, &opts.precond_par, queue);

  // Copy the system to the device
  magma_dmtransfer(A, &dA, Magma_CPU, Magma_DEV, queue);
  magma_dmtransfer(b, &db, Magma_CPU, Magma_DEV, queue);
  // initialize an initial guess for the iteration vector
  magma_dvinit(&dx, Magma_DEV, b.num_rows, b.num_cols, 0.0, queue);

  // Generate the preconditioner.
  magma_d_precondsetup(dA, db, &opts.solver_par, &opts.precond_par, queue);

  cudaDeviceSynchronize();

  // If we want to solve the problem, run:
  magma_d_solver(dA, db, &dx, &opts, queue);

  cudaDeviceSynchronize();

  // Then copy the solution back to the host...
  magma_dmtransfer(dx, &x, Magma_DEV, Magma_CPU, queue);
  // and back to the application code
  magma_dvcopy(x, &m, &n, sol.data(), queue);

  std::cout << "iterations: " << opts.solver_par.numiter << std::endl;
  std::cout << " res: " << opts.solver_par.iter_res << std::endl;

  // Free the allocated memory...
  magma_dmfree(&dx, queue);
  magma_dmfree(&db, queue);
  magma_dmfree(&dA, queue);
  magma_dmfree(&b, queue); // won't do anything as MAGMA does not own the data.
  magma_dmfree(&A, queue); // won't do anything as MAGMA does not own the data.
  // and finalize MAGMA.
  magma_queue_destroy(queue);
  magma_finalize();

  // output
  cout << sol << endl;
}
~~~


~~~shell:Makefile
#---------------------------------------------------------------
# User setting
#---------------------------------------------------------------
#path
MAGMADIR     ?= /usr/local/magma
CUDADIR      ?= /usr/local/cuda
OPENBLASDIR  ?= /usr/lib/x86_64-linux-gnu/openblas

#exe file
TARGET = magmacg.out

#set Eigen's path
INCLUDES = -I ../vendor/

#for magma
INCLUDES += -DADD_ -I$(MAGMADIR)/include -I$(MAGMADIR)/sparse/include -I$(CUDADIR)/include

#compiler
GCC = nvcc

#compiler option
CFLAGS = -std=c++14 -MMD -MP -O3 -DEIGEN_NO_DEBUG -DEIGEN_INITIALIZE_MATRICES_BY_ZERO

#link
LDFLAGS = -lm -L$(MAGMADIR)/lib -lmagma_sparse -lmagma -L$(CUDADIR)/lib64 -lcublas -lcudart -lcusparse -L$(OPENBLASDIR)/lib -lopenblas

#no make directory
NOMAKEDIR= .git% data% doc%

#object directory
OBJDIR = objs


#---------------------------------------------------------------
# Don't change the following
#---------------------------------------------------------------
#Run through the source directory with find command and list it down to the subdirectories.
SRCDIRS  := $(shell find . -type d)
#List all cpp files with the foreach command based on the SRCDIRS
CPPS   = $(foreach dir, $(SRCDIRS), $(wildcard $(dir)/*.cu))
#Remove directories not to make.
SRCS = $(filter-out $(NOMAKEDIR), $(CPPS))
#Get the name of the directory containing the SRCS.
DIRS = $(dir $(SRCS))
#Create a directory for binary with the same structure as the cpp files.
BINDIRS = $(addprefix $(OBJDIR)/, $(DIRS))
# Determine the object file name based on the SRCS
OBJS = $(addprefix $(OBJDIR)/, $(patsubst %.cu, %.o, $(SRCS)))
# Make the dependency file .d from the .o file
DEPS = $(OBJS:.o=.d)

default:
	@[ -d  $(OBJDIR)   ] || mkdir -p $(OBJDIR)
	@[ -d "$(BINDIRS)" ] || mkdir -p $(BINDIRS)
	@make all --no-print-directory
	@echo "make done!!"

all : $(OBJS) $(TARGET)

$(TARGET): $(OBJS) $(LIBS)
	$(GCC) -o $@ $^ $(LDFLAGS)

$(OBJDIR)/%.o: %.cu
	$(GCC) $(CFLAGS) $(INCLUDES) -o $@ -c $<

clean:
	@rm -rf $(TARGET) $(OBJDIR)

-include $(DEPS)
~~~

Eigenへのパスだけ編集すれば，makeできます．
大量の"warning"が出ますが動くのでよしとします．

~~~
make
./magmacg.out
~~~

対角行列なので反復計算はしません．

~~~
iterations: 1
 res: -0
20
20
20
20
20
20
20
20
20
20
~~~

<a id="#Chapter5"></a>
# 所感
これを自分が作っているFEMプログラムに差し込んでみました．

約２００万自由度のFEMの計算を行った結果です．

Eigen(CPU)のCG
![cpu.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1035656/0b20dad5-23df-c6d8-2454-214df1e623e4.png)


Magma(GPU)のCG

![gpu.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1035656/7625742f-dbea-1ad5-19db-5d741afb0a89.png)

精度に申し分はなく，計算時間は2/3程度になりました．
アセンブリングなどの時間込みなので，ソルバーの時間だけでみれば更に差がついていると考えられます．
多分行列のサイズを大きくすれば更に差がつくと思います．そのへんの検証はまた時間があるときにやります．
裏を返せば，１００万自由度超の計算でなければ恩恵は受けれないと思います．

<a id="#reference"></a>
# 参考文献
- [Seth Pendergrass氏のブログ](https://rgb.sh/blog/magma)
- [Matrix computations on the GPU, CUBLAS and MAGMA by example](https://www.researchgate.net/publication/280776423_Matrix_computations_on_the_GPU_CUBLAS_and_MAGMA_by_example)