元々Sunのサイトにあったんですが、Oracleに買収されてリンクが移ったのでメモ。What Every Computer Scientist Should Know About Floating-Point Arithmetic. PDF直リン１ PDF直リン２浮動小数点演算について PDF直リン 付録D

組み込み関数 CMPLX() は REAL(4) がデフォールト値なので、引数が倍精度でも第二引数で KIND を倍精度に指定しないと、単精度で変換されます。 実行結果 ソース・プログラム program Complex implicit none complex (kind = 8) :: c real (kind = 4) :: a r…

FORTRAN66 では DO LOOP 終了後の DO LOOP 変数の値は不定ですが、FORTRAN77 以降は終了時の値です。

昔Householder法とQR法のサブルーチンを外人からもらいましたが、EISPACKの tred と tql で F66 で書かれていました。のちに Numerical Recipes の F77 版に置き換えたけれども、基本的に LAPACK のものだと思います。微妙に変えてあって、tql は旧版で収束…

クリントン（雌の方）がForeign Policy誌10月11日号に発表した論文“America's Pacific Century”がネットで話題だというので暇つぶしにななめ読みしたが、日米安保破棄も何も書いてないｗ このおっさんの記事『日米安保破棄を真剣に検討し始めた米国』はなん…

子言衛霊公之無道也、康子曰、夫如是、奚而不喪。孔子曰、仲叔圉治賓客、祝鴕治宗廟、王孫賈治軍旅。夫如是、奚其喪。[通釈]孔子が衛の霊公の無道なことを語られたので、季康子が問うて曰うには、「お話のように無道なのに、どうして位を失わないのですか。…

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…

配列の割り当てに Fortran2008 の source オプションを使ってみました。 pivotを使うと精度が上がります。 実行結果 ソース・プログラム module m_mat implicit none integer, parameter :: kd = kind(0.0d0) contains subroutine lu2(a, b, x) real(kd), in…

プログラム中で１から１０までの整数列が欲しくなることがありますが[(i, i = 1, 10)]と配列生成子でつくるのも、捨て変数 i を宣言しないといけないので面倒です。Haskell のように 1..10 のように生成できるとうれしいのですが、Fortranならユーザー定義演…

基本的にガウス消去法と同じ方式になっています。 実行結果 LU分解では精度がやや悪くなるようです。 Eliza’s father=Alfred P.DoolittleStation master=AlfeeYoung Freddy=Sherlock Holmes ソース・プログラム 簡単に書けるのがうれしい。ほとんど実質３〜…

もともとの行列が上書きされてしまうので、中間サブルーチンを置いて作業行列などを確保します。枢軸（pivot)用のベクトルもここで用意します。 実行結果 500*500の行列の線型方程式を20回解いて残差の二乗平均を出力。おおむね数値誤差程度のまで求められて…

線型方程式 Ax=b をガウス消去法で解くプログラム。Fortran90 だと式をそのままプログラムに直せばいいので楽。でもPivotting を入れるとやや見苦しくなる。並列計算機では LU分解で Crout法より Gauss 消去法の方が有利になることもあるとか。ベクトル機全…

ローレンツ型の関数を‐1から1の範囲で、いくつかの分割数に対して台形公式で数値積分します。出てきた積分値をリチャードソン補外をして改善します。その結果は理屈上シンプソン法によって求めたものと一致します。 実行結果 まずもっとも細かい分割数での座…

半径１の１／４円中に入るか入らないかで円周率を求める古典的な例です。一千万回の試行をしているけれども、組み込み乱数がちゃんとしているかは不明ｗ 試行回数Nに対して、おおむね 1/sqrt(N) に比例して誤差が減るはずなので目安の数字も出力しています。…

[Fortran95]メモ帳 module m_stamp implicit none contains subroutine stamp(text) character (len = *), intent(in), optional :: text character (len = 40) :: buff character (len = 10) :: date, time real, save :: time0, time1 logical, save :: qf…

途中経過メモ 過去のルーチン流用で乱れ気味。 コンソールプログラムから呼び出せて、コンソールで入出力しながら簡単な二次元グラフを描いてゆける最小限度の簡潔なサブルーチン集を作りたいです。 ソース・プログラム Fortran2008 ではベッセル関数や誤差…

左クリックで翌月、右クリックで終了。１枚目に白紙が出ます。 ソース・プログラム 余り考えないで流用したので乱れている。 module m_sub use ifwina use kernel32 type :: t_wnd integer (HANDLE) :: hWnd integer (HANDLE) :: hDC end type type (t_wnd),…

重ね打ちを使っている絵。 ソース・プログラム module uho_win use ifwina use ifwinty use ifmt, only : RTL_CRITICAL_SECTION implicit none type :: t_wnd integer (HANDLE) :: hWnd integer (HANDLE) :: hDC integer (LPINT) :: hThread integer :: id i…

ラインプリンタの重ね打ちを使っているので、Windowsのグラフィックでラインプリンタ動作を再現するように作ってみた試し。もう少し完成度をあげたい。 ラインプリンタ出力＋α スヌーピーカレンダーデータは別プログラム。過去日記参照。 テスト 未完成↓ MOD…

OOP言語の常識に疎いので、コンストラクタ・デストラクタをどうするのが良いのかよく分かりません。デストラクタが二回呼ばれていますが、一回目はコンストラクタ内で仮生成したものが破棄されるところで出ているのかしらん？デストラクタはF2003規格のfinal…

sourcelair というオンライン実行環境が開発中のようで、Fortranも使えるようです。登録すると、ファイルの保存などもできます。今のところ設定ミスか Module が使えないようです。H23-11-08 メールしてみたところ、今のところ全言語で Module は使えないそ…

日本からもイケるのか分かりませんが、Stieve Lionel の webinar を見ると、Intel Fortran の更新ライセンスが2割引きで買えるようです。セミナーはここから。 http://t.co/JJJ5x8eJ 滑舌が悪いｗ TechXtendがパキスタンに聞こえた私は糞耳ｗ

素数をｎ個もとめる。ガウスの素数定理に基づいて必要なエラトステネスの篩の大きさを推測する。求めたおおむねｎ個の素数の配列から最初のｎ個を取りだす。オプションとして /assume:realloc_lhs が必要。求める素数の数が大きい場合にはオプション /heap-a…

ここまで来ると更に簡単に。尻からではなく、頭から sort してゆく。 module m_isort2 implicit none contains subroutine isort3y( x ) real, intent(in out) :: x(:) integer :: j, k do k = 2, size(x) do j = 1, k - 1 if ( x(k) < x(j) ) exit end do x…

少し改良。配列のサイクリックシフトで。 module m_isort2 implicit none contains subroutine isort2y( x ) real, intent(in out) :: x(:) integer :: j, k, n n = size(x) if (n <= 1) return do k = n - 1, 1, -1 do j = k + 1, n if ( x(k) < x(j) ) exi…

配列生成子を使ったものは、配列はみ出しの場合分けをしなくて済むが、死ぬほど遅いｗ 素朴な方がいい。 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(…

今度は $15,000 に引き上げ！リーマンショックの時も $6,000 くらいに引き上げて学生とか貧民からド顰蹙くらってたのに。 大貧民・大富豪を地で行く展開ﾜﾛﾀ

ベストセラーになっただけあって読みやすく面白い。 華々しいところしか報じないジャーナリストへの不満が軍人側から何度か表明されるのも良い。そのためか縁の下の力持ちとして地味な仕事をしている兵隊さんの姿も描いている。また腐敗しきった支那人の有様…