Maximum likelihood parameters optimization program for CSA added without examples.
[unres.git] / source / maxlik / src_CSA / minsumsl.f
diff --git a/source/maxlik/src_CSA/minsumsl.f b/source/maxlik/src_CSA/minsumsl.f
new file mode 100644 (file)
index 0000000..53105e5
--- /dev/null
@@ -0,0 +1,86 @@
+      subroutine minsumsl(nvar,x,minval)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      parameter (maxvar=maxene+3*nnbase)
+      parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2)) 
+*********************************************************************
+* OPTIMIZE sets up SUMSL or DFP and provides a simple interface for *
+* the calling subprogram.                                           *     
+* when d(i)=1.0, then v(35) is the length of the initial step,      *     
+* calculated in the usual pythagorean way.                          *     
+* absolute convergence occurs when the function is within v(31) of  *     
+* zero. unless you know the minimum value in advance, abs convg     *     
+* is probably not useful.                                           *     
+* relative convergence is when the model predicts that the function *   
+* will decrease by less than v(32)*abs(fun).                        *   
+*********************************************************************
+      dimension iv(liv)                                               
+      real*8 minval,x(nvar),d(maxvar),v(1:lv)                     
+      external funclik,grad,fdum
+      integer idum(1)
+      double precision rdum(1)
+      double precision urparm(maxT)
+      double precision g(maxvar)
+      call deflt(2,iv,liv,lv,v)                                         
+* 12 means fresh start, dont call deflt                                 
+      iv(1)=12                                                          
+* max num of fun calls                                                  
+      maxfun=1000
+      iv(17)=maxfun
+* max num of iterations                                                 
+      maxit=50
+      iv(18)=maxit
+* controls output                                                       
+      iv(19)=1                                                          
+* selects output unit                                                   
+      iv(21)=2
+* 1 means to print out result                                           
+      iv(22)=1                                                          
+* 1 means to print out summary stats                                    
+      iv(23)=1                                                          
+* 1 means to print initial x and d                                      
+      iv(24)=1                                                          
+* min val for v(radfac) default is 0.1                                  
+      v(24)=0.01D0                                                       
+* max val for v(radfac) default is 4.0                                  
+      v(25)=2.0D0                                                       
+c     v(25)=4.0D0                                                       
+* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)    
+* the sumsl default is 0.1                                              
+      v(26)=0.001D0
+* false conv if (act fnctn decrease) .lt. v(34)                         
+* the sumsl default is 100*machep                                       
+      v(34)=v(34)/100.0D0                                               
+* absolute convergence                                                  
+      tolf=1.0D-4
+      v(31)=tolf
+* relative convergence                                                  
+      rtolf=1.0D-12
+      v(32)=rtolf
+* controls initial step size                                            
+       v(35)=1.0D-6                                                    
+* large vals of d correspond to small components of step                
+      do 20 i=1,nvar
+         d(i)=1.0D0                                                     
+20    continue
+
+      nf=0
+      call funclik(nvar,x,nf,f,idum,urparm,fdum)  
+      write (2,'(a,1pe17.10)') 'Initial function value:',f
+      call grad(nvar,x,nf,g,idum,urparm,fdum)
+      write (2,*) "Initial gradient"
+      do i=1,nvar
+        write (2,'(i5,e15.5)') i,g(i)
+      enddo
+c minimize the log-likelihood function
+      print *,"iv1",iv(1)
+      call sumsl(nvar,d,x,funclik,grad,iv,liv,lv,v,idum,urparm,fdum)
+      minval=v(10)                                                      
+      write (2,*)
+      write (2,'(a,i4)') 'SUMSL return code:',iv(1)
+      write (2,'(a,1pe17.10)') 'Final function value:',minval
+c      print *,"exiting minsumsl"
+      return  
+      end  
+c---------------------------------------------------------------------
+