ADMB Documentation  11.1.2490
 All Classes Files Functions Variables Typedefs Friends Defines
jacobclc.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: jacobclc.cpp 2336 2014-09-13 08:12:31Z 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 __SUN__
00034   #include <iostream.h>
00035   #include <sys/stat.h>
00036   #include <sys/types.h>
00037   #include <unistd.h>
00038 #endif
00039 
00040 #ifdef _MSC_VER
00041   #define lseek _lseek
00042   #define  read _read
00043   #define write _write
00044   #define open _open
00045   #define close _close
00046 #else
00047   #include <iostream>
00048   using namespace std;
00049   #include <sys/stat.h>
00050   #include <sys/types.h>
00051   #include <unistd.h>
00052 #endif
00053 
00054 #if defined(__NDPX__ )
00055   extern "C" {
00056     int lseek(int, int, int);
00057     int read(int, char*, int);
00058   };
00059 #endif
00060 
00061 #if defined(__ZTC__)
00062   void _far * _cdecl _farptr_norm(void _far *);
00063   void _far * _cdecl _farptr_fromlong(unsigned long);
00064   long _cdecl _farptr_tolong(void _far *);
00065 #endif
00066 
00067 #include <math.h>
00068 
00073 void jacobcalc(int nvar, const dmatrix& jac)
00074 {
00075   gradient_structure::jacobcalc(nvar,jac);
00076 }
00077 
00082 void gradient_structure::jacobcalc(int nvar, const dmatrix& _jac)
00083 {
00084   ADUNCONST(dmatrix,jac)
00085 
00086   unsigned int i;
00087   off_t lpos;
00088   if(!instances)
00089   {
00090     jac.initialize();
00091     return;
00092   }
00093   if (jac.rowmin()!=1)
00094   {
00095     cerr << "Error in jacobcalc jacobian must have minimum valid"
00096             " index of 1" << endl;
00097     jac.initialize();
00098     return;
00099   }
00100   int depvar_count=DEPVARS_INFO->depvar_count;
00101   if (jac.rowmax()<depvar_count)
00102   {
00103     cerr << "Error in jacobcalc jacobian must have maximumvalid"
00104             " index >= the number of dependent variables " << endl
00105           << " which is " << depvar_count << endl;
00106     jac.initialize();
00107     return;
00108   }
00109 
00110   int& _GRADFILE_PTR=GRAD_STACK1->_GRADFILE_PTR;
00111   // check to see if anything has been written into the file
00112   off_t last_gpos=lseek(_GRADFILE_PTR,0L,SEEK_CUR);
00113 
00114   //save current contents of the buffer so we can get them later
00115   if (last_gpos)
00116   {
00117     GRAD_STACK1->write_grad_stack_buffer();
00118   }
00119 
00120   // check to see if anything has been written into the file
00121   off_t last_cpos=lseek(fp->file_ptr,0L,SEEK_CUR);
00122 
00123   //save current contents of the buffer so we can get them later
00124   if (last_cpos)
00125   {
00126     fp->write_cmpdif_stack_buffer();
00127   }
00128 
00129   for (i=jac.rowmin();i<=(unsigned int)jac.rowmax();i++)
00130   {
00131     if (jac(i).indexmin() !=1)
00132     {
00133       cerr  << "jacobian matrix minimum column index must equal 1"
00134                " for all rows" << endl;
00135       return;
00136     }
00137     if (jac(i).indexmax() < nvar)
00138     {
00139       cerr  << "jacobian matrix column size is less than the number of"
00140                " independent variables" << endl;
00141       return;
00142     }
00143   }
00144   // save variable values if desired
00145   if (save_var_flag)
00146   {
00147     save_arrays();
00148     save_variables();
00149   }
00150   // now evalueate the jacobian
00151   for (int ijac=1;ijac<=depvar_count;ijac++)
00152   {
00153     dvector& g=(dvector&)(jac(ijac));
00154     //max_num_dependent_variables=ndv;
00155     if (depvar_count>DEPVARS_INFO->max_num_dependent_variables)
00156     {
00157       cout << "maximum number of depdendent variables of "
00158          << DEPVARS_INFO->max_num_dependent_variables << " exceeded "
00159          << endl
00160          << "use gradient_structure::set_NUM_DEPENDENT_VARIABLES(int i);"
00161          << endl << "to increase the number of dependent variables"
00162          << endl;
00163       ad_exit(1);
00164     }
00165 
00166     fp->offset=DEPVARS_INFO->cmpdif_buffer_position(ijac);
00167     fp->toffset=fp->offset;
00168     _GRADFILE_PTR=DEPVARS_INFO->grad_file_count(ijac);
00169     fp->file_ptr=DEPVARS_INFO->cmpdif_file_count(ijac);
00170     lpos=DEPVARS_INFO->grad_file_position(ijac);
00171     // position the cmpdif file correctly;
00172     if (last_cpos)
00173     {
00174       off_t cmp_lpos=DEPVARS_INFO->cmpdif_file_position(ijac);
00175       lseek(fp->file_ptr,cmp_lpos,SEEK_SET);
00176       fp->read_cmpdif_stack_buffer(cmp_lpos);
00177     }
00178     GRAD_STACK1->_GRADFILE_PTR = GRAD_STACK1->gradfile_handle();
00179 
00180     if (last_gpos)
00181     {
00182       // just use the end of the buffer
00183       GRAD_STACK1->set_gbuffer_pointers();
00184 
00185       // check to sere if buffer was written into the beginning of
00186       // the next file
00187       if ( (GRAD_STACK1->_GRADFILE_PTR == GRAD_STACK1->_GRADFILE_PTR1)
00188          && (lpos == GRAD_STACK1->end_pos1) && (lpos>0) )
00189       {
00190         // get the next file
00191         GRAD_STACK1->increment_current_gradfile_ptr();
00192         lpos=0;
00193       }
00194       // and position the file to the begginig of the buffer image
00195       lseek(_GRADFILE_PTR,lpos,SEEK_SET);
00196       // now fill the buffer with the appropriate stuff
00197       GRAD_STACK1->read_grad_stack_buffer(lpos);
00198       // now reposition the grad_buffer pointer
00199     }
00200     GRAD_STACK1->ptr=
00201          (grad_stack_entry *)DEPVARS_INFO->grad_buffer_position(ijac);
00202 
00203     if(GRAD_STACK1->ptr <= GRAD_STACK1->ptr_first)
00204     {
00205 #ifdef SAFE_ALL
00206         cerr << "warning -- calling gradcalc when no calculations generating"
00207          << endl << "derivative information have occurred" << endl;
00208 #endif
00209       g.initialize();
00210       return;
00211     }    // current is one past the end so -- it
00212 
00213     gradient_structure::GRAD_STACK1->ptr--;
00214 
00215     for (i=0; i< GRAD_LIST->nlinks; i++)
00216     {
00217       * (double*) (GRAD_LIST->dlink_addresses[i]) = 0;
00218     }
00219 
00220     double_and_int* tmp =
00221       (double_and_int*)gradient_structure::ARRAY_MEMBLOCK_BASE;
00222 
00223     unsigned long int max_last_offset
00224                = gradient_structure::ARR_LIST1->get_max_last_offset();
00225 
00226     unsigned int size = sizeof(double_and_int);
00227 
00228     for (i = 0; i < (max_last_offset/size); i++)
00229     {
00230       tmp->x = 0;
00231 #if defined (__ZTC__)
00232       tmp = (double_and_int*)_farptr_norm((void*)(++tmp));
00233 #else
00234       tmp++;
00235 #endif
00236     }
00237 
00238     *gradient_structure::GRAD_STACK1->ptr->dep_addr  = 1;
00239 
00240     int break_flag=1;
00241     do
00242     {
00243       gradient_structure::GRAD_STACK1->ptr++;
00244       #ifdef FAST_ASSEMBLER
00245         gradloop();
00246       #else
00247         //int counter=0;
00248       while (gradient_structure::GRAD_STACK1->ptr-- >
00249              gradient_structure::GRAD_STACK1->ptr_first)
00250       {
00251         //grad_stack_entry* grad_ptr =
00252         //gradient_structure::GRAD_STACK1->ptr;
00253         {
00254           (* gradient_structure::GRAD_STACK1->ptr->func)();
00255         }
00256       }
00257       #endif
00258 
00259   // back up the file one buffer size and read forward
00260       lpos = lseek(gradient_structure::GRAD_STACK1->_GRADFILE_PTR,
00261         -((long int)(sizeof(grad_stack_entry)*gradient_structure::
00262         GRAD_STACK1->length)),SEEK_CUR);
00263 
00264        break_flag=gradient_structure::
00265                   GRAD_STACK1->read_grad_stack_buffer(lpos);
00266     }  while (break_flag); // do
00267 
00268     int mindx = g.indexmin();
00269     dvector & gg=(dvector&)(g);
00270     for (i=0; i<(unsigned int)nvar; i++)
00271     {
00272       gg[i+mindx] =  * gradient_structure::INDVAR_LIST->get_address(i);
00273       //g[i+mindx] =  * gradient_structure::INDVAR_LIST->get_address(i);
00274     }
00275     gradient_structure::GRAD_STACK1->ptr =
00276          gradient_structure::GRAD_STACK1->ptr_first;
00277   }// loop over dep vars
00278   DEPVARS_INFO->depvar_count=0;
00279   if (gradient_structure::save_var_flag)
00280   {
00281     gradient_structure::restore_arrays();
00282     gradient_structure::restore_variables();
00283   }
00284 }