Merge pull request #8 from hpc4cmb/local_pointing_periods

Local pointing periods

src/commonparam.f90

@@ -28,9 +28,6 @@ MODULE commonparam
 
   ! TOAST additions -RK
   integer(dp), allocatable :: chunk_offsets(:), chunk_sizes(:)
-  integer(i4b) :: n_chunkset, n_chunkset_tot, n_channel, n_chunk, n_chunk_tot
-  integer(i4b) :: n_chunk_max, n_sky_nonzero
-  integer(i4b) :: first_chunk, last_chunk ! starts from 1
   logical, allocatable :: use_toast_channel(:)
   character(len=256) :: telescope_name, sky_name
   logical :: flag_by_horn=.false., force_pol=.false.
@@ -134,7 +131,6 @@ MODULE commonparam
   ! Input files
   character(len=SLEN) :: file_param='', file_simulation='', &
        file_inmask='', file_spectrum='', file_gap='', &
-       file_pntperiod='', file_objectsize='', &
        file_fpdb_supplement=''
 
   ! Output files
@@ -177,12 +173,11 @@ MODULE commonparam
   integer(i4b), allocatable :: baselines_short_start(:)
   integer(i4b), allocatable :: baselines_short_stop(:)
   real(i8b), allocatable :: baselines_short_time(:)
-  integer(i4b), allocatable :: base_pntid_short(:)
 
   ! Number of pointing periods and their duration as a number of samples
-  integer(i4b) :: nopntperiods = -1
-  integer(i8b), allocatable :: pntperiods(:)
-  integer(i4b), allocatable :: pntperiod_id(:)
+  integer(i4b) :: ninterval = -1, ninterval_tot = -1
+  integer(i8b), allocatable :: intervals(:)
+  integer(i4b), allocatable :: interval_id(:)
 
   ! Baselines per pointing period (only used in case use_pntperiods=T)
   integer(i4b), allocatable :: noba_short_pp(:)

src/covmat.F90

@@ -47,7 +47,7 @@ MODULE covmat
   subroutine write_covmat(outroot)
     character(len=*) :: outroot
 
-    integer(i4b) :: idet, ichunk, ipsd, ipsd_det, pid, ierr
+    integer(i4b) :: idet, ichunk, ipsd, ipsd_det, ierr
     integer(i4b) :: imap, ipix, ibase, jmap, jpix, jbase, k, i, j, ip, kstart, noba
 
     real(dp) :: detweight, mm, ptf1
@@ -89,18 +89,16 @@ subroutine write_covmat(outroot)
 
        call tic(555)
 
-       loop_chunk : do ichunk = first_chunk, last_chunk ! global indices
+       loop_chunk : do ichunk = 1, ninterval ! global indices
           noba = noba_short_pp(ichunk) ! baselines on this pointing period
-          kstart = sum(noba_short_pp(first_chunk:ichunk-1)) ! first baseline, local index
+          kstart = sum(noba_short_pp(1:ichunk-1)) ! first baseline, local index
           ipsd = psd_index( idet, baselines_short_time(kstart+1) )
           ipsd_det = psd_index_det( idet, baselines_short_time(kstart+1) )
           if ( ipsd_det < 0 ) cycle loop_chunk
           ! ALWAYS use the optimal weights, even if the map has some other weighting scheme
           detweight = 1 / detectors(idet)%sigmas(ipsd_det)**2
           if (detweight == 0) cycle loop_chunk
 
-          pid = pntperiod_id(ichunk)
-
           call get_ptf( idet, kstart, noba, detweight )
 
           call get_hitmap( idet, kstart, noba )

src/inputparam.f90

@@ -338,10 +338,6 @@ SUBROUTINE read_line(line)
        file_spectrum = trim(value)
     case ('file_gap')
        file_gap = trim(value)
-    case ('file_pntperiod')
-       file_pntperiod = trim(value)
-    case ('file_objectsize')
-       file_objectsize = trim(value)
     case ('file_fpdb_supplement')
        file_fpdb_supplement = trim(value)
     case('binary_output')

src/noise_routines.f90

@@ -627,9 +627,9 @@ SUBROUTINE convolve_pp(y, x, fc, mode)
     !$OMP DO SCHEDULE(DYNAMIC,1)
     do idet = 1, nodetectors
 
-       do ichunk = first_chunk, last_chunk
+       do ichunk = 1, ninterval
           noba = noba_short_pp(ichunk)
-          kstart = sum(noba_short_pp(first_chunk:ichunk-1))
+          kstart = sum(noba_short_pp(1:ichunk-1))
           ipsd = psd_index(idet, baselines_short_time(kstart+1))
 
           if (ipsd == -1) then
@@ -704,7 +704,7 @@ END SUBROUTINE convolve_pp
   SUBROUTINE construct_preconditioner(nna)
 
     real(dp), intent(in)  :: nna(noba_short_max, nodetectors)
-    integer :: i, j, k, kstart, n, noba, pid, idet, ichunk, ipsd, ierr, nbandmin
+    integer :: i, j, k, kstart, n, noba, idet, ichunk, ipsd, ierr, nbandmin
     real(dp), allocatable :: invcov(:, :), blockm(:, :)
     logical :: neg
 
@@ -726,9 +726,9 @@ SUBROUTINE construct_preconditioner(nna)
 
        prec_diag = 0
        do idet = 1, nodetectors
-          do ichunk = first_chunk, last_chunk
+          do ichunk = 1, ninterval
              noba = noba_short_pp(ichunk)
-             kstart = sum(noba_short_pp(first_chunk:ichunk-1))
+             kstart = sum(noba_short_pp(1:ichunk-1))
              ipsd = psd_index_det(idet, baselines_short_time(kstart+1))
              if (ipsd < 0) then
                 ! No PSD available
@@ -774,39 +774,37 @@ SUBROUTINE construct_preconditioner(nna)
 
     !$OMP PARALLEL DO IF (nodetectors >= nthreads) &
     !$OMP     DEFAULT(NONE) &
-    !$OMP     SHARED(nodetectors, first_chunk, last_chunk, noba_short_pp, &
-    !$OMP         baselines_short_time, pntperiod_id, nband, invcov, id, &
+    !$OMP     SHARED(nodetectors, ninterval, noba_short_pp, &
+    !$OMP         baselines_short_time, nband, invcov, id, &
     !$OMP         sampletime, nof, baselines_short_start, &
     !$OMP         baselines_short_stop, bandprec, nna, prec_diag, &
     !$OMP         detectors, nthreads) &
-    !$OMP     PRIVATE(idet, ichunk, noba, kstart, ipsd, pid, blockm, &
+    !$OMP     PRIVATE(idet, ichunk, noba, kstart, ipsd, blockm, &
     !$OMP         i, j, k, n, neg, ierr, nbandmin)
     do idet = 1,nodetectors
 
        !$OMP PARALLEL DO IF (nodetectors < nthreads) &
        !$OMP     DEFAULT(NONE) &
-       !$OMP     SHARED(idet, nodetectors, first_chunk, last_chunk, &
-       !$OMP         noba_short_pp, baselines_short_time, pntperiod_id, nband, &
+       !$OMP     SHARED(idet, nodetectors, ninterval, &
+       !$OMP         noba_short_pp, baselines_short_time, nband, &
        !$OMP         invcov, id, sampletime, nof, baselines_short_start, &
        !$OMP         baselines_short_stop, bandprec, nna, prec_diag, &
        !$OMP         detectors) &
-       !$OMP     PRIVATE(ichunk, noba, kstart, ipsd, pid, blockm, &
+       !$OMP     PRIVATE(ichunk, noba, kstart, ipsd, blockm, &
        !$OMP         i, j, k, n, neg, ierr, nbandmin)
-       do ichunk = first_chunk, last_chunk
+       do ichunk = 1, ninterval
           noba = noba_short_pp(ichunk)
-          kstart = sum(noba_short_pp(first_chunk:ichunk-1))
+          kstart = sum(noba_short_pp(1:ichunk-1))
           ipsd = psd_index(idet, baselines_short_time(kstart+1))
 
-          pid = pntperiod_id(ichunk)
-
           if (noba < nband) then
              ! This chunk (pointing period) has fewer baselines than
              ! the width of the preconditioner.  This should be a very
              ! infrequent issue unless the user has chosen a too wide
              ! preconditioner.
              if (idet == 1 .and. omp_get_thread_num() == 0) then
                 write(*,*) id, ' : WARNING : noba < nband : ', noba, ' < ', &
-                     nband, ' : ichunk = ', ichunk, ', pid = ', pid, &
+                     nband, ' : ichunk = ', ichunk, &
                      ', tstart = ', &
                      sampletime(baselines_short_start(kstart+1)), &
                      ', length = ', &
@@ -920,27 +918,27 @@ SUBROUTINE preconditioning_band(z, r)
     else
        !$OMP PARALLEL DO IF (nodetectors >= nthreads) &
        !$OMP     DEFAULT(NONE) &
-       !$OMP     SHARED(noba_short_pp, first_chunk, nband, prec_diag, &
+       !$OMP     SHARED(noba_short_pp, ninterval, nband, prec_diag, &
        !$OMP         bandprec, baselines_short_time, pfx, pxx, z, r, &
-       !$OMP         nodetectors, last_chunk, id, nof, nofilter, notail, &
+       !$OMP         nodetectors, id, nof, nofilter, notail, &
        !$OMP         fcov, nshort, detectors, nthreads) &
        !$OMP     PRIVATE(idet, ichunk, kstart, noba, j, k, &
        !$OMP         ierr, ipsd, x0, m, no, xx, fx, id_thread, ipsddet, &
        !$OMP         nbandmin)
        do idet = 1, nodetectors
           !$OMP PARALLEL DO IF (nodetectors < nthreads) &
           !$OMP     DEFAULT(NONE) &
-          !$OMP     SHARED(idet, noba_short_pp, first_chunk, nband, prec_diag, &
+          !$OMP     SHARED(idet, noba_short_pp, ninterval, nband, prec_diag, &
           !$OMP         bandprec, baselines_short_time, pfx, pxx, z, &
-          !$OMP         r, nodetectors, last_chunk, id, nof, nofilter, notail, &
+          !$OMP         r, nodetectors, id, nof, nofilter, notail, &
           !$OMP         fcov, nshort, detectors) &
           !$OMP     PRIVATE(ichunk, kstart, noba, j, k, ierr, &
           !$OMP         ipsd, x0, m, no, xx, fx, id_thread, ipsddet, nbandmin)
-          do ichunk = first_chunk, last_chunk
+          do ichunk = 1, ninterval
 
              noba = noba_short_pp(ichunk)
              if (noba == 0) cycle
-             kstart = sum(noba_short_pp(first_chunk:ichunk-1))
+             kstart = sum(noba_short_pp(1:ichunk-1))
 
              if (noba < nband) then
                 nbandmin = noba
@@ -981,7 +979,7 @@ SUBROUTINE preconditioning_band(z, r)
                 call c_f_pointer(pxx(id_thread), xx, [nof])
                 !$OMP END CRITICAL
 
-                !Remove the PID average first.
+                !Remove the interval average first.
                 !This ensures that solution is independent of a constant
                 ! offset/pointing period
                 x0 = sum(r(kstart+1:kstart+noba, idet)) / noba