ADMB Documentation  11.1.1932
 All Classes Files Functions Variables Typedefs Friends Defines
jacob2.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: jacob2.cpp 1908 2014-04-17 22:45:07Z 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 #include <stdio.h>
00023 #include <stdlib.h>
00024 
00025 #ifdef __SUN__
00026   #include <iostream.h>
00027   #include <sys/stat.h>
00028   #include <sys/types.h>
00029   #include <unistd.h>
00030   #ifdef _MSC_VER
00031     #define lseek _lseek
00032     #define  read _read
00033     #define write _write
00034     #define open _open
00035     #define close _close
00036   #endif
00037 #endif
00038 
00039 #if !defined(_MSC_VER)
00040   #include <iostream>
00041   using namespace std;
00042   #include <sys/stat.h>
00043   #include <sys/types.h>
00044   #include <unistd.h>
00045 #endif
00046 
00047 #if defined(__NDPX__ )
00048   extern "C" {
00049     int lseek(int, int, int);
00050     int read(int, char*, int);
00051   };
00052 #endif
00053 
00054 #include <math.h>
00055 
00060 void jacobcalc(int nvar, const ofstream& ofs)
00061 {
00062   gradient_structure::jacobcalc(nvar,ofs);
00063 }
00064 
00069 void gradient_structure::jacobcalc(int nvar, const ofstream& _ofs)
00070 {
00071   ADUNCONST(ofstream,ofs)
00072   dvector jac(1,nvar);
00073   unsigned int i;
00074   my_off_t lpos;
00075   int depvar_count=DEPVARS_INFO->depvar_count;
00076 
00077   int& _GRADFILE_PTR=GRAD_STACK1->_GRADFILE_PTR;
00078   // check to see if anything has been written into the file
00079   my_off_t last_gpos=lseek(_GRADFILE_PTR,0L,SEEK_CUR);
00080 
00081   //save current contents of the buffer so we can get them later
00082   if (last_gpos)
00083   {
00084     GRAD_STACK1->write_grad_stack_buffer();
00085   }
00086 
00087   // check to see if anything has been written into the file
00088   my_off_t last_cpos=lseek(fp->file_ptr,0L,SEEK_CUR);
00089 
00090   //save current contents of the buffer so we can get them later
00091   if (last_cpos)
00092   {
00093     fp->write_cmpdif_stack_buffer();
00094   }
00095 
00096   // save variable values if desired
00097   if (save_var_flag)
00098   {
00099     save_arrays();
00100     save_variables();
00101   }
00102   // now evalueate the jacobian
00103   for (int ijac=1;ijac<=depvar_count;ijac++)
00104   {
00105     dvector& g=jac;
00106     //max_num_dependent_variables=ndv;
00107     if (depvar_count>DEPVARS_INFO->max_num_dependent_variables)
00108     {
00109       cout << "maximum number of depdendent variables of "
00110          << DEPVARS_INFO->max_num_dependent_variables << " exceeded "
00111          << endl
00112          << "use gradient_structure::set_NUM_DEPENDENT_VARIABLES(int i);"
00113          << endl << "to increase the number of dependent variables"
00114          << endl;
00115       ad_exit(1);
00116     }
00117 
00118     fp->offset=DEPVARS_INFO->cmpdif_buffer_position(ijac);
00119     fp->toffset=fp->offset;
00120     _GRADFILE_PTR=DEPVARS_INFO->grad_file_count(ijac);
00121     fp->file_ptr=DEPVARS_INFO->cmpdif_file_count(ijac);
00122     lpos=DEPVARS_INFO->grad_file_position(ijac);
00123     // position the cmpdif file correctly;
00124     if (last_cpos)
00125     {
00126       my_off_t cmp_lpos=DEPVARS_INFO->cmpdif_file_position(ijac);
00127       lseek(fp->file_ptr,cmp_lpos,SEEK_SET);
00128       fp->read_cmpdif_stack_buffer(cmp_lpos);
00129     }
00130     GRAD_STACK1->_GRADFILE_PTR = GRAD_STACK1->gradfile_handle();
00131 
00132     if (last_gpos)
00133     {
00134       // just use the end of the buffer
00135       GRAD_STACK1->set_gbuffer_pointers();
00136 
00137       // check to sere if buffer was written into the beginning of
00138       // the next file
00139       if ( (GRAD_STACK1->_GRADFILE_PTR == GRAD_STACK1->_GRADFILE_PTR1)
00140          && (lpos == GRAD_STACK1->end_pos1) && (lpos>0) )
00141       {
00142         // get the next file
00143         GRAD_STACK1->increment_current_gradfile_ptr();
00144         lpos=0;
00145       }
00146       // and position the file to the begginig of the buffer image
00147       lseek(_GRADFILE_PTR,lpos,SEEK_SET);
00148       // now fill the buffer with the appropriate stuff
00149       GRAD_STACK1->read_grad_stack_buffer(lpos);
00150       // now reposition the grad_buffer pointer
00151     }
00152     GRAD_STACK1->ptr=
00153          (grad_stack_entry *)DEPVARS_INFO->grad_buffer_position(ijac);
00154 
00155     if(GRAD_STACK1->ptr <= GRAD_STACK1->ptr_first)
00156     {
00157       #ifdef SAFE_ARRAYS
00158         cerr << "warning -- calling gradcalc when no calculations generating"
00159          << endl << "derivative information have occurred" << endl;
00160       #endif
00161       g.initialize();
00162       return;
00163     }    // current is one past the end so -- it
00164 
00165     gradient_structure::GRAD_STACK1->ptr--;
00166 
00167     for (i=0; i< GRAD_LIST->nlinks; i++)
00168     {
00169       * (double*) (GRAD_LIST->dlink_addresses[i]) = 0;
00170     }
00171 
00172     double_and_int* tmp =
00173       (double_and_int*)gradient_structure::ARRAY_MEMBLOCK_BASE;
00174 
00175     unsigned long int max_last_offset
00176                = gradient_structure::ARR_LIST1->get_max_last_offset();
00177 
00178     unsigned int size = sizeof(double_and_int );
00179 
00180     for (i = 0; i < (max_last_offset/size); i++)
00181     {
00182       tmp->x = 0;
00183 #if defined (__ZTC__)
00184       tmp = (double_and_int*)_farptr_norm((void*)(++tmp));
00185 #else
00186       tmp++;
00187 #endif
00188     }
00189 
00190     * gradient_structure::GRAD_STACK1->ptr->dep_addr  = 1;
00191     //double* zptr = gradient_structure::GRAD_STACK1->ptr->dep_addr;
00192 
00193     int break_flag=1;
00194 
00195     do
00196     {
00197       gradient_structure::GRAD_STACK1->ptr++;
00198       #ifdef FAST_ASSEMBLER
00199         gradloop();
00200       #else
00201         //int counter=0;
00202       while (gradient_structure::GRAD_STACK1->ptr-- >
00203              gradient_structure::GRAD_STACK1->ptr_first)
00204       {
00205         //grad_stack_entry* grad_ptr =
00206         //gradient_structure::GRAD_STACK1->ptr;
00207         {
00208           (* gradient_structure::GRAD_STACK1->ptr->func)();
00209         }
00210       }
00211       #endif
00212 
00213   // back up the file one buffer size and read forward
00214       lpos = lseek(gradient_structure::GRAD_STACK1->_GRADFILE_PTR,
00215         -((long int)(sizeof(grad_stack_entry)*gradient_structure::
00216         GRAD_STACK1->length)),SEEK_CUR);
00217 
00218        break_flag=gradient_structure::
00219                   GRAD_STACK1->read_grad_stack_buffer(lpos);
00220     }  while (break_flag); // do
00221 
00222     int mindx = g.indexmin();
00223     for (i=0; i<(unsigned int)nvar; i++)
00224     {
00225       g[i+mindx] =  * gradient_structure::INDVAR_LIST->get_address(i);
00226     }
00227     gradient_structure::GRAD_STACK1->ptr =
00228          gradient_structure::GRAD_STACK1->ptr_first;
00229     //ofs << setprecision(10) << g << endl;
00230     ofs.precision(10);
00231     ofs << g << endl;
00232   }// loop over dep vars
00233   DEPVARS_INFO->depvar_count=0;
00234   if (gradient_structure::save_var_flag)
00235   {
00236     gradient_structure::restore_arrays();
00237     gradient_structure::restore_variables();
00238   }
00239 }