いつもFotranを使っています.ただ今回共同研究者がC++でコードを書いてきたので,C++でコードを書かないといけない事案が発生しました.その時,コーディングのお約束が結構あって苦労したので次回から簡単に解決できるようにメモを残します.複数ファイルがあるような状態を想定しています.
以下でFortranはifort -extend-source
でコンパイルしています.C++はg++ -std=c++0x
でコンパイルしています.
namespace
なんでもかんでもグローバル変数にしないのがモダンなプログラミングの基本です.そのためにはnamespaceが役に立ちます.
Fortran
例えば以下のようにFotranで光速をcgs単位系で定義しているmoduleがあるとします(僕は長くとった固定幅で書く流儀です).
module cgs_units
implicit none
real(8),parameter:: c=2.997924562d+10
end module cgs_units
プログラム上では以下のように呼びますね.
program main
use cgs_units
implicit none
real(8):: kousoku
kousoku = c
end program
C++
これをC++で書くとどうなるか.ヘッダーファイルでは宣言(メモリを確保しない)だけして,定義(実際にメモリを確保する)は全部のプログラムで一回だけにしないといけません.なのでヘッダファイルではextern
を使って外部のファイルに変数が定義されていることを明示します.あとヘッダファイルは何回も読み込んではいけないので読み込みが一回きりになるように#ifndef _CGS_UNITS
を使っています.
# ifndef _CGS_UNITS
# define _CGS_UNITS
namespace cgs_units{
extern const double c;
};
# endif
プログラム上では以下のように呼びますね.まずどこかでc
を定義しておきます.Fortranでは数字の桁数をd10とか書いてましたが,e10に直しています.
# include <iostream>
# include "units.hpp"
namespace cgs_units{
const double c = 2.997924562e+10;
};
int main(){
using namespace std;
using namespace cgs_units;
double kousoku;
kousoku = c;
cout<< kousoku <<endl;
return 1;
}
多次元配列
C/C++の多次元配列の扱い方の情報は錯綜していて,なかなかベストなやり方が探せませんでした.コメント欄でのやりとりにあるように,C++では目的に応じていろんな多次元配列の取り扱いの仕方があります.目的にあったやり方をするのが良いでしょう.
Fortran的な直感的な多次元配列の取り扱いを望む場合,STLのvector classが便利です.既存の型ではなく,自分で作った型に対しても使えます.ただし,簡単な準備が必要です.用途によっては他のライブラリのほうが便利な場合もあります.2次元だったらEigen
とかが使いやすそうです.
この情報を得る前に書いた,何も準備がいらないやり方についても念のため残しておきます.ただし静的に配列を確保する場合に限ります.
Fortran
Fortranでの多次元配列の取り扱いは本当に簡単です.なのでこっちは動的に確保する場合を示しましょう.
module resolution
integer:: imax,jmax
end module resolution
module variables
implicit none
real(8),dimension(:,:),allocatable:: A
end module variables
配列が直感的に扱える様子が分かるかと思います.Fortranは基本的に参照渡しするので,渡した変数の値を関数(subroutine)の中で変えれば関数の外に出ても値は変更されています.
program main
use resolution
use variables
implicit none
real(8),dimension(:,:),allocatable:: B
imax=2
jmax=3
allocate(A(imax,jmax))
call SetA
allocate(B(imax,jmax))
B(:,:) = A (:,:)
call changeB(B)
end program
subroutine setA()
use resolution
use variables
implicit none
integer::i,j
do j=1,jmax
do i=1,imax
A(i,j) = 1.0d0
enddo
enddo
return
end subroutine setA
subroutine changeB(B)
use resolution
implicit none
real(8),dimension(imax,jmax),intent(inout):: B
integer::i,j
do j=1,jmax
do i=1,imax
B(i,j) = 2.0d0
enddo
enddo
return
end subroutine changeB
C++
初めは静的に確保する場合についてのみ書かれていたのですが,,コメント欄で情報をもらいまして,動的に確保もできるようにしたので,順番に説明します.
静的に確保する場合
こちらでは今回は静的に配列を確保する話を限ります.
以下のヘッダファイルでネームスペース内で配列を宣言しています.配列の宣言の仕方そのものは簡単ですが,コンパイル時に配列の個数が決まっていなければなりません.ヘッダファイルには宣言はせども定義はしないのが原則なのですが,constexpr
はグローバルな定数のようなものでヘッダファイルに書いたほうが便利なようです(コメント欄にあるように最新の書き方ではinline constexprというのがあって,この目的に適しています).
template
についてはあとで説明します.
# ifndef _EX2
# define _EX2
constexpr int imax=2;
constexpr int jmax=3;
namespace variables{
extern double A[jmax][imax];
};
void setA();
template <size_t imax, size_t jmax>
void changeB(double (&B)[jmax][imax]);
# endif
プログラム本体の方はもう少し高級なことをしています.imax
, jmax
はヘッダファイルの中に書かれています.
たぶんみなさんつまずくのは配列の渡し方かと思います.template
は関数を修飾していろんな型やいろんな配列サイズに対して同じ処理をすることを可能にします.size_t imax
はimax
は配列の数だということを宣言しています.
この場合のB
は関数へ参照渡しされておりますので,関数内で値を変えると関数の外に出ても値が変わります.
include <iostream>
# include "ex2.hpp"
namespace variables{
double A[jmax][imax];
};
int main(){
using namespace std;
using namespace variables;
setA();
cout<<"A[0][0]="<<A[0][0]<<endl;
double B[jmax][imax];
for(int j=0;j<jmax;j++){
for(int i=0;i<imax;i++){
B[j][i] =A[j][i];
}
}
cout<<"B[0][0]="<<B[0][0]<<endl;
changeB(B);
cout<<"B[0][0]="<<B[0][0]<<endl;
}
void setA(){
using namespace variables;
for(int j=0;j<jmax;j++){
for(int i=0;i<imax;i++){
A[j][i] =1.0e0;
}
}
}
template <size_t imax, size_t jmax>
void changeB(double (&B)[jmax][imax]){
for(int j=0;j<jmax;j++){
for(int i=0;i<imax;i++){
B[j][i] =2.0e0;
}
}
}
動的に確保する場合
以下のように2次元配列的なものを表すclass
のmatrix
を作ります.違う名前にしたい時はサンプルコードの中に現れているmatrix
をその名前に全部置換してください.データそのもの(data_
)と配列の数(rows_,cols_
)を変数として持っています.
# ifndef _EX3
# define _EX3
# include <iostream>
# include <vector>
# include <cassert>
template <typename T>
class matrix{
private:
std::vector <T> data_;
std::size_t rows_ = 0;
std::size_t cols_ = 0;
public:
// constructor
matrix () = default;
matrix (std::size_t cols, std::size_t rows) : data_(rows*cols), rows_{rows}, cols_{cols} {
// if you need initialize. write something here.
}
matrix(const matrix&) = default;
matrix(matrix&&) = default;
// =
matrix& operator=(const matrix&) = default;
matrix& operator=(matrix&&) = default;
// access
T& operator()(const std::size_t j,const std::size_t i) noexcept {
assert(j < cols_ && i < rows_); // check whether the index is in the range or not
return this-> data_[j*rows_ + i];
}
const T& operator()(const std::size_t j,const std::size_t i)const noexcept {
assert(j < cols_ && i < rows_); // check whether the index is in the range or not
return data_[j*rows_ + i];
}
// size
std::size_t size() const noexcept { return data_.size(); }
std::size_t rows() const noexcept { return rows_; }
std::size_t cols() const noexcept { return cols_; }
// resize, allocate
void allocate(const std::size_t newcols, const std::size_t newrows){
data_.resize(newrows*newcols);
rows_= newrows;
cols_= newcols;
}
};
namespace variables{
extern matrix<double> A;
};
void setA();
void changeB(matrix<double> (&B));
void showB(const matrix<double> (&B));
# endif
上記のコードの一部を抜き出して,何が行われているのかを説明していきます.実際に使う場合,以下の変数を初期化するコンストラクタが大事です.コロンを使ったコンストラクタは単に変数を代入することを表します.ここでは初期化時にゼロを代入しないバージョンで書いていますが,最初にゼロを入れたい場合はここに何か書いてください.ゼロの代入の仕方がdata_[...]=0
みたいにする場合もdata[...].zero()
みたいに初期化する関数を作る場合もあるかと思って書いていません.
matrix (std::size_t cols, std::size_t rows) : data_(rows*cols), rows_{rows}, cols_{cols} {
// if you need initialize. write something here.
}
続く以下のプログラムは成分にアクセスするオペレータを設定しています.初見だと呪文に見えるかもしれません.ここではv(j,i)
とか書いたときにj
列のi
番目の変数にアクセスするため()
を定義しています.データはこのclass
の中では1次元配列としてアクセスをしています.これは2次元配列を使うと変数がメモリ上で連続にならばないことを避けるためです.配列の引数の後ろに書かれている変数(v(j,i)
ならi
)に対して連続アクセスになるよう書かれています.
noexcept
とかassert
の行がないと確保している配列の外にアクセスできてバグのもとになりますから,そういう時にエラーを吐いて終了するようにしてあります.
要素を関数に渡すときもあるでしょう.func(const double &v1, double &v2)
みたいなときを考えます.関数内でv1
は変更せず,v2
は変更します.このv2
に使う場合には参照渡ししないといけないので一つ目の定義のようはT&
で始まり,T
型の参照を返しています.一方で v1
みたいにconst
と指定されている場合に対応するためにconst T&
を返す場合も記述してあります.このときconst std::size_t i)
のあとのconst
はメンバ変数であるdata_
,rows
,cols
を変更しないという意味なので重要です.
T& operator()(const std::size_t j,const std::size_t i) noexcept {
assert(j < cols_ && i < rows_); // check whether the index is in the range or not
return this-> data_[j*rows_ + i];
}
const T& operator()(const std::size_t j,const std::size_t i)const noexcept {
assert(j < cols_ && i < rows_); // check whether the index is in the range or not
return data_[j*rows_ + i];
}
以下の部分ではFortranのallocate
みたいに,定義時ではないときに配列の大きさを決めるため,最初に個数ゼロで変数を確保しておいて,それのサイズを変更しています.
void allocate(int newcols, int newrows){
data_.resize(newrows*newcols);
rows_= newrows;
cols_= newcols;
}
このヘッダファイルの使用例のプログラムは以下です.定義と同時に配列の個数を決める場合は,matrix<double> B(jmax,imax);
みたいにします.定義は別にしてあとで,後で配列の個数を決める場合はA.allocate(jmax,imax);
です.関数の引数にはmatrix<double> (&B)
のように渡します.
# include <iostream>
# include "ex3.hpp"
constexpr int imax=3;
constexpr int jmax=2;
namespace variables{
matrix<double> A;
};
int main(){
using namespace std;
using namespace variables;
A.allocate(jmax,imax);
setA();
cout<<"A[0][0]="<<A(0,0)<<endl;
matrix<double> B(jmax,imax);
for(int j=0;j<jmax;j++){
for(int i=0;i<imax;i++){
B(j,i) = A(j,i);
}
}
cout<<"B[0][0]="<<B(0,0)<<endl;
changeB(B);
cout<<"B[0][0]="<<B(0,0)<<endl;
}
void setA(){
using namespace variables;
for(int j=0;j<jmax;j++){
for(int i=0;i<imax;i++){
A(j,i) =1.0e0;
}
}
}
void changeB(matrix<double> (&B)){
for(int j=0;j<jmax;j++){
for(int i=0;i<imax;i++){
B(j,i) =2.0e0;
}
}
}
何度も呼ぶ関数で1回目と2回目以降の動作を変える
1回目は初期化と処理,2回目以降は処理だけしたい場合はよくあります.そんなときどうしたら良いか.
Fortran
変数にsave
という属性をつけると関数を抜けてもその値を保持してくれるので,1回目に呼ばれたのかそうでないのか区別ができます.
以下のサブルーチンでは1回目の呼び出しでC
の配列を確保して初期値を代入し,2回目以降の呼び出しでは'B=A+C'を繰り返します.
subroutine savetempvar(A,B)
use resolution, only: imax,jmax
implicit none
real(8),dimension(imax,jmax),intent(in):: A
real(8),dimension(imax,jmax),intent(out):: B
real(8),dimension(:,:),allocatable,save:: C
integer::i,j
logical,save::is_inited
data is_inited /.false./
if(.not. is_inited)then
allocate(C(imax,jmax))
do j=1,jmax
do i=1,imax
C(i,j) = 1.0d0
enddo
enddo
is_inited = .true.
endif
do j=1,jmax
do i=1,imax
B(i,j) =A(i,j) +C(i,j)
enddo
enddo
return
end subroutine savetempvar
C++
変数にstatic
という属性をつけると関数を抜けてもその値を保持してくれるので,1回目に呼ばれたのかそうでないのか区別ができます.以下の関数では static bool is_inited
を使って1回目はif(is_inited == false )
の中が実行されます.
関数の中で配列も使っています.static matrix<double> C(jmax,imax);
はこの関数の中だけで使われる使い捨ての配列なのですが,関数が呼ばれる度に確保して関数が終わるたびに解放しないほうが良いので,static
にしています.こういう大きな配列を確保したり解放したりする処理が遅いことはよくありますし,メモリの確保と解放を繰り返すのはメモリリーク(長くプログラムを実行すると使用メモリが増えていってしまう現象)にも繋がるので避けています.
# include "ex3.hpp"
void savetempvar(const matrix<double> (&A), matrix<double> (&B)){
static matrix<double> C(jmax,imax);
static bool is_inited = false;
if(is_inited == false ){
for(int j=0;j<jmax;j++){
for(int i=0;i<imax;i++){
C(j,i) = 1.0;
}
}
is_inited= true;
}
for(int j=0;j<jmax;j++){
for(int i=0;i<imax;i++){
B(j,i) = A(j,i) + C(j,i);
}
}
}