【MATLAB今昔】その処理能力の進化のほどを可視化してみた

## １．はじめに

　BASIC を使ってモタモタと数値計算していた Windows3.1 の時代に、突然、MATLAB の凄さを見せられたときの感動は今でも忘れません。それ以来、MATLAB にはずっとお世話になってきていますが、忘れられないのは感動だけではありません。その当時に独り善がりで会得した「利用上の心得」も記憶に残り、未だにその影響から抜け出せずにいます。その心得の一つが「<font color='#aaa'>実用上の行列サイズの上限は、精々250行×250列程度</font>」（風評被害防止のため淡色表示）というものです。誰に感化された訳でもなく、この程度のサイズの行列を扱う中で、なかなか計算結果が出てこないことに業を煮やして、身体で感じ取ったノウハウです。その後、あまり大きな行列を扱う機会は訪れず、この先入観を是正することができないまま現在に至ってしまいました。

　ところが、現行のバージョン（と言ってもチョット古い 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<br>
+ **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の<br>バージョン|同時代のパソコン|
|:-----------|:------|
|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));

<img width="400" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/932484/a2dfaaac-4047-9852-561c-ba1636e88fff.png">

　このようなシンプルな関係式であれば、解析的に簡単に導けるように思えますが、残念ながらまだ達成できておりません。なお、図１には、最適な $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 分割まで評価して何のメリットがあるのか？という疑問は残りますが、処理能力の把握のための試練だと割り切って強行しました。

　描かせる等高線の本数は処理時間に直接影響してくるはずですが、今回、その関係について深くは考えていません。特別な理由はなく、適度な本数として８レベルを選んでいます。

<img width="400" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/932484/4832e737-f318-6967-9d9d-84836f39e3f6.png">

**評価用のプログラム**

　基本的な処理内容は、コマンドラインからでも操作できる次のようなものです。

```
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 の分割数を大きくしていきながら繰り返します。行列演算の評価と同じく、全く同一の計算を５回にわたって実行させ、処理時間の平均値と最短・最長値を求めます。


## ３．評価結果

### ３.１　演算能力

　データ採取用プログラムから得られた情報を整理して、次の一連のグラフができました。**図３**の上は行列式、下は行列積、**図４**の上は逆行列、下は固有値について、それぞれ、行列サイズと処理時間の関係を示したものです。また、**図５**は等高線描画について、メッシュ分割数と描画時間の関係を示したものです。

<img width="800" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/932484/35de3d48-5d01-2cbf-3053-d10c7d71d9e3.png">

<img width="800" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/932484/a3f2e0ac-1ba2-bb1f-35c8-83f46bda6a45.png">

<img width="800" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/932484/dd7d77f4-f978-7a6b-928e-2f2e070c23ff.png">

　全体的には、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 関係の図はすべてこの方式によっています。画面上に表示されたものをそのまま最終成果物として利用できるという今の状態が、私にとっては究極の表現形態であると思っています。

<img width="400" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/932484/ee460456-0d5f-7537-a6f4-f73fc6221baf.png">

<img width="400" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/932484/43936435-7e7e-43c3-325d-83bbdda1083a.png">

<img width="400" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/932484/13e34da4-7274-d885-47ed-f8f657349931.png">

<img width="400" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/932484/07be9e30-0f41-7337-9c38-5e2b7c8a4af9.png">

<img width="400" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/932484/2a18330b-4485-23bf-e1c1-57dad140bef0.png">

<img width="400" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/932484/69cb477b-9d84-c7c8-2d1f-8517f9be40dc.png">

## ４．おわりに

　以上で新旧 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 の倍数でもある）にしているので、たまたま成立しているだけです。

## ５．プログラム

### ５.１　データ採取用プログラム（行列計算）

```MATLAB
% 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
```

### ５.２　データ採取用プログラム（等高線描画）

```MATLAB
% 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の変種」](https://qiita.com/tsubolabo/items/5b6149b7c7d6f71391bc)）

```MATLAB
% 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値のグラフ

　最終的な計測数値が明示されているので、念のため掲載します。

```MATLAB
% 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');
```