ページ

2016年8月30日火曜日

Qittaへ

最近はマークダウンで書くことが多くなりましたが,Bloggerでもできなくはないけど不便.
お世話になっているQittaに浮気しようと思います.

私の投稿プログラムに関することだし.

時々はBloggerの事も思い出してあげようと思います.

Hello Qitta, Goobye Blogger

2016年8月3日水曜日

ubuntu kernel update + nvidia = black screen になる場合

備忘録

nvidia のグラフィックボードを使っている場合,ubuntu のkernel updateの後にディスプレイは認識されているもののblack screenになる事がある.
nvidia 以外に原因は考えられないという事で対処法の備忘録


詳細はUbuntuでNVIDIAのディスプレイドライバが動作しない場合のチェック項目が非常に詳しく助かりました.

私の場合の症状は起動後メーカーロゴが表示,その後待てども暮らせどログイン画面が表示されないttyも表示されない,でもsshでログインできる状態.

主な原因はkernel updateでnvidiaのドライバーが重複して存在している事の様です.対処療法は以下の通り.

$ ls /lib/modules/カーネルバージョン-generic/updates/dkms/
nvidia_340.ko  nvidia_340_uvm.ko  nvidia.ko  nvidia-uvm.ko

と見るとと複数のnvidiaのドライバが混在しているので

$ modinfo nvidia
filename:       /lib/modules/カーネルバージョン-generic/updates/dkms/nvidia.ko
alias:          char-major-195-*
version:        340.96

よりバージョンを確認.不要なnvidia_340.koとnvidia_340_uvm.koを削除(もしくはホームディレクトリなどに退避)

lightdmを再起動しても認識されないので,rebootを掛ける.
待てども暮らせどshutdown されないので,電源長おしで強制終了を掛ける.

無事login画面が表示される.

うーん忙しい時にこの手のトラブルはご勘弁願いたい.
何も考えずにupgradeを実行するのはよくないのかなぁー

sphinx + ipython notebook (jupyter)

備忘録

ipython notebookで作ったファイルをipython nbconvert --to rst hoo.ipynb
すると
In[1]
Out[1]
などの情報が出力されずに困った.

http://nbsphinx.readthedocs.io/en/latest/usage.html

よりrstにconvertせずにipynbをそのままtoctreeで読み込む事によりこの問題を回避

sphinxの設定はデフォルトの場合,conf.pyがあるディレクトリで

python -m sphinx . _build/html

とすればよい

sphinx-autobuild が超便利

以上

2016年6月29日水曜日

4倍精度対角化

Lapackの前身にあたるeispackという一般複素行列の対角化を実行するfortrunのプログラムが有ります.3年ほど前に4倍精度に拡張して,別のプログラムに間接的にくっつけてgithubにあげていたんですが,需要があるらしい(?)ので単体でgithubにupしました.

(追記)gcc4倍精度に関するコーディングは昔のポストが参考になれば幸い

拡張と言っても宣言とリテラルを書き換えた位だった気がする.
もう昔のことで覚えていない.

今はmathematicaで多倍長対角化を実行していて,真面目にテストしてないのでバグが絶対混入していると思いますが,エルミートとユニタリーのクラスの対角化では,意図した計算結果を示すのでとりあえず良しとする.

将来的には任意精度まで拡張したいと考えていますが,速度がどうなることやら.
いやー,4倍精度演算がsoftwareレベルで実装されているので,計算速度が遅くなるのは仕方ないんだけど,やっぱりlapackめちゃくちゃ早くて素晴らしいな

2016年5月15日日曜日

クラウドストレージBoxをubuntu 16.04で使ってみる

Dropboxユーザー(無料)だが容量が足りない
それでも,ベータテスト参加した際の容量アップで9.8GB容量が有るのだが,
限界に近づいては圧縮を繰り返すという非生産的活動に時間を費やしたくない.

と言って有償版(2TB, 1200円/月)を使う程でもない.
50GB 300円/月とかあれば使うのだが・・・(昔はあった気がする)

実はCopy.comというDropboxと同等のサービスを使っていたのだが,
案の定(名前が悪いと思う.SEOの意味で)終了してしまったので,
boxなるサービスを使ってみる.

linuxクライアントが無いので若干のひと手間が必要.
一応dropboxと同様に使えるようになった.
いや正確ではない.このことは後述する
以下その備忘録

1. boxサイトからアカウントを作る

電話番号を要求されるので,若干の抵抗感がある
私の場合デフォルトで10GBまで無料
有料版は100GB で 1200円/月とDropboxより高めな印象

2. davfsをインストールする

$ sudo apt-get install davfs2
root以外でも使えるようにするか効かれるのでyesとする.
noを選択した場合
$ sudo dpkg-reconfigure davfs2
から再設定可能.

3. 何はともあれmountしてみる

適当にmount point
$ mkdir ~/Box
を作って,mountする
$ sudo mount -t davfs https://dav.box.com/dav /home/hanada/Box/
username (メールアドレス)とパスワードを入力してmountできたらとりあえず良し.

ただ,root権限でmountしておりユーザーは書き込みできないので修正する.
ここを参考にそのまんま設定すると上手くいったので詳細は省く.

最大ファイルサイズが250MBなのでデカイファイルのバックアップには向かない.
とりあえず,文献データ(5GB)を置いてみると,10MB/sec程度でコピーしてくれる.

課題と問題点

Dropboxよりえらく早くファイルコピーするなーと思って,よくよく考えたら
dropboxはローカルストレージとクラウドストレージ間の双方向同期でありoffline下でも使える.
更に異なるマシーン間の競合も上手く処理してくれる.

しかし上記方法はクラウドストレージのマウント方法なのでoffline下ではそもそもマウントできない.
この違いは大きい,というかDropboxのようには使えない.

offline環境下でも使えるようにするためには,別に母艦を用意しておいて
rsyncなど適当に設定するか?rsyncで双方向同期できるのか?unisonか?
マシーン毎に設定が必要か?
win mac は box syncなるアプリが有るらしいがwineでも試してみるか?
うーん,凝りだすと時間かかりそうで面倒やな・・・

2016年5月7日土曜日

ubuntu 16.04 へアップグレードした際の備忘録


linuxマシーンを何台かubuntu 16.04にする予定なので備忘録を残す.

なお,*マークが有る箇所はhomeを引き継ぎで実際には設定不要だった項目

前置き

ハードウェアーはthinkpad x1 carbon 3rd
ubuntu 15.10 install済みからubunt 16.04(日本語remix)のクリーンインストール
パーティション関係は
/
/home
swap
を分離

・homeをバックアップ
・usb bootable diskはstartup disk creatorから制作
・(NG): ubuntu 15.10を削除して16.04をインストール的なものを選択
   これをやると /のパーティション上にhomeがつくられて
   旧 /home パーティションが宙ぶらりんになっていた^^;
・インストール画面でその他を選びマウントポジションを指定して再インストール(OK)
無事終了

インストール自身は5分くらいで終わった.
昔はアップデート一日仕事だったけど3時間くらいに短縮されたなー
パーティション関係の失敗がなければ30分くらいで終わりそう
次のOS update 5年後になることを期待

以下インストール後の設定備忘録

いつもお世話になっているubuntu インストール直後に行う設定を参考に必要な設定を行う
wi-fiは非表示のネットワークから接続を構成しないと何故か繋がらない

システムの言語環境を英語に,特定のアプリを日本語環境で立ち上げる

ubuntuは日本語remix版をinstallしたが(オリジナル版はzipファイルに日本語名があると文字化けするらしい)が,ターミナルの出力が日本語だと,なんかエラーが出た際にstackoverflowなどでの検索に不都合.
またプログラムを書く際に,メニュー等やエラーを英語でどのように掛けばよいか悩む場合があるので,日頃から英語環境にしておく.

システム環境設定から言語サポート(?),使用する言語(English)を最上部に設定.

ただこれだと全てのアプリケーションが英語環境で立ち上がるので時に不便.

特定のアプリケーションを日本語でたちあげたい場合はlocale -a で ja_JP.utf8があることを確認して
$ locale -a | grep ja
ja_JP.utf8

$ env LANG=ja_JP.utf8 application_name
みたいな感じて起動すれば日本語環境で起動可能.メーラーとかeclipseとか設定ファイルが巨大なアプリは日本語で開いたほうがミスが少なくてすむ.

日常的に使うアプリはショートカットを自作し設定.
なお,現在の言語環境の設定はターミナルより確認可能.
$ locale
LANG=en_US.UTF-8
LANGUAGE=en
LC_CTYPE="en_US.UTF-8"
LC_NUMERIC=ja_JP.UTF-8
LC_TIME=ja_JP.UTF-8
LC_COLLATE="en_US.UTF-8"
LC_MONETARY=ja_JP.UTF-8
LC_MESSAGES="en_US.UTF-8"
LC_PAPER=ja_JP.UTF-8
LC_NAME=ja_JP.UTF-8
LC_ADDRESS=ja_JP.UTF-8
LC_TELEPHONE=ja_JP.UTF-8
LC_MEASUREMENT=ja_JP.UTF-8
LC_IDENTIFICATION=ja_JP.UTF-8
LC_ALL=

キーボード入力関係*

・macっぽい日本語/英語の切り替えにする
Mozc(fcitx)/入力メソットの設定/ホットキー
変換, カタカナひらがな:入力メソットオン
無変換:入力メソットオフ
全角半角: オフ

・句読点をコンマ,ピリオドに
  dash からmozcの設定を起動.基本設定から変更
・ここを参考にCuplock をctrlに

覚書:fcitx-mozcの辞書(ユーザー辞書)でコメント内容を見るためには変換候補を選んだ際に(Ctrl+Alt+H)でコメント内容を表示させる事が可能

touch padのon/offショートカット*

$ xinput list
⎡ Virtual core pointer                     id=2 [master pointer  (3)]
⎜   ↳ Virtual core XTEST pointer               id=4 [slave  pointer  (2)]
⎜   ↳ SynPS/2 Synaptics TouchPad               id=11 [slave  pointer  (2)]
⎜   ↳ TPPS/2 IBM TrackPoint                    id=12 [slave  pointer  (2)]
でidを調べて
$ xinput set-prop 11 "Device Enabled" 0 # disable
$ xinput set-prop 11 "Device Enabled" 1 # enable
あとはショートカットを設定

tex関係

ここを参考にパッケージをインストール

$ sudo apt-get install texlive texlive-lang-cjk xdvik-ja dvipsk-ja gv texlive-fonts-recommended texlive-fonts-extra
$ sudo apt-get install texlive-publishers
2番目はrevtexとのスタイルファイルが入っている.日本語及びrevtexは問題なく動く.

python 関係*

anaconda3をinstall

ノートのバッテリの電力制御

充電するタイミングを制御する為tlpをインストールする.
$ sudo add-apt-repository ppa:linrunner/tlp
$ sudo apt-get update    
$  sudo apt-get install tlp tlp-rdw  
thinkpad userなので別途
$  sudo apt-get install tp-smapi-dkms acpi-call-dkms 
をinstall.バッテリー情報を
$  sudo tlp-stat -b 
より確認.電源プラグを接続した際に充電するかどうかの制御は
$  sudo tlp setcharge 30 80 BAT0
のようにすれば良い.意味はバッテリーレベルが30%になると充電を開始し,バッテリーレベル80%まで充電する.

プリンター関係 (canon LBP5000を動かす)

Canonからlinux用のinstallerをダウンロードして install.shを実行
The install procedures are the following:
extract the downloaded tar file
  $ tar zxvf filename.tar.gz
The containing files in the extracted directory are
  $ ls
  64-bit_Driver  Documents  install.sh
In this directory, run
  $ sudo bash install.sh
then type [y] to answer some questions.
ここで
 OK: libglade2-0
 OK: libstdc++6:i386
 OK: libxml2:i386
 OK: libjpeg62:i386
 NG: libbeecrypt7:i386
 NG: libbeecrypt-dev:i386
を入れろと言われるのでOKはaptからinstall.NGはpackageが見つからず.

Then
 $ sudo lpadmin -p LBP5500 -P /usr/share/cups/model/CNCUPSLBPCSERIESLJ.ppd -v lpd://192.168.64.68/LBP5500 -E

That's all for printer setting. 

一応使えるので良しとする.

Evernote クライアント

nixnote  (evernote linux クライアント)を入れる.
debは動かなかったのでtarからinstall.sh を実行.
起動するとopencvが無いと怒られる
ibopencv-dev (もしかしたらlibopencv-core-devも必要かも)をインストール
syncからevernoteにloginすればok

Ricty fontのインストール

ここを参考

その他install した(する) packages

  • dropbox:公式の配布を利用
  • unity-tweak-tool: UIの拡張に
  • emacs, yatex, emacs-mozc
  • conky (gcal, lm-sensors, acpi)
  • gnuplot5-x11
  • duplicity :  backup tool
dejadup-monitorが暴走しなくなったので,システムbackupで十分かも.遅いけど
  • gparted: パーティションエディタ
  • zotero  : ここからPPAを利用した
  • gufw (+ python-gi) :ファイヤーフォール
  • eclipse: IDE
  • geany: IDE 
  • kate: 
editorは何が良いかいつも悩むんだけど,用途によって使い分けています. 
  1. 中,大規模のプログラミング: eclipse
  2. 小規模,使い捨てプログラミング,tex(スペルチェック込): emacs
  3. プレゼンスライド作成(reveal.js, html, 要スペルチェック): kate
geditではなくkateを使う理由はショートカットを設定が楽なので.emacsのキーバインド(Ctrl-a= Home button, Ctrl-e: End button)を恒常的に使用しているので必須. 
  • okular: pdf viewer
  • gfortran: fortrun コンパイラー
  • nkf: 文字コードの変換
  • libav-tools: (= ffmpeg) 動画の変換,作成に
  • pandoc: (ipython convertで必要)
  • gnome-sushi : nautilus でスペースキーを押すとサムネイルをpopup表示
(gloobus previewのほうが良い感じ).macのquick lookのような機能.epsだろうがsvgだろうが表示させることができる.
透明画像が暗転してしまうが・・・
  • gloobus preview: 複数枚のファイルをquick lookできる
透明画像が暗転するが,preview中に左右のキーを押せば次のファイルをpreviewするので非常に便利
参考: http://www.webupd8.org/2015/01/quick-file-previewer-gloobus-preview.html 

  • php:
php5からphp7になり,ブラウザからアクセスしてもソースコードが表示されるのでここを参考に再設定

 

Mathematica Benchimark

Benchmarkingパッケージを使ったMathematicaのベンチマーク

私物laptop
Intel(R) Core(TM) i7-5500U CPU @ 2.40GHz
# Core: 2
# Threading: 4
{"MachineName" -> "--", "System" -> "Linux x86 (64-bit)",
"BenchmarkName" -> "WolframMark", "FullVersionNumber" -> "10.3.1",
"Date" -> "January 31, 2016", "BenchmarkResult" -> 1.392,
"TotalTime" -> 9.945, "Results" -> {{"Data Fitting", 0.537},
{"Digits of Pi", 0.341}, {"Discrete Fourier Transform", 0.68},
{"Eigenvalues of a Matrix", 0.608}, {"Elementary Functions", 0.567},
{"Gamma Function", 0.467}, {"Large Integer Multiplication", 0.467},
{"Matrix Arithmetic", 0.334}, {"Matrix Multiplication", 0.754},
{"Matrix Transpose", 1.61}, {"Numerical Integration", 0.676},
{"Polynomial Expansion", 0.092}, {"Random Number Sort", 1.16},
{"Singular Value Decomposition", 0.765}, {"Solving a Linear System", 0.887}}}

------------------------------------------------------------------------------------------------
とあるサーバー1:
Intel(R) Xeon(R) CPU X5550 @ 2.67GHz
# Core: 4
# Threading: 8
{"MachineName" -> "--", "System" -> "Linux x86 (64-bit)",
"BenchmarkName" -> "WolframMark", "FullVersionNumber" -> "10.0.2",
"Date" -> "January 31, 2016", "BenchmarkResult" -> 1.332,
"TotalTime" -> 10.388, "Results" -> {{"Data Fitting", 1.059},
{"Digits of Pi", 0.507}, {"Discrete Fourier Transform", 0.545},
{"Eigenvalues of a Matrix", 1.001}, {"Elementary Functions", 0.389},
{"Gamma Function", 0.66}, {"Large Integer Multiplication", 0.668},
{"Matrix Arithmetic", 0.536}, {"Matrix Multiplication", 0.408},
{"Matrix Transpose", 1.022}, {"Numerical Integration", 1.002},
{"Polynomial Expansion", 0.136}, {"Random Number Sort", 1.159},
{"Singular Value Decomposition", 0.651}, {"Solving a Linear System", 0.645}}}

------------------------------------------------------------------------------------------------
とあるサーバー2:
Intel(R) Xeon(R) CPU E5607 @ 2.27GHz
# Core: 4?8?
# Threading: 8
{"MachineName" -> "--", "System" -> "Linux x86 (64-bit)",
"BenchmarkName" -> "WolframMark", "FullVersionNumber" -> "10.0.2",
"Date" -> "January 31, 2016", "BenchmarkResult" -> 1.056,
"TotalTime" -> 13.103, "Results" -> {{"Data Fitting", 0.962},
{"Digits of Pi", 0.648}, {"Discrete Fourier Transform", 0.6},
{"Eigenvalues of a Matrix", 1.076}, {"Elementary Functions", 0.481},
{"Gamma Function", 0.874}, {"Large Integer Multiplication", 0.916},
{"Matrix Arithmetic", 0.706}, {"Matrix Multiplication", 0.541},
{"Matrix Transpose", 1.485}, {"Numerical Integration", 1.342},
{"Polynomial Expansion", 0.14}, {"Random Number Sort", 1.542},
{"Singular Value Decomposition", 0.878}, {"Solving a Linear System", 0.912}}}

------------------------------------------------------------------------------------------------
Model Name: iMac
Processor Name: Intel Core i7
Processor Speed: 3.1 GHz
# Core: 4
# Threading: 8

{"MachineName" -> "--", "System" -> "Mac OS X x86 (64-bit)",
"BenchmarkName" -> "WolframMark", "FullVersionNumber" -> "10.3.1",
"Date" -> "January 31, 2016", "BenchmarkResult" -> 1.985,
"TotalTime" -> 6.975, "Results" -> {{"Data Fitting", 0.689},
{"Digits of Pi", 0.295}, {"Discrete Fourier Transform", 0.369},
{"Eigenvalues of a Matrix", 0.474}, {"Elementary Functions", 0.506},
{"Gamma Function", 0.353}, {"Large Integer Multiplication", 0.33},
{"Matrix Arithmetic", 0.632}, {"Matrix Multiplication", 0.259},
{"Matrix Transpose", 0.726}, {"Numerical Integration", 0.62},
{"Polynomial Expansion", 0.086}, {"Random Number Sort", 0.819},
{"Singular Value Decomposition", 0.435}, {"Solving a Linear System", 0.382}}}

------------------------------------------------------------------------------------------------
Model Name: iMac
Processor Name: Intel Core 2 Duo
Processor Speed: 3.06 GHz
# Core: 2
# Threading: 2(?)

{"MachineName" -> "--", "System" -> "Mac OS X x86 (64-bit)",
"BenchmarkName" -> "WolframMark", "FullVersionNumber" -> "10.0.2",
"Date" -> "January 31, 2016", "BenchmarkResult" -> 0.613,
"TotalTime" -> 22.586, "Results" -> {{"Data Fitting", 1.713},
{"Digits of Pi", 0.67}, {"Discrete Fourier Transform", 1.548},
{"Eigenvalues of a Matrix", 1.356}, {"Elementary Functions", 2.113},
{"Gamma Function", 0.768}, {"Large Integer Multiplication", 0.769},
{"Matrix Arithmetic", 2.632}, {"Matrix Multiplication", 1.646},
{"Matrix Transpose", 1.804}, {"Numerical Integration", 1.949},
{"Polynomial Expansion", 0.267}, {"Random Number Sort", 1.669},
{"Singular Value Decomposition", 1.992}, {"Solving a Linear System", 1.69}}}

-----------------------------------------------------------------------------------------------
raspberry pi 2
Coretex-A7
# Core: 4
# threading: 4(?)
{"MachineName" -> "--", "System" -> "Linux ARM (32-bit)",
"BenchmarkName" -> "WolframMark", "FullVersionNumber" -> "10.0.2",
"Date" -> "January 24, 2016", "BenchmarkResult" -> 0.029,
"TotalTime" -> 475.891,
"Results" -> {{"Data Fitting", 10.904}, {"Digits of Pi", 4.805}, {
"Discrete Fourier Transform", 30.977}, {
"Eigenvalues of a Matrix", 35.914}, {
"Elementary Functions", 25.479}, {"Gamma Function", 6.437}, {
"Large Integer Multiplication", 6.735}, {
"Matrix Arithmetic", 6.07}, {"Matrix Multiplication", 151.389}, {
"Matrix Transpose", 7.202}, {"Numerical Integration", 11.684}, {
"Polynomial Expansion", 1.086}, {"Random Number Sort", 8.376}, {
"Singular Value Decomposition", 75.031}, {
"Solving a Linear System", 93.802}}}

----------------------------------------------------------------------------------------------------
とあるサーバー

Intel(R) Xeon(R) CPU E5-1680 v3 @ 3.20GHz

{"MachineName" -> "--", "System" -> "Linux x86 (64-bit)",
 "BenchmarkName" -> "WolframMark", "FullVersionNumber" -> "10.4.1",
 "Date" -> "June 1, 2016", "BenchmarkResult" -> 2.06, "TotalTime" -> 6.721,
 "Results" -> {{"Data Fitting", 0.568}, {"Digits of Pi", 0.293},
   {"Discrete Fourier Transform", 0.274}, {"Eigenvalues of a Matrix",
    0.78}, {"Elementary Functions", 0.209}, {"Gamma Function", 0.395},
   {"Large Integer Multiplication", 0.389}, {"Matrix Arithmetic", 0.214},
   {"Matrix Multiplication", 0.282}, {"Matrix Transpose", 0.611},
   {"Numerical Integration", 0.705}, {"Polynomial Expansion", 0.098},
   {"Random Number Sort", 0.984}, {"Singular Value Decomposition", 0.589},
   {"Solving a Linear System", 0.33}}}

2016年5月6日金曜日

fortran 4倍精度で,計算速度はどの程度遅くなるか?

大学の環境で利用可能なコンパイラー
gfortran, ifort (intel), xl fortran (IBM), PGI fortran (NVIDIA)で
倍精度から4倍精度にした際にどの程度計算速度が遅くなるか比較する.
なお pgi fortranでは4倍精度のプログラムがコンパイルできなかったので,
実際には残りの3つの比較.

計算内容は割り算と掛け算をひたすら繰り返しただけ.

最適化(オプション)なしの比較.
----------------------- 結果1 (intel cpu)--------------------------------------

gfortran quadruple/double ratio
real: 9.6 倍遅延
complex : 19倍遅延

ifort quadruple/double ratio
real: 16倍 遅延
complex: 28倍遅延

efficiency (gfortran calc. time/ ifort calc. time)
gfortran -> ifort
real: 5.988sec/4.963sec -> 1.2倍 高速化
complex 22sec/16sec  ->1.3倍 高速化
(ifortの方がちょっとだけ早い)

--------------------------結果2 (ibm ワークステーション)------------------------

xlf90 quadruple/double ratio
real: 2.27倍遅延
complex: 5.83倍遅延

efficency  (gfortran -> xlf90)
gfortran 4.7.3で4倍精度コンパイルできなかったので比較できず.
(4.6から対応だと思っていたが・・・)

pgf90 4倍精度コンパイルできず...

-----------------------------------------------------------------------------
上記比較はオプションなしの場合.

ifortは倍精度は2-3倍オプションなしで早くなるが,
4倍精度では思ったより実力を発揮しなかった.
恐らく4倍精度はソフトウェア上での実装なので,最適化があまり意味がないのかと思う.
ただし,ifortの場合do文を自動並列化ができるので強制的に並列化してやれば,最大でcpu(thread)数倍くらい早くならないこともない.


ibmのコンパイラーはibmのマシーンに特化してる気もする.
他にもNAGや富士通(販売終了),absoftがあるらしいが,使ったこと無い.

価格2016年5月時点
IBM: $3380 
NAG: (63,000円) single user license 
Intel: 10-30万くらい (Parallel Studio XE Composer Edition)
PGI:5-10万くらい
absoft:6万-12万 
備考:intel compiler の非商用版(linuxのみ)は無償だったが,有償版だけになった模様.
学生に限りは無償使用のライセンスが在る模様.

2016年1月24日日曜日

Mathematica ベンチマーク @ Raspberry pi 2

Raspberry pi 2 を買ったのでMathematica(無償)のベンチマークを測ってみる。
mathematicaはnoobs から rasbian (jessei) installでデフォルトで導入されていた。


--------------- 参考 (手元のノート core-i7)-----------------------

In[1]:= Needs["Benchmarking`"]

In[2]:= Timing@Benchmark[]

Out[2]= {14.596, InputForm[{
  "MachineName" -> "----", "System" -> "Linux x86 (64-bit)", 
   "BenchmarkName" -> "WolframMark", "FullVersionNumber" -> "10.0.0", 
   "Date" -> "January 24, 2016", "BenchmarkResult" -> 1.197, 
   "TotalTime" -> 11.567, 
   "Results" -> {{"Data Fitting", 0.629}, {"Digits of Pi", 0.536}, {
     "Discrete Fourier Transform", 0.728}, {
     "Eigenvalues of a Matrix", 0.935}, {
     "Elementary Functions", 0.774}, {"Gamma Function", 0.456}, {
     "Large Integer Multiplication", 0.46}, {
     "Matrix Arithmetic", 0.416}, {"Matrix Multiplication", 0.983}, {
     "Matrix Transpose", 1.468}, {"Numerical Integration", 0.62}, {
     "Polynomial Expansion", 0.092}, {"Random Number Sort", 1.104}, {
     "Singular Value Decomposition", 1.268}, {
     "Solving a Linear System", 1.098}}}]}



-------------- Raspberry pi 2 (オーバクロック無し) ----------------------

In[1]:= Needs["Benchmarking`"]

In[2]:= Timing@Benchmark[]

Out[2]= {551.620000, InputForm[{
  "MachineName" -> "----", "System" -> "Linux ARM (32-bit)",
   "BenchmarkName" -> "WolframMark", "FullVersionNumber" -> "10.0.2",
   "Date" -> "January 24, 2016", "BenchmarkResult" -> 0.029,
   "TotalTime" -> 475.891,
   "Results" -> {{"Data Fitting", 10.904}, {"Digits of Pi", 4.805}, {
     "Discrete Fourier Transform", 30.977}, {
     "Eigenvalues of a Matrix", 35.914}, {
     "Elementary Functions", 25.479}, {"Gamma Function", 6.437}, {
     "Large Integer Multiplication", 6.735}, {
     "Matrix Arithmetic", 6.07}, {"Matrix Multiplication", 151.389}, {
     "Matrix Transpose", 7.202}, {"Numerical Integration", 11.684}, {
     "Polynomial Expansion", 1.086}, {"Random Number Sort", 8.376}, {
     "Singular Value Decomposition", 75.031}, {
     "Solving a Linear System", 93.802}}}]}


線形代数の差がすごいな。topで観察する限りElementary Functions の計算でcpu使用率が400%になる。

Raspberry pi 2 は4コアなので、フル性能使えているんだと思うけど、並列計算用karnelの起動できない。(LaunchKarnels[]の実行結果が空)
(ParallelDoかなんか実行た際に並列計算用のパッケージ?みたいなのをダウンロード後使えるようになった)

一方面白いことに、有料版Mathematicaはライセンスの制約上フロントエンドカーネルは2つまで(mathematicaを独立に2つまで)しか起動できないが、
Raspberry pi の場合はその制限が無いようだ。

他に気づいたことは環境設定が無いといこと、MathematicaScript、MatheKarnelが無いこと・・・(残念)

追記:代わりにwolframというコマンドがあるっぽい

Raspberry 版のMathematicaは何かしら機能制限がありそうだと思ってはいたけどちゃんと確かめてなかったので他にも制約がありそうです。
本当はラズパイ並列して、mathematicaタスク分散とか・・・、もう少し遊べそうかと思ったけど無理そうだな。

まー遊び道具としてもうちょっといじってみるか。

速度の問題はあるけど5000円程度でmathematicaが使えると思うと非常に安いのでは?
Mathematicaは日本人特別価格の設定でHome edition で7万くらい(日本以外は$300,日本から海外版購入不可)だからね.

-----------
raspberry piを安価pcとして使用した際に気づいたこと。

microsdの相性の問題で起動しなかったので、調べて買うこと。
wikiにまとまっているのでそれを買えば問題ないはず。


os がarmなのでubuntuのようにパッケージが豊富でないこと

例えばfirefox、やchromeがaptになく代わりにiceweaseelやmidoriがある。

chromium-browserは一つ前のversionのosではあったようだけど、最新版ではaptに無いため、debを拾ってきて インストールする必要がある。
Adobe flashが無いので困る人には困るかもしれない。

操作性は思ったより遅くはないけど機敏ではないかな
(El captanよりは応答が良いと思うのは私だけだろうか)
動画を見るのは結構大変かも?
オーバクロックのturboモードはめちゃくちゃ不安定だったけどpi2は割と安定してそう。

2015年6月4日木曜日

LaTeXでpdf画像が出力されない場合

日本語LaTeX (platex + dvipdfmx) 環境でpdf画像を\includegraphicsする際に,画像が出力されない場合があって困っていた.
具体的にはdviでは画像が出力されdvi -> pdfに変換後見えなくなる.

数年前に同じ経験をしていたのだが,解決方策を完全に忘れていたので,備忘録として残す.

私の環境下での解決方策は単純でpdfのversionが1-5だとdvi -> pdfに変換した際に画像が出力されなくなる.適当なテキストエディタ(emacsとか)でpdfの1行目にあるversionの情報
%pdf-1.5
を
%pdf-1.4
に書き直すだけで問題が解決された.副作用があるかどうかは知らない.

個人的な疑問だが,epsは原理的にテキストエディタで編集可能である事は知っていた.
その為nkf -g で文字コードを調べるとASKIIと帰ってくる.
しかしpdfはnkf -g で文字コードを調べるとBINARYである.
BINARYファイルを普通のテキストエディタで編集しても良いのだろうか?




Blog書くのひさしぶりだから(?),単文でも時間がかかるなぁ・・・.



2014年5月27日火曜日

MathemathicaScript の引数

MathematicaScriptの引数の引き渡しをいつも忘れてしまうので覚え書き.

MathematicaScriptの引数は

 $ScriptCommandLine
渡され,一番目はScriptのファイル名,2番目から引数となる.引数の制御は

If [Length[$ScriptCommandLine] !=2, 
 Print["usage: MathematicaScript -script" <> $ScriptCommandLine[[1]] <> " value1"];
 Exit[1] ]

とでもしておけば良い.引数はstringで受け取るので,実数とかにしたければ
ToExpression@$ScriptCommandLine[[2]]

とでもすればよい.

2014年4月29日火曜日

home page

home pageのコンテンツを一次非公開にしました.理由は内容が微妙だから.
また記事がまとまったら公開するかもしれません.

2014年1月22日水曜日

Mavericksのメモリー圧縮機能について

Mac OSX 10.9.1 (Mavericks) でカーネルパニックが発生したので記録する.

1.環境と問題

iMac
21.5-inch, Late 2013
プロセッサ  3.1 GHz Intel Core i7
メモリ  16 GB 1600 MHz DDR3


で数値計算を行った.メモリー使用量は

数百MB程度のメモリーを使用×6 + (瞬間的に)1GB程度メモリーを使用=最大で2GB程度

マシン性能から考えて余裕だろうとたかをくくっていたが,
結局2度強制終了するはめになった.

2. 現象

* 一定時間経過後,スワップアウトし始める
* CPUの使用率が100%にならない(10%)程度
* 最終的に応答不能状態に陥る->強制終了
アクティビティーモニターのメモリー観察
* kernel_task が10GB程度メモリーを占拠(何もしない状況では800MB程度)
* メモリー圧縮が8GB程度
* メモリープレッシャーが増大しつづける.






Mavericksから導入されたCompressed memory機能が怪しい!

3.解決

以下自己責任で
Disable compressed memory in Mac OS 10.9 Mavericks?
を参考にターミナルでメモリー圧縮modeを確認する(defaultは4)

$ sysctl -a vm.compressor_mode
vm.compressor_mode: 4

1でメモリー圧縮機能を切れるらしいので

$ sudo nvram boot-args="vm_compressor=1"

にして,再起動(なぜか再起動後も4だったので2回もやってしまった・・・)

$ sysctl -a vm.compressor_mode
vm.compressor_mode: 1

で vm.compressor_mode: 1 となっている事を確認する.

今のところ再度同じ計算をさせてみると,メモリープレッシャーも起こらないし,
kernel_taskがメモリーを馬鹿食いする事もなくなった.

メモリープレッシャーはマシン性能が悪い場合に威力を発揮するらしいが,
数値計算する上では邪魔な機能だなと思いました.

 おまけ

今まで知らなかったのだがターミナルで

$ sudo purge

を実行するとファイルキャッシュがクリアされる.

なんでも,duを行ってもファイルキャッシュが一掃されるらしいのだが,
効果は不明だし時間がかかるのでpurgeで十分な気がする
(少なくともmavericksでは).

2013年12月1日日曜日

Mathmaticaを用いた数値対角化の計算時間

Mathematicaを用いた対角化の計算速度をplotしてみました.

対角化した行列のclassは
ランダムではなくスパーズではない一般のエルミート行列(赤ドット)と
ランダムではなくスパースではない一般のユニタリー行列(青ドット).
また,行列は200桁精度で近似評価されています.

厳密評価のまま計算するとキットモット遅い.

行列の構成方法は気にしないと言う事でご勘弁.

実線は計算時間を3次の多項式でfittingしたもの.
最後に描画とfittingのプログラムを置いときます.
(もちろん計算時間のデータファイルがないと動きませんが・・・
欲しいなら上げますが,欲しいか?)




import numpy as np
import matplotlib.pyplot as plt

def f(x,para):
    return para[0]*x**3 + para[1]*x**2 + para[2]*x + para[3]
    
data = np.loadtxt("time_prec200.dat").transpose()
fig = plt.figure(figsize=(8,4))
ax0 = fig.add_subplot(1,2,1)
ax1 = fig.add_subplot(1,2,2)
colors=["red", "blue"]
for i, ax in enumerate([ax0,ax1]):
    ax.plot(data[0],data[1]/3600,'o', label='Hermitian matrix',markersize=10, color=colors[0])
    ax.plot(data[0],data[2]/3600,'o', label='Unitary matrix',markersize=10,color=colors[1])
#    ax.plot(data[0],data[3]/3600,'og')
    for j in range(1,3):
        x = np.arange(100, 10000,1)    
        para = np.polyfit(data[0],data[j],3)
        y = f(x,para)
        ax.plot(x,y/3600,'-',linewidth=3,color=colors[j-1])
    if i==1:
        ax.loglog()
        ax.set_xlim(50,10000)
    else:
        ax.set_ylim(0,5)
        ax.set_xlim(0,1000)
        ax.set_ylabel("calculation time (hours)",fontsize=16)            
    le = ax.legend(loc=2,numpoints=1)
    le.draw_frame(False)
    ax.grid()
    ax.set_xlabel("matrix dimension",fontsize=16)
sup = plt.suptitle("diagonalization time with Mathematica (precision=200)",fontsize=15)
sup.set_position((0.5,1.0))
plt.tight_layout()
plt.savefig("calc_time.png")
plt.show()

2013年11月21日木曜日

Maximaによる任意精度計算

Mathematicaはいつ使えなくなるか(ライセンス的な関係で)分らないので,Maximaを任意精度数値計算に使えないか試してみる.


精度(precision)のコントロールは

fpprec:50

演算による多倍長浮動小数点の有理数への変換を制御するため

bftorat: true

にする必要がある.デフォルトではfalse.
falseの場合多倍長浮動小数点は小さい有理数で近似される.

また,出力する桁のコントロールは

fpprintprec:50;

で可能.
maxima 5.30.0 だとload("fft")が出来なかったので,古いバージョンを使う.
参考:http://maxima.sourceforge.jp/maxima_5.html
$ maxima
Maxima 5.21post http://maxima.sourceforge.net
using Lisp SBCL 1.0.39
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) fpprec;
(%o1) 16
(%i2) fpprec:50;
(%o2) 50
(%i5) bftorat:true;
(%o5) true
(%i6) x:bfloat([0,1/3,0,0]);
(%o6) [0.0b0,3.3333333333333333333333333333333333333333333333333b−1,0.0b0,0.0b0]
(%i7) load("fft");
(%o7) "/Applications/Maxima.app/share/maxima/5.21post/share/numeric/fft.lisp"
(%i10) res:inverse_fft(fft(x));
(%o10) [0.0,.3333333333333333,0.0,0.0]
(%i11) bfloatp(res[1]);
(%o11) false

FFTは任意精度計算やってくれないっぽい.一方で行列mat=((i,pi),(pi,i))の対角化は任意精度でやってくれるっぽい.上のmatの固有値は(±pi, i)である.

(%i12) mat:bfloat(matrix([%i,%pi],[%pi,%i]));
(%o12) matrix([%i,3.1415926535897932384626433832795028841971693993751b0],[3.1415926535897932384626433832795028841971693993751b0,%i])
(%i13) load("eigen");
(%o13) "/Applications/Maxima.app/share/maxima/5.21post/share/matrix/eigen.mac"
(%i14) res:eigenvalues(mat);
`rat' replaced -9.8696044010893586188344909998761511353136994072408B0 by -41494057121218467653886768/4204226981644629723547549 = -9.8696044010893586188344909998761511353136994072408B0
(%o14) [[(4204226981644629723547549*%i−4*sqrt(10903152157933135742303986703335184117163956870727))/4204226981644629723547549,(4204226981644629723547549*%i+4*sqrt(10903152157933135742303986703335184117163956870727))/4204226981644629723547549],[1,1]]
(%i24) bfloat(realpart(res[1][1])); imagpart(res[1][1]);
(%o24) −3.1415926535897932384626433832795028841971693993751b0
(%o25) 1
(%i22) -bfloat(%pi);
(%o22) −3.1415926535897932384626433832795028841971693993751b0
(%i26) bfloat(realpart(res[1][2])); imagpart(res[1][2]);
(%o26) +3.1415926535897932384626433832795028841971693993751b0
(%o27) 1
---------- Mathmatica res ----------
In[2]:= N[Pi, 50]
Out[2]= 3.1415926535897932384626433832795028841971693993751

maxima対角化は時間が掛かりすぎるから現実的ではない・・・か
(%i1) f[i,j]:=random(1.0);
(%o1) f[i,j]:=random(1.0)
(%i6) showtime:true;
Evaluation took 0.0000 seconds (0.0000 elapsed) using 0 bytes.
(%o6) true
(%i12) mat:genmatrix(f,4,4);
Evaluation took 0.0000 seconds (0.0000 elapsed) using 0 bytes.
(%o12) matrix([.9138095996128959,.1951846177977887,.4823905248516196,.3554400571592862],[.4155627227214829,.6125738892538246,.8418126553359213,.5055894446438118],[0.0363933158999219,.6252975183418072,.9937314578694818,.04158924615444981],[.5931521815169765,.6546577661145074,.4707106740336644,.05192355366451751])
(%i13) eigenvalues(mat);
Evaluation took 9.2310 seconds (14.0640 elapsed) using 2764.567 MB.

2013年11月16日土曜日

Mathematicaの糖衣構文

Mathematicaの関数へのアクセスは糖衣構文によって色々な人に受け入れやすい事になっています.
多分,同じ数学概念でも棲んでいる世界によって解釈が異なったりするからだと思います.ただ,私には甘すぎては苦手です.

可読性って概念はMathematicaにはないのかね…
他人のプログラム読む時,絶望を覚えるんですが・・・

例えば以下の例は全てSqrt[3]を返します.

In[1]:= f[x_, y_: 2] := x^(1/y)

In[2]:= f[3]
f[3, 2]
3~f~2
3 // f
3 // f[#, 2] &
f @ 3
f @@ {3}
f @@ {3, 2}
f[#, 2] &  @ 3
Apply[f, {3}]
Apply[f, {3, 2}]
Apply[f, {##}] &  @@ {3, 2}


Out[2]= Sqrt[3]

Out[3]= Sqrt[3]

Out[4]= Sqrt[3]

Out[5]= Sqrt[3]

Out[6]= Sqrt[3]

Out[7]= Sqrt[3]

Out[8]= Sqrt[3]

Out[9]= Sqrt[3]

Out[10]= Sqrt[3]

Out[11]= Sqrt[3]

Out[12]= Sqrt[3]

Out[13]= Sqrt[3]

最近Mathematicaに関する事ばかりだな・・・

Mathematica Fourier変換について

バグを疑い始めたので,プログラムを純化させてみる.
In[1]:= IdentityTr[x_] := InverseFourier[Fourier[x]]

In[2]:= a = {0, 0, 0.33333333333333333, 0}
b = {0, 0, 0.333333333333333333, 0}
Precision[a[[3]]]
Precision[b[[3]]]

Out[2]= {0, 0, 0.333333, 0}

Out[3]= {0, 0, 0.33333333333333333, 0}

Out[4]= MachinePrecision

Out[5]= 17.5229

In[8]:= NestList[IdentityTr, a, 20]
NestList[IdentityTr, b, 20]

Out[8]= {{0, 0, 0.333333, 0}, {0., 0., 0.333333, 0.}, {0., 0., 
  0.333333, 0.}, {0., 0., 0.333333, 0.}, {0., 0., 0.333333, 0.}, {0., 
  0., 0.333333, 0.}, {0., 0., 0.333333, 0.}, {0., 0., 0.333333, 
  0.}, {0., 0., 0.333333, 0.}, {0., 0., 0.333333, 0.}, {0., 0., 
  0.333333, 0.}, {0., 0., 0.333333, 0.}, {0., 0., 0.333333, 0.}, {0., 
  0., 0.333333, 0.}, {0., 0., 0.333333, 0.}, {0., 0., 0.333333, 
  0.}, {0., 0., 0.333333, 0.}, {0., 0., 0.333333, 0.}, {0., 0., 
  0.333333, 0.}, {0., 0., 0.333333, 0.}, {0., 0., 0.333333, 0.}}

Out[9]= {{0, 0, 0.33333333333333333, 0}, {0.*10^-18, 
  0.*10^-18 + 0.*10^-18 I, 0.33333333333333333, 
  0.*10^-18 + 0.*10^-18 I}, {0.*10^-17, 0.*10^-18 + 0.*10^-18 I, 
  0.3333333333333333, 0.*10^-18 + 0.*10^-18 I}, {0.*10^-17, 
  0.*10^-17 + 0.*10^-17 I, 0.333333333333333, 
  0.*10^-17 + 0.*10^-17 I}, {0.*10^-16, 0.*10^-16 + 0.*10^-16 I, 
  0.333333333333333, 0.*10^-16 + 0.*10^-16 I}, {0.*10^-15, 
  0.*10^-15 + 0.*10^-15 I, 0.33333333333333, 
  0.*10^-15 + 0.*10^-15 I}, {0.*10^-14, 0.*10^-15 + 0.*10^-15 I, 
  0.3333333333333, 0.*10^-15 + 0.*10^-15 I}, {0.*10^-14, 
  0.*10^-14 + 0.*10^-14 I, 0.333333333333, 
  0.*10^-14 + 0.*10^-14 I}, {0.*10^-13, 0.*10^-13 + 0.*10^-13 I, 
  0.33333333333, 0.*10^-13 + 0.*10^-13 I}, {0.*10^-12, 
  0.*10^-12 + 0.*10^-12 I, 0.33333333333, 
  0.*10^-12 + 0.*10^-12 I}, {0.*10^-11, 0.*10^-11 + 0.*10^-11 I, 
  0.3333333333, 0.*10^-11 + 0.*10^-11 I}, {0.*10^-10, 
  0.*10^-11 + 0.*10^-11 I, 0.333333333, 
  0.*10^-11 + 0.*10^-11 I}, {0.*10^-10, 0.*10^-10 + 0.*10^-10 I, 
  0.33333333, 0.*10^-10 + 0.*10^-10 I}, {0.*10^-9, 
  0.*10^-9 + 0.*10^-9 I, 0.33333333, 
  0.*10^-9 + 0.*10^-9 I}, {0.*10^-8, 0.*10^-8 + 0.*10^-8 I, 0.3333333,
   0.*10^-8 + 0.*10^-8 I}, {0.*10^-7, 0.*10^-8 + 0.*10^-8 I, 0.333333,
   0.*10^-8 + 0.*10^-8 I}, {0.*10^-7, 0.*10^-7 + 0.*10^-7 I, 0.33333, 
  0.*10^-7 + 0.*10^-7 I}, {0.*10^-6, 0.*10^-6 + 0.*10^-6 I, 0.3333, 
  0.*10^-6 + 0.*10^-6 I}, {0.*10^-5, 0.*10^-5 + 0.*10^-5 I, 0.3333, 
  0.*10^-5 + 0.*10^-5 I}, {0.*10^-4, 0.*10^-4 + 0.*10^-4 I, 0.333, 
  0.*10^-4 + 0.*10^-4 I}, {0.*10^-3, 0.*10^-4 + 0.*10^-4 I, 0.33, 
  0.*10^-4 + 0.*10^-4 I}}

2013年11月9日土曜日

Mathematicaでの任意制度計算

Mathematicaが手に入ったので,いろいろテストしてみました.
その結果, 任意精度計算を行うにあたって精度の保持はユーザーが監視しないと行けない事が分りました.
また,監視を怠ると正しい結果が得られない事も分りました.

簡単な例として離散Fourier変換(FFTか?)と逆変換を繰り返す演算をさせます.

Wolfram Mathematica Document Foureirの説明にある通り,
Inputのlistが厳密な数値である場合,
N関数を適用し精度はMachinePrecisionに落とされます.
そのためFoureir関数は厳密精度計算を行ってくれません.

しかしながらユーザーがInputのlistにN関数を適用し精度をユーザー側でコントロールしたい状況が発生しますが,安易な方法では正しい結果が得られません.

その結果が In[3]です.

In[1]:= func1[y_, iter_, dps_: 10] := Module[{i, x},
  x = y;
  For[i = 0, i < iter, i++;
   If[ Mod[i, 2] != 0, x = SetPrecision[Fourier[x], dps], 
    x = SetPrecision[InverseFourier[x], dps]];
   ];
  Return[x]]
func0[y_, iter_] := Module[{i, x},
  x = y;
  For[i = 0, i < iter, i++;
   If[ Mod[i, 2] != 0, x = Fourier[x], x = InverseFourier[x]];
   ];
  Return[x]]


In[3]:= a = {0, 0, 1/3, 0}
res = func0[a, 100]
Precision[res[[1]]]
$MachinePrecision

b = N[{0, 0, 1/3, 0}, 10]
res = func0[b, 100]
Precision[res[[1]]]


Out[3]= {0, 0, 1/3, 0}

Out[4]= {0., 0., 0.333333, 0.}

Out[5]= MachinePrecision

Out[6]= 15.9546

Out[7]= {0, 0, 0.3333333333, 0}

Out[8]= {0.*10^5, 0.*10^4 + 0.*10^4 I, 0.*10^5, 0.*10^4 + 0.*10^4 I}

Out[9]= 0.
listに厳密な数値を代入した場合,正しい結果がMachinPrecisionで帰ってきます.
一方listにN関数を作用させて作った初期条件の場合,全ての結果のPrecisionが0になっています.

これでは使い物になりません.1回のFourier関数の演算でどうやら10^2-10^10程度の桁落ちが発生し最終的にゼロが帰ってくる様です.

上記問題を回避する為にはSetPrecision関数を用いて,強制的に精度を維持する事です
(func1関数を参照).
ただしこの方法が正当な回避方法かどうか私には分りませんが,少なくとも正しい答えが望んだ精度で帰ってくる様です(In[10]以降).
In[10]:= a = {0, 0, 1/3, 0}
res = func1[a, 100, 10]
Precision[res[[3]]]
b = N[{0, 0, 1/3, 0}, 10]
res = func1[b, 100, 10]
Precision[res[[3]]]

Out[10]= {0, 0, 1/3, 0}

Out[11]= {0, 0.*10^-21 + 0.*10^-21 I, 0.3333333333, 
 0.*10^-21 + 0.*10^-21 I}

Out[12]= 10.

Out[13]= {0, 0, 0.3333333333, 0}

Out[14]= {0, 0.*10^-21 + 0.*10^-21 I, 0.3333333333, 
 0.*10^-21 + 0.*10^-21 I}

Out[15]= 10.

In[34]:= Exit[]

2013年10月27日日曜日

mpmath 速度比べ

import mpmath
import numpy as np
import time

sample=500
iter=500

x = np.array([1.0/3.0]*sample)
init=x[0]
start = time.time()
for i in range(iter):
    if i<iter/2:
        x = x + 1.0/3.0
    else:
        x = x - 1.0/3.0
goal=time.time()
print("---- numpy array -----")
print("%.20f" % init)
print("%.20f" % x[0])
print("time:", goal-start)

mpmath.mp.dps = 18
x = np.array([mpmath.mpf("1.0")/mpmath.mpf("3.0")]*sample)
init=x[0]
start = time.time()
for i in range(iter):
    if i<iter/2:
        x = x + mpmath.mpf("1.0")/mpmath.mpf("3.0")
    else:
        x = x - mpmath.mpf("1.0")/mpmath.mpf("3.0")
goal=time.time()
print("---- mpmath with numpy array (dps=%d) -----" % mpmath.mp.dps)
print("%.20s" % init)
print("%.20s" % x[0])
print("time:", goal-start)

mpmath.mp.dps = 100
x = np.array([mpmath.mpf("1.0")/mpmath.mpf("3.0")]*sample)
init=x[0]
start = time.time()
for i in range(iter):
    if i<iter/2:
        x = x + mpmath.mpf("1.0")/mpmath.mpf("3.0")
    else:
        x = x - mpmath.mpf("1.0")/mpmath.mpf("3.0")
goal=time.time()
print("---- mpmath with numpy array (dps=%d) -----" % mpmath.mp.dps)
print("%.20s" % init)
print("%.20s" % x[0])
print("time:", goal-start)

mpmath.mp.dps = 18
x = mpmath.matrix([mpmath.mpf("1.0")/mpmath.mpf("3.0")]*sample)
init=x[0]
start = time.time()
for i in range(iter):
    if i<iter/2:
        x = x + mpmath.mpf("1.0")/mpmath.mpf("3.0")
    else:
        x = x - mpmath.mpf("1.0")/mpmath.mpf("3.0")
goal=time.time()
print("---- mpmath matrix (dps=%d)-----" % mpmath.mp.dps)
print("%.20s" % init)
print("%.20s" % x[0])
print("time:", goal-start)

mpmath.mp.dps = 100
x = mpmath.matrix([mpmath.mpf("1.0")/mpmath.mpf("3.0")]*sample)
init=x[0]
start = time.time()
for i in range(iter):
    if i<iter/2:
        x = x + mpmath.mpf("1.0")/mpmath.mpf("3.0")
    else:
        x = x - mpmath.mpf("1.0")/mpmath.mpf("3.0")
goal=time.time()
print("---- mpmath matrix (dps=%d)-----" % mpmath.mp.dps)
print("%.20s" % init)
print("%.20s" % x[0])
print("time:", goal-start)


x = [mpmath.mpf("1.0")/mpmath.mpf("3.0")]*sample
init=x[0]
start = time.time()
for i in range(iter):
    if i<iter/2:
        x = [xx + mpmath.mpf("1.0")/mpmath.mpf("3.0") for xx in x]
    else:
        x = [xx - mpmath.mpf("1.0")/mpmath.mpf("3.0") for xx in x]
goal=time.time()
print("---- list -----")
print("%.20f" % init)
print("%.20f" % x[0])
print("time:", goal-start)




実行結果
$ python speed_test_array.py 
---- numpy array -----
0.33333333333333331483
0.33333333333332632042
('time:', 0.0022330284118652344)
---- mpmath with numpy array (dps=18) -----
0.333333333333333333
0.333333333333333326
('time:', 1.7795300483703613)
---- mpmath with numpy array (dps=100) -----
0.333333333333333333
0.333333333333333333
('time:', 2.18923282623291)
---- mpmath matrix (dps=18)-----
0.333333333333333333
0.333333333333333326
('time:', 5.815176010131836)
---- mpmath matrix (dps=100)-----
0.333333333333333333
0.333333333333333333
('time:', 6.364114046096802)
---- list -----
0.33333333333333331483
0.33333333333333331483
('time:', 12.696530103683472)

2013年10月15日火曜日

gccの4倍精度ライブラリー

gcc 4.6 から4倍精度計算ができる様になりました.昔は4倍精度で計算できるマシンがあったみたいだけどね.

今回は、c, (c++はやってない), fortranの言語で実装するとしてどの言語が最も実装しやすいか比較する.
gcc のversionは4.8 を用いる.結論としてはfortranが一番楽だと思う.
cおよびc++は
http://homepage3.nifty.com/maikaze/float128.html
で解説されているようにc(, c++)では幾分実装が不自然な所もあるので,一応比較してみる.
計算する内容は1/3, 4*arctan(1)=pi,sqrt(2)である.

まずfortranから
program fort_quadmath
 implicit none
 real (kind(0D0)) dx, dy, dpi,dsqrt2
 real (kind(0Q0)) qx, qy, qpi,qsqrt2
 dx = 1.0D0;
 dy = 3.0D0;
 dsqrt2 = sqrt(2.0D0)
 dpi = 4.0D0*atan(1.0D0)
 qx = 1.0Q0;
 qy = 3.0Q0;
 qsqrt2 = sqrt(2.0Q0)
 qpi = 4.0Q0*atan(1.0Q0)
 print *, "---- double precision ----"
 print *, "1/3     =", dx/dy
 print *, "pi      =", dpi
 print *, "sqrt(2) =", dsqrt2
 print *, "----  quad  precision ----" 
 print *, "1/3     =", qx/qy
 print *, "pi      =", qpi
 print *, "sqrt(2) =", qsqrt2
end program fort_quadmath 
コンパイルは
 gfortran -o fort_quadmath fort_quadmath.f03 
であり,実行結果は
$ ./fort_quadmath 
 ---- double precision ----
 1/3     =  0.33333333333333331     
 pi      =   3.1415926535897931     
 sqrt(2) =   1.4142135623730951     
 ----  quad  precision ----
 1/3     =  0.333333333333333333333333333333333317      
 pi      =   3.14159265358979323846264338327950280      
 sqrt(2) =   1.41421356237309504880168872420969798   
 pi           3.1415926535897932384626433832795028841971693993751
 sqrt(2)      1.4142135623730950488016887242096980785696718753769
厳密解と比較して32桁くらいまでは遜色無く一致している.
実装に関して変数の型の宣言と、リテラル(0Q0)さえ気をつけていれば,関数の型を意識する必要がない(というか関数が代表名で与えられているので,変数の型によって関数が勝手に陰的に変わるんだと理解している.)

一方c言語に関して
#include <stdio.h>
#include <quadmath.h>
#include <math.h>

int main(void)
{
 char c1[80],c2[80],c3[80],c4[80],c5[80];
 double dx,dy,dpi;
 long double lx,ly,lpi;
 __float128 x,y, pi1,pi2,sqrt1,sqrt2;
 dx = 1.0;
 dy = 3.0;
 dpi = 4.0*atan(1.0);
 lx = 1.0L;
 ly = 3.0L;
 lpi = 4.0L*atan(1.0L);
 x = 1.0Q;
 y = 3.0Q;
 pi1 = 4.0Q*atan(1.0Q);
 pi2 = 4.0Q*atanq(1.0Q); 
 sqrt1 = sqrt(2.0Q);
 sqrt2 = sqrtq(2.0Q);  
 
 printf("---- double precision ----\n");
 printf("1/3 = %.30lf\n", dx/dy);
 printf("pi  = %lf\n", dpi); 
 printf("---- long double prescision ----\n");
 printf("1/3 = %.30Lf\n", lx/ly);
 printf("pi  = %.30Lf\n", lpi);
 printf("---- quad precision ----\n");
 quadmath_snprintf(c1,80,"%.40Qf", x/y);
 quadmath_snprintf(c2,80,"%.40Qf", pi1); 
 quadmath_snprintf(c3,80,"%.40Qf", pi2);  
 quadmath_snprintf(c4,80,"%.40Qf", sqrt1);  
 quadmath_snprintf(c5,80,"%.40Qf", sqrt2);    
 printf("1/3 = %s\n", c1);
 printf("pi1  = %s\n", c2);
 printf("pi2  = %s\n",c3);
 printf("sqrt1  = %s\n",c4);
 printf("sqrt2  = %s\n",c5);  
 return 0;
}
のコンパイルは
gcc -o c_quadmath c_quadmath.c -lquadmath
で行う.実行結果は
$ ./c_quadmath 
---- double precision ----
1/3 = 0.333333333333333314829616256247
pi  = 3.141593
---- long double prescision ----
1/3 = 0.333333333333333333342368351437
pi  = 3.141592653589793115997963468544
---- quad precision ----
1/3 = 0.3333333333333333333333333333333333172839
pi1  = 3.1415926535897931159979634685441851615906
pi2  = 3.1415926535897932384626433832795027974791
sqrt1  = 1.4142135623730951454746218587388284504414
sqrt2  = 1.4142135623730950488016887242096981769402
 
pi     3.1415926535897932384626433832795028841971693993751
sqrt(2)  1.4142135623730950488016887242096980785696718753769
の厳密解と比較して pi1とsqrt1がdouble程度の精度しかない. これは代入する変数の型と関数の返り値の型が一致していないため,関数の方の型で丸められてしまったせいだと考えられる.(c++でも?)
既存のプログラムの精度を倍精度から4倍精度に書き直すにはc言語では関数も書き換えなければならない. (gccがversion upしたらなんとかなるんか?) 精度の仕様変更に関しては圧倒的にfortranの方がいじりやすいと個人的に思う.
 (fortran はあんまり使ってないけどね・・・)