Fortran2003
Pivot選択のためにちょっと余分な計算をしなければなりません。 実行結果 ソース・プログラム module m_mat implicit none integer, parameter :: kd = kind(0.0d0) contains subroutine lu2(a, b, x) real(kd), intent(in) :: a(:, :), b(:) real(kd), inte…
実行結果 ソース・プログラム module m_mat implicit none integer, parameter :: kd = kind(0.0d0) contains subroutine lu2(a, b, x) real(kd), intent(in) :: a(:, :), b(:) real(kd), intent(out) :: x(:) real(kd), allocatable :: c(:, :) c = a ! F20…
プログラム中で１から１０までの整数列が欲しくなることがありますが[(i, i = 1, 10)]と配列生成子でつくるのも、捨て変数 i を宣言しないといけないので面倒です。Haskell のように 1..10 のように生成できるとうれしいのですが、Fortranならユーザー定義演…
もともとの行列が上書きされてしまうので、中間サブルーチンを置いて作業行列などを確保します。枢軸（pivot)用のベクトルもここで用意します。 実行結果 500*500の行列の線型方程式を20回解いて残差の二乗平均を出力。おおむね数値誤差程度のまで求められて…
ローレンツ型の関数を‐1から1の範囲で、いくつかの分割数に対して台形公式で数値積分します。出てきた積分値をリチャードソン補外をして改善します。その結果は理屈上シンプソン法によって求めたものと一致します。 実行結果 まずもっとも細かい分割数での座…
半径１の１／４円中に入るか入らないかで円周率を求める古典的な例です。一千万回の試行をしているけれども、組み込み乱数がちゃんとしているかは不明ｗ 試行回数Nに対して、おおむね 1/sqrt(N) に比例して誤差が減るはずなので目安の数字も出力しています。…
OOP言語の常識に疎いので、コンストラクタ・デストラクタをどうするのが良いのかよく分かりません。デストラクタが二回呼ばれていますが、一回目はコンストラクタ内で仮生成したものが破棄されるところで出ているのかしらん？デストラクタはF2003規格のfinal…
素数をｎ個もとめる。ガウスの素数定理に基づいて必要なエラトステネスの篩の大きさを推測する。求めたおおむねｎ個の素数の配列から最初のｎ個を取りだす。オプションとして /assume:realloc_lhs が必要。求める素数の数が大きい場合にはオプション /heap-a…
配列生成子を使ったものは、配列はみ出しの場合分けをしなくて済むが、死ぬほど遅いｗ 素朴な方がいい。 module m_isort2 contains subroutine isort2( x ) real, intent(in out) :: x(:) real :: tmp integer :: i, k, n n = size(x) if (n <= 1) return do…
再帰を使ったので、すぐスタックあふれで死ぬｗ コンパイラオプション /assume:realloc_lhs が必要。 module m_isort implicit none contains recursive function isort(arr) result(res) real, intent(in) :: arr(:) real, allocatable :: res(:) if (size(…
昨日の奴に PUBLIC / PRIVATE のアクセス指定を掛けてみたけど、まだよく挙動がわかりません。 ソース・コード MODULE m_locate ! Win32 API call USE ifwinty USE kernel32 IMPLICIT NONE CONTAINS ! SUBROUTINE locate(ix, iy) ! move cursor to (ix, iy) …
以前はコンパイラのバグで動かなかったポインタを用いたプログラム例。 ALLOCATE( r1(1), SOURCE = t_rectangle(20.0, 10.0, 35.0, 10.0) ) こういう記法ができるので、初期化・代入ルーチンを自分で書く必要がありません。リリース・ノートを見ると、多態型…
オブジェクト指向型の TYPE 中で、その TYPE に作用する副プログラム名に総称名を用いることが可能です。 実行結果 ソースコード MODULE m_TYPE IMPLICIT NONE TYPE :: t_vector REAL, ALLOCATABLE :: v(:) CONTAINS PROCEDURE, PASS :: init PROCEDURE, PAS…
3.2.1 Features from Fortran 2003 FINAL subroutines GENERIC keyword for type-bound procedures A generic interface may have the same name as a derived type Bounds specification and bounds remapping list on a pointer assignment 3.7 Fortran 20…
FINAL 属性のサブルーチンは、TYPE 解放のときに自動実行されます。Fortran2003 では、オブジェクト指向の機能が取り入れられたので、スカラー型の変数も ALLOCATE で動的に確保できるようになりました。ここでは動的に確保されたスカラー構造体 tt 内部に、…
日経サイエンス誌の数学コラムをやっていたイアン・スチュアート氏の『数学の秘密の本棚』という本があったので暇つぶしに買ってみました。古いパズルから最近流行の話題まで色々載っており、また問題がわりあい易し目なのが多いので、理解できずにがっくり…
Fortran2003では、Array操作がさらに柔軟になっています。ALLOCATABLEな配列を返り値とする関数と、ALLOCATABLE配列の自動再割付が、色々と危険なんですが、便利です。これらの新機能により、APPENDやCONSのような配列長の変わる操作が簡単に実現できます。M…
Haskellのプログラムを参考にCollatz問題のプログラムを書いてみました。 Fortran2003のREALLOCATE機能を使うことで、配列生成子を使うことで、配列をリストの代用とし、再帰を使って自動で伸ばしてゆきます。Collatz問題とは何かは、詳しくはここを見てくだ…
Fortranで、関数合成 f(g(x)) = f.g (x) のような記述ができないかと思ったんですが難しい。 .op.型の自己定義演算子は、関数引数の場合は出来ないし。関数ポインタの配列も出来ないのでまだいい考えが浮かびません。 適当にテストのメモ帳。 Fortranではポ…
CDには16bitのPCMの音声データが記録されているわけですが、それは一種の相対値であって再生時の絶対的な音量には任意性があります。それはCD再生機のボリュームつまみで音量を変えられることからわかります。MPEGの音声圧縮のためには、心理音響解析が必要…
前回は、Fortran2003の新機能TYPE,EXTENDSを無理に使おうとして、冗長になったので作り直しをば。 ■実行結果 入力データは、ネット上に落ちていた2kHzの信号ファイルに、info chunkを手動で入れてみたもの。 サンプリング・レートが44.1kHzなので、2kHzの信…
MPEG1/Layer2,3では1フレームあたりのデータ数が変わって、FFTにかけるデータ数も変わるので一般形にしておきます。Fortran2003には本来、パラメトリックTYPEというのがあって、TYPEの中の要素のKINDや定数をTYPEインスタンス生成のときに指定できるんですが…
昔書いたFFTをOOP風にしてみました。とりあえず、べきは固定になってますが、イニシャライズのところでパラメータをとるようにすればいいような気もしなくも無いです。まぁ将来直すということで。普通教科書に載っているFFTは2のべき乗のものですが、同じよ…
チラシ裏の下書きメモ帖です。テストしてません。 ■ソース・コード 配列生成子のところで、[CHARACTER(LEN=nn):: ....]と明示的に指定することで、中のデータが自動的にその長さに調節されます。Fortran95までは、宣言長にそろえて空白を入れておかないと駄…
CRCについてはNumerical Recipesの第二版以降にある、Less numerical Algorithmsの章にある原理の説明を読むのがよかろうかと思います。もう細部はすっかり忘れてしまいました。覚えているのは、１．割り算の余りを比較する。２．余りが違ければ元の数もちが…
16ビット整数のwav fileのデータを、倍精度実数にして吐き出すことを考えます。今はmpeg音声Layer 1を想定して、ファイル先頭の１チャンネルあたり384個のデータを読みこみ書き出します。 ■実行結果 ■ソース・プログラム F2003の新機能ASSOCIATEをあまり必然…
Intel Fortran v.11.1は、まだ完全にFortran2003に対応していないのですが、ある程度オブジェクト指向がらみの機能が実装されているので、練習してみます。ここでは、WAVファイルのヘッダー情報を読むことを考えます。特に、WAVファイルのうちでもCDをリップ…
Numerical Recipes 2nd ed. §5.2の計算法を使って連分数展開の式を求めました。Fortran2003では、ALLOCATABLEな配列に、最初に確保した大きさとは異なる配列代入を行うと、自動的に新しい大きさにreallocateされます。これはバグの元になりそうな機能ですが…
14桁ずつ数値が決まっていきます。 ■実行結果 収束の様子と真値との差を表示しています。最後に4*ATAN(1)で真値を書き出しています。 ■ソースコード MODULE m_pi IMPLICIT NONE INTEGER, PARAMETER :: kq = KIND(1.0q0) CONTAINS REAL(kq) ELEMENTAL FUNCTIO…
参照：http://crd.lbl.gov/~dhbailey/dhbpapers/pi-quest.pdf ■実行結果 ■ソースコード MODULE m_pi IMPLICIT NONE INTEGER, PARAMETER :: kq = KIND(1.0q0) CONTAINS ELEMENTAL REAL(kq) FUNCTION pi_Borwein4(n) !P.Borwein : quartic INTEGER, INTENT(IN)…
