ADMB Documentation  11.1.1927
 All Classes Files Functions Variables Typedefs Friends Defines
xgradclc.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: xgradclc.cpp 1905 2014-04-17 20:13:56Z 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   long int lpos;
00085   if(!gradient_structure::instances)
00086   {
00087     return;
00088   }
00089 
00090    gradient_structure::GRAD_STACK1->_GRADFILE_PTR =
00091               gradient_structure::GRAD_STACK1->gradfile_handle();
00092 
00093   int& _GRADFILE_PTR=gradient_structure::GRAD_STACK1->_GRADFILE_PTR;
00094 
00095   lpos = lseek(_GRADFILE_PTR,0L,SEEK_CUR);
00096 
00097   if(gradient_structure::GRAD_STACK1->ptr
00098        <= gradient_structure::GRAD_STACK1->ptr_first)
00099   {
00100     #ifdef SAFE_ARRAYS
00101       cerr <<
00102         "warning -- calling funnel_gradcalc when no calculations generating"
00103            << endl << "derivative information have occurred" << endl;
00104     #endif
00105     return;
00106   }    // current is one past the end so -- it
00107 
00108   //if (gradient_structure::save_var_flag)
00109   {
00110     gradient_structure::save_arrays();
00111     gradient_structure::save_variables();
00112   }
00113 
00114   gradient_structure::GRAD_STACK1->ptr--;
00115 
00116   for (i=0; i<gradient_structure::GRAD_LIST->nlinks; i++)
00117   {
00118     * (double*) (gradient_structure::GRAD_LIST->dlink_addresses[i]) = 0;
00119   }
00120 
00121   double_and_int* tmp =
00122     (double_and_int*)gradient_structure::ARRAY_MEMBLOCK_BASE;
00123 
00124      unsigned long int max_last_offset =
00125        gradient_structure::ARR_LIST1->get_max_last_offset();
00126 
00127      unsigned int size = sizeof(double_and_int );
00128 
00129   double * zptr;
00130 
00131    for (i = 0; i < (max_last_offset/size); i++)
00132    {
00133      tmp->x = 0;
00134 #if defined (__ZTC__)
00135      tmp = (double_and_int*)_farptr_norm((void*)(++tmp));
00136 #else
00137      tmp++;
00138 #endif
00139    }
00140 
00141     *gradient_structure::GRAD_STACK1->ptr->dep_addr = 1;
00142     zptr = gradient_structure::GRAD_STACK1->ptr->dep_addr;
00143 
00144 int break_flag=1;
00145 int funnel_flag=0;
00146 
00147 do
00148 {
00149   gradient_structure::GRAD_STACK1->ptr++;
00150   #ifdef FAST_ASSEMBLER
00151     gradloop();
00152   #else
00153     grad_stack_entry * grad_ptr_first=
00154       gradient_structure::GRAD_STACK1->ptr_first;
00155     while (gradient_structure::GRAD_STACK1->ptr-- >
00156            grad_ptr_first)
00157     {
00158       if (!gradient_structure::GRAD_STACK1->ptr->func)
00159       {
00160         funnel_flag=1;
00161         break;
00162       }
00163       else
00164         (*(gradient_structure::GRAD_STACK1->ptr->func))();
00165     }
00166 
00167   #endif
00168    if (funnel_flag) break;
00169 
00170   // back up the file one buffer size and read forward
00171   lpos = lseek(gradient_structure::GRAD_STACK1->_GRADFILE_PTR,
00172       -((long int)(sizeof(grad_stack_entry)*gradient_structure::
00173         GRAD_STACK1->length)),SEEK_CUR);
00174 
00175   break_flag=gradient_structure::GRAD_STACK1->read_grad_stack_buffer(lpos);
00176 }  while (break_flag); // do
00177 
00178  {
00179    if (lpos<0)
00180    {
00181      #ifdef GRAD_DIAG
00182       long int ttmp =
00183      #endif
00184       lseek(gradient_structure::GRAD_STACK1->_GRADFILE_PTR, 0,SEEK_CUR);
00185 
00186      #ifdef GRAD_DIAG
00187       cout << "Offset in file at end of gradcalc is " << ttmp
00188            << " bytes from the beginning\n";
00189      #endif
00190    }
00191  }
00192 
00193   //if (gradient_structure::save_var_flag)
00194   {
00195     long bytes_needed=min(gradient_structure::ARR_LIST1->get_last_offset()+1,
00196       gradient_structure::ARRAY_MEMBLOCK_SIZE);
00197     unsigned int dsize=bytes_needed/sizeof(double);
00198     //dvector dtmp(0,dsize-1);
00199     //memcpy((char*)&(dtmp(0)),(char*)gradient_structure::ARRAY_MEMBLOCK_BASE,
00200       //dsize*sizeof(double));
00201 
00202     double* dptr=(double*) gradient_structure::ARRAY_MEMBLOCK_BASE;
00203     dptr-=1;
00204     unsigned int ii=0;
00205     int nzero=0;
00206     int nnzero=0;
00207     int dcount=0;
00208     int zero_flag=0;
00209     ivector offset(0,dsize-1);
00210     save_identifier_string("ue");
00211     if (*(++dptr))
00212     {
00213       save_double_value(*dptr);
00214       dcount++;
00215       zero_flag=0;
00216       offset(ii++)=0;
00217       nnzero++;
00218     }
00219     else
00220     {
00221       zero_flag=1;
00222       nzero++;
00223     }
00224 
00225     for (unsigned int i1=1;i1<dsize;i1++)
00226     {
00227       if (*(++dptr))
00228       {
00229         save_double_value(*dptr);
00230         dcount++;
00231         nnzero++;
00232         if (zero_flag)
00233         {
00234           offset(ii++)=nzero;
00235           nzero=0;
00236           zero_flag=0;
00237         }
00238       }
00239       else
00240       {
00241         nzero++;
00242         if (!zero_flag)
00243         {
00244           offset(ii++)=nnzero;
00245           nnzero=0;
00246           zero_flag=1;
00247         }
00248       }
00249     }
00250     save_int_value(dcount);
00251 
00252     for (i=0;i<ii;i++)
00253     {
00254       save_int_value(offset(i));
00255     }
00256     save_int_value(ii);
00257 
00258     unsigned int ssize=gradient_structure::GRAD_LIST->nlinks;
00259     dvector stmp(0,ssize-1);
00260 
00261     for (i=0; i<gradient_structure::GRAD_LIST->nlinks; i++)
00262     {
00263       memcpy((char*)&(stmp(i)),
00264         gradient_structure::GRAD_LIST->dlink_addresses[i],sizeof(double));
00265     }
00266     //dtmp.save_dvector_value();
00267     //dtmp.save_dvector_position();
00268     stmp.save_dvector_value();
00269     stmp.save_dvector_position();
00270 
00271     // save the address of the dependent variable for the funnel
00272     int wsize=sizeof(double_and_int*);
00273     gradient_structure::get_fp()->fwrite(&zptr,size_t(wsize));
00274     save_identifier_string("ae");
00275 
00276     gradient_structure::GRAD_STACK1->set_gradient_stack(funnel_derivatives);
00277     gradient_structure::restore_arrays();
00278     gradient_structure::restore_variables();
00279   }
00280 }
00281 
00286 void funnel_derivatives(void)
00287 {
00288   verify_identifier_string("ae");
00289   prevariable_position deppos=restore_prevariable_position();
00290   dvector_position stmp_pos=restore_dvector_position();
00291   dvector stmp=restore_dvector_value(stmp_pos);
00292   //dvector_position dtmp_pos=restore_dvector_position();
00293   //dvector dtmp=restore_dvector_value(dtmp_pos);
00294   int ii=restore_int_value();
00295   int i;
00296   int ip=ii;
00297   if (!ip) ip=1;
00298   ivector offset(0,ip);
00299   offset(ip)=0;
00300   //ivector offset(0,ip-1);
00301   for (i=ii-1;i>=0;i--)
00302   {
00303     offset(i)=restore_int_value();
00304   }
00305   int dcount=restore_int_value();
00306   int dc=dcount;
00307   if (!dc) dc=1;
00308   dvector dx(0,dc-1);
00309   for (i=dcount-1;i>=0;i--)
00310   {
00311     dx(i)=restore_double_value();
00312   }
00313 
00314   verify_identifier_string("ue");
00315 
00316   double df = restore_prevariable_derivative(deppos);
00317   double * dptr= (double *) gradient_structure::ARRAY_MEMBLOCK_BASE;
00318 
00319   //double * dd = &(dx(1));
00320   ii=0;
00321   int ic=0;
00322   dptr+=offset(ii++);
00323   for (i=0;i<dcount;i++)
00324   {
00325     *(dptr++)+=dx(i)*df;
00326     if (++ic==offset(ii))
00327     {
00328       if (i==dcount-1)
00329       {
00330         break;
00331       }
00332       dptr+=offset(ii+1);
00333       ii+=2;
00334       ic=0;
00335     }
00336   }
00337 
00338   int smax=stmp.indexmax();
00339   for (i=0;i<smax;i++)
00340   {
00341     if (stmp(i))
00342     {
00343       *(double*)(gradient_structure::GRAD_LIST->dlink_addresses[i])
00344         +=stmp(i)*df;
00345     }
00346   }
00347 }
00348 
00353 dvariable& funnel_dvariable::operator=(const prevariable& t)
00354 {
00355   dvariable::operator = (t);
00356   funnel_gradcalc();
00357   return *this;
00358 }
00359 
00364 void ad_begin_funnel(void)
00365 {
00366   gradient_structure::GRAD_STACK1->set_gradient_stack(NULL);
00367 }
00368