スケジュールナースのブログ: 9月 2021

2021年9月30日木曜日

インストールの動画

正確には、大分前に、娘にアルバイトで作ってもらっていました。チュートリアル動画も発注してあったのですが、後期授業が始まってしまい、殆どできていません。基本的なところから、追々作っていきたいと思っています。こういう操作を知りたいというのがありましたら、お知らせください。時間はかかるかもしれませんが、追加していきたいと思います。

2021年9月29日水曜日

Windows11(Preview)にインストール

終夜運転用(世界記録更新狙い）のPCを自作しました。（パーツを集めて古いPCをリニューアル）新しいPCの動作確認も兼ね、Windows11での動作確認を行いました。

＜Windows10からの移行＞

PC正常性チェックアプリで確認したところ、最新CPU、最新マザーボードであるにも係わらず、移行できないとのメッセージ(TPM2.0..)が出ました。指摘に従って、BIOSのCPU設定を更新後にようやく、正常に移行できます、と出ました。これは、ちょっとしたバリアです。

＜Windows11にインストール＞

開くをクリックします。

Installをクリックします。

問題なく起動しました。

起動後は、以下のタスクバーから起動します。

確認に使用した筆者のハードウェアソフト環境です。

結論的には、特にWindows10と操作で変わる部分はないと思います。

2021年9月28日火曜日

疎行列密化コンパイル時間の削減

疎行列のSpMV演算時間本体は首尾よく短くなりましたが、今度は、密化するためのコンパイル時間の方が目立ってきており、本末転倒状態になってしまいました。そこで、コンパイル時間の削減に取り組みました。改善方法としては、ローカルキャッシュ的手法とマルチスレッド化です。

＜ローカルキャッシュ＞密化するためには、bit_vectorをKeyWordとしてKeyWordにできるだけ近いWORDをグループ化すれば、それがSIMD演算の対象になりえます。どれだけ、SIMD演算に持ち込めたのかは、

対象項目数/元のNonZeroエレメント数

で評価できます。実験的には、50%程度で、30%程度になる場合もありました。アルゴリズムが上手く行ったかどうかは、この指標で評価できます。アルゴリズムは、SIMD対象WORD毎にKEYWORDに対する距離演算結果をソートして、TOP　Kを持ってきます。これをSIMD対象BITS(SSE4.2で4bits,AVX2で8bits)毎に行うと大変時間がかかってしまいます。そこで、距離評価は、毎回行う部分と全体を行う部分に分け、全体を行う部分は、10-20回に一回程度としました。( この頻度を多くすれば、精度的には、理想に近づきますが、時間がかかります。逆に頻度を下げれば、コンパイル時間は短くなりますが、精度は低下します。) 毎回行う部分については、全体は見ずに、全体の上澄み部分だけを対象にします。いわば、ローカルキャッシュ手法です。

次に、マルチスレッド化です。対象行列は、二つあり、転置しない行列と転置した行列です。これらは独立なので、並列に作業できます。さらに、距離演算は、二つの評価項目があり、これも独立なので、並列に作業可能です。OpenMPで、異なる場所でのネストしたスレッドは、上手く記述できなかったので、std::threadで記述しました。

void compile() {
			 
			std::thread tha(&constraint_matrix_class::Tcompile, this);//転置行列
			std::thread thb(&constraint_matrix_class::nonTcompile, this);//転置しない行列
			
			
			tha.join();
			thb.join();
		}

上の転置行列中、さらに時間がかかる部分を並列化します。時間評価は、QueryPerformanceCounterで行い、粒度が細かすぎてスレッド起動の方が時間がかかる場合は、マルチスレッド化せずシーケンシャル駆動を行います。なお、threadの引数は、thisポインタ、const 参照にcref 、const参照でない参照にrefを使います。

void constraint_matrix_class::make_compiled_map(map>& hamming_map, map>& inclusion_map, const vector>& bitset_vec, const  set& entried, int prime_col)
{
    LARGE_INTEGER End1, End2;

    QueryPerformanceCounter(&End1);
    
    if (Tconstraint_matrix.cols() - entried.size() < 800) {
        make_hamming_map(hamming_map, bitset_vec, entried, prime_col);
        make_inclusion_map(inclusion_map, bitset_vec, entried, prime_col);
    } else {

        std::thread tha(&constraint_matrix_class::make_hamming_map, this, std::ref(hamming_map), std::cref(bitset_vec), std::cref(entried), prime_col);
        std::thread thb(&constraint_matrix_class::make_inclusion_map, this, std::ref(inclusion_map), std::cref(bitset_vec), std::cref(entried), prime_col);
        tha.join();
        thb.join();
    }
   
    QueryPerformanceCounter(&End2);
    Calc_sum4 += (double)(End2.QuadPart - End1.QuadPart) / freq.QuadPart;
}

2021年9月26日日曜日

行列の積分解

行列のリスケーリングは、効果的でした。そこで、これを前述のSpMV実装に適用したいのですが、一つ問題がありました。

\begin{align}(ABC)^{ \mathrm{ t } }=C^{ \mathrm{ t } }B^{ \mathrm{ t } }A^{ \mathrm{ t } } \end{align}

A,Cは、スケーリングの為の対角行列なので

\begin{align}(ABC)^{ \mathrm{ t } }=CB^{ \mathrm{ t } }A \end{align}

行列の積の結合法則により、

\begin{align}(ABC)^{ \mathrm{ t } }=C(B^{ \mathrm{ t } }A) \end{align}

\begin{align}\begin{pmatrix}y\end{pmatrix}=(ABC)^{ \mathrm{ t } }\begin{pmatrix}x\end{pmatrix}=C(B^{ \mathrm{ t } })\begin{pmatrix}Ax\end{pmatrix}\end{align}

と計算することが出来ます。ナーススケジューリングの場合、Bの殆どの係数は１です。そのままスケーリングしてしまうと、SpMVの積和演算をそのまま実行する必要がありますが、(3）のように計算すれば、その殆どを加算として計算できます。昨今のAVX2では、積和演算命令があり、加算のLatencyと変わらないし、むしろ積を分解することにより積回数は、増えてしまいます。その意味では得する部分はありません。しかし、行列の積係数をメモリしないでよいというアドバンテージが巨大インスタンスについては効いてくるはずです。

2021年9月22日水曜日

movemaskの逆命令はないの？

AVX512ならありますが、AVX2まではありません。AVX2では、速い代替手段がありますが、AVX以下なら、その手段も使えないので、いくつかの命令をシリーズで使うことになります。AVXで、折角VmaskLoad命令を備えたのにoperandが、__m128しかないおかげで余計なインストラクションを使わないといけません。（総じて、AVXのアーキテクチャは、急ごしらえの感が否めません。）そこで、FirstOrderサポート(超大なインスタンス対応、例えば、6ヶ月100人スタッフ）は、AVX2以降(2013年発売のHaswellアーキテクチャ以降）とすることにしました。机上計算では、これでEigenに対して、同じスレッド数下で、4倍以上のゲインを有するはずです。

2021年9月21日火曜日

Visual Studio のLoop　Unrollingは、制御不可

ループアンローリングは、勝手に入ったり入らなかったりします。コンパイルオプション、#pragma を試してみましたが、試した限り効力はありませんでした。

1重ループでは、何をやってもループアンローリングしてしまいます。しかし、2重ループ(run_code1A)では、ループアンローリングは入りませんでした。

while (address_cnt < n) {
			int col = address_ptr[address_cnt];
			float f = *(fp + col);
			__m128 B = _mm_load_ss(fp + address_ptr[address_cnt]);
			sum = _mm_add_ps(sum, B);
			++address_cnt;
		}

自前のループアンローリングは、次のように記述しています。

		}//Unrolling
#define USE_MY_UNROLLING
#ifdef USE_MY_UNROLLING
		while (address_cnt < n - 4) {

			__m128 B = _mm_load_ss(fp + address_ptr[address_cnt]);
			__m128 C = _mm_load_ss(fp + address_ptr[address_cnt + 1]);
			__m128 D = _mm_load_ss(fp + address_ptr[address_cnt + 2]);
			__m128 E = _mm_load_ss(fp + address_ptr[address_cnt + 3]);
			B = _mm_add_ps(B, C);
			D = _mm_add_ps(D, E);
			B = _mm_add_ps(B, D);
			sum = _mm_add_ps(sum, B);

			address_cnt += 4;
		}
#endif
		while (address_cnt < n) {
			int col = address_ptr[address_cnt];
			float f = *(fp + col);
			__m128 B = _mm_load_ss(fp + address_ptr[address_cnt]);
			sum = _mm_add_ps(sum, B);
			++address_cnt;
		}

ベクトル化手法を取らなくても、上の記述だけで、Eigenより速くなりました。(

Eigenは、SparseSpMVについて明示的アンローリングは行っていません。）

(instance19:マルチスレッド数4、FirstOrderIteration10万回　SpMVあたりのデータ)

■Eigen 　:3.96sec

■上記コード(ループアンローリングなし）:3.42sec

■ (ループアンローリングあり）:2.996sec

これをマルチスレッド化するには、OpenMPの ScheduleをDynamicにして、さらにchunk数を指定すると上手く動作しました。Eigenの記述をそのまま使いました。以上、ベクトル化前でも、ループアンローリング記述を追加するだけで、10-20%程度は、速度アップを期待できます。また、下のように2重ループで組んでも、データ的には、連続的に一重ループとして記述することが、重要なようです。キャッシュに収まり易いこともありますが、生成されるコードもシンプルになります。今回のようなクリティカルな関数では、一つのアセンブリコードの増減が速度に直結します。

void run_code1A(row_relay_struct* rrs, int rows, unsigned short* address_ptr, unsigned short* coeff_address_ptr, __m128* ifptr, float* ofptr, float* coeff_ptr) {


	/*LUT				 __m128
		rrs rows		 row_relay_struct rows

		one_address_ptr  unsigned short k
		address_ptr	     unsigned short n
		bits_ptr         unsigned char  n
		coff_address_ptr        unsigned short m
		coff_ptr				float          m




	*/
	int threads = Eigen::nbThreads();

#pragma omp parallel for schedule(dynamic,(rows+threads*4-1)/(threads*4)) num_threads(threads)
	
	for (auto row = 0;row < rows;++row) {
	
		__m128 sum = { 0,0,0,0 };



		int n = rrs[row].n;
		//assert(one_address_ptr + n == address_ptr);

		int address_cnt = rrs[row].address_cnt;
		if (row) {
			assert(rrs[row - 1].n == address_cnt);
		}
		float* fp = reinterpret_cast(ifptr);

		//Unrolling
#define USE_MY_UNROLLING
#ifdef USE_MY_UNROLLING
		while (address_cnt < n - 4) {

			__m128 B = _mm_load_ss(fp + address_ptr[address_cnt]);
			__m128 C = _mm_load_ss(fp + address_ptr[address_cnt + 1]);
			__m128 D = _mm_load_ss(fp + address_ptr[address_cnt + 2]);
			__m128 E = _mm_load_ss(fp + address_ptr[address_cnt + 3]);
			B = _mm_add_ps(B, C);
			D = _mm_add_ps(D, E);
			B = _mm_add_ps(B, D);
			sum = _mm_add_ps(sum, B);

			address_cnt += 4;
		}
#endif
		while (address_cnt < n) {
			int col = address_ptr[address_cnt];
			float f = *(fp + col);
			__m128 B = _mm_load_ss(fp + address_ptr[address_cnt]);
			sum = _mm_add_ps(sum, B);
			++address_cnt;
		}
		n = rrs[row].m;
		address_cnt = rrs[row].coff_cnt;
		if (row) {
			assert(rrs[row - 1].m == address_cnt);
		}
		while (address_cnt < n) {
			unsigned short col = coeff_address_ptr[address_cnt];

			float* fp = reinterpret_cast(ifptr);
			__m128 c = _mm_load_ss(fp + col);
			__m128 d = _mm_load_ss(coeff_ptr + address_cnt);
			c = _mm_mul_ss(c, d);
			sum = _mm_add_ps(sum, c);
			++address_cnt;
		}
		_mm_store_ss(ofptr + row, sum);
	}
	
}

まとめ１）VisualStudioは、ループアンローリングは、制御できない２）ループアンローリング効果は、10%-20%程度３）OpenMp for threadingは、staticよりもchunk付きdynamic scheduling が速い４）2重ループでもデータは、シーケンシャルで1重ループに。

2021年9月16日木曜日

疎行列の密化(SpMVのSIMD化）

部分的にでも密にできるところは、密にして並列度を高めよう、という発想です。サーベイした限り、この方法は見つかりませんでした。まずは、密化前の疎行列です。(instance19)

これを転置します。

y[m,1]=A[m,n]*X[n,1] 転置前

yt[1,m]=Xt[1,n]*At[n,m] 転置後

これを密化すると、つぎのようになります。

こうすると、SIMD化による性能向上が期待できます。AVXでは、単精度256bit,8個の並列演算まで可能です。

このインスタンスの場合、NonZero要素数が8万に対して、SIMD化しなければMAC演算数はそのまま8万、SIMD化によりMAC演算数は2万になるので、25%に演算数を削減できます。

（実際は、Load時Mask処理と後で述べる水平加算処理がプラスされます。）

2021年9月15日水曜日

疎行列の可視化

準備として、dynamic_bitsetをありがたく使用しました。c++ std::bitsetは、コンパイル時固定長でなければならず、面倒です。boostにもdynamic_bitsetはありますが、巨大なライブラリ郡をインストールするのも面倒なので、こちらを使うことにしました。

これで、疎行列部分のbitを対応させてbitmapにマッピングすると以下のような疎行列のbitmap(21MB)が生成されます。(instance19)

bitmapの生成には、こちらを使いました

一見、スパースのように見えないのですが、拡大すると、ちゃんとスパースに

なっています。

また、instance1は、下記です。画像のスケールが違うので、大きさの比較はできません。

しかし、instance1(7KB)の方が遥かに蜜(dense)行列であることは、イメージできます。

インスタンスが大きくなると、蜜度合いが薄くなる(疎度合いが大きくなる）ような傾向が感じられます。このように可視化すると改善の糸口が見えてきます。

2021年9月13日月曜日

Eigenスパース行列、既に最適化されていた

C++で置き換えてみたのですが、圧倒的に元のEigenの方が速かったです。それで、Eigen核心の部分のコードをDisassembleして見たのですが、次のように僅か5命令で、疎行列とベクタをロードし積和演算を行うことが達成されていました。見事なコードです。ソースはC++であり、どのように書いたらこのようなコードが生成されるのか、謎ですが、真似て書いても、中々、下のようなコードは生成されませんでした。(少なくとも二重ベクタループでは駄目でした。)

C++のvectorを使った二重ループの速度比較 - Qiita

00007FF6340034A0  movsxd      rax,dword ptr [r11+rcx*4]  
00007FF6340034A4  movsd       xmm0,mmword ptr [r9+rax*8]  
00007FF6340034AA  mulsd       xmm0,mmword ptr [r10+rcx*8]  
00007FF6340034B0  inc         rcx  
00007FF6340034B3  addsd       xmm1,xmm0  
00007FF6340034B7  cmp         rcx,rdx  
00007FF634003ま

また、マルチスレッド化は、OpenMPを用いた記述になっており、推定要素数に応じてマルチスレッドが駆動されます。実際、今回のインスタンスについては、ほぼ、スレッド数倍速になっており、既に大変よく効いていることが分かりました。

static void run(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const typename Res::Scalar& alpha)
  {
    LhsEval lhsEval(lhs);
    
    Index n = lhs.outerSize();
#ifdef EIGEN_HAS_OPENMP
    Eigen::initParallel();
    Index threads = Eigen::nbThreads();
#endif
    
    for(Index c=0; c1 && lhsEval.nonZerosEstimate() > 20000)
      {
        #pragma omp parallel for schedule(dynamic,(n+threads*4-1)/(threads*4)) num_threads(threads)
        for(Index i=0; i  
  
  
 疎行列とベクタの積は、SpMVというみたいです。いくつか論文がヒットし、実装を見てみました。 
  
[PDF] Autotuning Runtime Specialization For Sparse Matrix-Vector Multiplication by Buse Yilmaz, Baris Aktemur, M. Garzarán, Samuel N. Kamin, Mustafa Furkan Kıraç · 10.1145/2851500 (citationsy.com)
A sparse matrix‐vector multiplication method with low preprocessing cost - Aktemur - 2018 - Concurrency and Computation: Practice and Experience - Wiley Online Library
動的な並列実行機構を用いたSpMV実装の性能評価 | 文献情報 | J-GLOBAL 科学技術総合リンクセンター (jst.go.jp)
疎行列の問題は、疎であるが故に、メモリからのロードがネックになり最近のCPUの演算能力を生かせないことにあります。上記のような試みはなされていますが、本質的には、解決が難しいということです。Eigenのコードにしても、残る改善の余地は、ループアンローリングしかありません。しかし、これをC++で記述しても、Eigen程の精緻なコードを生成しないので、アセンブラで記述するしかないでしょう。
アンローリング｜行列積高速化#13｜やまもと｜note
アンローリングは、行うことにしても必要な性能(IPMの3倍速)は、得られそうにありません。そこで、さらなる改善としてNSPに特化した行列演算を試みることにしました。SIMDの命令セットを見たのですが、残念ながら効力があるだろうと思われるのは、AVXという命令セット以降です。最近のCPUならば、備えているとは思いますが、以降の検討は、AVXにすることにしました。（巨大インスタンスについては、それ以外に立ち向かう方法を思いつきません。）
記述は、アセンブラで行いますが、VisualStudioのアセンブラの使い方が良く分からず、JITアセンブラXbyakを使います。試してみたのですが、良い感じで使えそうでした。
また、現在の開発環境は古いCPUでAVXの評価が出来ないので、中古サーバを手配し開発環境を更新します。
にても、大規模NSP問題を解決するために、FirstOrderSolverを導入、FirstOrderSolverの問題を解決するために、疎行列演算問題を解く羽目に。疎行列問題を解決するために、SIMD
に手を出すことになるとは思ってもみませんでした。

2021年9月11日土曜日

たかが、FirstOrderMethod(一次法)　だから難しい　

untitled (orsj.or.jp)

大変参考になりました。合点のいく部分が多いです。

大元は、最急降下法なんですね。

「一次法は，二次法ほど汎用的な手法ではなく，ピタッとハマった問題に対してしか，よい性能を発揮できない．」また、「一次法は，数理最適化の知識だけでなく，モデリング，数値計算，アルゴリズム論，プログラミング技術を総動員しないと，よいものができない．」

もう少し頑張ってみます。

2021年9月10日金曜日

疎行列の計算

Eigenの動きを見ると,大規模インスタンスでは、勝手にマルチスレッドで動いています。しかしながら、バリアソルバと比べると未だバリアソルバの方が速いという結果です。予想外でした。

そもそもFirstOrderの利点としては、matrix-free (逆行列計算を必要としない）にあります。(Simplexにしてもバリアソルバにしても、疎行列の逆行列計算は、ほぼ必然といってよいと思います。)　利点部が逆にネックになっているというのは予想外でした。

良く調べていないのですが、ボトルネックは、規模から言ってベクトル演算部分ではなく、行列計算部であると思います。そこで、疎行列の積和演算に関しては、Eigenを使わずに自作することにしました。

まずは、疎行列演算規模の見積もりです。Branch＆Bound処理に行く直前の行列規模です。

nは、変数の数、mは、制約の数です。内ゼロでない係数が1.6-2.7%程度あるという結果になりました。やはり疎行列処理をすることが正しいという結果になりました。

instance21で、全体項数が625676　あり、行列xベクトル演算が2箇所、収束までIteration100000、全体での実行時間は、演算速度を1nsとすると

625676*2*100000*1e-9=125秒

となります。これがSingleThreadでの理論値になります。このままでは、バリアソルバでの実測値　75秒(Simplexでは数分かかるので論外）より大分劣ってしまいます。しかし、これから、SIMD化＋マルチスレッド化による性能向上は見込めるので、上手く実装できれば、バリアソルバの数倍、GPU実装では、さらなる効果が期待できます。(GPU実装は予定しません）

いずれにしても現在の状態を改善しない限り、1週間以上計算機を廻さないかぎり解が得られないので、高速化は必須です。

2021年9月7日火曜日

LP形式Conversion

しばしば、l <= a'x <= u (l !=u)というように両側で制約がある制約式に対して、スラックsを導入し、制約式に追加します。すると

`a'x - s = 0,

l <= s <= u

とConversionします。そうすると、両側制約式は、全て　RHS=0

と統一した行列で扱うことができます。スラック側に範囲制約が移っています。スラック変数の目的関数値への寄与係数ObjCoeffは、0でよいでしょう。目的関数値は変わりません。

2021年9月4日土曜日

Eigenで行列演算

C++テンプレートライブラリです。

でらうま倶楽部 : Eigen - C++で使える線形代数ライブラリ (livedoor.jp)

https://gist.github.com/AtsushiSakai/5227370

Numpy使いのためのEigenチートシート - ヤカンヒコウ (chachay.org)

FirstOrderSolverの記述で使います。記述してみると、確かにJulia並みに、行列演算が楽です。

FirstOrderは、http://www.orsj.or.jp/archive2/or64-6/or64_6_314.pdf

で述べられていますが、イメージプロセシング、機械学習の分野で成功を収めています。

一般的なOR問題をLpSolverとして扱う場合には、精度が不足し、収束が遅いという難点はあります。しかしながら、Simplexや、BarrierSolverで扱いにくい大規模な問題には、昨今の多コア、GPUとの親和性があるFirstOrderの方が歩があるのではないかと思います。今後、大規模LP問題解法として徐々にフォーカスがあたるのではないかと思います。

2021年9月3日金曜日

Barrier Solver Benchmarks updated

plato.asu.edu/ftp/lpbar.html

COPTが1位になりました。MDOPTも中国系です。BarrierSolverは、長くMOSEKが一位だったと思います。Highsが載っていないのですが、自分の評価では、Cplexよりやや落ちるくらいで、健闘していると思います。この表のなかでは、Highs>Tulip>CLPだと思います。

2021年9月2日木曜日

SAT Compretition 2021 Results

https://satcompetition.github.io/2021/results.html

BaseSolverが、Cadicalになって、CadicalのCベースバージョンKissatをベースにしたバージョンが上位を占めています。今年は、Cryptminisatは、参加していないようです。

2021年9月1日水曜日

大規模インスタンスで記録更新のためのPC性能

■CPU 3.3GHz→3.8GHz 20%弱速度UP

■スレッド数　４→16　2倍の速度アップ

■IPC　20%Up

なので、1.2*1.2*2=3倍

DDRのLatencyは、ほぼ変わらないのですが、キャッシュが効けば、マシンを置き換えるだけで、３倍の求解速度が期待できます。マシンを10日走らせるとき、3日で済むならば、投入する価値はあると考えます。メモリは、64GBを手配しました。

Ryzen 7 5800Xベンチマークレビュー：8コアCPUに6万円の価値はあるか？ | ちもろぐ (chimolog.co)

登録: 投稿 (Atom)

スケジュールナースのブログ