2021年10月31日日曜日

看護実習科目の年間スケジュール表の作成-GUI編

 前回、シフト連 成人1>成人2 を全スタッフについて制約するコアアルゴリズムをPythonで記述しました。今回は、任意のシフトの連順をGUIで操作できるように、GUIとPythonの記述追加を行います。

 下図は、シフト連の設定画面例です。スタッフ名1Gでは、成人1>成人2>老年>母性>地域>精神>小児>統合 の順でシフト連順を設定しています。その他のスタッフでは、成人1>成人2のみに設定する画面となっています。 

これにより、任意のスタッフに任意のシフト連順を設定することが出来ます。さて、スタッフプロパティでこのように設定するためには、スタッフプロパティ属性を記述しています。


つまり、

シフト連名=シフト名+'_'+連順(1,2,3,4,5,6,7,8)

としています。ルールです。このようにするのは、スタッフプロパティアイテム同士で、名前が衝突しないようにするためです。注意する点は、以下の2点です。

■シフト名と一致しないこと

■連順同士で区別できること

Python上で同じ名前になってしまうと意図しない動作となってしまい為です。(エラーは出ません。)

設定が多く面倒そうですが、テキストエディタで1行を記述しておいて、コピペ、置換、コピペ..で行うと、それほど面倒でもありません。



RUNした後、ソース全体を見ると、下のように成人1と、成人2には、全スタッフ[0,1,...]が設定されていることが分かります。また、老年は、スタッフ1Gである[0]だけが設定されていることが分かります。


よって、これらをPythonで拾ってくれば、GUIの設定を拾うことが出来ます。

<メンテナンス性を上げる>

Pythonソースを見ると特定のシフト名に対する記述がありません。今後のシフト名の変更、追加では、Pythonソースをいじらなくて良いようにしています。(プログラムは、書いたときは覚えていますが、1年後は、初見と同じです。)シフト名の変更追加があったときは、スタッフプロパティのみを変更すればOKです。覚えておくべきルールは、ルールだけです。

以下がPythonで書いた全ソースになります。



import sc3
import re #正規表現

def get_sequence_shift(person,seq):
    for s in globals():#全ての変数名をスキャンする
        rs=re.sub('_\d+', '', s)#正規表現 sから'_数字'を削除する ルール参照
        if s !=rs and rs in shiftdef:#削除した名前がシフト名と一致するものがあるなら
            ss=rs+'_'+str(seq)    
            if ss==s:#要求連順Noと等しいなら
                if person in eval(s):#文字列を評価して要求personと一致しているならば
                    #sc3.print('person='+str(person)+' '+rs+'\n')
                    return rs #シフト名を返す
    return ''#要求に一致するものはないことを報告する


def 順序制約(person,shift1,shift2):#person シフト連名 シフト連名+1の制約を行う。どちらも空白ではない

    first_day=今月[0]
    last_day=今月[len(今月)-1]

    st=staffdef[person]+' 順序制約 '+shift1+' > '+shift2 
    sc3.print(st+'  を制約します。\n')
    for day in 今月:
        vlist2=[]
        vlist1=[]
        for d in range(first_day,day+1):
            v2=sc3.GetShiftVar(person,d,shift2)
            vlist2.append(v2)
        for d in range(day+1,last_day+1):
            v1=sc3.GetShiftVar(person,d,shift1)
            vlist1.append(v1)
        V2=sc3.Or(vlist2)
        V1=sc3.Or(vlist1)
        st_day=st+' '+ str(day)
        sc3.AddHard(~(V2&V1),st_day) 
 
def 順序():
    n=len(shiftdef)
    for person in 全グループ:#全スタッフについて
        for seq in range(1,n+1):#全シフトについて
            shift1=get_sequence_shift(person,seq)# 連1
            shift2=get_sequence_shift(person,seq+1)#連2
            if len(shift1)>=1 and len(shift2)>=1:#どちらも空白でないときだけ
                順序制約(person,shift1,shift2)#制約する
 
順序()

上記ソースで、get_seq_shiftが、ルールに基づいたシフト名を返しているルーチンになります。ここが、一手にGUIとのインタフェースを引き受けています。正規表現とかglobals()とか、私もいちいち関数名を覚えているわけではなく、ググリながらコピペすることが多いです。

制約部については、前のコアアルゴリズムを改造しています。

以上です。

<短く書ける理由>

意外に短く書けているのは、Python自体のライブラリが充実していることと、既にGUIで記述した中身がPythonで記述されているからです。(ソース全体を見ると、結構長いソースになっています。)

<GUIは、まとめ制約にすぎない>

普段、なにげなく記述しているGUIも、裏では、今回記述したみたいに、展開すると沢山のprimitiveなAddHardや、AddSoftが動いていることが想像できると思います。Python上の制約も、Engineの中では、区別がありません。仮にC++で組み込んだとして比較しても性能上の優劣はありません。

その気になれば、殆どGUIを使わずにPython主体に組むことも可能ですが、GUIで設定仕切れない部分だけをPythonで記述することをお勧めします。


看護実習科目の年間スケジュール表の作成 - シフト連問題

 看護学部の実習先を割り当てる問題です。実習先を一つのシフトとして表現します。

シフト定義としては、

成人1

成人2

老年

地域

精神

小児

母性

統合

があります。期間60日内に上記シフトを割り当てます。シフトは、連続的であり、期間内に一回です。例えば、上の解では、成人1が3日、成人2が3日になっていますが、これらは、期間内に一回となっています。こういった連続シフト島をシフト連と呼称し、シフト連毎に順番の制約がある問題をシフト連問題と呼ぶことにします。

<シフト連問題は、何が問題か?>

シフト連毎に順番があることが問題です。看護師である妻に聞いたところ、たとえば、「成人1の前に、成人2は、有り得ない」、「統合は、まあ最後かな。 」という具合に、シフト連毎にハード制約とソフト制約が存在するようです。Day単位ならば、GUIでシフトパターンとして記述が可能ですが、不定長の連なので、GUIでは記述が困難です。


<難易度>

Pythonで制約を記述しますので、上級です。私も数十分悩んだのでアルゴリズムの難易度は高いです。加えて、GUIで操作インターフェースをどのようにPythonで記述するか?という問題もあります。やはり難易度が高いです。

<アルゴリズム>

便宜的にA>BをAは、Bに先行すると表記します。シフトA、BがあったときA>Bを実現するには、次のように、簡単な2日間から出発し、n日間では、どういう制約が必要になるかを考えてみます。


すると、期間をありえる全ての2つに分けて禁止パターンを記述していけばよいことに気づきます。ちなみに一つの連の長さは不定で、期間内に1個のみという仕様です。禁止パターンで記述しておけば、連が存在しないときでも悪さはしません。あくまで全てのB>Aを禁止するだけです。60日期間であれば、59個の禁止パターンを制約することになります。

<連の記述>

連は、期間内に1個なので、各Person,各DayでORを取れば十分です。

まずは、成人1>成人2 制約を記述してみます。

僅か24行です。

import sc3

def 順序():

    first_day=今月[0] #期間初日
    last_day=今月[len(今月)-1] #期間最終日
    for person in 全グループ:#全ての人について 
        st=staffdef[person]+' 順序制約'
        sc3.print(st+'します。\n')
        for day in 今月:#全期間あらゆる分割を行う
            vlist2=[]
            vlist1=[]
            for d in range(first_day,day+1):#初日からday日まで
                v2=sc3.GetShiftVar(person,d,'成人2')#person day shiftで状態取得 Binary 0/1
                vlist2.append(v2)
            for d in range(day+1,last_day+1):#day+1日から最終日まで
                v1=sc3.GetShiftVar(person,d,'成人1')#person day shiftで状態取得 Binary 0/1
                vlist1.append(v1)
            V2=sc3.Or(vlist2)#成人2連を求める
            V1=sc3.Or(vlist1)#成人1連を求める
            st_day=st+' '+ str(day)
            sc3.AddHard(~(V2&V1),st_day) #成人2>成人1 を禁止する。ハード制約する

順序()
これを打ち込んで、走らせてみます。最初の解で見られた、成人2>成人1 箇所は、無くなり全スタッフで、成人1>成人2 になりました。
成功です!

<Pythonを書くときの注意点>

■いつものように設定ボタンを押すのを忘れずに。また、前のソースに戻ってしまう現象が起こることがあります。RUNする前に保存も行ってください。

インデント自体に意味があります。インデントを揃えてください。下記は、インデントエラーの例です。赤いところをダブルクリックします。

すると、ソース全体に飛んでいきます。インデントが揃っていないことがわかります。
ソース全体は、ReadOnlyです。Editするには、ソースタブで元のソースを開いてください。


<一見宣言されていないような変数はどこにある?>
一回RUNして、ソース全体でご確認ください。例えば、”全グループ”は、ソース全体の最小の方で記述されています。配列になっていますが、これは、personを0から振った番号が収められています。for person in 全グループでは、person に0,1,2,3...とループします。動作が軽くなるように、person,day共、基数化しています。(最終的には、AddHardで、person,dayの数字を渡すので。)


<動作確認>
設計時は、頻繁にRUN-Editを繰り返すので、動作は軽い方が吉です。列制約をオフにして、ターゲット動作確認だけに注力するとよいでしょう。いきなりGUI部から記述するのではなく、上のように、アルゴリズム本質部の確認から行うのが吉です。

<汎用化>
成人1>成人2は、上記で確認出来ました。GUIで任意のシフト連の順序を記述できるようになると、より使い易いと思います。次回は、GUIインタフェースについて考えてみます。

2021年10月30日土曜日

Updated World Record on Instance21

Machine: E3-1226V3 3.3GHz

ScheduleNurse3 has found the Objective function value 21158. (LB=21133) Previously, it was LB=21124 UB=21159.





I reported it immediately, and it was approved.

schedulingbenchmarks.org


To check the objective function value, load the solution in AutoRoster, which will be displayed on the bottom right.  Please note "Feasible" is displayed.  It means the solution I submitted satisfies all constraints on Hard. 


The value in the lower-left corner of AutoRoster seems to indicate the size of the solution space, which is 10 to the 17367th power. I don't know how the software calculated, but I believe that even if we collect all the subatomic particles in the universe, they will not reach the 100th power of 10. It is also a remarkable coincidence that only one of the solutions was updated from 21159 to 22218 in the vast solution space.


This result is a prototype of the new Branching Algorithm. After we finish the basic algorithm, we expect to find the exact solution near LB.


It is still far from being ready for actual use, and there is still a lot of work to be done, so that it will take until the end of this year to update ScheduleNurse3.

Instance21世界記録更新

 マシン:E3-1226V3 3.3GHz

目的関数値21158を得ました。(LB=21133)

これまでは、LB=21124 UB=21159でした。



早速レポートして承認されました。

What's new? (schedulingbenchmarks.org)

目的関数値の確認は、解をAutoRosterで読み込ませて下右に表示されます。スケジュールナースのUBとAutoRosterのPenaltyが一致しているので間違いないでしょう。Feasibleというのが重要で、Not Feasibleが表示されるとハードエラーがあることを意味します。

AutoRoster左下の値は、解空間の広さを表しているようです。10の17367乗になっています。どういう計算でそういうなるのかは、謎ですが、確か全宇宙の素粒子を全部集めても10の100乗には届かなかったと思います。膨大な解空間中で、21159⇒22218へ1だけ更新というのもなんだか、凄い偶然です。

なお、本結果は、Branch Algorithmが未だ完成されたものではなく未だプロトタイプのものです。基本アルゴリズム完成後は、恐らく、LB近くで厳密解が発見できるものと期待しています。

実使用には程遠く、未だ多くの作業が残っているので年内一杯実リリースはかかりと思います。


2021年10月29日金曜日

Instance20速度改善

マシン:E3-1226V3 3.3GHz

改善前:23500sec

改善後:2023sec

今回の改善ポイントは、Instance13に同じですが、MyIPMSolverをVersionUpしています。

Benchmarks | What is Schedule NurseⅢ (nurse-scheduling-software.com)

によると、Cplexでは、Feasibleな解が8時間以内には得られず、Gurobiでは、6371sec となっています。NEOS SERVERとの違いはありますが、およそState Of Art Solver Gurobiに対して3倍の速度となっています。現リリースバージョンからは、10倍以上の速度アップとなります。



2021年10月28日木曜日

Instance13速度改善Ⅳ

 マシン:E3-1226V3 3.3GHz

改善前:1409sec

改善後:781sec

今回の改善ポイントは、

■MyIPMSolver使用

■Branch Algorithm改善 (Multithreaded Branching廃止)

■Algorithm1適用のタイミング改善

です。特に大きく効いたのは、Branch Algorithmです。

Linear SolverがSingle Threadで動くときは、この計算時間がボトルネックですので、マルチスレッドブランチングは意味がありました。が、Linear Solver自体をマルチスレッド化すると、他に廻せるリソース余裕はありません。マルチスレッド化する意味がないことになります。また、マルチスレッド化ブランチングは実装上難しいという問題も抱えていました。この際廃止することにしました。

これだけ速いとCplex/Gurobiと比べても大分速いと思います。少し古いデータですが、ルーベンカトリック大学のPieter Smetさんが計測結果では、12000秒となっています。 What's new? (schedulingbenchmarks.org)

今は、もう少し速くなっていると思いますので、後で再計測してみます。



2021年10月22日金曜日

MyIPMSolverの評価

 内点法による実装ソルバの評価を行いました。名前をMyIPMSolverとしました。

インスタンスは、SchedulingBenchmarksとINRCⅡから、規模の大きなインスタンスです。比較対象としては、オープンソースのCLP,HIGHS,商用ソルバ郡は、NEOS SERVERによる結果です。

結果としては、全てのインスタンスについて1位(王冠)となりました。黄色は2位オレンジは最下位です。また、参考として、最後にRyzen5800でのベンチマークを取っています。


グラフでは、背が低い程高性能を意味します。

<商用ソルバとの比較>

ただし、一般にIPMSolverは、IPMが結果が出た後、さらにSimplexのBasis形式に変換(Crossoverと呼ばれる)します。ユーザは、BranchAndBoundを行う場合に前の結果を使って(近い場所にある)演算を行います(WarmStart)。IPMの場合、殆ど場合、WarmStartが出来ない、もしくは行ったとしても効果がないことがあって、WarmStartが出来るSimplexSolverでResolveすることが多いのです。そのためには、Simplexで計算が出来るポイント解(Basis)を出力することが必要となります。

インスタンスによりますが、内点法で食う時間よりもCrossoverにかかる時間の方が大きい場合もあります。特に、巨大なインスタンスではその傾向が顕著になります。よって、上記結果は、MyIPMSolverがCrossoverを備えていませんのでフェアな比較になっていません。また、当然のことながら、NEOS SERVERとの比較ですので、なにがしかのハンデがあるのかもしれません。

<何故COPTが悪い成績なのか?>

商用IPMソルバは、インスタンスを見てSimplex/IPMを切り替えているのが普通です。大きくないインスタンスはSimplex,大きいインスタンスはIPMを使います。例えば、CPLEXのIPMがオンになるのは、Instance20から上です。ところがNEOS上にあるCOPTは、Simplex Onlyの模様です。結果、巨大インスタンスが来ると悲惨な結果となってしまいます。

<CROSSOVERはなくてもよいのか?>

確かにBranch&Boundを行う場合、Simplexの方が良さそうな気がするのですが、NSP巨大インスタンスの場合、少し離れただけで、SimplexではとてもReasonableな時間内に解が出てきませんでした。よってNSPに限れば、Primal/Dual解が得られれば、Simplexを行う必要はなく、そのままIPMを使い続けた方が得策と考えています。(FirstOrderと違って十分な精度1e-8がIPMでは得られているので、数値的な問題は考えにくい。)

<オープンソースソルバCLP,HIGHSとの比較>

Crossoverを引いたとしても明らかな優勢を見てとれます。CLP/HIGHSを使うと、Instance21以上で解を得るには、一週間以上計算機を廻さないといけないません。これらのInstance郡については、今後MyIPMSolverを使用して世界記録更新を狙います。また、INRCⅡの巨大インスタンスについても、今まで3日以上かかっていました。このソルバを使えば、Reasonableな時間内に解を示せるようになると思います。(n120w8インスタンスRyzen5800比で見ると今までより約30倍速) つまり、今までAlgorithm4で、実用的な時間に解を示せなかった巨大インスタンスも示せるようになる、ということです。

<E31226 4コア3.3GHzとRyzen 5800X 8コア3.8GHzの比較>

ターボブーストの影響もあり正確な性能比較ではないのですが、ほぼ1/2の時間で済んでいます。クロックUp・IPCの効果もありますが、やはりコア数が倍になっている影響が大きいと思います。ちなみにハイパースレッドをONにして論理コア数16にしてやってみたのですが、予想通り、改善されるどころか悪化してしまいました。これは、前に示したようにOneコアあたりのCPUパイプラインをほぼ埋めるようなコーディングをました。他に廻すCPUリソース余裕がない為、と推測されます。(なお、Windows11は、RyzenL3キャッシュ問題のFix済みVersionを使用しました。)

ともあれ、まがりなりにも、マルチスレッドの効果はある、と言えそうです。擬似論理コアベースではなく、実論理コアベースで。


2021年10月21日木曜日

コレスキー分解のマルチスレッド化 Intel OneTBBの導入

 修正コレスキー分解を普通にマルチスレッド化すると意外なほど性能がでませんでした。

恐らくキャッシングの問題だと思われます。この問題に対して、Lapackの実装は、キャッシュ階層を考慮した設計になっています。

また、EigenのLinearSolver郡も試しましたが、Dense行列ではLLTならマルチスレッド化されて速いです。(LDLTは、未だのようです。)しかしSparse行列を食わせるとシングルスレッドのままで、性能が出ませんでした。メモリも多々食います。

ここでのジレンマは、小行列にするとキャッシュ効果は出るけれども、マルチスレッド化時の性能は出にくくなる、というもので、奥が深いです。可能ならば、Lapackを持ってくるのが正解でしょう。


記述が再帰ですので、綺麗にスレッドに分解できません。そこで、スレッドよりも粒度が細かいタスクを導入するためにIntelOneTBBを導入することにしました。

本当は、OpenMPの omp taskがVisual Studioで使えればよいのですが、未だ対応していないので苦肉の策です。(こちらで計測の通り、OpenMPの方が性能が高いようです。)

Intel OneTBBの導入方法

1)こちらからダウンロード展開

2)CMAKEでソリューション・プロジェクトファイル郡を生成

3)タスクの導入用には、tbbプロジェクトだけで十分

4)モジュール定義が”Visual”で、何故か区切れてしまっているので、コンパイルエラーになってしまいます。defファイルの場所を指定してやるとエラーなく動作しました。また、リンカーのコマンドラインの追加部分にも残骸が残っているので削除します。

5)debug/release用に各々DLLファイルをコンパイル生成


oneapi::tbbのタスクの記述例

if ( nThis >= NSIZEH ){//シリアル実行の方が速い場合はマルチスレッド化しない
        oneapi::tbb::task_group tg;

        tg.run([thisStruct, aTri, nThis, aUnder, diagonal, work, nLeft2, iBlock, jBlock, numberBlocks,level]() {
            CholeskyTAK(thisStruct, aTri, nThis, aUnder, diagonal, work, nLeft2, iBlock, jBlock, numberBlocks, level+1);
         });
        tg.run([thisStruct, aTri, nThis, aUnder ,  diagonal, work, nLeft , nLeft2, iBlock, nb, jBlock, numberBlocks,level]() {
            CholeskyTAK(thisStruct, aTri, nThis, aUnder + number_entries(nb), diagonal, work, nLeft - nLeft2,
                iBlock + nb, jBlock, numberBlocks, level+1);
        });

        tg.wait();
}else ..

上の記述を適用した後のCPU Usageです。Windows10のタスクモニタ(左下)では、60%程度安定で100%になる様子は観測できませんが、 Visual Studioのモニタ(右下)の方が粒度が細かいようで、頻繁に変化する様子が観測され100%になることもあることが分かります。つまり、空いているCPUリソースを遊びがないようにタスクとして働かせていると考られます。



2021年10月20日水曜日

コレスキー分解の並列化

Eigen: Catalogue of dense decompositions

また、ベンチマークは、

https://eigen.tuxfamily.org/dox-devel/group__DenseDecompositionBenchmark.html

です。

現時点で、Eigenのコレスキー分解LDLTは、マルチスレッドに対応していません。またランク落ち(Rank-Revealing)についても、対応がありません。やはり自分で書くしかなさそうです。

コレスキー分解については、

1)コレスキー分解 - Wikipedia

2)コレスキー分解 - PukiWiki for PBCG Lab (tsukuba.ac.jp)

3)修正Cholesky分解 (coocan.jp)

4)main.dvi (kyushu-u.ac.jp)

を参照しました。

また、NSPの制約マトリクスは、基本的に疎行列です。正定対称化により、疎行列を維持できなくなる可能性が強いのですが、出来れば、そのままコレスキー分解をすることが望ましいです。従い、設計仕様は、

■SIMD化(Target AVX2)

■Multi thread化

■疎行列対応

■ランク落ち対応

を実現する修正コレスキー分解となります。

上の参照資料4)を見て思ったのですが、リニアソルバの核は、連立一次方程式の解法そのものです。取り掛かる前に見ておくべきでした。




2021年10月19日火曜日

コレスキー分解のSIMD化

 ボトルネックになっている部分をVisual Studio パフォーマンスプロファイラで調べ特定します。その部分のアセンブラコードを調べ、コンパイラが賢くAVX2命令を生成しているかどうかをチェックしました。意図したようにコンパイラが生成してくれればよいのですが、そうでない場合は、IntrinsicsCodeを書くことになります。(コンパイラは、ymmレジスタを使っていませんでした。)以下は、コンパイラの生成が思わしくないので自分で書いた唯一のコードです。

               for (i = 0; i < BLOCK; i += 4) {

               	   __m256d ymm_t0 = _mm256_load_pd(&aa[i + 0 + 0 * BLOCK]);
                   __m256d ymm_t1 = _mm256_load_pd(&aa[i + 0 + 1 * BLOCK]);
                   __m256d ymm_t2 = _mm256_load_pd(&aa[i + 0 + 2 * BLOCK]);
                   __m256d ymm_t3 = _mm256_load_pd(&aa[i + 0 + 3 * BLOCK]);
 
           

                    __m256d ymm_aUnderNow = _mm256_load_pd(aUnderNow);
                    
                    __m256d ymm_aboveNow0 = _mm256_broadcast_sd(  aboveNow);
                    __m256d ymm_aboveNow1 = _mm256_broadcast_sd(aboveNow + 1);
                    __m256d ymm_aboveNow2 = _mm256_broadcast_sd(aboveNow + 2);
                    __m256d ymm_aboveNow3 = _mm256_broadcast_sd(aboveNow + 3);

                    
                    __m256d ymm_tmp = _mm256_broadcast_sd(&(work[k]));
                    __m256d a=_mm256_mul_pd(ymm_aUnderNow, ymm_tmp);

                    ymm_t0 = _mm256_fnmadd_pd(a, ymm_aboveNow0,ymm_t0);//	c=-a*b+c; c-=ab;
                    ymm_t1 = _mm256_fnmadd_pd(a, ymm_aboveNow1, ymm_t1);//	c=-a*b+c; c-=ab;
                    ymm_t2 = _mm256_fnmadd_pd(a, ymm_aboveNow2, ymm_t2);//	c=-a*b+c; c-=ab;
                    ymm_t3 = _mm256_fnmadd_pd(a, ymm_aboveNow3, ymm_t3);//	c=-a*b+c; c-=ab;

                    aUnderNow += BLOCK;
                    aboveNow += BLOCK;
                }

                _mm256_store_pd(&aa[i + 0 + 0 * BLOCK], ymm_t0);

BLOCKは定数であり、コンパイラは、勝手にループ展開していましたが全く必要のないことです。(仮想的マシンの話ですが、同時命令実行数が4のとき、上記コードは、ベクタライズなしの同時命令実行数1に対して4x4=16倍の速度になります。) それはともかく、こんな風に書き換えるだけで10%程度速度が向上しました。

2021年10月18日月曜日

Ryzen Windows11

 大規模インスタンスで予定より全然性能が出ないので、おかしいなと思っていたら、

rigayaの日記兼メモ帳 【悲報】 AMD RyzenのL3 Cacheの大半を見失った件 【Windows11】 (fc2.com)

そういうことだったのですね。とりあえず、世界記録挑戦は、とりあえずパッチを待ってからにします。

2021年10月17日日曜日

Instance13速度改善Ⅲ

 開発環境をAVX2対応とした関係でマシンが変わっています。そこで、改善前のデータを取り直しました。

マシン:E3-1226V3 3.3GHz

改善前:1840sec

改善後:1409sec

今回の改善ポイントは、

■LinearSolverへのLoadProblem形式変更(constraints数の縮小)

PrimalSimplexの使用です。

■C++コンパイルオプションSSE2 Fastの使用

スケジュールナースは、既に64ビット版しかサポートしていないのでSSE2の使用をdefaultとしても問題ありません。

実行中の様子です。ほぼ100%CPUを使っています。


8月12日の初期報告では、7368secでしたので、マシンが変わったこともありますが、かなり改善出来たと思います。




2021年10月16日土曜日

BarrierSolverの製作

FirstOrderは、確かに速いのですが、収束の安定性に問題があることが分かりました。そこで、FirstOrderはDropし再びBarrierSolverを検討しようと思います。

 <Highsの高速化>

HighsのIPXは、マルチスレッド化がされていません。Julian Hall教授にどうしてなのか聞いてみました。現在の実装では、部分的にマルチスレッド化は、容易に出来る部分はあるけれども、アムダールの法則により効果が出にくい、とのことでした。特に変数が制約数を大幅に上回ることがなければ、難しいとのことです。(開発に資金提供する気があるか?とも聞かれました。。)

IPXの論文

これによれば、the crash basis (black bar at the bottom), for updating the basis during the interior point solve (dark grey), for running the linear solver (light grey) and for crossoverBasisに分類されるようです。なので恐らく、この全ての部分で並列化できないと効果が薄いという意味なのでしょう。

 C++ソースを見ると、Eigenも使っていませんし、AVX2によるベクタライズも見当たりません。

IPXに手を加えることは難しそうなので、別にBarrierSolverを作ることにしました。アルゴリズムは、同じMehrotra のプレディクタ・コレクタ法です。

untitled (orsj.or.jp)

IP4-PD-IIP.pdf (titech.ac.jp)

ボトルネックは、コレスキー分解で、

■解が近くなったときのランク落ち

■ベクタライズ(SIMD化)

■マルチスレッド化

がポイントになります。コレスキー分解については、Lapackの実装が最も速そうなのですが、ランク落ちについては、特にケアしていないので、やはり自分で書くしかなさそうです。(ちなみにLapackは、Fortranで書かれていますが、LLVMを使ってVisual Studioでもコンパイルすることが出来ました。良い時代になりました。)また、SuiteSparseについても見ましたが、対象とする超巨大NSPのインスタンス郡については、SuiteSparseを使用する意味が見いだせませんでした。まとめると、以下が設計仕様です。

■超巨大NSPインスタンス専用(それ以外は、Simplex使用)

■AVX2によるベクタライズ(それ以外はサポートしない)

■マルチスレッド

■クロスオーバは行わない

■Instance21で、10秒以下




2021年10月4日月曜日

物理コア数の検出

 前回結果より、スレッド数は、必ずしも多ければよい、という訳ではないということが分かったので、FirstOrderMethodを使用する場合には、次のAlgorithmで、スレッド数を決定しようと思います。

if 物理コア数が4以下 ⇒スレッド数=4

else ⇒スレッド数=物理コア数-1 

例えば、論理コア数が16、物理コア数が8のときは、7となります。

それで、物理コア数の検出が必要になるのですが、残念ながら物理コア数を得るLinux/Windowsで統一的に得る方法がない! らしいです。実際、以下のコードでコメントアウトしたのは、論理コア数を得るコードらしく、Ryzen5800Xでは、16が得られてしまいます。



//int cores = omp_get_num_procs();
//unsigned int cores = std::thread::hardware_concurrency();
      int cores = VCL_NAMESPACE::physicalProcessors();
      if (cores <= 4) {
            threads = cores;
        }else {
            threads = cores - 1;
        }
        cout << "threads=" << threads << endl;
}

これで、走らせてしまうと15スレッドが駆動されてしまうので、却ってパフォーマンスが低下してしまいます。(ベストな7スレッド駆動にしたい) そこで、StackOverFlowで紹介されたhttps://github.com/vectorclass/add-on/tree/master/physical_processors

を用いて、物理コア数を検出するようにしました。

2021年10月3日日曜日

The linear solver comparison

The graph below shows the processing time of each linear solver for three different instances, quite large(Instance21 from Scheduling Benchmarks), large(n120w8 from INRCⅡ), and nominal(Instance8 from Scheduling Benchmarks).



 When the instance size is quite large, the First Order Method has an overwhelming advantage. In contrast, the Simplex is dominant for light instances. In the middle graph, Interior Point Solver is dominant. As you can see, there is no single, all-purpose linear solver, so your best bet is to choose the best linear solver for your instance. In practical use, CLP will be the default linear solver since our nominal instance is less than 100 staffs per month. However, for long-term and large-scale applications such as 3 or 6 months, we should consider using Interior Point Solver or First Order Method.

2021年10月2日土曜日

Processing Times for First Order Method using SIMD and multiple threads

 Sparse Matrix Library "Eigen" is an excellent library for fast operations in general. However, if SpMV is the bottleneck of the processing time, there is room for improvement. Especially, AVX2 instruction sets make for faster computation time.

Here is a result.

The Eigen in the table is a bare SpMV description. The improved version of SIMD4 is a 128bits operation using SSE4.2 and LUT. 


SIMD8 is a 256bits operation using AVX2; instead of LUT, we use Load with mask instruction.


SIMD8 in the E3 processor is not so good as predicted. However, the latest CPU can expect the full power of AVX2. Thus we can expect almost three times performance gain for multi-threaded Eigen. 

In large scale operation, the first order solver has the advantage of light memory foot print for Simplex or interior point solver thus makes SIMD and multi- threading easily.