ADMB Documentation  11.1.2535
 All Classes Files Functions Variables Typedefs Friends Defines
sgradclc.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: sgradclc.cpp 2499 2014-10-25 00:19:55Z 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 offset = 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   int mindx = g.indexmin();
00206   for (i=0; i < (unsigned int)nvar; i++)
00207   {
00208     g[i + mindx] = *gradient_structure::INDVAR_LIST->get_address(i);
00209   }
00210 
00211   gradient_structure::GRAD_STACK1->ptr =
00212     gradient_structure::GRAD_STACK1->ptr_first;
00213 
00214   if (gradient_structure::save_var_flag)
00215   {
00216     gradient_structure::restore_arrays();
00217     gradient_structure::restore_variables();
00218   }
00219 }
00228 double gradcalc(int nvar, const dvector& _g, dvariable& f)
00229 {
00230   double v = value(f);
00231   gradcalc(nvar, _g);
00232   return v;
00233 }
00236 void gradient_structure::save_arrays()
00237 {
00238   void * temp_ptr;
00239   long bytes_needed=min(gradient_structure::ARR_LIST1->get_last_offset()+1,
00240     ARRAY_MEMBLOCK_SIZE);
00241   gradient_structure::save_var_file_flag=0;
00242 #ifdef __ZTC__
00243    if ( (temp_ptr = farmalloc(bytes_needed) ) == 0)
00244 #else
00245    //if ( (temp_ptr = malloc(bytes_needed) ) == 0)
00246    if ((temp_ptr = (void *)malloc(bytes_needed)) == 0)
00247   #define __USE_IOSTREAM__
00248 #endif
00249    {
00250      gradient_structure::save_var_file_flag=1;
00251      cerr << "insufficient memory to allocate space for ARRAY_MEMBLOCK"
00252           << " save buffer " << endl;
00253    }
00254    if (gradient_structure::save_var_file_flag==0)
00255    {
00256      ARRAY_MEMBLOCK_SAVE = temp_ptr;
00257 #if defined(DOS386)
00258   #ifndef USE_ASSEMBLER
00259          memcpy((char*)ARRAY_MEMBLOCK_SAVE,(char*)ARRAY_MEMBLOCK_BASE,
00260            bytes_needed);
00261   #else
00262          dw_block_move((double*)ARRAY_MEMBLOCK_SAVE,
00263            (double*)ARRAY_MEMBLOCK_BASE,bytes_needed/8);
00264   #endif
00265 #else
00266      unsigned long int max_move=50000;
00267      unsigned long int left_to_move=bytes_needed;
00268      humungous_pointer dest = ARRAY_MEMBLOCK_SAVE;
00269      humungous_pointer src = ARRAY_MEMBLOCK_BASE;
00270      while(left_to_move > max_move)
00271      {
00272        memcpy((char*)dest,(char*)src,max_move);
00273        left_to_move-=max_move;
00274        src+=max_move;
00275        dest+=max_move;
00276      }
00277      memcpy((char*)dest,(char*)src,left_to_move);
00278 #endif
00279   }
00280   else
00281   {
00282      humungous_pointer src = ARRAY_MEMBLOCK_BASE;
00283      lseek(gradient_structure::GRAD_STACK1->_VARSSAV_PTR,0L,SEEK_SET);
00284 #if defined(DOS386)
00285   #ifdef OPT_LIB
00286      write(gradient_structure::GRAD_STACK1->_VARSSAV_PTR,
00287        (char*)src, bytes_needed);
00288   #else
00289      ssize_t ret = write(gradient_structure::GRAD_STACK1->_VARSSAV_PTR,
00290        (char*)src, bytes_needed);
00291      assert(ret != -1);
00292   #endif
00293 #else
00294      unsigned long int max_move=500;
00295      unsigned long int left_to_move=bytes_needed;
00296      while(left_to_move > max_move)
00297      {
00298        write(_VARSSAV_PTR,(char*)src,max_move);
00299        left_to_move-=max_move;
00300        src+=max_move;
00301      }
00302      write(gradient_structure::GRAD_STACK1->_VARSSAV_PTR,(char*)src,
00303        left_to_move);
00304 #endif
00305   }
00306 }
00307 
00310 void gradient_structure::restore_arrays()
00311 {
00312   long bytes_needed=min(gradient_structure::ARR_LIST1->get_last_offset()+1,
00313     ARRAY_MEMBLOCK_SIZE);
00314   if (gradient_structure::save_var_file_flag==0)
00315   {
00316 #if defined(DOS386)
00317   #ifndef USE_ASSEMBLER
00318         memcpy((char*)ARRAY_MEMBLOCK_BASE,(char*)ARRAY_MEMBLOCK_SAVE,
00319           bytes_needed);
00320   #else
00321          dw_block_move((double*)ARRAY_MEMBLOCK_BASE,
00322            (double*)ARRAY_MEMBLOCK_SAVE,bytes_needed/8);
00323   #endif
00324 #else
00325      unsigned long max_move=50000;
00326 
00327      long int left_to_move=bytes_needed;
00328      humungous_pointer src = ARRAY_MEMBLOCK_SAVE;
00329      humungous_pointer dest = ARRAY_MEMBLOCK_BASE;
00330      while(left_to_move > max_move)
00331      {
00332        memcpy((char*)dest,(char*)src,max_move);
00333        left_to_move-=max_move;
00334        src+=max_move;
00335        dest+=max_move;
00336      }
00337      memcpy((char*)dest,(char*)src,left_to_move);
00338 #endif
00339     ARRAY_MEMBLOCK_SAVE.free();
00340   }
00341   else
00342   {
00343     humungous_pointer dest = ARRAY_MEMBLOCK_BASE;
00344     lseek(gradient_structure::GRAD_STACK1->_VARSSAV_PTR,0L,SEEK_SET);
00345 #if defined(DOS386)
00346   #ifdef OPT_LIB
00347     read(gradient_structure::GRAD_STACK1->_VARSSAV_PTR,
00348       (char*)dest,bytes_needed);
00349   #else
00350     ssize_t ret = read(gradient_structure::GRAD_STACK1->_VARSSAV_PTR,
00351       (char*)dest,bytes_needed);
00352     assert(ret != -1);
00353   #endif
00354 #else
00355      unsigned long int max_move=50000;
00356 
00357      long int left_to_move=bytes_needed;
00358      while(left_to_move > max_move)
00359      {
00360        read(gradient_structure::GRAD_STACK1->_VARSSAV_PTR,
00361          (char*)dest,max_move);
00362        left_to_move-=max_move;
00363        dest+=max_move;
00364      }
00365      read(gradient_structure::GRAD_STACK1->_VARSSAV_PTR,
00366        (char*)dest,left_to_move);
00367 #endif
00368   }
00369 }
00373 void gradient_structure::save_variables()
00374 {
00375   if ((variables_save = new double[gradient_structure::MAX_DLINKS])==NULL)
00376   {
00377     //_VARSSAV_PTR=open(var_store_file_name, O_RDWR | O_CREAT | O_TRUNC ,
00378     //               S_IREAD | S_IWRITE);
00379     cerr << "insufficient memory to allocate space for dvariables"
00380          << " save buffer " << endl;
00381     ad_exit(1);
00382   }
00383   for (unsigned int i=0; i<gradient_structure::GRAD_LIST->nlinks; i++)
00384   {
00385     //variables_save[i]= * (double*)
00386     //    (gradient_structure::GRAD_LIST->dlink_addresses[i]);
00387     memcpy(&(variables_save[i]),
00388        gradient_structure::GRAD_LIST->dlink_addresses[i],sizeof(double));
00389   }
00390 }
00394 void gradient_structure::restore_variables()
00395 {
00396   for (unsigned int i=0; i<gradient_structure::GRAD_LIST->nlinks; i++)
00397   {
00398     //* (double*)(gradient_structure::GRAD_LIST->dlink_addresses[i])
00399     //  = variables_save[i];
00400     memcpy(gradient_structure::GRAD_LIST->dlink_addresses[i],
00401        &(variables_save[i]),sizeof(double));
00402   }
00403   delete [] variables_save;
00404   variables_save = 0;
00405 }
00409 void reset_gradient_stack(void)
00410 {
00411   gradient_structure::GRAD_STACK1->ptr =
00412     gradient_structure::GRAD_STACK1->ptr_first;
00413 
00414   int& _GRADFILE_PTR=gradient_structure::GRAD_STACK1->_GRADFILE_PTR;
00415 
00416   lseek(_GRADFILE_PTR,0L,SEEK_SET);
00417 }
00427 void grad_stack::set_gradient_stack1(void (* func)(void),
00428   double* dep_addr, double* ind_addr1)
00429 {
00430 #ifdef NO_DERIVS
00431   if (!gradient_structure::no_derivatives)
00432   {
00433 #endif
00434     if (ptr > ptr_last)
00435     {
00436        // current buffer is full -- write it to disk and reset pointer
00437        // and counter
00438        this->write_grad_stack_buffer();
00439     }
00440     //test_the_pointer();
00441     ptr->func = func;
00442     ptr->dep_addr = dep_addr;
00443     ptr->ind_addr1 = ind_addr1;
00444     ptr++;
00445 #ifdef NO_DERIVS
00446   }
00447 #endif
00448 }
00449 void test_the_pointer(void)
00450 {
00451 /*
00452   static int inner_count = 0;
00453   static grad_stack_entry * pgse = (grad_stack_entry*) (0x1498fffc);
00454   inner_count++;
00455   if (inner_count == 404849)
00456   {
00457     cout << ptr << endl;
00458     cout << ptr->func << endl;
00459     cout << ptr->dep_addr << endl;
00460     cout << (int)(ptr->dep_addr)%8 << endl;
00461   }
00462   pgse->func = (void (*)())(100);
00463   pgse->dep_addr = (double*) 100;
00464   pgse->ind_addr1 = (double*) 100;
00465 */
00466 }