ADMB Documentation  11.1x.2711
 All Classes Files Functions Variables Typedefs Friends Defines
jacobclc.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: jacobclc.cpp 2643 2014-11-13 21:00:02Z 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   off_t lpos;
00087   if(!instances)
00088   {
00089     jac.initialize();
00090     return;
00091   }
00092   if (jac.rowmin()!=1)
00093   {
00094     cerr << "Error in jacobcalc jacobian must have minimum valid"
00095             " index of 1" << endl;
00096     jac.initialize();
00097     return;
00098   }
00099   int depvar_count=DEPVARS_INFO->depvar_count;
00100   if (jac.rowmax()<depvar_count)
00101   {
00102     cerr << "Error in jacobcalc jacobian must have maximumvalid"
00103             " index >= the number of dependent variables " << endl
00104           << " which is " << depvar_count << endl;
00105     jac.initialize();
00106     return;
00107   }
00108 
00109   int& _GRADFILE_PTR=GRAD_STACK1->_GRADFILE_PTR;
00110   // check to see if anything has been written into the file
00111   off_t last_gpos=lseek(_GRADFILE_PTR,0L,SEEK_CUR);
00112 
00113   //save current contents of the buffer so we can get them later
00114   if (last_gpos)
00115   {
00116     GRAD_STACK1->write_grad_stack_buffer();
00117   }
00118 
00119   // check to see if anything has been written into the file
00120   off_t last_cpos=lseek(fp->file_ptr,0L,SEEK_CUR);
00121 
00122   //save current contents of the buffer so we can get them later
00123   if (last_cpos)
00124   {
00125     fp->write_cmpdif_stack_buffer();
00126   }
00127 
00128   for (int i=jac.rowmin();i<=jac.rowmax();i++)
00129   {
00130     if (jac(i).indexmin() !=1)
00131     {
00132       cerr  << "jacobian matrix minimum column index must equal 1"
00133                " for all rows" << endl;
00134       return;
00135     }
00136     if (jac(i).indexmax() < nvar)
00137     {
00138       cerr  << "jacobian matrix column size is less than the number of"
00139                " independent variables" << endl;
00140       return;
00141     }
00142   }
00143   // save variable values if desired
00144   if (save_var_flag)
00145   {
00146     save_arrays();
00147     save_variables();
00148   }
00149   // now evalueate the jacobian
00150   for (int ijac=1;ijac<=depvar_count;ijac++)
00151   {
00152     dvector& g=(dvector&)(jac(ijac));
00153     //max_num_dependent_variables=ndv;
00154     if (depvar_count>DEPVARS_INFO->max_num_dependent_variables)
00155     {
00156       cout << "maximum number of depdendent variables of "
00157          << DEPVARS_INFO->max_num_dependent_variables << " exceeded "
00158          << endl
00159          << "use gradient_structure::set_NUM_DEPENDENT_VARIABLES(int i);"
00160          << endl << "to increase the number of dependent variables"
00161          << endl;
00162       ad_exit(1);
00163     }
00164 
00165     fp->offset=DEPVARS_INFO->cmpdif_buffer_position(ijac);
00166     fp->toffset=fp->offset;
00167     _GRADFILE_PTR=DEPVARS_INFO->grad_file_count(ijac);
00168     fp->file_ptr=DEPVARS_INFO->cmpdif_file_count(ijac);
00169     lpos=DEPVARS_INFO->grad_file_position(ijac);
00170     // position the cmpdif file correctly;
00171     if (last_cpos)
00172     {
00173       off_t cmp_lpos=DEPVARS_INFO->cmpdif_file_position(ijac);
00174       lseek(fp->file_ptr,cmp_lpos,SEEK_SET);
00175       fp->read_cmpdif_stack_buffer(cmp_lpos);
00176     }
00177     GRAD_STACK1->_GRADFILE_PTR = GRAD_STACK1->gradfile_handle();
00178 
00179     if (last_gpos)
00180     {
00181       // just use the end of the buffer
00182       GRAD_STACK1->set_gbuffer_pointers();
00183 
00184       // check to sere if buffer was written into the beginning of
00185       // the next file
00186       if ( (GRAD_STACK1->_GRADFILE_PTR == GRAD_STACK1->_GRADFILE_PTR1)
00187          && (lpos == GRAD_STACK1->end_pos1) && (lpos>0) )
00188       {
00189         // get the next file
00190         GRAD_STACK1->increment_current_gradfile_ptr();
00191         lpos=0;
00192       }
00193       // and position the file to the begginig of the buffer image
00194       lseek(_GRADFILE_PTR,lpos,SEEK_SET);
00195       // now fill the buffer with the appropriate stuff
00196       GRAD_STACK1->read_grad_stack_buffer(lpos);
00197       // now reposition the grad_buffer pointer
00198     }
00199     GRAD_STACK1->ptr=
00200          (grad_stack_entry *)DEPVARS_INFO->grad_buffer_position(ijac);
00201 
00202     if(GRAD_STACK1->ptr <= GRAD_STACK1->ptr_first)
00203     {
00204 #ifdef SAFE_ALL
00205         cerr << "warning -- calling gradcalc when no calculations generating"
00206          << endl << "derivative information have occurred" << endl;
00207 #endif
00208       g.initialize();
00209       return;
00210     }    // current is one past the end so -- it
00211 
00212     gradient_structure::GRAD_STACK1->ptr--;
00213 
00214     for (unsigned int i=0; i< GRAD_LIST->nlinks; i++)
00215     {
00216       * (double*) (GRAD_LIST->dlink_addresses[i]) = 0;
00217     }
00218 
00219     double_and_int* tmp =
00220       (double_and_int*)gradient_structure::ARRAY_MEMBLOCK_BASE;
00221 
00222     unsigned long int max_last_offset
00223                = gradient_structure::ARR_LIST1->get_max_last_offset();
00224 
00225     unsigned int size = sizeof(double_and_int);
00226 
00227     for (unsigned int i = 0; i < (max_last_offset/size); i++)
00228     {
00229       tmp->x = 0;
00230 #if defined (__ZTC__)
00231       tmp = (double_and_int*)_farptr_norm((void*)(++tmp));
00232 #else
00233       tmp++;
00234 #endif
00235     }
00236 
00237     *gradient_structure::GRAD_STACK1->ptr->dep_addr  = 1;
00238 
00239     int break_flag=1;
00240     do
00241     {
00242       gradient_structure::GRAD_STACK1->ptr++;
00243       #ifdef FAST_ASSEMBLER
00244         gradloop();
00245       #else
00246         //int counter=0;
00247       while (gradient_structure::GRAD_STACK1->ptr-- >
00248              gradient_structure::GRAD_STACK1->ptr_first)
00249       {
00250         //grad_stack_entry* grad_ptr =
00251         //gradient_structure::GRAD_STACK1->ptr;
00252         {
00253           (* gradient_structure::GRAD_STACK1->ptr->func)();
00254         }
00255       }
00256       #endif
00257       // \todo Need test
00258       // back up the file one buffer size and read forward
00259       off_t offset = (off_t)(sizeof(grad_stack_entry)
00260         * gradient_structure::GRAD_STACK1->length);
00261       lpos = lseek(gradient_structure::GRAD_STACK1->_GRADFILE_PTR,
00262         -offset, 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 (int i=0; i<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 }