To compile VERA and SAM, obtain Numerical Recipes software from www.nr.com

The specific routines used are those for conjugate gradient optimization  

frprmn.c linmin.c brent.c f1dim.c mnbrak.c

with variable declarations changed from float to double, and the utilities 

nrutil.c  nrutil.h     

Combine these and modify as explained below to create the file 
conjugate_gradient.c

Compilation is performed using the make utility and the makefile by typing

make

------------------------------------------------------------------------------

Protocol:

In frprmn.c linmin.c mnbrak.c brent.c f1dim.c

Change float to double
Change vector to dvector 
Change free_vector to free_dvector


We included these routines in a single file using unix tools:

cat frprmn.c linmin.c mnbrak.c brent.c f1dim.c > nr_conjugate_gradient.c   

and then performed the replacements globally using

sed 's/float/double/g' nr_conjugate_gradient.c > tempfile1
sed 's/vector/dvector/g' tempfile1 > tempfile2
mv tempfile2  nr_conjugate_gradient.c
rm tempfile1

Parameters and modifications in individual routines are given below

---------
frprmnc.c
 
#define ITMAX 10000

--------
limin.c  :

Include declaration

extern float LINMIN_TOL;

Replace 

TOL by LINMIN_TOL in call to brent


--------
mnbrak.c

The objective functions to be extremized are sometimes not defined in
a region of parameter space, for example when the determinant
of the covariance matrix is negative.  The bracketing function mnbrak
will sometimes attempt to choose values in this region. 
For such cases, mnbrak is told to revise its guess, typically by
going a shorter distance in along the extremization direction. 

Errors of this type are detected using the integer expression 
errno, declared in <errno.h>



1.Include in preamble

#include <errno.h>
#include <stdio.h>

2. Include declaration
  double utemp, scale; /* scale is a rescaling factor */

3. Following first instance of

	*fa=(*func)(*ax);
	*fb=(*func)(*bx);

Insert

	scale = 1.;
	while ( errno == EDOM ){
	  scale /= 3.;
	  errno = 0;
	  (*bx) *= scale;
	  *fb = (*func)(*bx);
	}  
	errno = 0;
  
4. Following


	*cx=(*bx)+GOLD*(*bx-*ax);
	*fc=(*func)(*cx);

Insert

	scale = GOLD;
	while ( errno == EDOM ){ 
	  scale /= GOLD;
	  errno = 0;
	  *cx=(*bx)+scale*(*bx-*ax);
	  *fc=(*func)(*cx);
	}
	errno = 0;

5. Following 

 			u=(*cx)+GOLD*(*cx-*bx);
			fu=(*func)(u);
Insert

			scale = GOLD;
			while ( errno == EDOM ){
			  scale /= GOLD;
			  errno = 0;
			  utemp=(*cx)+scale*(*cx-*bx);
			  fu=(*func)(utemp);
			}
			errno = 0;


6. Replace

		} else if ((*cx-u)*(u-ulim) > 0.0) {
			fu=(*func)(u);
			if (fu < *fc) {
				SHFT(*bx,*cx,u,*cx+GOLD*(*cx-*bx))
				SHFT(*fb,*fc,fu,(*func)(u))
			}

with


		} else if ((*cx-u)*(u-ulim) > 0.0) { 
		  fu=(*func)(u);
		  if ( errno == EDOM ){
		    u=(*cx)+GOLD*(*cx-*bx);
		    errno=0;
		    fu=(*func)(u);
		    scale = GOLD;
		    while ( errno == EDOM ){
		      scale /= GOLD;
		      errno = 0;
		      utemp=(*cx)+scale*(*cx-*bx);
		      fu=(*func)(utemp);
		    }
		    errno = 0;
		  }
		  if (fu < *fc) {
		    utemp =  *cx+GOLD*(*cx-*bx); 
		    fu=(*func)(u);
		    scale = GOLD;
		    while ( errno == EDOM ){
		      scale /= GOLD;
		      errno = 0;
		      utemp=(*cx)+scale*(*cx-*bx);
		      fu=(*func)(utemp);
		    }
		    errno = 0;
		    SHFT(*bx,*cx,u,utemp)
		    SHFT(*fb,*fc,fu,(*func)(utemp))
		  }

7. Following 

			u=ulim;
			fu=(*func)(u);

Insert

			scale = 1.;
			while ( errno == EDOM ){ 
			  scale /= 3.;
			  errno = 0;
			  u =(*bx)+ scale * GLIMIT*(*cx-*bx);
			  if ( scale* GLIMIT <= 1 ) printf ("trouble: u is less than cx\n"); 
			  fu = (*func)(u);
			}  
			errno = 0;


8. Following

			u=(*cx)+GOLD*(*cx-*bx);
			fu=(*func)(u);

Insert

			scale = GOLD;
			while ( errno == EDOM ){
			  scale /= GOLD;
			  errno = 0;
			  u=(*cx)+scale*(*cx-*bx);
			  fu=(*func)(u);
			}
			errno = 0;