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