cone.c
Go to the documentation of this file.
00001 #include "arith.h"
00002 
00003 cone(x,n,p,w)
00004 MATRIX x;
00005 int n,p;
00006 REAL w[MAXRANK][MAXEDGE];
00007 {
00008    MATRIX u,s,v,f;
00009    REAL fugo,fugo_calc();
00010    LINE linetab[MAXEDGE];
00011    register int i,j,en,ii,mi,m,mm,det,start,nn;
00012    register unsigned int planetab,planework;
00013    int ssvdc(),r[MAX];
00014 
00015    /* pre process, f= U1Drr */
00016    planetab=cone_pre(x,n,p,u,s,v,r,f);
00017 
00018    /* ask the edges of the initial cone */
00019    m=0;
00020    for(i=0; m<r[n]; i++){
00021      if((planetab & (1<<i)) != 0){
00022        planework=planetab & ~(1<<i);
00023        det=edge(f,n,planetab,planework,r[n],linetab,m);
00024        face_init(linetab,m,r[n]);
00025        m += 1;
00026      }
00027    }
00028 
00029    /* cut the initial cone by the remaining hyperplane */
00030    start=0;  mi=m;  mm=m;
00031    for(i=0; i<n; i++){
00032      /* check if hyperplane(i) is not used for the initial cone */
00033      if((planetab & (1<<i)) == 0){
00034        planetab= planetab | (1<<i);
00035        for(en=start; en<m; en=linetab[en].next){
00036          fugo= fugo_calc(f,linetab,i,en,r[n]);
00037          /* when edge(en) is outside the hyperplane(i), create new edges */
00038          if(fugo < -EPS){
00039            /* delete edge(en) from the edge table and update the pointers */
00040            start= update_edge(linetab,en,start);
00041            for(ii=0; linetab[en].face[ii] != EMPTY; ii++){
00042              /* nn is the edge number sharing a common face with edge(en) */
00043              nn=linetab[en].face[ii];
00044              /* when edge(nn) is not outside the cone, the index is TRUE */
00045              if(linetab[nn].index == TRUE){
00046                fugo= fugo_calc(f,linetab,i,nn,r[n]);
00047                /* when edge(nn) is not outside the cone, create a new edge */
00048                /* next test is (fugo >= 0.0) if no numerical error exists  */  
00049                if(fugo > -EPS){
00050                  /* planework is the common hyperplane between en and nn, */
00051                  /* and the new hyperplane(i)                             */
00052                  planework = linetab[en].plane & linetab[nn].plane;
00053                  planework = planework | (1<<i);
00054                  det=edge(f,n,planetab,planework,r[n],linetab,mm);
00055                  if(det == TRUE){ /* det is always TRUE in this case */
00056                    /* update the pointers between edge(nn) and the new edge */
00057                    for(j=0; linetab[nn].face[j] != en; j++)
00058                        ;
00059                    linetab[nn].face[j]=mm;
00060                    linetab[mm].face[0]=nn;
00061                    mm += 1;
00062                  }
00063                }
00064              }
00065            }
00066          }
00067        }
00068        m=mm;
00069        /* update the face pointers among new edges on the new hyperplane(i) */
00070        face(linetab,mi,m,r[n],n);
00071      }
00072      mi=m;
00073    }
00074 
00075    /* post process */
00076    en=cone_post(s,v,r,n,linetab,start,m,w);
00077    return(en);
00078 }


euslisp
Author(s): Toshihiro Matsui
autogenerated on Thu Mar 9 2017 04:57:49