fortran66のブログ

fortran について書きます。

2011-01-01から1年間の記事一覧

DAVID GOLDBERG 『浮動小数点について知っておきたいこと。』

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

CMPLX の KIND 指定

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

DO LOOP 終了後の DO LOOP 変数の値

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

Householder法による三重対角化

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

クリントン国務長官の論文

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

講談社学術文庫版 宇野哲人『論語新釈』432頁 憲問第十四の二十

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

LU分解 Crout法 部分Pivot付き

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…

LU分解 Crout法

実行結果 ソース・プログラム 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…

LU分解 pivot 付き Doolittle 法

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

メモ帳

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

LU分解 Doolittle 法

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

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

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

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

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

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

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

Monte Carlo で円周率

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

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…

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

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

左クリックで翌月、右クリックで終了。1枚目に白紙が出ます。 ソース・プログラム 余り考えないで流用したので乱れている。 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風エラトステネスの篩

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

オンライン Fortran 実行環境

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

Webセミナー見た人に Intel Fortran のライセンス割引販売あり

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

エラトステネスの篩 ふたたびw

素数をn個もとめる。ガウスの素数定理に基づいて必要なエラトステネスの篩の大きさを推測する。求めたおおむねn個の素数の配列から最初のn個を取りだす。オプションとして /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…

素朴挿入ソート

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

挿入ソート

再帰を使ったので、すぐスタックあふれで死ぬw コンパイラオプション /assume:realloc_lhs が必要。 module m_isort implicit none contains recursive function isort(arr) result(res) real, intent(in) :: arr(:) real, allocatable :: res(:) if (size(…

米 citi bank 預金維持手数料免除額引き上げ

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

火野葦平「麦と兵隊」

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