Go to Qiita Advent Calendar 2024 Top
LoginSignup
6
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

@dwatanabee(ワタナベ ダイキ)

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

Last updated at Posted at 2021-10-04

はじめに

GPUの性能が急速に伸びている中，数値計算界隈でもGPUでソルバーを解く機運が高まってきています．

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

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

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

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

という感じです．

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

を前提に書いていきます．
私は，こちらの方のサイトを参考にCUDAを入れました．

nvcc -V

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

目次

  1. MAGMAとは
  2. マシン情報
  3. インストール手順
  4. Eigen+MAGMAコード例
  5. 所感
  6. 参考文献

MAGMAとは

参考

マシン情報

ディストリビューション: 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

インストール手順

openblasのインストール

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

sudo apt install libopenblas-dev gfortran

だけです．

~10月7日追記~
ローカルに落としたい方はこちらが参考になります．

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/内でビルドがうまくいっているかテストすることができます．

Eigen+MAGMAコード例

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

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

所感

これを自分が作っているFEMプログラムに差し込んでみました．

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

Eigen(CPU)のCG
cpu.png

Magma(GPU)のCG

gpu.png

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

参考文献

6
1
1

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
6
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?