0

I am trying to parallelise a do-loop in Fortran 90 for the code that looks as follows:

!$OMP PARALLEL DEFAULT(private) &
!$OMP& SHARED(cmesh, S) ! user-defined types with allocatable members. 
!$OMP DO
do i=1,cmesh%nelements

   Ke=0.d0
   ! skip details. Ke is populated inside another subroutine.

   do iii=1,(size(cmesh%Edof(1,:))-1) 
       do jjj=1,(size(cmesh%Edof(1,:))-1)
      
           ! add matrix Ke( real*8 (6,6) ) to elements of sparse matrix S:
           !$OMP CRITICAL
           call add( S, &
                     cmesh%dof_2_dof_free( cmesh%Edof( i, 1+iii ) ), &
                     cmesh%dof_2_dof_free( cmesh%Edof( i, 1+jjj ) ), &
                     Ke( iii, jjj ) )
           !$OMP END CRITICAL

       enddo
   enddo
enddo !i=1,cmesh%nelements
!$OMP END DO
!$OMP END PARALLEL
                    
call PARDISO( pt, 1, 1, 1, 13, cmesh%ndof_free, &
              S%a,  ! real*8, allocatable(43656)
              S%ia, ! integer*4, allocatable(1691)
              S%ja, ! integer*4, allocatable(43656)
              perm, 1, iparm, 0, &
              p_Rp, p_sol, error )

This works without OpenMP. If I enable OpenMP, the call to PARDISO throws the error:

Thread 1 stopped in mkl_pds_lp64_pardiso_c with signal SIGSEGV (Segmentation fault).

Reason/Origin: address not mapped to object (attempt to access invalid address).

even with OMP_NUM_THREADS=1. What am I doing wrong?

I compile with ifort (IFORT) 19.0.5.281 20190815.

Vladimir F Героям слава
  • 57,977
  • 4
  • 76
  • 119
honey_badger
  • 472
  • 3
  • 10
  • 1
    Please use tag [tag:fortran] for all Fortran questions. We need a real minimal reproducible example and the complete output. – Vladimir F Героям слава Jan 03 '21 at 10:19
  • 2
    According to the documentation at https://pardiso-project.org/manual/manual.pdf you are missing an argument in the call to paradiso, namely dparam that comes after error. This could cause the above problem, but without a complete example it is impossible to tell. – Ian Bush Jan 03 '21 at 13:16
  • 2
    What's the declaration of `Ke`? Most Fortran compilers always put arrays declared `private` on the stack instead of placing the larger ones on the heap. This means the program will crash if the threads' stacks are too small. See here: https://stackoverflow.com/a/13266595/1374437 – Hristo Iliev Jan 04 '21 at 15:49
  • Thank you all for the suggestions. @IanBush, I use Intel's MKL Pardiso, which does not have that last argument: https://software.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-fortran/top/sparse-solver-routines/onemkl-pardiso-parallel-direct-sparse-solver-interface/pardiso.html#pardiso . But your comment made me go through all of the arguments again and it says about `pt` that "The entries must be set to zero prior to the first call to pardiso.". Introducing this fixes my problem. Oddly for me, it used to work without it without OpenMP. – honey_badger Jan 04 '21 at 20:08

0 Answers0