Merge pull request #12 from jacobwilliams/stack-overflow-fix

some changes to avoid stack overflow on windows for large arrays

src/numerical_differentiation_module.f90

@@ -222,6 +222,7 @@ module numerical_differentiation_module
 
         procedure,public :: failed !! to check if an exception was raised.
         procedure,public :: get_error_status !! the status of error condition
+        procedure,public :: set_dpert !! change the `dpert` values
 
         ! internal routines:
         procedure :: destroy_sparsity_pattern      !! destroy the sparsity pattern
@@ -1260,6 +1261,26 @@ subroutine set_numdiff_sparsity_bounds(me,xlow,xhigh)
     end subroutine set_numdiff_sparsity_bounds
 !*******************************************************************************
 
+!*******************************************************************************
+!>
+!  Change the `dpert` vector. Can be used after the class has been initialized
+!  to change the perturbation step sizes (e.g., after an iteration).
+
+    subroutine set_dpert(me,dpert)
+
+    class(numdiff_type),intent(inout) :: me
+    real(wp),dimension(:),intent(in) :: dpert !! perturbation vector for `x`
+
+    if (size(dpert)/=me%n) then
+        call me%raise_exception(29,'set_dpert',&
+            'incorrect size of dpert array')
+    else
+        me%dpert = abs(dpert) ! update
+    end if
+
+    end subroutine set_dpert
+!*******************************************************************************
+
 !*******************************************************************************
 !>
 !  Initialize a [[numdiff_type]] class. This must be called first.
@@ -1606,7 +1627,17 @@ subroutine columns_in_partition_group(me,igroup,n_cols,cols,nonzero_rows,indices
                                                 num_nonzero_elements_in_col
                 if (allocated(col_indices)) deallocate(col_indices)
                 allocate(col_indices(num_nonzero_elements_in_col))
-                col_indices = pack(me%indices,mask=me%icol==cols(i))
+                !col_indices = pack(me%indices,mask=me%icol==cols(i))
+                block
+                    integer :: j,n
+                    n = 0
+                    do j = 1, size(me%icol)
+                        if (me%icol(j)==cols(i)) then
+                            n = n + 1
+                            col_indices(n) = j
+                        end if
+                    end do
+                end block
                 if (allocated(nonzero_rows)) then
                     nonzero_rows = [nonzero_rows,me%irow(col_indices)]
                     indices = [indices,col_indices]
@@ -1652,7 +1683,10 @@ subroutine compute_indices(me)
     integer :: i !! counter
 
     allocate(me%indices(me%num_nonzero_elements))
-    me%indices = [(i,i=1,me%num_nonzero_elements)]
+    !me%indices = [(i,i=1,me%num_nonzero_elements)]
+    do i = 1, me%num_nonzero_elements
+        me%indices(i) = i
+    end do
 
     end subroutine compute_indices
 !*******************************************************************************
@@ -2780,7 +2814,18 @@ subroutine compute_jacobian_partitioned(me,x,dx,jac)
                         num_nonzero_elements_in_col = count(me%sparsity%icol==cols(i))
                         if (allocated(col_indices)) deallocate(col_indices)
                         allocate(col_indices(num_nonzero_elements_in_col))
-                        col_indices = pack(me%sparsity%indices,mask=me%sparsity%icol==cols(i))
+                        ! col_indices = pack(me%sparsity%indices,mask=me%sparsity%icol==cols(i))
+                        block
+                            integer :: j,n
+                            n = 0
+                            do j = 1, size(me%sparsity%icol)
+                                if (me%sparsity%icol(j)==cols(i)) then
+                                    n = n + 1
+                                    col_indices(n) = j
+                                end if
+                            end do
+                        end block
+
                         df(me%sparsity%irow(col_indices)) = df(me%sparsity%irow(col_indices)) / &
                                                             (me%meth(1)%df_den_factor*dx(cols(i)))
                     end do