matrix - Fortran Seg Fault when assigning Matrices -

[update] code , few sentences changed reflect realization explained in second comment. code should compile line below, however, have older gfortran , may not seeing errors might.

gfortran blu_implementation_copy.f90 -o blu_implementation_copy.x 

i'm getting incredibly inconsistent seg fault error @ runtime fortran 90.

the overall objective of code factor randomly generated, complex, symmetric, diagonally-dominant matrix in blocks can parallelized. stripped down version (except relevant functions) can found @ end.

i'm running error when breaking original matrix blocks (tiles) in function mat2tiles. line of code keeps failing:

o(i,j)=cell(x(a:min(i*nblock,m), b:min(j*nblock, n)))

this assigning block of input matrix x matrix o @ index i,j. accomplishes because each index of o matrix nblockxnblock through derived type:

type cell complex*16 :: block(nblock,nblock) end type cell

type(cell) :: o(m/nblock+rem,n/nblock+rem) 

for matrix sizes , tile sizes, works perfectly. larger sizes, error starts appear in tile sizes close matrix size , reach @ half size of matrix. example, smallest instance of happening have found 20x20 matrix tile size of 18. doesn't happen consistently. if try out smaller matrices run , build size, runs. if larger size seg faults , run 20x20 18, seg faults. smallest consistent parameters i've found have been 25x25 tile size of 23. yet size, if print out block fails after o(i,j) assignment line:

print*, o(2,1)%block

it runs through whole thing without problem. taking out print statement makes seg fault again. finally, annoying 1 [fixed - refer second comment me], if 2000x2000 matrix tile size of 1000 or more, seg fault after getting function (which means it's occurring during variable allocations). leads me believe issue may stem how matrix allocated, since i'm using derived type. when try diagnose size , contents of matrix, seems normal.

program blu_implementation     implicit none     integer :: rem     integer, parameter :: m=20, n=20, nblock=19     real*8 :: start, finish     complex*16 :: a(m,n)      type cell     complex*16 :: block(nblock,nblock)     end type cell  !determines if cell matrix doesn't need bigger     rem=1     if (modulo(m,nblock)==0)         rem=0     endif      call cpu_time(start)     call functioncalling(a, m, n, nblock, rem)     call cpu_time(finish)     print*, 'overall time:', finish-start, 'seconds'  contains !================================================================================================================== subroutine functioncalling(a, m, n, nblock, rem)     implicit none     integer :: ipiv, info, m, n, nblock, rem     real*8 :: start, finish     complex*16 :: a(m,n)      type(cell) :: c(m/nblock+rem,n/nblock+rem)      call cpu_time(start)      a= cspdmatrixfill(a,m,n)      call cpu_time(finish)     print*, 'matrix fill time:', finish-start, 'seconds'       call cpu_time(start)      c= mat2tiles(a,m,n,nblock,rem)       call cpu_time(finish)     print*, 'tiling time:', finish-start, 'seconds'  end subroutine  !=================================================================================================================== ! generates complex, symmetric, positive-definite matrix (based off of another's code)  function cspdmatrixfill(a, m, n) result (matrix)  ! initialization         implicit none       integer :: i, j       integer :: m, n       real*8 :: x, xi       complex*16 :: a(m, n), matrix(m, n), eye(m, n), mt(n, m)        eye=0       forall(j=1:m) eye(j,j)=1  ! execution        call random_seed        i=1, m         j=1, n           call random_number (x)             call random_number(xi)             matrix(i,j) = cmplx(x,xi)        end       end    ! construct symmetric matrix ( o(n^2) )       call mtranspose(matrix, m, n)       matrix = matrix+mt    ! make positive definite (diagonally dominant)         matrix = matrix + n*eye  end function cspdmatrixfill  !====================================================================================================== subroutine mtranspose(a, i, j) ! takes matrix , 2 parameters used make matrix: a(i,j) ! returns matrix switched indices: a(j,i)       implicit none        integer :: i, j       complex*16 :: a(i,j), mt(j,i)        mt=a(j,i)       return end subroutine mtranspose  !======================================================================================================= !mat2tiles - breaks array cell array of adjacent sub-arrays of equal sizes ! ! o=mat2tiles(x,m,n,nblock) ! !will produce cell array o containing adjacent chunks of array x(m,n) !with each chunk of dimensions nblockxnblock. if nblock !not divide evenly size(x,i), chunks @ upper boundary of !x along dimension have bogus values not affect factorization !in places matrix doesn't occupy. (according older versions. might have changed edits) ! function mat2tiles(x,m,n,nblock,rem) result(o)  ! initialization         implicit none       integer :: a, b, i, j, m, n, nblock, rem       complex*16 :: x(m,n)        type(cell) :: o(m/nblock+rem,n/nblock+rem)  ! diagnostic print statements         print*,size(o(1,1)%block), size(o(1,2)%block), size(o(2,1)%block), size(o(2,2)%block)         print*, 'got start'  ! turn matrix x cell matrix o        j=1, n/nblock+rem         if (j==1)           b=j         else            b=b+nblock         endif                i=1, m/nblock+rem           if (i==1)             a=i           else              a=a+nblock           endif ! diagnostic print statement           print*, 'writing o: i:', i, 'j:', j, 'i*nblock:', i*nblock, 'j*nblock:', j*nblock, 'min of i:', min(i*nblock, m), &                     'min of j:', min(j*nblock, n), 'a:', a, 'b:', b            o(i,j)=cell(x(a:min(i*nblock,m), b:min(j*nblock, n)))          enddo       enddo ! diagnostic print statement       print*, 'got end'       return  end function mat2tiles  !================================================================================================== end program 

after narrowing down issue using numbers vs. variables of same numbers, discovered fortran doesn't assigning matrix matrix of different dimensions if fit inside. weird because smaller values of m, n, , nblock got around issue perfectly.

the solution defining o(i,j)%block(1:nblock,1:nblock)=x(dim1:dim2,etc) instead of o(i,j)=cell(x(dim1:dim2,etc) i=1 , j=1 , having modified instance of each cases of i , j.

correct code (for cases of matrix size , tile size) looks like:

if (i==m/nblock+rem .and. j==n/nblock+rem .and. rem==1)      o(i,j)%block(1:m-(i-1)*nblock,1:n-(j-1)*nblock)=x(a:min(i*nblock,m), b:min(j*nblock, n)) else if (i==m/nblock+rem .and. j/=n/nblock+rem .and. rem==1)     o(i,j)%block(1:m-(i-1)*nblock,1:nblock)=x(a:min(i*nblock,m), b:min(j*nblock, n)) else if (i/=m/nblock+rem .and. j==n/nblock+rem .and. rem==1)     o(i,j)%block(1:nblock,1:n-(j-1)*nblock)=x(a:min(i*nblock,m), b:min(j*nblock, n)) else     o(i,j)%block(1:nblock,1:nblock)=x(a:min(i*nblock,m), b:min(j*nblock, n)) end if 

with correction, code runs on old , new versions of gfortran. interesting note, however, how mac os x version of gfortran never ran segfault, linux version.