-
Douglas Guptill authoredDouglas Guptill authored
volume.c 55.26 KiB
#include <stdlib.h>
#include <ctype.h>
#include <math.h>
#include <stdio.h>
/*--------------------------------------------------------------
This file contains six recursive C routines for calculating
the volume and derivatives of a convex polyhedron in n
dimensions. The routines volume, volumef, volumeb and dvda are
fortran callable interfaces to their C counterparts cvolume,
cvolumeb and cdvda.
ROUTINE COMMENT CALLS
cvolume calculates volume of region cvolume
which satisfies Ax <= b.
cvolumef calculates volume of region cvolumef
which satisfies Ax <= b.
(faster version, see below)
cvolumeb calculates volume and the cvolume
derivative d(vol)/db(0).
cdvda calculates derivatives cdvda, cvolumeb,
d(vol)/da(0,j) (j=1,n). & cvolumebj
cdvdaf calculates derivatives cdvdaf, cvolumeb,
d(vol)/da(0,j) (j=1,n). & cvolumebj
(faster version, see below)
cvolumebj calculates partial volumes
and derivatives d(vol)/d(b(j)
(j=1,n). (called by cdvda.)
Here b(j) is the jth element of vector b, and a(i,j) is the
(i,j)th coefficient of the matrix A corresponding to the
ith constraint and the jth dimension.
The two faster versions cvolumef and cdvdaf do not continue down
the recursive tree if b(j) = 0 for any j at any time. This avoids
extra work but restricts the applicability of the routines to
cases where the origin does not pass through any inconsistent
constraints.
Malcolm Sambridge, April 1995.
--------------------------------------------------------------*/
/*--------------------------------------------------------------
ROUTINE: cvolume
This routine calculates the `volume' in dimension n of the region
bounded by a set of m linear inequality constraints of the form
A x <= b, where a has m rows and n columns and is given by a(m,n),
b is the n-vector and is contained in b(m). The recursive formula
of Lasserre (1983) is used. Redundant constraints are allowed
If the inequality constraints are inconsistent then the volume
is returned as zero. If the original polyhedron is unbounded
then a warning is issued and the volume is return as -1.0
(even though it is undefined).
Jean Braun and Malcolm Sambridge, Jan 1995.
Lasserre, J. B., 1983. "An analytical expression and an Algorithm
for the Volume of a Convex Polyhedron in R^n.", JOTA, vol. 39,
no. 3, 363-377.
--------------------------------------------------------------*/
float cvolume (a,b,m,n,mmax,nmax)
int *n, *m, *mmax, *nmax;
float *a, *b;
{
float v,amax,pivot;
int i,j,t,k,l;
int jj,kk;
float *ai, *aj, *ajt, *apjj, *bi;
int kmm,tmm;
int nn,mm,nn_max,mm_max;
int nm1,mm1;
float *ap, *bp;
int firstmin,firstmax;
float bmin,bmax,bb;
float partialv;
nn= *n;
mm= *m;
nn_max= *nmax;
mm_max= *mmax;
/* one-dimensional case (full reduction) */
if (nn == 1)
{
firstmin=0;
firstmax=0;
for (l=0;l<mm;l++)
{
if ( *(a+l) > 0.)
{
bb= *(b+l)/ *(a+l);
if (firstmin==0) {firstmin=1;bmin=bb;}
else if (bb<bmin) bmin=bb;
}
else if ( *(a+l) < 0.)
{
bb= *(b+l)/ *(a+l);
if (firstmax==0) {firstmax=1;bmax=bb;}
else if (bb>bmax) bmax=bb;
}
else if ( *(b+l) < 0.) /* Constraints are inconsistent */
{
printf("Inconsistent constraints found after reduction to n = 1 \n");
return(0.);
}
}
v=0.;
if (firstmin*firstmax == 1) v=bmin-bmax;
else
{
/* printf("Volume is unbounded at 1-D; volume returned as -1\n");*/
return(-1.0);
}
if (v<0.) {v=0.;}
return(v);
}
nm1=nn-1;
mm1=mm-1;
v=0.;
for (i=0;i<mm;i++)
{ ai=a+i;
bi=b+i;
/* find largest pivot */
amax=0.;
for (j=0;j<nn;j++)
if (fabs( *(ai+j*mm_max)) >= amax) {amax= fabs( *(ai+j*mm_max)); t=j;}
tmm=t*mm_max;
pivot=*(ai+tmm);
/* finds contribution to v from this pivot (if not nil) */
if (amax == 0.)
{
/* Constraint is inconsistent */
if ( *(bi) < 0.)
{ printf("Constraint %d is inconsistent\n",i+1); return(0.); }
/* otherwise constraint is redundant */
printf("Constraint %d is redundant\n",i+1);
}
else
{
/* allocate memory */
ap = (float *) malloc((sizeof (float))*nm1*mm1);
bp = (float *) malloc((sizeof (float))*mm1);
/* reduce a and b into ap and bp eliminating variable t and constraint i */
jj=-1;
for (j=0;j<mm;j++)
if (j != i)
{
jj=jj+1;
aj=a+j;
ajt=aj+tmm;
*(bp+jj)= *(b+j) - *(bi) * *(ajt) / pivot;
apjj=ap+jj;
kk=-1;
for (k=0;k<nn;k++)
if (k != t)
{
kk=kk+1;
kmm=k*mm_max;
*(apjj+kk*mm1)= *(aj+kmm)- *(ajt) * *(ai+kmm)/ pivot;
}
}
/* add contribution to volume from volume calculated in smaller dimension */
partialv=cvolume(ap,bp, &mm1, &nm1, &mm1, &nm1);
if(partialv == -1.0)return(-1.0);
v=v+ *(bi)/amax*partialv/nn;
free(ap);
free(bp);
}
}
return(v);
}
/*--------------------------------------------------------------
ROUTINE: volume
A dummy routine used to call cvolume from a fortran routine
Jean Braun and Malcolm Sambridge, Jan 1995.
----------------------------------------------------------------*/
void volume_();
void volume();
void volume_(a,b,m,n,mmax,nmax,result)
int *n, *m, *mmax, *nmax;
float *a, *b;
float *result;
{
volume (a,b,m,n,mmax,nmax,result);
}
void volume (a,b,m,n,mmax,nmax,result)
int *n, *m, *mmax, *nmax;
float *a, *b;
float *result;
{
*result=cvolume(a,b,m,n,mmax,nmax);
}
/*--------------------------------------------------------------
ROUTINE: cvolumeb
This routine calculates the `volume' in dimension n of the region
bounded by a set of m linear inequality constraints of the form
A x <= b, where a has m rows and n columns and is given by a(m,n),
b is the n-vector and is contained in b(m). The recursive formula
of Lasserre (1983) is used. Redundant constraints are allowed and
a warning is issued if any are encountered. If the inequality
constraints are inconsistent then the volume is returned as zero.
If the original polyhedron is unbounded then a warning is issued
and the volume is return as zero (even though it is undefined).
This version also calculates the derivative of the volume with
respect to the parameter b(0) using the simple formula of
Lasserre (1983). If *opt == 2 then only the derivative is
calculated and not the volume.
This routine is a variation from the routine `cvolume'.
Calls are made to routine cvolume.
Malcolm Sambridge, March 1995.
--------------------------------------------------------------*/
float cvolumeb (a,b,m,n,mmax,nmax,opt,dvdb)
int *n, *m, *opt, *mmax, *nmax;
float *a, *b;
float *dvdb;
{
float v,amax,pivot;
int i,j,t,k,l;
int jj,kk;
float *ai, *aj, *ajt, *apjj, *bi;
int kmm,tmm;
int nn,mm,nn_max,mm_max;int nm1,mm1;
int lmin,lmax;
float *ap, *bp;
int firstmin,firstmax;
float bmin,bmax,bb;
float vol;
nn= *n;
mm= *m;
nn_max= *nmax;
mm_max= *mmax;
/* one-dimensional case (full reduction) */
if (nn == 1)
{
firstmin=0;
firstmax=0;
lmax=0;
lmin=0;
for (l=0;l<mm;l++)
{
if ( *(a+l) > 0.)
{
bb= *(b+l)/ *(a+l);
if (firstmin==0) {firstmin=1;bmin=bb;lmin=l;}
else if (bb<bmin) {bmin=bb;lmin=l;}
}
else if ( *(a+l) < 0.)
{
bb= *(b+l)/ *(a+l);
if (firstmax==0) {firstmax=1;bmax=bb;lmax=l;}
else if (bb>bmax) {bmax=bb;lmax=l;}
}
else if ( *(b+l) < 0.)
{
/* Constraints are inconsistent.
Set volume and derivative to zero. */
printf("Inconsistent constraints found after reduction to n = 1 \n");
*dvdb = 0.;
return(0.);
}
}
v=0.;
*dvdb = 0.;
if (firstmin*firstmax == 1) v=bmin-bmax;
else
{
printf("Volume is unbounded; volume returned as -1\n derivatives returned as zero\n");
return(-1.0);
}
if (v<0.) {v=0.;*dvdb=0.;}
else if (v>0. && lmin == 0) *dvdb = 1. / *a;
else if (v>0. && lmax == 0) *dvdb = -1. / *a;
return(v);
}
nm1=nn-1;
mm1=mm-1;
v=0.;
/* perform main loop over constraints */
for (i=0;i<mm;i++)
{
ai=a+i;
bi=b+i;
/* find largest pivot */
amax=0.;
for (j=0;j<nn;j++)
if (fabs( *(ai+j*mm_max)) >= amax) {amax= fabs( *(ai+j*mm_max)); t=j;}
tmm=t*mm_max;
pivot=*(ai+tmm);
/* finds contribution to v from this pivot (if not nil) */
if (amax == 0.)
{
/* Constraint is inconsistent */
if ( *(bi) < 0.)
{
printf("Constraint %d is inconsistent\n",i+1);
return(0.);
}
/* otherwise constraint is redundant */
printf("Constraint %d is redundant\n",i+1);
}
else
{
/* allocate memory */
ap = (float *) malloc((sizeof (float))*nm1*mm1);
bp = (float *) malloc((sizeof (float))*mm1);
/* reduce a and b into ap and bp eliminating variable t and constraint i */
jj=-1;
for (j=0;j<mm;j++)
if (j != i)
{
jj=jj+1;
aj=a+j;
ajt=aj+tmm;
*(bp+jj)= *(b+j) - *(bi) * *(ajt) / pivot;
apjj=ap+jj;
kk=-1;
for (k=0;k<nn;k++)
if (k != t)
{
kk=kk+1;
kmm=k*mm_max;
*(apjj+kk*mm1)= *(aj+kmm)- *(ajt) * *(ai+kmm)/ pivot;
}
}
/* add contribution to volume from volume calculated in smaller dimension */
vol=cvolume(ap,bp, &mm1, &nm1, &mm1, &nm1);
if(vol == -1.0)
{
*dvdb = 0.;
return(-1.0);
}
v=v+ *(bi)/amax*vol/nn;
free(ap);
free(bp);
/* calculate derivatives for first constraint only */
/* calculate volume and derivative
with respect to b_0 */ if (i == 0) *dvdb = vol/amax;
/* calculate derivative with respect
to b_0 but not volume */
if (*opt == 2) return (0.);
}
}
return(v);
}
/*--------------------------------------------------------------
ROUTINE: volumeb
A dummy routine used to call cvolumeb from a fortran routine
----------------------------------------------------------------*/
void volumeb_();
void volumeb();
void volumeb_(a,b,m,n,mmax,nmax,opt,volume,dvdb)
int *n, *m, *mmax, *nmax, *opt;
float *a, *b;
float *volume;
float *dvdb;
{
volumeb(a,b,m,n,mmax,nmax,opt,volume,dvdb);
}
void volumeb (a,b,m,n,mmax,nmax,opt,volume,dvdb)
int *n, *m, *mmax, *nmax, *opt;
float *a, *b;
float *volume;
float *dvdb;
{
if(*opt == 1) /* calculate volume and derivative with
respect to b(0) */
{
*volume=cvolumeb(a,b,m,n,mmax,nmax,opt,dvdb);
}
else if(*opt == 2) /* calculate derivative with respect to b(0)
but not volume */
{
*volume=cvolumeb(a,b,m,n,mmax,nmax,opt,dvdb);
}
else
{
printf(" Warning: routine volumeb called with an invalid option parameter\n");
printf(" valid options are 1 and 2, option given = %d \n",*opt);
}
}
/*--------------------------------------------------------------
ROUTINE: cvolumebj
This routine only calculates the j-th outer loop of the recursive
routine cvolumeb. It returns the partial volume for projection
onto the j-th constraint and the derivative of the total volume
with respect to parameter b(j) using the simple formula of
Lasserre(1983). (See cvolumeb for more details.)
The reason for this routine is so that the derivatives with respect to parameter b(j) can be calculated and passed back
to a fortran routine without creating and an extra array of size
m (because only one derivative is calculated per call). In
this way it also avoids recalculating the total volume m times
(which would be the case if we used a simple variation of routine
cvolumeb).
To calculate the total volume this routine must be called m times
and the partial volume summed
Constraint j is determined by the value of *con.
Calls are made to routine cvolume.
Malcolm Sambridge, March 1995.
--------------------------------------------------------------*/
float cvolumebj (a,b,m,n,mmax,nmax,con,dvdb)
int *n, *m, *mmax, *nmax, con;
float *a, *b;
float *dvdb;
{
float v,amax,pivot;
int i,j,t,k,l;
int jj,kk;
float *ai, *aj, *ajt, *apjj, *bi;
int kmm,tmm;
int nn,mm,nn_max,mm_max;
int nm1,mm1;
int lmin,lmax;
float *ap, *bp;
int firstmin,firstmax;
float bmin,bmax,bb;
float vol;
nn= *n;
mm= *m;
nn_max= *nmax;
mm_max= *mmax;
/* one-dimensional case (full reduction) */
if (nn == 1)
{
firstmin=0;
firstmax=0;
lmax=0;
lmin=0;
for (l=0;l<mm;l++)
{
if ( *(a+l) > 0.)
{
bb= *(b+l)/ *(a+l);
if (firstmin==0) {firstmin=1;bmin=bb;lmin=l;}
else if (bb<bmin) {bmin=bb;lmin=l;}
}
else if ( *(a+l) < 0.)
{
bb= *(b+l)/ *(a+l);
if (firstmax==0) {firstmax=1;bmax=bb;lmax=l;}
else if (bb>bmax) {bmax=bb;lmax=l;}
}
else if ( *(b+l) < 0.)
{
/* Constraints are inconsistent.
Set volume and derivative to zero. */
printf("Inconsistent constraints found after reduction to n = 1 \n");
*dvdb = 0.;
return(0.);
}
}
v=0.;
*dvdb = 0.;
if (firstmin*firstmax == 1) v=bmin-bmax;
else
{
printf("Volume is unbounded; volume returned as -1\n derivatives returned as zero\n");
return(-1.0);
}
if (v<0.) {v=0.;*dvdb=0.;}
else if (v>0. && lmin == con) *dvdb = 1. / *(a+lmin);
else if (v>0. && lmax == con) *dvdb = -1. / *(a+lmax);
return(v);
}
nm1=nn-1;
mm1=mm-1;
v=0.;
/* perform main loop over constraints */
/* for (i=0;i<mm;i++) */
i = con;
ai=a+i;
bi=b+i;
/* find largest pivot */
amax=0.;
for (j=0;j<nn;j++)
if (fabs( *(ai+j*mm_max)) >= amax) {amax= fabs( *(ai+j*mm_max)); t=j;}
tmm=t*mm_max;
pivot=*(ai+tmm);
/* finds contribution to v from this pivot (if not nil) */
if (amax == 0.)
{
/* Constraint is inconsistent */
if ( *(bi) < 0.)
{
printf("Constraint %d is inconsistent\n",i+1);
*dvdb = 0.;
return(0.);
}
/* otherwise constraint is redundant */
printf("Constraint %d is redundant\n",i+1);
*dvdb = 0.;
}
else
{
/* allocate memory */
ap = (float *) malloc((sizeof (float))*nm1*mm1);
bp = (float *) malloc((sizeof (float))*mm1);
/* reduce a and b into ap and bp eliminating variable t and constraint i */
jj=-1;
for (j=0;j<mm;j++)
if (j != i)
{
jj=jj+1;
aj=a+j;
ajt=aj+tmm;
*(bp+jj)= *(b+j) - *(bi) * *(ajt) / pivot;
apjj=ap+jj;
kk=-1;
for (k=0;k<nn;k++)
if (k != t)
{
kk=kk+1;
kmm=k*mm_max;
*(apjj+kk*mm1)= *(aj+kmm)- *(ajt) * *(ai+kmm)/ pivot;
}
}
/* calculate partial volume */
vol=cvolume(ap,bp, &mm1, &nm1, &mm1, &nm1);
if(vol == -1.0)
{
*dvdb = 0.;
return(-1.0);
}
v=*(bi)/amax*vol/nn;
/* calculate derivative of total volume
with respect to current constraint */
*dvdb = vol/amax;
free(ap);
free(bp);
}
return(v);
}
/*--------------------------------------------------------------
ROUTINE: volumebj
A dummy routine used to call cvolumebj from a fortran routine
----------------------------------------------------------------*/
void volumebj_();
void volumebj();
void volumebj_ (a,b,m,n,mmax,nmax,con,volume,dvdb)
int *n, *m, *mmax, *nmax,*con;
float *a, *b;
float *volume;
float *dvdb;
{
volumebj(a,b,m,n,mmax,nmax,con,volume,dvdb);
}
void volumebj (a,b,m,n,mmax,nmax,con,volume,dvdb)
int *n, *m, *mmax, *nmax,*con;
float *a, *b;
float *volume;
float *dvdb;
{
int tcon;
tcon = *con - 1;
/* calculate partial volume and
derivative with respect to b(con) */
*volume=cvolumebj(a,b,m,n,mmax,nmax,tcon,dvdb);
}
/*--------------------------------------------------------------
ROUTINE: cdvda
This routine calculates the derivative with respect to a(0,tdim)
of the `volume' in dimension n of the region bounded by a set
of m linear inequality constraints of the form A x <= b, where
a has m rows and n columns and is given by a(m,n), b is the
n-vector and is contained in b(m). The derivative expression
is recursive and derived from the formula of Lasserre (1983).
Redundant constraints are allowed and a warning is issued if any are
encountered. If the inequality constraints are inconsistent then the
derivative is returned as zero. If any constraint is orthogonal to
the component a(0,idim) then the reduction can only take place onto
variable idim. A special case is used to handle this which involves
no further recursive calls.
If the original polyhedron is unbounded then a warning is issued
and the derivative is return as zero.
Note: This code takes advantage of the fact that during recursive calls
constraint 0 does not change its position in the list of remaining
constraints (if it has not been eliminated), i.e. it is always the
first constraint. This would not be the case if the algorithm were
adapated to deal with other constraints, i.e. evaluate dvda_i,j
where i .ne. 0.
Calls itself cvolumeb, and cvolumebj.
Malcolm Sambridge, March 1995.
--------------------------------------------------------------*/
float cdvda (a,b,m,n,mmax,nmax,tdim,temp,jval,code)
int *n, *m, *nmax, *mmax, *tdim, *jval, *code;
float *a, *b, *temp;
{
float v,amax,pivot;
int i,j,t,k,l;
int jj,kk;
float *ai, *aj, *ajt, *apjj, *bi;
int kmm,tmm;
int nn,mm,nn_max,mm_max,ttdim;
int nm1,mm1;
int lmin,lmax, jjval, kval;
float *ap, *bp, *ttemp;
int firstmin,firstmax;
float bmin,bmax,bb;
float deriv, junk, vol, dvdb, dbda;
int special, opt;
nn= *n;
mm= *m;
nn_max= *nmax;
mm_max= *mmax;
/* one-dimensional case (full reduction) */
*code = 0;
if (nn == 1)
{
firstmin=0;
firstmax=0;
lmax=0;
lmin=0;
for (l=0;l<mm;l++)
{
if ( *(a+l) > 0.)
{
bb= *(b+l)/ *(a+l);
if (firstmin==0) {firstmin=1;bmin=bb;lmin=l;}
else if (bb<bmin) {bmin=bb;lmin=l;}
}
else if ( *(a+l) < 0.)
{
bb= *(b+l)/ *(a+l);
if (firstmax==0) {firstmax=1;bmax=bb;lmax=l;}
else if (bb>bmax) {bmax=bb;lmax=l;}
}
else if ( *(b+l) < 0.)
{
/* Constraint is inconsistent.
Set derivative to zero. */
printf("Inconsistent constraints found after reduction to n = 1 \n");
*code = 1;
return(0.);
}
}
v=0.;
if (firstmin*firstmax == 1) v=bmin-bmax;
else
{
printf("Volume is unbounded; derivative returned is zero\n");
*code = -1;
return(0.);
}
if (v<0.) return(0.);
if(*jval == 1) /* Constraint 0 has not yet been encountered */
{
if (lmin == 0) deriv = -bmin/ *a;
else if (lmax == 0) deriv = bmax/ *a;
else deriv = 0.;
return(deriv);
}
else if(*jval == 0) /* Constraint 0 has already been encountered */
{
deriv = ( *(temp+lmax) * (bmax/ *(a+lmax)) ) -
( *(temp+lmin) * (bmin/ *(a+lmin)) );
return(deriv);
}
}
nm1=nn-1;
mm1=mm-1;
v=0.;
/* perform main loop over constraints */
for (i=0;i<mm;i++)
{
ai=a+i;
bi=b+i;
ttdim = *tdim;
special = 0;
/* find largest pivot */
amax=0.;
t = 0;
for (j=0;j<nn;j++)
if (fabs( *(ai+j*mm_max)) >= amax && j != ttdim)
{amax= fabs( *(ai+j*mm_max)); t=j;}
/* finds contribution to v from
this pivot (if not nil) */
if (amax == 0.)
{
if(*(ai + ttdim * mm_max) == 0.0)
{
/* Constraint is inconsistent */
if ( *(bi) < 0.)
{
printf("Constraint %d is inconsistent\n",i+1);
*code = 1;
return(0.);
}
/* otherwise constraint is redundant */
printf("Constraint %d is redundant\n",i+1);
}
else
{
/* if projection can only be peformed
on dimension tdim then activate
special case */
special = 1;
t = ttdim;
amax = fabs(*(ai+t * mm_max));
}
}
tmm=t*mm_max;
pivot=*(ai+tmm);
if(t < ttdim) ttdim = ttdim -1;
if (amax != 0)
{
/* determine if constraint 0 has been encountered */
kval = 0;
if ( i == 0 && *jval == 1)
{
/* This is the first encounter of
constraint 0 on this path so we
allocate memory and store parameters
to be used when n = 1 */
if (special == 0)
{
ttemp = (float *) malloc((sizeof (float))*mm1);
for (j=0;j<mm1;j++) *(ttemp+j) = - *(a+j+1+tmm)/pivot;
kval = 1;
}
jjval = 0;
}
else if (*jval == 0) /* Constraint 0 has already been
encountered */ {
jjval = 0;
/* perform recursive update of component
derivative array temp. This eliminates
row i and copies into a new vector */
if(special == 0)
{
ttemp = (float *) malloc((sizeof (float))*mm1);
for (j=0;j<i;j++) *(ttemp+j) = *(temp+j)
-(*(temp+i) * *(a+j+tmm)/pivot);
for (j=i;j<mm1;j++) *(ttemp+j) = *(temp+j+1)
-(*(temp+i) * *(a+j+1+tmm)/pivot);
}
}
else
{ /* Constraint 0 has not yet been
encountered */
jjval = 1;
}
/* allocate memory */
ap = (float *) malloc((sizeof (float))*nm1*mm1);
bp = (float *) malloc((sizeof (float))*mm1);
/* reduce a and b into ap and bp eliminating variable t and constraint i */
jj=-1;
for (j=0;j<mm;j++)
if (j != i)
{
jj=jj+1;
aj=a+j;
ajt=aj+tmm;
*(bp+jj)= *(b+j) - *(bi) * *(ajt) / pivot;
apjj=ap+jj;
kk=-1;
for (k=0;k<nn;k++)
if (k != t)
{
kk=kk+1;
kmm=k*mm_max;
*(apjj+kk*mm1)= *(aj+kmm)- *(ajt) * *(ai+kmm)/ pivot;
}
}
/* add contribution to derivative from that calculated in smaller dimension */
/* Normal case method */
if(special == 0)
{
deriv=cdvda(ap,bp, &mm1, &nm1, &mm1, &nm1, &ttdim, ttemp, &jjval,code);
v=v+ *(bi)/amax*deriv/nn;
if (kval == 1 || *jval == 0) free(ttemp);
if(*code != 0)return (0.);
}
else /* Use special case method */
{
if( *jval == 1)
{
if(i == 0) /* This is constraint 0 */
{
deriv = 0.;
vol = 0.;
dvdb = 0.;
for (j=1;j<mm;j++)
{ k = j - 1;
junk=cvolumebj(ap,bp,&mm1,&nm1,&mm1,&nm1,k,&dvdb);
if(junk == -1.)
{
*code = -1;
return(0.);
}
deriv = deriv + dvdb * *(a + j + tmm) ;
vol = vol + junk;
}
if(nm1 == 1)vol = junk;
deriv = *(bi) * deriv/pivot;
deriv = (deriv - vol) /pivot;
v=v+ *(bi)/amax*deriv/nn;
}
else /* Constraint 0 not yet encountered */
{
opt = 2;
junk=cvolumeb(ap,bp,&mm1,&nm1,&mm1,&nm1,&opt,&deriv);
if(junk == -1.)
{
*code = -1;
return(0.);
}
deriv= -(deriv * *(bi)/pivot);
v=v+ *(bi)/amax*deriv/nn;
}
}
else if( *jval == 0) /* Constraint 0 already encountered */
{
vol = 0.;
deriv = 0.;
for (j=0;j<i;j++)
{
junk=cvolumebj(ap,bp,&mm1,&nm1,&mm1,&nm1,j,&dvdb);
if(junk == -1.)
{
*code = -1;
return(0.);
}
dbda = (*(a+j+tmm) * *(temp+i)/pivot) - *(temp+j);
deriv = deriv + dvdb*dbda;
vol = vol + junk;
}
for (j=i+1;j<mm;j++)
{
k = j - 1;
junk=cvolumebj(ap,bp,&mm1,&nm1,&mm1,&nm1,k,&dvdb);
if(junk == -1.)
{
*code = -1;
return(0.);
}
dbda = (*(a+j+tmm) * *(temp+i)/pivot) - *(temp+j);
deriv = deriv + dvdb*dbda;
vol = vol + junk;
}
if(nm1 == 1)vol = junk;
deriv = (*(bi) * deriv - *(temp+i)*vol)/pivot;
v=v+ *(bi)/amax*deriv/nn;
}
}
free(ap);
free(bp);
}
}
return(v);
}
/*--------------------------------------------------------------
ROUTINE: dvda
A dummy routine used to call cdvda from a fortran routine
----------------------------------------------------------------*/
/* Not used...won't compile...function name and variable name (were) identical */
void dvda_();
void dvda();
void dvda_ (a,b,m,n,mmax,nmax,idim,result,code)
int *n, *m, *mmax, *nmax, *idim, *code;
float *a, *b;
float *result;
{
dvda(a,b,m,n,mmax,nmax,idim,result,code);
}
void dvda (a,b,m,n,mmax,nmax,idim,result,code)
int *n, *m, *mmax, *nmax, *idim, *code;
float *a, *b;
float *result;
{
int j;
int jval, tdim;
float *temp;
jval = 1;
tdim = *idim - 1;
*result=cdvda(a,b,m,n,mmax,nmax,&tdim,temp,&jval,code);
free(temp);
}
/*--------------------------------------------------------------
ROUTINE: cvolumef
This routine calculates the `volume' in dimension n of the region
bounded by a set of m linear inequality constraints of the form
A x <= b, where a has m rows and n columns and is given by a(m,n),
b is the n-vector and is contained in b(m). The recursive formula
of Lasserre (1983) is used. Redundant constraints are allowed
If the inequality constraints are inconsistent then the volume
is returned as zero. If the original polyhedron is unbounded
then a warning is issued and the volume is return as -1.0
(even though it is undefined).
This is a faster version which does not continue down the
recursive tree if if encounters b(j) = 0 for any j. This
restricts its use to cases where the origin does not
pass through any inconsistent constraints. Also redundant
constraints that pass through the origin will not be detected.
Jean Braun and Malcolm Sambridge, Jan 1995.
Lasserre, J. B., 1983. "An analytical expression and an Algorithm
for the Volume of a Convex Polyhedron in R^n.", JOTA, vol. 39,
no. 3, 363-377.
--------------------------------------------------------------*/
float cvolumef (a,b,m,n,mmax,nmax)
int *n, *m, *mmax, *nmax;
float *a, *b;
{
float v,amax,pivot;
int i,j,t,k,l;
int jj,kk;
float *ai, *aj, *ajt, *apjj, *bi;
int kmm,tmm;
int nn,mm,nn_max,mm_max;
int nm1,mm1;
float *ap, *bp;
int firstmin,firstmax;
float bmin,bmax,bb;
float partialv;
nn= *n;
mm= *m;
nn_max= *nmax;
mm_max= *mmax;
/* one-dimensional case (full reduction) */
if (nn == 1)
{
firstmin=0;
firstmax=0;
for (l=0;l<mm;l++)
{
if ( *(a+l) > 0.)
{
bb= *(b+l)/ *(a+l);
if (firstmin==0) {firstmin=1;bmin=bb;}
else if (bb<bmin) bmin=bb;
}
else if ( *(a+l) < 0.)
{
bb= *(b+l)/ *(a+l);
if (firstmax==0) {firstmax=1;bmax=bb;}
else if (bb>bmax) bmax=bb;
}
else if ( *(b+l) < 0.) /* Constraints are inconsistent */
{
printf("Inconsistent constraints found after reduction to n = 1 \n");
return(0.);
}
}
v=0.;
if (firstmin*firstmax == 1) v=bmin-bmax;
else
{
printf("Volume is unbounded at 1-D; volume returned as -1\n");
return(-1.0);
}
if (v<0.) {v=0.;}
return(v);
}
nm1=nn-1;
mm1=mm-1;
v=0.;
for (i=0;i<mm;i++)
{
ai=a+i;
bi=b+i;
/* find largest pivot */
amax=0.;
for (j=0;j<nn;j++)
if (fabs( *(ai+j*mm_max)) >= amax) {amax= fabs( *(ai+j*mm_max)); t=j;}
tmm=t*mm_max;
pivot=*(ai+tmm);
/* finds contribution to v from this pivot (if not nil) */
if (amax == 0.)
{
/* Constraint is inconsistent */
if ( *(bi) < 0.)
{ printf("Constraint %d is inconsistent\n",i+1); return(0.); }
/* otherwise constraint is redundant */
printf("Constraint %d is redundant\n",i+1);
}
else
{
/* allocate memory */
ap = (float *) malloc((sizeof (float))*nm1*mm1);
bp = (float *) malloc((sizeof (float))*mm1);
/* reduce a and b into ap and bp eliminating variable t and constraint i */
jj=-1;
for (j=0;j<mm;j++)
if (j != i)
{
jj=jj+1;
aj=a+j;
ajt=aj+tmm;
*(bp+jj)= *(b+j) - *(bi) * *(ajt) / pivot;
apjj=ap+jj;
kk=-1;
for (k=0;k<nn;k++)
if (k != t)
{
kk=kk+1;
kmm=k*mm_max;
*(apjj+kk*mm1)= *(aj+kmm)- *(ajt) * *(ai+kmm)/ pivot;
}
}
/* add contribution to volume from volume calculated in smaller dimension */
partialv = 0.;
if (*(bi) != 0.) partialv=cvolumef(ap,bp, &mm1, &nm1, &mm1, &nm1);
if(partialv == -1.0)return(-1.0);
v=v+ *(bi)/amax*partialv/nn;
free(ap);
free(bp);
}
}
return(v);
}/*--------------------------------------------------------------
ROUTINE: volume
A dummy routine used to call cvolume from a fortran routine
----------------------------------------------------------------*/
void volumef_();
void volumef();
void volumef_ (a,b,m,n,mmax,nmax,volume)
int *n, *m, *mmax, *nmax;
float *a, *b;
float *volume;
{
volumef(a,b,m,n,mmax,nmax,volume);
}
void volumef (a,b,m,n,mmax,nmax,volume)
int *n, *m, *mmax, *nmax;
float *a, *b;
float *volume;
{
*volume=cvolumef(a,b,m,n,mmax,nmax);
}
/*--------------------------------------------------------------
ROUTINE: cdvdaf
This routine calculates the derivative with respect to a(0,tdim)
of the `volume' in dimension n of the region bounded by a set
of m linear inequality constraints of the form A x <= b, where
a has m rows and n columns and is given by a(m,n), b is the
n-vector and is contained in b(m). The derivative expression
is recursive and derived from the formula of Lasserre (1983).
Redundant constraints are allowed and a warning is issued if any are
encountered. If the inequality constraints are inconsistent then the
derivative is returned as zero. If any constraint is orthogonal to
the component a(0,idim) then the reduction can only take place onto
variable idim. A special case is used to handle this which involves
no further recursive calls.
If the original polyhedron is unbounded then a warning is issued
and the derivative is return as zero.
Note: This code takes advantage of the fact that during recursive calls
constraint 0 does not change its position in the list of remaining
constraints (if it has not been eliminated), i.e. it is always the
first constraint. This would not be the case if the algorithm were
adapated to deal with other constraints, i.e. evaluate dvda_i,j
where i .ne. 0.
This is a faster version which does not continue down the
recursive tree if if encounters b(j) = 0 for any j. This
restricts its use to cases where the origin does not
pass through any inconsistent constraints. Also redundant
constraints that pass through the origin will not be detected.
Calls itself cvolumeb, and cvolumebj.
Malcolm Sambridge, March 1995.
--------------------------------------------------------------*/
float cdvdaf (a,b,m,n,mmax,nmax,tdim,temp,jval,code)
int *n, *m, *nmax, *mmax, *tdim, *jval, *code;
float *a, *b, *temp;
{
float v,amax,pivot;
int i,j,t,k,l;
int jj,kk;
float *ai, *aj, *ajt, *apjj, *bi;
int kmm,tmm;
int nn,mm,nn_max,mm_max,ttdim;
int nm1,mm1;
int lmin,lmax, jjval, kval;
float *ap, *bp, *ttemp;
int firstmin,firstmax;
float bmin,bmax,bb;
float deriv, junk, vol, dvdb, dbda;
int special, opt;
nn= *n;
mm= *m;
nn_max= *nmax;
mm_max= *mmax;
/* one-dimensional case (full reduction) */
*code = 0;
if (nn == 1)
{
firstmin=0;
firstmax=0;
lmax=0;
lmin=0;
for (l=0;l<mm;l++)
{
if ( *(a+l) > 0.)
{
bb= *(b+l)/ *(a+l);
if (firstmin==0) {firstmin=1;bmin=bb;lmin=l;}
else if (bb<bmin) {bmin=bb;lmin=l;}
}
else if ( *(a+l) < 0.)
{
bb= *(b+l)/ *(a+l);
if (firstmax==0) {firstmax=1;bmax=bb;lmax=l;}
else if (bb>bmax) {bmax=bb;lmax=l;}
}
else if ( *(b+l) < 0.)
{
/* Constraint is inconsistent.
Set derivative to zero. */
printf("Inconsistent constraints found after reduction to n = 1 \n");
*code = 1;
return(0.);
}
}
v=0.;
if (firstmin*firstmax == 1) v=bmin-bmax;
else
{
printf("Volume is unbounded; derivative returned is zero\n");
*code = -1;
return(0.);
}
if (v<0.) return(0.);
if(*jval == 1) /* Constraint 0 has not yet been encountered */
{
if (lmin == 0) deriv = -bmin/ *a;
else if (lmax == 0) deriv = bmax/ *a;
else deriv = 0.;
return(deriv);
}
else if(*jval == 0) /* Constraint 0 has already been encountered */
{
deriv = ( *(temp+lmax) * (bmax/ *(a+lmax)) ) -
( *(temp+lmin) * (bmin/ *(a+lmin)) );
return(deriv);
}
}
nm1=nn-1;
mm1=mm-1;
v=0.;
/* perform main loop over constraints */
for (i=0;i<mm;i++)
{
ai=a+i;
bi=b+i;
ttdim = *tdim;
special = 0;
/* find largest pivot */
amax=0.;
t = 0;
for (j=0;j<nn;j++)
if (fabs( *(ai+j*mm_max)) >= amax && j != ttdim)
{amax= fabs( *(ai+j*mm_max)); t=j;}
/* finds contribution to v from
this pivot (if not nil) */
if (amax == 0.)
{
if(*(ai + ttdim * mm_max) == 0.0)
{
/* Constraint is inconsistent */
if ( *(bi) < 0.)
{
printf("Constraint %d is inconsistent\n",i+1);
*code = 1;
return(0.);
}
/* otherwise constraint is redundant */
printf("Constraint %d is redundant\n",i+1);
}
else
{
/* if projection can only be peformed
on dimension tdim then activate
special case */
special = 1;
t = ttdim;
amax = fabs(*(ai+t * mm_max));
}
}
tmm=t*mm_max; pivot=*(ai+tmm);
if(t < ttdim) ttdim = ttdim -1;
if (amax != 0)
{
/* determine if constraint 0 has been encountered */
kval = 0;
if ( i == 0 && *jval == 1)
{
/* This is the first encounter of
constraint 0 on this path so we
allocate memory and store parameters
to be used when n = 1 */
if (special == 0)
{
ttemp = (float *) malloc((sizeof (float))*mm1);
for (j=0;j<mm1;j++) *(ttemp+j) = - *(a+j+1+tmm)/pivot;
kval = 1;
}
jjval = 0;
}
else if (*jval == 0) /* Constraint 0 has already been
encountered */
{
jjval = 0;
/* perform recursive update of component
derivative array temp. This eliminates
row i and copies into a new vector */
if(special == 0)
{
ttemp = (float *) malloc((sizeof (float))*mm1);
for (j=0;j<i;j++) *(ttemp+j) = *(temp+j)
-(*(temp+i) * *(a+j+tmm)/pivot);
for (j=i;j<mm1;j++) *(ttemp+j) = *(temp+j+1)
-(*(temp+i) * *(a+j+1+tmm)/pivot);
}
}
else
{ /* Constraint 0 has not yet been
encountered */
jjval = 1;
}
/* allocate memory */
ap = (float *) malloc((sizeof (float))*nm1*mm1);
bp = (float *) malloc((sizeof (float))*mm1);
/* reduce a and b into ap and bp eliminating variable t and constraint i */
jj=-1;
for (j=0;j<mm;j++)
if (j != i)
{
jj=jj+1;
aj=a+j;
ajt=aj+tmm;
*(bp+jj)= *(b+j) - *(bi) * *(ajt) / pivot;
apjj=ap+jj;
kk=-1;
for (k=0;k<nn;k++)
if (k != t)
{ kk=kk+1;
kmm=k*mm_max;
*(apjj+kk*mm1)= *(aj+kmm)- *(ajt) * *(ai+kmm)/ pivot;
}
}
/* add contribution to derivative from that calculated in smaller dimension */
/* Normal case method */
if(special == 0)
{
deriv = 0.;
if (*(bi) != 0.)
{deriv=cdvdaf(ap,bp, &mm1, &nm1, &mm1, &nm1, &ttdim, ttemp, &jjval,code);}
v=v+ *(bi)/amax*deriv/nn;
if (kval == 1 || *jval == 0) free(ttemp);
if(*code != 0)return (0.);
}
else /* Use special case method */
{
if( *jval == 1)
{
if(i == 0) /* This is constraint 0 */
{
deriv = 0.;
vol = 0.;
dvdb = 0.;
for (j=1;j<mm;j++)
{
k = j - 1;
junk=cvolumebj(ap,bp,&mm1,&nm1,&mm1,&nm1,k,&dvdb);
if(junk == -1.)
{
*code = -1;
return(0.);
}
deriv = deriv + dvdb * *(a + j + tmm) ;
vol = vol + junk;
}
if(nm1 == 1)vol = junk;
deriv = *(bi) * deriv/pivot;
deriv = (deriv - vol) /pivot;
v=v+ *(bi)/amax*deriv/nn;
}
else /* Constraint 0 not yet encountered */
{
opt = 2;
junk=cvolumeb(ap,bp,&mm1,&nm1,&mm1,&nm1,&opt,&deriv);
if(junk == -1.)
{
*code = -1;
return(0.);
}
deriv= -(deriv * *(bi)/pivot);
v=v+ *(bi)/amax*deriv/nn;
}
}
else if( *jval == 0) /* Constraint 0 already encountered */
{
vol = 0.;
deriv = 0.;
for (j=0;j<i;j++)
{
junk=cvolumebj(ap,bp,&mm1,&nm1,&mm1,&nm1,j,&dvdb);
if(junk == -1.)
{
*code = -1;
return(0.); }
dbda = (*(a+j+tmm) * *(temp+i)/pivot) - *(temp+j);
deriv = deriv + dvdb*dbda;
vol = vol + junk;
}
for (j=i+1;j<mm;j++)
{
k = j - 1;
junk=cvolumebj(ap,bp,&mm1,&nm1,&mm1,&nm1,k,&dvdb);
if(junk == -1.)
{
*code = -1;
return(0.);
}
dbda = (*(a+j+tmm) * *(temp+i)/pivot) - *(temp+j);
deriv = deriv + dvdb*dbda;
vol = vol + junk;
}
if(nm1 == 1)vol = junk;
deriv = (*(bi) * deriv - *(temp+i)*vol)/pivot;
v=v+ *(bi)/amax*deriv/nn;
}
}
free(ap);
free(bp);
}
}
return(v);
}
/*--------------------------------------------------------------
ROUTINE: dvda
A dummy routine used to call cdvda from a fortran routine
----------------------------------------------------------------*/
void dvdaf_();
void dvdaf();
void dvdaf_ (a,b,m,n,mmax,nmax,idim,dvda,code)
int *n, *m, *mmax, *nmax, *idim, *code;
float *a, *b;
float *dvda;
{
dvdaf(a,b,m,n,mmax,nmax,idim,dvda,code);
}
void dvdaf (a,b,m,n,mmax,nmax,idim,dvda,code)
int *n, *m, *mmax, *nmax, *idim, *code;
float *a, *b;
float *dvda;
{
int j;
int jval, tdim;
float *temp;
jval = 1;
tdim = *idim - 1;
*dvda=cdvdaf(a,b,m,n,mmax,nmax,&tdim,temp,&jval,code);
free(temp);
}
/*--------------------------------------------------------------
ROUTINE: cdvda_debug
This routine calculates the derivative with respect to a(0,tdim)
of the `volume' in dimension n of the region bounded by a set
of m linear inequality constraints of the form A x <= b, where
a has m rows and n columns and is given by a(m,n), b is the
n-vector and is contained in b(m). The derivative expression
is recursive and derived from the formula of Lasserre (1983).
Redundant constraints are allowed and a warning is issued if any are
encountered. If the inequality constraints are inconsistent then the
derivative is returned as zero. If any constraint is orthogonal to
the component a(0,idim) then the reduction can only take place onto
variable idim. A special case is used to handle this which involves
no further recursive calls.
If the original polyhedron is unbounded then a warning is issued
and the derivative is return as zero.
Note: This code takes advantage of the fact that during recursive calls
constraint 0 does not change its position in the list of remaining
constraints (if it has not been eliminated), i.e. it is always the
first constraint. This would not be the case if the algorithm were
adapated to deal with other constraints, i.e. evaluate dvda_i,j
where i .ne. 0.
Calls itself cvolumeb, and cvolumebj.
This is a version of cdvda used for debugging and contains extra
code and print statements.
Malcolm Sambridge, March 1995.
--------------------------------------------------------------*/
float cdvda_debug (a,b,m,n,tdim,temp,jval,code)
int *n, *m, *tdim, *jval, *code;
float *a, *b, *temp;
{
float v,amax,pivot;
int i,j,t,k,l;
int jj,kk;
float *ai, *aj, *ajt, *apjj, *bi;
int kmm,tmm;
int nn,mm, ttdim;
int nm1,mm1;
int lmin,lmax, jjval, kval;
float *ap, *bp, *ttemp;
int firstmin,firstmax;
float bmin,bmax,bb;
float deriv, junk, vol, dvdb, dbda;
int special, opt;
float volp, volm, da, FD, FDa;
nn= *n;
mm= *m;
/* write out a and b (debug) */
printf("\n"); printf(" Current matrix A and vector b \n");
for (j=0;j<mm;j++)
{
aj=a+j;
for (k=0;k<nn;k++)
{ajt=aj+k * mm; printf(" a %d %d = %f",j,k,*(ajt));}
printf(" : b = %f \n", *(b+j));
}
/* one-dimensional case (full reduction) */
*code = 0;
if (nn == 1)
{
printf("One dimension left \n");
firstmin=0;
firstmax=0;
lmax=0;
lmin=0;
for (l=0;l<mm;l++)
{
if ( *(a+l) > 0.)
{
bb= *(b+l)/ *(a+l);
if (firstmin==0) {firstmin=1;bmin=bb;lmin=l;}
else if (bb<bmin) {bmin=bb;lmin=l;}
}
else if ( *(a+l) < 0.)
{
bb= *(b+l)/ *(a+l);
if (firstmax==0) {firstmax=1;bmax=bb;lmax=l;}
else if (bb>bmax) {bmax=bb;lmax=l;}
}
else if ( *(b+l) < 0.)
{
/* Constraint is inconsistent.
Set derivative to zero. */
printf("Inconsistent constraints found after reduction to n = 1 \n");
*code = 1;
return(0.);
}
}
v=0.;
if (firstmin*firstmax == 1) v=bmin-bmax;
else
{
printf("Volume is unbounded; derivative returned is zero\n");
*code = -1;
return(0.);
}
if (v<0.) return(0.);
if(*jval == 1) /* Constraint 0 has not yet been encountered */
{
if (lmin == 0) deriv = -bmin/ *a;
else if (lmax == 0) deriv = bmax/ *a;
else deriv = 0.;
printf(" No projection onto constraint 0 \n");
printf(" lmin = %d lmax = %d bmin = %f bmax = %f alpha = %f \n",
lmin,lmax,bmin,bmax,*a);
printf(" deriv contribution = %f \n",deriv);
return(deriv);
}
else if(*jval == 0) /* Constraint 0 has already been encountered */
{
deriv = ( *(temp+lmax) * (bmax/ *(a+lmax)) ) -
( *(temp+lmin) * (bmin/ *(a+lmin)) ); printf(" Projection onto constraint 0 has occurred \n");
printf(" lmin = %d bmin = %f alpha = %f temp = %f\n",
lmin,bmin,*(a+lmin),*(temp+lmin));
printf(" lmax = %d bmax = %f alpha = %f temp = %f\n",
lmax,bmax,*(a+lmax),*(temp+lmax));
printf(" deriv contribution = %f \n",deriv);
return(deriv);
}
}
nm1=nn-1;
mm1=mm-1;
v=0.;
printf(" About to start main loop n = %d m = %d \n",nn,mm);
/* perform main loop over constraints */
for (i=0;i<mm;i++)
{
ai=a+i;
bi=b+i;
ttdim = *tdim;
special = 0;
/* find largest pivot */
amax=0.;
t = 0;
for (j=0;j<nn;j++)
if (fabs( *(ai+j*mm)) >= amax && j != ttdim)
{amax= fabs( *(ai+j*mm)); t=j;}
/* finds contribution to v from
this pivot (if not nil) */
if (amax == 0.)
{
if(*(ai + ttdim * mm) == 0.0)
{
/* Constraint is inconsistent */
if ( *(bi) < 0.)
{
printf("Constraint %d is inconsistent\n",i+1);
*code = 1;
return(0.);
}
/* otherwise constraint is redundant */
printf("Constraint %d is redundant\n",i+1);
}
else
{
/* if projection can only be peformed
on dimension tdim then activate
special case */
printf("Constraint %d can only be projected onto dimension %d\n",i+1,*tdim);
printf("using special case method \n");
special = 1;
t = ttdim;
amax = fabs(*(ai+t * mm));
}
}
tmm=t*mm;
pivot=*(ai+tmm);
printf("\n Projection onto constraint %d variable %d k = %d\n",i,t,ttdim);
if(t < ttdim) ttdim = ttdim -1;
if (amax != 0)
{
/* determine if constraint 0 has been encountered */
kval = 0;
if ( i == 0 && *jval == 1)
{
/* This is the first encounter of
constraint 0 on this path so we
allocate memory and store parameters
to be used when n = 1 */
if (special == 0)
{
ttemp = (float *) malloc((sizeof (float))*mm);
for (j=0;j<mm1;j++) *(ttemp+j) = - *(a+j+1+tmm)/pivot;
kval = 1;
printf("\n Generating temp n = %d m = %d i = %d \n",nn,mm,i);
printf(" ttemp : \n");
for (j=0;j<mm1;j++) printf(" %f ",*(ttemp+j));
printf("\n");
}
jjval = 0;
}
else if (*jval == 0) /* Constraint 0 has already been
encountered */
{
jjval = 0;
/* perform recursive update of component
derivative array temp. This eliminates
row i and copies into a new vector */
if(special == 0)
{
ttemp = (float *) malloc((sizeof (float))*mm1);
for (j=0;j<i;j++) *(ttemp+j) = *(temp+j)
-(*(temp+i) * *(a+j+tmm)/pivot);
for (j=i;j<mm1;j++) *(ttemp+j) = *(temp+j+1)
-(*(temp+i) * *(a+j+1+tmm)/pivot);
printf(" Reduced ttemp : \n");
for (j=0;j<mm1;j++) printf(" %f ",*(ttemp+j));
printf("\n");
}
}
else
{ /* Constraint 0 has not yet been
encountered */
jjval = 1;
}
/* allocate memory */
ap = (float *) malloc((sizeof (float))*nm1*mm1);
bp = (float *) malloc((sizeof (float))*mm1);
/* reduce a and b into ap and bp eliminating variable t and constraint i */
jj=-1;
for (j=0;j<mm;j++)
if (j != i)
{
jj=jj+1;
aj=a+j;
ajt=aj+tmm; *(bp+jj)= *(b+j) - *(bi) * *(ajt) / pivot;
apjj=ap+jj;
kk=-1;
for (k=0;k<nn;k++)
if (k != t)
{
kk=kk+1;
kmm=k*mm;
*(apjj+kk*mm1)= *(aj+kmm)- *(ajt) * *(ai+kmm)/ pivot;
}
}
/* add contribution to derivative from that calculated in smaller dimension */
/* Normal case method */
if(special == 0)
{
/*
deriv = 0.;
if (*(bi) != 0.)
{deriv=cdvda(ap,bp, &mm1, &nm1, &mm1, &nm1, &ttdim, ttemp, &jjval,code);}
*/
deriv=cdvda(ap,bp, &mm1, &nm1, &mm1, &nm1, &ttdim, ttemp, &jjval,code);
v=v+ *(bi)/amax*deriv/nn;
if (kval == 1 || *jval == 0) free(ttemp);
if(*code != 0)return (0.);
}
else /* Use special case method */
{
if( *jval == 1)
{
if(i == 0) /* This is constraint 0 */
{
deriv = 0.;
vol = 0.;
dvdb = 0.;
printf("\n special case 2: reduction by 0 and variable k\n");
for (j=1;j<mm;j++)
{
k = j - 1;
junk=cvolumebj(ap,bp,&mm1,&nm1,&mm1,&nm1,k,&dvdb);
if(junk == -1.)
{
*code = -1;
return(0.);
}
printf("\n j %d k %d dvdb %f vol %f",j,k,dvdb,junk);
deriv = deriv + dvdb * *(a + j + tmm) ;
vol = vol + junk;
}
if(nm1 == 1)vol = junk;
deriv = *(bi) * deriv/pivot;
deriv = (deriv - vol) /pivot;
v=v+ *(bi)/amax*deriv/nn;
printf("dcont %f",deriv);
printf("\n total volume of face = %f",vol);
/* start debug section */
FDa = *(bi)/amax*deriv/nn;
da = 0.01* *(a + *tdim * mm);
if(da == 0.0)da = 0.01;
*(a + *tdim * mm) = *(a + *tdim * mm) + da;
jj=-1;
for (j=0;j<mm;j++)
if (j != i)
{
jj=jj+1; aj=a+j;
ajt=aj+tmm;
*(bp+jj)= *(b+j) - *(bi) * *(ajt) / pivot;
apjj=ap+jj;
kk=-1;
for (k=0;k<nn;k++)
if (k != t)
{
kk=kk+1;
kmm=k*mm;
*(apjj+kk*mm1)= *(aj+kmm)- *(ajt) * *(ai+kmm)/ pivot;
}
}
volp = cvolume(ap,bp,&mm1,&nm1);
amax = fabs(pivot+da);
volp = (*(bi)/(amax*nn)) * volp;
*(a + *tdim * mm) = *(a + *tdim * mm) - 2. * da;
jj=-1;
for (j=0;j<mm;j++)
if (j != i)
{
jj=jj+1;
aj=a+j;
ajt=aj+tmm;
*(bp+jj)= *(b+j) - *(bi) * *(ajt) / pivot;
apjj=ap+jj;
kk=-1;
for (k=0;k<nn;k++)
if (k != t)
{
kk=kk+1;
kmm=k*mm;
*(apjj+kk*mm1)= *(aj+kmm)- *(ajt) * *(ai+kmm)/ pivot;
}
}
volm = cvolume(ap,bp,&mm1,&nm1);
amax = fabs(pivot-da);
*(a + *tdim * mm) = *(a + *tdim * mm) + da;
volm = (*(bi)/(amax*nn)) * volm;
FD = (volp-volm)/(2. * da);
printf("\n k %d a = %f da = %f",*tdim,*(a + *tdim * mm),da);
printf("\n pivot = %f",pivot);
printf("\n volp = %f",volp);
printf("\n volm = %f",volm);
printf("\n FD estimate = %f",FD);
printf("\n analytical estimate = %f\n",FDa);
/* end debug section */
}
else /* Constraint 0 not yet encountered */
{
printf("\n special case 1: constraint 0 not yet reached\n");
opt = 2;
junk=cvolumeb(ap,bp,&mm1,&nm1,&mm1,&nm1,&opt,&deriv);
if(junk == -1.)
{
*code = -1;
return(0.);
}
deriv= -(deriv * *(bi)/pivot);
v=v+ *(bi)/amax*deriv/nn;
printf("\n spec: 0 not yet reached i %d dcont %f",i,deriv);
/* debug section: reduce system and repeat calculation with finite difference */
FDa = *(bi)/amax*deriv/nn;
da = 0.01* *(a + *tdim * mm);
if(da == 0.0)da = 0.01;
*(a + *tdim * mm) = *(a + *tdim * mm) + da;
/* reduce system */
jj=-1;
for (j=0;j<mm;j++)
if (j != i)
{
jj=jj+1;
aj=a+j;
ajt=aj+tmm;
*(bp+jj)= *(b+j) - *(bi) * *(ajt) / pivot;
apjj=ap+jj;
kk=-1;
for (k=0;k<nn;k++)
if (k != t)
{
kk=kk+1;
kmm=k*mm;
*(apjj+kk*mm1)= *(aj+kmm)- *(ajt) * *(ai+kmm)/ pivot;
}
}
volp = cvolume(ap,bp,&mm1,&nm1);
*(a + *tdim * mm) = *(a + *tdim * mm) - 2. * da;
/* reduce system */
jj=-1;
for (j=0;j<mm;j++)
if (j != i)
{
jj=jj+1;
aj=a+j;
ajt=aj+tmm;
*(bp+jj)= *(b+j) - *(bi) * *(ajt) / pivot;
apjj=ap+jj;
kk=-1;
for (k=0;k<nn;k++)
if (k != t)
{
kk=kk+1;
kmm=k*mm;
*(apjj+kk*mm1)= *(aj+kmm)- *(ajt) * *(ai+kmm)/ pivot;
}
}
volm = cvolume(ap,bp,&mm1,&nm1);
*(a + *tdim * mm) = *(a + *tdim * mm) + da;
FD = (*(bi)/(amax*nn)) * (volp-volm)/(2. * da);
printf("\n k %d a = %f da = %f",*tdim,*(a + *tdim * mm),da);
printf("\n volp = %f",volp);
printf("\n volm = %f",volm);
printf("\n FD estimate = %f",FD);
printf("\n analytical estimate = %f\n",FDa);
/* end debug section */
}
}
else if( *jval == 0) /* Constraint 0 already encountered */
{
printf("\n special case 3: constraint 0 already reduced\n");
vol = 0.;
deriv = 0.;
for (j=0;j<i;j++)
{
junk=cvolumebj(ap,bp,&mm1,&nm1,&mm1,&nm1,j,&dvdb);
if(junk == -1.)
{
*code = -1;
return(0.);
} dbda = (*(a+j+tmm) * *(temp+i)/pivot) - *(temp+j);
deriv = deriv + dvdb*dbda;
vol = vol + junk;
}
for (j=i+1;j<mm;j++)
{
k = j - 1;
junk=cvolumebj(ap,bp,&mm1,&nm1,&mm1,&nm1,k,&dvdb);
if(junk == -1.)
{
*code = -1;
return(0.);
}
dbda = (*(a+j+tmm) * *(temp+i)/pivot) - *(temp+j);
deriv = deriv + dvdb*dbda;
vol = vol + junk;
}
if(nm1 == 1)vol = junk;
deriv = (*(bi) * deriv - *(temp+i)*vol)/pivot;
v=v+ *(bi)/amax*deriv/nn;
printf("\n spec: 0 already encountered i = %d dcont %f",i,deriv);
}
}
free(ap);
free(bp);
}
}
return(v);
}