fvar.hpp

Go to the documentation of this file.

#if defined(__GNUDOS__)
#  pragma interface
#endif
//#define USE_DDOUBLE

#ifndef FVAR_HPP
#define FVAR_HPP
#define USE_VECTOR_SHAPE_POOL
//#define MYDEBUG
#if defined(USE_DDOUBLE)
#  include <qd/qd.h>
//#  include </usr/local/include/qd/qd.h>
# define  double dd_real
  dd_real operator ++ (dd_real& d) { return d+=1.0;}
  dd_real pow(const dd_real& d,int i);
  dd_real pow(const dd_real& d,const dd_real&);
#endif
#if defined(__TURBOC__) && defined(__linux__)
#define __STL_PTHREADS
#  include <k3stl.h>
#  include <pthread_alloc>
#endif
#define  __NUMBERVECTOR__
#define ADUNCONST(type,obj) type & obj = (type&) _##obj;

#define my_off_t long int
#define my_u_off_t unsigned long int

#define  MFCL2_CONSTRUCTORS

#if (defined(__GNUDOS__))
#  if (  (__GNUC__>3) || (__GNUC__==3 && __GNUC_MINOR__>=0)  )
#    define __GNU_NEWER__
#    define GCC3
#  endif
#endif

#if (defined(__MSVC32__))
#  if (__MSVC32__>=7)
#    define __MSC_NEWER__
#  endif
#endif

#if defined(__MINGW32__)
#  include <conio.h>
#endif
//#define USE_ADPVM 1
//#define AD_DEMO
#if defined(_ADEXEP)
#  define USE_EXECPTIONS
#endif
//#define CHK_ID_STRING
//#define SAFE_ALL
#define DOS386
//#define AD_FAST_ASSIGN
// file fvar.hpp
// try to just use BOR_CONST = const
#define __USE_IOSTREAM__
#if defined(__BORLANDC__)
#  if (__BORLANDC__  >= 0x0550) 
#    include <fcntl.h>
#  endif
#  if (__BORLANDC__  >= 0x0520) 
#    define BOR_CONST const
// #    define __USING_STD_NAMES__
#  else
#  define BOR_CONST
#endif
#else
#  define BOR_CONST
#endif
#define USE_HIGHER_ARRAYS
#if defined(__BORLANDC__) || defined (__MSVC32__) || defined(__WAT32__)
#  include <io.h>
#endif

// test
//#if defined (__MSVC32__) 
#  if defined (BOR_CONST) 
#    undef BOR_CONST
#  endif
#  define BOR_CONST const
//#endif

#if defined(linux) || defined (__GNUDOS__)
#define AD_LONG_INT long int
//#define AD_LONG_INT long long int
#else
#define AD_LONG_INT long int
#endif

#if defined(linux) || defined (__GNUDOS__)
#  include <unistd.h>
#endif

#if defined(linux) || defined(__ADSGI__) || defined(__linux__)
#include <fcntl.h> // to get fstreambase
#  if !defined(O_BINARY)
#    define O_BINARY 0
#  endif
  #define UNIXKLUDGE
#endif

#if !defined(linux) 
#  include <dos.h>
#  if !defined(__MSVC32__) && !defined(__BORLANDC__) && !defined (__GNUDOS__)
#    include <pc.h>
#  endif
#endif
#define AD_VERSION 7

#if defined(__GNUDOS__)
extern "C" {
      void clrscr(void);
}
  #define CGNU_DOUBLE double 
#else
#  if defined(__BORLANDC__)
#    if (__BORLANDC__  > 0x0520) 
#      define CGNU_DOUBLE double 
#    else
#      define CGNU_DOUBLE  double 
#    endif
#  else
#    define CGNU_DOUBLE  double 
#  endif
#endif

#ifdef __GNUDOS__
  #define DOS386 
  #define __GNU__
  #ifdef OLDGNU
  #include <gccmanip.h>
  #endif
//#define HUGE_VAL  1.7976931348623158E+308
#endif

// for compilers that have trouble with const
#define USE_CONST
#ifdef USE_CONST
  #define _CONST const
#else
  #define _CONST 
#endif
#ifdef __GNUDOS__
 // #include <admanip.h>
#endif

#ifdef _SNGL_ 
  #define double float
#endif

#ifndef NO_DERIVS
  #define  NO_DERIVS
#endif

#ifdef SAFE_ALL
//#define INCLUDE_BLOCKSIZE       // add size info to df1b2variable
  #define CHK_ID_STRING
  #define DSAFE_ALLOCATE
//  #define SAFE_INITIALIZE
  #define SAFE_ARRAYS
#endif 

#ifdef __NDP486__
#include "short.def"
#endif

#ifdef USE_ASSEMBLER
  #include "adass.h"
#endif

extern "C" {
typedef int (* fptr) (const char * format, ...) ;
extern fptr ad_printf;
typedef void (*exitptr)(int);
extern exitptr ad_exit;

void spdll_exit(int);
}

  class smart_counter
  {
    int * ncopies;
  public:
    int * get_ncopies(void);
    smart_counter(void);
    smart_counter(const smart_counter & sc);
    ~smart_counter(void);
  };

  class double_and_int
  {
  public:
    double x;
    inline double& xvalue(void){return x;}
  };

class banded_symmetric_dvar_matrix;
class banded_symmetric_dmatrix;
class banded_lower_triangular_dmatrix;
class banded_lower_triangular_dvar_matrix;
class random_number_generator;
class ad_integer;
class index_type;
class dlink;
class dvector;
class dmatrix;
class tdmatrix;
class imatrix;
class prevariable;
class dvariable;
class dvar_vector;
class ivector;
class lvector;
class dvar_matrix;
class uostream;
class uistream;
class arr_link;
class arr_list;
class d3_array;
class d4_array;
class d5_array;
class d6_array;
class d7_array;
class dvar3_array;
class dvar4_array;
class dvar5_array;
class dvar6_array;
class dvar7_array;
class dlist;
class indvar_offset_list;
class dvar_matrix_position;
class DF_FILE;
class dmatrix_position;
class d3_array_position;
class dvector_position;
class vector_shape;
class ivector_position;
class kkludge_object;
class dvar_vector_position;
dvector restore_dvar_vector_derivatives(void);
class gradient_structure;
class dependent_variables_information;
#if defined(__GNUC__) 
#  if (__GNUC__  >= 3) 
#    include <fstream>
#  else
#    include <fstream.h>
#  endif
#elif defined(__MSVC32__) 
#  if (__MSVC32__  >= 7)
#    include <fstream>
#  else
#    include <fstream.h>
#  endif
#else
#  include <fstream.h>
#endif
#include <stdio.h>
//#define __SUN__ 1
# if defined (__WAT32__)
#  if defined (max)
#    undef max
#  endif
#  if defined (min)
#    undef min
#  endif
#endif
#if defined (__TURBOC__) || defined (__WAT32__) 
#include <mem.h>  // to get memcpy and NULL
#include <fstream.h> // to get fstreambase
#include <iostream.h>
#include <dos.h>
#undef __ZTC__
#define __USE_IOSTREAM__
#undef __SUN__
#endif

#if ( defined(__GNUC__) && __GNUC__<3)
#if defined(linux) || defined(__CYGWIN__)
#   include <strstream.h>
#  else
#   include <strstrea.h>
#  endif
#endif

#if defined(__BORLANDC__)
#  if (__BORLANDC__  < 0x0550) 
#    include <fcntl.h> // to get fstreambase
#    include <strstrea.h>
#  else
#    include <strstream>
#  endif
#endif

#if defined(__MSVC32__ ) 
#  if (__MSVC32__  >= 7)
#    include <iostream>
#  else
#    include <iostream.h>
#    include <strstrea.h>
#  endif
#include <stddef.h>
#include <fcntl.h> // to get fstreambase
#include <dos.h>
#undef __ZTC__
#define __USE_IOSTREAM__
#undef __SUN__
#endif

#ifdef __ZTC__
#include <fstream.hpp> // to get fstreambase
#include <iostream.hpp>
#define __USE_IOSTREAM__
#include <dos.h>
#undef __SUN__
#endif

#ifdef __NDPX__
#include <fstream.h>
#include <iostream.h>
#define __USE_IOSTREAM__
#undef __SUN__
#endif

#ifdef __SUN__
#undef __NDPX__
#include <fstream.h>
#include <iostream.h>
#define __USE_IOSTREAM__

#ifndef _FPOS_T_DEFINED
#ifdef __GNUDOS__
#  if defined(linux) || defined(__ADSGI__)|| defined(__CYGWIN32__)
      typedef long int fpos_t;
#  else
      typedef unsigned long int fpos_t;
#  endif
# undef __SUN__
#elif defined (__WAT32__)
  typedef long fpos_t;
#else
  typedef int fpos_t;
#endif
  #define _FPOS_T_DEFINED
#endif

#if defined(__GNUC__)
extern "C" int getch(void);
#else
int getch(void);
#endif

#endif
#include <math.h>

#if defined(__BORLANDC__)
#  if (__BORLANDC__  > 0x0520) 
     using std::ofstream;
     using std::ifstream;
#  endif
#endif

#if ( (defined(__GNUC__) && __GNUC__>=3)   \
    || (defined(__MSVC32__) && __MSVC32__>=7)  \
    || (defined(__BORLANDC__) && __BORLANDC__>=0x0570  && defined(__linux__)) )
#  include <iostream>
#  include <iomanip>
#  include <sstream>
#  include <istream>
#  include <strstream>
     using std::ofstream;
     using std::ostream;
     using std::ifstream;
     using std::istream;
     using std::istrstream;
     using std::streampos;
     using std::streambuf;
     using std::setw;
     using std::setprecision;
     using std::istringstream;
     using std::ios;
     using std::cerr;
     using std::cin;
     using std::cout;
     using std::endl;
#else
#  include <iomanip.h>     
#endif     

#define BEGIN_MINIMIZATION(nvar, objective_function, ind_vars, gradient, cntrl) \
gradient_structure gs; \
while (cntrl.ireturn >= 0) \
{ \
  cntrl.fmin(objective_function,ind_vars,gradient ); \
  if (cntrl.ireturn > 0) \
  {

#define END_MINIMIZATION(nvar,gradient) \
    gradcalc(nvar, gradient); \
  } \
}

void default_evaluation(void);
void default_evaluation0(void);
void default_evaluation1(void);
void default_evaluation1m(void);
void default_evaluation2(void);
void default_evaluation3(void);
void default_evaluation4(void);
void default_evaluation4m(void);

void myheapcheck( char *);

void RETURN_ARRAYS_INCREMENT(void);
void RETURN_ARRAYS_DECREMENT(void);

void * farptr_norm(void *);
long int farptr_tolong(void *) ;
#if (!defined(__ZTC__) || defined(DOS386))
  long int _farptr_tolong(void *);
#endif
void denormalize_ptr(void * , unsigned int byte_offset);

class i3_array;

ostream& operator<<(BOR_CONST ostream& ostr,_CONST dmatrix& z);
istream& operator>>(BOR_CONST istream& istr,BOR_CONST dmatrix& z);
uostream& operator<<(BOR_CONST uostream& ostr,_CONST dmatrix& z);
uistream& operator>>(BOR_CONST uistream& istr,BOR_CONST dmatrix& z);

ostream& operator<<(BOR_CONST ostream& ostr,_CONST d3_array& z);
istream& operator>>(BOR_CONST istream& istr,BOR_CONST d3_array& z);
uostream& operator<<(BOR_CONST uostream& ostr,_CONST d3_array& z);
uistream& operator>>(BOR_CONST uistream& istr,BOR_CONST d3_array& z);

ostream& operator<<(BOR_CONST ostream& ostr,_CONST dvar3_array& z);
istream& operator>>(BOR_CONST istream& istr,BOR_CONST dvar3_array& z);
uostream& operator<<(BOR_CONST uostream& ostr,_CONST dvar3_array& z);
uistream& operator>>(BOR_CONST uistream& istr,BOR_CONST dvar3_array& z);

ostream& operator<<(BOR_CONST ostream& ostr,_CONST ivector& z);
istream& operator>>(BOR_CONST istream& istr,BOR_CONST ivector& z);
uostream& operator<<(BOR_CONST uostream& ostr,_CONST ivector& z);
uistream& operator>>(BOR_CONST uistream& istr,BOR_CONST ivector& z);

ostream& operator<<(BOR_CONST ostream& ostr,_CONST lvector& z);
istream& operator>>(BOR_CONST istream& istr,BOR_CONST lvector& z);
uostream& operator<<(BOR_CONST uostream& ostr,_CONST lvector& z);
uistream& operator>>(BOR_CONST uistream& istr,BOR_CONST lvector& z);

ostream& operator<<(BOR_CONST ostream& ostr,const dvector& z);
istream& operator>>(BOR_CONST istream& istr,BOR_CONST dvector& z);
uostream& operator<<(BOR_CONST uostream& ostr,_CONST dvector& z);
uistream& operator>>(BOR_CONST uistream& istr,BOR_CONST dvector& z);

ostream& operator<<(BOR_CONST ostream& ostr,_CONST dvar_vector& z);
istream& operator>>(BOR_CONST istream& istr,BOR_CONST dvar_vector& z);
uostream& operator<<(BOR_CONST uostream& ostr,_CONST dvar_vector& z);
uistream& operator>>(BOR_CONST uistream& istr,BOR_CONST dvar_vector& z);

ostream& operator<<(BOR_CONST ostream& ostr,_CONST dvar_matrix& z);
istream& operator>>(BOR_CONST istream& istr,BOR_CONST dvar_matrix& z);
uostream& operator<<(BOR_CONST uostream& ostr,_CONST dvar_matrix& z);
uistream& operator>>(BOR_CONST uistream& istr,BOR_CONST dvar_matrix& z);


ostream& operator<<(BOR_CONST ostream& ostr,_CONST prevariable& z);
istream& operator>>(BOR_CONST istream& istr,BOR_CONST prevariable& z);
uostream& operator<<(BOR_CONST uostream& ostr,_CONST prevariable& z);
uistream& operator>>(BOR_CONST uistream& istr,BOR_CONST prevariable& z);

ostream& setscientific(BOR_CONST ostream& s);

//ostream& setshowpoint(BOR_CONST ostream& s);
class preshowpoint
{
  int x;
};
preshowpoint  setshowpoint(void);
ostream& operator << (BOR_CONST ostream&,preshowpoint);

#if (__MSVC32__>=7 || __BORLANDC__  >= 0x0560) 
#   define setfixed() std::fixed
#else
  ostream& setfixed(BOR_CONST ostream& s);
  class prefixed
  {
    int x;
  };
  prefixed  setfixed(void);
  ostream& operator << (BOR_CONST ostream&,prefixed);
#endif
  
#if (__MSVC32__>=7 || __BORLANDC__  >= 0x0560) 
#   define setscientific() std::scientific
#else
  class prescientific
  {
    int x;
  };
  prescientific  setscientific(void);
  ostream& operator << (BOR_CONST ostream&,prescientific);
#endif

istream& operator>>(BOR_CONST istream& istr,BOR_CONST imatrix& z);
ostream& operator<<(BOR_CONST ostream& istr,_CONST imatrix& z);
istream& operator>>(BOR_CONST istream& istr,BOR_CONST i3_array& z);
ostream& operator<<(BOR_CONST ostream& istr,_CONST i3_array& z);
  class grad_stack;

  class kkludge_object
  {
    char u;
  };

      class humungous_pointer
      {
      public:
        #if defined(__BORLANDC__ ) && !defined(DOS386) && !defined(__WIN32__)
    char huge * ptr;
  #else
    char * ptr;
  #endif
          int adjustment;
        void free(void);
  ~humungous_pointer(){;}
        void adjust(int);
  humungous_pointer();
  humungous_pointer operator + (unsigned long int& offset);
  humungous_pointer& operator += (unsigned long int& offset);
  humungous_pointer& operator = (void * p);
  int operator == (void * p);
  int operator != (void * p); 
        operator char* (); 
        operator void* (); 
        operator double_and_int* ();
        operator double* ();
      };

#ifndef _VECTOR_SHAPE
  #define _VECTOR_SHAPE
#include <dfpool.h>
  class vector_shape_pool : public dfpool
  {
  public:
    vector_shape_pool(void);
    vector_shape_pool(int);
  };

  class vector_shape
  {
  public:
#if defined(USE_VECTOR_SHAPE_POOL)
    static vector_shape_pool * xpool;
    void * operator new(size_t);
    void operator delete(void * ptr,size_t n)
    {  xpool->free(ptr); }
#endif
    unsigned int ncopies;
    void shift(int min);
    int index_min;
    int index_max;
  private:
    friend class dvector;
    //friend class tdvector;
    friend class subdvector;
    friend class dvar_vector;
    friend class ivector;
    friend class lvector;
    friend class ptr_vector;
  public:
    int decr_ncopies(void) {return --ncopies;}
    int get_ncopies(void) {return ncopies;}
    int incr_ncopies(void) {return ++ncopies;}
    vector_shape(int lb, int lu) {index_min=lb;index_max=lu;ncopies=0;}
    int indexmin(){return index_min;}
    int indexmax(){return index_max;}
  };
#endif

  class ptr_vector
  {
    void ** v;
    vector_shape * shape;

  public:

    int indexmin() _CONST {return shape->index_min;} // returns the minimum allowable index
    int indexmax() _CONST  {return shape->index_max;} // returns the maximum allowable index
    int size() _CONST {return shape->index_max - shape->index_min+1;} // returns the maximum allowable index

    void shift(int min);
    ptr_vector( int ncl,  int ncu);
    ptr_vector();
    void allocate(int,int);
    //operator void ** ();
   ~ptr_vector();
    void*& operator[] (int i);
    void*& operator() (int i);
    //void*& elem(int i);
    void*& elem(int i) { return(*(v+i)); }
    int operator ! (void) _CONST {return (shape == NULL);}
    int operator () (void) _CONST {return (shape != NULL);}

    ptr_vector(BOR_CONST ptr_vector& t);
    ptr_vector& operator = (BOR_CONST ptr_vector& t);
    void initialize();
 }; // end of class ptr_vector

  class preivector
  {
    ivector * p;
    int lb;
    int ub;
    inline preivector(ivector * _p,int _lb,int _ub) {p=_p;lb=_lb,ub=_ub;}
    friend class ivector;
  };

  class vector_shapex;
  class ivector
  {
  protected:
    int index_min;
    int index_max;
    int * v;
    vector_shapex * shape;

  public:
    int allocated(void) _CONST {return (shape != NULL);}

    ivector& operator -- (void) 
    { 
      index_min--;index_max--;v++; return *this;
    }
    void reallocate(double size);

    ivector& operator ++ (void)  
    { 
      index_min++;index_max++;v--; return *this;
    }
    friend class vector_index;
    int operator ! (void) _CONST {return (shape == NULL);}

        int& elem(int i) { return(v[i]); }
#ifdef USE_CONST
        _CONST int& elem(int i) _CONST { return(v[i]); }
#endif
    int indexmin()_CONST  {return index_min;} // returns the minimum allowable index
    int indexmax()_CONST  {return index_max;} // returns the maximum allowable index
    int size() _CONST {return index_max - index_min+1;} // returns the maximum allowable index
    int * get_v(void) _CONST {return v;}
    void shift(int min);
    //ivector(unsigned int sz); //makes an array [0..sz]
    //inline allocated(void) {return ::allocated(this);}
    ivector(_CONST lvector&);
    ivector(void);

    ivector(_CONST dvector&);

    void fill(const char * s);
    void fill_seqadd(int, int);
    void fill_multinomial(BOR_CONST int& seed,_CONST dvector& p);
    void fill_multinomial(BOR_CONST random_number_generator& rng,_CONST dvector& p);
    ivector( int ncl,  int ncu);
    ivector(const ad_integer& ncl,const index_type& ncu);
    void allocate(const ad_integer& ncl,const index_type& ncu);
             // makes an array [ncl..ncu]

    ivector( unsigned int sz, long int * x);

//#ifdef HOME_VERSION
      void safe_allocate(int ncl,int ncu);
      void allocate(int ncl,int ncu);
      void allocate(_CONST ivector& iv);
      void allocate(_CONST dvector& dv);
      void allocate(const char *);
      void allocate(void);
      void deallocate(void);
      void safe_deallocate(void);
//#endif

    //operator int* () {return v;}

   ~ivector();

    void save_ivector_position(void) _CONST;

    void save_ivector_value(void) _CONST;

    void write_on(BOR_CONST ostream& s) _CONST;

    void read_from(BOR_CONST istream& s);

    void write_on(BOR_CONST uostream& s) _CONST;

    void read_from(BOR_CONST uistream& s);

 #ifdef OPT_LIB
    inline int& operator[] (int i) {return(*(v+i));}
    inline int& operator() (int i) {return(*(v+i));}
 #else
    int& operator[] (int i);
    int& operator() (int i);
 #endif

#ifdef USE_CONST
 #ifdef OPT_LIB
    inline _CONST int& operator[] (int i) _CONST {return(*(v+i));}
    inline _CONST int& operator() (int i) _CONST {return(*(v+i));}
 #else
    _CONST int& operator[] (int i) _CONST;
    _CONST int& operator() (int i) _CONST;
 #endif
#endif

    ivector sub(int lb,int ub)
    {
      return preivector(this,lb,ub);
    }
    ivector operator () (int lb,int ub)
    {
      return preivector(this,lb,ub);
    }
    ivector(const preivector& pd);

    ivector operator ()(_CONST ivector& u);

    ivector(_CONST ivector& t);
    void shallow_copy(_CONST ivector& t);

    ivector& operator = (_CONST ivector& t);
    ivector& operator = (int t);

    ivector& operator = (_CONST lvector& t);

    friend class imatrix;
    friend class i3_array;
    void initialize(void);
    ivector& operator += (_CONST ivector& t);
    ivector& operator += (int t);

 }; // end of class ivector

  ivector operator + (_CONST ivector& v,_CONST ivector& w);
  ivector operator - (_CONST ivector& v,_CONST ivector& w);
  ivector operator + (_CONST ivector& v,int w);
  ivector operator - (_CONST ivector& v,int w);
  ivector operator + (int v,_CONST ivector& w);
  ivector operator - (int v,_CONST ivector& w);

class lvector_position
{
  int min;
  int max;
  int * v;
public:
  lvector_position(void);
  lvector_position(_CONST lvector& v);
  lvector_position(BOR_CONST lvector_position& dvp);
};

  class lvector
  {
     AD_LONG_INT * v;
    vector_shape * shape;

  public:
    int operator ! (void) _CONST {return (shape == NULL);}

        inline AD_LONG_INT& elem(int i) { return(v[i]); }
#ifdef USE_CONST
        inline _CONST AD_LONG_INT& elem(int i) _CONST { return(v[i]); }
#endif
    int indexmin() _CONST  {return shape->index_min;} // returns the minimum allowable index
    int indexmax() _CONST {return shape->index_max;} // returns the maximum allowable index
    int size() _CONST {return shape->index_max - shape->index_min+1;} // returns the maximum allowable index
    void shift(int min);

    void fill(const char * s);
    void fill_seqadd( BOR_CONST AD_LONG_INT&,  BOR_CONST AD_LONG_INT&);
    void fill_multinomial(BOR_CONST int& seed,_CONST dvector& p);
    void fill_multinomial(BOR_CONST random_number_generator& rng,_CONST dvector& p);

    //lvector(unsigned int sz); //makes an array [0..sz]

    lvector(_CONST lvector&);

    lvector(_CONST dvector&);

    lvector(_CONST ivector&);

    lvector(void);

    lvector( int ncl,  int ncu);
    void allocate( int ncl,  int ncu);
    void allocate(_CONST lvector&);
    void allocate(void);
                                           // makes an array [ncl..ncu]

    lvector( unsigned int sz,  AD_LONG_INT * x);

    operator  AD_LONG_INT* ();

   ~lvector();

    void write_on(BOR_CONST ostream& s) _CONST;

    void read_from(BOR_CONST istream& s);

    void write_on(BOR_CONST uostream& s) _CONST;

    void read_from(BOR_CONST uistream& s);

#ifdef OPT_LIB
    inline  AD_LONG_INT& operator[] (int i) {return v[i];}
    inline  AD_LONG_INT& operator() (int i) {return v[i];}
#else
     AD_LONG_INT& operator[] (int i);
     AD_LONG_INT& operator() (int i);
#endif

#ifdef USE_CONST
#ifdef OPT_LIB
    inline _CONST  AD_LONG_INT& operator[] (int i) _CONST {return v[i];}
    inline _CONST  AD_LONG_INT& operator() (int i) _CONST {return v[i];}
#else
    _CONST  AD_LONG_INT& operator[] (int i) _CONST ;
    _CONST  AD_LONG_INT& operator() (int i) _CONST ;
#endif
#endif

    lvector operator ()(_CONST lvector& u);

    lvector& operator = (_CONST lvector& t);

    void initialize(void);
    friend class lmatrix;
 }; // end of class lvector

  AD_LONG_INT sum(_CONST lvector&);

  
  class dependent_variables_information
  {
    int max_num_dependent_variables;
    int depvar_count;
    ptr_vector grad_buffer_position;
    lvector cmpdif_buffer_position;
    lvector grad_file_position;
    lvector cmpdif_file_position;
    lvector grad_file_count;
    lvector cmpdif_file_count;
    dependent_variables_information(int ndv);
    friend class gradient_structure;
  };
  dvar_vector_position restore_dvar_vector_position(void);
  dvector restore_dvar_vector_value(BOR_CONST dvar_vector_position& tmp);
  void arr_free(double_and_int *);
  double_and_int * arr_new(unsigned int sz); 

  class gradient_structure
  {
      static char cmpdif_file_name[61];
      static DF_FILE * fp;
public:
      #if defined(NO_DERIVS)
        static int no_derivatives;
      #endif
private:
      static long int USE_FOR_HESSIAN;
      static long int NVAR;
      static int NUM_RETURN_ARRAYS;
      static dvariable** RETURN_ARRAYS;
      static int RETURN_ARRAYS_PTR;
      static dvariable ** RETURN_PTR_CONTAINER;
      static long int TOTAL_BYTES;
      static long int PREVIOUS_TOTAL_BYTES;
      static unsigned long ARRAY_MEMBLOCK_SIZE; //js
      static humungous_pointer ARRAY_MEMBLOCK_BASE;
      static humungous_pointer ARRAY_MEMBLOCK_BASEA;
      static humungous_pointer ARRAY_MEMBLOCK_SAVE;
public:    
      static double * get_ARRAY_MEMBLOCK_BASE() 
      { 
        return (double*) ARRAY_MEMBLOCK_BASE;
      }
private:
      static double * variables_save;
      static long int CMPDIF_BUFFER_SIZE;
      static long int GRADSTACK_BUFFER_SIZE;
      static unsigned int MAX_NVAR_OFFSET;
      static int save_var_file_flag;
      static int save_var_flag;

      static int MAX_DLINKS;
      static indvar_offset_list * INDVAR_LIST;
      static int NUM_DEPENDENT_VARIABLES;
      static dependent_variables_information * DEPVARS_INFO;

      // this needs to be a static member function so other static
      // member functions can call it
      static void check_set_error(const char* variable_name);

      static int instances;
      int   x;

    public:
      class arrmemblerr{}; // exception class
   
      static int Hybrid_bounded_flag;
      static double * hessian_ptr;
      static int get_USE_FOR_HESSIAN(void) {return USE_FOR_HESSIAN;}
      static void set_USE_FOR_HESSIAN(int i) {USE_FOR_HESSIAN=i;}
      friend class dfsdmat;
      gradient_structure(long int size=100000L);  // constructor
      ~gradient_structure(void); // destructor
      static void save_variables(void);
      static void restore_variables(void);
      static void save_arrays(void);
      static void restore_arrays(void);
      static long int totalbytes(void);
      friend dvector restore_dvar_vector_value(BOR_CONST dvar_vector_position& tmp);
      friend void cleanup_temporary_files();
      //friend dvector restore_dvar_vector_value(BOR_CONST dvar_vector_position&);
      friend dvector restore_dvar_vector_derivatives(void);
      friend dmatrix restore_dvar_matrix_derivatives(void);
      friend dmatrix restore_dvar_matrix_value(void);
      //friend dmatrix restore_derivatives(void);
      friend void gradfree(dlink * v);
      friend double_and_int * arr_new(unsigned int sz); //js
      friend void arr_free(double_and_int *);
      friend void RETURN_ARRAYS_DECREMENT(void);
      friend void RETURN_ARRAYS_INCREMENT(void);
      friend void make_indvar_list(_CONST dvar_vector& t);
      //friend void gradcalc( int , double *);
      friend void gradcalc(int nvar,BOR_CONST dvector& g);
      friend void slave_gradcalc(void);
      friend void funnel_gradcalc(void);
      friend void allocate_dvariable_space(void);
      friend void wide_funnel_gradcalc(void);
      friend dvar_vector_position restore_dvar_vector_position(void);
      static grad_stack  * GRAD_STACK1;
      friend double_and_int * gradnew();
      static dlist * GRAD_LIST;
      static int RETURN_ARRAYS_SIZE;
      //static int RETURN_INDEX;
      static dvariable * RETURN_PTR;
      static dvariable * MIN_RETURN;
      static dvariable * MAX_RETURN;
      static arr_list * ARR_LIST1;
      static arr_list * ARR_FREE_LIST1;
      //static void funnel_jacobcalc(void);
      static void jacobcalc(int nvar,BOR_CONST dmatrix& jac);
      static void jacobcalc(int nvar,BOR_CONST ofstream& jac);
      static void jacobcalc(int nvar,BOR_CONST uostream& jac);

      friend void default_evaluation(void);
      //access functions

      friend class DF_FILE;
      static DF_FILE * get_fp(void){return fp;}
      static void set_NUM_RETURN_ARRAYS(int i);
     #if defined(NO_DERIVS)
      static void set_NO_DERIVATIVES(void);
      static void set_YES_DERIVATIVES(void);
     #endif
      static void set_YES_SAVE_VARIABLES_VALUES();
      static void set_NO_SAVE_VARIABLES_VALUES();
      //static int _GRADFILE_PTR; // should be int gradfile_handle;
      //static int _GRADFILE_PTR1; // should be int gradfile_handle;
      //static int _GRADFILE_PTR2; // should be int gradfile_handle;
      //static int _VARSSAV_PTR; // should be int gradfile_handle;
      static void set_NUM_DEPENDENT_VARIABLES(int i);
      static void set_RETURN_ARRAYS_SIZE(int i);
      static void set_ARRAY_MEMBLOCK_SIZE(unsigned long i);
      static void set_CMPDIF_BUFFER_SIZE(long int i);
      static void set_GRADSTACK_BUFFER_SIZE(long int i);
      static void set_MAX_NVAR_OFFSET(unsigned int i);
      static void set_MAX_DLINKS(int i);
      static long int NUM_GRADSTACK_BYTES_WRITTEN(void);
      friend class dlist;
      friend class grad_stack;
      static void save_dependent_variable_position
        (_CONST prevariable& v1);
      static long int max_last_offset;
      friend class function_minimizer;
      friend void funnel_derivatives(void);
  };

  void jacobcalc(int nvar,BOR_CONST dmatrix& g);
  void jacobcalc(int nvar,BOR_CONST ofstream& ofs);
  void jacobcalc(int nvar,BOR_CONST uostream& ofs);

#if defined(__BORLANDC__ )
#  pragma interface
#endif
    class dvect_ptr_ptr
    {
      dvector ** m;
    };
    class dlink
    {
       double_and_int  di;
       dlink *       prev;
    public:  // comments
      dlink * previous();
      inline double_and_int * get_address()
      {
        return &di;
      } //access function

      //friend tempvar();
      //friend class prevariable;
      //friend class tempvar;
      friend class dlist;
      friend void gradcalc(int nvar,BOR_CONST dvector& g);
      friend void slave_gradcalc(void);
      friend void gradloop();
      friend double_and_int * gradnew();
      friend void allocate_dvariable_space(void);
    };

    class dlist
    {
      dlink * last;
      unsigned int last_offset;
      unsigned int nlinks;
      dlink ** dlink_addresses;
      friend double_and_int * gradnew();
      friend void df_check_derivative_values(void);
      friend void df_check_derivative_values_indexed(void);
      friend void df_check_derivative_values_indexed_break(void);

    public:
      dlist();  // constructor
      void check_list(void);  // check list integrity
     ~dlist();  // destructor
      dlink * create();     //create a new link 
      void free_remove(dlink * rem);
      dlink * append(dlink *);  // add a link
      dlink* last_remove();
      friend void funnel_gradcalc(void);
      friend void slave_gradcalc(void);
      friend void gradcalc(int nvar,BOR_CONST dvector& g);
      friend void gradloop();
      friend void gradient_structure::restore_variables();
      friend void gradient_structure::save_variables();
      friend void gradient_structure::jacobcalc(int nvar,BOR_CONST dmatrix& jac);
      friend void allocate_dvariable_space(void);
      //friend void gradient_structure::funnel_jacobcalc(void);
      friend void gradient_structure::jacobcalc(int nvar,BOR_CONST ofstream& jac);
      friend void gradient_structure::jacobcalc(int nvar,BOR_CONST uostream& jac);
      friend void funnel_derivatives(void);
    };

    class indvar_offset_list;
    class grad_stack_entry
    {
    public:
      void (* func) (void);
      double * dep_addr;
      double * ind_addr1;
      double  mult1;
      double * ind_addr2;
      double  mult2;
    public:
      friend void gradcalc(int nvar,BOR_CONST dvector& g);
      friend void slave_gradcalc(void);
      friend void gradloop();
      friend void default_evaluation(void);
      friend class grad_stack;
      friend void gradient_structure::jacobcalc(int nvar,BOR_CONST dmatrix& jac);
      //friend void gradient_structure::funnel_jacobcalc(void);
    };
    void default_evaluation3ind(void);

    class grad_stack
    {
      grad_stack_entry * true_ptr_first;
      grad_stack_entry * ptr_first;
      grad_stack_entry * ptr_last;
      long int length;
      long int true_length;
    public:
      grad_stack_entry * ptr;
    private:
      //lvector * table;
      // js
      int _GRADFILE_PTR; // should be int gradfile_handle;
      int _GRADFILE_PTR1; // should be int gradfile_handle;
      int _GRADFILE_PTR2; // should be int gradfile_handle;
      int _VARSSAV_PTR; // should be int gradfile_handle;
      char gradfile_name[61];
      char gradfile_name1[61];
      char gradfile_name2[61];
      char var_store_file_name[61];
      void create_gradfile();
      long end_pos;
      long end_pos1;
      long end_pos2;
      dmatrix * table;
    public:
      friend void gradcalc(int nvar,BOR_CONST dvector& g);
      friend void slave_gradcalc(void);
      friend void funnel_gradcalc(void);
      friend void default_evaluation(void);
      friend void default_evaluation3ind(void);
      friend void reset_gradient_stack(void);
      friend void default_evaluation4ind(void);
      friend void grad_chk(void);
      friend void gradloop();
      friend void cleanup_temporary_files();
      ostream& operator  <<(grad_stack);
     void print();
      grad_stack();
     ~grad_stack();
      void write_grad_stack_buffer(void);

       void set_gradient_stack(void (* func) (void),
   double * dep_addr,double * ind_addr1=NULL, double mult1=0,
   double * ind_addr2=NULL, double mult2=0);

       void set_gradient_stack(void (* func) (void),
   double * dep_addr,double * ind_addr1,double * ind_addr2);


      void set_gradient_stack0(void (* func)(void),
        double * dep_addr);

       void set_gradient_stack1(void (* func) (void),
   double * dep_addr,double * ind_addr1);

       void set_gradient_stack2(void (* func) (void),
   double * dep_addr,double * ind_addr1, double mult1);

       void set_gradient_stack4(void (* func) (void),
   double * dep_addr,double * ind_addr1, double * ind_addr2);

       void set_gradient_stack(void (* func)(void),
        double * dep_addr,double * ind_addr1, double mult1, double * ind_addr2,
        double mult2,double * ind_addr3, double mult3,double * ind_addr4,
        double mult4);
    
       void set_gradient_stack(void (* func)(void),
        double * dep_addr,double * ind_addr1, double mult1, double * ind_addr2,
        double mult2,double * ind_addr3, double mult3);
    
      int read_grad_stack_buffer(my_off_t& lpos);
      void set_gradient_stack(void (*ptr)(void));
      void set_gbuffer_pointers(void);
       //js
      void increment_current_gradfile_ptr(void);
      int decrement_current_gradfile_ptr(void);
      //void open_gradfile();
      //void close_gradfile();
      #ifndef __MSC__
  int& gradfile_handle();
      #else
  int gradfile_handle();
      #endif
      char* get_gradfile_name();
      friend class gradient_structure;
      //int get_ngradfiles();
    }; // end of grad_stack


#ifdef OPT_LIB
    
    inline void grad_stack::set_gradient_stack(void (* func)(void),
      double * dep_addr,double * ind_addr1, double mult1, double * ind_addr2,
      double mult2)
    {
      #ifdef NO_DERIVS
  if (!gradient_structure::no_derivatives)
  {
      #endif
#     if defined(MYDEBUG)
          int wrote_buffer=0;
#     endif
    if (ptr > ptr_last)
    {
       // current buffer is full -- write it to disk and reset pointer
       // and counter
       this->write_grad_stack_buffer();
#     if defined(MYDEBUG)
             wrote_buffer=1;
#     endif
    }
#     if defined(MYDEBUG)
          if (wrote_buffer==1)
          {
            cout << "WROTE BUFFER" << endl;
          }
#     endif
    ptr->func = func;
    ptr->dep_addr = dep_addr;
    ptr->ind_addr1 = ind_addr1;
    ptr->mult1=mult1;
    ptr->ind_addr2 = ind_addr2;
    ptr->mult2=mult2;
    ptr++;
      #ifdef NO_DERIVS
  }
      #endif
    }
#else
  //  void grad_stack::set_gradient_stack(void (* func)(void),
   //   double * dep_addr,double * ind_addr1, double mult1, double * ind_addr2,
    //  double mult2);
#endif

    inline void grad_stack::set_gradient_stack(void (* func)(void),
      double * dep_addr,double * ind_addr1, double mult1, double * ind_addr2,
      double mult2,double * ind_addr3, double mult3)
    {
      #ifdef NO_DERIVS
  if (!gradient_structure::no_derivatives)
  {
      #endif
    if (ptr > ptr_last)
    {
       // current buffer is full -- write it to disk and reset pointer
       // and counter
       this->write_grad_stack_buffer();
    }
    ptr->func = NULL;
    ptr->dep_addr = dep_addr;
    ptr->ind_addr1 = ind_addr1;
    ptr->mult1=mult1;
    ptr->ind_addr2 = ind_addr2;
    ptr->mult2=mult2;
    ptr++;
    if (ptr > ptr_last)
    {
       // current buffer is full -- write it to disk and reset pointer
       // and counter
       this->write_grad_stack_buffer();
    }
    ptr->func = func;
    ptr->ind_addr1 = ind_addr3;
    ptr->mult1=mult3;
    ptr++;
      #ifdef NO_DERIVS
  }
      #endif
    }

    inline void grad_stack::set_gradient_stack(void (* func)(void),
      double * dep_addr,double * ind_addr1, double mult1, double * ind_addr2,
      double mult2,double * ind_addr3, double mult3,double * ind_addr4,
      double mult4)
    {
      #ifdef NO_DERIVS
  if (!gradient_structure::no_derivatives)
  {
      #endif
    if (ptr > ptr_last)
    {
       // current buffer is full -- write it to disk and reset pointer
       // and counter
       this->write_grad_stack_buffer();
    }
    ptr->func = NULL;
    ptr->dep_addr = dep_addr;
    ptr->ind_addr1 = ind_addr1;
    ptr->mult1=mult1;
    ptr->ind_addr2 = ind_addr2;
    ptr->mult2=mult2;
    ptr++;
    if (ptr > ptr_last)
    {
       // current buffer is full -- write it to disk and reset pointer
       // and counter
       this->write_grad_stack_buffer();
    }
    ptr->func = func;
    ptr->ind_addr1 = ind_addr3;
    ptr->mult1=mult3;
    ptr->ind_addr2 = ind_addr4;
    ptr->mult2=mult4;
    ptr++;
      #ifdef NO_DERIVS
  }
      #endif
    }

    inline void grad_stack::set_gradient_stack(void (* func)(void),
      double * dep_addr,double * ind_addr1,double * ind_addr2)
    {
      #ifdef NO_DERIVS
  if (!gradient_structure::no_derivatives)
  {
      #endif
    if (ptr > ptr_last)
    {
       // current buffer is full -- write it to disk and reset pointer
       // and counter
       this->write_grad_stack_buffer();
    }
    ptr->func = func;
    ptr->dep_addr = dep_addr;
    ptr->ind_addr1 = ind_addr1;
    ptr->ind_addr2 = ind_addr2;
    ptr++;
      #ifdef NO_DERIVS
  }
      #endif
    }

    inline void grad_stack::set_gradient_stack2(void (* func)(void),
      double * dep_addr,double * ind_addr1, double mult1)
    {
      #ifdef NO_DERIVS
        if (!gradient_structure::no_derivatives)
        {
      #endif
          if (ptr > ptr_last)
          {
             // current buffer is full -- write it to disk and reset pointer
             // and counter
             this->write_grad_stack_buffer();
          }
          ptr->func = func;
          ptr->dep_addr = dep_addr;
          ptr->ind_addr1 = ind_addr1;
          ptr->mult1=mult1;
          ptr++;
      #ifdef NO_DERIVS
        }
      #endif
    }

    inline void grad_stack::set_gradient_stack4(void (* func)(void),
      double * dep_addr,double * ind_addr1, double * ind_addr2)
    {
      #ifdef NO_DERIVS
  if (!gradient_structure::no_derivatives)
  {
      #endif
    if (ptr > ptr_last)
    {
       // current buffer is full -- write it to disk and reset pointer
       // and counter
       this->write_grad_stack_buffer();
    }
    ptr->func = func;
    ptr->dep_addr = dep_addr;
    ptr->ind_addr1 = ind_addr1;
    ptr->ind_addr2 = ind_addr2;
    ptr++;
      #ifdef NO_DERIVS
  }
      #endif
    }



   /*
    inline void grad_stack::set_gradient_stack(void (* func)(void),
      double * dep_addr,double * ind_addr1, double mult1)
    {
      #ifdef NO_DERIVS
  if (!gradient_structure::no_derivatives)
  {
      #endif
    if (ptr > ptr_last)
    {
       // current buffer is full -- write it to disk and reset pointer
       // and counter
       this->write_grad_stack_buffer();
    }
    ptr->func = func;
    ptr->dep_addr = dep_addr;
    ptr->ind_addr1 = ind_addr1;
    ptr->mult1=mult1;
          ptr++;
      #ifdef NO_DERIVS
        }
      #endif
    }
    */
   /*
    void grad_stack::set_gradient_stack1(void (* func)(void),
      double * dep_addr,double * ind_addr1)
    {
      #ifdef NO_DERIVS
        if (!gradient_structure::no_derivatives)
        {
      #endif
          if (ptr > ptr_last)
          {
             // current buffer is full -- write it to disk and reset pointer
             // and counter
             this->write_grad_stack_buffer();
          }
          ptr->func = func;
          ptr->dep_addr = dep_addr;
          ptr->ind_addr1 = ind_addr1;
          ptr++;
      #ifdef NO_DERIVS
        }
      #endif
    }
  */

   /*
    inline void grad_stack::set_gradient_stack(void (* func)(void),
      double * dep_addr,double * ind_addr1)
    {
      #ifdef NO_DERIVS
        if (!gradient_structure::no_derivatives)
        {
      #endif
          if (ptr > ptr_last)
          {
             // current buffer is full -- write it to disk and reset pointer
             // and counter
             this->write_grad_stack_buffer();
          }
          ptr->func = func;
          ptr->dep_addr = dep_addr;
          ptr->ind_addr1 = ind_addr1;
          ptr++;
      #ifdef NO_DERIVS
        }
      #endif
    }
   */

   /*
    inline void grad_stack::set_gradient_stack0(void (* func)(void),
      double * dep_addr)
    {
      #ifdef NO_DERIVS
        if (!gradient_structure::no_derivatives)
        {
      #endif
          if (ptr > ptr_last)
          {
             // current buffer is full -- write it to disk and reset pointer
             // and counter
             this->write_grad_stack_buffer();
          }
          ptr->func = func;
          ptr->dep_addr = dep_addr;
          ptr++;
      #ifdef NO_DERIVS
        }
      #endif
    }
   */

   /*
    inline void grad_stack::set_gradient_stack(void (* func)(void),
      double * dep_addr)
    {
      #ifdef NO_DERIVS
        if (!gradient_structure::no_derivatives)
        {
      #endif
          if (ptr > ptr_last)
          {
             // current buffer is full -- write it to disk and reset pointer
             // and counter
             this->write_grad_stack_buffer();
          }
          ptr->func = func;
          ptr->dep_addr = dep_addr;
          ptr++;
      #ifdef NO_DERIVS
        }
      #endif
    }
*/
    inline void grad_stack::set_gradient_stack(void (*func)(void))
    {
      #ifdef NO_DERIVS
        if (!gradient_structure::no_derivatives)
        {
      #endif
          if (ptr > ptr_last)
          {
            // current buffer is full -- write it to disk and reset pointer
            // and counter
            this->write_grad_stack_buffer();
          }

          ptr->dep_addr = NULL;
          ptr->func = func;
          ptr->ind_addr2 = NULL ;// want to put a long int                                            //into the memory space of a double
          ptr->mult2=0;
          ptr++;
      #ifdef NO_DERIVS
        }
      #endif
    }

    class indvar_offset_list
    {
      int nvar; // The number of independent variables
      double **  address;
  
    public:
      friend class gradient_structure;
      inline double * get_address(BOR_CONST int& i) { return(address[i]);}
      void put_address(unsigned int& i, double * iaddress) 
      { 
        address[i] = iaddress;
        //  cerr << "In put_address i = " << i << "\n";
      }
    };

    void gradfree(dlink *);

  class prevariable_position;

  class prevariable
  {
  protected:
#ifndef __SUN__
    prevariable(void) { }
#endif
#ifndef __NDPX__
    prevariable(double_and_int * u) {v=u;}
#endif

  public:
    double_and_int * v;
    friend class dvar_vector_iterator;
    friend class dvar_vector;
    friend class dvar_matrix;
    friend class dvar3_array;
    //shclass sc;
    friend class indvar_offset_list;
    friend class gradient_structure;
    friend double_and_int * gradnew();
    friend void make_indvar_list(int ,dvariable *);

    friend class banded_symmetric_dvar_matrix;
    friend class banded_lower_triangular_dvar_matrix;
    friend double& value(_CONST prevariable& v1);

    friend double * address(_CONST prevariable& v1);

    //void gradfree(dlink * v)

    friend prevariable& operator *(_CONST prevariable& v1,_CONST prevariable& v2);

    friend prevariable& operator *( double v1,_CONST prevariable& v2);

    friend prevariable& operator *(_CONST prevariable& v1, double v2);

    friend prevariable& operator / (_CONST prevariable& t1,_CONST prevariable& t2);

    friend prevariable& operator / ( double t1,_CONST prevariable& t2);

    friend prevariable& operator / (_CONST prevariable& t1, double t2);

    friend prevariable& sin(_CONST prevariable& t1);

    friend prevariable& fabs(_CONST prevariable& t1);
    friend prevariable& sigmoid(_CONST prevariable& t1);

    friend prevariable& sfabs(_CONST prevariable& t1); //"smoothed absolute value function

    friend prevariable& sqrt(_CONST prevariable& t1);
    friend prevariable& sqr(_CONST prevariable& t1);

    friend prevariable& exp(_CONST prevariable& t1);

    friend prevariable& atan(_CONST prevariable& t1);

    friend prevariable& tan(_CONST prevariable& t1);
    friend prevariable& tanh(_CONST prevariable& t1);

    friend prevariable& atan2(_CONST prevariable& t1,_CONST prevariable& t2);
    friend prevariable& atan2(_CONST prevariable& t1, double t2);
    friend prevariable& atan2( double t1,_CONST prevariable& t2);

    friend prevariable& acos(_CONST prevariable& t1);

    friend prevariable& asin(_CONST prevariable& t1);

    friend prevariable& cos(_CONST prevariable& t1);
    friend prevariable& sinh(_CONST prevariable& t1);
    friend prevariable& cosh(_CONST prevariable& t1);

    friend prevariable& log(_CONST prevariable& t1);
    friend prevariable& log10(_CONST prevariable& t1);

    friend prevariable& ldexp(_CONST prevariable&, BOR_CONST int&);


  public:

   void save_prevariable_position(void) _CONST;
   prevariable_position restore_prevariable_position(void);
   void save_prevariable_value(void) _CONST;
   double restore_prevariable_value(void);
   double restore_prevariable_derivative(void);


   inline double * xadr() { return( &(v->x) ); }
   inline double& xval() { return( (v->x) ); }

   inline double_and_int * & get_v()
   {
     return v;
   }
   #ifdef USE_CONST
     inline double_and_int *  get_v() _CONST 
     {
       return v;
     }
   #endif

   prevariable& operator=(_CONST prevariable&);
   prevariable& operator=( double );
#if (__BORLANDC__  >= 0x0540) 
   prevariable& operator=(_CONST prevariable&)_CONST ;
   prevariable& operator=( double )_CONST ;
#endif

   int operator ==(_CONST prevariable& v1) _CONST;
   int operator <=(_CONST prevariable& v1) _CONST;
   int operator >=(_CONST prevariable& v1) _CONST;
   int operator >(_CONST prevariable& v1) _CONST;
   int operator <(_CONST prevariable& v1) _CONST;
   int operator !=(_CONST prevariable& v1) _CONST;

   int operator ==( double v1) _CONST;
   int operator <=( double v1) _CONST;
   int operator >=( double v1) _CONST;
   int operator >( double v1) _CONST;
   int operator <( double v1) _CONST;
   int operator !=( double v1) _CONST;
#if defined(USE_DDOUBLE)
   int operator ==( int v1) _CONST;
   int operator <=( int v1) _CONST;
   int operator >=( int v1) _CONST;
   int operator >( int v1) _CONST;
   int operator <( int v1) _CONST;
   int operator !=( int v1) _CONST;
#endif

  public:
#ifdef __SUN__
    prevariable(void) { }
#endif
#ifdef __NDPX__
    prevariable(double_and_int * u) {v=u;}
#endif

    void initialize(void);

    friend char* fform(const char*,_CONST prevariable&);

    void operator +=(_CONST prevariable& t1);
    void operator +=( double t1);

    void operator -=(_CONST prevariable& t1);
    void operator -=( double t1);

    void operator /=(_CONST prevariable& v1);
    void operator /=( double v1);

    void operator *=(_CONST prevariable& v1);
    void operator *=( double v1);

    friend prevariable& operator +(_CONST prevariable& v1,_CONST prevariable& v2);

    friend prevariable& operator +( double v1,_CONST prevariable& v2);

    friend prevariable& operator +(_CONST prevariable& v1, double v2);

    friend prevariable& operator -(_CONST prevariable& v1);

    friend prevariable& operator -(_CONST prevariable& v1,_CONST prevariable& v2);

    friend prevariable& operator -( double v1,_CONST prevariable& v2);

    friend prevariable& operator -(_CONST prevariable& v1, double v2);

    friend prevariable& pow(_CONST prevariable& t1,_CONST prevariable& t2);

    friend prevariable& pow(_CONST prevariable& t1, double t2);

    friend prevariable& pow( double t1,_CONST prevariable& t2);

  }; // end of class prevariable

  inline double& value(_CONST prevariable& v1) { return(v1.v->x); }
  inline double * address(_CONST prevariable& v1) { return(&(v1.v->x)); }


  prevariable& operator << (BOR_CONST prevariable& v1,_CONST prevariable& v2);
  dvar_vector& operator << (BOR_CONST dvar_vector& v1,_CONST dvar_vector& v2);
  dvar_matrix& operator << (BOR_CONST dvar_matrix& v1,_CONST dvar_matrix& v2);

  class df1_one_variable;

  class dvariable : public prevariable
  {
  public:
     dvariable();
     ~dvariable();
     dvariable( double t);
     dvariable(BOR_CONST int& t);
     dvariable(kkludge_object);
     dvariable(_CONST prevariable&);
     dvariable& operator=(_CONST prevariable&);
     dvariable& operator = (const df1_one_variable& v);
     dvariable& operator=( double);
#if defined(USE_DDOUBLE)
#  undef double
     dvariable& operator=( double);
#  define double dd_real
#endif
     dvariable(_CONST dvariable&);
//#  if (__BORLANDC__  > 0x0520) 
//     dvariable& operator+=(_CONST prevariable&);
//#  endif

  }; // end of class dvariable

#if defined(max)
#undef max
#endif
#if defined(min)
#undef min
#endif

  class funnel_dvariable : public dvariable
  {
  public:
     dvariable& operator=(_CONST prevariable&);
  };

 // spreadsheet like element wise operations

  dvar_vector operator * (_CONST dvar_vector& t1, double x);
  dvar_vector operator / ( double x,_CONST dvar_vector& t1);
  dvar_vector operator / (_CONST dvar_vector& t1, double x);
  dvar_vector operator + ( double x,_CONST dvar_vector& t1);
  dvar_vector operator + (_CONST dvar_vector& t1, double x);
  dvar_vector operator - ( double x,_CONST dvar_vector& t1);
  dvar_vector operator - (_CONST dvar_vector& t1, double x);
  dvar_vector operator - (_CONST dvar_vector& t1);
  dvar_vector operator * (_CONST dvar_vector& t1,_CONST prevariable& x);
  dvar_vector operator / (_CONST prevariable& x,_CONST dvar_vector& t1);
  dvar_vector operator / (_CONST dvar_vector& t1,_CONST prevariable& x);
  dvar_vector operator + (_CONST prevariable& x,_CONST dvar_vector& t1);
  dvar_vector operator + (_CONST dvar_vector& t1,_CONST prevariable& x);
  dvar_vector operator - (_CONST prevariable& x,_CONST dvar_vector& t1);
  dvar_vector operator - (_CONST dvar_vector& t1,_CONST prevariable& x);
  dvar_vector operator - (_CONST dvector& t1,_CONST prevariable& x);
  dvar_vector operator * (_CONST dvector& t1,_CONST prevariable& x);
  dvar_vector operator * (_CONST prevariable& x,_CONST dvector& t1);

  dvector operator * (_CONST dvector& t1, double x);
  dvector operator / ( double x,_CONST dvector& t1);
  dvector operator / (_CONST dvector& t1, double x);
  dvector operator + ( double x,_CONST dvector& t1);
  dvector operator + (_CONST dvector& t1, double x);
  dvector operator - ( double x,_CONST dvector& t1);
  dvector operator - (_CONST  dvector& t1, double x);
  dvector operator - (_CONST dvector& t1);

  double max(_CONST dmatrix&);
  int max(const imatrix&);
  double max(_CONST dvector&);
  dvariable max(_CONST dvar_vector&);
  dvariable min(_CONST dvar_vector&);

  dmatrix symmetrize(_CONST dmatrix& m1);
  dvector eigenvalues(_CONST dmatrix& m1);
  dmatrix eigenvectors(_CONST dmatrix& m1);
  dmatrix eigenvectors(_CONST dmatrix& m1,BOR_CONST dvector& eigenvalues);

  dvar_matrix symmetrize(_CONST dvar_matrix& m1);
  dvar_vector eigenvalues(_CONST dvar_matrix& m1);
  dvar_matrix eigenvectors(_CONST dvar_matrix& m1);

  dmatrix outer_prod(_CONST dvector& t1 ,_CONST dvector& t2 );
  dvar_matrix outer_prod(_CONST dvar_vector& t1 ,_CONST dvar_vector& t2 );
  dvar_matrix outer_prod(_CONST dvector& t1 ,_CONST dvar_vector& t2 );
  dvar_matrix outer_prod(_CONST dvar_vector& t1 ,_CONST dvector& t2 );
  //double trace(_CONST dmatrix&);
  //dvariable trace(_CONST dvar_matrix&);
  dmatrix operator * ( double x,_CONST dmatrix& m);
  dmatrix operator * (_CONST dmatrix& m, double d);
  dmatrix operator / (_CONST dmatrix& m, double d);
  dmatrix operator / ( double d,_CONST dmatrix& m);
  dmatrix operator + ( double x,_CONST dmatrix& m);
  dvar_matrix operator + ( const dvariable& x,const dmatrix& m);
  dvar_matrix operator - ( const dvariable& x,const dmatrix& m);
  dmatrix operator + (_CONST dmatrix& m, double d);
  dmatrix operator - ( double x,_CONST dmatrix& m);
  dmatrix operator - (_CONST dmatrix& m, double d);
  dvar_matrix operator / (_CONST dvar_matrix& m,_CONST prevariable& x);
  dvar_matrix operator / (_CONST dmatrix& m,_CONST prevariable& x);
  dvar_matrix operator / (_CONST dvar_matrix& m, double x);
  dvar_matrix operator / ( double x,_CONST dvar_matrix& m);
  dvar_matrix operator / (_CONST prevariable& x,_CONST dvar_matrix& m);

  dvar_matrix operator * (_CONST prevariable& x,_CONST dmatrix& m);
  dvar_matrix operator * (_CONST dvar_matrix& m,_CONST prevariable& x);
  dvar_matrix operator * (_CONST prevariable& x,_CONST dvar_matrix& m);
  dvar_matrix operator * ( double x,_CONST dvar_matrix& m);

  dvector operator & (_CONST dvector& t1,_CONST dvector& t2);
  dvar_vector operator & (_CONST dvar_vector& t1,_CONST dvar_vector& t2);


  ivector column(_CONST imatrix& m,int i);
  dvector extract_column(_CONST dmatrix& m,int i);
  dvector column(_CONST dmatrix& m,int i);
  dvector extract_row(_CONST dmatrix& m, int j);
  dvector row(_CONST dmatrix& m, int j);
  dvar_vector extract_column(_CONST dvar_matrix& m, int i);
  dvar_vector column(_CONST dvar_matrix& m, int i);
  dvector column_value(const dvar_matrix& m, int i);
  dvar_vector extract_row(_CONST dvar_matrix& m, int j);
  dvar_vector row(_CONST dvar_matrix& m, int j);

// dvector mathematical functions

     dvector sin(_CONST dvector& t1);
     dvector sqrt(_CONST dvector& t1);
     dvector sqr(_CONST dvector& t1);
     dvector exp(_CONST dvector& t1);
     dvector mfexp(_CONST dvector& t1);
     dvector mfexp(_CONST dvector& t1, double d);
     dvector atan(_CONST dvector& t1);
     dvector tan(_CONST dvector& t1);
     dvector tanh(_CONST dvector& t1);
     dvector atan2(_CONST dvector& t1,_CONST dvector& t2);
     dvector atan2(_CONST dvector& t1, double t2);
     dvector atan2( double t1,_CONST dvector& t2);
     dvector acos(_CONST dvector& t1);
     dvector asin(_CONST dvector& t1);
     dvector cos(_CONST dvector& t1);
     dvector sinh(_CONST dvector& t1);
     dvector cosh(_CONST dvector& t1);
     dvector log(_CONST dvector& t1);
     dvector log10(_CONST dvector& t1);
     dvector pow(_CONST dvector& t1, double);
     dvector pow(_CONST dvector& t1,int);
     dvector pow( double,_CONST dvector& t1);
     ivector pow(const ivector& v1, int x); 
     ivector pow(int x,const ivector& v1); 

// end of dvector mathematical functions

// dvar_vector mathematical functions

     dvar_vector sin(_CONST dvar_vector& t1);
     dvar_vector sqrt(_CONST dvar_vector& t1);
     dvar_vector sqr(_CONST dvar_vector& t1);
     dvar_vector exp(_CONST dvar_vector& t1);
     dvar_vector mfexp(_CONST dvar_vector& t1);
     dvar_vector mfexp(_CONST dvar_vector& t1, double d);
     dvar_vector atan(_CONST dvar_vector& t1);
     dvar_vector tan(_CONST dvar_vector& t1);
     dvar_vector tanh(_CONST dvar_vector& t1);
     dvar_vector atan2(_CONST dvar_vector& t1,_CONST dvar_vector& t2);
     dvar_vector atan2(_CONST dvar_vector& t1, double t2);
     dvar_vector atan2( double t1,_CONST dvar_vector& t2);
     dvar_vector acos(_CONST dvar_vector& t1);
     dvar_vector asin(_CONST dvar_vector& t1);
     dvar_vector cos(_CONST dvar_vector& t1);
     dvar_vector sinh(_CONST dvar_vector& t1);
     dvar_vector cosh(_CONST dvar_vector& t1);
     dvar_vector log(_CONST dvar_vector& t1);
     dvar_vector log10(_CONST dvar_vector& t1);
     dvar_vector pow(_CONST dvar_vector&,_CONST dvar_vector& t1);
     dvar_vector pow(_CONST dvar_vector&,_CONST dvector& t1);
     dvar_vector pow(_CONST dvector&,_CONST dvar_vector& t1);
     dvar_vector pow(_CONST dvar_vector& t1, double);
     dvar_vector pow(_CONST dvar_vector& t1,int);
     dvar_vector pow(_CONST dvar_vector& t1,_CONST prevariable&);
     dvar_vector pow(_CONST dvector& t1,_CONST prevariable&);
     dvar_vector pow(_CONST prevariable&,_CONST dvar_vector& t1);

// end of dvar_vector mathematical functions

// dmatrix mathematical functions
    dmatrix exp(_CONST dmatrix& m);
    dmatrix mfexp(_CONST dmatrix& m);
    dmatrix mfexp(_CONST dmatrix& m, double d);
    dmatrix sqrt(_CONST dmatrix& m);
    dmatrix sqr(_CONST dmatrix& m);
    dmatrix pow(_CONST dmatrix& m, double e);
    dmatrix pow(_CONST dmatrix& m,int e);
    dmatrix log(_CONST dmatrix& m);
    dmatrix sin(_CONST dmatrix& m);
    dmatrix cos(_CONST dmatrix& m);
    dmatrix tan(_CONST dmatrix& m);
    dmatrix elem_div(_CONST dmatrix& m,_CONST dmatrix& m2);
    dmatrix elem_prod(_CONST dmatrix& m,_CONST dmatrix& m2);
// end of dmatrix mathematical functions

//  dvar_matrix mathematical functions
    dvar_matrix exp(_CONST dvar_matrix& m);
    dvar_matrix mfexp(_CONST dvar_matrix& m);
    dvar_matrix mfexp(_CONST dvar_matrix& m, double d);
    dvar_matrix sqrt(_CONST dvar_matrix& m);
    dvar_matrix sqr(_CONST dvar_matrix& m);
    dvar_matrix log(_CONST dvar_matrix& m);
    dvar_matrix sin(_CONST dvar_matrix& m);
    dvar_matrix cos(_CONST dvar_matrix& m);
    dvar_matrix tan(_CONST dvar_matrix& m);
    dvar_matrix pow(_CONST dvar_matrix& m, double e);
    dvar_matrix pow(_CONST dvar_matrix& m,_CONST prevariable& e);
    dvar_matrix pow(_CONST dmatrix& m,_CONST prevariable& e);
    dvar_matrix pow(_CONST dvar_matrix& m,int e);
    dvar_matrix elem_prod(_CONST dvar_matrix& m,_CONST dvar_matrix& m2);
    dvar_matrix elem_prod(_CONST dvar_matrix& m,_CONST dmatrix& m2);
    dvar_matrix elem_prod(_CONST dmatrix& m,_CONST dvar_matrix& m2);
    dvar_matrix elem_div(_CONST dvar_matrix& m,_CONST dvar_matrix& m2);
    dvar_matrix elem_div(_CONST dvar_matrix& m,_CONST dmatrix& m2);
    dvar_matrix elem_div(_CONST dmatrix& m,_CONST dvar_matrix& m2);
//  end of dvar_matrix mathematical functions

  int min(_CONST ivector& t1);
  int max(_CONST ivector& t1);
  int Max(_CONST ivector& t1);

     double mfexp(double);
     double mfexp(double,double bound);
     dvariable mfexp(const prevariable&);
     dvariable mfexp(const prevariable&,double bound);

#ifndef DVEC_ARR_HPP
#define DVEC_ARR_HPP

  class vector_shapex 
  {
  public:
    void * trueptr;
    vector_shapex(int lb,int ub,void * p) : index_min(lb), 
      index_max(ub), ncopies(0), trueptr(p) {}
    void * get_truepointer(void){ return trueptr; }
    friend class dvector;
    friend class ivector;
    //friend class tdvector;
    friend class dvar_vector;
#if defined(USE_VECTOR_SHAPE_POOL)
    static vector_shape_pool * xpool;
    void * operator new(size_t);
    void operator delete(void * ptr,size_t n)
    {  xpool->free(ptr); }
#endif
    unsigned int ncopies;
    void shift(int min);
    int index_min;
    int index_max;
  private:
    friend class subdvector;
    friend class lvector;
    friend class ptr_vector;
  public:
    int decr_ncopies(void) {return --ncopies;}
    int get_ncopies(void) {return ncopies;}
    int incr_ncopies(void) {return ++ncopies;}
    int indexmin(){return index_min;}
    int indexmax(){return index_max;}
  };
  class predvector
  {
    dvector * p;
    int lb;
    int ub;
    inline predvector(dvector * _p,int _lb,int _ub) {p=_p;lb=_lb,ub=_ub;}
    friend class dvector;
  };

  class predvar_vector
  {
    dvar_vector * p;
    int lb;
    int ub;
    inline predvar_vector(dvar_vector * _p,int _lb,int _ub) {p=_p;lb=_lb,ub=_ub;}
    friend class dvar_vector;
  };
  class dvector 
  {
  protected:
    double * v;
    int index_min;
    int index_max;
    vector_shapex * shape;

  public:
    dvector operator -();
    int operator ! (void) _CONST {return (shape == NULL);}
    int allocated(void) _CONST {return (shape != NULL);}
    dvector& operator -- (void) 
    { 
      index_min--;index_max--;v++; return *this;
    }
    dvector& operator ++ (void)  
    { 
      index_min++;index_max++;v--; return *this;
    }
    void reallocate(double size);
  dvector sub(int lb,int ub)
  {
    return predvector(this,lb,ub);
  }
  dvector operator () (int lb,int ub)
  {
    return predvector(this,lb,ub);
  }
  dvector(const predvector& pd);
  //virtual void write_message(void) { cout << " This is a dvector" << endl; }
  
  void fill_randpoisson(double lambda, BOR_CONST random_number_generator& rng);
  void fill_randnegbinomial(double lambda,double tau,
   const random_number_generator& rng);
  void shallow_copy(const dvector&);

    friend class banded_symmetric_dvar_matrix;
    friend class banded_lower_triangular_dvar_matrix;
    friend class banded_symmetric_dmatrix;
    friend class banded_lower_triangular_dmatrix;
    
//#ifdef HOME_VERSION
      void allocate(int ncl,int ncu);
      void allocate(_CONST dvector& dv);
      void allocatec(_CONST dvector& dv);
      void allocate(_CONST dvar_vector&);
      void allocate(const char *);
      void allocate(void);
      void deallocate(void);
      void safe_deallocate(void);
      void safe_allocate(int,int);
//#endif

    _CONST double * address() _CONST { return (v); }
    double*& get_v(void) { return (v); }
    double& elem(int i) { return(v[i]); }
    double* initpointer(void) { return (v+indexmin()); }
    #ifdef USE_CONST
      _CONST double* initpointer(void) _CONST { return (v+indexmin()); }
      double* get_v(void) _CONST { return (v); }
      _CONST double& elem(int i) _CONST { return(v[i]); }
    #endif

    void fill(const char * s);
    void fill_randu_ni(long int& n);
    void fill_randn_ni(long int& n);
    void fill_randbi_ni(long int& n,double);

    void fill_randu(long int& n);
    void fill_randn(long int& n);
    void fill_randbi(long int& n,double);


    void fill_randbi(double p,BOR_CONST random_number_generator& rng);
    void fill_randu(BOR_CONST random_number_generator& rng);
    void fill_randn(BOR_CONST random_number_generator& rng);
    void fill_randcau(BOR_CONST random_number_generator& rng);

    void fill_seqadd( double, double);
    void fill_multinomial(BOR_CONST int& seed,_CONST dvector& p);
    void fill_multinomial(BOR_CONST random_number_generator& rng,_CONST dvector& p);
    void initialize(void);

    int& testmin()  {return shape->index_min;} // returns the minimum allowable index

    int& testmax() {return shape->index_max;} // returns the minimum allowable index

    int indexmin() _CONST {return index_min;} // returns the minimum allowable index

    int indexmax() _CONST {return index_max;} // returns the maximum allowable index
    int size() _CONST {return (index_max - index_min+1);} // returns the number of elements
    dvector& shift(int min);

  
    dvector(BOR_CONST dvar_vector_position& dvp,BOR_CONST kkludge_object&);

    dvector(const ad_integer&,const index_type&);
    void allocate(const ad_integer&,const index_type&);
    dvector(void);

    dvector(_CONST dvector&);

    /* Not used.
    //dvector(_CONST dvector&,int lb,int ub);
    */

    dvector(_CONST ivector&);

    dvector(_CONST lvector&);

    dvector(const char *);

    dvector( int ncl,  int ncu);
                                           // makes an array [ncl..ncu]

    dvector( unsigned int sz, double * x);

    dvector(char * filename, const int& column);

    //operator double* () { return v;}

   ~dvector();

    void save_dvector_position(void) _CONST;
    void save_dvector_derivatives(BOR_CONST dvar_vector_position& pos) _CONST;
    void save_dvector_derivatives_na(BOR_CONST dvar_vector_position& pos) _CONST;
    void save_dvector_derivatives(void) _CONST;
    void save_dvector_value(void) _CONST;

    //dvector operator()(int,int);
    dvector operator()(_CONST lvector&);
    dvector operator ()(_CONST ivector& u);

    dvector& operator += (_CONST dvector& v1);
    dvector& operator -= (_CONST dvector& v1);
    dvector& operator += ( double v1);
    dvector& operator /= ( double v1);
    dvector& operator *= ( double v1);
    dvector& operator -= ( double v1);

    void read_from(BOR_CONST uistream&);

    friend class sdmatrix;
    friend double norm(_CONST dvector&);
    friend double norm2(_CONST dvector&);
    friend class dvar_vector;
    friend class dmatrix;
    friend class d3_array;
    friend char* fform(const char*,_CONST dvector&);

    void write_on(BOR_CONST ostream&) _CONST;
    void write_on(BOR_CONST uostream&) _CONST;
    void read_from(BOR_CONST istream&);
    friend double operator * (_CONST dvector& ,_CONST dvector& );

    friend dvariable operator * (_CONST dvector& ,_CONST dvar_vector& );

    friend dvariable operator * (_CONST dvar_vector& ,_CONST dvector& );

    friend dvar_vector operator * (_CONST prevariable& ,_CONST dvector& );

    friend dvector operator * ( double ,_CONST dvector& );

    friend dvector operator + (_CONST dvector& ,_CONST dvector& );

    friend dvector elem_prod(_CONST dvector& ,_CONST dvector& );

    friend dvector first_difference(_CONST dvector&);
    friend dvector second_difference(_CONST dvector&);

    friend dvector elem_div(_CONST dvector& ,_CONST dvector& );

    friend dvar_vector elem_prod(_CONST dvector& ,_CONST dvar_vector& );

    friend dvar_vector elem_div(_CONST dvector&,_CONST dvar_vector&);

    friend dvar_vector elem_prod(_CONST dvar_vector&,_CONST dvector& );

    friend dvar_vector elem_div(_CONST dvar_vector& ,_CONST dvector& );

    friend dvar_vector operator + (_CONST dvar_vector& ,_CONST dvector& );

    friend dvar_vector operator + (_CONST dvector& ,_CONST dvar_vector& );

    friend dvector operator - (_CONST dvector& ,_CONST dvector& );

    friend dvar_vector operator - (_CONST dvar_vector& ,_CONST dvector& );

    friend dvar_vector operator - (_CONST dvector& ,_CONST dvar_vector& );

    friend dvector  operator * (_CONST dvector& x,_CONST dmatrix& m);

    friend dvector  operator * (_CONST dmatrix& x,_CONST dvector& m);

    friend dvar_vector  operator * (_CONST dvector& x,_CONST dvar_matrix& m);

    friend dvar_vector  operator * (_CONST dvar_matrix& x,_CONST dvector& m);


#ifdef OPT_LIB
      inline double& operator[] (register int i) { return(*(v+i)); }
      inline double& operator() (register int i) { return(*(v+i)); }
#else
      double& operator[] (int i);
      double& operator() (int i);
#endif

#ifdef USE_CONST
  #ifdef OPT_LIB
    inline _CONST double& operator[] (register int i) _CONST { return(*(v+i)); }
    inline _CONST double& operator() (register int i) _CONST { return(*(v+i)); }
  #else
      _CONST double& operator[] (int i) _CONST;
      _CONST double& operator() (int i) _CONST;
  #endif
#endif


    dvector& operator = (_CONST dvector& t);

    dvector& operator = ( double x);

    friend dvector exp(_CONST dvector&);

    friend dvector log(_CONST dvector&);

    friend dvector fabs(_CONST dvector&);

    friend double max(_CONST dvector&);

    friend double min(_CONST dvector&);

 }; // end of class dvector

class independent_variables : public dvector 
{
  public:

    independent_variables(BOR_CONST independent_variables& v) : dvector(v) {}

    independent_variables( int ncl,  int ncu) :
      dvector(ncl,ncu) {}
                                           // makes an array [ncl..ncu]

    independent_variables( unsigned int sz, double * x) :
      dvector(sz, x) {}

    independent_variables& operator = (_CONST dvector& t);

};
     

#endif

dvariable dfatan1( dvariable , double , double ,BOR_CONST prevariable&fpen);

double boundp( double x, double fmin, double fmax,const double& fpen);
double boundp( double x, double fmin, double fmax);

dvariable boundp(const prevariable & x, double fmin, double fmax,const prevariable& fpen);
dvariable boundp(const prevariable & x, double fmin, double fmax,const prevariable& fpen,double s);

double nd2fboundp(double x,double minb,double maxb,const double& pen);
double boundpin(double x, double fmin, double fmax);
double boundpin(const prevariable& x, double fmin, double fmax);
double boundpin(const prevariable& x, double fmin, double fmax,double s);

double dmin(double, double );

double dmax(double i, double j);




#include <stdlib.h>
#ifdef __TURBOC__
#  if !defined(__linux__)
#    include <alloc.h>
#  endif
#endif


   double sigmoid( double t1);

  class mat_shape 
  {
    unsigned int ncopies;
    unsigned int nrows;
    unsigned int ncols; 
    int row_min;
    int row_max;
    int col_min;
    int col_max;
    mat_shape(int rl,int ru,int cl=0,int cu=-1);
    mat_shape(){};
    void colshift(int min);
    void rowshift(int min);

    //friend class const_dmatrix;
    friend class dmatrix;
    friend class sdmatrix;
    friend class dvar_matrix;
    friend class imatrix;
    friend class lmatrix;
    friend class i3_array;
  };

  class mat_shapex 
  {
  public:
    void * trueptr;
    unsigned int ncopies;
    mat_shapex(_CONST void * m) { trueptr=(void*)m;ncopies=0;}
    mat_shapex(){trueptr=NULL;ncopies=0;};

    void * get_pointer (void) {return trueptr;}
    friend class dmatrix;
    friend class sdmatrix;
    friend class dvar_matrix;
    friend class imatrix;
    friend class i3_array;
    friend class lmatrix;
  };

    class arr_link;
    class arr_list
    {
      arr_link * last;
      arr_link * free_last;
      unsigned long int last_offset;
      unsigned long int max_last_offset;
    public:
      long int number_arr_links;
      friend class arr_link;
  
    public:
    
      arr_list(void)
      {
        last = 0;
        free_last = 0;
        last_offset = 0;
        max_last_offset = 0;
        number_arr_links = 0;
      }
      unsigned long int get_last_offset() { return last_offset;}
      unsigned long int get_number_arr_links(){ return(number_arr_links);}
      unsigned long int get_max_last_offset() { return (max_last_offset);}
      void reset_max_last_offset() {max_last_offset=0;}
      friend double_and_int * arr_new(unsigned int);
      friend void arr_free(double_and_int *);
      friend void arr_remove(arr_link **);
      friend void arr_free_list_remove(arr_link **);
      friend void arr_free_add(arr_link *);
      friend void arr_free_remove(arr_link *);
    };

    class arr_link
    {
#if defined(USE_VECTOR_SHAPE_POOL)
    static vector_shape_pool * xpool;
    void * operator new(size_t);
    void operator delete(void * ptr,size_t n) 
    {  xpool->free(ptr); }
#endif
       arr_link *          prev;
       arr_link *          next;
       arr_link *          free_prev;
       arr_link *          free_next;
       unsigned int        status;
       // unsigned int     free_list_status;
       unsigned int        size;
       unsigned long int   offset;
     public:
       arr_link();
    
       friend double_and_int * arr_new(unsigned int);
       friend void arr_free(double_and_int *);
       friend void arr_remove(arr_link **);
       friend void arr_free_remove(arr_link *);
       friend void arr_free_add(arr_link *);
    };

#if defined(__NUMBERVECTOR__)
  class param_init_number_vector;
  class param_init_bounded_number_vector;
  class param_init_vector_vector;
  class param_init_bounded_vector_vector;
#endif

 class dvar_vector
 {
 public:
   double_and_int * va;
   int index_min;
   int index_max;
   arr_link * link_ptr;
   vector_shapex * shape;
 public:
    dvar_vector operator -();

    dvar_vector& operator -- (void) 
    { 
      index_min--;index_max--;va++; return *this;
    }
    dvar_vector& operator ++ (void)  
    { 
      index_min++;index_max++;va--; return *this;
    }
  dvar_vector sub(int lb,int ub)
  {
    return predvar_vector(this,lb,ub);
  }
  dvar_vector operator () (int lb,int ub)
  {
    return predvar_vector(this,lb,ub);
  }
  void shallow_copy(const dvar_vector&);
 int operator ! (void) _CONST {return (shape == NULL);}
 friend class dvar_matrix;
 friend class dvar3_array;
    friend class banded_symmetric_dvar_matrix;
    friend class banded_lower_triangular_dvar_matrix;
    friend class banded_symmetric_dmatrix;
    friend class banded_lower_triangular_dmatrix;

 void fill_randpoisson(double lambda, BOR_CONST random_number_generator& rng);
 void fill_randnegbinomial(double lambda,double tau,
  const random_number_generator& rng);
 prevariable elem(int i) { return (va+i); }

 double& elem_value(int i) { return (va[i].x); }

 double_and_int * get_va(){return va;}

#ifdef USE_CONST
  prevariable elem(int i) _CONST { return (va+i); }

  double& elem_value(int i) _CONST { return (va[i].x); }

  double_and_int * get_va() _CONST {return va;}
#endif


    friend dvar_matrix  operator * (_CONST dvar_matrix& m1,_CONST dmatrix& m2 );

    void deallocate();
    dvar_vector(_CONST dvar_vector&);
    dvar_vector(const predvar_vector&);
    dvar_vector();
    dvar_vector(int ncl,int ncu); // makes an array [ncl..ncu]
    dvar_vector(int ncl,int ncu,kkludge_object);

    //dvar_vector(const ad_integer&,const ad_integer&);
    dvar_vector( unsigned int sz, double * x);
    dvar_vector(BOR_CONST independent_variables&);
    friend char* fform(const char*,_CONST dvar_vector& );
#   if defined(__NUMBERVECTOR__)
    dvar_vector(const param_init_number_vector&);
    dvar_vector(const param_init_bounded_number_vector&);
#   endif
    dvar_vector(_CONST dvector&);
    dvar_vector(const char *);
   ~dvar_vector();
    void allocate(int,int);
    void allocate(void);
    void allocate(_CONST dvector&);
    void allocatec(_CONST dvar_vector&);
    void allocate(_CONST dvar_vector&);

    void allocate(const ad_integer&,const ad_integer&);
    void initialize(_CONST dvector& ww);
    void initialize(void);
    void save_dvar_vector_position(void) _CONST;
    void save_dvar_vector_value(void) _CONST;
    void write_on(BOR_CONST ostream&) _CONST;
    void write_on(BOR_CONST uostream&) _CONST;
    void read_from(BOR_CONST istream&);
    void read_from(BOR_CONST uistream&);
    int indexmin() _CONST {return index_min;} // returns the minimum allowable index
    int indexmax() _CONST {return index_max;} // returns the maximum allowable index
    int size() _CONST {return (index_max - index_min+1);} // returns the number of elements
    dvar_vector& shift(int min);

#ifdef OPT_LIB
 #if defined(__NDPX__) || defined(__SUN__)
  inline prevariable operator() (register int i) { return (va+i); }
  inline prevariable operator[] (register int i) { return (va+i); }
      #ifdef USE_CONST
  inline _CONST prevariable operator() (int i) _CONST { return (va+i); }
  inline _CONST prevariable operator[] (int i) _CONST { return (va+i); }
      #endif
 #else
  inline prevariable operator() (int i) { return (va+i); }
  inline prevariable operator[] (int i) { return (va+i); }
  #ifdef USE_CONST
  inline _CONST prevariable operator() (int i) _CONST { return (va+i); }
  inline _CONST prevariable operator[] (int i) _CONST { return (va+i); }
  #endif
 #endif
#else
  prevariable operator[] (int i);
  prevariable operator() (int i);
  #ifdef USE_CONST
  _CONST prevariable operator[] (int i) _CONST;
  _CONST prevariable operator() (int i) _CONST;
  #endif
#endif

    double* initpointer(void) { return ((double*)(va+indexmin())); }
  #ifdef USE_CONST
    _CONST double* initpointer(void) _CONST { return ((double *)(va+indexmin())); }
#endif
    dvar_vector operator()(_CONST lvector&);
    //dvar_vector operator()(int,int);
    dvar_vector operator ()(_CONST ivector& u);
    dvar_vector& operator+= (_CONST prevariable& d);
    dvar_vector& operator+= (double d);
    dvar_vector& operator/= (_CONST prevariable& d);
    //dvar_vector& operator*= (_CONST dvariable& d);
    dvar_vector& operator*= (_CONST prevariable& d);
    dvar_vector& operator*= ( double d);
    dvar_vector& operator/= ( double d);
    dvar_vector& operator-= (_CONST prevariable& d);
    dvar_vector& operator-= (double d);
    dvar_vector& operator+= (_CONST dvector& v1);
    dvar_vector& operator-= (_CONST dvector& v1);
    dvar_vector& operator+= (_CONST dvar_vector& v1);
    dvar_vector& operator-= (_CONST dvar_vector& v1);
    dvar_vector& operator = (_CONST dvar_vector& t);
    dvar_vector& operator = (_CONST dvector& t);
    dvar_vector& operator = ( double t);
    dvar_vector& operator = (_CONST prevariable& t);
    void fill(const char *);
    void fill_randu(long int& n);
    void fill_randn(long int& n);
    void fill_randbi(long int& n,double);

    void fill_randu_ni(long int& n);
    void fill_randn_ni(long int& n);
    void fill_randbi_ni(long int& n,double);

    void fill_seqadd( double, double);
    void fill_multinomial(BOR_CONST int& seed,_CONST dvector& p);
    void fill_multinomial(BOR_CONST random_number_generator& rng,_CONST dvector& p);
    friend dvar_vector operator + (_CONST dvar_vector& ,_CONST dvar_vector&);
    friend dvar_vector operator + (_CONST dvar_vector& ,_CONST dvector&);
    friend dvar_vector operator + (_CONST dvector& ,_CONST dvar_vector&);
    friend dvar_vector operator - (_CONST dvar_vector& ,_CONST dvar_vector&);

    friend dvar_vector operator - (_CONST dvector& ,_CONST dvar_vector&);

    friend dvar_vector operator - (_CONST dvar_vector& ,_CONST dvector&);

    friend dvar_vector sigmoid(_CONST dvar_vector& t1);

    friend dvariable operator * (_CONST dvar_vector& ,_CONST dvar_vector&);

    friend dvar_vector elem_div(_CONST dvar_vector& ,_CONST dvar_vector&);

    friend dvariable operator * (_CONST dvector& ,_CONST dvar_vector&);

    friend dvariable operator * (_CONST dvar_vector& ,_CONST dvector&);

    friend dvar_vector operator * (_CONST prevariable& ,_CONST dvar_vector&);

    friend dvar_vector operator * (_CONST prevariable& ,_CONST dvector& );

    friend dvar_vector operator * ( double ,_CONST dvar_vector&);

    friend dvar_vector operator * (_CONST dvar_vector& ,_CONST dmatrix& );

    friend dvar_vector operator * (_CONST dmatrix& ,_CONST dvar_vector& );

    friend dvar_vector operator * (_CONST dvar_vector& ,_CONST dvar_matrix& );

    friend dvar_vector operator * (_CONST dvar_matrix& ,_CONST dvar_vector& );

    friend dvar_matrix operator * (_CONST dvar_matrix& ,_CONST dvar_matrix& );

    friend dvar_matrix operator * (_CONST dmatrix& , _CONST dvar_matrix& );

    friend dvar_vector elem_prod(_CONST dvar_vector&,_CONST dvar_vector& );

    friend dvar_vector first_difference(_CONST dvar_vector&);
    friend dvar_vector second_difference(_CONST dvar_vector&);

    //friend dvar_vector elem_div(_CONST dvar_vector& ,_CONST dvar_vector& ); // js, see above

    friend dvar_vector elem_prod(_CONST dvector& ,_CONST dvar_vector& );

    friend dvar_vector elem_div(_CONST dvector&,_CONST dvar_vector&);

    friend dvar_vector elem_prod(_CONST dvar_vector&,_CONST dvector& );

    friend dvar_vector elem_div(_CONST dvar_vector& ,_CONST dvector& );

    friend dvariable norm(_CONST dvar_vector&);
    friend dvariable norm2(_CONST dvar_vector&);


    friend void copy_status(BOR_CONST ostream& s,_CONST dvar_vector& v);

    friend dvar_vector exp(_CONST dvar_vector&);

    friend dvar_vector log(_CONST dvar_vector&);

    friend dvar_vector sin(_CONST dvar_vector&);

    friend dvar_vector fabs(_CONST dvar_vector&);

    friend dvector value(_CONST dvar_vector& v1);

    friend dvar_vector sfabs(_CONST dvar_vector&);

    friend void make_indvar_list(_CONST dvar_vector&);
    friend class array_size;

 }; // end of class dvar_vector

 /*
 class funnel_dvar_vector : public dvar_vector
 {
 public:
    funnel_dvar_vector(int l,int u);
    dvar_vector& operator=(_CONST dvar_vector&);
 };
 */
   class fvar_ptr
   {
     dvar_vector * p;
   };


 class dvar_matrix
 {
    int index_min;
    int index_max;
    dvar_vector * m;
    mat_shapex * shape;

 public:
    dvar_matrix& operator -- (void) 
    { 
      index_min--;index_max--;m++; return *this;
    }
    dvar_matrix& operator ++ (void)  
    { 
      index_min++;index_max++;m--; return *this;
    }
  
    int operator ! (void) _CONST {return (shape == NULL);}
  inline dvar_vector& elem(int i) { return (m[i]); }
  inline prevariable elem(int i, int j) { return (elem(i).elem(j) ) ; }
  #ifdef USE_CONST
    inline dvar_vector& elem(int i) const { return (m[i]); }
    inline prevariable elem(int i, int j) const {return (elem(i).elem(j) ) ; }
  #endif

    friend class banded_symmetric_dvar_matrix;
    friend class banded_lower_triangular_dvar_matrix;
    friend class banded_symmetric_dmatrix;
    friend class banded_lower_triangular_dmatrix;
friend class dvar3_array;
//#ifdef HOME_VERSION
  void shallow_copy(const dvar_matrix&);
  dvar_matrix();
  void allocate(int nrl,int nrh,int ncl,int nch);
  void allocate(int nrl,int nrh);
  void allocate(ad_integer nrl,ad_integer nrh);
  void allocate(_CONST dmatrix& m1);
  void allocate(_CONST dvar_matrix& m1);
  void allocate(int nrl,int nrh,_CONST ivector& ncl,_CONST ivector& nch);
  void allocate(int nrl,int nrh,int ncl,_CONST ivector& nch);
  void allocate(int nrl,int nrh,_CONST ivector& ncl,int nch);
  void allocate(void);
  void deallocate();
//#endif
  dvar_matrix(_CONST banded_symmetric_dvar_matrix&v);
  dvar_matrix(_CONST banded_lower_triangular_dvar_matrix&v);
# if defined(__NUMBERVECTOR__)
    dvar_matrix(const param_init_vector_vector& );
    dvar_matrix(const param_init_bounded_vector_vector&);
# endif
  dvar_matrix sub(int,int);

    double fill_seqadd( double, double);

    int colmin(void) _CONST { return((*this)(indexmin()).indexmin());}
    int colmax(void) _CONST { return((*this)(indexmin()).indexmax());}
    int rowmin(void) _CONST { return(index_min);}
    int rowmax(void) _CONST { return(index_max);}
    int indexmin(void) _CONST { return(index_min);}
    int indexmax(void) _CONST { return(index_max);}
    int rowsize() _CONST {return (rowmax()-rowmin()+1);} // returns the number of rows
    int colsize() _CONST {return (colmax()-colmin()+1);} // returns the number of columns
    void colshift(int min);
    void rowshift(int min);

    friend char* fform(const char*,_CONST dvar_matrix&);

    friend class dvar_vector;

    dvar_matrix(const ad_integer& nrl,const ad_integer& nrh,
      const index_type& ncl,const index_type& nch);

    void allocate(const ad_integer& nrl,const ad_integer& nrh,
      const index_type& ncl,const index_type& nch);

    dvar_matrix( int,  int,  int,  int );
    dvar_matrix( int,  int);
    dvar_matrix( int,  int, kkludge_object kk);
             // makes a matrix [nrl..nrh][ncl..nch]

    dvar_matrix(int,int,_CONST ivector&,_CONST ivector&);
             // makes a ragged dvar_matrix [nrl..nrh][ncl..nch]

    dvar_matrix(int,int,int,_CONST ivector&);
             // makes a ragged dvar_matrix [nrl..nrh][ncl..nch]

    dvar_matrix(_CONST dvar_matrix&);
             // copy initializer
    void initialize(void);

    dvar_matrix(_CONST dmatrix&);

    //dvar_matrix(char *);

   ~dvar_matrix();

    void save_dvar_matrix_position(void) _CONST;
    void save_dvar_matrix_value(void) _CONST;

    void fill(const char *);
    //void colfill(BOR_CONST int&n,...);
    //void rowfill(BOR_CONST int&n,...);

    void colfill_randu(BOR_CONST int&j,long int&n);
    void rowfill_randu(BOR_CONST int& i,long int& n);
    void colfill_randn(BOR_CONST int&j,long int&n);
    void rowfill_randn(BOR_CONST int& i,long int& n);
    void fill_randn(long int&n);
    void fill_randu(long int&n);

    void colfill_seqadd_ni(BOR_CONST int&, double, double);
    void colfill_randu_ni(BOR_CONST int&j,long int&n);
    void rowfill_randu_ni(BOR_CONST int& i,long int& n);
    void colfill_randn_ni(BOR_CONST int&j,long int&n);
    void rowfill_randn_ni(BOR_CONST int& i,long int& n);
    void fill_randn_ni(long int&n);
    void fill_randu_ni(long int&n);

    void colfill_seqadd(BOR_CONST int&, double, double);
    void rowfill_seqadd(BOR_CONST int&, double, double);
    void colfill(int j,_CONST dvar_vector& v);
    void rowfill(int j,_CONST dvar_vector& v);

    void write_on(BOR_CONST ostream&) _CONST;
    void write_on(BOR_CONST uostream&) _CONST;
    void read_from(BOR_CONST istream&);
    void read_from(BOR_CONST uistream&);

#ifdef OPT_LIB
      inline dvar_vector& operator[] (register int i) { return (m[i]); }
      inline dvar_vector& operator() (register int i) { return (m[i]); }
#else
      dvar_vector& operator () (int i);
      dvar_vector& operator[] (int) ;
#endif

#ifdef USE_CONST
  #ifdef OPT_LIB
      inline _CONST dvar_vector& operator[] (register int i) _CONST { return (m[i]); }
      inline _CONST dvar_vector& operator() (register int i) _CONST { return (m[i]); }
  #else
      _CONST dvar_vector& operator () (int i) _CONST;
      _CONST dvar_vector& operator[] (int) _CONST ;
  #endif
#endif

  #ifdef OPT_LIB
    #ifdef __NDPX__
  prevariable operator () (register int i, register int j)
    { return ( prevariable((m[i]).va+j) ); }
    #else
  inline prevariable operator () (register int i, register int j)
    { return ( (m[i]).va+j ); }
    #endif
  #else
      prevariable operator () (int i, int j);
  #endif

  inline double& elem_value (register int i, register int j) 
    { return * (double*) ((m[i]).va+j); }

#ifdef USE_CONST
  inline _CONST double& elem_value (register int i, register int j) _CONST
    { return * (double*) ((m[i]).va+j); }
  #ifdef OPT_LIB
    #ifdef __NDPX__
  prevariable operator () (register int i, register int j) _CONST
    { return ( prevariable((m[i]).va+j) ); }
    #else
  inline prevariable operator () (register int i, register int j) _CONST
    { return ( (m[i]).va+j ); }
    #endif
  #else
      _CONST prevariable operator () (int i, int j) _CONST;
  #endif
#endif

    dvar_matrix& operator+= (_CONST dvar_matrix& x);
    dvar_matrix& operator-= (_CONST dvar_matrix& x);
    dvar_matrix& operator+= (_CONST dmatrix& x);
    dvar_matrix& operator-= (_CONST dmatrix& x);


    dvar_matrix& operator = (_CONST dvar_matrix& );

    dvar_matrix& operator = (_CONST dmatrix& );
    dvar_matrix& operator = ( double t);
    dvar_matrix& operator = (_CONST prevariable& t);

    dvar_matrix& operator *= (_CONST prevariable& t);
    dvar_matrix& operator *= ( double t);
    dvar_matrix& operator /= (_CONST prevariable& t);
    dvar_matrix& operator /= ( double t);

    friend dvar_vector operator * (_CONST dvar_vector& ,_CONST dvar_matrix& );

    friend dvar_vector operator * (_CONST dvar_matrix& ,_CONST dvar_vector& );

    friend dvar_vector operator * (_CONST dvector& ,_CONST dvar_matrix& );

    friend dvar_vector operator * (_CONST dvar_matrix& ,_CONST dvector& );

    friend dvar_matrix operator * (_CONST dvar_matrix& , _CONST dvar_matrix& );

    friend dvar_matrix operator * (_CONST dvar_matrix& , _CONST dmatrix& );

    friend dvar_matrix operator * (_CONST dmatrix& , _CONST dvar_matrix& );

    friend dvar_matrix operator + (_CONST dvar_matrix& ,_CONST dvar_matrix& );
    friend dvar_matrix operator + (_CONST dvar_matrix& ,_CONST dmatrix& );
    friend dvar_matrix operator + (_CONST dmatrix& ,_CONST dvar_matrix& );

    friend dvar_matrix operator + ( double ,_CONST dvar_matrix& );
    friend dvar_matrix operator + (_CONST dvar_matrix& , double );
    friend dvar_matrix operator - ( double ,_CONST dvar_matrix& );
    friend dvar_matrix operator - (_CONST dvar_matrix& , double );

    friend dvar_matrix operator + (_CONST dvariable& ,_CONST dvar_matrix& );
    friend dvar_matrix operator + (_CONST dvar_matrix& ,_CONST dvariable& );
    friend dvar_matrix operator - (_CONST dvariable& ,_CONST dvar_matrix& );
    friend dvar_matrix operator - (_CONST dvar_matrix& ,_CONST dvariable& );

    friend dvar_matrix operator - (_CONST dvar_matrix& ,_CONST dvar_matrix& );
    friend dvar_matrix operator - (_CONST dvar_matrix& ,_CONST dmatrix& );
    friend dvar_matrix operator - (_CONST dmatrix& ,_CONST dvar_matrix& );

    friend dvar_matrix inv(_CONST dvar_matrix&);

    friend dvariable det(_CONST dvar_matrix&);
    friend dvariable ln_det(const dvar_matrix&,const int& sgn);

    //friend dvar_matrix testsub(dvar_matrix);

    friend  dvar_matrix trans(_CONST dvar_matrix&);

    friend dvariable norm(_CONST dvar_matrix&);
    friend dvariable norm2(_CONST dvar_matrix&);

    friend void copy_status(BOR_CONST ostream& s,_CONST dvar_matrix& m1);

 }; // end of class dvar_matrix

dvariable ln_det(const dvar_matrix&);
dvar_matrix operator * (const dvar_matrix& t1, double x);
dmatrix value(_CONST dvar_matrix& m);
d3_array value(_CONST dvar3_array& a);
dvar_vector sort(_CONST dvar_vector&,int NSTACK=60);
dvector sort(_CONST dvector&,int NSTACK=60);
ivector sort(_CONST ivector&,int NSTACK=60);
dvector sort(_CONST dvector&,BOR_CONST ivector& index,int NSTACK=60);
ivector sort(_CONST ivector&,BOR_CONST ivector& index,int NSTACK=60);
dmatrix sort(_CONST dmatrix&,int column,int NSTACK=60);
imatrix sort(_CONST imatrix&,int column,int NSTACK=60);

 //void gradcalc( int , double *);
 void gradcalc(int nvar,BOR_CONST dvector& g);
 void slave_gradcalc(void);

 class dmatrix
 {
 protected:
   int index_min;
   int index_max;
    dvector * m;
    mat_shapex * shape;
    friend char* fform(const char*,_CONST dmatrix&);
    friend class dvar_matrix;
 public:
 
    dmatrix& operator -- (void) 
    { 
      index_min--;index_max--;m++; return *this;
    }
    dmatrix& operator ++ (void)  
    { 
      index_min++;index_max++;m--; return *this;
    }
    void shallow_copy(const dmatrix&);
    int operator ! (void) _CONST {return (shape == NULL);}

    dmatrix sub(int,int);
//#ifdef HOME_VERSION
    dmatrix(void); 
    dmatrix(int,int,kkludge_object); 
    dmatrix(int,int); 
    void allocate(void);
    void allocate(_CONST dmatrix& dm);
    void allocate(_CONST dvar_matrix&);
    void allocate(int nrl,int nrh,int ncl,int nch);
    void allocate(int nrl,int nrh);
    void allocate(ad_integer nrl,ad_integer nrh);
    void allocate(int nrl,int nrh,int ncl,const ivector& nch);
    //void allocate(int nrl,int nrh,
    // const index_type& ncl,const index_type& nch);
    void allocate(int nrl,int nrh,const ivector& ncl,int nch);
    void deallocate();
    void allocate(const ad_integer& nrl,const ad_integer& nrh,
      const index_type& ncl,const index_type& nch);
    void allocate(int nrl,int nrh,const ivector& ncl,const ivector& nch);
//#endif
    friend class banded_symmetric_dmatrix;
    friend class banded_lower_triangular_dmatrix;

    dmatrix( int,  int,  int,  int );
             // makes a matrix [nrl..nrh][ncl..nch]

    dmatrix(const ad_integer& nrl,const ad_integer& nrh,
      const index_type& ncl,const index_type& nch);

    dmatrix(int,int,_CONST ivector& coll,_CONST ivector& colh);
             // makes a ragged dmatrix[nrl..nrh][ncl..nch]

    dmatrix(int,int,int coll,_CONST ivector& colh);
             // makes a ragged dmatrix[nrl..nrh][ncl..nch]

    dmatrix(BOR_CONST dvar_matrix_position&);

    dmatrix(BOR_CONST dmatrix_position&);

    dmatrix(_CONST dmatrix&);
    dmatrix(_CONST banded_symmetric_dmatrix&);
    dmatrix(_CONST banded_lower_triangular_dmatrix&);
    dmatrix(char *);
    void fill(const char *);
    double fill_seqadd( double, double);
    void initialize(void);
             // copy initializer

   ~dmatrix();
    void save_dmatrix_derivatives(BOR_CONST dvar_matrix_position& pos) _CONST;
    void save_dmatrix_derivatives_na(BOR_CONST dvar_matrix_position& pos) _CONST;
    void save_dmatrix_value(void) _CONST;
    void save_dmatrix_position(void) _CONST;
    //void save_dmatrix_derivatives(void);

    int indexmin(void) _CONST { return index_min;}
    int indexmax(void) _CONST { return index_max;}
    int rowmin(void) _CONST { return index_min;}
    int rowmax(void) _CONST { return index_max;}
    int colmin(void) _CONST { return((*this)(indexmin()).indexmin());}
    int colmax(void) _CONST { return((*this)(indexmin()).indexmax());}
    int rowsize() _CONST {return (rowmax()-rowmin()+1);} // returns the number of rows
    int colsize() _CONST {return (colmax()-colmin()+1);} // returns the number of columns
    void rowshift(int min);
    void colshift(int min);

    void write_on(BOR_CONST ostream&) _CONST;
    void write_on(BOR_CONST uostream&) _CONST;
    void read_from(BOR_CONST istream&);
    void read_from(BOR_CONST uistream&);

    //void colfill(BOR_CONST int&n,...);
    //void rowfill(BOR_CONST int&n,...);

    void colfill_randu(BOR_CONST int&j,long int&n);
    void rowfill_randu(BOR_CONST int& i,long int& n);
    void colfill_randn(BOR_CONST int&j,long int&n);
    void fill_randn(long int&n);
    void fill_randu(long int&n);
    void rowfill_randn(BOR_CONST int& i,long int& n);



    void colfill_randu(BOR_CONST int&j,BOR_CONST random_number_generator& rng);
    void rowfill_randu(BOR_CONST int& i,BOR_CONST random_number_generator& rng);
    void fill_randn(BOR_CONST random_number_generator& rng);
    void fill_randcau(BOR_CONST random_number_generator& rng);
    void fill_randu(BOR_CONST random_number_generator& rng);
    void colfill_randn(BOR_CONST int&j,BOR_CONST random_number_generator& rng);
    void rowfill_randn(BOR_CONST int& i,BOR_CONST random_number_generator& rng);

    void colfill_randu_ni(BOR_CONST int&j,long int&n);
    void rowfill_randu_ni(BOR_CONST int& i,long int& n);
    void colfill_randn_ni(BOR_CONST int&j,long int&n);
    void fill_randn_ni(long int&n);
    void fill_randu_ni(long int&n);
    void rowfill_randn_ni(BOR_CONST int& i,long int& n);



    void colfill_seqadd(BOR_CONST int&, double, double);
    void rowfill_seqadd(BOR_CONST int&, double, double);
    void colfill(int j,_CONST dvector& v);
    void rowfill(int j,_CONST dvector& v);

#ifdef OPT_LIB
      inline dvector& operator() (register int i) { return m[i]; }
      inline double& operator() (register int i, register int j) {return(*(m[i].v+j));}
      inline dvector& operator[] (register int i) { return m[i]; }
#else
      double& operator() (int i, int j);
      dvector& operator() (int i);
      dvector& operator[] (int) ;
#endif

#ifdef USE_CONST
  #ifdef OPT_LIB
      inline _CONST dvector& operator() (register int i) _CONST { return m[i]; }
      inline _CONST double& operator() (register int i, register int j) _CONST 
      {return(*(m[i].v+j));}
      inline _CONST dvector& operator[] (register int i) _CONST  { return m[i]; }
#else
      _CONST double& operator() (int i, int j) _CONST;
      _CONST dvector& operator() (int i) _CONST;
      _CONST dvector& operator[] (int)  _CONST;
  #endif
#endif

    inline dvector& elem(int i) { return( *(m+i) ); }
    inline double& elem(int i,int j) { return( *((*(m+i)).v+j) ); }
    #ifdef USE_CONST
      inline _CONST dvector& elem(int i) _CONST { return( *(m+i) ); }
      inline _CONST double& elem(int i,int j) _CONST 
        { return( *((*(m+i)).v+j) ); }
    #endif
    friend class d3_array;
    friend dvector operator * (_CONST dvector& ,_CONST dmatrix& );

    friend dvector operator * (_CONST dmatrix& ,_CONST dvector& );

    friend dvar_vector operator * (_CONST dvar_vector& ,_CONST dmatrix& );

    friend dvar_vector operator * (_CONST dmatrix& ,_CONST dvar_vector& );

    friend dmatrix operator * (_CONST dmatrix& , _CONST dmatrix& );

    friend dvar_matrix operator * (_CONST dvar_matrix& ,_CONST dmatrix& );

    friend dvar_matrix operator * (_CONST dmatrix& , _CONST  dvar_matrix& );

    friend dvar_matrix::dvar_matrix(_CONST dmatrix&);

    friend dmatrix operator - (_CONST dmatrix& ,_CONST dmatrix& );
    friend dmatrix operator + (_CONST dmatrix& ,_CONST dmatrix& );

    friend dvar_matrix operator + (_CONST dvar_matrix& ,_CONST dmatrix& );

    friend dvar_matrix operator + (_CONST dmatrix& ,_CONST dvar_matrix& );

    friend dmatrix  trans(_CONST dmatrix& m1);

    friend dmatrix inv(_CONST dmatrix&);
    friend dmatrix inv(_CONST dmatrix& m1,const double& _ln_det, const int& _sgn);

    friend double det(_CONST dmatrix&);
    friend double ln_det(_CONST dmatrix& m1,BOR_CONST int& sgn);

    friend double norm(_CONST dmatrix&);
    friend double norm2(_CONST dmatrix&);

    dmatrix& operator += (_CONST dmatrix& t);
    dmatrix& operator -= (_CONST dmatrix& t);

    dmatrix& operator = (const dmatrix& t);
    dmatrix& operator = ( double t);