LU 分解で連立一次方程式を解く（Fortran バージョン）
どうも，筆者です．前回に続き，Fortran で LU 分解を行うプログラムを作成してみた．LU 分解については，以前の記事を参照すること．
workspacememory.hatenablog.com
LU 分解のプログラム
需要あるのだろうか．これは，「luDecompSolver.f」というファイル名で保存している．
!----------------------------------------------------------------------! !---------------------- LU Decomposition ----------------------! !----------------------------------------------------------------------! !23456789|123456789|123456789|123456789|123456789|123456789|123456789|12 ! 1 2 3 4 5 6 7 subroutine luDecomp(n,pivot,matrix) implicit none ! === arguments === integer*4 :: n integer*4 :: pivot(n) double precision :: matrix(n,n) ! n: rank, scalar ! pivot: pivot vector, pivot(row) ! matrix: coefficient matrix, matrix(row,col) ! --- local variables --- integer*4 :: row,col,k,ip,iptmp double precision :: tmp,valmax,weight,mineps mineps=1e-10 do k=1,n pivot(k)=k end do do k=1,n-1 valmax=abs(matrix(k,k)) ip=k do row=k+1,n tmp=abs(matrix(row,k)) if (valmax.lt.tmp) then valmax=tmp ip=row end if end do ! the diagonal component is less than mineps if (valmax.lt.mineps) then write(*,100) 'diagonal component is less than',mineps 100 format(a,1X,G15.7) return else if (ip.ne.k) then do col=k,n tmp=matrix(ip,col) matrix(ip,col)=matrix(k,col) matrix(k,col)=tmp end do iptmp=pivot(ip) pivot(ip)=pivot(k) pivot(k)=iptmp do col=1,k-1 tmp=matrix(k,col) matrix(k,col)=matrix(ip,col) matrix(ip,col)=tmp end do end if do row=k+1,n weight=matrix(row,k)/matrix(k,k) matrix(row,k)=weight do col=k+1,n matrix(row,col) * =matrix(row,col)-weight*matrix(k,col) end do end do end do return end !----------------- last line of LU Decomposition ------------------! ! 1 2 3 4 5 6 7 !23456789|123456789|123456789|123456789|123456789|123456789|123456789|12 !----------------------------------------------------------------------! !---------------------- LU Solver ----------------------! !----------------------------------------------------------------------! !23456789|123456789|123456789|123456789|123456789|123456789|123456789|12 ! 1 2 3 4 5 6 7 subroutine luSolver(n,pivot,matrix,b) implicit none ! === arguments === integer*4 :: n integer*4 :: pivot(n) double precision :: matrix(n,n),b(n) ! n: rank, scalar ! pivot: pivot vector, pivot(row) ! matrix: coefficient matrix, matrix(row,col) ! b: right side vector, b(row) ! --- local variables --- integer*4 :: row,id double precision :: tmp,y(n) ! ===================== ! forkward substitution ! ===================== y(1)=b(pivot(1)) do row=2,n id=row-1 y(row)=b(pivot(row))-sum(matrix(row,1:id)*y(1:id)) end do ! ===================== ! backward substitution ! ===================== b(n)=y(n)/matrix(n,n) do row=n-1,1,-1 id=row+1 tmp=y(row)-sum(matrix(row,id:n)*b(id:n)) b(row)=tmp/matrix(row,row) end do end !----------------- last line of LU Solver ------------------! ! 1 2 3 4 5 6 7 !23456789|123456789|123456789|123456789|123456789|123456789|123456789|12
main routine
いつも通り，main routine のコードも載せておく．ちょっと調べて，動的配列を使ってみた．この動的配列は Fortran90 から導入されたらしい．エラー処理等を入れたので，goto 文とかがあって読みにくい．「main_lu.f」として保存した．
!----------------------------------------------------------------------! !--------------------------- main code ---------------------------! !----------------------------------------------------------------------! !23456789|123456789|123456789|123456789|123456789|123456789|123456789|12 ! 1 2 3 4 5 6 7 implicit none character :: fmtstr*32,argv*16,fname*16 integer*4 :: n,i,j integer*4,allocatable :: pivot(:) double precision,allocatable :: matrix(:,:),b(:) ! === check arguments === if (iargc().ne.1) then call getarg(0,argv) write(*,'("Usage: ", a, "[data file name]")') argv stop ! === read arguments === else call getarg(1,argv) read(argv,'(a16)') fname end if ! === read data from file === open(10,file=fname,status='old',err=1) goto 100 1 write(*,'("Error: ", a, "cannot open")') fname stop 100 continue read(10,*)n allocate(pivot(n),matrix(n,n),b(n)) do i=1,n read(10,*,end=200) (matrix(i,j),j=1,n),b(i) end do 200 continue close(unit=10) ! =========================== write(fmtstr,'("("I0"(G15.7,1X),G15.7)")') n do i=1,n write(*,fmtstr) (matrix(i,j),j=1,n),b(i) end do ! execute LU decomposition call luDecomp(n,pivot,matrix) ! solve call luSolver(n,pivot,matrix,b) write(*,'(/ "answer:")') write(*,'(G15.7)') (b(i),i=1,n) deallocate(pivot,matrix,b) stop end
読み込むファイルは，
[次元数（n: rank）] [係数行列 (1,1) 成分] [係数行列 (1,2) 成分] ... [係数行列 (1,n) 成分] [右辺ベクトル (1)] [係数行列 (2,1) 成分] [係数行列 (2,2) 成分] ... [係数行列 (2,n) 成分] [右辺ベクトル (2)] ... ... ... ... ... ... ... ... [係数行列 (n,1) 成分] [係数行列 (n,2) 成分] ... [係数行列 (n,n) 成分] [右辺ベクトル (n)]
という形式にする．前回の例だと，
3 2.0 1.0 -2.0 1.0 1.0 1.0 -1.0 4.0 1.0 -2.0 3.0 -1.0
となる．
Makefile
前回と同様に，Makefile も載せておく．
FC := gfortran PROG := solveLU SRCS := luDecompSolver.f main_lu.f OBJS := $(SRCS:%.f=%.o) all: $(PROG) $(PROG): $(OBJS) $(FC) -o $@ $^ %.o: %.f $(FC) -c -O2 $< .PHONY: clean clean: rm -f $(OBJS) $(PROG)
実行方法
Fortran のプログラムを上から順に，「luDecompSolver.f」，「main_lu.f」として保存したとする．また，Makefile を「Makefile」として保存したとする．その段階で，
[hoge@local]$ ls Makefile luDecompSolver.f main_lu.f
となっているはずである．まず，make コマンドを利用してコンパイルする．
[hoge@local]$ make gfortran -c -O2 luDecompSolver.f gfortran -c -O2 main_lu.f gfortran -o solveLU luDecompSolver.o main_lu.o
「solveLU」という実行ファイルが出来上がる．もし，make コマンドがなければ，上の実行結果（gfortran から始まる部分）をコピペして端末で実行すればよい．後は，データファイル input.dat を用意して，
[hoge@local]$ ./solveLU input.dat
とすれば，計算結果が出力される．また，実行ファイル，オブジェクトファイルの削除は，
[hoge@local]$ make clean rm -f luDecompSolver.o main_lu.o solveLU
とすればよい．
ガウスの消去法（Fortran バージョン）
お久しぶりです．筆者です．
アクセス解析のデータをみた
久しぶりに，ブログに戻ってきた．最初に，アクセス解析のデータをみた． すると，多くの人がガウスの消去法関連でアクセスしている事が分かった．その中でも，Fortran 言語でガウスの消去法を探している人が多い事に気付いた．
ガウスの消去法を Fortran 言語に書き直す
筆者は，ほとんど C 言語を利用するが，Fortran も少しは触れた事があるので，色々調べながら，以前書いた記事「ガウスの消去法のプログラム」を Fortran77 の固定形式で書き直す事にした． 何故，Fortran77 の固定形式かというと，筆者がそれしか読んだ事がないためである．この平成の時代に，Fortran77 の固定形式を勉強するとは思わなかったが，昔の方が書いたプログラムを使う以上，仕方ない事である．
workspacememory.hatenablog.com
書き直したもの
という事で，色々調べながら，C 言語のプログラムを Fortran で書き直した．読んだ事あるのは固定形式だが，書いた事はないので，コンパイルが通る程度にしか書けない．90/95 の表現もあるかもしれない．
!----------------------------------------------------------------------! !---------------------- Gaussian elimination ----------------------! !----------------------------------------------------------------------! !23456789|123456789|123456789|123456789|123456789|123456789|123456789|12 subroutine gaussSolver(n,matrix,b) implicit none ! === arguments === integer*4 :: n double precision :: matrix(n,n),b(n) ! n: rank, scalar ! matrix: coefficient matrix, matrix(row,col) ! b: right side vector, b(row) ! --- local variables --- integer*4 :: i,row,col,pivot double precision :: tmp,valmax,diag,scl,mineps mineps=1e-10 do i=1,n-1 ! execute partial pivoting pivot=i valmax=abs(matrix(i,i)) ! search maximum value in a column do row=i+1,n tmp=abs(matrix(row,i)) if (valmax.lt.tmp) then pivot=row valmax=tmp end if end do ! replace if (pivot.ne.i) then do col=1,n tmp=matrix(i,col) matrix(i,col)=matrix(pivot,col) matrix(pivot,col)=tmp end do tmp=b(i) b(i)=b(pivot) b(pivot)=tmp end if diag=matrix(i,i) if (abs(diag).lt.mineps) then return end if ! =================== ! forward elimination ! =================== do row=i+1,n scl=matrix(row,i)/diag do col=i,n matrix(row,col)=matrix(row,col)-matrix(i,col)*scl end do b(row)=b(row)-b(i)*scl end do end do ! ===================== ! backward substitution ! ===================== b(n)=b(n)/matrix(n,n) do row=n-1,1,-1 tmp=b(row)-sum(matrix(row,(row+1):n)*b((row+1):n)) b(row)=tmp/matrix(row,row) !------------------------------------------------------------ ! do row=n,1,-1 ! valsum=b(row) ! ! do col=row+1,n ! valsum=valsum-matrix(row,col)*b(col) ! end do ! b(row)=valsum/matrix(row,row) ! end do !------------------------------------------------------------ return end
main routine
ついでに，main routine の部分も以下に示しておく．
!----------------------------------------------------------------------! !--------------------------- main code ---------------------------! !----------------------------------------------------------------------! !23456789|123456789|123456789|123456789|123456789|123456789|123456789|12 implicit none character :: fmtstr*32 integer*4 :: n,i,j parameter(n=3) double precision :: matrix(n,n),b(n) matrix(1,1)=2.D0 matrix(1,2)=1.D0 matrix(1,3)=-2.D0 matrix(2,1)=1.D0 matrix(2,2)=1.D0 matrix(2,3)=-1.D0 matrix(3,1)=1.D0 matrix(3,2)=-2.D0 matrix(3,3)=3.D0 b(1)=1.D0 b(2)=4.D0 b(3)=-1.D0 write(fmtstr,'("("I0"(G15.7,1X),G15.7)")') n do i=1,n write(*,fmtstr) (matrix(i,j),j=1,n),b(i) end do ! output matrix and right side vector open(10,file='input.txt',status='replace') do i=1,n write(10,fmtstr) (matrix(i,j),j=1,n),b(i) matrix(i,:)=0.D0 b(i)=0.D0 end do close(unit=10) ! read data from file open(11,file='input.txt',status='old') do i=1,n read(11,*,end=100) (matrix(i,j),j=1,n),b(i) end do 100 continue close(unit=11) ! execute Gaussian elimination call gaussSolver(n,matrix,b) write(*,'(/ "answer:")') write(*,'(G15.7)') (b(i),i=1,n) stop end
いつも通り，英語が出来ないのはスルーして欲しい．無駄にファイル出力して，ファイル読み込みしているが，これは，単純に，Fortran で読み書きが出来るかどうかを確認しているだけである．たいした意味は無い．
Makefile
上のプログラムは自由に使ってもらってかまわないが，コンパイルに困る可能性もある．そのため，Makefile を作成した．Ubuntu 14.04 LTS でインストールできる gfortran で動作を確認している． Makefile を使った事がない人は，1 つのファイルに全て書いて，それをコンパイルしてもらえればよいと思う． また，この gfortran は，Fortran95 用のコンパイラなので，フリーダムな書き方をしていても文法のミスさえなければ通るみたいだ．以下に，Makefile を示す．
FC := gfortran PROG := solveGauss SRCS := gaussSolver.f prog_main.f OBJS := $(SRCS:%.f=%.o) all: $(PROG) $(PROG): $(OBJS) $(FC) -o $@ $^ %.o: %.f $(FC) -c -O2 $< .PHONY: clean clean: rm -f $(OBJS) $(PROG)
sub routine 「gaussSolver」を「gaussSolver.f」として，main routine 「main code」を「prog_main.f」として保存し，make とすると，実行ファイル「solveGauss」が出来上がる．後は，実行すると解が得られる． make clean とすると，.o ファイルと solveGauss が削除される．
[hoge@local]$ ls Makefile gaussSolver.f prog_main.f [hoge@local]$ make gfortran -c -O2 gaussSolver.f gfortran -c -O2 prog_main.f gfortran -o solveGauss gaussSolver.o prog_main.o [hoge@local]$ make clean rm -f gaussSolver.o prog_main.o solveGauss
【更新】レーベンバーグ・マーカート法と準ニュートン法（C 言語）
どうも，筆者です．
以前，同一のタイトルでレーベンバーグ・マーカート法と準ニュートン法の使い方を説明したが，あれは，Intel コンパイラがないと動作しないという欠点があった（また，3 次元のフィッティングが限界である）．
workspacememory.hatenablog.com
他に良いものが無いか探していたところ，GNUPLOT の fit コマンドが精度良くフィッティングを行ってくれるという事が分った．ただ，これは 4 変数までしか使用出来なかった． 仕方ないので，ソースコードを軽く読んで，グローバル変数等を使用しない C 言語プログラムに落とし込むことを考えた．また，準ニュートン法に関しては，Numerical Recipes in C（Second Edition，1992，小学生の頃だ）が，公開されていたので，そこにあるソースコードを自分好みに変更して使用するようにした．このため，ソースコードの殆どの部分にコメントが入っていない．
ここでは，main 関数だけ紹介しておく．
プログラムとして，何が使いやすいか分らないが，main 関数にユーザー定義関数（フィッティング関数）を定義しておき，最適化のときにこの関数を利用するようにプログラムしてある．なので，ほぼ全ての関数にユーザー関数を引数として投げる形になっている．作ったときはこれで良いと思ったが，今考えると修正が一切出来ない不便なものになってしまっている．
ここにプログラムをおいておく．パスワードは【gm4tj6dv】である．英語が全く出来ていないのは勘弁して欲しい．苦手なので，ある程度調べてはいるが，それでも（変数名とか関数名とか）色々なところに間違いがある．
https://drive.google.com/folderview?id=0BxLebOn2HFHoRVJYMTdGTzNaY3c&usp=sharing
【改訂版】 LaTeX に複数のソースコードを載せる
どうも，お久しぶりです．筆者です．
更新理由
以前，LaTeX の記事 を書いたが，プリアンブルの記述が汚かったので，もう少し綺麗にかけないものかと思い，色々調べた．
方法
方法という程しっかりしたものではないが，以前作成したものを少しいじればよい． 具体的には，言語を指定しないスタイルを定義（lstdefinestyle）し，後で，これを言語毎に呼び出せばよい．
比較的簡単ではあるが，以前ここまで出来なかったのは，単純に，複数の言語を LaTeX に取り込めるという段階で満足してしまっていたためである．お恥ずかしい．
コード
さて，いつも通り，コードを載せておく．今回は，アルファベット順に並び替え，色の定義等も行っておいたので，適宜変更すれば自分好みのものに変更出来る．
また，ソースコードの一番下に「includeCode」というコマンドも定義しておいた．これは，ファイルから呼び出すときに便利であると思う． 例えば，「./src/main.c」にある C コードを LaTeX に貼りたい場合は，
\includeCode[C]{./src/main.c}{main}
としてもらえればよい．最初の C は「CustomC」の C である．という事は，「./sample.sh」の shell script を LaTeX に貼りたい場合は，
\includeCode[Bash]{./sample.sh}{title}
とすれば良いのである．また，3 つ目の引数は，ソースコードの caption として表示される項目である．上記の場合，それぞれ「main」と「title」と表示されるはずである．
一言
今まで gist にコードを書いて，スクリプトを持ってきていたが，はてなブログでも，シンタックス・ハイライトに対応しているんだな．最近知ったことである．今後どうやってコードを書いていこうか，悩みどころだ． まぁ，今回みたいに，短いものはここに書いて，長いものは gist だろうなぁ．
Ubuntu 14.04 LTE の入力メソッド
どうも筆者です．
最近，仮想 OS に Ubuntu 14.04 LTE をインストールした． しかし，デフォルトの Unity では動作が重いため，Xubuntu を導入することにした．
Xubuntu の導入
Xubuntu 自体は，コマンド叩けば導入が出来た．
ところが，いつも通り日本語を入力しようとすると，直接入力と日本語入力が，高速に切り替わる状況となってしまった．すなわち，日本語を入力しようとしても，タイミングが悪いと直接入力になってしまい，直接入力に変更しようとしても，タイミングが悪いと日本語入力のままという状況になった．このままでは，使いづらいので，何とかする必要があった．
入力メソッドについて
まずは，入力メソッドについて調べた．すると，IBus や Fcitx，uim というものがあることが分かった．最近では，Fcitx が主に利用されているようである．
入力メソッドの確認
ここでは，自分の計算機が Fcitx になっているか調べるため，「設定」→「言語サポート」から，「キーボード入力に使う IM システム」の項目を調べ，fcitx となっていることを確認した． 次に，「設定」→「fcitx 設定」から，fcitx の設定を確認した．
入力メソッドの設定内容
入力メソッドには，「日本語キーボード」と「Mozc」が設定されており，全体の設定では，Zenkakuhankaku と Ctrl + Space で，入力メソッドのオンオフが出来るようになっていた． 試しに，Zenkakuhankaku のキーを ESC キーで無効にしたところ，問題はなくなったが，Ctrl + Space で入力切替を行わないといけない状況になってしまった．コレでは不便である． そこで，どこかに設定ファイルがあると思い，探してみると，以下の場所にあることが分かった．
入力メソッドの設定ファイルの変更
設定ファイルの場所が分かったので，これを編集して設定を行う．まずは，ESC キーで切り替えが行えるようにした．
これで，荒ぶることなく切り替えが行えるようになったが，Vim を使用する筆者にとって，これでは少し不便であった．設定ファイルをもう少し見ていくと，インプットメソッドをオンにする時と，オフにする時のキーが別々で設定できることが分かった．
そこで，インプットメソッドをオンにするときは，半角/全角キーを，オフにするときは ESC キーを使うように設定を行った．
これで，荒ぶることなく，Vim に支障なく日本語と直接入力の切り替えが行えるようになった．
