ADMB Documentation  11.1.2490
 All Classes Files Functions Variables Typedefs Friends Defines
sgradclc.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: sgradclc.cpp 2402 2014-09-22 22:28:18Z johnoel $
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   #define lseek _lseek
00039   #define  read _read
00040   #define write _write
00041   #define ssize_t int
00042 #else
00043   #include <iostream>
00044   using namespace std;
00045   #include <sys/stat.h>
00046   #include <sys/types.h>
00047   #include <unistd.h>
00048 #endif
00049 
00050 #ifdef __SUN__
00051   #include <iostream.h>
00052   #include <sys/stat.h>
00053   #include <sys/types.h>
00054   #include <unistd.h>
00055 #endif
00056 
00057 #if defined(__NDPX__ )
00058   extern "C" {
00059     int lseek(int, int, int);
00060     int read(int, char*, int);
00061   };
00062 #endif
00063 
00064 #include <math.h>
00065 
00066 #if defined(__ZTC__)
00067   void _far * _cdecl _farptr_norm(void _far *);
00068   void _far * _cdecl _farptr_fromlong(unsigned long);
00069   long _cdecl _farptr_tolong(void _far *);
00070 #endif
00071 
00072 #ifndef OPT_LIB
00073   #include <cassert>
00074 #endif
00075 
00082 void gradcalc(int nvar, const dvector& _g)
00083 {
00084   if (nvar!=0)
00085   {
00086     if (nvar != gradient_structure::NVAR)
00087     {
00088       cerr << "nvar != gradient_structure::NVAR in gradcalc" << endl;
00089       cerr << "  nvar = " << nvar << endl;
00090       cerr << "  gradient_structure::NVAR = " << gradient_structure::NVAR
00091            << endl;
00092       cerr << "  in " __FILE__ << endl;
00093       ad_exit(1);
00094     }
00095   }
00096   dvector& g= (dvector&) _g;
00097   gradient_structure::TOTAL_BYTES = 0;
00098   gradient_structure::PREVIOUS_TOTAL_BYTES=0;
00099   unsigned int i;
00100   if(!gradient_structure::instances)
00101   {
00102     g.initialize();
00103     return;
00104   }
00105 
00106   if (g.size() < nvar)
00107   {
00108     cerr  << "gradient vector size is less than the number of variables.\n";
00109     ad_exit(1);
00110   }
00111 
00112    gradient_structure::GRAD_STACK1->_GRADFILE_PTR =
00113      gradient_structure::GRAD_STACK1->gradfile_handle();
00114 
00115   int& _GRADFILE_PTR=gradient_structure::GRAD_STACK1->_GRADFILE_PTR;
00116 
00117   lseek(_GRADFILE_PTR,0L,SEEK_CUR);
00118 
00119   if (gradient_structure::GRAD_STACK1->ptr <=
00120         gradient_structure::GRAD_STACK1->ptr_first)
00121   {
00122 #ifdef SAFE_ALL
00123     cerr << "warning -- calling gradcalc when no calculations generating"
00124          << endl << "derivative information have occurred" << endl;
00125 #endif
00126     g.initialize();
00127     return;
00128   }    // current is one past the end so -- it
00129 
00130   if (gradient_structure::save_var_flag)
00131   {
00132     gradient_structure::save_arrays();
00133     gradient_structure::save_variables();
00134   }
00135 
00136   gradient_structure::GRAD_STACK1->ptr--;
00137 
00138   for (i = 0; i < gradient_structure::GRAD_LIST->nlinks; i++)
00139   {
00140     *(double*)(gradient_structure::GRAD_LIST->dlink_addresses[i]) = 0;
00141   }
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   unsigned int size = sizeof(double_and_int);
00149   for (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 lpos = lseek(gradient_structure::GRAD_STACK1->_GRADFILE_PTR,
00187       -((long int)(sizeof(grad_stack_entry)*gradient_structure::
00188         GRAD_STACK1->length)),SEEK_CUR);
00189 
00190     break_flag=gradient_structure::GRAD_STACK1->read_grad_stack_buffer(lpos);
00191   } while (break_flag);
00192 
00193   {
00194 #ifdef GRAD_DIAG
00195     long int ttmp =
00196 #endif
00197     lseek(gradient_structure::GRAD_STACK1->_GRADFILE_PTR, 0,SEEK_CUR);
00198 #ifdef GRAD_DIAG
00199     cout << "Offset in file at end of gradcalc is " << ttmp
00200          << " bytes from the beginning\n";
00201 #endif
00202   }
00203 
00204   int mindx = g.indexmin();
00205   for (i=0; i < (unsigned int)nvar; i++)
00206   {
00207     g[i + mindx] = *gradient_structure::INDVAR_LIST->get_address(i);
00208   }
00209 
00210   gradient_structure::GRAD_STACK1->ptr =
00211     gradient_structure::GRAD_STACK1->ptr_first;
00212 
00213   if (gradient_structure::save_var_flag)
00214   {
00215     gradient_structure::restore_arrays();
00216     gradient_structure::restore_variables();
00217   }
00218 }
00227 double gradcalc(int nvar, const dvector& _g, dvariable& f)
00228 {
00229   double v = value(f);
00230   gradcalc(nvar, _g);
00231   return v;
00232 }
00235 void gradient_structure::save_arrays()
00236 {
00237   void * temp_ptr;
00238   long bytes_needed=min(gradient_structure::ARR_LIST1->get_last_offset()+1,
00239     ARRAY_MEMBLOCK_SIZE);
00240   gradient_structure::save_var_file_flag=0;
00241 #ifdef __ZTC__
00242    if ( (temp_ptr = farmalloc(bytes_needed) ) == 0)
00243 #else
00244    //if ( (temp_ptr = malloc(bytes_needed) ) == 0)
00245    if ((temp_ptr = (void *)malloc(bytes_needed)) == 0)
00246   #define __USE_IOSTREAM__
00247 #endif
00248    {
00249      gradient_structure::save_var_file_flag=1;
00250      cerr << "insufficient memory to allocate space for ARRAY_MEMBLOCK"
00251           << " save buffer " << endl;
00252    }
00253    if (gradient_structure::save_var_file_flag==0)
00254    {
00255      ARRAY_MEMBLOCK_SAVE = temp_ptr;
00256 #if defined(DOS386)
00257   #ifndef USE_ASSEMBLER
00258          memcpy((char*)ARRAY_MEMBLOCK_SAVE,(char*)ARRAY_MEMBLOCK_BASE,
00259            bytes_needed);
00260   #else
00261          dw_block_move((double*)ARRAY_MEMBLOCK_SAVE,
00262            (double*)ARRAY_MEMBLOCK_BASE,bytes_needed/8);
00263   #endif
00264 #else
00265      unsigned long int max_move=50000;
00266      unsigned long int left_to_move=bytes_needed;
00267      humungous_pointer dest = ARRAY_MEMBLOCK_SAVE;
00268      humungous_pointer src = ARRAY_MEMBLOCK_BASE;
00269      while(left_to_move > max_move)
00270      {
00271        memcpy((char*)dest,(char*)src,max_move);
00272        left_to_move-=max_move;
00273        src+=max_move;
00274        dest+=max_move;
00275      }
00276      memcpy((char*)dest,(char*)src,left_to_move);
00277 #endif
00278   }
00279   else
00280   {
00281      humungous_pointer src = ARRAY_MEMBLOCK_BASE;
00282      lseek(gradient_structure::GRAD_STACK1->_VARSSAV_PTR,0L,SEEK_SET);
00283 #if defined(DOS386)
00284   #ifdef OPT_LIB
00285      write(gradient_structure::GRAD_STACK1->_VARSSAV_PTR,
00286        (char*)src, bytes_needed);
00287   #else
00288      ssize_t ret = write(gradient_structure::GRAD_STACK1->_VARSSAV_PTR,
00289        (char*)src, bytes_needed);
00290      assert(ret != -1);
00291   #endif
00292 #else
00293      unsigned long int max_move=500;
00294      unsigned long int left_to_move=bytes_needed;
00295      while(left_to_move > max_move)
00296      {
00297        write(_VARSSAV_PTR,(char*)src,max_move);
00298        left_to_move-=max_move;
00299        src+=max_move;
00300      }
00301      write(gradient_structure::GRAD_STACK1->_VARSSAV_PTR,(char*)src,
00302        left_to_move);
00303 #endif
00304   }
00305 }
00306 
00309 void gradient_structure::restore_arrays()
00310 {
00311   long bytes_needed=min(gradient_structure::ARR_LIST1->get_last_offset()+1,
00312     ARRAY_MEMBLOCK_SIZE);
00313   if (gradient_structure::save_var_file_flag==0)
00314   {
00315 #if defined(DOS386)
00316   #ifndef USE_ASSEMBLER
00317         memcpy((char*)ARRAY_MEMBLOCK_BASE,(char*)ARRAY_MEMBLOCK_SAVE,
00318           bytes_needed);
00319   #else
00320          dw_block_move((double*)ARRAY_MEMBLOCK_BASE,
00321            (double*)ARRAY_MEMBLOCK_SAVE,bytes_needed/8);
00322   #endif
00323 #else
00324      unsigned long max_move=50000;
00325 
00326      long int left_to_move=bytes_needed;
00327      humungous_pointer src = ARRAY_MEMBLOCK_SAVE;
00328      humungous_pointer dest = ARRAY_MEMBLOCK_BASE;
00329      while(left_to_move > max_move)
00330      {
00331        memcpy((char*)dest,(char*)src,max_move);
00332        left_to_move-=max_move;
00333        src+=max_move;
00334        dest+=max_move;
00335      }
00336      memcpy((char*)dest,(char*)src,left_to_move);
00337 #endif
00338     ARRAY_MEMBLOCK_SAVE.free();
00339   }
00340   else
00341   {
00342     humungous_pointer dest = ARRAY_MEMBLOCK_BASE;
00343     lseek(gradient_structure::GRAD_STACK1->_VARSSAV_PTR,0L,SEEK_SET);
00344 #if defined(DOS386)
00345   #ifdef OPT_LIB
00346     read(gradient_structure::GRAD_STACK1->_VARSSAV_PTR,
00347       (char*)dest,bytes_needed);
00348   #else
00349     ssize_t ret = read(gradient_structure::GRAD_STACK1->_VARSSAV_PTR,
00350       (char*)dest,bytes_needed);
00351     assert(ret != -1);
00352   #endif
00353 #else
00354      unsigned long int max_move=50000;
00355 
00356      long int left_to_move=bytes_needed;
00357      while(left_to_move > max_move)
00358      {
00359        read(gradient_structure::GRAD_STACK1->_VARSSAV_PTR,
00360          (char*)dest,max_move);
00361        left_to_move-=max_move;
00362        dest+=max_move;
00363      }
00364      read(gradient_structure::GRAD_STACK1->_VARSSAV_PTR,
00365        (char*)dest,left_to_move);
00366 #endif
00367   }
00368 }
00372 void gradient_structure::save_variables()
00373 {
00374   if ((variables_save = new double[gradient_structure::MAX_DLINKS])==NULL)
00375   {
00376     //_VARSSAV_PTR=open(var_store_file_name, O_RDWR | O_CREAT | O_TRUNC ,
00377     //               S_IREAD | S_IWRITE);
00378     cerr << "insufficient memory to allocate space for dvariables"
00379          << " save buffer " << endl;
00380     ad_exit(1);
00381   }
00382   for (unsigned int i=0; i<gradient_structure::GRAD_LIST->nlinks; i++)
00383   {
00384     //variables_save[i]= * (double*)
00385     //    (gradient_structure::GRAD_LIST->dlink_addresses[i]);
00386     memcpy(&(variables_save[i]),
00387        gradient_structure::GRAD_LIST->dlink_addresses[i],sizeof(double));
00388   }
00389 }
00393 void gradient_structure::restore_variables()
00394 {
00395   for (unsigned int i=0; i<gradient_structure::GRAD_LIST->nlinks; i++)
00396   {
00397     //* (double*)(gradient_structure::GRAD_LIST->dlink_addresses[i])
00398     //  = variables_save[i];
00399     memcpy(gradient_structure::GRAD_LIST->dlink_addresses[i],
00400        &(variables_save[i]),sizeof(double));
00401   }
00402   delete [] variables_save;
00403   variables_save = 0;
00404 }
00408 void reset_gradient_stack(void)
00409 {
00410   gradient_structure::GRAD_STACK1->ptr =
00411     gradient_structure::GRAD_STACK1->ptr_first;
00412 
00413   int& _GRADFILE_PTR=gradient_structure::GRAD_STACK1->_GRADFILE_PTR;
00414 
00415   lseek(_GRADFILE_PTR,0L,SEEK_SET);
00416 }
00426 void grad_stack::set_gradient_stack1(void (* func)(void),
00427   double* dep_addr, double* ind_addr1)
00428 {
00429 #ifdef NO_DERIVS
00430   if (!gradient_structure::no_derivatives)
00431   {
00432 #endif
00433     if (ptr > ptr_last)
00434     {
00435        // current buffer is full -- write it to disk and reset pointer
00436        // and counter
00437        this->write_grad_stack_buffer();
00438     }
00439     //test_the_pointer();
00440     ptr->func = func;
00441     ptr->dep_addr = dep_addr;
00442     ptr->ind_addr1 = ind_addr1;
00443     ptr++;
00444 #ifdef NO_DERIVS
00445   }
00446 #endif
00447 }
00448 void test_the_pointer(void)
00449 {
00450 /*
00451   static int inner_count = 0;
00452   static grad_stack_entry * pgse = (grad_stack_entry*) (0x1498fffc);
00453   inner_count++;
00454   if (inner_count == 404849)
00455   {
00456     cout << ptr << endl;
00457     cout << ptr->func << endl;
00458     cout << ptr->dep_addr << endl;
00459     cout << (int)(ptr->dep_addr)%8 << endl;
00460   }
00461   pgse->func = (void (*)())(100);
00462   pgse->dep_addr = (double*) 100;
00463   pgse->ind_addr1 = (double*) 100;
00464 */
00465 }