fortranのライブラリ

Debian busterで古いプログラムを動かそうと思ったら、libg2c.so.0が無いと言われた。調べてみたら、これはg77のライブラリらしい。しかし、今のdebianではg77はgfortranに置き換えられてしまっている。昔のubuntuのdebを使うと動かせるらしいので、試してみた。http://old-releases.ubuntu.com/ubuntu/pool/universe/g/gcc-3.4/から必要なdebファイルを取ってきて、以下のようにインストールする。

dpkg -i --force-depends cpp-3.4_3.4.6-6ubuntu5_i386.deb gcc-3.4_3.4.6-6ubuntu5_i386.deb gcc-3.4-base_3.4.6-6ubuntu5_i386.deb libg2c0_3.4.6-6ubuntu5_i386.deb libg2c0-dev_3.4.6-6ubuntu5_i386.deb g77-3.4_3.4.6-6ubuntu5_i386.deb
ln -s /usr/bin/g77-3.4 /usr/bin/g77

必要で無い操作もある気がするけど、とりあえず動いたので、それなりうまく行っている。しかし、upgradeするときには、文句を言われるので、以下のように、一旦uninstallしてからupgradeする必要がある。

dpkg -r cpp-3.4:i386 gcc-3.4:i386 gcc-3.4-base:i386 libg2c0:i386 libg2c0-dev:i386 g77-3.4:i386

Rとrubyの比較

LANを介してデータのやり取りをする必要があったので、rubyのsocketを使って、必要なプログラムを書いた。通常は一つのサーバーに複数のクライアントが何回も接続することが想定されるので、threadを使うことが多い。しかし、今回は一つのクライアントが一度だけ接続するというケースだったので、threadが必要が無いことに気がついたので、以下のように単純化した。

serv=TCPServer.open(2051)
cl=serv.accept
str=cl.readpartial(10)
cl.close
serv.close

acceptでclientが接続するまで待ってくれる。クライアント側は、TCPSocket.open("localhost",2051)で接続して、読み書きすれば良い。

ふと、これをRでできないかを考えてみた。クライアント側は、以下のように、socketを開けて、読み書きすれば良い。

tcp<-socketConnection("localhost",2051)
write("123",tcp)
cat("123",file=tcp)
readLines(tcp,n=1)
close(tcp)

サーバーもやってみたら簡単だった。

tcp<-socketConnection(port=2051,server=TRUE)
close(tcp)

定義するときに、クライアントが接続するまで待ってくれるので、クライアント側とほとんど同じ書き方で、サーバー側のプログラムも書けてしまう。ただ、複数のクライアントが接続する場合は面倒そうだ。

Rでrubyと同様にsocketを用いた通信ができそうだが、違いがある。rubyではreadpartialでは、データが来るまで待ってくれるが、RだとreadBinとかでは、データがなかったら読まずに終わってしまう。特定の長さのデータを読み込みたいときには、rubyではread(length)とすれば良いが、Rでは長さを指定しても、それよりも短いデータしか読み込まない可能性があるので、読み込んだデータを数えて、特定の長さになるまで繰り返さないといけない。これは面倒だが、何か楽に固定長を読み込む方法は無いのかな。

平方根の計算

球面調和関数を取り扱っていると、Condon-Shortley係数とかGaunt係数というものを使うことがある。この値は表になっているのだが、maximaを使って行列の形に表すプログラムを書いてみた。

load(clebsch_gordan);
ck(k,l1,m1,l2,m2):= rootscontract((-1)^m1*sqrt(2*l1+1)*sqrt(2*l2+1)*wigner_3j(l1,l2,k,0,0,0)*wigner_3j(l1,l2,k,-m1,m2,m1-m2))$
mck(k,l):=genmatrix(lambda([a,b],ck(k,l,l+1-a,l,l+1-b)),2*l+1,2*l+1)$

苦労したのは、平方根を単純化する部分である。sqrt(2)*sqrt(6)を2*sqrt(3)と変形するやり方が分からなくて、自乗してからsqrtしたりする関数を書いたりしていたが、rootscontractを使うと良いことがやっとわかった。

first trial
AMDのRyzenはコストパフォーマンスがよいということで、興味を持っていたが、Linuxのインストールが難しいという噂を耳にしていた。最近、メインPCがいきなり電源が落ちることがあって、最初は電源ケーブルが抜けかかっているだけかと思っていたのだが、何度か再発するので、PCを更新することにした。そろそろRyzenにもLinuxのインストールが簡単にできるようになっているだろうということで、新しいPCはRyzen5にした。 今回は、CPUがryzen5になったことと、SSDとHDD両方載せていることが、今までのPCとは大きく異なる点である。パーティションを切るときに、SSDをルートに、HDDをhomeにして、インストールしようとしたが、EFIが無いと言われて進まない。よく分からなかったので、SSDのみに自動でパーティションを設定して、EFIを自動生成して、そこからSWAPを除いてルートに割り当てたあと、HDDにhomeを割り当てたら、進むことが出来た。 インストール自体は、うまく行ったようだが、再起動しても、画面が真っ黒なままである。今日は時間が無いので、後日に再挑戦しよう。

2020/6/1追記 とりあえず、前回はDebian10.3をUSBブートで試したが、最新は10.4だったので、それで試すことにした。しかし、同じような症状だったので、bootの設定がおかしいのかと思って、BIOSからいじってみた。Debianが二つあるように見えるけど、single user modeなのかと思っていたが、それだと画面が真っ黒なことの説明ができない。二つある方のどちらかを選ぶべきなのかと思って、いろいろ試していたら、grubが立ち上がるようになった。webの情報だと、ここからいろいろ設定しないといけないという話だったが、その後は問題無く動いているようである。まあ、しばらく使って試してみよう。

古いPCにlubuntu
コロナウイルスの影響で，人と会う機会を減らすためにオンラインで様々なことをすることが求められるようになっている．いろいろな出先でネット環境を手軽に使えるようにするために，古いノートPCを発掘したので，これにlinuxを入れてみた．入れたOSはlubuntu18.04である．もう少ししたら20.04が出るはずなのだが．さらに，古いPCなので，amd64では無くi386である．インストール自体はうまくいったのだが，wifiとSDカードが認識できかった．前者に対しては，bmc-kernel-sourceを入れたらうまく行った．後者に関しては，/etc/modprobe.d/optionsにoptions sdhci debug_quirks=1を記述したら認識するようになった．

win7とwinXPの違い

測定用のPCにおいては、WinXPやWin7などのサポートの切れたOSを使っているものが多くある。これらのPCからデータを取り出すときに、どのようにするかが問題である。直接USBメモリなどを挿すと、そこからウイルスが感染する可能性があるし、ネットワークに繋ぐと同様である。そこで、以前構築したのが、閉じたネットワークにLinuxPCを参加させて、そのPCからUSBメモリで取り出す仕組みである。Linuxにウイルスが感染してそれがWindowsに伝染ったらどうしようも無いが、その可能性は非常に低くなる。

十年ぐらい前に作った時には、/etc/fstabには以下のように記述していた。

//ip_win7/DATA  /home/user/dir1/  smbfs iocharset=utf8,password=,defaults,ro,user 0 0
//ip_winxp/DATA  /home/user/dir2/  smbfs iocharset=utf8,password=,defaults,ro,user 0 0

これで、起動時にwindowsのディレクトリがLinuxにマウントされる。今回、PCも古くなったので、一新することになった。Debian10を入れて、sambaを使うためにcifs-utilsを入れる。smbfsとcifsで使い方はほとんど一緒だと思ったのだが、微妙なところが違って苦労した。結局、以下のように記述した。

//ip_win7/DATA /home/user/dir1 cifs iocharset=utf8,password=,user,dir_mode=0755,file_mode=0755 0 0
//ip_winxp/Data /home/user/dir2 cifs iocharset=utf8,password=,user,dir_mode=0755,file_mode=0755,vers=1.0 0 0

まず、重要なのが、ディレクトリの最後のスラッシュの有無である。次に問題なのが、マウントしたときのmodeの指定である。最後に悩んだのが、win7とwinXPの違いである。XPではSMB1.0なのだが、vers=1.0を指定しないとアクセスできない。7はSMB2.0にも対応しているので、指定しないでも大丈夫なようだ。

BLAS/LAPACKとparallel

行列を使って対角化などの計算をやっていたのだが、高速化ができないかを試してみた。

Rの通常の計算はそれほど速くないが、Rの行列計算は、BLAS/LAPACKを使っているので、それなりに速い。sessionInfo()とすると、どのBLAS/LAPACKを使っているかが表示される。標準のもの以外にも、openblasやatlasやintel-mklなどがあり、これらを使うことにより、多少速度が変わるらしい。debianの場合は、libopenblas-base, libatlas3-base, intel-mklをaptで入れて、

sudo update-alternatives --config libblas.so.3-x86_64-linux-gnu
sudo update-alternatives --config liblapack.so.3-x86_64-linux-gnu

とすると、これらを変更できる。私の環境では、atlasが若干速かったので、これを使うようにした。

繰り返し計算などで有効なのが、並列計算である。近年のCPUの多くが、複数のコアを持つので、これらを使って同時に計算をすることができれば、原理的にはコアの数の分だけ速くなるはずである。Rでは、parallelというlibraryを使うと、比較的簡単に並列計算ができるらしい。

library(parallel)
cores <- detectCores()
cl<-makeCluster(cores,type="PSOCK")

として、コアの数だけ同時に計算できるようにして、clusterApply()やparSapply()を使って、リストやベクトルの計算をすれば良いのだが、それぞれのコアライブラリを組み込むためにclusterEvalQ()で実行したり、変数を認識させるためにclusterExport()を使ったりする必要がある。計算の準備などをするために、実際にはコアの数分の一にはならないが、それなりには速くなる。また、最後にはstopCluster(cl)としておく。

単純な計算の並列化には成功したが、関数の中でmulti-coreを使って計算しようとしたら、環境の部分などでトラブって、なかなかうまく行かない。Rを使うと比較的簡単に並列計算ができることが分かったが、いろいろと癖がありそうである。

Rで外積

Rにベクトルの計算をさせていたら、通常の外積が無いことに気がついた。とりあえず、適当に定義してみた。

"%v%"<-function(x,y)c(x[2]*y[3]-x[3]*y[2],x[3]*y[1]-x[1]*y[3],x[1]*y[2]-x[2]*y[1])

これでa%v%bとすれば、ベクトル外積が計算できる。名前は平行六面体の体積を計算するときに使うので、vにしてみた。oもxもすでに使われているし。

"%v%"<-function(x,y)x[c(2,3,1)]*y[c(3,1,2)]-y[c(2,3,1)]*x[c(3,1,2)]

の方がRチックで良いかな。

arduinoとLUFA それなりに昔のことになるが，書き込み機が簡単に自作できるということを知って，マイコンとしてAVRを使うようになった．当初はtinyを使っていたが，ピンの数が足りなくなって，megaも使うようになった．ハンダ付けや配線と書き込みが面倒だと感じて居たのだが，arduino自作用の基盤を使うと，それらが楽になることに気が付いて，しばらく使っていた．そして，中国製のarduinoが安価に手に入るようになったので，arduinoも使うようになった．arduinoとしては，UNOを最初に使ったが，nanoはそれとほぼ同等の性能で小型なので，小さい部分に使う時に使うようになった． USBと通信のできるnanoよりも小さくて16個以上のピンのものが必要になった．nanoの基板の片側を切り取れば，希望の大きさになり，その部分には重要な回路は無いので一応動作することが分かったが，面倒だし美しくない．いろいろと調べたら，arduino pro microというのが，それなりに小さくてピンの数も大丈夫だということが分かった．さらに，これはAVR本体でUSBと通信しているので，USBをシリアル以外としても使うことができる． arduino pro microは，arduino Leonardoを小型にしたもののようで，AVRとしては，atmega32u4を使っている．通常のarduinoとは違って，USB-serial変換ICが無いので，書き込みの癖が強いらしい．実際にやってみたが，うまく行ったり行かなかったりで，調子は良くない．このAVRについて調べていたら，LUFAというものを発見した．arduinoでは無いのだが，独自のbootloaderを作っており，今回の目的に一致するように感じたので，少し使ってみた． LUFAをdownloadしてunzipすると，いろいろなfileができるが，Bootloaders/MassStorageにあるBootloaderMassStorageを使うことにした．これは，AVRをharddiskのような形で認識させて，そこにファイルのような形でプログラムを転送して，AVRに実行させるものである．まず，makefileを編集して，以下のようにparameterを変更する． MCU = atmega32u4 BOARD = LEONARDO F_CPU = 16000000 FLASH_SIZE_KB = 32 BOOT_SECTION_SIZE_KB = 4 そして，BootloaderMassStorage.cを必要に応じて変更する．bootloaderが実行されると，watchdogリセットがかかってApplication_Jump_Check()が実行され，ここからプログラムに飛ぶ．このsubroutineのLEONARDOのところのピンの条件を変更すると，その条件が満たされたときには，プログラムに飛ばずにプログラム転送modeになる．makeすることによってhexファイルを作成し，ISPでpro microに書き込む．arduino as ISPを使う場合には，以下のコマンドである． avrdude -p m32u4 -c avrisp -P /dev/ttyUSB0 -b 19200 -B 4 -U flash:w:"BootloaderMassStorage.hex" -v ちなみに，fuseはE:CB, H:D8, L:FFとなっていた．BOOTRSTは0となっているが，1にすると直接プログラムに飛ぶのだと思う．必要に応じて，fuseも書き換えれば良い．プログラム転送modeになったときには，二つのLEDが交互に光る．この状態でAVRはmediaとして認識されるので，binファイルをddで転送する．cpでcopyしたらダメだった．sdbとして認識している場合には，以下のコマンドでファイル転送できる．もし，sdbが別のdiskを示している場合には，それを上書きしてしまわないように注意が必要である．書き込んだ後にumountすることもお忘れ無く． sudo dd if=output.bin of=/dev/sdb seek=4 arduinoを使って書けるようなプログラムは，その方が簡単に作ることができる．転送するプログラムのbinをarduinoで作る方法を説明しよう．まず，arduino IDEはaptで入れるのでは無く，本家から新しいものをdownloadしてinstallする必要がある．File - Preferences - Additional Boards Manager URLs:のところにhttps://raw.githubusercontent.com/sparkfun/Arduino_Boards/master/IDE_Board_Manager/package_sparkfun_index.jsonを指定して，Tools - Boards Managerからを選んで，pro microを使えるようにする．ToolsでBoardはSparkFunProMicroとして，適切な電圧とクロックを選択する．プログラムができたら，compileして，Sketch-Export compiled Binaryとすると，hexファイルができるので，以下のコマンドでbinファイルを作ると，これを上記の用に転送して再起動すればプログラムを動かすことができる． avr-objcopy -I ihex -O binary -R .

maximaで角運動量演算子
Rで角運動量演算子の行列を計算するプログラムは、以前に書いて、時々使っているが、Rでは変数を含んだ計算をすることができない。そこで、maximaを使ってこの行列を計算するプログラムを書いてみた。maximaについてはど素人なので、無駄な部分が多いと思うが、これを使うと変数を含んだ変形が簡単になるはずだ。ただし、行列の掛け算は.を、行列のn乗は^^を使うことに注意しないといけない。最後の行は、核四重極相互作用の行列を例として示した。

mjz(j):=genmatrix(lambda([a,b],if a=b then j+1-a else 0),2*j+1,2*j+1);
mjp(j):=genmatrix(lambda([a,b],if a=b-1 then sqrt(a*(2*j+1-a)) else 0),2*j+1,2*j+1);
mjm(j):=transpose(mjp(j));
mjx(j):=(mjp(j)+mjm(j))/2;
mjy(j):=(mjp(j)-mjm(j))/2/%i;
me(j):=ident(2*j+1);
3*mjz(j)^^2-j*(j+1)*me(j)+n/2*(mjp(j)^^2+mjm(j)^^2);

2020/4/15追記 久々に使おうとしたら、代入がコロンだということを忘れていて、なかなか動かなくて苦労した。