SPModel: libDelaunay/src/volume.c Source File
VPAC - Computational Software Development
Main | SPModel | StGermain FrameWork |
Main Page | Alphabetical List | Class List | Directories | File List | Class Members | File Members

volume.c

Go to the documentation of this file.
00001 #include <stdlib.h>
00002 #include <ctype.h>
00003 #include <math.h>
00004 #include <stdio.h>
00005 /*--------------------------------------------------------------
00006 
00007 This file contains six recursive C routines for calculating
00008 the volume and derivatives of a convex polyhedron in n 
00009 dimensions. The routines volume, volumef, volumeb and dvda are 
00010 fortran callable interfaces to their C counterparts cvolume,
00011 cvolumeb and cdvda.
00012 
00013 ROUTINE      COMMENT            CALLS 
00014 
00015 cvolume  calculates volume of region    cvolume 
00016 which satisfies Ax <= b.
00017 
00018 cvolumef     calculates volume of region    cvolumef 
00019 which satisfies Ax <= b.
00020 (faster version, see below)
00021 
00022 cvolumeb     calculates volume and the      cvolume
00023 derivative d(vol)/db(0).
00024 
00025 cdvda        calculates derivatives     cdvda, cvolumeb, 
00026 d(vol)/da(0,j) (j=1,n).        & cvolumebj
00027 
00028 cdvdaf        calculates derivatives    cdvdaf, cvolumeb, 
00029 d(vol)/da(0,j) (j=1,n).        & cvolumebj
00030 (faster version, see below)
00031 
00032 cvolumebj    calculates partial volumes
00033 and derivatives d(vol)/d(b(j)
00034 (j=1,n). (called by cdvda.)
00035 
00036 Here b(j) is the jth element of vector b, and a(i,j) is the
00037 (i,j)th coefficient of the matrix A corresponding to the
00038 ith constraint and the jth dimension.
00039 
00040 The two faster versions cvolumef and cdvdaf do not continue down
00041 the recursive tree if b(j) = 0 for any j at any time. This avoids
00042 extra work but restricts the applicability of the routines to
00043 cases where the origin does not pass through any inconsistent
00044 constraints.
00045 
00046 Malcolm Sambridge, April 1995.
00047 
00048 --------------------------------------------------------------*/
00049 
00050 /*--------------------------------------------------------------
00051 
00052 ROUTINE: cvolume
00053 
00054 This routine calculates the `volume' in dimension n of the region
00055 bounded by a set of m linear inequality constraints of the form
00056 A x <= b, where a has m rows and n columns and is given by a(m,n),
00057 b is the n-vector and is contained in b(m). The recursive formula
00058 of Lasserre (1983) is used. Redundant constraints are allowed 
00059 If the inequality constraints are inconsistent then the volume 
00060 is returned as zero.  If the original polyhedron is unbounded 
00061 then a warning is issued and the volume is return as -1.0 
00062 (even though it is undefined). 
00063 
00064 Jean Braun and Malcolm Sambridge, Jan 1995.
00065 
00066 Lasserre, J. B., 1983. "An analytical expression and an Algorithm
00067 for the Volume of a Convex Polyhedron in R^n.", JOTA, vol. 39, 
00068 no. 3, 363-377.
00069 
00070 --------------------------------------------------------------*/
00071 
00072 float cvolume (a,b,m,n,mmax,nmax)
00073 
00074 int   *n, *m, *mmax, *nmax;
00075 float *a, *b; 
00076 
00077 {
00078     
00079     float v,amax,pivot;
00080     int   i,j,t,k,l;
00081     int   jj,kk;
00082     float   *ai, *aj, *ajt, *apjj, *bi;
00083     int   kmm,tmm;
00084     int   nn,mm,nn_max,mm_max;
00085     int   nm1,mm1;
00086     float  *ap, *bp;
00087     int   firstmin,firstmax;
00088     float bmin,bmax,bb;
00089     float partialv;
00090     FILE *file_ptr;
00091     
00092     
00093     nn= *n;
00094     mm= *m;
00095     nn_max= *nmax;
00096     mm_max= *mmax;
00097     file_ptr = NULL;
00098     
00099     //file_ptr = fopen("volume.out", "w");
00100     //printf ("nn = %d\tmm = %d\tnn_max = %d\tmm_max = %d\n", nn, mm, nn_max, mm_max);
00101     
00102     /* one-dimensional case (full reduction) */
00103     
00104     if (nn == 1) 
00105     {
00106         firstmin=0;
00107         firstmax=0;
00108         for (l=0;l<mm;l++)
00109         {
00110             //printf ("\n*a = %f\nl = %d\n*(a+l) = %f\n", *a, l, *(a+l));
00111             //printf ("*b = %f\nl = %d\n*(b+l) = %f\n", *b, l, *(b+l));
00112             
00113             if ( *(a+l) > 0.) 
00114             {
00115                 //  printf ("*a = %f\nl = %d\n*(a+l) = %f\n", *a, l, *(a+l));
00116                 //  printf ("*b = %f\nl = %d\n*(b+l) = %f\n", *b, l, *(b+l));
00117                 bb= *(b+l)/ *(a+l);
00118             if (firstmin==0) {firstmin=1;bmin=bb;}
00119                 else if (bb<bmin) bmin=bb;
00120             }
00121             else if ( *(a+l) < 0.)
00122             {
00123                 //  printf ("*a = %f\nl = %d\n*(a+l) = %f\n", *a, l, *(a+l));
00124                 //  printf ("*b = %f\nl = %d\n*(b+l) = %f\n", *b, l, *(b+l));
00125                 bb= *(b+l)/ *(a+l);
00126             if (firstmax==0) {firstmax=1;bmax=bb;}
00127                 else if (bb>bmax) bmax=bb;
00128             }
00129             else if ( *(b+l) < 0.)         /* Constraints are inconsistent */
00130             { 
00131                 //  printf ("*a = %f\nl = %d\n*(a+l) = %f\n", *a, l, *(a+l));
00132                 //  printf ("*b = %f\nl = %d\n*(b+l) = %f\n", *b, l, *(b+l));
00133                 printf ("HEY Inconsistent constraints found after reduction to n = 1 \n\n");
00134                 return(0.); 
00135             }
00136         }
00137         v=0.;
00138         if (firstmin*firstmax == 1) v=bmin-bmax;
00139         else 
00140         {
00141             /*     printf("Volume is unbounded at 1-D; volume returned as -1\n"); */
00142             return(-1.0);
00143         }
00144         
00145     if (v<0.) {v=0.;} 
00146         
00147         return(v);
00148     }
00149     
00150     nm1=nn-1;
00151     mm1=mm-1;
00152     v=0.;
00153     
00154     for (i=0;i<mm;i++)
00155     {
00156         ai=a+i;
00157         bi=b+i;
00158         
00159         /* find largest pivot */
00160         
00161         amax=0.;
00162         for (j=0;j<nn;j++) 
00163     if (fabs( *(ai+j*mm_max)) >= amax) {amax= fabs( *(ai+j*mm_max)); t=j;}
00164         tmm=t*mm_max;
00165         pivot=*(ai+tmm);
00166         
00167         /* finds contribution to v from this pivot (if not nil) */
00168         
00169         if (amax == 0.)
00170         {
00171             /* Constraint is inconsistent */
00172             
00173             if ( *(bi) < 0.) 
00174         { printf("Constraint %d is inconsistent\n",i+1); return(0.); }
00175             
00176             /* otherwise constraint is redundant */
00177             
00178             printf("Constraint %d is redundant\n",i+1); 
00179         }
00180         else
00181         {
00182             
00183             /* allocate memory */
00184             
00185             ap = (float *) malloc(4*nm1*mm1);
00186             bp = (float *) malloc(4*mm1);
00187             
00188             /* reduce a and b into ap and bp eliminating variable t and constraint i */
00189             
00190             jj=-1;
00191             for (j=0;j<mm;j++)
00192             if (j != i)
00193             {
00194                 jj=jj+1;
00195                 aj=a+j;
00196                 ajt=aj+tmm;
00197                 *(bp+jj)= *(b+j) - *(bi) * *(ajt) / pivot;
00198                 apjj=ap+jj;
00199                 kk=-1;
00200                 for (k=0;k<nn;k++)
00201                 if (k != t)
00202                 {
00203                     kk=kk+1;
00204                     kmm=k*mm_max;
00205                     *(apjj+kk*mm1)= *(aj+kmm)- *(ajt) * *(ai+kmm)/ pivot;
00206                 }
00207             }
00208             
00209             /* add contribution to volume from volume calculated in smaller dimension */
00210             
00211             partialv=cvolume(ap,bp, &mm1, &nm1, &mm1, &nm1);
00212             if(partialv == -1.0)return(-1.0);
00213             v=v+ *(bi)/amax*partialv/nn; 
00214             
00215             free(ap);
00216             free(bp);
00217         }
00218     }
00219     
00220     return(v);
00221     
00222 }
00223 
00224 /*--------------------------------------------------------------
00225 
00226 ROUTINE: volume
00227 
00228 A dummy routine used to call cvolume from a fortran routine 
00229 
00230 
00231 Jean Braun and Malcolm Sambridge, Jan 1995.
00232 
00233 ----------------------------------------------------------------*/
00234 
00235 void volume_ (a,b,m,n,mmax,nmax,volume)
00236 
00237 int   *n, *m, *mmax, *nmax;
00238 float *a, *b;
00239 float *volume; 
00240 
00241 {
00242     
00243     *volume=cvolume(a,b,m,n,mmax,nmax);
00244     
00245 }
00246 
00247 /*--------------------------------------------------------------
00248 
00249 ROUTINE: cvolumeb
00250 
00251 This routine calculates the `volume' in dimension n of the region
00252 bounded by a set of m linear inequality constraints of the form
00253 A x <= b, where a has m rows and n columns and is given by a(m,n),
00254 b is the n-vector and is contained in b(m). The recursive formula
00255 of Lasserre (1983) is used. Redundant constraints are allowed and
00256 a warning is issued if any are encountered. If the inequality 
00257 constraints are inconsistent then the volume is returned as zero.
00258 If the original polyhedron is unbounded then a warning is issued
00259 and the volume is return as zero (even though it is undefined). 
00260 
00261 This version also calculates the derivative of the volume with
00262 respect to the parameter b(0) using the simple formula of
00263 Lasserre (1983). If *opt == 2 then only the derivative is 
00264 calculated and not the volume.
00265 
00266 This routine is a variation from the routine `cvolume'.
00267 
00268 Calls are made to routine cvolume.
00269 
00270 Malcolm Sambridge, March 1995.
00271 
00272 
00273 --------------------------------------------------------------*/
00274 
00275 float cvolumeb (a,b,m,n,mmax,nmax,opt,dvdb)
00276 
00277 int   *n, *m, *opt, *mmax, *nmax;
00278 float *a, *b; 
00279 float *dvdb; 
00280 
00281 {
00282     
00283     float v,amax,pivot;
00284     int   i,j,t,k,l;
00285     int   jj,kk;
00286     float   *ai, *aj, *ajt, *apjj, *bi;
00287     int   kmm,tmm;
00288     int   nn,mm,nn_max,mm_max;
00289     int   nm1,mm1;
00290     int   lmin,lmax;
00291     float  *ap, *bp;
00292     int   firstmin,firstmax;
00293     float bmin,bmax,bb;
00294     float vol;
00295     
00296     nn= *n;
00297     mm= *m;
00298     nn_max= *nmax;
00299     mm_max= *mmax;
00300     
00301     /* one-dimensional case (full reduction) */
00302     
00303     if (nn == 1) 
00304     {
00305         firstmin=0;
00306         firstmax=0;
00307         lmax=0;
00308         lmin=0;
00309         for (l=0;l<mm;l++)
00310         {
00311             if ( *(a+l) > 0.) 
00312             {
00313                 bb= *(b+l)/ *(a+l);
00314             if (firstmin==0) {firstmin=1;bmin=bb;lmin=l;}
00315             else if (bb<bmin) {bmin=bb;lmin=l;}
00316             }
00317             else if ( *(a+l) < 0.)
00318             {
00319                 bb= *(b+l)/ *(a+l);
00320             if (firstmax==0) {firstmax=1;bmax=bb;lmax=l;}
00321             else if (bb>bmax) {bmax=bb;lmax=l;}
00322             }
00323             else if ( *(b+l) < 0.) 
00324             { 
00325                 /* Constraints are inconsistent.  
00326                 Set volume and derivative to zero. */
00327                 
00328                 printf("Inconsistent constraints found after reduction to n = 1 \n");
00329                 *dvdb = 0.; 
00330                 return(0.);
00331             }
00332         }
00333         v=0.;
00334         *dvdb = 0.; 
00335         if (firstmin*firstmax == 1) v=bmin-bmax;
00336         else 
00337         {
00338             /*     printf("Volume is unbounded; volume returned as -1\n derivatives returned as zero\n"); */
00339             return(-1.0);
00340         }
00341         
00342     if (v<0.) {v=0.;*dvdb=0.;} 
00343         else if (v>0. && lmin == 0) *dvdb = 1. / *a;
00344         else if (v>0. && lmax == 0) *dvdb = -1. / *a;
00345         return(v);
00346     }
00347     
00348     nm1=nn-1;
00349     mm1=mm-1;
00350     v=0.;
00351     
00352     /*  perform main loop over constraints */
00353     
00354     for (i=0;i<mm;i++)
00355     {
00356         ai=a+i;
00357         bi=b+i;
00358         
00359         /* find largest pivot */
00360         
00361         amax=0.;
00362         for (j=0;j<nn;j++) 
00363     if (fabs( *(ai+j*mm_max)) >= amax) {amax= fabs( *(ai+j*mm_max)); t=j;}
00364         tmm=t*mm_max;
00365         pivot=*(ai+tmm);
00366         
00367         /* finds contribution to v from this pivot (if not nil) */
00368         
00369         if (amax == 0.)
00370         {
00371             /* Constraint is inconsistent */
00372             
00373             if ( *(bi) < 0.) 
00374             {
00375                 printf("Constraint %d is inconsistent\n",i+1); 
00376                 return(0.);
00377             } 
00378             
00379             /* otherwise constraint is redundant */
00380             
00381             printf("Constraint %d is redundant\n",i+1); 
00382         }
00383         else
00384         {
00385             
00386             /* allocate memory */
00387             
00388             ap = (float *) malloc(4*nm1*mm1);
00389             bp = (float *) malloc(4*mm1);
00390             
00391             /* reduce a and b into ap and bp eliminating variable t and constraint i */
00392             
00393             jj=-1;
00394             for (j=0;j<mm;j++)
00395             if (j != i)
00396             {
00397                 jj=jj+1;
00398                 aj=a+j;
00399                 ajt=aj+tmm;
00400                 *(bp+jj)= *(b+j) - *(bi) * *(ajt) / pivot;
00401                 apjj=ap+jj;
00402                 kk=-1;
00403                 for (k=0;k<nn;k++)
00404                 if (k != t)
00405                 {
00406                     kk=kk+1;
00407                     kmm=k*mm_max;
00408                     *(apjj+kk*mm1)= *(aj+kmm)- *(ajt) * *(ai+kmm)/ pivot;
00409                 }
00410             }
00411             
00412             /* add contribution to volume from volume calculated in smaller dimension */
00413             
00414             vol=cvolume(ap,bp, &mm1, &nm1, &mm1, &nm1);
00415             if(vol == -1.0)
00416             {
00417                 *dvdb = 0.;
00418                 return(-1.0);
00419             }
00420             v=v+ *(bi)/amax*vol/nn;
00421             
00422             free(ap);
00423             free(bp);
00424             
00425             /* calculate derivatives for first constraint only */
00426             
00427             /* calculate volume and derivative 
00428             with respect to b_0 */
00429             if (i == 0) *dvdb = vol/amax; 
00430             /* calculate derivative with respect 
00431             to b_0 but not volume */
00432             if (*opt == 2) return (0.); 
00433         }
00434     }
00435     
00436     return(v);
00437     
00438 }
00439 
00440 /*--------------------------------------------------------------
00441 
00442 ROUTINE: volumeb
00443 
00444 A dummy routine used to call cvolumeb from a fortran routine 
00445 
00446 ----------------------------------------------------------------*/
00447 
00448 void volumeb_ (a,b,m,n,mmax,nmax,opt,volume,dvdb)
00449 
00450 int   *n, *m, *mmax, *nmax, *opt;
00451 float *a, *b;
00452 float *volume; 
00453 float *dvdb; 
00454 
00455 {
00456     
00457     if(*opt == 1)                      /* calculate volume and derivative with
00458     respect to b(0) */
00459     {
00460         *volume=cvolumeb(a,b,m,n,mmax,nmax,opt,dvdb);
00461     }
00462     else if(*opt == 2)                /* calculate derivative with respect to b(0) 
00463     but not volume */
00464     {
00465         *volume=cvolumeb(a,b,m,n,mmax,nmax,opt,dvdb);
00466     }
00467     else 
00468     {
00469         printf(" Warning: routine volumeb called with an invalid option parameter\n");
00470         printf("          valid options are 1 and 2, option given = %d \n",*opt);
00471     } 
00472     
00473 }
00474 
00475 /*--------------------------------------------------------------
00476 
00477 ROUTINE: cvolumebj
00478 
00479 This routine only calculates the j-th outer loop of the recursive
00480 routine cvolumeb. It returns the partial volume for projection 
00481 onto the j-th constraint and the derivative of the total volume 
00482 with respect to parameter b(j) using the simple formula of 
00483 Lasserre(1983). (See cvolumeb for more details.)
00484 
00485 The reason for this routine is so that the derivatives
00486 with respect to parameter b(j) can be calculated and passed back
00487 to a fortran routine without creating and an extra array of size
00488 m (because only one derivative is calculated per call). In
00489 this way it also avoids recalculating the total volume m times
00490 (which would be the case if we used a simple variation of routine
00491 cvolumeb).
00492 
00493 To calculate the total volume this routine must be called m times
00494 and the partial volume summed
00495 
00496 Constraint j is determined by the value of *con.
00497 
00498 Calls are made to routine cvolume.
00499 
00500 Malcolm Sambridge, March 1995.
00501 
00502 --------------------------------------------------------------*/
00503 
00504 float cvolumebj (a,b,m,n,mmax,nmax,con,dvdb)
00505 
00506 int   *n, *m, *mmax, *nmax, con;
00507 float *a, *b; 
00508 float *dvdb; 
00509 
00510 {
00511     
00512     float v,amax,pivot;
00513     int   i,j,t,k,l;
00514     int   jj,kk;
00515     float   *ai, *aj, *ajt, *apjj, *bi;
00516     int   kmm,tmm;
00517     int   nn,mm,nn_max,mm_max;
00518     int   nm1,mm1;
00519     int   lmin,lmax;
00520     float  *ap, *bp;
00521     int   firstmin,firstmax;
00522     float bmin,bmax,bb;
00523     float vol;
00524     
00525     nn= *n;
00526     mm= *m;
00527     nn_max= *nmax;
00528     mm_max= *mmax;
00529     
00530     /* one-dimensional case (full reduction) */
00531     
00532     if (nn == 1) 
00533     {
00534         firstmin=0;
00535         firstmax=0;
00536         lmax=0;
00537         lmin=0;
00538         for (l=0;l<mm;l++) 
00539         {
00540             if ( *(a+l) > 0.) 
00541             {
00542                 bb= *(b+l)/ *(a+l);
00543             if (firstmin==0) {firstmin=1;bmin=bb;lmin=l;}
00544             else if (bb<bmin) {bmin=bb;lmin=l;}
00545             }
00546             else if ( *(a+l) < 0.)
00547             {
00548                 bb= *(b+l)/ *(a+l);
00549             if (firstmax==0) {firstmax=1;bmax=bb;lmax=l;}
00550             else if (bb>bmax) {bmax=bb;lmax=l;}
00551             }
00552             else if ( *(b+l) < 0.) 
00553             { 
00554                 /* Constraints are inconsistent.  
00555                 Set volume and derivative to zero. */
00556                 
00557                 printf("Inconsistent constraints found after reduction to n = 1 \n");
00558                 *dvdb = 0.; 
00559                 return(0.);
00560             }
00561         }
00562         v=0.;
00563         *dvdb = 0.; 
00564         if (firstmin*firstmax == 1) v=bmin-bmax;
00565         else 
00566         {
00567             /*     printf("Volume is unbounded; volume returned as -1\n derivatives returned as zero\n"); */
00568             return(-1.0);
00569         }
00570         
00571     if (v<0.) {v=0.;*dvdb=0.;} 
00572         else if (v>0. && lmin == con) *dvdb = 1. / *(a+lmin);
00573         else if (v>0. && lmax == con) *dvdb = -1. / *(a+lmax);
00574         return(v);
00575     }
00576     
00577     nm1=nn-1;
00578     mm1=mm-1;
00579     v=0.;
00580     
00581     /*  perform main loop over constraints */
00582     
00583     /* for (i=0;i<mm;i++) */
00584     
00585     i = con;
00586     
00587     ai=a+i;
00588     bi=b+i;
00589     
00590     /* find largest pivot */
00591     
00592     amax=0.;
00593     for (j=0;j<nn;j++) 
00594 if (fabs( *(ai+j*mm_max)) >= amax) {amax= fabs( *(ai+j*mm_max)); t=j;}
00595     tmm=t*mm_max;
00596     pivot=*(ai+tmm);
00597     
00598     /* finds contribution to v from this pivot (if not nil) */
00599     
00600     if (amax == 0.)
00601     {
00602         /* Constraint is inconsistent */
00603         
00604         if ( *(bi) < 0.) 
00605         {
00606             printf("Constraint %d is inconsistent\n",i+1); 
00607             *dvdb = 0.; 
00608             return(0.);
00609         } 
00610         
00611         /* otherwise constraint is redundant */
00612         
00613         printf("Constraint %d is redundant\n",i+1); 
00614         *dvdb = 0.; 
00615     }
00616     else
00617     {
00618         
00619         /* allocate memory */
00620         
00621         ap = (float *) malloc(4*nm1*mm1);
00622         bp = (float *) malloc(4*mm1);
00623         
00624         /* reduce a and b into ap and bp eliminating variable t and constraint i */
00625         
00626         jj=-1;
00627         for (j=0;j<mm;j++)
00628         if (j != i)
00629         {
00630             jj=jj+1;
00631             aj=a+j;
00632             ajt=aj+tmm;
00633             *(bp+jj)= *(b+j) - *(bi) * *(ajt) / pivot;
00634             apjj=ap+jj;
00635             kk=-1;
00636             for (k=0;k<nn;k++)
00637             if (k != t)
00638             {
00639                 kk=kk+1;
00640                 kmm=k*mm_max;
00641                 *(apjj+kk*mm1)= *(aj+kmm)- *(ajt) * *(ai+kmm)/ pivot;
00642             }
00643         }
00644         
00645         /* calculate partial volume */ 
00646         
00647         vol=cvolume(ap,bp, &mm1, &nm1, &mm1, &nm1);
00648         if(vol == -1.0)
00649         {
00650             *dvdb = 0.;
00651             return(-1.0);
00652         }
00653         v=*(bi)/amax*vol/nn;
00654         
00655         /* calculate derivative of total volume
00656         with respect to current constraint */
00657         *dvdb = vol/amax; 
00658         
00659         free(ap);
00660         free(bp);
00661     }
00662     
00663     return(v);
00664     
00665 }
00666 
00667 /*--------------------------------------------------------------
00668 
00669 ROUTINE: volumebj
00670 
00671 A dummy routine used to call cvolumebj from a fortran routine
00672 
00673 ----------------------------------------------------------------*/
00674 
00675 void volumebj_ (a,b,m,n,mmax,nmax,con,volume,dvdb)
00676 
00677 int   *n, *m, *mmax, *nmax,*con;
00678 float *a, *b;
00679 float *volume; 
00680 float *dvdb; 
00681 
00682 {
00683     
00684     int tcon;
00685     tcon = *con - 1;
00686     /* calculate partial volume and 
00687     derivative with respect to b(con) */
00688     
00689     *volume=cvolumebj(a,b,m,n,mmax,nmax,tcon,dvdb);
00690     
00691 }
00692 /*--------------------------------------------------------------
00693 
00694 ROUTINE: cdvda
00695 
00696 This routine calculates the derivative with respect to a(0,tdim)
00697 of the `volume' in dimension n of the region bounded by a set 
00698 of m linear inequality constraints of the form A x <= b, where 
00699 a has m rows and n columns and is given by a(m,n), b is the 
00700 n-vector and is contained in b(m). The derivative expression
00701 is recursive and derived from the formula of Lasserre (1983). 
00702 
00703 Redundant constraints are allowed and a warning is issued if any are 
00704 encountered. If the inequality constraints are inconsistent then the 
00705 derivative is returned as zero.  If any constraint is orthogonal to 
00706 the component a(0,idim) then the reduction can only take place onto 
00707 variable idim. A special case is used to handle this which involves
00708 no further recursive calls.
00709 
00710 If the original polyhedron is unbounded then a warning is issued
00711 and the derivative is return as zero. 
00712 
00713 Note: This code takes advantage of the fact that during recursive calls
00714 constraint 0 does not change its position in the list of remaining 
00715 constraints (if it has not been eliminated), i.e. it is always the 
00716 first constraint.  This would not be the case if the algorithm were 
00717 adapated to deal with other constraints, i.e. evaluate dvda_i,j 
00718 where i .ne. 0.
00719 
00720 Calls itself cvolumeb, and cvolumebj.
00721 
00722 Malcolm Sambridge, March 1995.
00723 
00724 --------------------------------------------------------------*/
00725 
00726 float cdvda (a,b,m,n,mmax,nmax,tdim,temp,jval,code)
00727 
00728 int   *n, *m, *nmax, *mmax, *tdim, *jval, *code;
00729 float *a, *b, *temp; 
00730 
00731 {
00732     
00733     float v,amax,pivot;
00734     int   i,j,t,k,l;
00735     int   jj,kk;
00736     float   *ai, *aj, *ajt, *apjj, *bi;
00737     int   kmm,tmm;
00738     int   nn,mm,nn_max,mm_max,ttdim;
00739     int   nm1,mm1;
00740     int   lmin,lmax, jjval, kval;
00741     float  *ap, *bp, *ttemp;
00742     int   firstmin,firstmax;
00743     float bmin,bmax,bb;
00744     float deriv, junk, vol, dvdb, dbda;
00745     int   special, opt;
00746     
00747     nn= *n;
00748     mm= *m;
00749     nn_max= *nmax;
00750     mm_max= *mmax;
00751     
00752     /* one-dimensional case (full reduction) */
00753     
00754     *code = 0;
00755     
00756     if (nn == 1) 
00757     {
00758         firstmin=0;
00759         firstmax=0;
00760         lmax=0;
00761         lmin=0;
00762         for (l=0;l<mm;l++)
00763         {
00764             if ( *(a+l) > 0.) 
00765             {
00766                 bb= *(b+l)/ *(a+l);
00767             if (firstmin==0) {firstmin=1;bmin=bb;lmin=l;}
00768             else if (bb<bmin) {bmin=bb;lmin=l;}
00769             }
00770             else if ( *(a+l) < 0.)
00771             {
00772                 bb= *(b+l)/ *(a+l);
00773             if (firstmax==0) {firstmax=1;bmax=bb;lmax=l;}
00774             else if (bb>bmax) {bmax=bb;lmax=l;}
00775             }
00776             else if ( *(b+l) < 0.) 
00777             {
00778                 /* Constraint is inconsistent.
00779                 Set derivative to zero. */
00780                 
00781                 printf("Inconsistent constraints found after reduction to n = 1 \n");
00782                 *code = 1;
00783                 return(0.);
00784             }
00785         }
00786         v=0.;
00787         if (firstmin*firstmax == 1) v=bmin-bmax;
00788         else 
00789         {
00790             /*     printf("Volume is unbounded; derivative returned is zero\n"); */
00791             *code = -1;
00792             return(0.);
00793         }
00794         
00795         if (v<0.) return(0.);       
00796         
00797         if(*jval == 1)           /* Constraint 0 has not yet been encountered */
00798         {
00799             if (lmin == 0) deriv = -bmin/ *a;
00800             else if (lmax == 0) deriv = bmax/ *a;
00801             else deriv = 0.;
00802             return(deriv); 
00803         }
00804         else if(*jval == 0)     /* Constraint 0  has already been encountered */
00805         {
00806             deriv =  ( *(temp+lmax) * (bmax/ *(a+lmax)) ) -
00807             ( *(temp+lmin) * (bmin/ *(a+lmin)) );
00808             return(deriv); 
00809         }
00810     }
00811     nm1=nn-1;
00812     mm1=mm-1;
00813     v=0.;
00814     
00815     /*  perform main loop over constraints */
00816     
00817     for (i=0;i<mm;i++)
00818     {
00819         ai=a+i;
00820         bi=b+i;
00821         ttdim = *tdim;
00822         special = 0;
00823         
00824         /* find largest pivot */
00825         
00826         amax=0.;
00827         t = 0;
00828         for (j=0;j<nn;j++) 
00829         if (fabs( *(ai+j*mm_max)) >= amax && j != ttdim) 
00830     {amax= fabs( *(ai+j*mm_max)); t=j;}
00831         
00832         /* finds contribution to v from 
00833         this pivot (if not nil) */
00834         
00835         if (amax == 0.)
00836         {
00837             if(*(ai + ttdim * mm_max) == 0.0) 
00838             { 
00839                 /* Constraint is inconsistent */
00840                 if ( *(bi) < 0.)
00841                 {
00842                     printf("Constraint %d is inconsistent\n",i+1);
00843                     *code = 1;
00844                     return(0.);
00845                 }
00846                 
00847                 /* otherwise constraint is redundant */
00848                 
00849                 printf("Constraint %d is redundant\n",i+1); 
00850             }
00851             else
00852             {
00853                 /* if projection can only be peformed
00854                 on dimension tdim then activate 
00855                 special case */ 
00856                 
00857                 special = 1;
00858                 t = ttdim;
00859                 amax = fabs(*(ai+t * mm_max));
00860                 
00861             }
00862         }
00863         
00864         tmm=t*mm_max;
00865         pivot=*(ai+tmm);
00866         
00867         if(t < ttdim) ttdim = ttdim -1;
00868         
00869         
00870         if (amax != 0)
00871         {
00872             
00873             /* determine if constraint 0 has been encountered */
00874             
00875             kval = 0;
00876             if ( i == 0 && *jval == 1)
00877             {
00878                 /* This is the first encounter of
00879                 constraint 0 on this path so we 
00880                 allocate memory and store parameters 
00881                 to be used when n = 1 */
00882                 
00883                 if (special == 0) 
00884                 {
00885                     ttemp = (float *) malloc(4*mm1);
00886                     for (j=0;j<mm1;j++) *(ttemp+j) = - *(a+j+1+tmm)/pivot;
00887                     kval = 1;
00888                 }
00889                 jjval = 0;
00890             }
00891             else if (*jval == 0)             /* Constraint 0 has already been
00892             encountered */
00893             {
00894                 jjval = 0;  
00895                 /* perform recursive update of component
00896                 derivative array temp. This eliminates 
00897                 row i and copies into a new vector */
00898                 
00899                 if(special == 0)
00900                 {
00901                     ttemp = (float *) malloc(4*mm1);
00902                     for (j=0;j<i;j++) *(ttemp+j) = *(temp+j) 
00903                     -(*(temp+i) * *(a+j+tmm)/pivot);
00904                     for (j=i;j<mm1;j++) *(ttemp+j) = *(temp+j+1)
00905                     -(*(temp+i) * *(a+j+1+tmm)/pivot);
00906                 }
00907             }
00908             else 
00909             {                             /* Constraint 0 has not yet been 
00910                 encountered */
00911                 jjval = 1;                
00912             }
00913             
00914             
00915             /* allocate memory */
00916             
00917             ap = (float *) malloc(4*nm1*mm1);
00918             bp = (float *) malloc(4*mm1);
00919             
00920             /* reduce a and b into ap and bp eliminating variable t and constraint i */
00921             
00922             jj=-1;
00923             for (j=0;j<mm;j++)
00924             if (j != i)
00925             {
00926                 jj=jj+1;
00927                 aj=a+j;
00928                 ajt=aj+tmm;
00929                 *(bp+jj)= *(b+j) - *(bi) * *(ajt) / pivot;
00930                 apjj=ap+jj;
00931                 kk=-1;
00932                 for (k=0;k<nn;k++)
00933                 if (k != t)
00934                 {
00935                     kk=kk+1;
00936                     kmm=k*mm_max;
00937                     *(apjj+kk*mm1)= *(aj+kmm)- *(ajt) * *(ai+kmm)/ pivot;
00938                 }
00939             }
00940             
00941             /* add contribution to derivative from that calculated in smaller dimension */
00942             
00943             /* Normal case method */
00944             if(special == 0)
00945             {
00946                 deriv=cdvda(ap,bp, &mm1, &nm1, &mm1, &nm1, &ttdim, ttemp, &jjval,code);
00947                 v=v+ *(bi)/amax*deriv/nn;
00948                 if (kval == 1 || *jval == 0) free(ttemp);
00949                 if(*code != 0)return (0.);
00950                 
00951             }
00952             else                    /* Use special case method */
00953             {
00954                 if( *jval == 1)                    
00955                 {                     
00956                     if(i == 0)                      /* This is constraint 0 */
00957                     {
00958                         deriv = 0.; 
00959                         vol = 0.; 
00960                         dvdb = 0.;
00961                         for (j=1;j<mm;j++) 
00962                         {
00963                             k = j - 1;
00964                             junk=cvolumebj(ap,bp,&mm1,&nm1,&mm1,&nm1,k,&dvdb);
00965                             if(junk == -1.)
00966                             {
00967                                 *code = -1;
00968                                 return(0.);
00969                             }
00970                             deriv = deriv + dvdb * *(a + j + tmm) ;
00971                             vol = vol + junk;
00972                         }
00973                         if(nm1 == 1)vol = junk;
00974                         deriv = *(bi) * deriv/pivot;
00975                         deriv = (deriv - vol) /pivot; 
00976                         v=v+ *(bi)/amax*deriv/nn;
00977                         
00978                     }
00979                     else                /* Constraint 0 not yet encountered */
00980                     {
00981                         opt = 2;
00982                         junk=cvolumeb(ap,bp,&mm1,&nm1,&mm1,&nm1,&opt,&deriv);
00983                         if(junk == -1.)
00984                         {
00985                             *code = -1;
00986                             return(0.);
00987                         }
00988                         deriv= -(deriv *  *(bi)/pivot);
00989                         v=v+ *(bi)/amax*deriv/nn;
00990                     }
00991                 }
00992                 else if( *jval == 0)               /* Constraint 0 already encountered */
00993                 {
00994                     vol = 0.; 
00995                     deriv = 0.; 
00996                     for (j=0;j<i;j++) 
00997                     {
00998                         junk=cvolumebj(ap,bp,&mm1,&nm1,&mm1,&nm1,j,&dvdb);
00999                         if(junk == -1.)
01000                         {
01001                             *code = -1;
01002                             return(0.);
01003                         }
01004                         dbda = (*(a+j+tmm) * *(temp+i)/pivot) - *(temp+j);
01005                         deriv = deriv + dvdb*dbda;
01006                         vol = vol + junk;
01007                     }
01008                     for (j=i+1;j<mm;j++) 
01009                     {
01010                         k = j - 1;
01011                         junk=cvolumebj(ap,bp,&mm1,&nm1,&mm1,&nm1,k,&dvdb);
01012                         if(junk == -1.)
01013                         {
01014                             *code = -1;
01015                             return(0.);
01016                         }
01017                         dbda = (*(a+j+tmm) * *(temp+i)/pivot) - *(temp+j);
01018                         deriv = deriv + dvdb*dbda;
01019                         vol = vol + junk;
01020                     }
01021                     if(nm1 == 1)vol = junk;
01022                     deriv = (*(bi) * deriv - *(temp+i)*vol)/pivot;
01023                     v=v+ *(bi)/amax*deriv/nn;
01024                 }
01025             }
01026             
01027             free(ap);
01028             free(bp);
01029         }
01030     }
01031     
01032     return(v);
01033     
01034 }
01035 
01036 /*--------------------------------------------------------------
01037 
01038 ROUTINE: dvda
01039 
01040 A dummy routine used to call cdvda from a fortran routine 
01041 
01042 ----------------------------------------------------------------*/
01043 
01044 
01045 void dvda_ (a,b,m,n,mmax,nmax,idim,dvda,code)
01046 
01047 int   *n, *m, *mmax, *nmax, *idim, *code;
01048 float *a, *b;
01049 float *dvda; 
01050 
01051 {
01052     
01053     int    jval, tdim;
01054     // tentative try, Feng VPAC Oct. 3, 2002, initialize temp as zero (?)
01055     float *temp=0;
01056     
01057     
01058     jval = 1;
01059     tdim = *idim - 1;
01060     
01061     *dvda=cdvda(a,b,m,n,mmax,nmax,&tdim,temp,&jval,code);
01062     
01063     free(temp);
01064 }
01065 
01066 /*--------------------------------------------------------------
01067 
01068 ROUTINE: cvolumef
01069 
01070 This routine calculates the `volume' in dimension n of the region
01071 bounded by a set of m linear inequality constraints of the form
01072 A x <= b, where a has m rows and n columns and is given by a(m,n),
01073 b is the n-vector and is contained in b(m). The recursive formula
01074 of Lasserre (1983) is used. Redundant constraints are allowed 
01075 If the inequality constraints are inconsistent then the volume 
01076 is returned as zero.  If the original polyhedron is unbounded 
01077 then a warning is issued and the volume is return as -1.0 
01078 (even though it is undefined). 
01079 
01080 This is a faster version which does not continue down the
01081 recursive tree if if encounters b(j) = 0 for any j. This
01082 restricts its use to cases where the origin does not
01083 pass through any inconsistent constraints. Also redundant
01084 constraints that pass through the origin will not be detected.
01085 
01086 Jean Braun and Malcolm Sambridge, Jan 1995.
01087 
01088 
01089 Lasserre, J. B., 1983. "An analytical expression and an Algorithm
01090 for the Volume of a Convex Polyhedron in R^n.", JOTA, vol. 39, 
01091 no. 3, 363-377.
01092 
01093 --------------------------------------------------------------*/
01094 
01095 float cvolumef (a,b,m,n,mmax,nmax)
01096 
01097 int   *n, *m, *mmax, *nmax;
01098 float *a, *b; 
01099 
01100 {
01101     
01102     float v,amax,pivot;
01103     int   i,j,t,k,l;
01104     int   jj,kk;
01105     float   *ai, *aj, *ajt, *apjj, *bi;
01106     int   kmm,tmm;
01107     int   nn,mm,nn_max,mm_max;
01108     int   nm1,mm1;
01109     float  *ap, *bp;
01110     int   firstmin,firstmax;
01111     float bmin,bmax,bb;
01112     float partialv;
01113     
01114     nn= *n;
01115     mm= *m;
01116     nn_max= *nmax;
01117     mm_max= *mmax;
01118     
01119     /* one-dimensional case (full reduction) */
01120     
01121     if (nn == 1) 
01122     {
01123         firstmin=0;
01124         firstmax=0;
01125         for (l=0;l<mm;l++)
01126         {
01127             if ( *(a+l) > 0.) 
01128             {
01129                 bb= *(b+l)/ *(a+l);
01130             if (firstmin==0) {firstmin=1;bmin=bb;}
01131                 else if (bb<bmin) bmin=bb;
01132             }
01133             else if ( *(a+l) < 0.)
01134             {
01135                 bb= *(b+l)/ *(a+l);
01136             if (firstmax==0) {firstmax=1;bmax=bb;}
01137                 else if (bb>bmax) bmax=bb;
01138             }
01139             else if ( *(b+l) < 0.)         /* Constraints are inconsistent */
01140             { 
01141                 printf("Inconsistent constraints found after reduction to n = 1 \n");
01142                 return(0.); 
01143             }
01144         }
01145         v=0.;
01146         if (firstmin*firstmax == 1) v=bmin-bmax;
01147         else 
01148         {
01149             /*     printf("Volume is unbounded at 1-D; volume returned as -1\n"); */
01150             return(-1.0);
01151         }
01152         
01153     if (v<0.) {v=0.;} 
01154         
01155         return(v);
01156     }
01157     
01158     nm1=nn-1;
01159     mm1=mm-1;
01160     v=0.;
01161     
01162     for (i=0;i<mm;i++)
01163     {
01164         ai=a+i;
01165         bi=b+i;
01166         
01167         /* find largest pivot */
01168         
01169         amax=0.;
01170         for (j=0;j<nn;j++) 
01171     if (fabs( *(ai+j*mm_max)) >= amax) {amax= fabs( *(ai+j*mm_max)); t=j;}
01172         tmm=t*mm_max;
01173         pivot=*(ai+tmm);
01174         
01175         /* finds contribution to v from this pivot (if not nil) */
01176         
01177         if (amax == 0.)
01178         {
01179             /* Constraint is inconsistent */
01180             
01181             if ( *(bi) < 0.) 
01182         { printf("Constraint %d is inconsistent\n",i+1); return(0.); }
01183             
01184             /* otherwise constraint is redundant */
01185             
01186             printf("Constraint %d is redundant\n",i+1); 
01187         }
01188         else
01189         {
01190             
01191             /* allocate memory */
01192             
01193             ap = (float *) malloc(4*nm1*mm1);
01194             bp = (float *) malloc(4*mm1);
01195             
01196             /* reduce a and b into ap and bp eliminating variable t and constraint i */
01197             
01198             jj=-1;
01199             for (j=0;j<mm;j++)
01200             if (j != i)
01201             {
01202                 jj=jj+1;
01203                 aj=a+j;
01204                 ajt=aj+tmm;
01205                 *(bp+jj)= *(b+j) - *(bi) * *(ajt) / pivot;
01206                 apjj=ap+jj;
01207                 kk=-1;
01208                 for (k=0;k<nn;k++)
01209                 if (k != t)
01210                 {
01211                     kk=kk+1;
01212                     kmm=k*mm_max;
01213                     *(apjj+kk*mm1)= *(aj+kmm)- *(ajt) * *(ai+kmm)/ pivot;
01214                 }
01215             }
01216             
01217             /* add contribution to volume from volume calculated in smaller dimension */
01218             
01219             partialv = 0.;
01220             if (*(bi) != 0.) partialv=cvolumef(ap,bp, &mm1, &nm1, &mm1, &nm1);
01221             if(partialv == -1.0)return(-1.0);
01222             v=v+ *(bi)/amax*partialv/nn; 
01223             
01224             free(ap);
01225             free(bp);
01226         }
01227     }
01228     
01229     return(v);
01230     
01231 }
01232 /*--------------------------------------------------------------
01233 
01234 ROUTINE: volume
01235 
01236 A dummy routine used to call cvolume from a fortran routine 
01237 
01238 ----------------------------------------------------------------*/
01239 
01240 void volumef_ (a,b,m,n,mmax,nmax,volume)
01241 
01242 int   *n, *m, *mmax, *nmax;
01243 float *a, *b;
01244 float *volume; 
01245 
01246 {
01247     
01248     *volume=cvolumef(a,b,m,n,mmax,nmax);
01249     
01250 }
01251 /*--------------------------------------------------------------
01252 
01253 ROUTINE: cdvdaf
01254 
01255 This routine calculates the derivative with respect to a(0,tdim)
01256 of the `volume' in dimension n of the region bounded by a set 
01257 of m linear inequality constraints of the form A x <= b, where 
01258 a has m rows and n columns and is given by a(m,n), b is the 
01259 n-vector and is contained in b(m). The derivative expression
01260 is recursive and derived from the formula of Lasserre (1983). 
01261 
01262 Redundant constraints are allowed and a warning is issued if any are 
01263 encountered. If the inequality constraints are inconsistent then the 
01264 derivative is returned as zero.  If any constraint is orthogonal to 
01265 the component a(0,idim) then the reduction can only take place onto 
01266 variable idim. A special case is used to handle this which involves
01267 no further recursive calls.
01268 
01269 If the original polyhedron is unbounded then a warning is issued
01270 and the derivative is return as zero. 
01271 
01272 Note: This code takes advantage of the fact that during recursive calls
01273 constraint 0 does not change its position in the list of remaining 
01274 constraints (if it has not been eliminated), i.e. it is always the 
01275 first constraint.  This would not be the case if the algorithm were 
01276 adapated to deal with other constraints, i.e. evaluate dvda_i,j 
01277 where i .ne. 0.
01278 
01279 This is a faster version which does not continue down the
01280 recursive tree if if encounters b(j) = 0 for any j. This
01281 restricts its use to cases where the origin does not
01282 pass through any inconsistent constraints. Also redundant
01283 constraints that pass through the origin will not be detected.
01284 
01285 Calls itself cvolumeb, and cvolumebj.
01286 
01287 Malcolm Sambridge, March 1995.
01288 
01289 --------------------------------------------------------------*/
01290 
01291 float cdvdaf (a,b,m,n,mmax,nmax,tdim,temp,jval,code)
01292 
01293 int   *n, *m, *nmax, *mmax, *tdim, *jval, *code;
01294 float *a, *b, *temp; 
01295 
01296 {
01297     
01298     float v,amax,pivot;
01299     int   i,j,t,k,l;
01300     int   jj,kk;
01301     float   *ai, *aj, *ajt, *apjj, *bi;
01302     int   kmm,tmm;
01303     int   nn,mm,nn_max,mm_max,ttdim;
01304     int   nm1,mm1;
01305     int   lmin,lmax, jjval, kval;
01306     float  *ap, *bp, *ttemp;
01307     int   firstmin,firstmax;
01308     float bmin,bmax,bb;
01309     float deriv, junk, vol, dvdb, dbda;
01310     int   special, opt;
01311     
01312     nn= *n;
01313     mm= *m;
01314     nn_max= *nmax;
01315     mm_max= *mmax;
01316     
01317     /* one-dimensional case (full reduction) */
01318     
01319     *code = 0;
01320     
01321     if (nn == 1) 
01322     {
01323         firstmin=0;
01324         firstmax=0;
01325         lmax=0;
01326         lmin=0;
01327         for (l=0;l<mm;l++)
01328         {
01329             if ( *(a+l) > 0.) 
01330             {
01331                 bb= *(b+l)/ *(a+l);
01332             if (firstmin==0) {firstmin=1;bmin=bb;lmin=l;}
01333             else if (bb<bmin) {bmin=bb;lmin=l;}
01334             }
01335             else if ( *(a+l) < 0.)
01336             {
01337                 bb= *(b+l)/ *(a+l);
01338             if (firstmax==0) {firstmax=1;bmax=bb;lmax=l;}
01339             else if (bb>bmax) {bmax=bb;lmax=l;}
01340             }
01341             else if ( *(b+l) < 0.) 
01342             {
01343                 /* Constraint is inconsistent.
01344                 Set derivative to zero. */
01345                 
01346                 printf("Inconsistent constraints found after reduction to n = 1 \n");
01347                 *code = 1;
01348                 return(0.);
01349             }
01350         }
01351         v=0.;
01352         if (firstmin*firstmax == 1) v=bmin-bmax;
01353         else 
01354         {
01355             /*     printf("Volume is unbounded; derivative returned is zero\n"); */
01356             *code = -1;
01357             return(0.);
01358         }
01359         
01360         if (v<0.) return(0.);       
01361         
01362         if(*jval == 1)           /* Constraint 0 has not yet been encountered */
01363         {
01364             if (lmin == 0) deriv = -bmin/ *a;
01365             else if (lmax == 0) deriv = bmax/ *a;
01366             else deriv = 0.;
01367             return(deriv); 
01368         }
01369         else if(*jval == 0)     /* Constraint 0  has already been encountered */
01370         {
01371             deriv =  ( *(temp+lmax) * (bmax/ *(a+lmax)) ) -
01372             ( *(temp+lmin) * (bmin/ *(a+lmin)) );
01373             return(deriv); 
01374         }
01375     }
01376     nm1=nn-1;
01377     mm1=mm-1;
01378     v=0.;
01379     
01380     /*  perform main loop over constraints */
01381     
01382     for (i=0;i<mm;i++)
01383     {
01384         ai=a+i;
01385         bi=b+i;
01386         ttdim = *tdim;
01387         special = 0;
01388         
01389         /* find largest pivot */
01390         
01391         amax=0.;
01392         t = 0;
01393         for (j=0;j<nn;j++) 
01394         if (fabs( *(ai+j*mm_max)) >= amax && j != ttdim) 
01395     {amax= fabs( *(ai+j*mm_max)); t=j;}
01396         
01397         /* finds contribution to v from 
01398         this pivot (if not nil) */
01399         
01400         if (amax == 0.)
01401         {
01402             if(*(ai + ttdim * mm_max) == 0.0) 
01403             { 
01404                 /* Constraint is inconsistent */
01405                 if ( *(bi) < 0.)
01406                 {
01407                     printf("Constraint %d is inconsistent\n",i+1);
01408                     *code = 1;
01409                     return(0.);
01410                 }
01411                 
01412                 /* otherwise constraint is redundant */
01413                 
01414                 printf("Constraint %d is redundant\n",i+1); 
01415             }
01416             else
01417             {
01418                 /* if projection can only be peformed
01419                 on dimension tdim then activate 
01420                 special case */ 
01421                 
01422                 special = 1;
01423                 t = ttdim;
01424                 amax = fabs(*(ai+t * mm_max));
01425                 
01426             }
01427         }
01428         
01429         tmm=t*mm_max;
01430         pivot=*(ai+tmm);
01431         
01432         if(t < ttdim) ttdim = ttdim -1;
01433         
01434         
01435         if (amax != 0)
01436         {
01437             
01438             /* determine if constraint 0 has been encountered */
01439             
01440             kval = 0;
01441             if ( i == 0 && *jval == 1)
01442             {
01443                 /* This is the first encounter of
01444                 constraint 0 on this path so we 
01445                 allocate memory and store parameters 
01446                 to be used when n = 1 */
01447                 
01448                 if (special == 0) 
01449                 {
01450                     ttemp = (float *) malloc(4*mm1);
01451                     for (j=0;j<mm1;j++) *(ttemp+j) = - *(a+j+1+tmm)/pivot;
01452                     kval = 1;
01453                 }
01454                 jjval = 0;
01455             }
01456             else if (*jval == 0)             /* Constraint 0 has already been
01457             encountered */
01458             {
01459                 jjval = 0;  
01460                 /* perform recursive update of component
01461                 derivative array temp. This eliminates 
01462                 row i and copies into a new vector */
01463                 
01464                 if(special == 0)
01465                 {
01466                     ttemp = (float *) malloc(4*mm1);
01467                     for (j=0;j<i;j++) *(ttemp+j) = *(temp+j) 
01468                     -(*(temp+i) * *(a+j+tmm)/pivot);
01469                     for (j=i;j<mm1;j++) *(ttemp+j) = *(temp+j+1)
01470                     -(*(temp+i) * *(a+j+1+tmm)/pivot);
01471                 }
01472             }
01473             else 
01474             {                             /* Constraint 0 has not yet been 
01475                 encountered */
01476                 jjval = 1;                
01477             }
01478             
01479             
01480             /* allocate memory */
01481             
01482             ap = (float *) malloc(4*nm1*mm1);
01483             bp = (float *) malloc(4*mm1);
01484             
01485             /* reduce a and b into ap and bp eliminating variable t and constraint i */
01486             
01487             jj=-1;
01488             for (j=0;j<mm;j++)
01489             if (j != i)
01490             {
01491                 jj=jj+1;
01492                 aj=a+j;
01493                 ajt=aj+tmm;
01494                 *(bp+jj)= *(b+j) - *(bi) * *(ajt) / pivot;
01495                 apjj=ap+jj;
01496                 kk=-1;
01497                 for (k=0;k<nn;k++)
01498                 if (k != t)
01499                 {
01500                     kk=kk+1;
01501                     kmm=k*mm_max;
01502                     *(apjj+kk*mm1)= *(aj+kmm)- *(ajt) * *(ai+kmm)/ pivot;
01503                 }
01504             }
01505             
01506             /* add contribution to derivative from that calculated in smaller dimension */
01507             
01508             /* Normal case method */
01509             if(special == 0)
01510             {
01511                 deriv = 0.;
01512                 if (*(bi) != 0.) 
01513             {deriv=cdvdaf(ap,bp, &mm1, &nm1, &mm1, &nm1, &ttdim, ttemp, &jjval,code);}
01514                 v=v+ *(bi)/amax*deriv/nn;
01515                 if (kval == 1 || *jval == 0) free(ttemp);
01516                 if(*code != 0)return (0.);
01517                 
01518             }
01519             else                    /* Use special case method */
01520             {
01521                 if( *jval == 1)                    
01522                 {                     
01523                     if(i == 0)                      /* This is constraint 0 */
01524                     {
01525                         deriv = 0.; 
01526                         vol = 0.; 
01527                         dvdb = 0.;
01528                         for (j=1;j<mm;j++) 
01529                         {
01530                             k = j - 1;
01531                             junk=cvolumebj(ap,bp,&mm1,&nm1,&mm1,&nm1,k,&dvdb);
01532                             if(junk == -1.)
01533                             {
01534                                 *code = -1;
01535                                 return(0.);
01536                             }
01537                             deriv = deriv + dvdb * *(a + j + tmm) ;
01538                             vol = vol + junk;
01539                         }
01540                         if(nm1 == 1)vol = junk;
01541                         deriv = *(bi) * deriv/pivot;
01542                         deriv = (deriv - vol) /pivot; 
01543                         v=v+ *(bi)/amax*deriv/nn;
01544                         
01545                     }
01546                     else                /* Constraint 0 not yet encountered */
01547                     {
01548                         opt = 2;
01549                         junk=cvolumeb(ap,bp,&mm1,&nm1,&mm1,&nm1,&opt,&deriv);
01550                         if(junk == -1.)
01551                         {
01552                             *code = -1;
01553                             return(0.);
01554                         }
01555                         deriv= -(deriv *  *(bi)/pivot);
01556                         v=v+ *(bi)/amax*deriv/nn;
01557                     }
01558                 }
01559                 else if( *jval == 0)               /* Constraint 0 already encountered */
01560                 {
01561                     vol = 0.; 
01562                     deriv = 0.; 
01563                     for (j=0;j<i;j++) 
01564                     {
01565                         junk=cvolumebj(ap,bp,&mm1,&nm1,&mm1,&nm1,j,&dvdb);
01566                         if(junk == -1.)
01567                         {
01568                             *code = -1;
01569                             return(0.);
01570                         }
01571                         dbda = (*(a+j+tmm) * *(temp+i)/pivot) - *(temp+j);
01572                         deriv = deriv + dvdb*dbda;
01573                         vol = vol + junk;
01574                     }
01575                     for (j=i+1;j<mm;j++) 
01576                     {
01577                         k = j - 1;
01578                         junk=cvolumebj(ap,bp,&mm1,&nm1,&mm1,&nm1,k,&dvdb);
01579                         if(junk == -1.)
01580                         {
01581                             *code = -1;
01582                             return(0.);
01583                         }
01584                         dbda = (*(a+j+tmm) * *(temp+i)/pivot) - *(temp+j);
01585                         deriv = deriv + dvdb*dbda;
01586                         vol = vol + junk;
01587                     }
01588                     if(nm1 == 1)vol = junk;
01589                     deriv = (*(bi) * deriv - *(temp+i)*vol)/pivot;
01590                     v=v+ *(bi)/amax*deriv/nn;
01591                 }
01592             }
01593             
01594             free(ap);
01595             free(bp);
01596         }
01597     }
01598     
01599     return(v);
01600     
01601 }
01602 
01603 /*--------------------------------------------------------------
01604 
01605 ROUTINE: dvda
01606 
01607 A dummy routine used to call cdvda from a fortran routine 
01608 
01609 ----------------------------------------------------------------*/
01610 
01611 
01612 void dvdaf_ (a,b,m,n,mmax,nmax,idim,dvda,code)
01613 
01614 int   *n, *m, *mmax, *nmax, *idim, *code;
01615 float *a, *b;
01616 float *dvda; 
01617 
01618 {
01619     int    jval, tdim;
01620     // Feng VPAC, Oct. 3, 2002---initialize temp to zero(?)
01621     
01622     float *temp=0;
01623     
01624     
01625     jval = 1;
01626     tdim = *idim - 1;
01627     
01628     *dvda=cdvdaf(a,b,m,n,mmax,nmax,&tdim,temp,&jval,code);
01629     
01630     free(temp);
01631 }
01632 
01633 
01634 /*--------------------------------------------------------------
01635 
01636 ROUTINE: cdvda_debug
01637 
01638 This routine calculates the derivative with respect to a(0,tdim)
01639 of the `volume' in dimension n of the region bounded by a set 
01640 of m linear inequality constraints of the form A x <= b, where 
01641 a has m rows and n columns and is given by a(m,n), b is the 
01642 n-vector and is contained in b(m). The derivative expression
01643 is recursive and derived from the formula of Lasserre (1983). 
01644 
01645 Redundant constraints are allowed and a warning is issued if any are 
01646 encountered. If the inequality constraints are inconsistent then the 
01647 derivative is returned as zero.  If any constraint is orthogonal to 
01648 the component a(0,idim) then the reduction can only take place onto 
01649 variable idim. A special case is used to handle this which involves
01650 no further recursive calls.
01651 
01652 If the original polyhedron is unbounded then a warning is issued
01653 and the derivative is return as zero. 
01654 
01655 Note: This code takes advantage of the fact that during recursive calls
01656 constraint 0 does not change its position in the list of remaining 
01657 constraints (if it has not been eliminated), i.e. it is always the 
01658 first constraint.  This would not be the case if the algorithm were 
01659 adapated to deal with other constraints, i.e. evaluate dvda_i,j 
01660 where i .ne. 0.
01661 
01662 Calls itself cvolumeb, and cvolumebj.
01663 
01664 This is a version of cdvda used for debugging and contains extra
01665 code and print statements.
01666 
01667 Malcolm Sambridge, March 1995.
01668 
01669 --------------------------------------------------------------*/
01670 
01671 float cdvda_debug (a,b,m,n,tdim,temp,jval,code)
01672 
01673 int   *n, *m, *tdim, *jval, *code;
01674 float *a, *b, *temp; 
01675 
01676 {
01677     
01678     float v,amax,pivot;
01679     int   i,j,t,k,l;
01680     int   jj,kk;
01681     float   *ai, *aj, *ajt, *apjj, *bi;
01682     int   kmm,tmm;
01683     int   nn,mm, ttdim;
01684     int   nm1,mm1;
01685     int   lmin,lmax, jjval, kval;
01686     float  *ap, *bp, *ttemp;
01687     int   firstmin,firstmax;
01688     float bmin,bmax,bb;
01689     float deriv, junk, vol, dvdb, dbda;
01690     int   special, opt;
01691     float volp, volm, da, FD, FDa;
01692     
01693     
01694     nn= *n;
01695     mm= *m;
01696     
01697     /* write out a and b (debug) */
01698     
01699     printf("\n"); 
01700     printf(" Current matrix A and vector b \n");
01701     for (j=0;j<mm;j++)
01702     {
01703         aj=a+j;
01704         for (k=0;k<nn;k++) 
01705     {ajt=aj+k * mm; printf(" a %d %d = %f",j,k,*(ajt));}
01706         printf(" : b = %f \n", *(b+j));
01707     }
01708     
01709     /* one-dimensional case (full reduction) */
01710     
01711     *code = 0;
01712     
01713     if (nn == 1) 
01714     {
01715         printf("One dimension left \n");
01716         firstmin=0;
01717         firstmax=0;
01718         lmax=0;
01719         lmin=0;
01720         for (l=0;l<mm;l++)
01721         {
01722             if ( *(a+l) > 0.) 
01723             {
01724                 bb= *(b+l)/ *(a+l);
01725             if (firstmin==0) {firstmin=1;bmin=bb;lmin=l;}
01726             else if (bb<bmin) {bmin=bb;lmin=l;}
01727             }
01728             else if ( *(a+l) < 0.)
01729             {
01730                 bb= *(b+l)/ *(a+l);
01731             if (firstmax==0) {firstmax=1;bmax=bb;lmax=l;}
01732             else if (bb>bmax) {bmax=bb;lmax=l;}
01733             }
01734             else if ( *(b+l) < 0.) 
01735             {
01736                 /* Constraint is inconsistent.
01737                 Set derivative to zero. */
01738                 
01739                 printf("Inconsistent constraints found after reduction to n = 1 \n");
01740                 *code = 1;
01741                 return(0.);
01742             }
01743         }
01744         v=0.;
01745         if (firstmin*firstmax == 1) v=bmin-bmax;
01746         else 
01747         {
01748             /*     printf("Volume is unbounded; derivative returned is zero\n"); */
01749             *code = -1;
01750             return(0.);
01751         }
01752         
01753         if (v<0.) return(0.);       
01754         
01755         if(*jval == 1)           /* Constraint 0 has not yet been encountered */
01756         {
01757             if (lmin == 0) deriv = -bmin/ *a;
01758             else if (lmax == 0) deriv = bmax/ *a;
01759             else deriv = 0.;
01760             printf(" No projection onto constraint 0 \n");
01761             printf(" lmin = %d lmax = %d bmin = %f bmax = %f alpha = %f \n",
01762             lmin,lmax,bmin,bmax,*a);
01763             printf(" deriv contribution = %f \n",deriv);
01764             return(deriv); 
01765         }
01766         else if(*jval == 0)     /* Constraint 0  has already been encountered */
01767         {
01768             deriv =  ( *(temp+lmax) * (bmax/ *(a+lmax)) ) -
01769             ( *(temp+lmin) * (bmin/ *(a+lmin)) );
01770             printf(" Projection onto constraint 0 has occurred \n");
01771             printf(" lmin = %d bmin = %f alpha = %f temp = %f\n",
01772             lmin,bmin,*(a+lmin),*(temp+lmin));
01773             printf(" lmax = %d bmax = %f alpha = %f temp = %f\n",
01774             lmax,bmax,*(a+lmax),*(temp+lmax));
01775             printf(" deriv contribution = %f \n",deriv);
01776             return(deriv); 
01777         }
01778     }
01779     nm1=nn-1;
01780     mm1=mm-1;
01781     v=0.;
01782     
01783     printf(" About to start main loop n = %d m = %d \n",nn,mm);
01784     
01785     /*  perform main loop over constraints */
01786     
01787     for (i=0;i<mm;i++)
01788     {
01789         ai=a+i;
01790         bi=b+i;
01791         ttdim = *tdim;
01792         special = 0;
01793         
01794         /* find largest pivot */
01795         
01796         amax=0.;
01797         t = 0;
01798         for (j=0;j<nn;j++) 
01799         if (fabs( *(ai+j*mm)) >= amax && j != ttdim) 
01800     {amax= fabs( *(ai+j*mm)); t=j;}
01801         
01802         /* finds contribution to v from 
01803         this pivot (if not nil) */
01804         
01805         if (amax == 0.)
01806         {
01807             if(*(ai + ttdim * mm) == 0.0) 
01808             { 
01809                 /* Constraint is inconsistent */
01810                 if ( *(bi) < 0.)
01811                 {
01812                     printf("Constraint %d is inconsistent\n",i+1);
01813                     *code = 1;
01814                     return(0.);
01815                 }
01816                 
01817                 /* otherwise constraint is redundant */
01818                 
01819                 printf("Constraint %d is redundant\n",i+1); 
01820             }
01821             else
01822             {
01823                 /* if projection can only be peformed
01824                 on dimension tdim then activate 
01825                 special case */ 
01826                 
01827                 printf("Constraint %d can only be projected onto dimension %d\n",i+1,*tdim);
01828                 printf("using special case method \n");
01829                 special = 1;
01830                 t = ttdim;
01831                 amax = fabs(*(ai+t * mm));
01832                 
01833             }
01834         }
01835         
01836         tmm=t*mm;
01837         pivot=*(ai+tmm);
01838         
01839         printf("\n Projection onto constraint %d variable %d k = %d\n",i,t,ttdim);
01840         
01841         if(t < ttdim) ttdim = ttdim -1;
01842         
01843         
01844         if (amax != 0)
01845         {
01846             
01847             /* determine if constraint 0 has been encountered */
01848             
01849             kval = 0;
01850             if ( i == 0 && *jval == 1)
01851             {
01852                 /* This is the first encounter of
01853                 constraint 0 on this path so we 
01854                 allocate memory and store parameters 
01855                 to be used when n = 1 */
01856                 
01857                 if (special == 0) 
01858                 {
01859                     ttemp = (float *) malloc(4*mm);
01860                     for (j=0;j<mm1;j++) *(ttemp+j) = - *(a+j+1+tmm)/pivot;
01861                     kval = 1;
01862                     printf("\n Generating temp n = %d m = %d i = %d \n",nn,mm,i); 
01863                     printf(" ttemp : \n"); 
01864                     for (j=0;j<mm1;j++) printf(" %f ",*(ttemp+j));
01865                     printf("\n"); 
01866                 }
01867                 jjval = 0;
01868             }
01869             else if (*jval == 0)             /* Constraint 0 has already been
01870             encountered */
01871             {
01872                 jjval = 0;  
01873                 /* perform recursive update of component
01874                 derivative array temp. This eliminates 
01875                 row i and copies into a new vector */
01876                 
01877                 if(special == 0)
01878                 {
01879                     ttemp = (float *) malloc(4*mm1);
01880                     for (j=0;j<i;j++) *(ttemp+j) = *(temp+j) 
01881                     -(*(temp+i) * *(a+j+tmm)/pivot);
01882                     for (j=i;j<mm1;j++) *(ttemp+j) = *(temp+j+1)
01883                     -(*(temp+i) * *(a+j+1+tmm)/pivot);
01884                     printf(" Reduced ttemp : \n"); 
01885                     for (j=0;j<mm1;j++) printf(" %f ",*(ttemp+j));
01886                     printf("\n"); 
01887                 }
01888             }
01889             else 
01890             {                             /* Constraint 0 has not yet been 
01891                 encountered */
01892                 jjval = 1;                
01893             }
01894             
01895             
01896             /* allocate memory */
01897             
01898             ap = (float *) malloc(4*nm1*mm1);
01899             bp = (float *) malloc(4*mm1);
01900             
01901             /* reduce a and b into ap and bp eliminating variable t and constraint i */
01902             
01903             jj=-1;
01904             for (j=0;j<mm;j++)
01905             if (j != i)
01906             {
01907                 jj=jj+1;
01908                 aj=a+j;
01909                 ajt=aj+tmm;
01910                 *(bp+jj)= *(b+j) - *(bi) * *(ajt) / pivot;
01911                 apjj=ap+jj;
01912                 kk=-1;
01913                 for (k=0;k<nn;k++)
01914                 if (k != t)
01915                 {
01916                     kk=kk+1;
01917                     kmm=k*mm;
01918                     *(apjj+kk*mm1)= *(aj+kmm)- *(ajt) * *(ai+kmm)/ pivot;
01919                 }
01920             }
01921             
01922             /* add contribution to derivative from that calculated in smaller dimension */
01923             
01924             /* Normal case method */
01925             if(special == 0)
01926             {
01927                 /*
01928                 deriv = 0.;
01929                 if (*(bi) != 0.) 
01930             {deriv=cdvda(ap,bp, &mm1, &nm1, &mm1, &nm1, &ttdim, ttemp, &jjval,code);}
01931                 */
01932                 deriv=cdvda(ap,bp, &mm1, &nm1, &mm1, &nm1, &ttdim, ttemp, &jjval,code);
01933                 v=v+ *(bi)/amax*deriv/nn;
01934                 if (kval == 1 || *jval == 0) free(ttemp);
01935                 if(*code != 0)return (0.);
01936                 
01937             }
01938             else                    /* Use special case method */
01939             {
01940                 if( *jval == 1)                    
01941                 {                     
01942                     if(i == 0)                      /* This is constraint 0 */
01943                     {
01944                         deriv = 0.; 
01945                         vol = 0.; 
01946                         dvdb = 0.;
01947                         printf("\n special case 2: reduction by 0 and variable k\n"); 
01948                         for (j=1;j<mm;j++) 
01949                         {
01950                             k = j - 1;
01951                             junk=cvolumebj(ap,bp,&mm1,&nm1,&mm1,&nm1,k,&dvdb);
01952                             if(junk == -1.)
01953                             {
01954                                 *code = -1;
01955                                 return(0.);
01956                             }
01957                             printf("\n j %d k %d dvdb %f vol %f",j,k,dvdb,junk); 
01958                             deriv = deriv + dvdb * *(a + j + tmm) ;
01959                             vol = vol + junk;
01960                         }
01961                         if(nm1 == 1)vol = junk;
01962                         deriv = *(bi) * deriv/pivot;
01963                         deriv = (deriv - vol) /pivot; 
01964                         v=v+ *(bi)/amax*deriv/nn;
01965                         printf("dcont %f",deriv); 
01966                         printf("\n total volume of face = %f",vol); 
01967                         
01968                         /* start debug section */
01969                         
01970                         FDa = *(bi)/amax*deriv/nn;
01971                         da = 0.01* *(a + *tdim * mm);
01972                         if(da == 0.0)da = 0.01;
01973                         *(a + *tdim * mm) = *(a + *tdim * mm) + da;
01974                         
01975                         jj=-1;
01976                         for (j=0;j<mm;j++)
01977                         if (j != i)
01978                         {
01979                             jj=jj+1;
01980                             aj=a+j;
01981                             ajt=aj+tmm;
01982                             *(bp+jj)= *(b+j) - *(bi) * *(ajt) / pivot;
01983                             apjj=ap+jj;
01984                             kk=-1;
01985                             for (k=0;k<nn;k++)
01986                             if (k != t)
01987                             {
01988                                 kk=kk+1;
01989                                 kmm=k*mm;
01990                                 *(apjj+kk*mm1)= *(aj+kmm)- *(ajt) * *(ai+kmm)/ pivot;
01991                             }
01992                         }
01993                         volp = cvolume(ap,bp,&mm1,&nm1,&mm1,&nm1); 
01994                         amax = fabs(pivot+da);
01995                         volp = (*(bi)/(amax*nn)) * volp;
01996                         *(a + *tdim * mm) = *(a + *tdim * mm) - 2. * da;
01997                         jj=-1;
01998                         for (j=0;j<mm;j++)
01999                         if (j != i)
02000                         {
02001                             jj=jj+1;
02002                             aj=a+j;
02003                             ajt=aj+tmm;
02004                             *(bp+jj)= *(b+j) - *(bi) * *(ajt) / pivot;
02005                             apjj=ap+jj;
02006                             kk=-1;
02007                             for (k=0;k<nn;k++)
02008                             if (k != t)
02009                             {
02010                                 kk=kk+1;
02011                                 kmm=k*mm;
02012                                 *(apjj+kk*mm1)= *(aj+kmm)- *(ajt) * *(ai+kmm)/ pivot;
02013                             }
02014                         }
02015                         volm = cvolume(ap,bp,&mm1,&nm1,&mm1,&nm1); 
02016                         amax = fabs(pivot-da);
02017                         *(a + *tdim * mm) = *(a + *tdim * mm) + da;
02018                         volm = (*(bi)/(amax*nn)) * volm;
02019                         FD = (volp-volm)/(2. * da);
02020                         printf("\n k  %d a = %f da = %f",*tdim,*(a + *tdim * mm),da); 
02021                         printf("\n pivot = %f",pivot); 
02022                         printf("\n volp = %f",volp); 
02023                         printf("\n volm = %f",volm); 
02024                         printf("\n FD estimate = %f",FD); 
02025                         printf("\n analytical estimate = %f\n",FDa); 
02026                         
02027                         /* end debug section */
02028                         
02029                     }
02030                     else                /* Constraint 0 not yet encountered */
02031                     {
02032                         printf("\n special case 1: constraint 0 not yet reached\n"); 
02033                         opt = 2;
02034                         junk=cvolumeb(ap,bp,&mm1,&nm1,&mm1,&nm1,&opt,&deriv);
02035                         if(junk == -1.)
02036                         {
02037                             *code = -1;
02038                             return(0.);
02039                         }
02040                         deriv= -(deriv *  *(bi)/pivot);
02041                         v=v+ *(bi)/amax*deriv/nn;
02042                         printf("\n spec: 0 not yet reached i %d dcont %f",i,deriv); 
02043                         
02044                         /* debug section: reduce system and repeat calculation with finite difference */
02045                         
02046                         FDa = *(bi)/amax*deriv/nn;
02047                         da = 0.01* *(a + *tdim * mm);
02048                         if(da == 0.0)da = 0.01;
02049                         *(a + *tdim * mm) = *(a + *tdim * mm) + da;
02050                         
02051                         /* reduce system */
02052                         
02053                         jj=-1;
02054                         for (j=0;j<mm;j++)
02055                         if (j != i)
02056                         {
02057                             jj=jj+1;
02058                             aj=a+j;
02059                             ajt=aj+tmm;
02060                             *(bp+jj)= *(b+j) - *(bi) * *(ajt) / pivot;
02061                             apjj=ap+jj;
02062                             kk=-1;
02063                             for (k=0;k<nn;k++)
02064                             if (k != t)
02065                             {
02066                                 kk=kk+1;
02067                                 kmm=k*mm;
02068                                 *(apjj+kk*mm1)= *(aj+kmm)- *(ajt) * *(ai+kmm)/ pivot;
02069                             }
02070                         }
02071                         volp = cvolume(ap,bp,&mm1,&nm1,&mm1,&nm1); 
02072                         *(a + *tdim * mm) = *(a + *tdim * mm) - 2. * da;
02073                         
02074                         /* reduce system */
02075                         
02076                         jj=-1;
02077                         for (j=0;j<mm;j++)
02078                         if (j != i)
02079                         {
02080                             jj=jj+1;
02081                             aj=a+j;
02082                             ajt=aj+tmm;
02083                             *(bp+jj)= *(b+j) - *(bi) * *(ajt) / pivot;
02084                             apjj=ap+jj;
02085                             kk=-1;
02086                             for (k=0;k<nn;k++)
02087                             if (k != t)
02088                             {
02089                                 kk=kk+1;
02090                                 kmm=k*mm;
02091                                 *(apjj+kk*mm1)= *(aj+kmm)- *(ajt) * *(ai+kmm)/ pivot;
02092                             }
02093                         }
02094                         volm = cvolume(ap,bp,&mm1,&nm1,&mm1,&nm1); 
02095                         *(a + *tdim * mm) = *(a + *tdim * mm) + da;
02096                         FD = (*(bi)/(amax*nn)) * (volp-volm)/(2. * da);
02097                         printf("\n k  %d a = %f da = %f",*tdim,*(a + *tdim * mm),da); 
02098                         printf("\n volp = %f",volp); 
02099                         printf("\n volm = %f",volm); 
02100                         printf("\n FD estimate = %f",FD); 
02101                         printf("\n analytical estimate = %f\n",FDa); 
02102                         
02103                         /* end debug section */
02104                         
02105                     }
02106                 }
02107                 else if( *jval == 0)               /* Constraint 0 already encountered */
02108                 {
02109                     printf("\n special case 3: constraint 0 already reduced\n"); 
02110                     vol = 0.; 
02111                     deriv = 0.; 
02112                     for (j=0;j<i;j++) 
02113                     {
02114                         junk=cvolumebj(ap,bp,&mm1,&nm1,&mm1,&nm1,j,&dvdb);
02115                         if(junk == -1.)
02116                         {
02117                             *code = -1;
02118                             return(0.);
02119                         }
02120                         dbda = (*(a+j+tmm) * *(temp+i)/pivot) - *(temp+j);
02121                         deriv = deriv + dvdb*dbda;
02122                         vol = vol + junk;
02123                     }
02124                     for (j=i+1;j<mm;j++) 
02125                     {
02126                         k = j - 1;
02127                         junk=cvolumebj(ap,bp,&mm1,&nm1,&mm1,&nm1,k,&dvdb);
02128                         if(junk == -1.)
02129                         {
02130                             *code = -1;
02131                             return(0.);
02132                         }
02133                         dbda = (*(a+j+tmm) * *(temp+i)/pivot) - *(temp+j);
02134                         deriv = deriv + dvdb*dbda;
02135                         vol = vol + junk;
02136                     }
02137                     if(nm1 == 1)vol = junk;
02138                     deriv = (*(bi) * deriv - *(temp+i)*vol)/pivot;
02139                     v=v+ *(bi)/amax*deriv/nn;
02140                     printf("\n spec: 0 already encountered i = %d dcont %f",i,deriv); 
02141                 }
02142             }
02143             
02144             free(ap);
02145             free(bp);
02146         }
02147     }
02148     
02149     return(v);
02150     
02151 }
02152