volume.c

           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,&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,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,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);

}