読者です 読者をやめる 読者になる 読者になる

くれなゐの雑記

身の回りの人や自分が困っていたことを記事にしています

rustをインストールから行列で連立方程式を解くところまで

Motivation

なぜrust?

環境

rustのインストール

このコマンドで簡単にインストールする

curl https://sh.rustup.rs -sSf | sh

以下のようなメッセージができる. 初心者なのでdefaultのstableをinsatllする.

   default host triple: x86_64-unknown-linux-gnu
     default toolchain: stable
  modify PATH variable: yes

1) Proceed with installation (default)
2) Customize installation
3) Cancel installation

インストールが終わると以下のメッセージが出てくるので

To configure your current shell run source $HOME/.cargo/env

.bashrcかなんかに

source $HOME/.cargo/env

を追加する

rustc --version

を入力して,それらしいものがでてきたらとりあえずOK.

References

www.rustup.rs

qiita.com

crateの作成

Rust では create (木箱) 単位でプログラムを開発するらしい
今回は簡単な行列計算をやってみるので,mat_calcという名前でcrateを作る

cargo new mat_calc --bin

なんかできた

output: 
Created binary (application) `mat_calc` project

ディレクトリを覗いてみると

$ ls mat_calc
.  ..  Cargo.toml  .git  .gitignore  src

どうやら.gitも生成されているらしいぱっと見た感じ.gitignoreもいい感じに定義されてた

実行テスト

src/main.rsにはすでにHello Worldを出力するプログラムが入っているのでとりあえず実行してみる とりあえずディレクトリに移動して

$ cd mat_calc
$ cargo run
   Compiling mat_calc v0.1.0 (file:///home/kurenaif/Desktop/mat_calc)
    Finished dev [unoptimized + debuginfo] target(s) in 0.53 secs
     Running `target/debug/mat_calc`
Hello, world!

無事出力に成功した

References

qiita.com

dependenciesの設定

行列を扱いたいけどFortranでもないしrustにはそんなライブラリは標準搭載されていない. rustで行列が扱えるライブラリを探したら以下の3つが見つかった

一番使いやすそうで,commitが盛んなnalgebraを使うことにした. 速度面はよくわからなかった.

ここで, Cargo.tomlを編集し,nalgebraを使用できるようにする(pip installに近いイメージ)

[dependencies]
nalgebra = "0.11.0"

versionは2017-05-05の時点での最新バージョンを使用した. githubなりなんなりをみてうまいこと指定する. githubのURLでの指定も可能だけど https://crates.io に登録されているのでこんな感じでURLなしでもcrateを入れることが出来る.

References

qiita.com

連立方程式を解いてみる

行列の表示

以下のソースコードsrc/main.rsに実装し,実行してみる 今後の流体計算を考慮して普通のMatrixよりも実行時に長さを決められるDMatrixのほうが都合がよさそうなのでDMatrixを使用する. 使い方その他もろもろはソースコードみて察していただきたい.
個々で気をつけなければならないのが行列の数字はすべて実数で表記しないとちゃんとprintlnしてくれない.
つまり 1,2,3… のような整数を書く場合でも, 1., 2., 3. と書いたほうが後々楽になる.
rustは整数と小数,また同じ整数どうしや小数どうしてもbit幅が違えばコンパイルエラーを出すほど厳密.

extern crate nalgebra as na;

fn main(){
    
    let mat = na::DMatrix::from_iterator(3,3,[
        1.,2.,3.,
        4.,5.,6.,
        7.,8.,9.
            ].iter().cloned());
    println!("{}",mat);

    let b = na::DVector::from_iterator(3,[
        1.,2.,3.
            ].iter().cloned());
    println!("{}",b);

}
$ cargo run
  ┌                   ┐
  │ 1.000 4.000 7.000 │
  │ 2.000 5.000 8.000 │
  │ 3.000 6.000 9.000 │
  └                   ┘

  ┌       ┐
  │ 1.000 │
  │ 2.000 │
  │ 3.000 │
  └       ┘
  • めっちゃいい感じの出力をしてくれる.(DMatrixじゃなくて普通のMatrixはこうはいかない)
  • どうやらソースコードに書いたものに比べて転置したものが出てくるらしい 縦ベクトルが3つ横に並んでるイメージ.

逆行列

方程式を解くためには逆行列が必要不可欠 さっき用意したmatは実は逆行列が存在しない.
面白そうなのでさっきのmat逆行列を計算してみる.
nalgebraでは,.try_inverse()を用いて逆行列を計算する.
無理だったら無理なりの挙動をしてくれる.

try_inverseを出力するプログラム.
println!("{}")println!("{:?}")では呼びされるメソッドが微妙に異なっており,後者の方がちょっと詳しいメッセージが出てくることが多い(デバッグ用出力?).
今回,前者を使うとエラーを吐くので後者を使っている.

extern crate nalgebra as na;

fn main(){
    
    let mat = na::DMatrix::from_iterator(3,3,[
        1.,2.,3.,
        4.,5.,6.,
        7.,8.,9.
            ].iter().cloned());
    println!("{:?}",mat.try_inverse());
}
$ cargo run
None

Noneが出力された.
リファレンスに書いているが,try_inverse()をしてなかった時はNoneが返ってくるらしい

では逆行列が存在する時はどんな出力をするのだろうか?
今回はtry_inverseの結果を"{}""{:?}"の二つを出力している.
最初のいもすさんの記事でも書いてある通り,rustでは基本所有権がmoveしちゃうので最初のprintln!matの所有権が切れる.
そこで,clone()することで所有権を継続させている.

extern crate nalgebra as na;

fn main(){
    
    let mat = na::DMatrix::from_iterator(3,3,[
        1.,3.,5.,
        2.,3.,8.,
        2.,1.,5.
            ].iter().cloned());
    println!("{:?}",mat.clone().try_inverse());
    println!("{}",mat.try_inverse().unwrap());
}
Some(Matrix { data: MatrixVec { data: [1.4, -2, 1.8, 1.2, -1, 0.4, -0.8, 1, -0.6], nrows: Dynamic { value: 3 }, ncols: Dynamic { value: 3 } }, _phantoms: PhantomData })
  ┌                      ┐
  │  1.400  1.200 -0.800 │
  │ -2.000 -1.000  1.000 │
  │  1.800  0.400 -0.600 │
  └                      ┘

一つ目の出力がtry_inverse()の中身を詳細に出力したものである.
Some()がついていて,実はこのままでは{}の出力ができない.
ここでSome()は値がNullになるかもしれない的なやつらしい.
Some()を取るためには,unwrap()をすると良い.

詳しくはここを参照する:

qiita.com

これの結果があってるかわからないので元のmatにかけて単位行列が出ること確認する
rustでは演算子オーバーロードができる.
このライブラリでは行列の掛け算もちゃんと定義されているので,以下のソースコードみたいな感じで表記できる.

extern crate nalgebra as na;

fn main(){
    
    let mat = na::DMatrix::from_iterator(3,3,[
        1.,3.,5.,
        2.,3.,8.,
        2.,1.,5.
            ].iter().cloned());
    println!("mat\n{}",mat.clone());
    println!("inv(mat)\n{}",mat.clone().try_inverse().unwrap());
    println!("mat*inv(mat)\n{}",mat.clone().try_inverse().unwrap()*mat);
}
$ cargo run
mat
  ┌                   ┐
  │ 1.000 2.000 2.000 │
  │ 3.000 3.000 1.000 │
  │ 5.000 8.000 5.000 │
  └                   ┘

inv(mat)
  ┌                      ┐
  │  1.400  1.200 -0.800 │
  │ -2.000 -1.000  1.000 │
  │  1.800  0.400 -0.600 │
  └                      ┘

mat*inv(mat)
  ┌                      ┐
  │  1.000 -0.000  0.000 │
  │  0.000  1.000  0.000 │
  │  0.000  0.000  1.000 │
  └                      ┘

ちゃんと逆行列が計算されているっぽい

連立方程式を解いてみる

いよいよ本題

とりあつかう問題はこれの例1: 連立方程式の解き方

連立方程式を $$ x = A^{-1} b$$ の形にしてinvして解くだけ

extern crate nalgebra as na;

fn main(){
    let a = na::DMatrix::from_iterator(3,3,[
        1.,1.,0.,
        2.,-1.,1.,
        3.,2.,1.
    ].iter().cloned());
    println!("A\n{}",a);

    let b = na::DVector::from_iterator(3,[
        0.,3.,-1.
    ].iter().cloned());
    println!("B\n{}",b);

    let invA = a.try_inverse().unwrap();
    println!("invA\n{}",invA);

    println!("x\n{}", invA*b)
}
A
  ┌                      ┐
  │  1.000  2.000  3.000 │
  │  1.000 -1.000  2.000 │
  │  0.000  1.000  1.000 │
  └                      ┘

B
  ┌        ┐
  │  0.000 │
  │  3.000 │
  │ -1.000 │
  └        ┘

invA
  ┌                      ┐
  │  1.500 -0.500 -3.500 │
  │  0.500 -0.500 -0.500 │
  │ -0.500  0.500  1.500 │
  └                      ┘

x
  ┌        ┐
  │  2.000 │
  │ -1.000 │
  │  0.000 │
  └        ┘

ページの答えと一致してるしちゃんと解けてるっぽい.

おまけ

今回色々やったリポジトリをここにおいておきます github.com

次回は簡単な熱拡散問題とか解いてみるかも