Making a new release that has the decoupled matrix treatment optional…

… rather than default

configure.ac

@@ -4,7 +4,7 @@ dnl
 dnl +------------------------
 dnl | Initialize package info
 dnl +------------------------
-AC_INIT([libmadam], [1.0.0], [github:hpc4cmb/libmadam])
+AC_INIT([libmadam], [1.0.1], [github:hpc4cmb/libmadam])
 AC_CONFIG_SRCDIR([Makefile.am])
 AM_INIT_AUTOMAKE([foreign])
 AC_CONFIG_HEADERS(config.h)

src/commonparam.f90

@@ -103,7 +103,7 @@ MODULE commonparam
   logical :: bin_subsets = .false.
   logical :: mcmode = .false., cached = .false.
 
-  real(dp) :: good_baseline_fraction=0 ! default acceps all baselines
+  real(dp) :: good_baseline_fraction=0 ! default accepts all baselines
   ! monte Carlo mode
   integer(idp) :: mc_increment=1e7, mc_loops=1, mc_id=0, rng_base=0
   logical :: incomplete_matrices = .false.
@@ -125,6 +125,7 @@ MODULE commonparam
 
   integer :: pixmode_map=2, pixmode_cross=2
   real(dp) :: pixlim_map=1e-6, pixlim_cross=1e-3
+  logical :: allow_decoupling = .false.
 
   real(dp) :: dnshort=-1
   integer :: nlong=-1, nshort=-1

src/inputparam.f90

@@ -236,6 +236,10 @@ SUBROUTINE read_line(line)
        read(value, *, iostat=ierr) pixlim_map
     case ('pixlim_cross')
        read(value, *, iostat=ierr) pixlim_cross
+    case('incomplete_matrices')
+       read(value, *, iostat=ierr) incomplete_matrices
+    case('allow_decoupling')
+       read(value, *, iostat=ierr) allow_decoupling
     case ('kfirst')
        read(value, *, iostat=ierr) kfirst
     case ('basis_func')
@@ -354,8 +358,6 @@ SUBROUTINE read_line(line)
        read(value, *, iostat=ierr) bin_subsets
     case('mcmode')
        read(value, *, iostat=ierr) mcmode
-    case('incomplete_matrices')
-       read(value, *, iostat=ierr) incomplete_matrices
 
     case default
        call abort_mpi('Unknown parameter: ' // trim(key))

src/matrix.f90

@@ -3,7 +3,7 @@ MODULE matrix
   ! Routines for handling of symmetric matrices
   !
 
-  use commonparam, only : incomplete_matrices
+  use commonparam, only : incomplete_matrices, allow_decoupling
 
   implicit none
   private
@@ -23,25 +23,29 @@ subroutine test_decoupling(cc, n, istart)
     integer, intent(out) :: istart
     integer :: i
 
-    !istart = 2
-    !do i = 1, n
-    !   if (cc(i, i) > 1e-30) then
-    !      if (any(abs(cc(1, 2:)) / cc(i, i) > 1e-10)) then
-    !         istart = 1
-    !      end if
-    !   end if
-    !end do
-    if (any(abs(cc(2:, 1)) > 1e-10)) then
+    if (.not. allow_decoupling) then
        istart = 1
     else
-       ! Decoupled case
        istart = 2
-    end if
+       do i = 1, n
+          if (cc(i, i) > 1e-30) then
+             if (any(abs(cc(1, 2:)) / cc(i, i) > 1e-10)) then
+                istart = 1
+             end if
+          end if
+       end do
+       if (any(abs(cc(2:, 1)) > 1e-10)) then
+          istart = 1
+       else
+          ! Decoupled case
+          istart = 2
+       end if
 
-    if (istart == 2 .and. cc(1, 1) > 1e-30) then
-       cc(1, 1) = 1 / cc(1, 1)
-       cc(2:, 1) = 0
-       cc(1, 2:) = 0
+       if (istart == 2 .and. cc(1, 1) > 1e-30) then
+          cc(1, 1) = 1 / cc(1, 1)
+          cc(2:, 1) = 0
+          cc(1, 2:) = 0
+       end if
     end if
 
   end subroutine test_decoupling