Fortran で巨大配列の計算を行うときの注意

概要

• Fortran で巨大配列の計算を行うと、セグメンテーション違反することがある。
• マシンの stack size 制限に注意。

情報

• 実行環境
  • CentOS 6.4
  • fortran 90
  • fortran のソースコードは Cから呼び出し
• lapack は yum でインストールしたもの。普通に使うには古いかも。。

$ yum install lapack lapack-devel

実行

• C で、任意サイズの行列を作成し、fortran で計算させる。ここでは SVD。

$ ./a.out
Matrix size = ( 5, 6 )
Malloc memory.
Calculate SVD ...

• 巨大な配列になると、セグメンテーション違反する！！

$ ./a.out 1000 1000
Matrix size = ( 1000, 1000 )
Malloc memory.
セグメンテーション違反です (コアダンプ)

• fortran のメモリ領域は、マシンの stack size 分しか取らないらしい。
  • デフォルト ( 10240 ) のまま fortran のプログラムを動かしたのが原因！！
• ulimit コマンドを使って、マシンの stack size を調整すると正常に動く。

$ ulimit -s unlimited
$ ./a.out 1000 1000
Matrix size = ( 1000, 1000 )
Malloc memory.
Calculate SVD ...

サンプルソース

• main.c

#include <stdio.h>
#include <stdlib.h>

#define min(x,y) (((x) < (y)) ? (x) : (y))
#define max(x,y) (((x) < (y)) ? (y) : (x))

/**
 * Calculate SVD
 * matrix SRC, matrix U, array s, matrix Vt, matrix row, matrix col
 */
extern void svd_( double *, double *, double *, double *, int *, int * );

/**
 * Main
 */
int main( int argc, char *argv[] )
{
    int row = 5, col = 6;

    if ( argc >= 2 ) {
        row = atoi( argv[1] );
    }
    if ( argc >= 3 ) {
        col = atoi( argv[2] );
    }
    printf( "Matrix size = ( %ld, %ld )\n", row, col );

    size_t s_size = max( min( row, col ), 1 );
    size_t x, y;

    printf( "Malloc memory.\n" );
    double Src[ row * col ], U[ row * row ], s[ s_size ], Vt[ col * col ];

    int random_seed = 1;
    int threshold   = 60;
    for ( y = 0; y < row; y++ ) {
        for ( x = 0; x < col; x++ ) {
            Src[ x + y * col ] = ( rand() % 100 < threshold ) ? 0 : 1;
        }
    }

    printf( "Calculate SVD ... \n" );
    svd_( Src, U, s, Vt, &row, &col );

    return 0;
}

• svd.f90

subroutine svd( src, u, s, vt, row, col )
  implicit none

  integer row, col
  double precision src( row * col )
  double precision u( row * row ), s( max( row, col ) ), vt( col * col )

  double precision w_src( row, col )
  double precision w_u( row, row ), w_vt( col, col )

  integer lwork
  integer lda, ldu, ldvt

  double precision :: work( max( 3 * min( row, col ) + max( row, col ), 5 * min( row, col ) ) )

  integer i, j, info

  lwork = max( 3 * min( row, col ) + max( row, col ), 5 * min( row, col ) )
  lda   = max( row,   1 )
  ldu   = max( row,   1 )
  ldvt  = max(   1, col )

  ! Array => Matrix
  do i = 1, row
     do j = 1, col
        w_src( i, j ) = src( j + ( i - 1 ) * col )
     end do
  end do

  ! SVD
  call dgesvd( 'S', 'S', row, col, w_src, lda, s, w_u, ldu, w_vt, ldvt, work, lwork, info )

  ! Matrix => Array
  do i = 1, row
     do j = 1, row
        u( j + ( i - 1 ) * row ) = w_u( i, j )
     end do
  end do

  do i = 1, col
     do j = 1, col
        vt( j + ( i - 1 ) * col ) = w_vt( i, j )
     end do
  end do

  return
end subroutine svd

• Makefile

CC  = gcc
FTN = gfortran

a.out : svd.o main.o
	${FTN} -I /usr/local/lib64 -llapack $^

svd.o : svd.f90
	${FTN} -c $<

main.o : main.c
	${CC} -c $<


clean :
	rm -rf *.out
	rm -rf *.o