rework nthr logic replacing #431: no cap on user o.nthreads override; ifndef _OPENMP hard nthr=1 override; warnings at the finufft level only; binsort controlled by spopts.sort_threads matching docs

Author: ahbarnett

CHANGELOG

@@ -1,8 +1,9 @@
 List of features / changes made / release notes, in reverse chronological order.
 If not stated, FINUFFT is assumed (cuFINUFFT <=1.3 is listed separately).
 
-* CPU plan stage prevents now caps # threads at omp_get_max_threads (being 1
-  for single-thread build); warns if this cap was activated (PR 431)
+* CPU plan stage allows any # threads, warns if > omp_get_max_threads(); or
+  if single-threaded fixes nthr=1 and warns opts.nthreads>1 attempt.
+  Sort now respects spread_opts.sort_threads not nthreads. Supercedes PR 431.
 * new docs troubleshooting accuracy limitations due to condition number of the
   NUFFT problem.
 * new sanity check on nj and nk (<0 or too big); new err code, tester, doc.

src/finufft.cpp

@@ -591,15 +591,22 @@ int FINUFFT_MAKEPLAN(int type, int dim, BIGINT* n_modes, int iflag,
   p->fftSign = (iflag>=0) ? 1 : -1;         // clean up flag input
 
   // choose overall # threads...
-  int maxnthr = MY_OMP_GET_MAX_THREADS();
-  int nthr = maxnthr;                       // use as many as OMP gives us
+#ifdef _OPENMP
+  int ompmaxnthr = MY_OMP_GET_MAX_THREADS();
+  int nthr = ompmaxnthr;                    // default: use as many as OMP gives us
+  // (the above could be set, or suggested set, to 1 for small enough problems...)
   if (p->opts.nthreads>0) {
-    nthr = min(maxnthr,p->opts.nthreads);   // user override up to max avail
-    if (p->opts.nthreads > maxnthr)         // if no OMP, maxnthr=1
-      fprintf(stderr,"%s warning: user requested %d threads, but only %d threads available; enforcing nthreads=%d.\n",__func__,p->opts.nthreads,maxnthr,nthr);
+    nthr = p->opts.nthreads;                // user override, now without limit
+    if (p->opts.showwarn && (nthr > ompmaxnthr))
+      fprintf(stderr,"%s warning: using opts.nthreads=%d, more than the %d OpenMP claims available; note large nthreads can be slower.\n",__func__,nthr,ompmaxnthr);
   }
+#else
+  int nthr = 1;                             // always 1 thread (avoid segfault)
+  if (p->opts.nthreads>1)
+    fprintf(stderr,"%s warning: opts.nthreads=%d but library is single-threaded; ignoring!\n",__func__,p->opts.nthreads);
+#endif
   p->opts.nthreads = nthr;                  // store actual # thr planned for
-  // (this sets all downstream spread/interp, 1dkernel, and FFT thread counts)
+  // (this sets/limits all downstream spread/interp, 1dkernel, and FFT thread counts...)
 
   // choose batchSize for types 1,2 or 3... (uses int ceil(b/a)=1+(b-1)/a trick)
   if (p->opts.maxbatchsize==0) {            // logic to auto-set best batchsize

src/spreadinterp.cpp

@@ -267,14 +267,18 @@ int indexSort(BIGINT* sort_indices, BIGINT N1, BIGINT N2, BIGINT N3, BIGINT M,
   timer.start();                 // if needed, sort all the NU pts...
   int did_sort=0;
   int maxnthr = MY_OMP_GET_MAX_THREADS();
-  if (opts.nthreads>0)           // user override up to max avail
-    maxnthr = min(maxnthr,opts.nthreads);
-  
+  if (opts.sort_threads>0)           // user override, now without limit
+    maxnthr = opts.sort_threads;     // maxnthr = the max threads sorting could use
+  // (we don't print warning here, since: no showwarn in spread_opts, and finufft
+  // already warned about it. spreadinterp-only advanced users will miss a warning)
   if (opts.sort==1 || (opts.sort==2 && better_to_sort)) {
     // store a good permutation ordering of all NU pts (dim=1,2 or 3)
     int sort_debug = (opts.debug>=2);    // show timing output?
-    int sort_nthr = opts.sort_threads;   // choose # threads for sorting
-    if (sort_nthr==0)   // use auto choice: when N>>M, one thread is better!
+    int sort_nthr = opts.sort_threads;   // 0, or proposed max # threads for sorting
+#ifndef _OPENMP
+    sort_nthr = 1;                   // if single-threaded lib, override user
+#endif
+    if (sort_nthr==0)   // multithreaded auto choice: when N>>M, one thread is better!
       sort_nthr = (10*M>N) ? maxnthr : 1;      // heuristic
     if (sort_nthr==1)
       bin_sort_singlethread(sort_indices,M,kx,ky,kz,N1,N2,N3,opts.pirange,bin_size_x,bin_size_y,bin_size_z,sort_debug);
@@ -323,9 +327,12 @@ int spreadSorted(BIGINT* sort_indices,BIGINT N1, BIGINT N2, BIGINT N3,
   int ndims = ndims_from_Ns(N1,N2,N3);
   BIGINT N=N1*N2*N3;            // output array size
   int ns=opts.nspread;          // abbrev. for w, kernel width
-  int nthr = MY_OMP_GET_MAX_THREADS();  // # threads to use to spread
+  int nthr = MY_OMP_GET_MAX_THREADS();  // guess # threads to use to spread
   if (opts.nthreads>0)
-    nthr = min(nthr,opts.nthreads);     // user override up to max avail
+    nthr = opts.nthreads;       // user override, now without limit
+#ifndef _OPENMP
+  nthr = 1;                   // if single-threaded lib, override user
+#endif
   if (opts.debug)
     printf("\tspread %dD (M=%lld; N1=%lld,N2=%lld,N3=%lld; pir=%d), nthr=%d\n",ndims,(long long)M,(long long)N1,(long long)N2,(long long)N3,opts.pirange,nthr);
 
@@ -445,9 +452,12 @@ int interpSorted(BIGINT* sort_indices,BIGINT N1, BIGINT N2, BIGINT N3,
   int ndims = ndims_from_Ns(N1,N2,N3);
   int ns=opts.nspread;          // abbrev. for w, kernel width
   FLT ns2 = (FLT)ns/2;          // half spread width, used as stencil shift
-  int nthr = MY_OMP_GET_MAX_THREADS();   // # threads to use to interp
+  int nthr = MY_OMP_GET_MAX_THREADS();   // guess # threads to use to interp
   if (opts.nthreads>0)
-    nthr = min(nthr,opts.nthreads);      // user override up to max avail
+    nthr = opts.nthreads;       // user override, now without limit
+#ifndef _OPENMP
+  nthr = 1;                   // if single-threaded lib, override user
+#endif
   if (opts.debug)
     printf("\tinterp %dD (M=%lld; N1=%lld,N2=%lld,N3=%lld; pir=%d), nthr=%d\n",ndims,(long long)M,(long long)N1,(long long)N2,(long long)N3,opts.pirange,nthr);
 
@@ -1292,7 +1302,7 @@ void bin_sort_multithread(BIGINT *ret, BIGINT M, FLT *kx, FLT *ky, FLT *kz,
   nbins2 = isky ? N2/bin_size_y+1 : 1;
   nbins3 = iskz ? N3/bin_size_z+1 : 1;
   BIGINT nbins = nbins1*nbins2*nbins3;
-  if (nthr==0)
+  if (nthr==0)                      // should never happen in spreadinterp use
     fprintf(stderr,"[%s] nthr (%d) must be positive!\n",__func__,nthr);
   int nt = min(M,(BIGINT)nthr);     // handle case of less points than threads
   std::vector<BIGINT> brk(nt+1);    // list of start NU pt indices per thread