Custom Query (986 matches)

Andy reported a bug in the do loop computing the Hessian eigenvectors in routine cgradient_lanczos.h. The loop counter is backwards.

In getting ROMS to run in OpenMP mode on my PC with Gfortran, I struck some anomalies in function my_numthreads (file ROMS/Utility/mp_routines.F). Here is the relevant code

#elif defined _OPENMP
      integer :: omp_get_max_threads
!!    integer :: omp_get_num_threads

      my_numthreads=omp_get_max_threads()
!!    my_numthreads=omp_get_num_threads()
...

Obviously someone has tried both omp_get_num_threads and omp_get_max_threads, and settled on the latter. But it didn't work for me with Gfortran, as ROMS always found just one thread and didn't pay any attention to environment variable OMP_NUM_THREADS. So, if all else fails, read the specification:

​http://www.openmp.org/drupal/mp-documents/spec25.pdf

I *think* I've worked it out. Function omp_get_max_threads should be called *before* you enter a parallel region to determine the number of threads available and omp_get_num_threads should be called *inside* a parallel region to determine the number of threads you have. When you call omp_get_max_threads inside a parallel region or omp_get_num_threads outside a parallel region, the results are undefined.

In ROMS the number of threads is determined once inside subroutine initialize_parallel, file ROMS/Modules/mod_parallel.F. So I think omp_get_max_threads is the correct function to use, but it is currently called inside a parallel region:

!$OMP PARALLEL SHARED(numthreads)
          numthreads=my_numthreads()
!$OMP END PARALLEL

I think this is wrong and it appears to give the wrong results with Gfortran. I removed the parallel declarations and it appears to work OK (tested on Gfortran and Intel Fortra). Modified file is attached.

By the way, the following F90 test program illustrates how these functions are supposed to work

program hello
  use omp_lib
  integer :: id, mthreads, nthreads
  mthreads = omp_get_max_threads()
  write (*,*) 'Maximum is', mthreads, 'threads'
  !$omp parallel private(id)
  id = omp_get_thread_num()
  write (*,*) 'Hello World from thread', id
  !$omp barrier
  if ( id .eq. 0 ) then
    nthreads = omp_get_num_threads()
    write (*,*) 'There are', nthreads, 'threads'
  end if
  !$omp end parallel
end program         

We need to have the capability of specifying the number of Hessian approximation eigenvectors to use during preconditioning. Currently, the conjugate gradient cgradient.h uses all the eigenpairs computed from the Lanczos algorithm.