１．はじめに
BASIC を使ってモタモタと数値計算していた Windows3.1 の時代に、突然、MATLAB の凄さを見せられたときの感動は今でも忘れません。それ以来、MATLAB にはずっとお世話になってきていますが、忘れられないのは感動だけではありません。その当時に独り善がりで会得した「利用上の心得」も記憶に残り、未だにその影響から抜け出せずにいます。その心得の一つが「実用上の行列サイズの上限は、精々250行×250列程度」（風評被害防止のため淡色表示）というものです。誰に感化された訳でもなく、この程度のサイズの行列を扱う中で、なかなか計算結果が出てこないことに業を煮やして、身体で感じ取ったノウハウです。その後、あまり大きな行列を扱う機会は訪れず、この先入観を是正することができないまま現在に至ってしまいました。
ところが、現行のバージョン（と言ってもチョット古い R2019a）で、コマンドラインから下記を実行させてみたところ、3000行×3000列の大行列の演算を2秒も掛からないうちに完了し、$|A^{-1}||A|{=}1$ を実証する高精度で妥当な答を返してきました。肝心な計算は手抜きして、公式集を参考にして答えただけではないか？と疑うほどの速さです。しかし、微妙な誤差を含んでいることから見ると、一つ一つの要素を大切にし、几帳面に計算した結果であることは明らかです。いつからこんなにも処理能力が向上したのでしょうか？
>> tic; A=rand(3000)/10; [det(inv(A))*det(A) toc]
ans =
1.000000000002983e+00 1.415134600000000e+00
・・・ということで、新旧バージョンの処理能力（実行時間と処理可能サイズなど）を比較してみようと思い立ち、色々とやってみました。
２．評価方法
２.１ 使用するMATLABのバージョン
次の４種類を使用しました。Windows3.1 の時代の MATLAB も仲間に入れたいところですが、当時の勤務先の会社が購入したものなので手元にはありません。インストール用のメディアは何十枚もの3.5インチのフロッピーディスクでしたから、残っていたとしても、昔の環境の再現にはかなり苦労したことでしょう。
- R2019a Home
-
R2007b Education
以上を走らせるパソコン(FMVU75B3B)は
Windows 10 Home, 64bit, Core i5-8250U, 1.6GHz, RAM 4GB, SSD
LCD 13.3型ワイド 1920×1080 拡大率150%で使用（1280×720相当）, 1.0kg
- ver7.1(R14SP3) Education
-
ver5.3(R11) Student
以上を走らせるパソコン(FMVNB70ET)は
Windows XP Home Edition, 32bit, pentium4-M, 2.2GHz, RAM 768MB, HDD
LCD 15型XGA 1024×768, 4.0kg(含:ACアダプタ、0.2kg分解能体重計による)
上記の Windows10 パソコンは、性能よりも可搬性を重視して購入したものです。スタバで格好良く使って、若者に負けていないことを誇示したいという邪心からです。MATLAB からは、「このレベルのパソコンによるものを、いかにも性能限界であるかのように吹聴されては迷惑だ！」とお叱りを受けそうです。しかし、一泡沫ユーザーによる私的で我流な評価に過ぎませんので何卒ご容赦ください。
後者の XP パソコンは当家の４世代前のビンテージものです。埃を被っていたものを引っ張り出してきました。液晶画面は一部の表示が欠け、キーボードの文字も消えかかっていましたが、評価に必要な機能は何とか確保できました。内蔵タイマーがリセットされたままで MATLAB が起動できなかったり、評価中、HDD から壊れそうな音が聞こえて不安に陥ったりと、予想しなかった苦労も色々とありましたが・・・。
上述のように、MATLAB の処理能力はバージョンだけで決まる訳ではなく、パソコンの仕様にも大きく左右されます。その観点から、MATLAB の各バージョンの発表時期を調べてみると、当時、最も可能性があったパソコンとの組み合わせは、次のようなものと思われます。
| MATLABの バージョン |
同時代のパソコン |
|---|---|
| R2019a(2019) | Windows10(2015/7) |
| R2007b(2007) | WindowsVista(2006/11) or XP(2001/9) |
| ver7.1(2005) | WindowsXP(2001/9) |
| ver5.3(1999) | Windows98(1998/7) or 95(1995/11) |
| ？？？ | Windows3.1(1993/5) |
これによると、今回使用した MATLAB のうちの R2007b と ver5.3 については、組合せているパソコンが当時のものよりも新型になっています。これでは、当時実感できた処理能力よりも良好な評価結果が出てくる可能性があります。したがって、これらの結果については、少々割り引いて評価するのが良さそうです。しかし、同時代のパソコンであっても、価格に応じて性能はピンからキリまで色々です。細かなことを言い始めると収拾がつかなくなりそうなのでこの程度でやめておきます。
２.２ 評価対象にする演算処理
計算規模が大きく時間がかかりそうな処理として、各種の行列計算と meshgrid 応用計算について調べることにします。
２.２.１ 行列計算
演算対象の行列の生成
行列計算としては、行列式 (det)・掛け算 (*, mtimes)・逆行列 (inv)・固有値 (eig) の演算能力について調査します。具体的には、行列サイズを徐々に大きくしていって、これらの処理に要する時間を測定することになります。
そのためには、まず、演算対象となる巨大な行列を生成しておく必要があります。これには乱数を使うのが最も手軽ですが、一つ解決しておかなければならない問題があります。「はじめに」で示したコマンドラインからの入力では、A=rand(3000) で、0～1 の範囲の乱数を要素にした 3000×3000 のサイズの行列を生成しています。しかし、それだけでは行列式がオーバーフローして ±Inf になってしまうため、全要素の値を 1/10 に減じてこの問題を回避しています。
ところが、この除数の値は常に 10 にしておけば良い訳ではありせん。かなりの自由度はありますが、行列のサイズに応じて正しく調整しないと、答がオーバーフローやアンダーフローして評価どころではなくなります。そこで、まずは、最適な除数を決めるための実験をしました。「最適」の定義は自分で勝手に決めました。行列式の値が、オーバーフローやアンダーフローから最も離れた状態、すなわち、ほぼ e+00 オーダーになることであるとします。
行列式の値を何故そこまで重要視するのか疑問を持たれるかもしれません。その理由は次のとおりです。この実験では、せっかく処理時間を計測しても、その演算結果が誤っているのでは意味がありません。まずは、演算結果が正しいことを確認しなければなりません。しかし、演算結果の巨大サイズの行列要素を一つずつ確認するのは無理です。そこで、行列式が持つ性質を利用して、巨視的に正否の判定をします。したがって、行列式の値が全ての判断基準になります。これが正しく表現できなければ演算結果の評価もできません。
コマンドラインから下記を入力し、行列サイズ $n$ の値を設定したうえ、除数 $K$ の値を微妙に変えながら、出力された $|A|$ の値がほぼ e+00 オーダーになる点を捜します。
K=10; n=3000; disp(datestr(now,'dd-mmm-yyyy HH:MM:SS')), pause(0.5); A=rand(n)/K; disp(det(A))
この結果から得られたのが図１の青〇のプロットです。このプロットが、行列サイズに応じた最適な除数の関係を示している訳ですが、なんと、$\pmb{K{=}0.175\sqrt{n}}$ の赤色の線上に綺麗に載っていることが分かります。そこで、当実験では、次の式で操作対象の行列を生成することにします。
A=rand(n)/(0.175*sqrt(n));
このようなシンプルな関係式であれば、解析的に簡単に導けるように思えますが、残念ながらまだ達成できておりません。なお、図１には、最適な $K$（Kbest） だけでなく、アンダーフロー直前の Kmax と オーバーフロー直前の Kmin の曲線も併記しています。Kmin ～ Kmax で囲まれた範囲が、正常な $|A|$ を出力できる領域になります。かなりの自由度がありますが、n が非常に大きくなると、K の値の許容範囲がどんどん狭くなっていきます。
ver5.3 の Kmax の曲線は、最適K線を漸近線にして素直に伸び上がっています。しかし、これよりも新しい他のバージョン（図では R2019a で代表している）のものは、なにやら縮こまった不自然な形になっているのが気になります。両者では行列式の計算方法や乱数生成のアルゴリズムが異なるのでしょうか？。これが原因で、最適K線を基準にしている当実験では、演算可能な行列のサイズは、新しいバージョンの方が低くなる傾向にあります。
ところで、これらのデータの採取には途方もない時間がかかって閉口しました。いくら暇だらけの老人とはいえ、このような作業はもうコリゴリです。
演算結果の正否の判断
下表の左欄の各演算結果について、右欄の式を満足できていることで正答であることを確認します。正答でなかったケースについては、処理時間を計測したところで意味がないので計測対象から外します。
| 判定対象 | 正否判定基準 |
|---|---|
| 行列式 $|A|$ | $|A^{\mathrm{T}}|/|A|{=}1$ |
| 行列積 $AA$ | $|AA|/|A|^2{=}1$ |
| 逆行列 $A^{-1}$ | $|A^{-1}||A|{=}1$ |
| 固有値 $\lambda_i$ | $(\prod_{i=1}^{n}\lambda_i)/|A|{=}1$ |
誤答のときには、正否判定値が 1 から懸け離れた値（ ±Inf, 0, NaN, 4.5e-123, etc ）になるので、判断に迷うようなことはありません。多くは 1.000000 のレベルで 1 を満たしますが、一部、1.000099 となるものもありました。これでも正答扱いにしています（ ver7.1 の n=4000 における行列積）。
処理可能な行列サイズの見積もり
特に、古いバージョンでは、一連のデータを採取するのに途方もなく長い時間がかかります。データ採取プログラムを適当に走らせるだけでは、途中でフリーズして強制終了にでもなると最初からやり直しになってしまいます。・・・ということで、各バージョンごとに、処理可能な行列サイズの範囲を事前に見定めておく必要があります。そのために、コマンドラインで直接操作する下記のコマンド文字列を使いました。データ採取プログラムの中で指示している「処理すべき行列サイズ」の範囲は、これらの結果から決めています。
% 性能限界の事前確認（行列式計算）
n=1000; disp(datestr(now,'dd-mmm-yyyy HH:MM:SS')), pause(0.5); A=rand(n)/(0.175*sqrt(n)); detA=det(A); B=A'; tic; ansA=det(B); ta=toc; jg=ansA/detA; tb=toc; disp([ta jg tb])
% 性能限界の事前確認（行列積計算）
n=1000; disp(datestr(now,'dd-mmm-yyyy HH:MM:SS')), pause(0.5); A=rand(n)/(0.175*sqrt(n)); detA=det(A); tic; ansA=A*A; ta=toc; jg=det(ansA)/detA^2; tb=toc; disp([ta jg tb])
% 性能限界の事前確認（逆行列計算）
n=1000; disp(datestr(now,'dd-mmm-yyyy HH:MM:SS')), pause(0.5); A=rand(n)/(0.175*sqrt(n)); detA=det(A); tic; ansA=inv(A); ta=toc; jg=det(ansA)*detA; tb=toc; disp([ta jg tb])
% 性能限界の事前確認（固有値計算）
n=1000; disp(datestr(now,'dd-mmm-yyyy HH:MM:SS')), pause(0.5); A=rand(n)/(0.175*sqrt(n)); detA=det(A); tic; ansA=eig(A); ta=toc; jg=abs(prod(ansA)/detA); tb=toc; disp([ta jg tb])
評価用のプログラム
以上を踏まえ、末尾の「プログラム」欄に添付しているデータ採取用プログラムを作りました。使用するどの MATLAB のバージョンでも動作できるように調整してあります。古い MATLAB では受け付けてくれないコマンドも多いので、昔を思い出しながら共用化するのには結構苦労しました。ネット上では旧バージョンの情報はなかなか得られないので、コマンドウィンドウ上だけで、そのバージョンについての使用方法が分かる help コマンドは大いに役立ちます。
古いバージョンで特に不便を感じたのは、R2007b を含め、これ以前のバージョンでは、スクリプト中でローカル関数を定義できない点です。処理能力の向上を評価するためには、処理できるデータ規模や処理時間以外にも、このような便利機能を考慮に入れることも大切だと感じました。
さて、計算の都度、実行時間は微妙に変化します。そのバラつきの傾向も把握するために、全く同一の計算を５回にわたって実行させます。これらから、平均時間と最短・最長時間を求めてグラフ化することにします。
なお、データ採取用プログラムを走らせているときには、外乱を避けるために他のプログラムは並走させないように注意しました。しかし、便宜上、待機状態のエクスプローラと秀丸エディタのウィンドウは開いたままでも良いことにしました。また、ver5.3 では、プログラム実行中にコマンドウィンドウ上をクリックすると「反応なし」となり、タスクマネージャーによる強制終了以外の操作は受付けなくなります。結果が出るのが待ち遠しくてついついやってしまいがちですが、プログラムが終了するまでひたすら我慢です。
２.２.２ meshgrid 応用計算
等高線の描画
処理能力の評価のために、contour コマンドを使って図２のような等高線を描かせます。等高線はメッシュを細かくすればするほど滑らかになりますが、一般には、一辺を 100 ～ 200 分割程度にすれば十分な滑らかさが得られます。それ以上にしても見分けはつきません。5000 分割まで評価して何のメリットがあるのか？という疑問は残りますが、処理能力の把握のための試練だと割り切って強行しました。
描かせる等高線の本数は処理時間に直接影響してくるはずですが、今回、その関係について深くは考えていません。特別な理由はなく、適度な本数として８レベルを選んでいます。
評価用のプログラム
基本的な処理内容は、コマンドラインからでも操作できる次のようなものです。
n=1000; x=linspace(-1,1,n+1); [X,Y]=meshgrid(x,x); Z=X.^2.*(1-X.^2)-Y.^2; tic; [c,h]=contour(X,Y,Z,[-0.5:0.1:0.2]); toc
上記の [c,h]=contour(...) の部分を単に contour(...) としても等高線は描けますが、toc から得られる処理時間があまりにも短いように感じました。どうも、contour だけでは返り値が存在しないので、メインのルーチンが contour に引数を渡しただけの段階で処理完了を宣言してしまっているのではないかと疑念を持ちました。返り値の取得を強制させるように [c,h]= を追加してみたところ、処理時間は妥当と思える長さまで延びたので、こちらの方法を採用しました。
評価用プログラムでは、上記と同様の操作を meshgrid の分割数を大きくしていきながら繰り返します。行列演算の評価と同じく、全く同一の計算を５回にわたって実行させ、処理時間の平均値と最短・最長値を求めます。
３．評価結果
３.１ 演算能力
データ採取用プログラムから得られた情報を整理して、次の一連のグラフができました。図３の上は行列式、下は行列積、図４の上は逆行列、下は固有値について、それぞれ、行列サイズと処理時間の関係を示したものです。また、図５は等高線描画について、メッシュ分割数と描画時間の関係を示したものです。
全体的には、MATLAB のバージョンが新しくなるほど、処理時間が短くなっているのが良く分かります。しかも、対数目盛でなければ表現できないほどの大きな進化です。ものによっては３デカード（1000 倍）近くの高速化が達成されています。Windows3.1 時代の私の MATLAB への評価が妥当なものであったかどうかは不明ですが、漠然としていた新旧の性能差を明確にすることができました。
しかし、処理できる行列のサイズは、最も古い ver5.3 の方が大きかったり、等高線の処理時間は一つ古いバージョンの方が短かったりするのは意外でした。また、行列のサイズ $n$ と処理時間 $t$ の関係はほぼ $t=an^3$ 、領域の分割数 $n$ と等高線描画時間 $t$ の関係はほぼ $t=an^2+b$ で表せそうなことも分かりました。
同一条件で５回繰り返した結果を示していますが、薄い色で塗り潰した領域はそれほど広くはなく、殆どのケースでバラツキは意外に少ないように見えます。しかし、時期を改めて同一のデータを採ってみると、プロットの形はかなり変化するのが実状です。処理速度の順番が入れ替わるほどのことはありませんが、ここに示したものが不動のデータであるとの認識は禁物です。
３.２ 可視化能力
処理能力の一部として扱うのが適切なのかどうかは分かりませんが、一連の作業の中で気づいたのが、演算結果を可視化する figure の表現力のバージョンによる違いです。
図６～１１は MATLAB の各バージョンの figure をスクリーンショットし、そのままの画素数で png 化したものです。なお、R2019a と R2007b については、ディスプレイの拡大率を指定できるパソコンで動作させているので、それぞれを、拡大率 100% と 150% の２つの条件で描かせています。また、R2019a では図２のものが規定色ですが、旧バージョンと条件を揃えるように設定変更しています。
gcf から得られる Position 情報によれば、どの figure も幅 560px 高さ 420px でバージョン間で差はありません。また gca の Position 情報でも figure 内の axis の配置はすべて共通です。要するに、figure が作図に利用できる画素数は、バージョン間で差がありません。しかし、画像の表現力には差が見られます。
旧バージョンによる図８～１１のグラフはどれも似たり寄ったりで、ビットマップ感が強くて美しさに欠けます。しかし、R2019a による図６,７では、線や文字が非常に緻密に描かれており、一気に表現力が上がっています。中でも、拡大率 150% で描かせた図６が秀逸です。実質的に $1.5^2$ 倍に増えた画素数を作図に有効に利用しているようです。なお、R2007b には、拡大率を表現力の改善に活かすような機能はないようで、100% と 150% の図に大きな違いは見られません。
さて、ここで認識しておいて欲しいのは、旧バージョンでも、print コマンドで紙に印刷したり、各種画像のフォーマットに変換したりすれば、画面上の図よりも解像度の高い美しい図が得られるということです。画面上の表現力だけですべてを評価しようとするのは適切ではありません。しかし、印刷した図が、画面上からイメージしたものと微妙に異なって納得できず、その修正のためにプログラムの手直しを繰り返すこともありました。やはり、画面上の図と印刷結果が WYSIWYG の関係にあるのが理想です。
その点で、R2019a の figure の表現力の向上は大いに評価すべきものです。最近では、最終成果物として、印刷したものよりもスクリーンコピーしたものをそのまま利用する機会も多くなっています。私の Qiita への投稿も、MATLAB 関係の図はすべてこの方式によっています。画面上に表示されたものをそのまま最終成果物として利用できるという今の状態が、私にとっては究極の表現形態であると思っています。
４．おわりに
以上で新旧 MATLAB の能力の比較は終わりましたが、私の環境では確認できない次の性能についても興味が湧いてきました。面白い情報をお持ちの方が居られればお教えください。
- 最新の R2023a と高性能パソコンによる処理能力
- スーパーコンビュータ上で動作している MATLAB の処理能力
最後に、この記事を書いている過程で出くわした不可解現象について紹介しておきます。これは、どのように解釈すれば辻褄が合うのでしょうか？
既述の行列積の評価では、A*A なる演算を対象にして 4000×4000 のサイズの行列まで正答が得られたことになっています。しかし、当初は、これを一捻りしたつもりで A*A' としていたのですが、正答できる行列のサイズが 1500×1500 までしか上がりませんでした。ところが、 A*A'でも、ver5.3 だけはこれ以上のサイズでも大丈夫でした。また、A*flipud(A)、A*fliplr(A) とした場合には、どのバージョンでも大丈夫でした。確かに、A*A' では、同じ行列の行同士で内積をとる形になって特殊条件にはなりますが、乱数が一様であれば、条件に関係なく正答が得られても良いはずです。図１の Kmax のグラフにおいて、 ver5.3 以外は素直な形ではなかった点と何か関係はあるのでしょうか？
注：
|A*flipud(A)|=|A|^2、|A*fliplr(A)|=|A|^2 は、行列 A のサイズ n×n が n=4N or n=4N+1 (N=1,2,3,...)の場合に限り成立する関係です。この記事では n を 100 の倍数（ 4 の倍数でもある）にしているので、たまたま成立しているだけです。
５．プログラム
５.１ データ採取用プログラム（行列計算）
% testMatrix03.m
% MATLAB の行列計算能力の測定
% 「データ採取用プログラム」
% 評価対象関数： det, *(mtimes), inv, eig
clear
close all
testrep=5; % テストの繰返し回数
thisfile='testMatrix03.m';
% ■■■ 要選択 ■■■ ==========================
% テストするMATLABのバージョンを番号で選択
mat_ver=2; % 1:R2019a, 2:R2007b, 3:MATLAB7.1, 4:MATLAB5.3
% ===============================================
switch mat_ver
case 1
mat_name='R2019a';
case 2
mat_name='R2007b';
case 3
mat_name='ver7.1';
case 4
mat_name='ver5.3';
end
% 評価する行列サイズの一覧
% R2019a用
nm11=[1 2 4 6 8 [10:5:40]]*100; % det用 n=4500では誤答
nm12=[1 2 4 6 8 [10:5:40]]*100; % 積*用 n=4500では誤答
nm13=[1 2 4 6 8 [10:5:40]]*100; % inv用 n=4500では誤答
nm14=[1 2 4 6 8 [10:5:35]]*100; % eig用 n=4000では誤答
% R2007b用
nm21=[1 2 4 6 8 [10:5:40]]*100; % det用 n=4500では誤答
nm22=[1 2 4 6 8 [10:5:40]]*100; % 積*用 n=4500では誤答
nm23=[1 2 4 6 8 [10:5:40]]*100; % inv用 n=4500では誤答
nm24=[1 2 4 6 8 [10:5:35]]*100; % eig用 n=4000では誤答
% MATLAB7.1用
nm31=[1 2 4 6 8 [10:5:40]]*100; % det用 n=4500では誤答
nm32=[1 2 4 6 8 [10:5:40]]*100; % 積*用 n=4500では誤答
nm33=[1 2 4 6 8 [10:5:40]]*100; % inv用 n=4500では誤答
nm34=[1 2 4 6 8 [10:5:35]]*100; % eig用 n=4000では誤答
% MATLAB5.3用（ n>5000 は調査せず。未知）
nm41=[1 2 4 6 8 [10:5:50]]*100; % det用 n=5000でもOK。
nm42=[1 2 4 6 8 [10:5:45]]*100; % 積*用 n=5000で
% 「メモリが足りません」
nm43=[1 2 4 6 8 [10:5:50]]*100; % inv用 n=5000でもOK。
nm44=[1 2 4 6 8 [10:5:30]]*100; % eig用 n=3500でも正答
% はするものの、
% 1.5時間もかかる。
% n=3000で打ち切る。
nM={ nm11 nm12 nm13 nm14 ; ...
nm21 nm22 nm23 nm24 ; ...
nm31 nm32 nm33 nm34 ; ...
nm41 nm42 nm43 nm44 };
% データ記録用の変数
size_vs_time=[];
% 途中でフリーズしてもデータが残るように、ファイルにも逐次記録。
% そのためのファイルの準備。
Date=datestr(now,'dd-mmm-yyyy HH:MM:SS');
Date=strrep(strrep(strrep(Date,'-','_'),':','_'),' ','_');
fid=fopen(['xx_size_vs_time_' Date '.txt'],'w+');
fprintf(fid,['xx_size_vs_time_' Date '.txt']);
fprintf(fid,'\n');
fprintf(fid,'Ver Com Size Time_ave Time_min Time_max');
fprintf(fid,'\n');
fclose(fid);
com=[1 2 3 4]; % テストするコマンドの種類一覧
% 1:det, 2:mtimes, 3:inv, 4:eig
for nc=com; % 各種のコマンドについて繰り返す。
switch nc
case 1
com_name='行列式'; % テストする機能名
nmat=nM{mat_ver,1}; % 計算すべき行列サイズの一覧
command='B=transpose(A); tic; ansA=det(B); te=toc;';
% 実行時間計測用コマンド
% R2007bでは、文字列""はダメで''のみ。
% ゆえに文字列中で転置を'とはできない。
% それゆえtransposeを使う。
judge='Jg=ansA/detA;'; % 正否判定値（1.0で正）
case 2
com_name='行列積';
nmat=nM{mat_ver,2};
command='tic; ansA=A*A; te=toc;';
judge='Jg=det(ansA)/(detA^2);';
case 3
com_name='逆行列';
nmat=nM{mat_ver,3};
command='tic; ansA=inv(A); te=toc;';
judge='Jg=det(ansA)*detA;';
case 4
com_name='固有値';
nmat=nM{mat_ver,4};
command='tic; ansA=eig(A); te=toc;';
judge='Jg=abs(prod(ansA)/detA);';
end
for nm=nmat; % 各種の行列サイズについて繰り返す。
A=rand(nm)/(0.175*sqrt(nm)); % 操作対象の行列の生成。
detA=det(A); % 計算が暴走していないかの検証用の基準
timedata=[]; % 毎回の演算実行時間を記録する行ベクトル
for ne=1:testrep % バラツキの平均化のため繰り返し処理
eval(command); % det, *, inv, eig, tic, toc などを含む
% コマンド文字列の実行
if ne==1; % 1回目のみ実行。2回目以降は同値のはず
% なので、時間の節約のためバイパス。
eval(judge); % 正否判定値を生成する
% コマンド文字列の実行（Jgを生成）
end
% 正否判定値がほぼ1.000で正しいことと、実行時間の傾向を
% 目視監視するため、コマンドラインに毎回印字する。
disp([mat_name ' ' com_name ' ' ...
sprintf('%4s',num2str(nm)) '×' ...
sprintf('%4s',num2str(nm)) ...
' 正否判定 ' num2str(Jg,'%8.6f') ...
' 処理時間 ' num2str(te) ' 秒'])
timedata=[timedata te]; % 毎回の実行時間をすべて記録
end
% 繰返し実行の結果から、平均・最短・最長時間を計算・表示
t_ave=sum(timedata)/testrep; % 平均時間
t_min=min(timedata); % 最短時間
t_max=max(timedata); % 最長時間
disp(' ')
disp([' 平均時間：' num2str(t_ave,'%11.5f') ...
' 最短時間：' num2str(t_min,'%11.5f') ...
' 最長時間：' num2str(t_max,'%11.5f')]);
disp(' ')
% データ記録用の変数に計測結果を追記。
% （MATLABバージョン・コマンド・行列サイズ・諸時間）
size_vs_time=[size_vs_time; ...
mat_ver nc nm t_ave t_min t_max];
% 途中でフリーズしてもデータが残るように、ファイルにも逐次記録。
% ver5.3 では、曖昧に'%d'では不可。桁数の明示をしないとエラー。
fid=fopen(['xx_size_vs_time_' Date '.txt'],'a');
fprintf(fid,'%3s',num2str(mat_ver,'%1d'));
fprintf(fid,'%3s',num2str(nc,'%1d'));
fprintf(fid,'%6s',num2str(nm,'%4d'));
fprintf(fid,'%12s',num2str(t_ave,'%11.5f'));
fprintf(fid,'%12s',num2str(t_min,'%11.5f'));
fprintf(fid,'%12s',num2str(t_max,'%11.5f'));
fprintf(fid,'\n');
fclose(fid);
end
end
% 計測結果をコマンドラインに印字
disp(datestr(now,'dd-mmm-yyyy HH:MM:SS'));
disp( [' created by ' thisfile ' on ' mat_name]);
disp('Ver Com Size Time_ave Time_min Time_max');
B=size_vs_time;
for n=1:size(B,1)
disp([sprintf('%3s',num2str(B(n,1),'%1d')) ...
sprintf('%3s',num2str(B(n,2),'%1d')) ...
sprintf('%6s',num2str(B(n,3),'%4d')) ...
sprintf('%12s',num2str(B(n,4),'%11.5f')) ...
sprintf('%12s',num2str(B(n,5),'%11.5f')) ...
sprintf('%12s',num2str(B(n,6),'%11.5f'))])
end
５.２ データ採取用プログラム（等高線描画）
% testMatrix06.m
% MATLAB の等高線描画能力の測定
% 「データ採取用プログラム」
clear
close all
testrep=5; % テストの繰返し回数
thisfile='testMatrix06.m';
% ■■■ 要選択 ■■■ ==========================
% テストするMATLABのバージョンを番号で選択
mat_ver=1; % 1:R2019a, 2:R2007b, 3:MATLAB7.1, 4:MATLAB5.3
% ===============================================
switch mat_ver
case 1
mat_name='R2019a';
case 2
mat_name='R2007b';
case 3
mat_name='ver7.1';
case 4
mat_name='ver5.3';
end
% 評価する分割数の一覧
% R2019a用
nd1=[1 2 4 6 8 [10:5:50]]*100; % n=5500 は未確認
% R2007b用
nd2=[1 2 4 6 8 [10:5:50]]*100; % n=5500 は未確認
% MATLAB7.1用
nd3=[1 2 4 6 8 [10:5:35]]*100; % n=4000 では異常動作
% MATLAB5.3用
nd4=[1 2 4 6 8 [10:5:35]]*100; % n=4000 では異常動作
nD={ nd1 nd2 nd3 nd4 };
% データ記録用の変数
div_vs_time=[];
% 途中でフリーズしてもデータが残るように、ファイルにも逐次記録。
% そのためのファイルの準備。
Date=datestr(now,'dd-mmm-yyyy HH:MM:SS');
Date=strrep(strrep(strrep(Date,'-','_'),':','_'),' ','_');
fid=fopen(['xx_div_vs_time_' Date '.txt'],'w+');
fprintf(fid,['xx_div_vs_time_' Date '.txt']);
fprintf(fid,'\n');
fprintf(fid,'Ver Com Div Time_ave Time_min Time_max');
fprintf(fid,'\n');
fclose(fid);
com=5; % テストするコマンドの種類
% 5:contour
com_name='等高線';
ndiv=nD{mat_ver}; % 分割数の一覧
for nd=ndiv; % 各種の分割数について繰り返す。
x=linspace(-1,1,nd+1); % -1～1の間をnd分割する。
y=x;
[X,Y]=meshgrid(x,y); % x,y平面をnd×nd分割したメッシュ
Z=X.^2.*(1-X.^2)-Y.^2; % 曲面 Z(X,Y)
timedata=[]; % 毎回の演算実行時間を記録する行ベクトル
for ne=1:testrep % バラツキの平均化のため繰り返し処理
tic; % 計時開始
[dc,hc]=contour(X,Y,Z,[-0.5:0.1:0.2]); % 等高線の計算
te=toc; % 所要時間記録
% 実行時間の傾向を目視監視するため、コマンドラインに毎回印字する。
disp([mat_name ' ' com_name ' ' ...
sprintf('%4s',num2str(nd)) '×' ...
sprintf('%4s',num2str(nd)) ...
' 処理時間 ' num2str(te) ' 秒'])
timedata=[timedata te]; % 毎回の実行時間をすべて記録
drawnow; % 描画の催促
pause(2); % 描画結果を確認のため2秒間だけ表示
close all; % 1回目描画との条件合わせのため
end
% 繰返し実行の結果から、平均・最短・最長時間を計算・表示
t_ave=sum(timedata)/testrep; % 平均時間
t_min=min(timedata); % 最短時間
t_max=max(timedata); % 最長時間
disp(' ')
disp([' 平均時間：' num2str(t_ave,'%11.5f') ...
' 最短時間：' num2str(t_min,'%11.5f') ...
' 最長時間：' num2str(t_max,'%11.5f')]);
disp(' ')
% データ記録用の変数に計測結果を追記。
% （MATLABバージョン・コマンド・分割数・諸時間）
div_vs_time=[div_vs_time; ...
mat_ver com nd t_ave t_min t_max];
% 途中でフリーズしてもデータが残るように、ファイルにも逐次記録。
% ver5.3 では、曖昧に'%d'では不可。桁数の明示をしないとエラー。
fid=fopen(['xx_div_vs_time_' Date '.txt'],'a');
fprintf(fid,'%3s',num2str(mat_ver,'%1d'));
fprintf(fid,'%3s',num2str(com,'%1d'));
fprintf(fid,'%6s',num2str(nd,'%4d'));
fprintf(fid,'%12s',num2str(t_ave,'%11.5f'));
fprintf(fid,'%12s',num2str(t_min,'%11.5f'));
fprintf(fid,'%12s',num2str(t_max,'%11.5f'));
fprintf(fid,'\n');
fclose(fid);
end
% 計測結果をコマンドラインに印字
disp(datestr(now,'dd-mmm-yyyy HH:MM:SS'));
disp( [' created by ' thisfile ' on ' mat_name]);
disp('Ver Com Div Time_ave Time_min Time_max');
B=div_vs_time;
for n=1:size(B,1)
disp([sprintf('%3s',num2str(B(n,1),'%1d')) ...
sprintf('%3s',num2str(B(n,2),'%1d')) ...
sprintf('%6s',num2str(B(n,3),'%4d')) ...
sprintf('%12s',num2str(B(n,4),'%11.5f')) ...
sprintf('%12s',num2str(B(n,5),'%11.5f')) ...
sprintf('%12s',num2str(B(n,6),'%11.5f'))])
end
５.３ データのグラフ化プログラム
最終的に採用した計測数値が明示されているので、念のため掲載します。
ユーザー定義関数 make_axes_tidily を使っています。
（【MATLABユーザー定義関数】自由度が高い「subplotの変種」）
% testMatrix07.m
% MATLAB の行列計算能力、等高線計算能力の測定の結果をグラフ化
clear
close all
% 測定結果（testMatric03.m, testMatric06.m による）
% １列目：1:R2019a, 2:R2007b, 3:ver7.1, 4:ver5.3
% ２列目：1:行列式, 2:行列積, 3:逆行列, 4:固有値, 5:等高線
% ３列目：行列サイズ, 分割数
% ４列目：平均処理時間 [秒]
% ５列目：最短処理時間 [秒]
% ６列目：最長処理時間 [秒]
% MATLAB R2019a の測定結果
data01=[
1 1 100 0.00229 0.00054 0.00480
1 1 200 0.00091 0.00078 0.00118
1 1 400 0.00240 0.00184 0.00386
1 1 600 0.00741 0.00340 0.01272
1 1 800 0.01265 0.00887 0.01765
1 1 1000 0.02539 0.01786 0.03267
1 1 1500 0.06241 0.05843 0.07142
1 1 2000 0.09998 0.09678 0.10558
1 1 2500 0.16163 0.15030 0.18691
1 1 3000 0.29388 0.25771 0.33068
1 1 3500 0.41446 0.37204 0.44449
1 1 4000 0.61044 0.55276 0.69935
1 2 100 0.00130 0.00016 0.00344
1 2 200 0.00052 0.00034 0.00070
1 2 400 0.00303 0.00224 0.00368
1 2 600 0.00885 0.00681 0.01207
1 2 800 0.01956 0.01260 0.02503
1 2 1000 0.03461 0.02382 0.03907
1 2 1500 0.08758 0.08034 0.09270
1 2 2000 0.19124 0.15106 0.22087
1 2 2500 0.31125 0.28061 0.35784
1 2 3000 0.54782 0.45339 0.62738
1 2 3500 0.86390 0.76462 0.94933
1 2 4000 1.20136 1.12973 1.37760
1 3 100 0.00660 0.00052 0.02969
1 3 200 0.00216 0.00095 0.00622
1 3 400 0.00796 0.00570 0.01080
1 3 600 0.01694 0.01346 0.01912
1 3 800 0.02830 0.02311 0.03417
1 3 1000 0.05378 0.05067 0.06332
1 3 1500 0.13477 0.12015 0.16201
1 3 2000 0.27424 0.26614 0.28695
1 3 2500 0.50273 0.48888 0.51644
1 3 3000 0.80817 0.76344 0.83770
1 3 3500 1.35569 1.27317 1.44354
1 3 4000 1.87069 1.83551 1.92074
1 4 100 0.01440 0.00402 0.05400
1 4 200 0.02349 0.02117 0.02777
1 4 400 0.07628 0.07195 0.09006
1 4 600 0.22723 0.21886 0.24064
1 4 800 0.35808 0.34878 0.37077
1 4 1000 0.57684 0.55290 0.58788
1 4 1500 1.55741 1.53836 1.57451
1 4 2000 3.45890 3.43195 3.53665
1 4 2500 6.59489 6.56158 6.61371
1 4 3000 11.66058 11.53541 11.75744
1 4 3500 17.74509 17.57694 18.11189
1 5 100 0.09018 0.06668 0.11260
1 5 200 0.10157 0.09928 0.10434
1 5 400 0.14529 0.14058 0.15164
1 5 600 0.21915 0.21437 0.22513
1 5 800 0.29538 0.28526 0.30494
1 5 1000 0.41234 0.40947 0.41646
1 5 1500 0.82213 0.81714 0.82565
1 5 2000 1.46583 1.43961 1.52399
1 5 2500 2.28179 2.23877 2.33725
1 5 3000 3.25969 3.22792 3.31430
1 5 3500 4.40113 4.34763 4.47549
1 5 4000 5.81969 5.77652 5.92862
1 5 4500 7.32599 7.16071 7.62689
1 5 5000 9.38206 8.91221 9.80008];
% MATLAB R2007b の測定結果
data02=[
2 1 100 0.00021 0.00016 0.00027
2 1 200 0.00081 0.00075 0.00092
2 1 400 0.00725 0.00558 0.00943
2 1 600 0.01737 0.01550 0.01909
2 1 800 0.03563 0.03304 0.03692
2 1 1000 0.06679 0.06561 0.06792
2 1 1500 0.21936 0.21639 0.22419
2 1 2000 0.50645 0.50120 0.51059
2 1 2500 0.96557 0.95803 0.97042
2 1 3000 1.63351 1.62234 1.64124
2 1 3500 2.55805 2.54849 2.57340
2 1 4000 3.76943 3.75093 3.78148
2 2 100 0.00190 0.00020 0.00796
2 2 200 0.00191 0.00145 0.00242
2 2 400 0.01235 0.01109 0.01380
2 2 600 0.03692 0.03560 0.03789
2 2 800 0.08474 0.08313 0.08558
2 2 1000 0.16265 0.16159 0.16379
2 2 1500 0.59976 0.52985 0.63757
2 2 2000 1.37499 1.32621 1.43914
2 2 2500 2.73006 2.64790 2.81836
2 2 3000 4.24834 4.15515 4.53684
2 2 3500 6.66432 6.59504 6.77264
2 2 4000 9.84562 9.82265 9.88438
2 3 100 0.00649 0.00039 0.03086
2 3 200 0.00287 0.00238 0.00385
2 3 400 0.01665 0.01475 0.01827
2 3 600 0.04610 0.04334 0.04745
2 3 800 0.10424 0.09864 0.10840
2 3 1000 0.19463 0.19100 0.19821
2 3 1500 0.63645 0.63124 0.64269
2 3 2000 1.47813 1.47590 1.48018
2 3 2500 2.90117 2.85418 3.04439
2 3 3000 4.90081 4.88524 4.91227
2 3 3500 7.70092 7.68362 7.72392
2 3 4000 11.41932 11.39077 11.44221
2 4 100 0.01035 0.00381 0.03483
2 4 200 0.02477 0.02341 0.02847
2 4 400 0.11796 0.11756 0.11858
2 4 600 0.38650 0.38588 0.38752
2 4 800 0.68896 0.67834 0.70601
2 4 1000 1.11384 1.10799 1.12115
2 4 1500 3.11351 3.10732 3.13447
2 4 2000 6.76006 6.73526 6.77717
2 4 2500 12.74787 12.69403 12.92670
2 4 3000 20.21341 20.13866 20.28293
2 4 3500 30.33170 30.24185 30.38879
2 5 100 0.08335 0.04799 0.11860
2 5 200 0.11303 0.09013 0.17929
2 5 400 0.12802 0.11861 0.13658
2 5 600 0.17055 0.14985 0.18171
2 5 800 0.22761 0.19060 0.28158
2 5 1000 0.29125 0.27005 0.30682
2 5 1500 0.52285 0.49355 0.54245
2 5 2000 0.84889 0.83426 0.87938
2 5 2500 1.26317 1.22916 1.31744
2 5 3000 1.74134 1.72806 1.75393
2 5 3500 2.34496 2.31857 2.36714
2 5 4000 3.03385 3.01381 3.05601
2 5 4500 3.84646 3.79442 3.96145
2 5 5000 4.73073 4.65457 4.78454];
% MATLAB ver7.1 の測定結果
data03=[
3 1 100 0.00127 0.00121 0.00141
3 1 200 0.00643 0.00616 0.00719
3 1 400 0.04904 0.04460 0.06494
3 1 600 0.12576 0.12271 0.12849
3 1 800 0.25747 0.25283 0.26485
3 1 1000 0.47301 0.46722 0.48090
3 1 1500 1.46980 1.42981 1.53355
3 1 2000 3.20180 3.16059 3.23868
3 1 2500 5.95840 5.88513 6.02364
3 1 3000 9.73500 9.67580 9.79549
3 1 3500 17.47936 17.22252 17.68559
3 1 4000 22.15679 21.25188 24.53420
3 2 100 0.09726 0.00108 0.48192
3 2 200 0.01445 0.00691 0.04423
3 2 400 0.07857 0.04869 0.09058
3 2 600 0.18564 0.15302 0.22108
3 2 800 0.37529 0.35730 0.38196
3 2 1000 0.66245 0.64896 0.68863
3 2 1500 2.15878 2.08352 2.25796
3 2 2000 4.88904 4.83850 4.94647
3 2 2500 9.76765 9.64881 9.87365
3 2 3000 16.92031 16.30056 17.70760
3 2 3500 29.39530 27.54045 32.75164
3 2 4000 47.33367 46.02548 48.61029
3 3 100 0.21919 0.00286 1.08436
3 3 200 0.20068 0.01369 0.93799
3 3 400 0.28501 0.13389 0.78561
3 3 600 0.36959 0.33940 0.42750
3 3 800 0.73782 0.71540 0.75283
3 3 1000 1.31694 1.28087 1.37239
3 3 1500 3.92831 3.85538 4.06811
3 3 2000 8.72406 8.40220 9.17784
3 3 2500 15.95907 15.91925 16.00671
3 3 3000 26.45950 26.42847 26.52666
3 3 3500 41.50427 41.40494 41.69957
3 3 4000 59.58553 59.39745 59.94225
3 4 100 0.06946 0.01901 0.25213
3 4 200 0.16241 0.12322 0.19546
3 4 400 1.32765 1.29709 1.34607
3 4 600 4.61600 4.59640 4.64813
3 4 800 10.68325 10.61664 10.75417
3 4 1000 20.35686 20.31969 20.42164
3 4 1500 73.34735 70.19733 76.82302
3 4 2000 166.99092 166.53052 167.43886
3 4 2500 407.40939 360.27411 420.39288
3 4 3000 721.54401 720.45538 722.05789
3 4 3500 1418.72562 1407.31936 1439.20624
3 5 100 1.06960 0.44859 2.34288
3 5 200 0.86968 0.55586 1.20123
3 5 400 1.17377 0.87546 1.88831
3 5 600 1.23008 1.18047 1.31700
3 5 800 1.68130 1.62462 1.76356
3 5 1000 2.26599 2.15029 2.45652
3 5 1500 4.16328 4.12798 4.20845
3 5 2000 6.91324 6.80712 6.98531
3 5 2500 10.28977 10.19799 10.39823
3 5 3000 14.76414 14.70922 14.86422
3 5 3500 36.66370 21.48734 51.33949];
% MATLAB ver5.3 の測定結果
data04=[
4 1 100 0.00200 0.00000 0.01000
4 1 200 0.01000 0.01000 0.01000
4 1 400 0.40260 0.39000 0.41100
4 1 600 1.45200 1.40200 1.48200
4 1 800 3.70340 3.55500 3.79600
4 1 1000 6.91800 6.88000 6.97000
4 1 1500 23.64820 23.26400 24.29500
4 1 2000 54.83040 54.60800 55.03900
4 1 2500 107.97120 105.55200 116.42700
4 1 3000 212.24080 184.10500 269.91800
4 1 3500 289.54240 287.56300 291.86000
4 1 4000 429.14680 427.20400 432.72200
4 1 4500 619.80320 609.27600 626.07000
4 1 5000 892.76560 848.74000 1052.77400
4 2 100 0.00800 0.00000 0.01000
4 2 200 0.02000 0.02000 0.02000
4 2 400 0.74540 0.74100 0.76100
4 2 600 2.38520 2.36300 2.43300
4 2 800 5.87440 5.83800 5.91800
4 2 1000 10.73340 10.65500 10.79500
4 2 1500 34.69000 34.63000 34.74000
4 2 2000 80.72420 80.63600 80.83700
4 2 2500 157.97300 157.85700 158.22700
4 2 3000 270.51480 270.32800 270.81000
4 2 3500 427.68500 424.63000 434.55500
4 2 4000 659.19580 647.12000 677.38400
4 2 4500 906.63760 905.34200 909.18700
4 3 100 0.01200 0.00000 0.03000
4 3 200 0.03400 0.03000 0.04000
4 3 400 0.76520 0.75100 0.78100
4 3 600 3.31480 3.25500 3.42500
4 3 800 7.89300 7.78100 7.99100
4 3 1000 16.32760 15.76300 17.17400
4 3 1500 52.63580 51.49400 53.75800
4 3 2000 128.40840 122.53600 134.69300
4 3 2500 239.24800 238.04200 240.15500
4 3 3000 406.95520 406.13400 407.90600
4 3 3500 655.38040 645.58900 677.36400
4 3 4000 968.09220 950.82800 982.79300
4 3 4500 1387.02240 1371.95300 1401.71500
4 3 5000 1952.38780 1867.45500 2097.83700
4 4 100 0.04400 0.03000 0.07000
4 4 200 0.25040 0.25000 0.25100
4 4 400 3.22440 3.17500 3.26500
4 4 600 13.08700 12.91900 13.17900
4 4 800 32.50880 32.10600 33.15800
4 4 1000 64.41860 62.57000 69.97100
4 4 1500 227.06860 219.56600 239.10400
4 4 2000 537.15420 526.91700 552.50400
4 4 2500 1460.32360 1438.17800 1493.33700
4 4 3000 2308.87600 2297.12300 2333.61600
4 5 100 0.63500 0.53100 0.70100
4 5 200 0.64280 0.61100 0.69100
4 5 400 0.94140 0.82200 1.03100
4 5 600 1.40000 1.32200 1.50200
4 5 800 1.93060 1.89300 1.99300
4 5 1000 2.65980 2.60300 2.78400
4 5 1500 5.01920 4.97700 5.09700
4 5 2000 8.12000 8.04200 8.20200
4 5 2500 12.34780 12.20700 12.68800
4 5 3000 17.82580 17.60500 18.40700
4 5 3500 24.31860 24.08400 24.71500];
DD={data01,data02,data03,data04};
for nc=[1:5] % 各コマンドについて繰り返す。
if nc==1||nc==3||nc==5
% 描画用紙の準備
[hax,hsp,h_number] = make_axes_tidily ...
('A4','land',[2 1],[250 75],[30 24 30 10]);
delete(h_number); % 不要な注記を消す。
axes(hsp(1))
if nc==1||nc==3
text(0,10,'MATLAB の処理能力の進化（行列サイズと処理時間）','FontSize',16, ...
'HorizontalAlign','center')
else
text(0,10,'MATLAB の処理能力の進化（分割数と処理時間）','FontSize',16, ...
'HorizontalAlign','center')
end
if nc==1
text(-140,10,'図３','FontSize',20)
elseif nc==3
text(-140,10,'図４','FontSize',20)
elseif nc==5
text(-140,10,'図５','FontSize',20)
end
axes(hax(1))
ax=gca;
lincol=ax.ColorOrder; % 標準線色（７色）の取得
end
for nv=[1:4] % MATLAB の各バージョンについて繰り返す。
D=DD{nv};
data=D(D(:,2)==nc,:);
data(1,:)=[]; % 100×100は誤差が大きすぎるので削除
col=lincol(nv,:);
if nc==1||nc==2
hf=hax(nc);
elseif nc==3||nc==4
hf=hax(nc-2);
else
hf=hax(nc-4);
end
[hl(nv)]=plot_data(data,col,hf); % 線の描画（ローカル関数呼び出し）
end
xlim([0.0001 10000])
ylim([80 6000])
grid on
if nc==5
ylabel('分割数（ n×n の n ）')
else
ylabel('行列サイズ（ n×n の n ）')
end
aa=gca;
aa.XAxis.FontSize=9.5;
aa.YAxis.FontSize=9.5;
aa.GridColor='k';
aa.GridAlpha=1;
aa.MinorGridColor='k';
aa.MinorGridAlpha=0.3;
aa.MinorGridLineStyle='-';
aa.LabelFontSizeMultiplier=1.3;
legend(hl,'MATLAB R2019a','MATLAB R2007b', ...
'MATLAB ver7.1','MATLAB ver5.3','Location','southeast')
if nc==1||nc==3
xticklabels('')
else
xlabel('処理時間 [秒]')
end
switch nc
case 1
text(1.5e-4,3.5e3,'行列式','FontSize',16);
case 2
text(1.5e-4,3.5e3,'行列積','FontSize',16);
case 3
text(1.5e-4,3.5e3,'逆行列','FontSize',16);
case 4
text(1.5e-4,3.5e3,'固有値','FontSize',16);
case 5
text(1.5e-4,3.5e3,'等高線（8レベル）','FontSize',16);
end
if nc==5
axes(hax(2));
axis off;
end
end
% =================================================
function [hl]=plot_data(data,col,hf)
nn=data(:,3); % サイズ（列ベクトル）
av=data(:,4); % 平均時間
mn=data(:,5); % 最短時間
mx=data(:,6); % 最長時間
axes(hf);
% 最短時間と最長時間のグラフ
loglog(mn,nn,'.-','Color',col,'LineWidth',0.5);
hold on
loglog(mx,nn,'.-','Color',col,'LineWidth',0.5);
% 最短時間と最長時間に囲まれた領域を淡く塗り潰す。
fill([mn;flipud(mx)],[nn;flipud(nn)],col,'FaceAlpha',0.3,'EdgeColor',col);
% 平均時間のグラフ
hl=loglog(av,nn,'o-','Color',col,'LineWidth',2);
end
５.４ 適正K値のグラフ
最終的な計測数値が明示されているので、念のため掲載します。
% testMatrix05.m
% お試し計算用のn×nの巨大行列をA=rand(n)/Kで作るとき、
% 行列サイズnと下記のKの値の関係をグラフ化する。
% データは、コマンドラインから下記を打ち込み、Kやnの値を変えながら、
% 根気よく取得したものをそのまま使用。
% K=10; n=3000; disp(datestr(now,'dd-mmm-yyyy HH:MM:SS')),...
% pause(0.5); A=rand(n)/K; disp(det(A))
%
% Kの最適値 Kbest : |A|の値がほぼe+00オーダになるようなKの値
% Kの許容下限値 Kmin : |A|の値が絶対に±InfにならないKの下限値
% Kの許容上限値 Kmax : |A|の値が絶対にアンダーフロー0にならないKの上限値
clear all
close all
% 計測結果のデータ（貴重で重要。データ採りには相当な時間を使った。）
% A=rand(n)/K としたときの最適なK(Kbest)、
% 下限・上限の K(Kmin,Kmax ...MATLAB R2019a他、Kmin5,Kmax5 ...MATLAB ver5.3)。
n= [100 200 300 400 500 700 1000 1500 2000 3000 4000 5000];
Kbest=[1.78 2.49 3.04 3.51 3.94 4.64 5.55 6.79 7.83 9.59 11.08 NaN];
Kmin =[0.00153 0.0733 0.290 0.602 0.958 1.692 2.74 4.24 5.51 7.59 9.29 NaN];
Kmax =[3039 102 36.1 22.4 17.3 13.41 11.65 11.15 11.51 11.54 11.84 NaN];
Kmin5=[NaN NaN NaN NaN NaN NaN 2.74 4.24 5.51 7.59 9.29 10.76];
Kmax5=[NaN NaN NaN NaN NaN NaN 11.65 11.15 11.34 12.27 13.33 14.36];
figure(1)
h1=loglog(n,Kbest,'o','LineWidth',2);
hold on
grid on
xlim([50 12000])
ylim([0.9 15])
h2=loglog([1 10000],[0.1 10]*1.75,'LineWidth',2);
h3=loglog(n,Kmin,'o-','LineWidth',1);
h4=loglog(n,Kmax,'o-','LineWidth',1);
h5=loglog(n,Kmin5,'.-','LineWidth',1,'MarkerSize',14);
h6=loglog(n,Kmax5,'.-','LineWidth',1,'MarkerSize',14);
aa=gca;
aa.GridColor='k';
aa.GridAlpha=1;
aa.MinorGridColor='k';
aa.MinorGridAlpha=0.3;
aa.MinorGridLineStyle='-';
text(32,15.5,'図１','FontSize',20);
text(60,10^0.95,'Kbest .... |A| の指数部がほぼ e+00 となるK')
text(60,10^0.90,'Kmin .... |A| が ±Inf とならないK')
text(60,10^0.85,'Kmax ... |A| が 0 (underflow) とならないK')
title([' 乱数行列 A=rand(n)/K について、' ...
'|A| をパラメータにした n と K の関係'],'FontSize',11);
xlabel('行列サイズ n （正方行列の一辺）');
ylabel('K');
legend([h1 h2 h3 h4 h5 h6], ...
{'Kbest' 'K=0.175*sqrt(n)' 'Kmin(R2019a)' 'Kmax(R2019a)'...
'Kmin(ver5.3)' 'Kmax(ver5.3)'},'Location','southeast');