ADMB Documentation  11.4.2891
 All Classes Files Functions Variables Typedefs Friends Defines
sgradclc.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id$
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #if !defined(DOS386)
00012   #define DOS386
00013 #endif
00014 
00015 #include "fvar.hpp"
00016 
00017 #include <sys/stat.h>
00018 #include <fcntl.h>
00019 #include <string.h>
00020 
00021 #ifdef __TURBOC__
00022   #pragma hdrstop
00023   #include <iostream.h>
00024 #endif
00025 
00026 #ifdef __ZTC__
00027   #include <iostream.hpp>
00028 #endif
00029 
00030 #if defined (__WAT32__)
00031 #  include <io.h>
00032 #endif
00033 
00034 #include <stdio.h>
00035 #include <stdlib.h>
00036 
00037 #ifdef _MSC_VER
00038   #ifdef _M_X64
00039   typedef __int64 ssize_t;
00040   #else
00041   typedef int ssize_t;
00042   #endif
00043   #define lseek _lseek
00044   #define  read _read
00045   #define write _write
00046 #else
00047   #include <iostream>
00048   using namespace std;
00049   #include <sys/stat.h>
00050   #include <sys/types.h>
00051   #include <unistd.h>
00052 #endif
00053 
00054 #ifdef __SUN__
00055   #include <iostream.h>
00056   #include <sys/stat.h>
00057   #include <sys/types.h>
00058   #include <unistd.h>
00059 #endif
00060 
00061 #if defined(__NDPX__ )
00062   extern "C" {
00063     int lseek(int, int, int);
00064     int read(int, char*, int);
00065   };
00066 #endif
00067 
00068 #include <math.h>
00069 
00070 #if defined(__ZTC__)
00071   void _far * _cdecl _farptr_norm(void _far *);
00072   void _far * _cdecl _farptr_fromlong(unsigned long);
00073   long _cdecl _farptr_tolong(void _far *);
00074 #endif
00075 
00076 #if !defined(OPT_LIB) || defined(_MSC_VER)
00077   #include <cassert>
00078 #endif
00079 
00086 void gradcalc(int nvar, const dvector& _g)
00087 {
00088   if (nvar!=0)
00089   {
00090     if (nvar != gradient_structure::NVAR)
00091     {
00092       cerr << "nvar != gradient_structure::NVAR in gradcalc" << endl;
00093       cerr << "  nvar = " << nvar << endl;
00094       cerr << "  gradient_structure::NVAR = " << gradient_structure::NVAR
00095            << endl;
00096       cerr << "  in " __FILE__ << endl;
00097       ad_exit(1);
00098     }
00099   }
00100   dvector& g= (dvector&) _g;
00101   gradient_structure::TOTAL_BYTES = 0;
00102   gradient_structure::PREVIOUS_TOTAL_BYTES=0;
00103   if(!gradient_structure::instances)
00104   {
00105     g.initialize();
00106     return;
00107   }
00108 
00109   if (g.size() < nvar)
00110   {
00111     cerr  << "gradient vector size is less than the number of variables.\n";
00112     ad_exit(1);
00113   }
00114 
00115    gradient_structure::GRAD_STACK1->_GRADFILE_PTR =
00116      gradient_structure::GRAD_STACK1->gradfile_handle();
00117 
00118   int& _GRADFILE_PTR=gradient_structure::GRAD_STACK1->_GRADFILE_PTR;
00119 
00120   lseek(_GRADFILE_PTR,0L,SEEK_CUR);
00121 
00122   if (gradient_structure::GRAD_STACK1->ptr <=
00123         gradient_structure::GRAD_STACK1->ptr_first)
00124   {
00125 #ifdef SAFE_ALL
00126     cerr << "warning -- calling gradcalc when no calculations generating"
00127          << endl << "derivative information have occurred" << endl;
00128 #endif
00129     g.initialize();
00130     return;
00131   }    // current is one past the end so -- it
00132 
00133   if (gradient_structure::save_var_flag)
00134   {
00135     gradient_structure::save_arrays();
00136     gradient_structure::save_variables();
00137   }
00138 
00139   gradient_structure::GRAD_STACK1->ptr--;
00140 
00141   gradient_structure::GRAD_LIST->initialize();
00142 
00143   double_and_int* tmp =
00144     (double_and_int*)gradient_structure::ARRAY_MEMBLOCK_BASE;
00145 
00146   unsigned long int max_last_offset =
00147     gradient_structure::ARR_LIST1->get_max_last_offset();
00148   size_t size = sizeof(double_and_int);
00149   for (unsigned int i = 0; i < (max_last_offset/size); i++)
00150   {
00151      tmp->x = 0;
00152 #if defined (__ZTC__)
00153      tmp = (double_and_int*)_farptr_norm((void*)(++tmp));
00154 #else
00155      tmp++;
00156 #endif
00157   }
00158 
00159   *gradient_structure::GRAD_STACK1->ptr->dep_addr = 1;
00160 
00161   //int icount=0;
00162   int break_flag=1;
00163   do
00164   {
00165     gradient_structure::GRAD_STACK1->ptr++;
00166 #ifdef FAST_ASSEMBLER
00167     gradloop();
00168 #else
00169     grad_stack_entry* grad_ptr_first =
00170       gradient_structure::GRAD_STACK1->ptr_first;
00171     while (gradient_structure::GRAD_STACK1->ptr-- >
00172              grad_ptr_first)
00173     {
00174       (*(gradient_structure::GRAD_STACK1->ptr->func))();
00175 /*
00176       icount++;
00177       if (icount%1000==0)
00178       {
00179         //cout << "icount = " << icount << endl;
00180       }
00181 */
00182     }
00183 #endif
00184 
00185    // back up the file one buffer size and read forward
00186    off_t offset = (off_t)(sizeof(grad_stack_entry)
00187                   * gradient_structure::GRAD_STACK1->length);
00188    off_t lpos = lseek(gradient_structure::GRAD_STACK1->_GRADFILE_PTR,
00189       -offset,SEEK_CUR);
00190 
00191     break_flag=gradient_structure::GRAD_STACK1->read_grad_stack_buffer(lpos);
00192   } while (break_flag);
00193 
00194   {
00195 #ifdef GRAD_DIAG
00196     long int ttmp =
00197 #endif
00198     lseek(gradient_structure::GRAD_STACK1->_GRADFILE_PTR, 0,SEEK_CUR);
00199 #ifdef GRAD_DIAG
00200     cout << "Offset in file at end of gradcalc is " << ttmp
00201          << " bytes from the beginning\n";
00202 #endif
00203   }
00204 
00205 #ifndef OPT_LIB
00206   assert(g.indexmin() > 0);
00207 #endif
00208   int mindx = g.indexmin();
00209   for (int i=0; i < nvar; i++)
00210   {
00211     g[i + mindx] = *gradient_structure::INDVAR_LIST->get_address(i);
00212   }
00213 
00214   gradient_structure::GRAD_STACK1->ptr =
00215     gradient_structure::GRAD_STACK1->ptr_first;
00216 
00217   if (gradient_structure::save_var_flag)
00218   {
00219     gradient_structure::restore_arrays();
00220     gradient_structure::restore_variables();
00221   }
00222 }
00231 double gradcalc(int nvar, const dvector& _g, dvariable& f)
00232 {
00233   double v = value(f);
00234   gradcalc(nvar, _g);
00235   return v;
00236 }
00239 void gradient_structure::save_arrays()
00240 {
00241   void * temp_ptr;
00242   unsigned long bytes_needed =
00243     min(gradient_structure::ARR_LIST1->get_last_offset() + 1,
00244         ARRAY_MEMBLOCK_SIZE);
00245   gradient_structure::save_var_file_flag=0;
00246 #ifdef __ZTC__
00247    if ( (temp_ptr = farmalloc(bytes_needed) ) == 0)
00248 #else
00249    //if ( (temp_ptr = malloc(bytes_needed) ) == 0)
00250    if ((temp_ptr = (void *)malloc(bytes_needed)) == 0)
00251   #define __USE_IOSTREAM__
00252 #endif
00253    {
00254      gradient_structure::save_var_file_flag=1;
00255      cerr << "insufficient memory to allocate space for ARRAY_MEMBLOCK"
00256           << " save buffer " << endl;
00257    }
00258    if (gradient_structure::save_var_file_flag==0)
00259    {
00260      ARRAY_MEMBLOCK_SAVE = temp_ptr;
00261 #if defined(DOS386)
00262   #ifndef USE_ASSEMBLER
00263          memcpy((char*)ARRAY_MEMBLOCK_SAVE,(char*)ARRAY_MEMBLOCK_BASE,
00264            bytes_needed);
00265   #else
00266          dw_block_move((double*)ARRAY_MEMBLOCK_SAVE,
00267            (double*)ARRAY_MEMBLOCK_BASE,bytes_needed/8);
00268   #endif
00269 #else
00270      unsigned long int max_move=50000;
00271      unsigned long int left_to_move=bytes_needed;
00272      humungous_pointer dest = ARRAY_MEMBLOCK_SAVE;
00273      humungous_pointer src = ARRAY_MEMBLOCK_BASE;
00274      while(left_to_move > max_move)
00275      {
00276        memcpy((char*)dest,(char*)src,max_move);
00277        left_to_move-=max_move;
00278        src+=max_move;
00279        dest+=max_move;
00280      }
00281      memcpy((char*)dest,(char*)src,left_to_move);
00282 #endif
00283   }
00284   else
00285   {
00286      humungous_pointer src = ARRAY_MEMBLOCK_BASE;
00287      lseek(gradient_structure::GRAD_STACK1->_VARSSAV_PTR,0L,SEEK_SET);
00288 #if defined(DOS386)
00289   #ifdef OPT_LIB
00290      write(gradient_structure::GRAD_STACK1->_VARSSAV_PTR,
00291        (char*)src, bytes_needed);
00292   #else
00293      ssize_t ret = write(gradient_structure::GRAD_STACK1->_VARSSAV_PTR,
00294        (char*)src, bytes_needed);
00295      assert(ret != -1);
00296   #endif
00297 #else
00298      unsigned long int max_move=500;
00299      unsigned long int left_to_move=bytes_needed;
00300      while(left_to_move > max_move)
00301      {
00302        write(_VARSSAV_PTR,(char*)src,max_move);
00303        left_to_move-=max_move;
00304        src+=max_move;
00305      }
00306      write(gradient_structure::GRAD_STACK1->_VARSSAV_PTR,(char*)src,
00307        left_to_move);
00308 #endif
00309   }
00310 }
00311 
00314 void gradient_structure::restore_arrays()
00315 {
00316   unsigned long bytes_needed =
00317     min(gradient_structure::ARR_LIST1->get_last_offset() + 1,
00318         ARRAY_MEMBLOCK_SIZE);
00319   if (gradient_structure::save_var_file_flag==0)
00320   {
00321 #if defined(DOS386)
00322   #ifndef USE_ASSEMBLER
00323         memcpy((char*)ARRAY_MEMBLOCK_BASE,(char*)ARRAY_MEMBLOCK_SAVE,
00324           bytes_needed);
00325   #else
00326          dw_block_move((double*)ARRAY_MEMBLOCK_BASE,
00327            (double*)ARRAY_MEMBLOCK_SAVE,bytes_needed/8);
00328   #endif
00329 #else
00330      unsigned long max_move=50000;
00331 
00332      long int left_to_move=bytes_needed;
00333      humungous_pointer src = ARRAY_MEMBLOCK_SAVE;
00334      humungous_pointer dest = ARRAY_MEMBLOCK_BASE;
00335      while(left_to_move > max_move)
00336      {
00337        memcpy((char*)dest,(char*)src,max_move);
00338        left_to_move-=max_move;
00339        src+=max_move;
00340        dest+=max_move;
00341      }
00342      memcpy((char*)dest,(char*)src,left_to_move);
00343 #endif
00344     ARRAY_MEMBLOCK_SAVE.free();
00345   }
00346   else
00347   {
00348     humungous_pointer dest = ARRAY_MEMBLOCK_BASE;
00349     lseek(gradient_structure::GRAD_STACK1->_VARSSAV_PTR,0L,SEEK_SET);
00350 #if defined(DOS386)
00351   #if defined(OPT_LIB) && !defined(_MSC_VER)
00352     read(gradient_structure::GRAD_STACK1->_VARSSAV_PTR,
00353       (char*)dest,bytes_needed);
00354   #else
00355     ssize_t ret = read(gradient_structure::GRAD_STACK1->_VARSSAV_PTR,
00356       (char*)dest,bytes_needed);
00357     assert(ret != -1);
00358   #endif
00359 #else
00360      unsigned long int max_move=50000;
00361 
00362      long int left_to_move=bytes_needed;
00363      while(left_to_move > max_move)
00364      {
00365        read(gradient_structure::GRAD_STACK1->_VARSSAV_PTR,
00366          (char*)dest,max_move);
00367        left_to_move-=max_move;
00368        dest+=max_move;
00369      }
00370      read(gradient_structure::GRAD_STACK1->_VARSSAV_PTR,
00371        (char*)dest,left_to_move);
00372 #endif
00373   }
00374 }
00378 void gradient_structure::save_variables()
00379 {
00380   GRAD_LIST->save_variables();
00381 }
00385 void gradient_structure::restore_variables()
00386 {
00387   GRAD_LIST->restore_variables();
00388 }
00392 void reset_gradient_stack(void)
00393 {
00394   gradient_structure::GRAD_STACK1->ptr =
00395     gradient_structure::GRAD_STACK1->ptr_first;
00396 
00397   int& _GRADFILE_PTR=gradient_structure::GRAD_STACK1->_GRADFILE_PTR;
00398 
00399   lseek(_GRADFILE_PTR,0L,SEEK_SET);
00400 }
00410 void grad_stack::set_gradient_stack1(void (* func)(void),
00411   double* dep_addr, double* ind_addr1)
00412 {
00413 #ifdef NO_DERIVS
00414   if (!gradient_structure::no_derivatives)
00415   {
00416 #endif
00417     if (ptr > ptr_last)
00418     {
00419        // current buffer is full -- write it to disk and reset pointer
00420        // and counter
00421        this->write_grad_stack_buffer();
00422     }
00423     //test_the_pointer();
00424     ptr->func = func;
00425     ptr->dep_addr = dep_addr;
00426     ptr->ind_addr1 = ind_addr1;
00427     ptr++;
00428 #ifdef NO_DERIVS
00429   }
00430 #endif
00431 }
00432 void test_the_pointer(void)
00433 {
00434 /*
00435   static int inner_count = 0;
00436   static grad_stack_entry * pgse = (grad_stack_entry*) (0x1498fffc);
00437   inner_count++;
00438   if (inner_count == 404849)
00439   {
00440     cout << ptr << endl;
00441     cout << ptr->func << endl;
00442     cout << ptr->dep_addr << endl;
00443     cout << (int)(ptr->dep_addr)%8 << endl;
00444   }
00445   pgse->func = (void (*)())(100);
00446   pgse->dep_addr = (double*) 100;
00447   pgse->ind_addr1 = (double*) 100;
00448 */
00449 }