Fortran90で構造解析をしよう！

last update on Dec. 19, 2007

はじめに

このページの内容は個人的な備忘録として作ったものですが，ひょっとしたら誰かの参考になるかもしれないと思って公開しているものです。なるべく間違いはないように気をつけてはいますが，（言い訳ですが）講義の資料や教科書用に作ったわけではないので間違いや勘違い（特に用語の使い方など）を含んだままの可能性が高いものです。ここを参考にしたせいで損害が生じたり恥をかいたりしても責任をもてませんので悪しからずご了承下さい。

目指すのは骨組の地震応答解析とか，有限要素法を用いた部材の非線形解析とか，マトリックス法による構造解析を行うことです。構造解析そのものについては別のところで勉強するものとして（英語ですがW. McGuireらの「MATRIX STRUCTURAL ANALYSIS」（WILEY）なんかわかりやすいと思います。），その環境を構築する方法を書いていきます。

環境構築の方針としては，次のように設定することにしました。

なるべくただで（パソコンはしょうがないけど高いソフトを買うなどは論外。パソコンもできれば今は使っていないもの再利用）
なるべく汎用性のある（特殊な環境ではなく，PCから大型計算機まで対応できそうな）

要するに研究室で誰も使わなくなったような古めのPCに環境を構築して，多少時間がかかってもいいからひたすら計算をしているようなマシンを作ってしまおうという感じです。解析してみていい結果がでそうなあたりがついたら，それをもとにドカンと研究費の申請でもしましょう。

このページの先頭

last update on Dec. 20, 2007

なぜFortran90？

さて，構造解析をするとしてプログラミング言語は何がよいか・・・。マトリックス法による構造解析なのでマトリックス（配列）が扱いやすいものがよい，となるとやはりMATLABです。極端にいうと連立一次方程式の解法とか，逆行列や固有値を求めたりとか，そのための数値計算法やそれをプログラミング化するための知識がなくても用意された関数を利用するだけで解析できてしまいます。マトリックスを用いた構造解析について勉強している段階や，解析のアルゴリズム（フロー）を検討している段階では余計なことに煩わされなくてよいでしょう。

しかし，MATLABはなかなか高価なソフトです。アメリカの大学なんかでは基本的にサイトライセンスで大学内では自由に使えたりするそうですが，日本ではそんなこともなく使いたかったらライセンスを購入しなければいけません。我々の研究室では使える環境になっていますが，少ないライセンスを長時間占有するわけにもいきません。また，基本的にインタプリタ式のプログラミング言語のため，実行するためにはMATLABの（もしくは互換の）環境が必要で，どこへでも持っていって実行するわけにもいきません。

MATLAB互換のフリーな環境としては，OctaveやScilab，FreeMATがあります。どれも使い込んでいないので詳しい論評は避けますが，いずれも非線形な地震応答解析などを本格的にやろうとすると遅くてやる気がなくなります。

それと，MATLABはインタプリタ式のプログラミング言語のため，やはりFortranやC言語などのコンパイラで実行体を作るものに比べると上記の互換環境ではなくMATLAB上で実行しても実行速度は落ちるような気がします。（MATLABのソースからC言語のソースを作成してCコンパイラでコンパイルするMATLAB Compilerというものや，MATLABからFortranやCのサブルーチンを呼ぶ方法もあります。）

ということで，ここから先は基本的にFortranでいくことにします。なぜCやC++ではなくFortranか？それは私が学生時代からFortranしか使ったことがないからです(笑)。両方使ってみて数値計算におけるFortranの優越性を確信したからではありません。ただし，Fortran90（もしくはFortran95）ではFortran77に比べると飛躍的に扱いやすくなっており，基本はFortran90にすることにします。

ただそれしか知らないからといってしまっては身も蓋もないので，数値計算を行う上でFortranが優れていると一般的に言われている点も含めて，Fortranから抜け出してCやC++にいかなかった理由を列挙しておきます。

配列の取り扱いが強力である。配列の添字に負も使えるし，Fortran90になってからは配列演算が強化されています。（配列の部分配列や豊富な配列関数が使えるようになりました。）
豊富な科学技術計算ライブラリがある。線形代数の計算を行うルーチンをまとめたLAPACKや高速フーリエ変換ライブラリのFFTWなど，比較的簡単にMATLABの関数の代わりが見つかります。（MATLABのマトリックスの関数の多くは内部でLAPACKのルーチンを用いているそうです。）
先人の残した多くの資産はFortranで書かれていることが多い。研究室の先輩方が作った多くのプログラムもFortranです。
Fortran90では配列の動的割り当てや構造型データの新設，DOループの構文化(文番号が不要になった)，自由形式コーディングの導入などFortran77までではプログラミング上不便だった点が大幅に改善された。
相変わらずオブジェクト指向プログラミングには不便だが，そもそも構造解析のプログラム自体あまりオブジェクト指向プログラミングに向いていない気がする。

手元に置いている参考書はL. Nyhoff, S. Leestma「入門Fortran90」（ピアソン・エデュケーション）ですが（この本の世間での評判は知りません），その他にFortran90を使ったプログラミングには以下のサイトをよく参考にしています。

Ｆｏｒｔｒａｎ９０プログラミング（京都大学冨田博之先生）
Fortran 90 骨の髄まで（いぐちさん"縦々横々～じゅーじゅーおーおー"内）
Fortran 標準コーディングルール（気象庁 ※公式なものではないと思う。）

このページの先頭

last update on Dec. 20, 2007

Linux(CentOS 5)のインストール

計算用のPCをFortranを使えるようにセットアップします。

まずはOSの選定です。フリーなFortranコンパイラが使える環境で，安定していて，コストがかからず，PCのスペックが低めでも動いて，ということでLinuxにしました。最近のPCはWindowsがプリインストールされた状態で購入しているためOSをWindowsにしても新たにコストがかかるわけではないですが，なかなかフリーのFortran90に対応したコンパイラがありません。

Windowsでも個人利用であれば，登録は必要ですがSilverfrost FTN95などが無料で使えるようですが，ここはこの先の数値計算のライブラリの導入などもUnix系の方が選択肢も情報も豊富だろうということでLinuxにします。LinuxではGNUのgccに含まれるgfortran（gccのver.4からg77に代わって標準のFortranコンパイラ）や非商用利用ならIntelのIntel FortranやMath Kernel Library（数値解析ライブラリ）が利用できます。また，Windows上にUnix系の環境を構築する手段としてCygwinやcoLinuxもありがますが，計算専用のマシンなのでここは素直にOSをUnix系にします。

OSはLinuxに決めたといっても，まだ星の数ほどあるディストリビューションの中からどれかを選ばなければいけません。Linuxについてたいした知識があるわけでもなくこだわりもないのですが，今回はCentOSにします。CentOSは商用パッケージであるRed Hat Enterprise Linuxと互換することを目指したフリーのLinuxディストリビューション（いわゆるRHELクローン）なので，安定性，情報量などの面でよいかな，と。最新バージョンは 5.1 （'07年12月現在）です。

いよいよインストールです。参考にしたのは下記サイトです。ただ，いちいちCDを焼くのも面倒なのでインストールはネットワークインストールにしました。

以下に設定の備忘録も兼ねてインストール手順を示しておきます。

配布ミラーのひとつ理研のFTPサイトからboot.isoをダウンロードしてCDに焼きます。
作成した起動用CDでPCを起動（BIOSでCDからブートできるようにするのを忘れない）します。
起動したらEnterキーを押してGUIによるインストールを始めます。
しばらくはCUI（マウスは使えない）だし，英語だけどどんどん進みます。
- choose a language: Japanese
- keyboard tyle: jp106
- installation method: FTP
インストール方法にFTPを選ぶとネットワークの設定が始まります。研究室ではNATを利用しておりこの計算用PCを接続するプライベートなネットワーク内ではDHCPも使えますが，IPアドレスが固定されていないと後々何かと不便なのでここでネットワークの設定もします。面倒ならここはDHCPにしておいて，インストール後ネットワークの設定をしても問題ありません。
- Enable IPv4 support: チェック
- Manual configuration: チェック（DHCPのチェックがはずれる）
- Enable IPv6 support: チェックしない
- IPv4 address: 192.168.0.100 / 255.255.255.0（/の前がマシンのIPアドレス，後ろがサブネットマスク：環境にあわせて）
- Gateway: 192.168.0.1（デフォルトゲートウェイ：環境にあわせて）
- Name Server: 192.168.0.1（1番目のネームサーバ：環境にあわせて）
続いてインストールに利用するFTPの設定です。日本国内だと理研のミラーが速くてよいでしょう。
- FTP site name: ftp.riken.jp
- CentOS directory: /Linux/centos/5/os/i386
CentOSのGUIインストールがいよいよ始まります。ここからはマウスが利用できますし，説明も日本語です。まずはディスクパーティションの設定方法を選びます。
- Disk Druidを使用して手動パーティション設定を選ぶ
パーティションを設定していきます。LVMはよくわからないので今回はパス。パーティションは4つに分けてそれぞれは次の通りにしました。/bootはわざわざ分けなくてもいいかもしれません。/homeを独立した基本領域（Primary Partition）にしておくと後日システムを再構築する際などにデータがある/homeをそのまま引き継げます。
1. /boot(128MB)
2. /(10GB)　※ このメモの通りインストールした場合の使用量は5GB強
3. swap(1024MB)
4. /home(残り)
iからiiiはまとめて拡張領域(Extended Partition)内に論理領域（Logical Partition）として作るのもありかと。システム再構築ではどうせまとめてフォーマットするんだし。
ブートローダーの設定になります。デフォルトのままで問題ありません。何かの事情でWindows XPとデュアルブートのマシンを構築する場合もGRUBが勝手にWindows XPがインストールされたパーティションを探してきてブートセレクタに取り込んでくれるはずです。
それでもWindowsとのデュアルブートにおいてブートセレクタにWindowsのNTLDRを用いたい場合は以下の手順で設定します。
1. 「高度なブートローダーのオプション」をチェックしてオプションの設定に進む。
2. インストール先はMBRではなくブートパーティション(/boot)の最初のセクタを選択する。
3. CentOSのインストール終了後PCを再起動するとWindows XPが起動する。（CentOSは起動できない）
4. BootpartをダウンロードしC:\に解凍する。
5. コマンドプロンプトでC:\に移動して'bootpart'と実行する。
6. 表示されたパーティション一覧から /boot があるパーティションを選んで
7. 'bootpart #(#:該当するパーティションの番号) c:\bootsect.lnx "CentOS 5" 'と実行する。これでブートセクタのファイル化，boot.iniへの登録ができる。
8. PCを再起動するとOSの選択画面が出るようになります。この画面の表示時間など詳細を変更したい場合はWindows上でboot.iniを直接編集するか，システムのプロパティ-詳細設定-起動と回復で設定できます。
続いてネットワークの設定です。5.で設定していればIPアドレスなどはすでに設定されています。
- ホスト名：drunk00.drunks（環境に応じて適宜）
- 2番目のDNS：（必要に応じて）
タイムゾーンの設定はデフォルト（アジア／東京）でよいでしょう。
rootのパスワードを設定します。
インストールするソフトウェアのパッケージ等を選択をします。（必要なソフトウェアは後でネットワークからインストールすることもできます。）
- 初期状態でインストールパッケージ: Desktop-Gnomeのみチェックをいれる
- 追加リポジトリ: CentOS Extrasにチェックをいれる
- 今すぐカスタマイズを選んで次に進む
インストールするソフトウェアのカスタマイズです。当面必要になりそうないくつかを追加と不要なものの削除をします。
- アプリケーション / Emacs: 追加 → デフォルト
- アプリケーション / グラフィカルインターネット: 追加 → firefoxのみチェック
- アプリケーション / 上記以外は不要
- 開発 / レガシーなソフトウェア開発: 追加 → デフォルト
- 開発 / 開発ツール: 追加 → デフォルト
- サーバー / Windowsファイルサーバー: 追加 → デフォルト
- サーバー / サーバー設定ツール: samba，securitylevel，servicesあたりをチェック
- サーバー / ネットワークサーバー，レガシーなネットワークサーバー: デフォルトのまま
- サーバー / 上記以外は不要
- ベースシステム / ダイヤルアップネットワークサポート: 不要
以上で準備は終了です。必要なパッケージをFTPで自動でダウンロードしながらインストールが進んでいきます。インストールが終わったらPCを再起動（起動用CDを取り出すのを忘れずに）します。再起動したらそのまま基本設定を行っていきます。
ファイアウォールの設定は，ファイアウォールに守られたプライベートなネットワークにのみ接続するのでとりあえず無効にします。
SELinuxも同様に安全なプライベートネットワークにのみ接続するのでとりあえず無効にします。有効にすればセキュリティ的に強固になるりますが，設定が複雑になるのでやめておきます。
日付と時刻の設定をします。ntp（ネットワークプロトコル）を有効にしておきます。
root以外のユーザーを作成します。
- ユーザー名: sake（管理用のユーザー）
- フルネーム: （適宜）
- パスワード: （適宜）
サウンドカードの設定は特に何もせず次へ。
追加のCDも，特にCDから追加するものはないのでそのまま終了します。
PCを再起動して全て終了です。

このページの先頭

last update on Dec. 20, 2007

CentOS 5の設定

まずはネットワークに接続して運用する上で最低限必要な設定と，運用していく上で設定しておくと便利な設定をしていきます。なお，以下の設定方針はあまり凝ったことはせず小さい労力でとりあえず動くマシンを作ることなので，セキュリティに関しては特に大きな注意は払っていません。

以下で，ターミナルからの入力のうち，"#"に続いて記述されたものは root 権限で実行するもの，"$"に続いて記述されたものは一般ユーザーの権限で実行するものです。

sshの設定

デフォルトでは ssh でいきなり root としてログインできてしまいますが，さすがにこれはセキュリティ上好ましくないのでこれを無効にします。root になって /etc/ssh/sshd_config を開き

#PermitRootLogin no
↓
PermitRootLogin no

sshd を再起動しておきます。

# /sbin/service sshd restart

rootになれるユーザーを管理者のみにする

さらに su コマンドで root になれるユーザーを管理用ユーザー（sake）のみにしておきます。まず，管理者ユーザーをwheelグループに追加します。

# /usr/sbin/usermod -G wheel sake

/etc/login.defs 内に以下の一行を追加します。

SU_WHEEL_ONLY yes

/etc/pam.d/suを編集します。

#auth required pam_wheel.so use_uid
↓
auth required pam_wheel.so use_uid

パッケージ管理ソフト yum のリポジトリの追加

パッケージ管理，アップデートソフトウェア yum で利用するリポジトリを追加します。

CentOSPlusの追加

/etc/yum.repos.d/CentOS-Base.repo を次のように編集します。

[centosplus]
enabled=0
↓
enabled=1

Fedora EPEL リポジトリを設定

RHEL および CentOS を含むそれに互換なディストリビューションを補うためのFedora EPEL リポジトリを，次の通り rpm をダウンロードしてインストールすることで設定します。

# wget http://download.fedora.redhat.com/pub/epel/5/i386/epel-release-5-2.noarch.rpm
# rpm -ivh epel-release-5-2.noarch.rpm

fedora-extrasを使えるように設定する

RHEL5やCentOS5用のリポジトリに登録されていないソフトウェアはでも，FC6用のリポジトリには登録されていてそこからインストールすると動作する場合もあります。ただし，必ず動くという保証はないし，ここからインストールすると整合性が失われる可能性もありますが，有用なものがある場合もあるので設定をしておきます。必ずデフォルトをオフにしておきます。

/etc/yum.repos.d/fedora.repoというファイルを作って次の内容を記述します。

[fedora-extras]
name=Fedora Extras 6 - $basearch
#baseurl=http://download.fedora.redhat.com/pub/fedora/linux/extras/6/$basearch/
mirrorlist=http://mirrors.fedoraproject.org/mirrorlist?repo=extras-6&arch=$basearch
#mirrorlist=file:///etc/yum.repos.d/local-extras
enabled=0
gpgkey=http://ftp.riken.jp/Linux/fedora/extras/RPM-GPG-KEY-Fedora-Extras
gpgcheck=1

なお，このリポジトリを利用するには yum を実行時に"--enablerepo=fedora-extras"を指定します。

yumの設定とソフトウェアのアップデート

yumの接続先を自動で最も速い場所にします。

# yum install yum-fastestmirror

サービスの設定でyumをオンにしてソフトウェアの自動アップデートが夜間に行われるようにしておきます。

# /sbin/service yum-updatesd start
# /sbin/chkconfig yum-updatesd on

ソフトウェアをアップデートします。（初回はかなり時間がかかります。）

# yum -y update

Windowsとのファイル共有（Samba）の設定

計算用PCはLinuxにしたとはいっても，私の場合，普段使うメインマシンはWindows XPです。Windowsとのファイルのやりとりが便利なように Samba を動かすことにします。

CentOSのインストール時にインストールするパッケージに"Windowsファイルサーバー"を選んであれば既に Samba はインストール済みです。後からインストールするには yum を使って次のようにします。

# yum install samba

Samba の設定ファイルを編集します。あまり凝った運用はしない（できない）のでデフォルトの設定ファイルに最小限の変更を加えるのみです。/etc/samba/smb.conf を開いて該当箇所を修正します。他はそのままにしておきます。

[global]（全般の設定）
workgroup = drunks（Windowsのワークグループ名を指定）
※ドメイン管理はしないのでドメイン名と同じにする必要はありません。
load printers = no（Sambaでプリンタを共有しない場合）

[homes]（ホームディレクトリの設定）
comment = Home Directories
browseable = no
writable = yes
※ [homes]以外の共有の設定はすべて";"または"#"でコメントアウト

Samba を起動します。

# /sbin/service smb start

マシンの起動時にSambaが自動で立ち上がるように設定します。

# /sbin/chkconfig smb on

手動インストール用のシンボリックリンクの設定

正式なLinuxの運用ルールは知りませんが，システムが管理しているものとは別に手でソフトをインストール場合，/usr/local の下にインストール場合と /opt の下にインストールする場合があるようです。2箇所にわかれていると管理する上で何かと不便なのでどちらかに統一したいところですが，なるべくデフォルトの設定のままインストールを行いたいような場合は，次のようにシンボリックリンクを作っておくと便利です。（すでに /opt 以下に何かインストールされた状態ではしないように。）

# ln -s /usr/local /opt

一般ユーザーの登録

セットアップしたPCを計算に使う前に一般のユーザーを登録します。1人で使うPCであっても，管理用ユーザー（今回はsake）と一般の作業を行うユーザーは分けておいたほうがよいでしょう。

CUIでユーザーを追加するには次のようにします。指定しているオプション"-m"はこのユーザーのホームディレクトリを作成するためのものです。ユーザーIDを指定したい場合は'-u [uid]'で指定します。

# useradd -m wine（wineという名のユーザーを作成）
# passwd wine（wineのパスワード設定）
Changing password for user wine.
New password:
Retype new password:
passwd: all authentication tokens updated successfully.

続いてユーザーwine を Samba のユーザーに登録します。

# pdbedit -a wine（wineをSambaのユーザーに追加）
new password:
retype new password:

なお，上のSambaの設定ではシステムとSambaのパスワードの同期をとるようになっていないので，システムのパスワードは 'passwd' で，Sambaのパスワードは 'smbpasswd' でそれぞれ変更する必要があります。

このページの先頭

last update on Dec. 21, 2007

Fortranコンパイラのインストール

まずはFortran90に対応したコンパイラを用意しなければいけません。幸いなことに GNU gcc の gfortran がFortran95をサポートしていますので基本的にこれを使うことにします。

また，プログラムの開発時にバグをつぶしていく段階では複数のコンパイラを使ってみて問題ないかのチェックもしたほうがいいでしょう。そのために gfortran に加えて他のコンパイラの場合もあわせてインストールしておくとかもしれません。

gcc-gfortranのインストール

gfortranはCentOSのインストール時にインストールするパッケージに"開発ツール"を選択していれば既にインストールされているはずです。後から gfortran をインストールするには yum を利用して以下のように行います。（Fortran以外のコンパイラなどもインストールされます。）

# yum install gcc* compat-gcc* compat-glibc* compat-lib*

G95のインストール

gcc ではありませんがスタンドアローンのFortran95 コンパイラとしてg95があります。（これは元はといえばgccのgfortranから分岐したものだそうです。）

g95のインストールは次のように行います。といってもg95のダウンロードのページからLinux x86用のバイナリパッケージをダウンロードして展開するだけです。/usr/local の下にインストール（下はソースを/usr/local/srcにダウンロードした場合）することにします。

# cd /usr/local
# tar zxvf src/g95-x86-linux.tgz

使いやすいように /usr/local/bin （ここなら一般ユーザーもデフォルトでPATHが通っている）にシンボリックリンクを作っておきます。

# ln -s /usr/local/g95-install/bin/*g95* /usr/local/bin/g95

Intel Fortran Compilerのインストール

また，詳細は省略しますが（いずれ追記するかもしれませんが），非商用であればLinux版の Intel Fortran Compiler を利用することもできます。Intel Software Evaluation Center内にあるNon-Commercial Software Developmentからダウンロードしてインストールします。メールアドレスなどの登録が必要です。

last update on Dec. 21, 2007

数値計算ツールのインストール

もちろん，Fortranで構造解析のプログラムを作るにあたって数値計算上必要なサブルーチンを自分で作ってもよいですが，精度が高く，しかも速度が出るものを作るのはなかなかたいへんです。もちろん，解法の原理を知るために一度ぐらいは自分でプログラムを作るのも大事かもしれませんが，それよりは定評があるライブラリを使うほうが少ない労力で確実な結果が得られます。なので，今回もなるべく既存の無償で公開されているライブラリを使うことを考えます。

LAPACKのインストール

LAPACK(Linear Algebra PACKage)は netlib で公開されている有名な線形代数ライブラリの1つです。連立1次方程式を解いたり，固有値問題を解いたり，LU分解やコレスキー分解するためのルーチンが用意されています。これを使えるようにインストールすることにします。

yum を使って以下のようにインストールするだけです。

# yum install blas lapack

実はLAPACKはベクトル演算や行列演算など計算における核の部分で BLAS (Basic Linear Algebra Subprograms) ライブラリを呼び出します。このBLASも用意しなければいけないわけですが，yum を使ったインストールの場合，もし blas を指定し忘れても自動で依存関係を調べ BLAS がインストールされていなければ一緒にインストールしてくれます。

コンパイル時にライブラリを指定しやすいようにシンボリックリンクを作成しておきます。

# ln -s /usr/lib/libblas.so.3 /usr/lib/libblas.so
# ln -s /usr/lib/liblapack.so.3 /usr/lib/liblapack.so

今回のように自分で管理しているCentOSのPCにインストールする場合は上記のようにするだけですみますが，使用しているLinuxのディストリビューションのパッケージ管理のリストにLAPACKが含まれていない場合や，管理者でない環境でLAPACKを使用する場合ソースからコンパイルする必要があります。また，ライブラリを静的にリンクしたい場合も同様です。

netlibのLAPACKのページからlapack.tgzを適当なディレクトリ（管理者で全員が利用できるようにインストールするなら /usr/loca/src がいいと思います）にダウンロードして展開し，できた lapack-3.1.1 に移動します。

$ cd /usr/local/src
$ tar zxvf lapack.tgz
$ cd lapack-3.1.1

make が参照する make.inc の雛形をコピーしてきます。

$ cp INSTALL/make.inc.LINUX ./make.inc

Frotranのコンパイラにgfortranを使う場合は make.inc を次のように編集します。

FORTRAN  = gfortran
OPTS     = -funroll-all-loops -O3
DRVOPTS  = $(OPTS)
NOOPT    =
LOADER   = gfortran
LOADOPTS =

ifort（Intel Fortran Compiler）の場合は make.inc はこんな感じです。

FORTRAN  = ifort
OPTS     = -O3
DRVOPTS  = $(OPTS)
NOOPT    =
LOADER   = ifort
LOADOPTS =

次の順にライブラリを作成します。tmglibは時間計測用のライブラリのようなので必要なければ作成しなくてもよいでしょう。

$ make blaslib
$ make lapacklib
$ make tmglib

うまくいけば blas_LINUX.a，lapack_LINUX.a，tmg_LINUX.a の3つのファイル（ライブラリ）が完成しています。これを任意のディレクトリに置いておきます。管理者なら /usr/local/lib，自分用のライブラリ専用ディレクトリを作るなら ~/lib あたりがよいでしょう。

これでライブラリは完成ですが，このままではコンパイルのたびにライブラリを指定するのにファイル名を全て入力しなければなりません。GNUの開発ツールでは，ライブラリには lib で始まるファイル名をつけることになっているとのことなので，ファイル名を変えておきます（シンボリックリンクを作るのでもよいでしょう）。ライブラリを置いてあるディレクトリに移動して次のように実行します。

$ mv blas_LINUX.a libblas.a
$ mv lapack_LINUX.a liblapack.a
$ mv tmg_LINUX.a libtmg.a

LAPACK95のインストール

以上で線形代数の計算にはLAPACKを利用できますが，LAPACKはFortran 77からも利用できるように書かれているので？，せっかくFortran 90を使っているのにLAPACKを呼び出して使おうとすると引数が多かったりといまいちに感じるかもしれません。そこでLAPACKをFortran95の様式で利用できるようにしたLAPACK95もインストールすることにします。

残念ながらLAPACK95はyumでインストールすることができませんので，netlibからソースをとってきてインストールします。netlibのLAPACK95のページからlapack95.tgzを適当なディレクトリにダウンロードして展開し，できた LAPACK95 に移動します。

$ tar zxvf lapack95.tgz
$ cd LAPACK95

いくつか準備をします。

$ mv SRC/makefile SRC/Makefile
$ mkdir lapack95_modules

make が参照する make.inc を次のように編集します。

FC = gfortran
FC1 = gfortran

OPTS0 = -O3
MODLIB = -I./../lapack95_modules
OPTS1 = -c $(OPTS0)
OPTS3 = $(OPTS1) $(MODLIB)
OPTL = -o
OPTLIB =

LAPACK_PATH = /usr/lib

LAPACK95 = ../lapack95.a
LAPACK77 = $(LAPACK_PATH)/liblapack.so
TMG77 = 
BLAS = $(LAPACK_PATH)/libblas.so

Intel Fortranを使う場合は上記のうちコンパイラを指定している部分を以下のようにします。

FC = ifort
FC1 = ifort

SRCに移動して make を実行します。

$ cd SRC
$ make single_double_complex_dcomplex

複素数型のサブルーチンを使わない場合は下のようにします。

$ make single_double

コンパイルが終了すると LAPACK95 ディレクトリに lapack95.a ができているはずです。また，LAPACK95/lapack95_modules には，関連する mod ファイルができています。これらのファイルは，任意のディレクトリに置いておけば使えますが，今回はライブラリは /usr/local/lib に，mod ファイルは /usr/local/include に置くことにします。また，ライブラリは指定しやすいように名前を変えておきます。

# cp /usr/local/src/LAPACK95/lapack95.a /usr/local/lib/liblapack95.a
# cp /usr/local/src/LAPACK95/lapack95_modules/*.mod /usr/local/include

管理者ではなく，個人用ライブラリとして作成することも可能です。上記の作業を任意のディレクトリで行い，できたライブラリと mod ファイルは ~/lib と ~/include にでも配置すればよいでしょう。

FFTW3のインストール

地震動の解析などではフーリエ変換が必要になることもあるので，DFT のためのライブラリをインストールします。

FFTW は1次元または多次元の任意のデータサイズについて DFT を行うことのできる C のライブラリですが，Fortranからも利用できます。また，単体で高速なだけでなく Multi Thread や MPI にも対応可能だそうです。

これも yum でインストールすることができます。（CentOSでも4.xまではyumでインストールできるのはFFTWのver2の系列なので下に示す方法でソースからコンパイルします。）

# yum install fftw3

このFFTWをFortranで利用するにはソースの最初にincludeするヘッダファイルが必要なのですが，これがどこにあるかわからず・・・。しかたないのでソースをダウンロードしてインストールすることにします。ダウンロードページから最新版であるfftw-3.1.2.tar.gzをダウンロードして展開後，configure，makeと進みます。

$ tar zxvf fftw-3.1.2.tar.gz
$ cd fftw-3.1.2
$ ./configure
$ make

コンパイル終了後，最後に以下のようにすることで /usr/local/include にヘッダが，/usr/local/libにライブラリがインストールされます。ライブラリのほうは yum でインストールしたものを使うなら削除してしまってもかまいません。

# make install

Octaveのインストール

基本的にプログラムはFortranで作ることにしていますが，ちょっとした計算などにはやはりMATLABライクな環境があると便利です。ということで，Octaveをインストールします。上にも書いたようにプログラムの実行速度はかなり遅いですが，X window（GUI環境）を使っているのであれば計算結果のグラフ描画などもちゃんとできます。（gnuplot を利用しているようです。）

これも yum でインストールできるので簡単です。Octave に必要な LAPACK などは既に導入済みですが，その他にもスパース行列用のライブラリなどは勝手に依存関係を調べて不足があればあわせてインストールされます。

# yum install octave

このページの先頭

last update on Jan. 31, 2008

LAPACKを使ったマトリックスの計算

LAPACKを使ったプログラムのコンパイル

LAPACKのライブラリを呼び出すプログラムは次のようにしてコンパイルします。使用するライブラリは'-l'（小文字のエル）オプションによって指定します。指定する順序は必ず下の例の順序でないといけません。'-o'オプションで作成されるファイル名の指定をしなければ同じディレクトリに'a.out'という名前の実行ファイルができます。

$ gfortran test.f90 -llapack -lblas

できたプログラムの実行は次のようにして行います。

$ ./a.out

'/usr/lib/'，'/usr/local/lib/'以外にライブラリが置いてある場合は'-L'オプションでライブラリの置いてあるディレクトリを指定します。例えば自分用のライブラリとして作成し，'~/lib'にライブラリを置いてあるような場合は次のようにしてコンパイルします。

$ gfortran test.f90 -L~/lib -llapack -lblas

LAPACKのルーチンの呼び出しにLAPACK95を利用している場合はソースファイルの最初でLAPACK95のコンパイル時にあわせて作成されたモジュールを使う宣言をしなければなりません。

program test95

  use f95_lapack

...

LAPACK95を使う場合は，LAPACKだけの場合に加えてコンパイル時にlapack95のライブラリの指定と，使用するモジュールファイル（.mod）が置いてある場所の'-I'（大文字のアイ）オプションによる指定が必要です。lapack95のライブラリはlapackのライブラリの前に指定します。

$ gfortran test95.f90 -I/usr/local/include -llapack95 -llapack -lblas

ライブラリやモジュールファイル（.mod）が他の場所にある場合は，ライブラリの場所はLAPACKだけの場合と同様に'-L'で指定し，'-I'による指定も変更します。ライブラリは'~/lib'に，モジュールファイル（.mod）は'~/include'にある場合は次のようになります。

$ gfortran test95.f90 -I~/include -L~/lib -llapack95 -llapack -lblas

なお，Intel Fortranを使う場合も同様の方法で指定できます。

Intel Math Kernel Library (MKL)

計算ツールのインストールの項では触れませんでしたが，Intel Fortran Compiler と同様 MKL も非商用であれば無償で使用できます。この MKL に含まれるLAPACKを使用する場合，次のように指定します。（MKL が /opt/intel/mkl/9 にインストールされている場合）

$ gfortran test.f90 -L/opt/intel/mkl/9/lib/32 -lmkl_lapack -lmkl

LAPACK95のインターフェイスも用意されていますが，私はまだ試していません。

連立一次方程式を解く

建築の構造分野で扱う構造解析の問題の代表的なものの1つに，{P}=[K]{d} の形をした連立方程式（{P}が外力，[K]が系の剛性マトリックス，{d}が変位）を解いて，外力{P}に対する未知の変位{d}を求めるというものがあります。

簡単な例題ですが，右図のような3層のせん断系に水平力が作用する場合の各層の変位を求めてみましょう。なぜ剛性マトリックスが以下のソースファイルの例のようになるかは構造解析をきちんと勉強して下さい。

MATLAB（Octave）で解く

まずはMATLAB（Octave）で解くにはどうするか？これはとても簡単です。m-ファイルにするとこんな感じになります。

（m-ファイルの例：leqsol.m）

% 連立一次方程式の例
K = [ 2, -1,  0;...
     -1,  2, -1;...
      0, -1,  1]

P = [ 1;...
      2;...
      3]

d = K \ P

Octaveで実行してみます。

octave:1> leqsol
K =

   2  -1   0
  -1   2  -1
   0  -1   1

P =

   1
   2
   3

d =

    6.0000
   11.0000
   14.0000

LAPACK95を使う

ではLAPACKを使ってFortranでコーディングするにはどうするか？まずはLAPACK95を使った例です。この例のようにフレームの構造解析や有限要素法では，境界条件も導入した最終的な全体剛性マトリックスは正定値対称行列になります（非線形解析ではこの限りではありません）。そこで，LAPACKのルーチンのうちコレスキー分解を利用するものを使うことにします。ただ必要なサブルーチンを呼び出すだけです。十分シンプルに書けると思いますがいかがでしょう？

（LAPACK95を用いたソースファイルの例：leqsol95.f90）

program leqsol95

! LAPACK95のコンパイル時にあわせて作成されたモジュールを使います。
  use f95_lapack

  implicit none

  integer, parameter :: N = 3
  integer :: i
  real(8) :: K(1:N,1:N), P(1:N)

  K(1:3,1) = (/  2.0d0, -1.0d0,  0.0d0 /)
  K(1:3,2) = (/ -1.0d0,  2.0d0, -1.0d0 /)
  K(1:3,3) = (/  0.0d0, -1.0d0,  1.0d0 /)

  write(*,'(A)') 'K ='
  write(*,'(3F12.4)') ( K(i,:), i=1,N )

  P(1:3)   = (/  1.0d0,  2.0d0,  3.0d0 /)

  write(*,'(A)') 'P ='
  write(*,'(F12.4)') P

! LAPACK95（ラッパー）を介してLAPACKのルーチンを呼び出します。
! K に K を上三角行列にコレスキー分解した結果が，
! P に連立一次方程式の解が格納されます。
! K や P を書き換えられたくない場合は注意しましょう。
  call LA_POSV( K, P )

  write(*,'(A)') 'd ='
  write(*,'(F12.4)') P

end program leqsol95

まずコンパイルします。

$ gfortran leqsol95.f90 -I/usr/local/include -llapack95 -llapack -lblas

実行してみます。

$ ./a.out
K =
      2.0000     -1.0000      0.0000
     -1.0000      2.0000     -1.0000
      0.0000     -1.0000      1.0000
P =
      1.0000
      2.0000
      3.0000
d =
      6.0000
     11.0000
     14.0000

LAPACKを直接呼び出す

LAPACK95を使うにはそのプログラムのコンパイルに使うコンパイラで作成されたモジュールを用意しなければならないし，実際に計算をするルーチンを呼び出すまでに内部で踏むステップは1ステップ多くなってしまいます。そのため，直接LAPACKのルーチンを呼び出したいようなときもあるかもしれません。

LAPACKを直接呼び出そうと思うとこんな感じになります。同じくただサブルーチンを呼び出すだけですみますが，やや引き数が多くなってしまいます。

（LAPACKを直接呼び出すソースファイルの例：leqsol.f90）

program leqsol

  implicit none

  integer, parameter :: N = 3
  integer :: i
  real(8) :: K(1:N,1:N), P(1:N)

  integer :: INFO

  K(1:3,1) = (/  2.0d0, -1.0d0,  0.0d0 /)
  K(1:3,2) = (/ -1.0d0,  2.0d0, -1.0d0 /)
  K(1:3,3) = (/  0.0d0, -1.0d0,  1.0d0 /)

  write(*,'(A)') 'K ='
  write(*,'(3F12.4)') ( K(i,:), i=1,N )

  P(1:3)   = (/  1.0d0,  2.0d0,  3.0d0 /)

  write(*,'(A)') 'P ='
  write(*,'(F12.4)') P

! LAPACKのルーチンを呼び出します。
! K に K を上三角行列にコレスキー分解した結果が，
! P に連立一次方程式の解が格納されます。
! K や P を書き換えられたくない場合は注意しましょう。
  call DPOSV('U', N, 1, K, N, P, N, INFO )

  write(*,'(A)') 'd ='
  write(*,'(F12.4)') P

end program leqsol

コンパイルしてみます。

$ gfortran leqsol.f90 -llapack -lblas

実行してみます。とうぜん同じ結果になります。

$ ./a.out
K =
      2.0000     -1.0000      0.0000
     -1.0000      2.0000     -1.0000
      0.0000     -1.0000      1.0000
P =
      1.0000
      2.0000
      3.0000
d =
      6.0000
     11.0000
     14.0000

固有値解析を行う

構造物の地震応答解析など振動の問題を扱う場合，対象とする系の固有値解析を行い得られた固有値から固有周期を求めたり，固有ベクトルから得られる系の振動モードを解析に用いたりします。

この場合，系の質量マトリックスを[M]，剛性マトリックスを[K]とし，系の固有角振動数ωの2乗をλとおけば，[K]{x}=λ[M]{x}の形をした一般化固有値問題を解くことになります。

では，右のような簡単なせん断系の固有値解析を行ってみます。

MATLAB（Octave）で解く

まずはMATLAB（Octave）で解くにはどうするか？これはとても簡単です。m-ファイルにするとこんな感じになります。ただ，Octaveは一般化固有値問題の形のままでは解くことができないので，標準固有値問題の形に変形します。

（m-ファイルの例：eigsol.m）

% 固有値解析の例

  K = [ 2, -1,  0;...
       -1,  2, -1;...
        0, -1,  1]

  M = [ 1,  0,  0;...
        0,  1,  0;...
        0,  0,  1]

% octaveは一般化固有値問題を解くことができないので標準固有値問題に変換します。
% 最も簡単な方法はこのように K を M で左割りしてしまう方法です。
% U に規格化された固有ベクトルが，D の対角要素に固有値が昇順に格納されます。
% [U, D] = eig( M \ K )

% 上の方法では一般に M\K は対称ではなくなるため必ずしもよい方法とはいえません。
% (octaveやMATLABは解いてくれますが精度は必ずしも確保されません)
% 一般には M のコレスキー分解を利用することが多いようです。

  R = chol(M);                  % M=R'*R の形にコレスキー分解します
  [U, D] = eig( R'\ K / R );    % KK=R'\K/R, V=R*U として D*KK = D*V を解きます
  D                             % 固有値は元の一般化固有値問題と等しくなります
  U = R \ U                     % 元の一般化固有値問題の固有ベクトルに変換します

% MATLABは一般化固有値問題をそのまま解くことができます。
% [U, D] = eig( K, M )

Octaveで実行してみます。

octave:1> eigsol
K =

   2  -1   0
  -1   2  -1
   0  -1   1

U =

  -0.32799   0.73698  -0.59101
  -0.59101   0.32799   0.73698
  -0.73698  -0.59101  -0.32799

D =

   0.19806   0.00000   0.00000
   0.00000   1.55496   0.00000
   0.00000   0.00000   3.24698

LAPACK95を使う

次はLAPACK95を使ったプログラムです。質量マトリックスも剛性マトリックスも対称行列になるので，実対称行列用のルーチンを使うことにします。厳密には，このルーチンを使うためには質量マトリックスは正定値でなければなりません（途中で質量マトリックスをコレスキー分解するため）。通常の骨組解析の解析のように質量を節点への集中質量として与える場合は問題ありませんが，有限要素法において質量を要素内の分布質量とした場合，その与え方によっては正定値でなくなることもあるそうなので注意が必要です。

この場合も1行サブルーチンを呼び出すだけですみます。

（LAPACK95を用いたソースファイルの例：eigsol95.f90）

program eigsol95

! LAPACK95のコンパイル時にあわせて作成されたモジュールを使います。
  use f95_lapack

  implicit none

  integer, parameter :: N = 3
  integer :: i
  real(8) :: K(1:N,1:N), M(1:N,1:N), D(1:N)

  K(1:3,1) = (/  2.0d0, -1.0d0,  0.0d0 /)
  K(1:3,2) = (/ -1.0d0,  2.0d0, -1.0d0 /)
  K(1:3,3) = (/  0.0d0, -1.0d0,  1.0d0 /)

  write(*,'(A)') 'K ='
  write(*,'(3F12.4)') ( K(i,:), i=1,N )

  M(1:3,1) = (/  1.0d0,  0.0d0,  0.0d0 /)
  M(1:3,2) = (/  0.0d0,  1.0d0,  0.0d0 /)
  M(1:3,3) = (/  0.0d0,  0.0d0,  1.0d0 /)

  write(*,'(A)') 'M ='
  write(*,'(3F12.4)') ( M(i,:), i=1,N )

! LAPACK95（ラッパー）を介してLAPACKのルーチンを呼び出します。
! K(:,i)に i 番目の固有ベクトルが，D に固有値が昇順に格納されます。
! M には M を上三角行列にコレスキー分解した結果が格納されています。
! K や M を書き換えられたくない場合は注意しましょう。
  call LA_SYGV( K, M, D, 1,'V' )

! 固有ベクトルが必要ない場合は5番目の引き数を 'N' とします。
! もしくは4番目以降の引き数は省略してもかまいません。
! この場合にも K の対角要素を含む上三角要素は書き換えられており
! M には M を上三角行列にコレスキー分解した結果が格納されます。
! K や M を書き換えられたくない場合は注意しましょう。
! call LA_SYGV( K, M, D, 1, 'N' )
! or
! call LA_SYGV( K, M, D )

  write(*,'(A)') 'U ='
  write(*,'(3F12.4)') ( K(i,:), i=1,N )

  write(*,'(A)') 'D ='
  write(*,'(3F12.4)') D

end program eigsol95

コンパイルして実行してみます。

$ gfortran eigsol95.f90 -I/usr/local/include -llapack95 -llapack -lblas
$ ./a.out
K =
      2.0000     -1.0000      0.0000
     -1.0000      2.0000     -1.0000
      0.0000     -1.0000      1.0000
M =
      1.0000      0.0000      0.0000
      0.0000      1.0000      0.0000
      0.0000      0.0000      1.0000
U =
     -0.3280      0.7370     -0.5910
     -0.5910      0.3280      0.7370
     -0.7370     -0.5910     -0.3280
D =
      0.1981      1.5550      3.2470

LAPACKを直接呼び出す

LAPACKを直接呼び出して同じ問題を解こうとするとこんな感じになります。この場合，サブルーチンを1回呼び出すだけではなく，いくつかの処理が必要になります。また，扱える行列についての注意点はLAPACK95の場合と同じです。（LAPACK95ではラッパーがここでしているような処理をしているだけで，同じLAPACKのルーチンを使ってます。）

（LAPACKを直接呼び出すソースファイルの例：eigsol.f90）

program eigsol

  implicit none

  integer, parameter :: N = 3
  integer :: i
  real(8) :: K(1:N,1:N), M(1:N,1:N), D(1:N)

  integer              :: LWORK, INFO
  real(8)              :: LWORK0
  real(8), allocatable :: WORK(:)

  K(1:3,1) = (/  2.0d0, -1.0d0,  0.0d0 /)
  K(1:3,2) = (/ -1.0d0,  2.0d0, -1.0d0 /)
  K(1:3,3) = (/  0.0d0, -1.0d0,  1.0d0 /)

  write(*,'(A)') 'K ='
  write(*,'(3F12.4)') ( K(i,:), i=1,N )

  M(1:3,1) = (/  1.0d0,  0.0d0,  0.0d0 /)
  M(1:3,2) = (/  0.0d0,  1.0d0,  0.0d0 /)
  M(1:3,3) = (/  0.0d0,  0.0d0,  1.0d0 /)

  write(*,'(A)') 'M ='
  write(*,'(3F12.4)') ( M(i,:), i=1,N )


! まず最適なワークスペースの大きさを求めます。
! 固有ベクトルが必要ない場合は2番目の引き数を 'N' とします。
  call DSYGV( 1, 'V', 'U', N, K, N, M, N, D, LWORK0, -1, INFO )

! ワークスペース用の配列を確保します。
  LWORK = INT(LWORK0)
  allocate( WORK(1:LWORK) )

! LAPACKのルーチンを呼び出します。
! K(:,i)に i 番目の固有ベクトルが，D に固有値が昇順に格納されます。
! M には M を上三角行列にコレスキー分解した結果が格納されています。
! K や M を書き換えられたくない場合は注意しましょう。
! 固有ベクトルが必要ない場合は2番目の引き数を 'N' とします。
  call DSYGV( 1, 'V', 'U', N, K, N, M, N, D, WORK, LWORK, INFO )

! ワークスペース用の配列を解放します。
  deallocate( WORK )

  write(*,'(A)') 'U ='
  write(*,'(3F12.4)') ( K(i,:), i=1,N )

  write(*,'(A)') 'D ='
  write(*,'(3F12.4)') D

end program eigsol

コンパイルして実行してみます。とうぜん同じ結果になります。

$ gfortran eigsol.f90 -llapack -lblas
$ ./a.out
K =
      2.0000     -1.0000      0.0000
     -1.0000      2.0000     -1.0000
      0.0000     -1.0000      1.0000
M =
      1.0000      0.0000      0.0000
      0.0000      1.0000      0.0000
      0.0000      0.0000      1.0000
U =
     -0.3280      0.7370     -0.5910
     -0.5910      0.3280      0.7370
     -0.7370     -0.5910     -0.3280
D =
      0.1981      1.5550      3.2470

逆行列を求める

[A]{x}={b}の形の連立一次方程式を解く場合，上の例のように数値解法によって直接解{x}を求めるので，逆行列を求めること（[A]の逆行列を求めて{x}=[A]**-1 {b}のようにする）は通常あまりしません。

だからといって，構造解析において逆行列を求める必要がないかというとそういうわけではないので，逆行列を求める例を示しておきます。（上の例で使った剛性マトリックスの逆行列を求めてみます。）

MATLAB（Octave）で求める

まずはMATLAB（Octave）で求めてみます。これはとても簡単です。m-ファイルにするとこんな感じになります。

（m-ファイルの例：matinv.m）

% 逆行列を求めるの例

  K = [ 2, -1,  0;...
       -1,  2, -1;...
        0, -1,  1]

  F = inv(K)

Octaveで実行してみます。

octave:1> matinv
K =

   2  -1   0
  -1   2  -1
   0  -1   1

F =

  1.00000  1.00000  1.00000
  1.00000  2.00000  2.00000
  1.00000  2.00000  3.00000

LAPACK95を使って求める

LAPACKには逆行列を直接求めるルーチンはありません。また，連立方程式を解いたときのように対称性を利用してコレスキー分解を使いところですが，LAPACK95にはコレスキー分解をするルーチンはありますが，分解した上（下）三角行列を使って一発で逆行列を求めるルーチンも用意されていません。そこで，一般の場合と同様にLU分解をし，それ使って逆行列を求めることにします。

（LAPACK95で逆行列を求める例：matinv95.f90）

program matinv95

! LAPACK95のコンパイル時にあわせて作成されたモジュールを使います。
  use f95_lapack

  implicit none

  integer, parameter :: N = 3
  integer :: i
  real(8) :: K(1:N,1:N)

  integer :: IPIV(1:N)

  K(1:3,1) = (/  2.0d0, -1.0d0,  0.0d0 /)
  K(1:3,2) = (/ -1.0d0,  2.0d0, -1.0d0 /)
  K(1:3,3) = (/  0.0d0, -1.0d0,  1.0d0 /)

  write(*,'(A)') 'K ='
  write(*,'(3F12.4)') ( K(i,:), i=1,N )

! LAPACK95（ラッパー）を介してLAPACKのルーチンを呼び出して行列をLU分解します。
! K には分解した L と U があわせて格納されています。
! （Lの対角要素（全て1）は格納されません。）
! K を書き換えられたくない場合は注意が必要です。
! IPIVにはLU分解をに用いたピボットの情報が格納されます。
  call LA_GETRF( K, IPIV )

! GETRFでLU分解した結果を用いてもとの行列の逆行列を求めます。
! もとの行列 K の逆行列が K に格納されます。
  call LA_GETRI( K, IPIV )

  write(*,'(A)') 'F ='
  write(*,'(3F12.4)') ( K(i,:), i=1,N )

end program matinv95

コンパイルして実行してみます。

$ gfortran matinv95.f90 -I/usr/local/include -llapack95 -llapack -lblas
$ ./a.out
K =
      2.0000     -1.0000      0.0000
     -1.0000      2.0000     -1.0000
      0.0000     -1.0000      1.0000
F =
      1.0000      1.0000      1.0000
      1.0000      2.0000      2.0000
      1.0000      2.0000      3.0000

LAPACKで求める

LAPACKにはコレスキー分解した結果を用いて逆行列を求めるルーチンも用意されているので，まずはコレスキー分解を用いる方法です。ただし，結果は（対称なので）上三角要素の範囲のみに格納されています。

（LAPACKのコレスキー分解のルーチンを用いて逆行列を求める例：matinv_ch.f90）

program matinv_ch

  implicit none

  integer, parameter :: N = 3
  integer :: i
  real(8) :: K(1:N,1:N)

  integer :: INFO

  K(1:3,1) = (/  2.0d0, -1.0d0,  0.0d0 /)
  K(1:3,2) = (/ -1.0d0,  2.0d0, -1.0d0 /)
  K(1:3,3) = (/  0.0d0, -1.0d0,  1.0d0 /)

  write(*,'(A)') 'K ='
  write(*,'(3F12.4)') ( K(i,:), i=1,N )

! LAPACKのルーチンを呼び出して行列をコレスキー分解します。
! K には上三角行列にコレスキー分解した結果が格納されます。
! K を書き換えられたくない場合は注意が必要です。
  call DPOTRF( 'U', N, K, N, INFO )

! LAPACKのルーチンを直接呼び出して
! POTRFでコレスキー分解した結果を用いてもとの行列の逆行列を求めます。
! もとの行列 K の逆行列が上三角要素に格納されています。
  call DPOTRI( 'U', N, K, N, INFO )

! 下半分の要素を埋めます。
  do i=1,N
    K(i+1:N,i) = K(i,i+1:N)
  end do

  write(*,'(A)') 'F ='
  write(*,'(3F12.4)') ( K(i,:), i=1,N )

end program matinv_ch

コンパイルして実行してみます。とうぜん同じ結果になります。

$ gfortran matinv_ch.f90 -llapack -lblas
$ ./a.out
K =
      2.0000     -1.0000      0.0000
     -1.0000      2.0000     -1.0000
      0.0000     -1.0000      1.0000
F =
      1.0000      1.0000      1.0000
      1.0000      2.0000      2.0000
      1.0000      2.0000      3.0000

LAPACK95のときと同じようにLU分解をして逆行列を求めようとするとこうなります。

（LAPACKのLU分解のルーチンを用いて逆行列を求める例：matinv_lu.f90）

program matinv_lu

  implicit none

  integer, parameter :: N = 3
  integer :: i
  real(8) :: K(1:N,1:N)

  integer              :: IPIV(1:N), LWORK, INFO
  real(8)              :: LWORK0
  real(8), allocatable :: WORK(:)

  K(1:3,1) = (/  2.0d0, -1.0d0,  0.0d0 /)
  K(1:3,2) = (/ -1.0d0,  2.0d0, -1.0d0 /)
  K(1:3,3) = (/  0.0d0, -1.0d0,  1.0d0 /)

  write(*,'(A)') 'K ='
  write(*,'(3F12.4)') ( K(i,:), i=1,N )

! LAPACKのルーチンを呼び出して行列をLU分解します。
! LAPACK95（ラッパー）を介してLAPACKのルーチンを呼び出して行列をLU分解します。
! K には分解した L と U があわせて格納されています。
! （Lの対角要素（全て1）は格納されません。）
! K を書き換えられたくない場合は注意が必要です。
! IPIVにはLU分解をに用いたピボットの情報が格納されます。
  call DGETRF( N, N, K, N, IPIV, INFO )

! LU分解された行列の逆行列を求めるのに最適なワークスペースの大きさを求めます。
  call DGETRI(N, K, N, IPIV, LWORK0, -1, INFO)

! ワークスペース用の配列を確保します。
  LWORK = INT(LWORK0)
  allocate( WORK(1:LWORK) )

! GETRFでLU分解した結果を用いてもとの行列の逆行列を求めます。
! もとの行列 K の逆行列が K に格納されます。
  call DGETRI(N, K, N, IPIV, WORK, LWORK, INFO)

! ワークスペース用の配列を解放します。
  deallocate( WORK )

  write(*,'(A)') 'F ='
  write(*,'(3F12.4)') ( K(i,:), i=1,N )

end program matinv_lu

コンパイルして実行してみます。とうぜんこれも同じ結果になります。

$ gfortran matinv_lu.f90 -llapack -lblas
$ ./a.out
K =
      2.0000     -1.0000      0.0000
     -1.0000      2.0000     -1.0000
      0.0000     -1.0000      1.0000
F =
      1.0000      1.0000      1.0000
      1.0000      2.0000      2.0000
      1.0000      2.0000      3.0000

このページの先頭

last update on Dec. 24, 2007

FFTW3を使ったフーリエ変換

FFTWは自称西側（古い言い方ですが）で最速のFFTだそうです。使う上で注意すべきことは，FFTWによるDFTではフーリエ変換（forward）時もフーリエ逆変換（backward）時も結果は規格化されていないということです（規格化定数は乗じられていない）。従って forward 処理により得られたフーリエ係数をそのまま backward 処理（逆変換）すると，元のデータの N 倍（N：データ数）になります。

FFTW3を使ったプログラムのコンパイル

FFTW3のライブラリを呼び出す場合は，ソースファイル中（変数宣言部の次）に必要な変数などを定義しているヘッダーファイルを組み込まなければなりません。

program test

<変数宣言部>

  include 'fftw3.f'

...

FFTW3のライブラリを呼び出すプログラムは次のようにしてコンパイルします。（ライブラリが'/usr/lib/'または'/usr/local/lib/'にある場合）使用するライブラリは'-l'（小文字のエル）オプションによって指定します。また，使用するヘッダーファイル（fftw3.f）が置いてある場所を'-I'（大文字のアイ）オプションにより指定します。'-o'オプションを指定していなければ'a.out'という実行ファイルができます。

$ gfortran test.f90 -I/usr/local/include -lfftw3

できたプログラムの実行は次のようにして行います。

$ ./a.out

ライブラリやヘッダーファイルが他の場所にある場合は，ライブラリの場所をLAPACKの場合と同様に'-L'で指定し，'-I'による指定も変更します。ライブラリは'~/lib'に，ヘッダーファイルは'~/include'にある場合は次のようになります。

$ gfortran test.f90 -I~/include -L~/lib -lfftw3

なお，Intel Fortranを使う場合も同様の方法で指定できます。また，MKLにもFFTを行うためのライブラリがあるようですが，私は使ったことがないので詳細はわかりません。

加速度時刻歴をフーリエ変換する

地震計や常時微動計測などで得られた加速度時刻歴をフーリエ変換することを考えます。フーリエ変換する目的はさまざまですが，ここではフーリエ振幅スペクトルとフーリエ位相スペクトルを求めるまでを行います。加速度時刻歴は次のように'acc.dat'というファイに縦一列（幅12桁）に記録されているとします。（例で用いるのは私が作成した模擬地震動です。）

（例で用いる加速度のデータ：acc.dat）

また，MATLAB（Octave）の関数（実はFFTWを利用している），FFTW3ともフーリエ変換するデータの長さは2の累乗ではなくてもかまわないのですが，以下の例ではデータの長さが2の累乗になるように（この場合が変換が速いので）加速度時刻歴の末尾にゼロのデータを付け足しています。

MATLAB（Octave）で求める

まずはMATLAB（Octave）で求めてみます。MATLAB（Octave）の関数 fft は，フーリエ変換（FFT）時の規格化定数は 1，指数関数部の符号は負です。

m-ファイルにするとこんな感じになります。フーリエ変換するデータの長さ N を求めた後，関数 fft によりフーリエ変換し，得られた複素フーリエ係数 C からフーリエ振幅 F とフーリエ位相 P を求めています。

（離散フーリエ変換の例：dft.m）

% 離散フーリエ変換（DFT）の例

% 加速度データを読み込む
  acc = load('acc.dat');

% フーリエ変換するデータの長さ NN を 2 の累乗になるように決める
  N = 2^ceil( log2( length(acc) ) );

% FFTを実行
  C = fft(acc, N);

% フーリエ振幅スペクトル，フーリエ位相スペクトルを求める
  F = abs(C);
  P = angle(C);

% フーリエ振幅スペクトル，フーリエ位相スペクトルの書き出し
% fid = fopen('fas.out','w');
% fprintf(fid, '%12.4e\n', F(1:N/2+1) )
% fclose(fid);

% fid = fopen('fps.out','w');
% fprintf(fid, '%12.5f\n', P(1:N/2+1) )
% fclose(fid);

Octaveで実行してみます。

octave:1> dft
octave:2> F
ans =

   6.8249e+03
   8.8953e+03
   9.4913e+03
   8.7184e+03
   8.7203e+03
   1.0009e+04
   7.7715e+03
   1.0737e+04
   1.2446e+04
   5.9300e+03
...

octave:3> P
ans =

   3.14159
   0.58705
  -1.68769
   2.14171
   0.39895
  -1.28268
   2.57205
   0.60957
  -0.88621
  -1.44486
...

FFTW3を用いて求める

FFTW3を用いて求めてみます。FFTW3には1次元の実数データをフーリエ変換する方法が用意されているので，それを使うことにします。実数データのフーリエ変換では得られう複素フーリエ係数はナイキスト振動数をはさんで共役になるため，このルーチンではFFTの結果はナイキスト振動数までの複素フーリエ係数が格納されます。

前述のようにFFTWでのFFTでは，結果の規格化はされていませんが，一般の離散フーリエ変換（forward処理）では規格化定数は1とすることが多いため，この例では規格化の処理をしていません。

以下の例では，'acc.dat'から加速度データを読み込み，FFTにより求めたフーリエ振幅とフーリエ位相をそれぞれ'fas.out'と'fps.out'に書き出します。

（FFTWによる離散フーリエ変換の例：dft.f90）

program dft

  implicit none

  integer                 :: N, ND, i, ios
  real(8)   , allocatable :: acc(:), F(:), P(:)
  complex(8), allocatable :: C(:)

  character(7), parameter :: inp ='acc.dat', out1='fas.out', out2='fps.out'

  integer :: plan(8)

! FFTW3を呼び出すのに必要なヘッダーファイルを include する
  include 'fftw3.f'

! ファイル（データ）の長さ ND を調べる
  open(51,FILE=inp,ACTION='read')
  ND = 0
  do
    read(51,'(F12.0)',IOSTAT=ios)
    if (ios<0) exit	! ファイルの末尾に来たらループを抜ける
    ND = ND + 1
  end do
  close(51)

! フーリエ変換するデータの長さ N を 2 の累乗になるように決めて
! 加速度データ，複素フーリエ係数，フーリエ振幅，フーリエ位相の配列を確保
  N = 2**int( log( dble(ND) ) / log( 2.0d0 ) + 0.5d0 )
  allocate( acc(1:N), C(0:N/2), F(0:N/2), P(0:N/2) )

! 加速度データを読み込む
  open(51,FILE=inp,ACTION='read')
  read(51,'(F12.3)') acc(1:ND)
  close(51)

! 加速度データの末尾（不足分）は 0 にする
  acc(ND+1:N) = 0.0d0

! FFTW3 を呼び出してフーリエ変換を行う
! 1) plan を作成する（1次元実数データ用）
  call dfftw_plan_dft_r2c_1d( plan, N, acc, C, FFTW_ESTIMATE )

! 2) FFT を実行する
  call dfftw_execute(plan)

! 3) plan を破棄する
  call dfftw_destroy_plan(plan)

! フーリエ振幅スペクトル，フーリエ位相スペクトルを求める
  F = ABS( C )
  P = ATAN2( AIMAG(C), DBLE(C) )

! フーリエ振幅スペクトル，フーリエ位相スペクトルの書き出し
  open(61,FILE=out1,STATUS='replace',ACTION='write')
  write(61,'(E12.4)') F
  close(61)

  open(61,FILE=out2,STATUS='replace',ACTION='write')
  write(61,'(F12.5)') P
  close(61)

end program dft

コンパイルして実行してみます。

$ gfortran dft.f90 -I/usr/local/include -lfftw3
$ ./a.out

結果をターミナルに表示してみます。MATLAB（Octave）と同じ結果です。

$ cat fas.out
  0.6825E+04
  0.8895E+04
  0.9491E+04
  0.8718E+04
  0.8720E+04
  0.1001E+05
  0.7771E+04
  0.1074E+05
  0.1245E+05
  0.5930E+04
...

$ cat fps.out
     3.14159
     0.58705
    -1.68769
     2.14171
     0.39895
    -1.28268
     2.57205
     0.60957
    -0.88621
    -1.44486
...

複素フーリエ係数をフーリエ逆変換する

例えば，得られたデータ列に対して周波数領域でフィルタ処理（バンドパスフィルタなど）をしたデータ列を得たい場合，フィルタを乗じたフーリエ係数から時間領域のデータ列をフーリエ逆変換により求める必要があります。

ここでは，単に先ほど得たフーリエ振幅とフーリエ位相から加速度時刻歴を復元する例を示します。

MATLAB（Octave）で求める

まずはMATLAB（Octave）で求めてみます。MATLAB（Octave）のフーリエ逆変換を行う関数 ifft は，変換時の規格化定数は 1/N，指数関数部の符号は正です。

m-ファイルにするとこんな感じになります。'fas.out'と'fps.out'から読み込んだフーリエ振幅とフーリエ位相から複素フーリエ係数を作り，ifft によりフーリエ逆変換をして加速度時刻歴を得ます。但し，ifft の結果は複素数（虚部は 0）のためその実部のみを取り出しています。

（離散フーリエ逆変換の例：idft.m）

% 離散フーリエ逆変換（IDFT）の例

% フーリエ振幅スペクトル，フーリエ位相スペクトルを読み込む
  F = load('fas.out');
  P = load('fps.out');

% 複素フーリエ係数を求める
  N = (length(F)-1) * 2;

  C = F .* complex( cos(P), sin(P) );
  C(N/2+2:N) = conj( C(N/2:-1:2) );

% FFTを実行し，実部を取り出す
  acc = real( ifft(C) );

% 加速度時刻歴の書き出し
% fid1 = fopen('acc.out','w');
% fprintf(fid1,'%12.5f\n',acc);
% fclose(fid1);

Octaveで実行してみます。

octave:1> idft
octave:2> acc
ans =

   0.118271
   0.110924
   0.101463
   0.097692
   0.108344
   0.105536
   0.097079
   0.136624
   0.155325
   0.208218
...

FFTW3を用いて求める

では，FFTW3を用いて求めてみます。1次元のフーリエ変換（逆変換）するルーチンを用います。また，FFTWでのフーリエ逆変換では結果の規格化はされないので，得られた結果に対して規格化定数1/Nを乗じています。

フーリエ振幅とフーリエ位相をそれぞれ'fas.out'と'fps.out'から読み込み，求めた加速度時刻歴を'acc.out'に書き出すプログラムの例を示します。

（FFTWによる離散フーリエ逆変換の例：idft.f90）

program idft

  implicit none

  integer                 :: N, N1, i, ios
  real(8)   , allocatable :: acc(:), F(:), P(:)
  complex(8), allocatable :: C(:), cacc(:)

  character(7), parameter :: inp1='fas.out', inp2='fps.out', out ='acc.out'

  integer :: plan(8)

! FFTW3を呼び出すのに必要なヘッダーファイルを include する
  include 'fftw3.f'

! ファイル（データ）の長さ N を調べる
  open(51,FILE=inp1,ACTION='read')
  N = 0
  do
    read(51,'(F12.0)',IOSTAT=ios)
    if (ios<0) exit	! ファイルの末尾に来たらループを抜ける
    N = N + 1
  end do
  close(51)

! フーリエ逆変換により作成される時刻歴データの長さ N を求め
! 加速度データ，複素フーリエ係数，フーリエ振幅，フーリエ位相の配列を確保
  N = (N-1) * 2
  allocate( acc(1:N), cacc(1:N), C(0:N-1), F(0:N/2), P(0:N/2) )

! フーリエ振幅スペクトル，フーリエ位相スペクトルを読み込む
  open(51,FILE=inp1,ACTION='read')
  read(51,'(E12.4)') F(0:N/2)
  close(51)

  open(51,FILE=inp2,ACTION='read')
  read(51,'(F12.5)') P(0:N/2)
  close(51)

! 複素フーリエ係数を求める
  C = F * CMPLX( COS(P), SIN(P) );
  C(N/2+1:N-1) = CONJG( C(N/2-1:1:-1) );

! FFTW3 を呼び出してフーリエ逆変換を行う
! 1) plan を作成する
  call dfftw_plan_dft_1d( plan, N, C, cacc, FFTW_BACKWARD, FFTW_ESTIMATE )

! 2) FFT を実行する
  call dfftw_execute(plan)

! 3) plan を破棄する
  call dfftw_destroy_plan(plan)

! フーリエ逆変換で得られたもののを規格化し，実部を取り出す
  acc = DBLE( cacc / DBLE(N) )

! 加速度時刻歴の書き出し
  open(61,FILE=out,STATUS='replace',ACTION='write')
  write(61,'(F12.5)') acc
  close(61)

end program idft

コンパイルして実行してみます。

$ gfortran idft.f90 -I/usr/local/include -lfftw3
$ ./a.out

結果をターミナルに表示してみます。MATLAB（Octave）と同じ結果になります。元の'acc.dat'とはわずかに異なりますが，これは数値計算上の誤差や，特にフーリエ振幅とフーリエ位相を書き出すときの丸め誤差のせいです。

$ cat acc.out
     0.11827
     0.11092
     0.10146
     0.09769
     0.10834
     0.10554
     0.09708
     0.13662
     0.15533
     0.20822
...

このページの先頭

研究メモ

Fortran90で構造解析をしよう！

はじめに

なぜFortran90？

Linux(CentOS 5)のインストール

CentOS 5の設定

sshの設定

rootになれるユーザーを管理者のみにする

パッケージ管理ソフト yum のリポジトリの追加

CentOSPlusの追加

Fedora EPEL リポジトリを設定

fedora-extrasを使えるように設定する

yumの設定とソフトウェアのアップデート

Windowsとのファイル共有（Samba）の設定

手動インストール用のシンボリックリンクの設定

一般ユーザーの登録

Fortranコンパイラのインストール

gcc-gfortranのインストール

G95のインストール

Intel Fortran Compilerのインストール

数値計算ツールのインストール

LAPACKのインストール

LAPACK95のインストール

FFTW3のインストール

Octaveのインストール

LAPACKを使ったマトリックスの計算

LAPACKを使ったプログラムのコンパイル

Intel Math Kernel Library (MKL)

連立一次方程式を解く

MATLAB（Octave）で解く

LAPACK95を使う

LAPACKを直接呼び出す

固有値解析を行う

MATLAB（Octave）で解く

LAPACK95を使う

LAPACKを直接呼び出す

逆行列を求める

MATLAB（Octave）で求める

LAPACK95を使って求める

LAPACKで求める

FFTW3を使ったフーリエ変換

FFTW3を使ったプログラムのコンパイル

加速度時刻歴をフーリエ変換する

MATLAB（Octave）で求める

FFTW3を用いて求める

複素フーリエ係数をフーリエ逆変換する

MATLAB（Octave）で求める

FFTW3を用いて求める