ADMB Documentation  11.1.2495
 All Classes Files Functions Variables Typedefs Friends Defines
xgradclc.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: xgradclc.cpp 2340 2014-09-13 10:22:50Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include "fvar.hpp"
00012 
00013 #include <sys/stat.h>
00014 #include <fcntl.h>
00015 #include <string.h>
00016 
00017 #ifdef __TURBOC__
00018   #pragma hdrstop
00019   #include <iostream.h>
00020 #endif
00021 
00022 #ifdef __ZTC__
00023   #include <iostream.hpp>
00024 #endif
00025 
00026 #if defined (__WAT32__)
00027 #  include <io.h>
00028 #endif
00029 
00030 #include <stdio.h>
00031 #include <stdlib.h>
00032 
00033 #ifdef _MSC_VER
00034   #define lseek _lseek
00035   #define  read _read
00036   #define write _write
00037 #else
00038   #include <iostream>
00039   using namespace std;
00040   #include <sys/stat.h>
00041   #include <sys/types.h>
00042   #include <unistd.h>
00043 #endif
00044 
00045 #ifdef __SUN__
00046   #include <iostream.h>
00047   #include <sys/stat.h>
00048   #include <sys/types.h>
00049   #include <unistd.h>
00050 #endif
00051 
00052 #if defined(__NDPX__ )
00053   extern "C" {
00054     int lseek(int, int, int);
00055     int read(int, char*, int);
00056   };
00057 #endif
00058 
00059 #include <math.h>
00060 
00061 void funnel_derivatives(void);
00062 
00063 #if defined(__ZTC__)
00064   void _far * _cdecl _farptr_norm(void _far *);
00065   void _far * _cdecl _farptr_fromlong(unsigned long);
00066   long _cdecl _farptr_tolong(void _far *);
00067 #endif
00068 
00073 void funnel_gradcalc(void)
00074 {
00075 #  ifdef NO_DERIVS
00076      if (gradient_structure::no_derivatives)
00077      {
00078        return;
00079      }
00080 #  endif
00081   gradient_structure::TOTAL_BYTES = 0;
00082   gradient_structure::PREVIOUS_TOTAL_BYTES=0;
00083   unsigned int i;
00084   if(!gradient_structure::instances)
00085   {
00086     return;
00087   }
00088 
00089    gradient_structure::GRAD_STACK1->_GRADFILE_PTR =
00090               gradient_structure::GRAD_STACK1->gradfile_handle();
00091 
00092   int& _GRADFILE_PTR=gradient_structure::GRAD_STACK1->_GRADFILE_PTR;
00093 
00094   off_t lpos = lseek(_GRADFILE_PTR,0L,SEEK_CUR);
00095 
00096   if(gradient_structure::GRAD_STACK1->ptr
00097        <= gradient_structure::GRAD_STACK1->ptr_first)
00098   {
00099 #ifdef SAFE_ALL
00100       cerr <<
00101         "warning -- calling funnel_gradcalc when no calculations generating"
00102            << endl << "derivative information have occurred" << endl;
00103 #endif
00104     return;
00105   }    // current is one past the end so -- it
00106 
00107   //if (gradient_structure::save_var_flag)
00108   {
00109     gradient_structure::save_arrays();
00110     gradient_structure::save_variables();
00111   }
00112 
00113   gradient_structure::GRAD_STACK1->ptr--;
00114 
00115   for (i=0; i<gradient_structure::GRAD_LIST->nlinks; i++)
00116   {
00117     * (double*) (gradient_structure::GRAD_LIST->dlink_addresses[i]) = 0;
00118   }
00119 
00120   double_and_int* tmp =
00121     (double_and_int*)gradient_structure::ARRAY_MEMBLOCK_BASE;
00122 
00123      unsigned long int max_last_offset =
00124        gradient_structure::ARR_LIST1->get_max_last_offset();
00125 
00126      unsigned int size = sizeof(double_and_int );
00127 
00128   double * zptr;
00129 
00130    for (i = 0; i < (max_last_offset/size); i++)
00131    {
00132      tmp->x = 0;
00133 #if defined (__ZTC__)
00134      tmp = (double_and_int*)_farptr_norm((void*)(++tmp));
00135 #else
00136      tmp++;
00137 #endif
00138    }
00139 
00140     *gradient_structure::GRAD_STACK1->ptr->dep_addr = 1;
00141     zptr = gradient_structure::GRAD_STACK1->ptr->dep_addr;
00142 
00143 int break_flag=1;
00144 int funnel_flag=0;
00145 
00146 do
00147 {
00148   gradient_structure::GRAD_STACK1->ptr++;
00149   #ifdef FAST_ASSEMBLER
00150     gradloop();
00151   #else
00152     grad_stack_entry * grad_ptr_first=
00153       gradient_structure::GRAD_STACK1->ptr_first;
00154     while (gradient_structure::GRAD_STACK1->ptr-- >
00155            grad_ptr_first)
00156     {
00157       if (!gradient_structure::GRAD_STACK1->ptr->func)
00158       {
00159         funnel_flag=1;
00160         break;
00161       }
00162       else
00163         (*(gradient_structure::GRAD_STACK1->ptr->func))();
00164     }
00165 
00166   #endif
00167    if (funnel_flag) break;
00168 
00169   // back up the file one buffer size and read forward
00170   lpos = lseek(gradient_structure::GRAD_STACK1->_GRADFILE_PTR,
00171       -((long int)(sizeof(grad_stack_entry)*gradient_structure::
00172         GRAD_STACK1->length)),SEEK_CUR);
00173 
00174   break_flag=gradient_structure::GRAD_STACK1->read_grad_stack_buffer(lpos);
00175 }  while (break_flag); // do
00176 
00177  {
00178    if (lpos<0)
00179    {
00180      #ifdef GRAD_DIAG
00181       long int ttmp =
00182      #endif
00183       lseek(gradient_structure::GRAD_STACK1->_GRADFILE_PTR, 0,SEEK_CUR);
00184 
00185      #ifdef GRAD_DIAG
00186       cout << "Offset in file at end of gradcalc is " << ttmp
00187            << " bytes from the beginning\n";
00188      #endif
00189    }
00190  }
00191 
00192   //if (gradient_structure::save_var_flag)
00193   {
00194     long bytes_needed=min(gradient_structure::ARR_LIST1->get_last_offset()+1,
00195       gradient_structure::ARRAY_MEMBLOCK_SIZE);
00196     unsigned int dsize=bytes_needed/sizeof(double);
00197     //dvector dtmp(0,dsize-1);
00198     //memcpy((char*)&(dtmp(0)),(char*)gradient_structure::ARRAY_MEMBLOCK_BASE,
00199       //dsize*sizeof(double));
00200 
00201     double* dptr=(double*) gradient_structure::ARRAY_MEMBLOCK_BASE;
00202     dptr-=1;
00203     unsigned int ii=0;
00204     int nzero=0;
00205     int nnzero=0;
00206     int dcount=0;
00207     int zero_flag=0;
00208     ivector offset(0,dsize-1);
00209     save_identifier_string("ue");
00210     if (*(++dptr))
00211     {
00212       save_double_value(*dptr);
00213       dcount++;
00214       zero_flag=0;
00215       offset(ii++)=0;
00216       nnzero++;
00217     }
00218     else
00219     {
00220       zero_flag=1;
00221       nzero++;
00222     }
00223 
00224     for (unsigned int i1=1;i1<dsize;i1++)
00225     {
00226       if (*(++dptr))
00227       {
00228         save_double_value(*dptr);
00229         dcount++;
00230         nnzero++;
00231         if (zero_flag)
00232         {
00233           offset(ii++)=nzero;
00234           nzero=0;
00235           zero_flag=0;
00236         }
00237       }
00238       else
00239       {
00240         nzero++;
00241         if (!zero_flag)
00242         {
00243           offset(ii++)=nnzero;
00244           nnzero=0;
00245           zero_flag=1;
00246         }
00247       }
00248     }
00249     save_int_value(dcount);
00250 
00251     for (i=0;i<ii;i++)
00252     {
00253       save_int_value(offset(i));
00254     }
00255     save_int_value(ii);
00256 
00257     unsigned int ssize=gradient_structure::GRAD_LIST->nlinks;
00258     dvector stmp(0,ssize-1);
00259 
00260     for (i=0; i<gradient_structure::GRAD_LIST->nlinks; i++)
00261     {
00262       memcpy((char*)&(stmp(i)),
00263         gradient_structure::GRAD_LIST->dlink_addresses[i],sizeof(double));
00264     }
00265     //dtmp.save_dvector_value();
00266     //dtmp.save_dvector_position();
00267     stmp.save_dvector_value();
00268     stmp.save_dvector_position();
00269 
00270     // save the address of the dependent variable for the funnel
00271     int wsize=sizeof(double_and_int*);
00272     gradient_structure::get_fp()->fwrite(&zptr,size_t(wsize));
00273     save_identifier_string("ae");
00274 
00275     gradient_structure::GRAD_STACK1->set_gradient_stack(funnel_derivatives);
00276     gradient_structure::restore_arrays();
00277     gradient_structure::restore_variables();
00278   }
00279 }
00280 
00285 void funnel_derivatives(void)
00286 {
00287   verify_identifier_string("ae");
00288   prevariable_position deppos=restore_prevariable_position();
00289   dvector_position stmp_pos=restore_dvector_position();
00290   dvector stmp=restore_dvector_value(stmp_pos);
00291   //dvector_position dtmp_pos=restore_dvector_position();
00292   //dvector dtmp=restore_dvector_value(dtmp_pos);
00293   int ii=restore_int_value();
00294   int i;
00295   int ip=ii;
00296   if (!ip) ip=1;
00297   ivector offset(0,ip);
00298   offset(ip)=0;
00299   //ivector offset(0,ip-1);
00300   for (i=ii-1;i>=0;i--)
00301   {
00302     offset(i)=restore_int_value();
00303   }
00304   int dcount=restore_int_value();
00305   int dc=dcount;
00306   if (!dc) dc=1;
00307   dvector dx(0,dc-1);
00308   for (i=dcount-1;i>=0;i--)
00309   {
00310     dx(i)=restore_double_value();
00311   }
00312 
00313   verify_identifier_string("ue");
00314 
00315   double df = restore_prevariable_derivative(deppos);
00316   double * dptr= (double *) gradient_structure::ARRAY_MEMBLOCK_BASE;
00317 
00318   //double * dd = &(dx(1));
00319   ii=0;
00320   int ic=0;
00321   dptr+=offset(ii++);
00322   for (i=0;i<dcount;i++)
00323   {
00324     *(dptr++)+=dx(i)*df;
00325     if (++ic==offset(ii))
00326     {
00327       if (i==dcount-1)
00328       {
00329         break;
00330       }
00331       dptr+=offset(ii+1);
00332       ii+=2;
00333       ic=0;
00334     }
00335   }
00336 
00337   int smax=stmp.indexmax();
00338   for (i=0;i<smax;i++)
00339   {
00340     if (stmp(i))
00341     {
00342       *(double*)(gradient_structure::GRAD_LIST->dlink_addresses[i])
00343         +=stmp(i)*df;
00344     }
00345   }
00346 }
00347 
00352 dvariable& funnel_dvariable::operator=(const prevariable& t)
00353 {
00354   dvariable::operator = (t);
00355   funnel_gradcalc();
00356   return *this;
00357 }
00358 
00363 void ad_begin_funnel(void)
00364 {
00365   gradient_structure::GRAD_STACK1->set_gradient_stack(NULL);
00366 }
00367