ADMB Documentation  11.1.2274
 All Classes Files Functions Variables Typedefs Friends Defines
jacobclc.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: jacobclc.cpp 1946 2014-04-29 05:22: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 __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   long int 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   long int 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   long int 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       long int cmp_lpos=DEPVARS_INFO->cmpdif_file_position(ijac);
00175       lseek(fp->file_ptr,cmp_lpos,SEEK_SET);
00176       #ifndef __MSC__
00177         fp->read_cmpdif_stack_buffer(cmp_lpos);
00178       #else
00179         fp->read_cmpdif_stack_buffer((long int&)cmp_lpos);
00180       #endif
00181     }
00182     GRAD_STACK1->_GRADFILE_PTR = GRAD_STACK1->gradfile_handle();
00183 
00184     if (last_gpos)
00185     {
00186       // just use the end of the buffer
00187       GRAD_STACK1->set_gbuffer_pointers();
00188 
00189       // check to sere if buffer was written into the beginning of
00190       // the next file
00191       if ( (GRAD_STACK1->_GRADFILE_PTR == GRAD_STACK1->_GRADFILE_PTR1)
00192          && (lpos == GRAD_STACK1->end_pos1) && (lpos>0) )
00193       {
00194         // get the next file
00195         GRAD_STACK1->increment_current_gradfile_ptr();
00196         lpos=0;
00197       }
00198       // and position the file to the begginig of the buffer image
00199       lseek(_GRADFILE_PTR,lpos,SEEK_SET);
00200       // now fill the buffer with the appropriate stuff
00201       GRAD_STACK1->read_grad_stack_buffer(lpos);
00202       // now reposition the grad_buffer pointer
00203     }
00204     GRAD_STACK1->ptr=
00205          (grad_stack_entry *)DEPVARS_INFO->grad_buffer_position(ijac);
00206 
00207     if(GRAD_STACK1->ptr <= GRAD_STACK1->ptr_first)
00208     {
00209 #ifdef SAFE_ALL
00210         cerr << "warning -- calling gradcalc when no calculations generating"
00211          << endl << "derivative information have occurred" << endl;
00212 #endif
00213       g.initialize();
00214       return;
00215     }    // current is one past the end so -- it
00216 
00217     gradient_structure::GRAD_STACK1->ptr--;
00218 
00219     for (i=0; i< GRAD_LIST->nlinks; i++)
00220     {
00221       * (double*) (GRAD_LIST->dlink_addresses[i]) = 0;
00222     }
00223 
00224     double_and_int* tmp =
00225       (double_and_int*)gradient_structure::ARRAY_MEMBLOCK_BASE;
00226 
00227     unsigned long int max_last_offset
00228                = gradient_structure::ARR_LIST1->get_max_last_offset();
00229 
00230     unsigned int size = sizeof(double_and_int);
00231 
00232     for (i = 0; i < (max_last_offset/size); i++)
00233     {
00234       tmp->x = 0;
00235 #if defined (__ZTC__)
00236       tmp = (double_and_int*)_farptr_norm((void*)(++tmp));
00237 #else
00238       tmp++;
00239 #endif
00240     }
00241 
00242     *gradient_structure::GRAD_STACK1->ptr->dep_addr  = 1;
00243 
00244     int break_flag=1;
00245     do
00246     {
00247       gradient_structure::GRAD_STACK1->ptr++;
00248       #ifdef FAST_ASSEMBLER
00249         gradloop();
00250       #else
00251         //int counter=0;
00252       while (gradient_structure::GRAD_STACK1->ptr-- >
00253              gradient_structure::GRAD_STACK1->ptr_first)
00254       {
00255         //grad_stack_entry* grad_ptr =
00256         //gradient_structure::GRAD_STACK1->ptr;
00257         {
00258           (* gradient_structure::GRAD_STACK1->ptr->func)();
00259         }
00260       }
00261       #endif
00262 
00263   // back up the file one buffer size and read forward
00264       lpos = lseek(gradient_structure::GRAD_STACK1->_GRADFILE_PTR,
00265         -((long int)(sizeof(grad_stack_entry)*gradient_structure::
00266         GRAD_STACK1->length)),SEEK_CUR);
00267 
00268        break_flag=gradient_structure::
00269                   GRAD_STACK1->read_grad_stack_buffer(lpos);
00270     }  while (break_flag); // do
00271 
00272     int mindx = g.indexmin();
00273     dvector & gg=(dvector&)(g);
00274     for (i=0; i<(unsigned int)nvar; i++)
00275     {
00276       gg[i+mindx] =  * gradient_structure::INDVAR_LIST->get_address(i);
00277       //g[i+mindx] =  * gradient_structure::INDVAR_LIST->get_address(i);
00278     }
00279     gradient_structure::GRAD_STACK1->ptr =
00280          gradient_structure::GRAD_STACK1->ptr_first;
00281   }// loop over dep vars
00282   DEPVARS_INFO->depvar_count=0;
00283   if (gradient_structure::save_var_flag)
00284   {
00285     gradient_structure::restore_arrays();
00286     gradient_structure::restore_variables();
00287   }
00288 }