いつもFotranを使っています．ただ今回共同研究者がC++でコードを書いてきたので，C++でコードを書かないといけない事案が発生しました．その時，コーディングのお約束が結構あって苦労したので次回から簡単に解決できるようにメモを残します．複数ファイルがあるような状態を想定しています．

以下でFortranはifort -extend-sourceでコンパイルしています．C++はg++ -std=c++0xでコンパイルしています．

namespace

なんでもかんでもグローバル変数にしないのがモダンなプログラミングの基本です．そのためにはnamespaceが役に立ちます．

Fortran

例えば以下のようにFotranで光速をcgs単位系で定義しているmoduleがあるとします（僕は長くとった固定幅で書く流儀です）．

modules.f

      module cgs_units
      implicit none
      real(8),parameter:: c=2.997924562d+10
      end module cgs_units

プログラム上では以下のように呼びますね．

main.f

      program main
      use cgs_units
      implicit none
      real(8):: kousoku
      kousoku = c
      end program

C++

これをC++で書くとどうなるか．ヘッダーファイルでは宣言（メモリを確保しない）だけして，定義（実際にメモリを確保する）は全部のプログラムで一回だけにしないといけません．なのでヘッダファイルではexternを使って外部のファイルに変数が定義されていることを明示します．あとヘッダファイルは何回も読み込んではいけないので読み込みが一回きりになるように#ifndef _CGS_UNITSを使っています．

units.hpp

# ifndef _CGS_UNITS
# define _CGS_UNITS
namespace cgs_units{
  extern const double c;
};
# endif

プログラム上では以下のように呼びますね．まずどこかでcを定義しておきます．Fortranでは数字の桁数をd10とか書いてましたが，e10に直しています．

main.cpp

# 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での多次元配列の取り扱いは本当に簡単です．なのでこっちは動的に確保する場合を示しましょう．

modules.f

      module resolution
      integer:: imax,jmax
      end module resolution

      module variables
      implicit none
      real(8),dimension(:,:),allocatable:: A
      end module variables

配列が直感的に扱える様子が分かるかと思います．Fortranは基本的に参照渡しするので，渡した変数の値を関数（subroutine）の中で変えれば関数の外に出ても値は変更されています．

main.f

      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についてはあとで説明します．

ex2.hpp

# 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は関数へ参照渡しされておりますので，関数内で値を変えると関数の外に出ても値が変わります．

ex2.cpp

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_)を変数として持っています．

ex3.hpp

# 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)のように渡します．

ex3.cpp

# 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);
  }
  }

}

参考資料

Physics Lab. 2009 計算機班 :C++ 初歩