ドレスデン午前４時

いろいろとややこしいんですが、ポスドクは基本的にTV-L-13という分類になるらしい？です。サラリーはいろんな税金さっ引いて(元の48%くらいなくなる？詳しく知らない確認したことないw)、だいたい1800-1900ユーロくらいです。年で取れる有休は基本的に２９日です。使わないとこれまたいろいろややこしいことになります。 西と東でこれまた若干の違いがあります。 それでいて労働時間は９時５時で、土日は完全にフリーです。 平日の５時以降はもっといたけりゃそれでも構いませんが、 ７時過ぎるとあんまり人がいません。 特に冬は４時に空が暗くなるとちらほら人が帰り始めます。 家賃はそれこそピンキリですけど、私のフラットはちょっと例外的で格安(200ユーロで水代・電気代・ガス代・ネット代フルカバー)です。学生と共同ですけど。ひとりで住みたいとなるとドレスデンだと大体400ユーロ？/月で５０平米。 昼食代だけは少し高く、週５日で30ユーロ（ちなみに学食にビールサーバーがあって真っ昼間にビールを飲むことも可能です。自分はやらないけどｗ）。夜は自炊しているので週15ユーロです。 ドイツの場合は（というか欧米全般）、修士以上の学生はポスドクの半分の給料？（無償の奨学金？奨励金？）をもらってるので親からの仕送りとかないと思う（ドイツ人に聞いたことない）。

手取りの給料を見ると安っと思うかもしれませんが、労働環境や充実した（有名無実でない）有給、安い酒と食材、てきどにかんそうした気候を考えると（あと、言語の難しさをのぞくと）いい国だと思います。ただパーマネントのポジションを取るのはかなり熾烈だと思います。もともと永住するつもりがサラサラないうえに、若い人はほとんど英語が通じるので全くドイツ語を話そうという気が起らず今にいたってます。

１年半で４報の論文を手がけ（１報は表紙を飾りました。残り２報はfirst author, 残りは共著）ることができました。

条件にもよるけど誤差が１％増加するだけで１０００倍速くなるコードを作ることができた。びっくりしたと同時にうれしかった。この原因が分かれば方法論の論文も書ける。

計算機の資源が豊富だから一気に計算をやってしまえるのではあるが、その分解析も忙しくなる。Perlをちゃんと勉強しようかな。
ここで整理整頓の能力が問われる。忙しいとかいいつつ、7時頃帰宅です。それでも施設では2番目くらいに遅いです。
それでいてみんなは残すべき結果は残しています（自分も：D）。

明日は祭日。何の祭日かは不明。

偏微分方程式を数値的に解くパッケージMUDPACK(泥パックではありません)。偏微分方程式を解く数値計算のソフトは応用範囲が広く
需要がそこそこあると思うんですが、ネットで調べても日本語のページは全く見当たらなかったので（検索の仕方が悪い？^ ^;）、
超々初心者ど素人ながら防備録も兼ねて使い方をメモしたいと思います。
日本語へたくそで、間違いもあるかもしれないですけど、そこはどうか。。。

そもそもこのパッケージを使おうと思ったのは、
Poisson方程式を数値的に解いて、
電位の空間分布を調べる方法を習得したかったからです。

Poisson方程式は、∇²φ = -ρ/ε₀であらわされます。∇はナブラと読み、∇ = (∂/∂x +∂/∂y +∂/∂z )です。
φは電位で、φ=φ(x,y,z)のように空間座標の関数です。
ρは電荷密度でρ=ρ(x,y,z)のように空間座標の関数です。。
Poisson方程式は、与えられた電荷密度の分布ρと境界条件からφ電位の分布を求める際に用いられる方程式で、
これを（数値的にあるいは分析的に）解くことにより空間における電位の分布（並びにその空間微分である電界）を求めることができます。

MUDPACKのパッケージ群に関する情報は
http://www.cisl.ucar.edu/css/software/mudpack/　にまとめられています。
ライセンスに同意して下記のフォームを入力すると
http://www.cisl.ucar.edu/css/software/mudpack/mud.form.html
ソフトを入手することができます。

この記事で問題とする方程式はPoisson方程式であり、
∇²φ = ∂²φ/∂x² +∂²φ/∂y² +∂²φ/∂z² = -ρ/ε₀の形をしています。
えーと、だから3次元で、分離可能（x、y、zの各次元で独立してる。）で、なおかつ電位φと電化密度ρは
実数なので、用いるサブルーチン（サブプログラム）は2nd order/real 3D separable　(mud3p)か
4th order/real 3D separable (mud34p)となります。

4th orderと2nd orderの違いがわからないので（たぶん精度の違いだと思う）、とりあえずは
シンプルであろう2nd order(mud3p)を用いて話を進めたいと思います。
サブルーチンmud3sp.f　の仕様は　下記のURLで詳しく解説されています。
http://www.cisl.ucar.edu/css/software/mudpack/mud3sp.txt
サブルーチンmud3sp.f　を用いたサンプルコード(test code)は　下記のURLで見ることができます。
http://www.cisl.ucar.edu/css/software/mudpack/tmud3sp.txt
また、なんのこっちゃわからないんですけど予備のコード(residual code, )は　下記のURLで見ることができます。
http://www.cisl.ucar.edu/css/software/mudpack/resm3sp.txt

それではサブルーチンmud3sp.f　の仕様について見ていきましょう。
これがまた小難しく、見た瞬間に意識が飛びそうになりますが、辛抱強く一つづつ。。。。

(1)　2行目はこのファイル名です。　.dはドキュメントを表す拡張子です。
(2)　10行くらいコピーライトの記述が続きます。
(3)　サブルーチンmud3sp.f　はインプット（引数）としてiparm, fparm, work, cofx, cofy, cofz, bndyc, rhs, phi, mgopt, ierror が必要です。
(4)　サンプルコードは”tmud3sp.f”です。と言ってます。
(5)　このパッケージ群を開発・コーディングしたジョン・C・アダムスさんの署名です。
(6)　サブルーチンmud3sp.fを実行する際に必要となる親プログラムはmudcom.f　とあります。
(7)　目的：　サブルーチンmud3sp.fは自動的に空間を区分化し、2ndオーダーの有限差分法を、ある境界条件のもと（ジリクレ境界条件、周期的境界条件など）、
箱形の空間における分離可能な楕円型偏微分方程式に対して適用します。

このサブルーチンでは、以下のような楕円型偏微分方程式を想定しています。
∇²φ =
c_xx∂²φ/∂x²  + c_x∂φ/∂x + c_exφ(x, y, z)+
c_yy∂²φ/∂y² + c_y∂φ/∂y + c_eyφ(x, y, z)+
c_zz∂²φ/∂z² +  c_z∂φ/∂z + c_ezφ(x, y, z) = ρ(x, y, z) の形をしています。

ここではc_xx、c_x、c_ex、c_yy、c_y、c_ey、c_zz、c_z、c_ez、はそれぞれ既知の実数の係数です。
φ(x, y, z) は未知です。この方程式の解です。つまりこれから求めようとする関数です。
ρ(x, y, z) は既知の関数です。
c_xxとc_yyとc_zzは実数の正の値でなくてはなりません。

サブルーチンmud3sp.fは分離可能な偏微分方程式の数値計算を前提としていますが、
この「分離可能(separable)」とは何かというと簡単に言うと
方程式に∂²φ/∂y∂z　とか ∂²φ/∂x∂y　とかいう項がなくて変数分離できるってことです（たぶん）。
そのような項が紛れ込む際は、muh3.fやmud3.fを使うとあります。

(8)　メッシュの記述について
長方形型の箱を想定し、xyz各方向にnx、ny、nz個のグリッドを設けます。
長方形のBoxは

 x_a<= x <= x_b,  y_c <= y <= y_d,  z_e <= z <= z_f　で定義されます。

グリッド幅は
dlx = (x_b - x_a) /(n_x-1),  dly = (y_d - y_c)/(n_y-1),  dlz = (z_f - z_e)/(n_z-1)で与えられます。

よって各グリッドは
x_i =  x_a + (i-1)*dlx,  i = 1, 2, ..., n_x
y_j =  y_c + (j-1)*dly,  j = 1, 2, ..., n_y
z_k =  y_e + (k-1)*dlz, k = 1, 2, ...,  n_z
で与えられます。

ちなみに言わずもがなですが、定義域？は、(n_x-1)*(n_y-1)*(n_z-1)個の小さな長方形で埋め尽くされることになります。

(9)　プログラム言語についてです。フォートラン90(f90)あるいはフォートラン77(f77)　とあります。
(10)　ポータビリティ（移植可能性）について　いろんなプラットフォームでf90とf77でコンパイルしてきた　とあります。
(11)　方法論については文献[1]と[9]を、パフォーマンスについては[1]と[2]と[9]を参照にするべし　とあります。
(12) 　区分化と解について。偏微分方程式と境界条件は、すべてのグリッドレベルにおいて、二次の有限差分の式を用いて、自動的に区分化される。粗いグリッドレ ベルにおける対角のdominanceは、ゼロでない一次の項存在下において、二次の係数を調整するによって保持される（必要な場合は）。
(13)　結果的に得られるtri-diagonalなlinearなブロック行列の系はマルチグリッドの繰り返しを使って近似される。これについては文献[10, 11, 13, 15, 16, 18]を参照。
(14)　Mudpack バージョン5.0.1は完全にウェイトを置いたresidualな制限のみを用いる。
(15)　デフォルトはcubicな延長とw(2,1)を含む。　w(2,1)ってなんだ？
(16)　これらは　選択されたマルチグリッドのオプションにしたがい上書きされる。（詳しくは"mgopt"オプションを参照。）
(17)　maximum relative differencesに基づくerrorのコントロールは可能で、(a)full multigrid cycling(fmg)、(b)cycling beginning、(c)が選択可能である。　
(18)　Relaxationの方法に関しては、Gauss-Seidal点、line relaxation, planar relaxation(3次元の異方性をもつ問題に対して)が提供されている。
(19)　上記(18)のすべての方法は、cray vectorizationと改善された収束率について分岐点（赤/黒）、線、あるいは面に基づく並べ替え（順序？）を用いている。参考文献は[14]
(20) 　高次の解について：　もしマルチグリッドのサイクルが2nd orderのestimateで終わると（すなわち区分化レベルerrorに達したとき？）、解はさらに４次のorderのestimate へ、"Deferred　correction"のテクニックを用いて、改善することができる。
(21)　解を入れていた配列の値は、切り捨て誤差に対する４次のorderの近似を生み出すために使われる。
(22)
(23)
(24)
(25)　参考文献の列挙

ここからが実質的に必要となる知識
(26)　引数の記述
(27)　iparm　２２の整数からなるベクトル値。効率的に整数の引数を渡すために使われる。iparmははじめにmud3spにおいてsetされる。定義は以下のよう...

(27-1)　intl  = iparm(1):
(27-2)　nxa = iparm(2): x = x_a におけるxy平面の境界条件。
　　　　0なら関数p(x,y,z)はx方向について周期的境界条件に従う。つまりすべてのx,y,zについてp(x,y,z) = p(x+x_b-x_a,y,z)
           1なら関数p(x_a,y,z)は特定される。　（これはphi(1,j,k)を介したインプットでなくてはならない。）
　　　　２ならx = x_a において混合型の境界条件に従う。（"bndyc"の項のkbdy = 1を参照しましょう。）
(27-3)　nxb = iparm(3): x = x_b におけるxy平面の境界条件
           0なら関数p(x,y,z)は周期的境界条件に従う。つまりすべてのx,y,zについてp(x,y,z) = p(x+x_b-x_a,y,z)
           1なら関数p(x_a,y,z)は特定される。　（これはphi(n_x,j,k)を介したインプットでなくてはならない。）
　　　　２ならx = x_b において混合型の境界条件に従う。（"bndyc"の項のkbdy = 2を参照しましょう。）
(27-4)　nyc = iparm(4): y = y_c におけるxz平面の境界条件
           0なら関数p(x,y,z)はy方向について周期的境界条件に従う。つまりすべてのx,y,zについてp(x,y,z) = p(x,y+y_d-y_c,z)
           1なら関数p(x,y_c,z)は特定される。　（これはphi(i,1,k)を介したインプットでなくてはならない。）
　　　　２ならy= y_c において混合型の境界条件に従う。（"bndyc"の項のkbdy = 3を参照しましょう。）
(27-5)　nyd = iparm(5): y = y_d におけるxz平面の境界条件
           0なら関数p(x,y,z)はy方向について周期的境界条件に従う。つまりすべてのx,y,zについてp(x,y,z) = p(x,y+y_d-y_c,z)
           1なら関数p(x,y_d,z)は特定される。　（これはphi(i,n_y,k)を介したインプットでなくてはならない。）
　　　　２ならy= y_d において混合型の境界条件に従う。（"bndyc"の項のkbdy = 4を参照しましょう。）
(27-6)　nze = iparm(6): z = z_e におけるxy平面の境界条件
           0なら関数p(x,y,z)はz方向について周期的境界条件に従う。つまりすべてのx,y,zについてp(x,y,z) = p(x,y,z_f-z_e)
           1なら関数p(x,y,z_e)は特定される。　（これはphi(i,j,1)を介したインプットでなくてはならない。）
　　　　２ならz= z_e において混合型の境界条件に従う。（"bndyc"の項のkbdy = 5を参照しましょう。）　
(27-7)　nzf = iparm(7): z = z_f におけるxy平面の境界条件
           0なら関数p(x,y,z)はz方向について周期的境界条件に従う。つまりすべてのx,y,zについてp(x,y,z) = p(x,y,z_f-z_e)
           1なら関数p(x,y,z_e)は特定される。　（これはphi(i,j,z_f)を介したインプットでなくてはならない。）
　　　　２ならz= z_f において混合型の境界条件に従う。（"bndyc"の項のkbdy = 6を参照しましょう。）　　
/******/ グリッドサイズの引数 /******/
(27-8)　ixp = iparm(8): 1より大きな整数、それはx方向のグリッドの数を定義する際に使われる。（詳しくはnx = iparm(14)を参照のこと）。”ixp+1”はマルチグリッドのサイクルの最中のもっとも粗いxに関するグリッド数になる。ixpはできるだけ小さ くなければならない。推奨される値は素数である2か3あるいは（おそらく）5。大きな値は、マルチグリッドの収束率をひどく鈍らせうる、とりわけもしx方 向に関するline relaxationが使われない場合は。もしixpが2より大きい場合、それは2または小さな奇数であるべきである、なぜならixpの因数2のパワー は、nx=iparm(14)を変えずにiex=iparm(11)を増やすことで取り除くことができるからである。
(27-9)　jyq = iparm(9): 1より大きい整数で、それはy方向のグリッドの数を定義する際に使われる。（詳しくはny = iparm(15)を参照のこと）。”jyq+1”はマルチグリッドのサイクルの最中のもっとも粗いxに関するグリッド数になる。jyqはできるだけ小さ くなければならない。推奨される値は素数である2か3あるいは（おそらく）5。大きな値は、マルチグリッドの収束率をひどく鈍らせうる、とりわけもしy方 向に関するline relaxationが使われない場合は。もしjyqが2より大きい場合、それは2または小さな奇数であるべきである、なぜならjyqの因数2のパワー は、ny=iparm(15)を変えずにjey=iparm(12)を増やすことで取り除くことができるからである。
(27-10)　kzr = iparm(10): 1より大きい整数で、それはz方向のグリッドの数を定義する際に使われる。（詳しくはnz = iparm(16)を参照のこと）。”kzr+1”はマルチグリッドのサイクルの最中のもっとも粗いzに関するグリッド数になる。kzrはできるだけ小さ くなければならない。推奨される値は素数である2か3あるいは（おそらく）5。大きな値は、マルチグリッドの収束率をひどく鈍らせうる、とりわけもしz方 向に関するline relaxationが使われない場合は。もしjyqが2より大きい場合、それは2または小さな奇数であるべきである、なぜならkzrの因数2のパワー は、nz=iparm(16)を変えずにkez=iparm(13)を増やすことで取り除くことができるからである。
(27-11)　iex = iparm(11): 正の整数で2を基底とする指数関数の指数、これはx方向のグリッドの数を定義する際に使われる。（詳しくはnx = iparm(14)を参照のこと）。iexが50よりも小さいことが要求される。マルチグリッディングのサイクルの効率を高めるため、iexはできるだけ 大きいことが求められる、そしてnx=iparm(14)を定義したらグリッドサイズの制限内でixp=iparm(8)はできるだけ小さいことが求めら れる。
(27-12)　jey = iparm(12): 正の整数で2を基底とする指数関数の指数、これはy方向のグリッドの数を定義する際に使われる。（詳しくはny = iparm(15)を参照のこと）。iexが50よりも小さいことが要求される。マルチグリッディングのサイクルの効率を高めるため、jeyはできるだけ 大きいことが求められる、そしてny=iparm(15)を定義したらグリッドサイズの制限内でjyq=iparm(9)はできるだけ小さいことが求めら れる。
(27-13)　kez = iparm(13): 正の整数で2を基底とする指数関数の指数、これはz方向のグリッドの数を定義する際に使われる。（詳しくはnz = iparm(16)を参照のこと）。kezが50よりも小さいことが要求される。マルチグリッディングのサイクルの効率を高めるため、kezはできるだけ 大きいことが求められる、そしてnz=iparm(16)を定義したらグリッドサイズの制限内でkzr=iparm(10)はできるだけ小さいことが求めら れる。
(27-14)　nx = iparm(14):区間 x_a≦ x ≦x_b（境界も含む）における均等に区切られたグリッドの数。nxは次の形式を満たす。
nx = ixp*2^(iex-1)+1. ここでixp = iparm(8), iex = iparm(11)
(27-15)　ny = iparm(15):区間 y_c≦ y ≦y_d（境界も含む）における均等に区切られたグリッドの数。nyは次の形式を満たす。
ny = jyq*2^(jey-1)+1. ここでjyq = iparm(9), jey = iparm(12)
(27-16)　nz = iparm(16):区間 z_e≦ z ≦z_f（境界も含む）における均等に区切られたグリッドの数。nzは次の形式を満たす。
nz = kzr*2^(kez-1)+1. ここでkzr = iparm(10), iex = iparm(13)
/******/　例　/******/
33 x 65 x 97 のグリッドを持つ系の解を求めたいとする。その時は、ixp=2, jyq=4, kzr=6, iex=jey=kez=5 とおける。
よりよい定数の選択はixp=jyq=2, kzr=3, iex=5, jey=kez=6である。
/******/　注意　/******/
Gがn_x x n_y x n_zの細かいグリッドを持ち、nを n=max₀(iex, jey, kez)とすると、Mudpackでは従順のマルチグリッドサイクルは以下のように実装される。

G(1) < G(2) < ... < G(k) < ... < G(n) = G.
それぞれのG(k)　( k = 1, 2, ... , n)は m_x x m_y x m_z　のグリッド点を持ち、グリッド点は、
m_x(k) = ixp * 2^{(max0(iex+k-n,1)-1)}+ 1
m_y(k) = jyq * 2^{(max0(jey+k-n,1)-1)}+ 1
m_z(k) = kzr * 2^{(max0(kez+k-n,1)-1)}+ 1
(27-17)　iguess = iparm(17)
０：　周期的境界条件のいかなる初期値も与えられていない、かつ（あるいは）、fullなマルチグリッドサイクルをもっとも粗いレベルから始めたいとき。
１：　周期的境界条件の初期値がもっとも細かいレベルでphiとして(phiについては下記参照)与えられていているとき。この場合、サイクルの始まり、または再開は最も細かいグリッドで開始される。
/********/　iguessについてのコメント　/********/
iguess=0 の設定はfullなマルチグリッディングあるいは”fmg”のサイクルを強いる。phiはすべてのグリッド点で開始？されていなくてはならない。もしよい 値がないならばすべての非ジリクレのグリッド点で０にセットすることもできる。iguess=0の時のphiに設定されたその値は渡され、低次元化され、 もっとも粗いグリッドレベル（サイクルの始まる）での周期的境界条件の初期値を作るのにつかわれる。

もしiguess=1ならば、inputの値phiはサイクルがはじまる最も細かいグリッドレベルでの周期的境界条件の初期値となる。
このオプションはもしとてもよい初期値が入手可能なときのみ使われる。（たとえば前のiguess=0のcallから再開する時）
/********/　時間依存の問題　/********/
楕円型偏微分方程式を各時間のステップでmarching問題の形式で解きたいとする。
l(p(t)) = r(t)
ここで、微分演算子l　は時間に依存しないとする。"p(t)"は解である。"r(t)"は時間tにおける右側の項である。
dtを時間ステップの最小の幅とする。するとp(t)は、l(p(t+dt)) = r(t+dt)を解くときに、intl = 1 のもとでp(t+dt)への初期値として使える。
はじめの2ステップののち、（つづける場合は？）、"correction"項を定めるのがよかろう。つまり、
e(t,dt) = p(t+dt) - p(t)
これは明らかに次の式を満たす。
l(e(t,dt)) = r(t+dt) - r(t).
これはiguess=0、intl=1のもとでで解かれるべきである。e(t,dt)の境界条件は、時刻t+dtで得られた値から求めた時刻tにおける値をsubtraceして、そこから得られたp(t)の境界条件から求められる。
たとえばもし、
あるxの境界で　d(p(t))/dx = f(t), d(p(t+dt))/dx = f(t+dt);　
そして　e(t,dt)が
d(e(t,dt))/dx = f(t+dt) -f(t)　を満たす場合。
e(t,dt)は0.0（at non-Dirchlet points）あるいは（もしp(t-dt)が保存されているなら）p(t)-p(t-dt)　であらわされる。
iguess = 0のもとで、これらの値は、もっとも粗いグリッドレベルでのe(t,dt)の初期値として働くであろう。
こ のアプローチは以下の場合に有効です。それは、もっとも粗いレベルから始まるマルチグリッディング　サイクルのfull sequenceが、e(t,dt)を解くときにいつも呼び出される場合。e(t,dt)の数桁の高精度な値（普通はp(t)よりもずっと小さい）は、最 終的な近似においてさらにいくつかの桁のよい精度をもたらすであろう。
p(t+dt) = p(t) + e(t, dt).
時間について積分する際、このアプローチを使えば、より高い精度を与え、すべての時間ステップについてp(t)をp(t+dt)の初期値としてつかうことができるでしょう。それは付加的なストレージを必要とします。
もし微分演算子"l"が時間依存性を持つ場合（偏微分方程式の係数を介してあるいは微分形式の境界条件の係数を介して）、
l(t+dt)(p(t+dt)) = r(t+dt) をintl=0とともにすべての時間ステップについて解くとき、p(t)をp(t+dt)の初期値として使う。(区分化は毎回新しい"t"ごとに繰り返されなければならない。)
なおiguess=0（そのときp(t)は計算の始まるもっとも粗いグリッドレベルでの初期値になる。）またはiguess=1（このときp(t)は計算サイクルの固定されたもっとも細かいグリッドレベルでの初期値になる。）の場合でも共に実行することができる。

(27-18) 　maxcy = iparm(18): maxcyはtolmax=fparm(7)=0.0の時（no error control）の最も粗いグリッドレベル（ (ixp+1) x (jyq+1) x (kzr+1) ）からもっとも細かいグリッドレベル（ n_x x n_y x n_z） まで実行されるサイクルの数（整数）。tolmax=fparm(7) > 0.0 がインプット( error control)のとき、maxcyは最も細かいグリッドのレベルともっとも粗いグリッドのレベルのサイクルの数の限界である。いかなるケースにおいて も、せいぜい最大でも、maxcy*(iprer+ipost)回のrelaxation sweepが最もグリッドの細かいレベルで(下記のiprer=mgopt(2), ipost=mgopt(3)を見よ。)行われる。マルチグリッドの繰り返しが正しく働いている時、ほんの数回のサイクルで収束する。大きなmaxcyは 必要とされるべきではない。

(27-19)　method = iparm(19): 0の時、red/black orderingでgauss-seidel法（pointwise relaxation）が用いられる。これはmud3spにおいて提供される唯一の緩和法である。Lineあるいはplanar relaxationはワークスペースを節約するというmud3spの利点を失わせる。　もしlineあるいはplanar relaxationが必要ならばmuh3あるいはmud3を用いるべし。　引数methodはただmud3spの目的に注意を払うために一つしか種類の ない引数として使われる。

(27-20)　length = iparm(20)　はワークスペースの長さ（ベクトルの長さ）。これはmuh3やmud3といった非分離型のソルバーで必要とされるワークスペースと比べてひどく小さい。　length = 7*(n_x+2)*(n_y+2)*(n_z+2)/2であるがいつもこれで十分であるわけではない。正確な最小のlengthの値はグリッドサイズの引数に依存する。それは現在の入力引数に対してiparm(20)=0と設定してiparm(21)にプリントするようにmud3spを呼ぶことで既定の値にできる。またはf90の場合は後に続くmud3spのコールにおいてiparm(21)の値を使って動的にwork spaceを確保することでそれができる。

/********/　fparmについて　/********/
fparmは浮動少数からなるベクトルでその長さは８である。効率的に浮動小数の引数を渡すために使われる。fparmはmud3spにおいて以下のように定義される。
(28-1) x_a = fparm(1), x_b = fparm(2)はx方向の範囲。x_aはx_bよりも小さくなくてはならない。

(28-2) y_c = fparm(3), y_d = fparm(4)はy方向の範囲。y_cはy_dよりも小さくなくてはならない。

(28-3) z_e = fparm(5), z_f = fparm(6)はz方向の範囲。z_eはz_fよりも小さくなくてはならない。

(28-4) tolmax = fparm(7):

Advances in Organic Photoconductor Technology Chem. Rev. ASAP

du -hs

kompare