LU分解 pivot 付き Doolittle 法

Fortran2008

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

2011-11-29

メモ帳

Fortran2003

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

2011-11-27

LU分解 Doolittle 法

Fortran95

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

2011-11-26

ガウス消去法　枢軸交換付き

Fortran2003

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

2011-11-25

ガウス消去法　枢軸変換なし

Fortran90

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

2011-11-23

台形公式をリチャードソン補外でシンプソン法に。

Fortran2003

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

2011-11-22

Monte Carlo で円周率

Fortran2003

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

2011-11-21

STAMP CPUTIME

[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…

2011-11-17

■

uhoplot

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

2011-11-15

DOS窓経由ではなく window だけで。

Win32

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

2011-11-14

少し改良

uhoplot

重ね打ちを使っている絵。ソース・プログラム 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…