ADMB Documentation  11.2.2828
 All Classes Files Functions Variables Typedefs Friends Defines
sgradclc.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: sgradclc.cpp 2735 2014-11-26 23:47:25Z 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   #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   for (unsigned int i = 0; i < gradient_structure::GRAD_LIST->nlinks; i++)
00142   {
00143     *(double*)(gradient_structure::GRAD_LIST->dlink_addresses[i]) = 0;
00144   }
00145 
00146   double_and_int* tmp =
00147     (double_and_int*)gradient_structure::ARRAY_MEMBLOCK_BASE;
00148 
00149   unsigned long int max_last_offset =
00150     gradient_structure::ARR_LIST1->get_max_last_offset();
00151   unsigned int size = sizeof(double_and_int);
00152   for (unsigned int i = 0; i < (max_last_offset/size); i++)
00153   {
00154      tmp->x = 0;
00155 #if defined (__ZTC__)
00156      tmp = (double_and_int*)_farptr_norm((void*)(++tmp));
00157 #else
00158      tmp++;
00159 #endif
00160   }
00161 
00162   *gradient_structure::GRAD_STACK1->ptr->dep_addr = 1;
00163 
00164   //int icount=0;
00165   int break_flag=1;
00166   do
00167   {
00168     gradient_structure::GRAD_STACK1->ptr++;
00169 #ifdef FAST_ASSEMBLER
00170     gradloop();
00171 #else
00172     grad_stack_entry* grad_ptr_first =
00173       gradient_structure::GRAD_STACK1->ptr_first;
00174     while (gradient_structure::GRAD_STACK1->ptr-- >
00175              grad_ptr_first)
00176     {
00177       (*(gradient_structure::GRAD_STACK1->ptr->func))();
00178 /*
00179       icount++;
00180       if (icount%1000==0)
00181       {
00182         //cout << "icount = " << icount << endl;
00183       }
00184 */
00185     }
00186 #endif
00187 
00188    // back up the file one buffer size and read forward
00189    off_t offset = (off_t)(sizeof(grad_stack_entry)
00190                   * gradient_structure::GRAD_STACK1->length);
00191    off_t lpos = lseek(gradient_structure::GRAD_STACK1->_GRADFILE_PTR,
00192       -offset,SEEK_CUR);
00193 
00194     break_flag=gradient_structure::GRAD_STACK1->read_grad_stack_buffer(lpos);
00195   } while (break_flag);
00196 
00197   {
00198 #ifdef GRAD_DIAG
00199     long int ttmp =
00200 #endif
00201     lseek(gradient_structure::GRAD_STACK1->_GRADFILE_PTR, 0,SEEK_CUR);
00202 #ifdef GRAD_DIAG
00203     cout << "Offset in file at end of gradcalc is " << ttmp
00204          << " bytes from the beginning\n";
00205 #endif
00206   }
00207 
00208 #ifndef OPT_LIB
00209   assert(g.indexmin() > 0);
00210 #endif
00211   int mindx = g.indexmin();
00212   for (int i=0; i < nvar; i++)
00213   {
00214     g[i + mindx] = *gradient_structure::INDVAR_LIST->get_address(i);
00215   }
00216 
00217   gradient_structure::GRAD_STACK1->ptr =
00218     gradient_structure::GRAD_STACK1->ptr_first;
00219 
00220   if (gradient_structure::save_var_flag)
00221   {
00222     gradient_structure::restore_arrays();
00223     gradient_structure::restore_variables();
00224   }
00225 }
00234 double gradcalc(int nvar, const dvector& _g, dvariable& f)
00235 {
00236   double v = value(f);
00237   gradcalc(nvar, _g);
00238   return v;
00239 }
00242 void gradient_structure::save_arrays()
00243 {
00244   void * temp_ptr;
00245   unsigned long bytes_needed =
00246     min(gradient_structure::ARR_LIST1->get_last_offset() + 1,
00247         ARRAY_MEMBLOCK_SIZE);
00248   gradient_structure::save_var_file_flag=0;
00249 #ifdef __ZTC__
00250    if ( (temp_ptr = farmalloc(bytes_needed) ) == 0)
00251 #else
00252    //if ( (temp_ptr = malloc(bytes_needed) ) == 0)
00253    if ((temp_ptr = (void *)malloc(bytes_needed)) == 0)
00254   #define __USE_IOSTREAM__
00255 #endif
00256    {
00257      gradient_structure::save_var_file_flag=1;
00258      cerr << "insufficient memory to allocate space for ARRAY_MEMBLOCK"
00259           << " save buffer " << endl;
00260    }
00261    if (gradient_structure::save_var_file_flag==0)
00262    {
00263      ARRAY_MEMBLOCK_SAVE = temp_ptr;
00264 #if defined(DOS386)
00265   #ifndef USE_ASSEMBLER
00266          memcpy((char*)ARRAY_MEMBLOCK_SAVE,(char*)ARRAY_MEMBLOCK_BASE,
00267            bytes_needed);
00268   #else
00269          dw_block_move((double*)ARRAY_MEMBLOCK_SAVE,
00270            (double*)ARRAY_MEMBLOCK_BASE,bytes_needed/8);
00271   #endif
00272 #else
00273      unsigned long int max_move=50000;
00274      unsigned long int left_to_move=bytes_needed;
00275      humungous_pointer dest = ARRAY_MEMBLOCK_SAVE;
00276      humungous_pointer src = ARRAY_MEMBLOCK_BASE;
00277      while(left_to_move > max_move)
00278      {
00279        memcpy((char*)dest,(char*)src,max_move);
00280        left_to_move-=max_move;
00281        src+=max_move;
00282        dest+=max_move;
00283      }
00284      memcpy((char*)dest,(char*)src,left_to_move);
00285 #endif
00286   }
00287   else
00288   {
00289      humungous_pointer src = ARRAY_MEMBLOCK_BASE;
00290      lseek(gradient_structure::GRAD_STACK1->_VARSSAV_PTR,0L,SEEK_SET);
00291 #if defined(DOS386)
00292   #ifdef OPT_LIB
00293      write(gradient_structure::GRAD_STACK1->_VARSSAV_PTR,
00294        (char*)src, bytes_needed);
00295   #else
00296      ssize_t ret = write(gradient_structure::GRAD_STACK1->_VARSSAV_PTR,
00297        (char*)src, bytes_needed);
00298      assert(ret != -1);
00299   #endif
00300 #else
00301      unsigned long int max_move=500;
00302      unsigned long int left_to_move=bytes_needed;
00303      while(left_to_move > max_move)
00304      {
00305        write(_VARSSAV_PTR,(char*)src,max_move);
00306        left_to_move-=max_move;
00307        src+=max_move;
00308      }
00309      write(gradient_structure::GRAD_STACK1->_VARSSAV_PTR,(char*)src,
00310        left_to_move);
00311 #endif
00312   }
00313 }
00314 
00317 void gradient_structure::restore_arrays()
00318 {
00319   unsigned long bytes_needed =
00320     min(gradient_structure::ARR_LIST1->get_last_offset() + 1,
00321         ARRAY_MEMBLOCK_SIZE);
00322   if (gradient_structure::save_var_file_flag==0)
00323   {
00324 #if defined(DOS386)
00325   #ifndef USE_ASSEMBLER
00326         memcpy((char*)ARRAY_MEMBLOCK_BASE,(char*)ARRAY_MEMBLOCK_SAVE,
00327           bytes_needed);
00328   #else
00329          dw_block_move((double*)ARRAY_MEMBLOCK_BASE,
00330            (double*)ARRAY_MEMBLOCK_SAVE,bytes_needed/8);
00331   #endif
00332 #else
00333      unsigned long max_move=50000;
00334 
00335      long int left_to_move=bytes_needed;
00336      humungous_pointer src = ARRAY_MEMBLOCK_SAVE;
00337      humungous_pointer dest = ARRAY_MEMBLOCK_BASE;
00338      while(left_to_move > max_move)
00339      {
00340        memcpy((char*)dest,(char*)src,max_move);
00341        left_to_move-=max_move;
00342        src+=max_move;
00343        dest+=max_move;
00344      }
00345      memcpy((char*)dest,(char*)src,left_to_move);
00346 #endif
00347     ARRAY_MEMBLOCK_SAVE.free();
00348   }
00349   else
00350   {
00351     humungous_pointer dest = ARRAY_MEMBLOCK_BASE;
00352     lseek(gradient_structure::GRAD_STACK1->_VARSSAV_PTR,0L,SEEK_SET);
00353 #if defined(DOS386)
00354   #if defined(OPT_LIB) && !defined(_MSC_VER)
00355     read(gradient_structure::GRAD_STACK1->_VARSSAV_PTR,
00356       (char*)dest,bytes_needed);
00357   #else
00358     ssize_t ret = read(gradient_structure::GRAD_STACK1->_VARSSAV_PTR,
00359       (char*)dest,bytes_needed);
00360     assert(ret != -1);
00361   #endif
00362 #else
00363      unsigned long int max_move=50000;
00364 
00365      long int left_to_move=bytes_needed;
00366      while(left_to_move > max_move)
00367      {
00368        read(gradient_structure::GRAD_STACK1->_VARSSAV_PTR,
00369          (char*)dest,max_move);
00370        left_to_move-=max_move;
00371        dest+=max_move;
00372      }
00373      read(gradient_structure::GRAD_STACK1->_VARSSAV_PTR,
00374        (char*)dest,left_to_move);
00375 #endif
00376   }
00377 }
00381 void gradient_structure::save_variables()
00382 {
00383   if ((variables_save = new double[gradient_structure::MAX_DLINKS])==NULL)
00384   {
00385     //_VARSSAV_PTR=open(var_store_file_name, O_RDWR | O_CREAT | O_TRUNC ,
00386     //               S_IREAD | S_IWRITE);
00387     cerr << "insufficient memory to allocate space for dvariables"
00388          << " save buffer " << endl;
00389     ad_exit(1);
00390   }
00391   for (unsigned int i=0; i<gradient_structure::GRAD_LIST->nlinks; i++)
00392   {
00393     //variables_save[i]= * (double*)
00394     //    (gradient_structure::GRAD_LIST->dlink_addresses[i]);
00395     memcpy(&(variables_save[i]),
00396        gradient_structure::GRAD_LIST->dlink_addresses[i],sizeof(double));
00397   }
00398 }
00402 void gradient_structure::restore_variables()
00403 {
00404   for (unsigned int i=0; i<gradient_structure::GRAD_LIST->nlinks; i++)
00405   {
00406     //* (double*)(gradient_structure::GRAD_LIST->dlink_addresses[i])
00407     //  = variables_save[i];
00408     memcpy(gradient_structure::GRAD_LIST->dlink_addresses[i],
00409        &(variables_save[i]),sizeof(double));
00410   }
00411   delete [] variables_save;
00412   variables_save = 0;
00413 }
00417 void reset_gradient_stack(void)
00418 {
00419   gradient_structure::GRAD_STACK1->ptr =
00420     gradient_structure::GRAD_STACK1->ptr_first;
00421 
00422   int& _GRADFILE_PTR=gradient_structure::GRAD_STACK1->_GRADFILE_PTR;
00423 
00424   lseek(_GRADFILE_PTR,0L,SEEK_SET);
00425 }
00435 void grad_stack::set_gradient_stack1(void (* func)(void),
00436   double* dep_addr, double* ind_addr1)
00437 {
00438 #ifdef NO_DERIVS
00439   if (!gradient_structure::no_derivatives)
00440   {
00441 #endif
00442     if (ptr > ptr_last)
00443     {
00444        // current buffer is full -- write it to disk and reset pointer
00445        // and counter
00446        this->write_grad_stack_buffer();
00447     }
00448     //test_the_pointer();
00449     ptr->func = func;
00450     ptr->dep_addr = dep_addr;
00451     ptr->ind_addr1 = ind_addr1;
00452     ptr++;
00453 #ifdef NO_DERIVS
00454   }
00455 #endif
00456 }
00457 void test_the_pointer(void)
00458 {
00459 /*
00460   static int inner_count = 0;
00461   static grad_stack_entry * pgse = (grad_stack_entry*) (0x1498fffc);
00462   inner_count++;
00463   if (inner_count == 404849)
00464   {
00465     cout << ptr << endl;
00466     cout << ptr->func << endl;
00467     cout << ptr->dep_addr << endl;
00468     cout << (int)(ptr->dep_addr)%8 << endl;
00469   }
00470   pgse->func = (void (*)())(100);
00471   pgse->dep_addr = (double*) 100;
00472   pgse->ind_addr1 = (double*) 100;
00473 */
00474 }