Line data    Source code

       1             : /**
       2             :    Bojan Nikolic <bojan@bnikolic.co.uk> 
       3             :    Initial version 2003
       4             : 
       5             :    This file is part of BNMin1 and is licensed under GNU General
       6             :    Public License version 2.
       7             : 
       8             :    \file lmmin.cxx
       9             :    Renamed to lmmin.cc 2023.
      10             : 
      11             : */
      12             : 
      13             : #include "lmmin.h"
      14             : 
      15             : #include "pda.h"
      16             : 
      17             : 
      18             : 
      19             : namespace Minim {
      20             : 
      21             :   static Minim::LMMin* LMMin_this;
      22             :   static int LMMin_lock (0);
      23             : 
      24           0 :   extern "C" void LMMin_helper (int * /*m*/,
      25             :                                 int * /*n*/,
      26             :                                 double *x,
      27             :                                 double *fvec,
      28             :                                 int * /*iflag*/)
      29             :   {
      30           0 :     LMMin_this->copytopars(x);
      31           0 :     LMMin_this->ResEval();
      32           0 :     LMMin_this->copyres(fvec);
      33           0 :   }
      34             :   
      35             :   
      36           0 :   LMMin::LMMin( Minimisable &pm ): 
      37             :     Minimiser(pm),
      38           0 :     ftol(1e-3), 
      39           0 :     xtol(1e-3),
      40           0 :     gtol(1e-3),
      41           0 :     maxfev(1000),
      42           0 :     epsfcn(0.0),
      43           0 :     fjac(NParam() * NRes()),
      44           0 :     ldfjac(NRes()),
      45           0 :     ipvt(NParam())
      46             :   {   
      47           0 :     MonitorChi_stride=NParam();
      48           0 :   }
      49             :   
      50             :   
      51           0 :   void LMMin::solve(void)
      52             :   {
      53             :     
      54           0 :     if ( LMMin_lock ) 
      55             :       {
      56           0 :         return;
      57             :       }
      58             :     else 
      59             :       {
      60           0 :         LMMin_lock  =1;
      61             :       }
      62             :     
      63             :     
      64             :     
      65             :     // this is needed to communicate with the 
      66           0 :     LMMin_this = this; 
      67             :     
      68           0 :     InitRes();
      69             :     
      70           0 :     int n= NParam(); 
      71           0 :     int m= NRes();
      72             :     
      73           0 :     std::vector< double >   x0(  n ); 
      74           0 :     copyfrompars(&x0[0]);
      75             :     
      76           0 :     std::vector< double>    fvec( m );
      77             :     
      78           0 :     fjac.resize( n * m);
      79           0 :     ldfjac = m;
      80           0 :     ipvt.resize(n);
      81             :     
      82             :     
      83             : #ifndef BNMIN1_NO_PDA
      84             : 
      85             :     std::vector< double>    diag( n);
      86             :     int mode =1;
      87             :     double factor = 100.0;
      88             :     int nprint = 0;
      89             :     int info=0;
      90             :     int nfev=0;
      91             :     
      92             :     std::vector<double> qtf ( n);
      93             :     
      94             :     std::vector<double> wa1( n) , wa2( n), wa3( n) , wa4( m);
      95             :     
      96             :     pda_lmdif__( LMMin_helper,
      97             :                  &m,
      98             :                  &n,
      99             :                  &x0[0],
     100             :                  &fvec[0],
     101             :                  &ftol, &xtol , &gtol , &maxfev, &epsfcn,
     102             :                  &diag[0],
     103             :                  &mode ,
     104             :                  &factor,
     105             :                  &nprint, 
     106             :                  &info,
     107             :                  &nfev,
     108             :                  &fjac[0],
     109             :                  &ldfjac,
     110             :                  &ipvt[0],
     111             :                  &qtf[0],
     112             :                  &wa1[0],&wa2[0],&wa3[0],&wa4[0]
     113             :                  );
     114             : #else
     115           0 :     throw "No PDA algorithms so don't know how to do Lavenberg-Marquarndt";
     116             : #endif
     117             :     
     118             :     LMMin_lock  =0;
     119           0 :   }
     120             :   
     121             :   
     122             : }