仮眠プログラマーのつぶやき

自分がプログラムやっていて、思いついたことをつぶやいていきます。

自作ゲームやツール、ソースなどを公開しております。
①ポンコツ自動車シュライシュラー
DOWNLOAD
②流体力学ソース付き
汚いほうDOWNLOAD
綺麗なほうDOWNLOAD
③ミニスタヲズ
DOWNLOAD
④地下鉄でGO
DOWNLOAD
⑤ババドン
DOWNLOAD
⑥圧縮拳(ツール)
DOWNLOAD
⑦複写拳
DOWNLOAD
⑧布シミュレーション
DOWNLOAD
⑨minecraft巨大電卓地形データ
DOWNLOAD
⑩フリュードランダー
デジゲー博頒布α版
DOWNLOAD
⑪パズドラルート解析GPGPU版
DOWNLOAD
⑫ゲーム「流体de月面着陸」
DOWNLOAD

HSP

HSPCL ver1.10リリース

http://sourceforge.jp/downloads/users/2/2285/HSPCL1.1.rar/

とりあえず加算サンプル。
まずデバイスを取得するところまでが鬼門だった気がする。 そもそも対応デバイスがそれなりに古くないグラボ(オンボード不可)だから相当限られるだろう

ローカルメモリ(共有メモリ)使えて、double計算できるだけでも十分だけど

HSPSHADプラグイン開発の続き

HSPSHADはHLSLを併用することによってGPGPUを擬似的にやるプラグインであった

だがこの前の記事のように、HLSLの記述に慣れなくてはいけなかったり、GPGPUには必要不可欠な指定ドット(任意の番地)への書き込みが不可能だったり、結局不便なことが多く、プラグイン開発計画は途中で凍結した。

そして、ここにきてopenCLを使ったGPGPUをはじめてみたところ、とても便利に感じてきたのでHSPプラグインにしてしまおうと考えた。


新しいプラグインの名前はHSPCL
openCLではGPUだけでなくCellやCPU上でも実行できる非常にオールマイティーな言語なので、使い方によってはCPUのマルチスレッド動作のためにこのプラグインを使うこともできる。

ただひとつ問題なのは、動作する環境がそれなりに新しいCPU、GPUでないといけないということだ

CPUではSandy Bridgeでは動作しなくIvy やSandy Bridge EからopenCLに対応、GPUではRadeonかnVidiaのGPUが搭載されているか、オンボードならintel でも最低でもIntel HD Graphics 4000つまりIvy Bridge世代からでしか使うことができない。つまりSandy Bridge 世代以前のノーパや、デスクトップでもオンボードのものでSandy以前ならこのプラグインは使えないということになる


まぁこの問題は時間が解決してくれるということで・・・

ところでHSPは32bitアプリケーションなのでメモリ確保は2Gまでしかできない
これはどうしようもない制限だが、このプラグインで確保できるグラフィックボード上のメモリは2Gではない。
グラボの最大メモリ容量まで可能なはずだ。
よって6Gbyteなど大容量メモリを搭載しているグラボがあればかなり大規模な計算も可能になってくる。

とりあえずそんなわけでHSPCLの開発を頑張ってみることにする


追記
ver1.0UP
とりあえずお使いの環境でopenCLが実行できるかの確認だけできる
http://sourceforge.jp/downloads/users/2/2280/HSPCL.zip/
ver 1.1
http://sourceforge.jp/downloads/users/2/2285/HSPCL1.1.rar/

話は変わってHGIMG4の話題
これはopenGLを使ったシェーダ機能も盛りこんであり、HSPdish上でも使えるかもしれないというもの!
ということはスマホ端末で擬似GPGPUが簡単にできるかもしれないというわけだ!
スマホを使った分散コンピューティングなんてできたら、非常に面白いかもしれない。
これは新しい情報待ちということで 

HSPでFFT(高速フーリエ変換)モジュール作った、あと多倍長乗算も

--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

#const beki 3;2^bekiこの配列が作られます、例beki=23の時 約400万配列×400万配列が800万配列に格納されます
#const kisuu 10000;1~32768まで、浮動小数配列ひとつに格納できる上限数=基数、何進数かという話、2^beki*log10底のkisuu が乗算後の値の桁数に値する

#module
#const beki beki@
#deffunc fft初期設定
n=1
repeat beki;桁数計算
n*=2
loop

ddim qcosq,n/2+1;sin cos予め計算格納
ddim qsinq,n/2+1
n2=n/2
repeat n/2+1
itjrad=3.1415926535897932384626*cnt/n2
qcosq.cnt=cos(itjrad)
qsinq.cnt=sin(itjrad)
loop
qcosq.n2=-1.0
qsinq.n2=0.0

itjn=1;バタフライインデックス作成
dim inind1,itjn
dim inind2,itjn*2
repeat beki
repeat itjn:inind2.cnt*=2:loop
dim inind1,itjn
memcpy inind1,inind2,4*itjn,0,0
dim inind2,itjn*2
memcpy inind2,inind1,4*itjn,0,0
repeat itjn:inind1.cnt+:loop
memcpy inind2,inind1,4*itjn,4*itjn,0
itjn*=2
loop
dim inind1,1;2の方に結果が
itjn=0

ddim chasr,n;キャッシュ
ddim chasi,n;キャッシュ
return

#deffunc fft array reitjo,array imitjo
q0q=1
q1q=n/2

repeat beki-1
qct0q=0
qct1q=q0q
qcn2q=0
repeat n/2
dup reitjoqccnt1,reitjo(inind2.qct0q)
dup imitjoqccnt1,imitjo(inind2.qct0q)
dup reitjoqccnt2,reitjo(inind2.qct1q)
dup imitjoqccnt2,imitjo(inind2.qct1q)
ita=reitjoqccnt2*qcosq.qcn2q-imitjoqccnt2*qsinq.qcn2q
itb=imitjoqccnt2*qcosq.qcn2q+reitjoqccnt2*qsinq.qcn2q
reitjoqccnt2=reitjoqccnt1-ita
imitjoqccnt2=imitjoqccnt1-itb
reitjoqccnt1+=ita
imitjoqccnt1+=itb
qct0q+
qct1q+
qcn2q+=q1q
if qct0q\q0q=0:qct0q+=q0q:qct1q+=q0q:qcn2q=0
loop
q0q*=2
q1q/=2
loop

qct0q=0
qct1q=q0q
qcn2q=0
repeat n/2
dup reitjoqccnt1,reitjo(inind2.cnt)
dup imitjoqccnt1,imitjo(inind2.cnt)
dup reitjoqccnt2,reitjo(inind2.qct1q)
dup imitjoqccnt2,imitjo(inind2.qct1q)
ita=reitjoqccnt2*qcosq.cnt-imitjoqccnt2*qsinq.cnt
itb=imitjoqccnt2*qcosq.cnt+reitjoqccnt2*qsinq.cnt
chasr.cnt=reitjoqccnt1+ita
chasi.cnt=imitjoqccnt1+itb
chasr.qct1q=reitjoqccnt1-ita
chasi.qct1q=imitjoqccnt1-itb
qct1q+
loop

memcpy reitjo,chasr,8*n,0,0
memcpy imitjo,chasi,8*n,0,0
return



#deffunc ifft array reitjo,array imitjo
q0q=1
q1q=n/2

repeat beki-1
qct0q=0
qct1q=q0q
qcn2q=n2
repeat n/2
dup reitjoqccnt1,reitjo(inind2.qct0q)
dup imitjoqccnt1,imitjo(inind2.qct0q)
dup reitjoqccnt2,reitjo(inind2.qct1q)
dup imitjoqccnt2,imitjo(inind2.qct1q)
ita=reitjoqccnt2*qcosq.qcn2q-imitjoqccnt2*qsinq.qcn2q
itb=imitjoqccnt2*qcosq.qcn2q+reitjoqccnt2*qsinq.qcn2q
reitjoqccnt2=reitjoqccnt1+ita
imitjoqccnt2=imitjoqccnt1+itb
reitjoqccnt1-=ita
imitjoqccnt1-=itb
qct0q+
qct1q+
qcn2q-=q1q
if qct0q\q0q=0:qct0q+=q0q:qct1q+=q0q:qcn2q=n2
loop
q0q*=2
q1q/=2
loop

qct0q=0
qct1q=q0q
qcn2q=n2
repeat n/2
dup reitjoqccnt1,reitjo(inind2.cnt)
dup imitjoqccnt1,imitjo(inind2.cnt)
dup reitjoqccnt2,reitjo(inind2.qct1q)
dup imitjoqccnt2,imitjo(inind2.qct1q)
ita=reitjoqccnt2*qcosq.qcn2q-imitjoqccnt2*qsinq.qcn2q
itb=imitjoqccnt2*qcosq.qcn2q+reitjoqccnt2*qsinq.qcn2q
chasr.cnt=(reitjoqccnt1-ita)/n
chasi.cnt=(imitjoqccnt1-itb)/n
chasr.qct1q=(reitjoqccnt1+ita)/n
chasi.qct1q=(imitjoqccnt1+itb)/n
qct1q+
qcn2q-
loop


memcpy reitjo,chasr,8*n,0,0
memcpy imitjo,chasi,8*n,0,0
return

#global





screen 0,480,700
fft初期設定
randomize
n=1
repeat beki
n*=2
loop


ddim rex,n
ddim rey,n
ddim reitjo,n;実数代入先
ddim imitjo,n;虚数代入先
ddim chasr,n;キャッシュ
ddim chasi,n;キャッシュ

; gosub*indexsakusei
repeat n/2;「A*B=c」のAとBにランダムな値を格納
rex.cnt=1.0*rnd(kisuu);2.7*sin(3.1415926535*490*cnt/n)+2.0*sin(3.1415926535*190.0*cnt/n)
loop
repeat n/2
rey.cnt=1.0*rnd(kisuu);1.0*sin(3.1415926535*50.0*cnt/n)+20.0*sin(3.1415926535*90.0*cnt/n)
loop

; aatim=gettime(4)*3600000+gettime(5)*60000+gettime(6)*1000+gettime(7)
; mes aatim
memcpy reitjo,rex,8*n,0,0
memcpy imitjo,rey,8*n,0,0

fft reitjo,imitjo

memcpy chasr,reitjo,8*n,0,0
memcpy chasi,imitjo,8*n,0,0


ReTm=0.0
imTm=0.0
ReTN=0.0
imTN=0.0
ReTm=chasr
ImTm=chasi
ReTN=chasr
ImTN=chasi
reitjo=(Retm*ImTm+ReTN*ImTN)/2.0
imitjo=( ((ReTN+ImTN)*(ReTN-ImTN)) + ( (ImTm+ReTm)*(ImTm-ReTm) ) )/4.0
itjn=n-1
repeat n-1,1
ReTm=chasr.cnt
ImTm=chasi.cnt
ReTN=chasr.itjn
ImTN=chasi.itjn
reitjo.cnt=0.5*(Retm*ImTm+ReTN*ImTN)
imitjo.cnt=0.25*( (ReTN+ImTN)*(ReTN-ImTN) + (ImTm+ReTm)*(ImTm-ReTm) )
itjn-
loop


ddim chasr,1;キャッシュ
ddim chasi,1;キャッシュ

ifft reitjo,imitjo


; mes gettime(4)*3600000+gettime(5)*60000+gettime(6)*1000+gettime(7)
; mes gettime(4)*3600000+gettime(5)*60000+gettime(6)*1000+gettime(7)-aatim

/*
repeat n
if absf(reitjo.cnt-int(reitjo.cnt))<0.7:if absf(reitjo.cnt-int(reitjo.cnt))>0.3:mes reitjo.cnt
; if absf(imitjo.cnt-int(imitjo.cnt))<0.7:if absf(imitjo.cnt-int(imitjo.cnt))>0.3:mes imitjo.cnt
loop
*/
;: bsave "reitjo",reitjo
; bsave "imitjo",imitjo
; stop
dim kekka,n;計算結果入るint型変数
gosub*繰り上げ計算
goto*hyouji2



*繰り上げ計算
w=n-1
w1=n-2
repeat n-1
kekka.w+=int(reitjo.w+0.5)
if kekka.w>=kisuu:kekka.w1+=kekka.w/kisuu:kekka.w\kisuu
w-
w1-
loop
kekka.w+=int(reitjo.w+0.5)
return

*hyouji;波形表示
color 0,0,0
repeat n
a=reitjo.cnt
if cnt=0:pos 0,int(240.0-a)
line cnt*640/n,int(240.0-a)
loop

color 255,0,0
repeat n
a=imitjo.cnt
if cnt=0:pos 0,int(240.0-a)
line cnt*640/n,int(240.0-a)
loop

color 0,255,0
repeat n
a=(imitjo.cnt*imitjo.cnt+reitjo.cnt*reitjo.cnt)
if cnt=0:pos 0,int(240.0-a)
line cnt*640/n,int(240.0-a)
loop
stop


*hyouji2
font "MS ゴシック",10,1
pos 0,0
mes "掛け算結果配列表示"
repeat limit(n,1,300)
if strlen(str(kekka.cnt))=1:mes "000"+str(kekka.cnt)
if strlen(str(kekka.cnt))=2:mes "00"+str(kekka.cnt)
if strlen(str(kekka.cnt))=3:mes "0"+str(kekka.cnt)
if strlen(str(kekka.cnt))>=4:mes str(kekka.cnt)
loop

pos 150,0
mes "掛ける数A配列表示"
repeat limit(n,1,300)
if strlen(str(int(rex.cnt)))=1:mes "000"+str(int(rex.cnt))
if strlen(str(int(rex.cnt)))=2:mes "00"+str(int(rex.cnt))
if strlen(str(int(rex.cnt)))=3:mes "0"+str(int(rex.cnt))
if strlen(str(int(rex.cnt)))>=4:mes str(int(rex.cnt))
loop

pos 300,0
mes "掛ける数B配列表示"
repeat limit(n,1,300)
if strlen(str(int(rey.cnt)))=1:mes "000"+str(int(rey.cnt))
if strlen(str(int(rey.cnt)))=2:mes "00"+str(int(rey.cnt))
if strlen(str(int(rey.cnt)))=3:mes "0"+str(int(rey.cnt))
if strlen(str(int(rey.cnt)))>=4:mes str(int(rey.cnt))
loop
stop

--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
















これはhsp 3.31で動きます。他のverでの動作は保証できません

ソースの補足として、
命令として「fft」「ifft」「fft初期設定」の3つを登録した。それぞれ離散フーリエ変換、逆フーリエ、バタフライ演算のインデックス作成、となっている。
「fft」 「ifft」ともに、とる引数は実数配列(double),虚数配列(double)で、変換結果がそれに返される。
配列の要素数は、2^beki個であることが条件、それ以上でも以下でも、1でもずれたらうまく計算されない。

実行すると、
実行画面の左には、A*B=CのCの配列が表示されている
実行画面の真ん中には、A*B=CのAの配列が表示されている
実行画面の右には、A*B=CのBの配列が表示されている
bekiの値を増やすほど、多くの桁数の乗算が行える
仮にbekiを23にすれば、約400万配列×400万配列が800万配列に格納される
bekiが23でkisuuが10000の場合、約1600万桁×1600万桁=3200万桁の乗算となる


fftモジュールはもちろん乗算以外の用途にも使える。 

HSPで高速並列演算

(※並列にできるのは整数型、float型の四則演算のみです。命令や関数の並列処理は難or不可です。)


 前にブログに書いたようにhgimg3やHSPDXfixなどのプラグインは標準命令と比べて描画がとても早い。
 標準命令ではスペックに限界があり、例えば自分の作りたいシューティングゲームなどがとても60FPS出せない場合はこれらのプラグインを使うと解決されることが多い。

 たいていこういうプラグインは、ハードウェアに処理をさせているため高速に描画できます、とか書いてある。

o0380008511359222298


 その「ハードウェア」ってのはつまるところGPU(広義のグラフィックボードやビデオカード)のことで、端的に言えばこいつもCPUみたく四則演算ができる演算装置だ。ただCPUみたく高度な処理はきなくて、全画面のピクセルの書き換えとか、同時に何万もの同じ処理をしたいときなどに向いているとされている。

 GPUは本来グラフィック関連の処理に特化されたチップなので、その用途以外で使われることは考えられていなかったようなのだが、GPGPUと言って最近このGPUがCPUにかわって様々な処理をするようになってきた。 動画変換や計算科学(流体力学やブラックホールシミュレーションなど)といった、グラフィックに関係ないところに使われるようになってきて、それもまたCPUで処理するより何十、何百倍も早いそうなのだ。

 http://dic.nicovideo.jp/a/gpgpu (ニコニコ大百科がまとまってて分かりやすい!)
 後にも書くが、流体力学の計算をCPUとGPUでやらせたときの処理時間比が私のパソコンでは最大180倍になるという結果になった!
これは2^n×2^nの計算領域でCPUとGPUで流体計算させたときの1ループの処理時間(秒)

o0703046311359223918



 あまりにもGPUの計算時間が早くて、いまいちすごさがパッとしないので処理時間の比をとってみた


o0725041611359223917


 横軸が若干省略されてしまっているが、
512^2と2048^2の間にある1024^2の計算格子領域で一番GPU/CPU比が高く181倍も高速化ができた。
大岡山に存在するとされている某スパコンは、主にこのGPUで計算させているのでCPUのみのスパコンと比べすごく演算が早いとのこと。(もしCPUだけで同じ性能にしようとすると莫大なお金がかかる)
周波数が高いわけではなく、何万もの計算を並列に行なえるから結果的に早いのだそうだ。


・・・さて長い前置きは置いておいて、 私はHSPコンテストで分子シミュレーター
o0640048011359227067


とか
重力シミュレーター
o0431032111359227068

とか


まぁ普通にくだらないクソソフトを毎年量産しては応募しているのだが・・・
何がクソかというと、なんといっても処理が重いのだ!

これはプログラムが非効率なのもあるが、多数の粒子や相互作用などを計算すること自体がそもそも物理的に膨大な計算量を必要とする。

非力な私のノーパソのCPUでは、一億桁×一億桁の計算に何日もかかるというのに、分子のシミュレーションをリアルタイムでやろうとするなんぞ分不相応というか背伸びのしすぎであった・・・

今よりもCPUの性能が上がればいいのだが、生憎CPUのクロック数はもう頭打ちと言われていてせいぜい3.6Ghzが一般CPUの限界だから、今の1.5倍にしかならない。
HSPはマルチスレッド対応ではないため、私のやりたいシミュレーション系のプログラムはもうCPUだけでは限界なのだ。 多重起動の擬似マルチスレッド化とかして工夫はしてきたものやはりCPUに頼っているだけでは根本解決にはならない
そこで、本格的にGPGPUプログラミングに手を出そうと考えていた。

調べていくうちにGPGPUプログラミングのためにはnVIDIA製のGPU搭載デスクトップPCを買い、CUDAを勉強するのが一番手っ取り早いということがわかったのだが・・・

1 買う金がない
2 デスクトップ置く場所がない
3 CUDAを勉強する時間がない

の理由から諦めかけていた。
HSPでGPGPUプログラムができれば全て解決んだけど・・と思って何気なくHSPのサンプルを眺めていたら・・ これはもしかしたら、 HSPでGPGPUぽいことができるんじゃないのかと、そんな確信に近いある予感がひらめいた!!

きっかけはeasy3Dのサンプルのe3dhsp3_CbCr.hspを見た時だった
実行してみると、左には普通の3Dモデルのサンプルがリアルタイムで動いていて、右にはそのモデルの色だけがセピア色に変換されて同じく動いている

o0230017311356902500


 そうセピア変換とは、昔私がブログに無駄に長いソースを貼っつけたアレである。

まず画面の全ドットのrgb成分をyuvに逐一変換してuvに一定値を代入してまたrgbに変換・・・というやつである。
これをeasy3Dの中で実現しているということは、GPUにrgbをyuvに変換する並列処理をさせているということに他ならない!(CPUでリアルタイムなんかまず無理であろうし)

ということはソースかdllのプログラムのどこかにこのrgb成分をyuv成分に変換する、掛け、割り、足し、引きの計算式が書かれているはずだ!
そこのソースを改変すれば、自分の好きな手順、順番で四則演算ができるということではないのか!? この仮説が正しければeasy3Dのプラグインを使って高速並列計算が可能ということになる!
高速シミュレーションの道を開くために、執念でdllを解析しようかとか考えていたら、思いのほかもっと近くにその四則演算部分のプログラムを発見した
e3dhsp3.dllのと同じフォルダにある「E3D_HLSL」だ!
o0480024011361295228


そのフォルダにposteffect.fxが入っていて、こいつがセピア変換の四則演算を記述してあるソースだった。
調べたらこの.fxという第二のソースはHLSLといってHigh Level Shader Languageというもうひとつの言語らしい
おや、fxファイルをメモ帳で開くとプログラムソースが見れる!
ということは編集もできるではないか!?

o0730063411361295227



このファイルに自分で式を打って正常に動けば、自分の思い通りの並列処理が可能ってこと!?
試しに計算式を消したり追加したり実験してみたところ、コンパイルが正常に行なわれた

驚くべきことに、HSPで並列処理ができることがこれでわかった!!(正確にはHLSLでだけど・・・)


さらに驚くべきことに、easy3Dのソースを調べていくと律儀にも、私みたいなユーザーがposteffect.fxを改変できるようにと、なんとすでにカスタム領域が備わっていたのだ!
あとからおちゃっこさんのブログを見てわかったのだが「カスタムHLSL機能の需要は?」という名目で2年以上も前にこの自由に定義できるカスタムHLSLの機能を作っていたらしい。
つまりHSP側で並列計算させたいデータをeasy3Dの命令でGPUに送り、GPU側の処理の操作はHLSL(posteffect.fx)のカスタム領域で自分の好きな処理をさせられ、処理結果をまたHSP側で取り出すといった一連の操作が、我々ユーザーのためにすでに可能にされていたということだ。

HLSLのカスタムを可能にしてくれたおちゃっこさんには本当に感謝しないといけない!!
おちゃっこさんは需要なさそうと気にしていらっしゃったようだけど、これぞまさに私の求めていたもの!
この機能に出会えて本当に良かった!




HSPで並列演算ができるとわかった今、シミュレーションプログラムの幅が大きく広がること間違いなしである!!
これでもう数値解析やシミュレーションのデバッグに何日もかけないで済む!
HSPコンテストに糞ソフトと呼ばれるような生半可なシミュレーションソフトは出さないぞ。(呼んでいるのは主に自分)
さらにさらに驚くべきことに自分の非力なノーパにもそれなりのGPU・・オンボードだが・・がそなわっていたらしい!
o0399042211361313932

今の今まで気づかなかった・・そもそもGPUなんて付いてないものかと・・
最大単精度浮動小数点バッファまで対応してるみたいだし、シェーダモデル4.0で、なんだか考えられる中で最高の環境が気づかぬ間に手元にそろっていたようだ。
で、 具体的な並列処理の方法だが、自分の手順をさかのぼってみて、例をおりまぜてまとめてみた。





まずHSP3.2は入っている前提で

・easy3dのサイトからEasy3DForHSP3_ver5233をダウンロードする
・解凍してreadme読む ・readmeの言われたとおりにする
・次に解凍した中のE3D_HLSLのフォルダもhsp3フォルダにコピー
・下のリンクからダウンロードしてきたposteffect.fxを、hsp3フォルダのE3D_HLSLの中にあるposteffect.fxに上書き

いきなりHLSLを自分でカスタムするのは難しいので、他人のカスタムしたものを参考に少しづつ変えて勉強していくのが一番いいだろう。またHLSLのサイトにいけば命令とか演算子が全部見られる。
サンプル 上のサンプルは、左の2枚の画像を加算、減算、割り算、掛け算しその結果を右に出力するというものだ。1~8のキー操作を受け付けて8種類の処理が見れるようになっている。

o0777040011361295229

話を戻して、例えば足し算の並列演算をしたい場合なら、

・#include "e3dhsp3.as"して、E3DCreateRenderTargetTextureで、足す方を格納する2枚のテクスチャと、足した結果を出力するためのテクスチャをひとつ、合計3つ、同じ大きさで作る。
 ようはc=a+bという計算をさせたい場合、aというテクスチャとbというテクスチャとcというテクスチャを作ってください、ということだ。
 テクスチャの大きさは、並列に処理したい要素数によって変わってくる。65536要素の計算を並列に処理したいときは256*256の大きさのテクスチャをつくればいい。
・E3DShaderConstUserTexで足す2枚のテクスチャをGPUに送る
・E3DCallUserShaderでshadernoを0、passnoを1にして、足し算実行。これでE3DCallUserShaderの描画先のスワップチェインID(テクスチャ)に足し算結果が出力される。(ここらへんはサンプルの中にソースあり)





といった流れが、足し算の並列処理の方法だ。
passnoを変えれば割り算などを指定できる。
当然posteffect.fxを書き換えれば3つ4つ連続の足し算や、割った余りを求めるものなど、無限の処理が可能となる。
出力された計算結果(テクスチャ)はE3DShaderConstUserTexで指定しE3DCallUserShaderでつぎの処理に使いまわしていく事ができる。

上のスクリーンショットにも書いたが、E3DCreateRenderTargetTextureで作成するテクスチャのbit数を16bitや32bitにすれば、浮動小数点同士の足し算掛け算を並列処理できることとなる。

問題は、並列処理するときのデータが全てテクスチャ形式で行なわれているという点である。

はっきりいってHSP側で作ったDOUBLE型の配列変数を浮動小数点テクスチャに変換するというのはかなり至難の業である。(整数テクスチャならば話は簡単で、bmp形式でデータのやりとりを行なえばよい。)
というのもE3DCreateRenderTargetTextureで作ったテクスチャに1ドットずつ特定の値を書き込む命令がない。(というかeasy3Dはおそらくこの使いかたは想定していなかった)。だから何らかの方法で工夫しないといけないのだ。
それに、作れるテクスチャの精度が最大32bitなので標準命令のdoble型(64bit)と比べ精度が落ちる。グラフィックチップによっては16bitが限界だったりするから、CPUでのdoble型の計算をGPUで精度を落とさずやるには、16bitを4回計算みたいな感じでやるなどさらなる工夫が必要となる。

そして浮動小数点テクスチャから配列変数に戻すのも同様に大変である。
そこらへんを考慮した設計でないと、予想と結果が大きく異なってしまう可能性がある。


ということは、HSP側で配列変数を作る必要は本当にあるのかという疑問がわく。
並列演算出力結果を例えば視覚的に見るだけで良い場合(流体解析や分子動力学など)なら、posteffect.fx内で可視化プログラムを作ってしまえば、easy3Dのテクスチャレンダリング命令でスクリーンにいっきに可視化できる。わざわざHSP側で作った配列変数に値を戻さなくても、だ。
理想はそのような形で、最初から最後までテクスチャとして数値を保存しておく方が楽で良いのだが、厳密な数値を文字列として出力しなければいけないなど、どうしてもHSP側の変数に数値を戻す作業が欠かせない場合がでてくるだろう。 浮動小数点バッファを小数のまま保存できたりHSPのメモリに転送できる命令があればいいが、ない以上、1枚の浮動小数点バッファを何枚ものbmpに情報を分解して保存するアルゴリズムを構築するしかないのかもしれない。少なくとも私のお思いつく方法は今のとこそれしかない・・・
それだけがHSPでの並列処理においての大きな問題である。逆にもし単精度浮動小数点テクスチャと標準命令のdouble変数が両方向に対して変換可能になれば、まぁ64→32bitで精度が悪くなるのは仕方ないとして、タイトル通りの「HSPで高速並列演算」が可能になるというわけだ。これからは浮動小数点の構造の勉強をしてなんとか数値の取り込み、取り出しの問題を解決していきたいと思う・・・。








最後に、標準命令だけで作った流体力学プログラムを、HLSLに移植して並列処理化してみた。 計算は32bit浮動小数点バッファを使用、格子数は256*256で、CPUの80倍以上高速に処理できた。

o0512025611361324083

(一応解説すると、上から流体が流れてきてて、画面上の方にあるみにくいけど黒い動かないのが障害物。画面右は圧力を可視化した物。緑=高圧、赤=低圧)

これは、出力がスクリーンのみなので、posteffect.fxのほうに可視化プログラムが記述されており、可視化したテクスチャをeasy3DのE3DRenderSpriteでスクリーンにコピーしている、といったところ。
また入力のほうは、障害物の場所情報とか、0(なし) 1(あり) の2進数で表せるような情報しかないため、HSP側でbmpの生成→easy3Dでbmpの読み込み、bmpをスプライト(整数)に貼りつけ→32bit浮動小数点テクスチャにスプライト貼りつけ、・・と、数値の劣化なくGPUに送ることが出来ている。
並列処理を利用した数値解析シミュレーション特に流体力学はHSPコンテスト2011に向けて温めているもののひとつなのでネタバレはこのくらいにしておくとしよう。


次回:流体力学のプログラムを作りたい!その1













追記
HSP→浮動小数点バッファ への32bit浮動小数のデータ転送は
userFL4_0 ~ userFL4_9 のユーザー定義定数を使って正確に 行なえることが分かった。

また32bit浮動小数点を、4つの8bit*3(r,g,b)バッファに分解して転送して、HLSL側で数値を復元する方法も試みたが、
仮数部の最小1bit部分でたまに誤差が生まれる模様・・・


また浮動小数データの分割&整数値化による GPUバッファ→HSP へのデータ転送は全くうまくいかず・・・orz

HGIMG3で軸を中心としたカメラの回転

フライトシミュレータ・・・・具体的にいうとエースコンバットや零式艦上戦闘機などのゲームをHSPで真似て作ってみたいと思うのは多分私だけではないはず・・



そこでHGIMG3を使ってやろうとするのだが、どうもsetangで指定するカメラの回転角度の計算でいつもつまずく。



早い話、視線を軸としてカメラを回転させたいのだ。

キー操作で説明すると、→←キーでローリングし、↓キーでピッチが上がり↑でピッチが下がるようなプログラムにしたいという訳である。


余談・・・

(実はHGIMG3のサンプルの中に飛行機ゲームがあるのだが、残念ながらそれは私の考えていたキーと回転の対応をしていなかった。逆にd3mのモジュールの中にこれだ、というサンプルを見つけた。しかしd3mとHGIMG3はまったく命令が違う上に、d3mの視点の回転命令が1つの命令にまとめられてしまって肝心な計算式がブラックボックスの中になっていた・・)



結果からいうと、setangの3つの角度だけではこの問題は解決できない!!


ではどうするかというと、新たに視点方向用の3つの方向ベクトルを作り、それで主な視点の計算を行うという方法をとらなければならない。


3つのベクトルというのは自分基準のx、y、z軸ベクトル、例えば飛行機でいうところの

 ・パイロットが見るような視線ベクトルと

 ・その法線である翼に平行なベクトル

 ・飛行機を真下から真上に見たときのベクトル

であり、この3つはそれぞれ直角同士でなければならないという制約がある。

補足・・・この3つのベクトルが視線方向を決めるのだから、逆に視線が回転すれば3つのベクトル全ての値が変わるというのは言うまでもない



主役を降ろされたsetangの3つの角度の使いどころは、本当に最後画面に出力するときにsetangでカメラ向きをセットするときだけ、なのである。


次に3つのベクトルでピッチやヨーやロールを取り扱う計算には、絶対クォータニオン(四元数)を知らなくては話が始まらない。

クオータニオンとは任意の軸に対してベクトルを回転させることができる計算方法である。

(計算には行列などがでてくる)



まぁクオータニオンの詳しい説明はこのブログでは端折らせてもらうとして、まとめると





ピッチ、ヨー、ロールのアルゴリズムの大まかな流れ・・・以下の例はロール

まず上記の3つのベクトルを作る

視線ベクトルを軸にして翼に平行なベクトルと機体に垂直なベクトルを回転する(ここでクオータニオンを使う)。

3つのベクトルをsetangの3つの値に変換(この変換が面倒)

setangでカメラ設定


という感じである。



ピッチのアルゴリズムは、翼に平行なベクトルを軸にして視線ベクトルと機体に垂直なベクトルを回転してやればいい。同じように機体に垂直ベクトルを軸にして視線ベクトルと翼ベクトルを回転させればヨーのアルゴリズムができる。



クオータニオンや座標変換については説明するのがすごく大変なので「いなえの鉛筆」さんなどのもっとわかりやすいサイトを見てもらって、このブログでは主に完成系のプログラムを公開することにする。





↓は上の3つの方向ベクトルsetangの3つの値stx,sty,stzに変換するプログラム(サブルーチン)

まず前提として3つの視点方向用ベクトルを

px py pz   (パイロットが見る視線ベクトル)

jx jy jz    (翼に平行なベクトル)

tx ty tz    (縦のベクトル)

とする。(これらは単位ベクトルでかつそれぞれが垂直である必要がある)


※ベクトル用の変数以外に勝手に使ってる変数は

a,c,d,x,y,z,w,ra,hx,hy,hz,stz1,stz2,q.0~q.3,r.0~r.3,rp.0~rp.3,rpq.0~rpq.3,π

配列は皆実数型でπ=3.1415・・・とする



*main

<クオータニオンの計算など>

・・・・・

・・・・・

gosub*henkan

setang HGOBJ_CAMERA,stx,sty,stz

・・・・・

・・・・・

goto*main




*henkan
if tx=0.0{
a=0.0
}else{
if (ty*jx-tx*jy)!0.0:a=(px*ty-py*tx)/(ty*jx-tx*jy)*jx/tx-px/tx:else:a=0.0
}
c=sqrt(1.0/(1.0+a*a))
d=c*a
hx=c*px+d*tx
if jx>0.0{
if hx<0.0:c=-c:d=-d
}else{
if hx>0.0:c=-c:d=-d
}
stx=atan(-d,c)
hz=c*pz+d*tz
hx=c*px+d*tx
hy=c*py+d*ty
sty=0.0
if (jz*hx-jx*hz)!0.0:d=hx/(jz*hx-jx*hz)
if hx!0.0:c=-d*jx/hx
sty=-atan(-d,c)
x=tx:y=ty:z=tz
ra=-stx/2.0
q=cos(ra):q.1=jx*sin(ra):q.2=jy*sin(ra):q.3=jz*sin(ra):r=cos(ra):r.1=-jx*sin(ra):r.2=-jy*sin(ra):r.3=-jz*sin(ra)
gosub*ks
hx=x:hy=y:hz=z

x=jx:y=jy:z=jz
ra=-sty/2.0
q=cos(ra):q.1=hx*sin(ra):q.2=hy*sin(ra):q.3=hz*sin(ra):r=cos(ra):r.1=-hx*sin(ra):r.2=-hy*sin(ra):r.3=-hz*sin(ra)
gosub*ks
c=hy
d=sqrt(1.0-c*c)
stz1=atan(d,c)
c=y
d=sqrt(1.0-c*c)
stz2=atan(d,c)
if absf(stz2-stz1-π/2.0)<0.0005:stz=-stz1
if absf(stz2+stz1-π*1.5)<0.0005:stz=-stz1
if absf(stz1-stz2-π/2.0)<0.0005:stz=-2.0*π+stz1
if absf(stz2+stz1-π/2.0)<0.0005:stz=-2.0*π+stz1
return





*ks
ddim rp,4;r*p
rp=-r.1*x-r.2*y-r.3*z
rp.1=r*x+r.2*z-r.3*y
rp.2=r*y+r.3*x-r.1*z
rp.3=r*z+r.1*y-r.2*x
ddim rpq,4;r*p*q
rpq=rp*q-q.1*rp.1-q.2*rp.2-q.3*rp.3
rpq.1=rp*q.1+q*rp.1+rp.2*q.3-rp.3*q.2
rpq.2=rp*q.2+q*rp.2+rp.3*q.1-rp.1*q.3
rpq.3=rp*q.3+q*rp.3+rp.1*q.2-rp.2*q.1
x=rpq.1
y=rpq.2
z=rpq.3
w=sqrt(x*x+y*y+z*z)
x=x/w
y=y/w
z=z/w
return





この、ベクトル→setang回転角変換プログラムは私がオリジナルで考えて、なおかつ途中から適当にやったら偶然できたプログラムなので実は自分でも良くわかってない。

とても見にくくて、スパゲッティになっているが、あしからず・・




次にベクトルの回転のみのプログラム(=クオータニオン)

※このプログラムだけではsetangの3つの値に変化はない

→←キーでローリング、↑↓キーでピッチ上下するというもの。


説明・・・・・回転させたいベクトルの値を x、 y、 z に代入し、軸となるベクトルは配列変数 q と r に sin(rad) を掛けて代入する。 配列の最初には cos(rad) をそれぞれ入れておくこと。 この状態で ↑ の *ks へ gosub すると回転計算結果が x、 y、 z に出力される。

またradにはベクトル回転角度を初期設定しておく


*main

stick key,$3ff
hgdraw
hgsync 16

if key&4:rad=0.02:gosub*k1
if key&1:rad=-0.02:gosub*k1
if key&8:rad=0.02:gosub*k2
if key&2:rad=-0.02::gosub*k2

<中略・・>

・・・・

・・・・

・・・・

goto*main



*k1
repeat 2
if cnt=0:x=jx:y=jy:z=jz
if cnt=1:x=tx:y=ty:z=tz
q=cos(rad):q.1=px*sin(rad):q.2=py*sin(rad):q.3=pz*sin(rad):r=cos(rad):r.1=-px*sin(rad):r.2=-py*sin(rad):r.3=-pz*sin(rad)
gosub*ks
if cnt=0:jx=x:jy=y:jz=z
if cnt=1:tx=x:ty=y:tz=z
loop
return

*k2
repeat 2
if cnt=0:x=px:y=py:z=pz
if cnt=1:x=tx:y=ty:z=tz
q=cos(rad):q.1=jx*sin(rad):q.2=jy*sin(rad):q.3=jz*sin(rad):r=cos(rad):r.1=-jx*sin(rad):r.2=-jy*sin(rad):r.3=-jz*sin(rad)
gosub*ks
if cnt=0:px=x:py=y:pz=z
if cnt=1:tx=x:ty=y:tz=z
loop

return


*ks

<上と同じなので中略>

return









一応この2つのプログラムをうまく駆使すれば、フライトシミュレータっぽくHGIMG3でも任意の軸回転ができるハズ!


ちなみに・・

初期設定で設定すべき

px py pz   (パイロットが見るような視線ベクトル)

jx jy jz    (翼に平行なベクトル)

tx ty tz    (機体に垂直なベクトル)

は、適当に決めていいのだが・・・

かならずベクトルは単位ベクトルになるように、また3つそれぞれが直角になるように設定お願いします。

そうじゃないと私のプログラムでは動きません・・・


また発見したバグは

px=0.0
py=0.0
pz=1.0

jx=1.0
jy=0.0
jz=0.0

tx=0.0
ty=1.0
tz=0.0
で初期設定すると、最初ピッチを上げようとすると画面がバグる

というものですが、すぐにローリング(左右キーを押す)すればバグは消えます。


初期設定では

px=sin(1)
py=0.0
pz=cos(1)

jx=cos(1)
jy=0.0
jz=-sin(1)

tx=0.0
ty=1.0
tz=0.0
をお勧めします。

プロフィール

toropippi

記事検索
アクセスカウンター

    QRコード
    QRコード
    • ライブドアブログ