use allocatable to move large local variables to heap

Author: helske

ChangeLog

@@ -9,6 +9,8 @@ Changes from version 1.5.1 to 1.6.0
      `n.ahead` to define constant (series-specific) `u` for future predictions, 
      thus avoiding the need for the `newdata` if model is otherwise 
      time-invariant.
+   * Adjustments to Fortran codes to deal with optimizations of the flang 
+     compiler which lead to importance sampling to crash due to stack overflow.
 
 Changes from version 1.5.0 to 1.5.1
    * Added explicit alias for KFAS-package due to changes in roxygen2.

src/covmeanw.f90

@@ -29,27 +29,31 @@ subroutine covmeanwprotect(x,w,m,n,k,meanx,covx)
 
     implicit none
     integer, intent(in) :: m, n, k
-    integer :: t,i
+    integer :: t,i,j,l
     double precision, intent(in), dimension(m,n,k) :: x
     double precision, intent(in), dimension(k) :: w
     double precision, intent(inout), dimension(m,n) :: meanx
     double precision, intent(inout), dimension(m,m,n) :: covx
-    double precision, dimension(m,n,k) :: x2
+    double precision, dimension(:,:), allocatable :: x2
 
     external dgemm
-
-
-    do i = 1, k
-        meanx = meanx + x(:,:,i)*w(i)
-    end do
-    do i = 1, k
-        x2(:,:,i) = sqrt(w(i))*(x(:,:,i) - meanx)
+     do i = 1, k
+        do j = 1, m
+          do l = 1,n
+              meanx(j, l) = meanx(j,l)+x(j,l,i)*w(i)
+          end do
+        end do
     end do
-
+    allocate(x2(m,k))
     do t = 1, n
-        call dgemm('n','t',m,m,k,1.0d0,x2(:,t,:),m,x2(:,t,:),m,0.0d0,covx(:,:,t),m)
+        do i = 1, k
+          do j = 1, m
+              x2(j, i) = sqrt(w(i))*(x(j,t,i) - meanx(j,t))
+          end do
+        end do
+        call dgemm('n','t',m,m,k,1.0d0,x2,m,x2,m,0.0d0,covx(:,:,t),m)
     end do
-
+    deallocate(x2)
 end subroutine covmeanwprotect
 
 subroutine varmeanw(x,w,m,n,k,meanx,varx,var)

src/isample.f90

@@ -30,13 +30,13 @@ subroutine isample(yt, ymiss, timevar, zt, tt, rtv, qt, a1, p1,p1inf, u, dist, &
     double precision, intent(inout), dimension(n,p) :: theta
     double precision, dimension(p,p,n) :: ht
     double precision, intent(inout), dimension(simdim,n,3 * nsim * antithetics + nsim) :: sim
-    double precision, dimension(p,(3 * nsim * antithetics + nsim)*(5-simwhat)) :: tsim
     double precision, dimension(n,p) :: ytilde
     double precision, dimension(n) :: tmp
-    double precision, dimension(3 * nsim * antithetics + nsim) :: w
+    double precision, intent(inout), dimension(3 * nsim * antithetics + nsim) :: w
     double precision :: diff
     double precision, external :: ddot
     double precision :: lik
+    double precision, dimension(:,:), allocatable :: tsim
 
     external approx, simgaussian
 
@@ -112,6 +112,8 @@ subroutine isample(yt, ymiss, timevar, zt, tt, rtv, qt, a1, p1,p1inf, u, dist, &
         end do
 
     else
+    
+        allocate(tsim(p,(3 * nsim * antithetics + nsim)*(5-simwhat)))
         do j=1,p
             select case(dist(j))
                 case(2)    !poisson
@@ -165,9 +167,7 @@ subroutine isample(yt, ymiss, timevar, zt, tt, rtv, qt, a1, p1,p1inf, u, dist, &
                     end do
             end select
         end do
-
-
+        deallocate(tsim)
     end if
-
 
 end subroutine isample

src/isamplefilter.f90

@@ -30,21 +30,21 @@ subroutine isamplefilter(yt, ymiss, timevar, zt, tt, rtv, qt, a1, p1,p1inf, u, d
     double precision, intent(inout), dimension(n,p) :: theta
     double precision, dimension(p,p,n) :: ht
     double precision, intent(inout), dimension(simdim,n,3 * nsim * antithetics + nsim) :: sim
-    double precision, dimension(p,(3 * nsim * antithetics + nsim)*(5-simwhat)) :: tsim
+    double precision, intent(inout), dimension(n,3 * nsim * antithetics + nsim) :: w
     double precision, dimension(n,p) :: ytilde
     double precision, dimension(n) :: tmp
-    double precision, dimension(n,3 * nsim * antithetics + nsim) :: w
     double precision, external :: ddot
 
     integer, dimension(n,p) :: ymiss2
     double precision, dimension(p,n,nsim) :: epsplus2
     double precision, dimension(r,n,nsim) :: etaplus2
     double precision, dimension(m,nsim) :: aplus12
-    double precision, dimension(simdim,n,3 * nsim * antithetics + nsim) :: sim2
     double precision :: diff
     double precision :: lik
-
+    double precision, dimension(:,:), allocatable :: tsim
+    double precision, dimension(:,:,:), allocatable :: sim2
     external approx, simgaussian
+    allocate(sim2(simdim,n,3 * nsim * antithetics + nsim))
 
     ht=0.0d0
     ytilde=0.0d0
@@ -155,6 +155,7 @@ subroutine isamplefilter(yt, ymiss, timevar, zt, tt, rtv, qt, a1, p1,p1inf, u, d
             end do
 
         else
+            allocate(tsim(p,(3 * nsim * antithetics + nsim)*(5-simwhat)))
             do j=1,p
                 select case(dist(j))
                     case(2)    !poisson
@@ -208,9 +209,10 @@ subroutine isamplefilter(yt, ymiss, timevar, zt, tt, rtv, qt, a1, p1,p1inf, u, d
                         end do
                 end select
             end do
+            deallocate(tsim)
         end if
         sim(:,i+1,:) = sim2(:,i+1,:)
     end do
     maxiter=maxitermax
-   
+    deallocate(sim2)
 end subroutine isamplefilter

src/ngfilter.f90

@@ -35,13 +35,16 @@ subroutine ngfilter(yt, ymiss, timevar, zt, tt, rtv, qt, a1, p1,p1inf, u, theta,
     double precision, intent(inout), dimension((p-1)*smoothy+1,(p-1)*smoothy+1,(n-1)*smoothy+1) :: yvar
     double precision, intent(inout), dimension((p-1)*smooths+1,(n-1)*smooths+1) :: thetahat
     double precision, intent(inout), dimension((p-1)*smooths+1,(p-1)*smooths+1,(n-1)*smooths+1) :: thetavar
-    double precision, dimension(smootha*m+(1-smootha)*p,n,4*nsim) :: sim
-    double precision, dimension(smootha*p+(1-smootha),n,4*nsim*smootha+(1-smootha)) :: osim
-    double precision, dimension(n,4*nsim) :: w
     double precision, dimension(p,m) :: pm
-    double precision, external :: ddot
+    double precision, dimension(:,:,:), allocatable :: sim
+    double precision, dimension(:,:,:), allocatable :: osim
+    double precision, dimension(:,:), allocatable :: w
 
+    double precision, external :: ddot
     external isamplefilter, covmeanwprotect, dgemv, dsymm, dgemm, covmeanw
+    
+    allocate(w(n,4*nsim))
+    allocate(sim(smootha*m+(1-smootha)*p,n,4*nsim))
 
     if(smootha.EQ.1) then
 
@@ -67,6 +70,7 @@ subroutine ngfilter(yt, ymiss, timevar, zt, tt, rtv, qt, a1, p1,p1inf, u, theta,
         end if
 
         if(smoothy.EQ.1) then
+            allocate(osim(smootha*p+(1-smootha),n,4*nsim*smootha+(1-smootha)))
             do t = 1, n
                 do j = 1,p
                     call dgemv('t',m,4*nsim,1.0d0,sim(:,t,:),m,zt(j,:,(t-1)*timevar(1)+1),1,0.0d0,osim(j,t,:),1)
@@ -90,6 +94,7 @@ subroutine ngfilter(yt, ymiss, timevar, zt, tt, rtv, qt, a1, p1,p1inf, u, theta,
             do t=1,n
                 call covmeanw(osim(:,t,:),w(t,:),p,1,4*nsim,yhat(:,t),yvar(:,:,t))
             end do
+            deallocate(osim)
         end if
     else
         call isamplefilter(yt, ymiss, timevar, zt, tt, rtv, qt, a1, p1,p1inf, u, dist, &
@@ -129,6 +134,7 @@ subroutine ngfilter(yt, ymiss, timevar, zt, tt, rtv, qt, a1, p1,p1inf, u, theta,
                 call covmeanw(sim(:,t,:),w(t,:),p,1,4*nsim,yhat(:,t),yvar(:,:,t))
             end do
         end if
-
     end if
+    deallocate(sim)
+    deallocate(w)
 end subroutine ngfilter

src/ngloglik.f90

@@ -31,10 +31,10 @@ subroutine ngloglik(yt, ymiss, timevar, zt, tt, rtv, qt, a1, p1,p1inf, p,m, &
     double precision, intent(inout), dimension(p,n,nsim) :: epsplus
     double precision, intent(inout), dimension(r,n,nsim) :: etaplus
     double precision, intent(inout) :: lik
-    double precision, dimension(p,n,nsim2) :: tsim
     double precision, dimension(n) :: tmp
-    double precision, dimension(nsim2) :: w
     double precision, intent(inout) :: diff
+    double precision, dimension(:), allocatable :: w
+    double precision, dimension(:,:,:), allocatable :: tsim
     double precision, external :: ddot
 
     external approx, marginalxx, dpoisf, dnormf, dbinomf, dgammaf, dnbinomf, simgaussian
@@ -93,7 +93,8 @@ subroutine ngloglik(yt, ymiss, timevar, zt, tt, rtv, qt, a1, p1,p1inf, p,m, &
     end do
 
     if(sim .EQ. 1) then
-
+        allocate(w(nsim))
+        allocate(tsim(p, n, nsim2))
         w=1.0d0
         info2 = 0
 
@@ -151,7 +152,8 @@ subroutine ngloglik(yt, ymiss, timevar, zt, tt, rtv, qt, a1, p1,p1inf, p,m, &
             info = info2
             return
         end if
-
+        deallocate(tsim)
+        deallocate(w)
     end if
 
 end subroutine ngloglik

src/ngsmooth.f90

@@ -35,13 +35,15 @@ subroutine ngsmooth(yt, ymiss, timevar, zt, tt, rtv, qt, a1, p1,p1inf, u, theta,
     double precision, intent(inout), dimension((p-1)*smoothy+1,(p-1)*smoothy+1,(n-1)*smoothy+1) :: yvar
     double precision, intent(inout), dimension((p-1)*smooths+1,(n-1)*smooths+1) :: thetahat
     double precision, intent(inout), dimension((p-1)*smooths+1,(p-1)*smooths+1,(n-1)*smooths+1) :: thetavar
-    double precision, dimension(smootha*m+(1-smootha)*p,n,4*nsim) :: sim
-    double precision, dimension(smootha*p+(1-smootha),n,4*nsim*smootha+(1-smootha)) :: osim
-    double precision, dimension(4*nsim) :: w
     double precision, dimension(p,m) :: pm
+    double precision, dimension(:,:,:), allocatable :: sim
+    double precision, dimension(:,:,:), allocatable :: osim
+    double precision, dimension(:), allocatable :: w
     double precision, external :: ddot
-
     external isample, covmeanwprotect, dgemv, dsymm, dgemm, covmeanw
+    
+    allocate(sim(smootha*m+(1-smootha)*p,n,4*nsim))
+    allocate(w(4*nsim))
 
     if(smootha.EQ.1) then
 
@@ -66,6 +68,7 @@ subroutine ngsmooth(yt, ymiss, timevar, zt, tt, rtv, qt, a1, p1,p1inf, u, theta,
         end if
 
         if(smoothy.EQ.1) then
+            allocate(osim(smootha*p+(1-smootha),n,4*nsim*smootha+(1-smootha)))
             do t = 1, n
                 do j = 1,p
                     call dgemv('t',m,4*nsim,1.0d0,sim(:,t,:),m,zt(j,:,(t-1)*timevar(1)+1),1,0.0d0,osim(j,t,:),1)
@@ -86,6 +89,8 @@ subroutine ngsmooth(yt, ymiss, timevar, zt, tt, rtv, qt, a1, p1,p1inf, u, theta,
                 end select
             end do
             call covmeanw(osim,w,p,n,4*nsim,yhat,yvar)
+            
+           deallocate(osim)
         end if
     else
         call isample(yt, ymiss, timevar, zt, tt, rtv, qt, a1, p1,p1inf, u, dist, &
@@ -120,5 +125,6 @@ subroutine ngsmooth(yt, ymiss, timevar, zt, tt, rtv, qt, a1, p1,p1inf, u, theta,
         end if
 
     end if
-
+    deallocate(sim)
+    deallocate(w)
 end subroutine ngsmooth