source: icGREP/icgrep-devel/cudd-2.5.1/nanotrav/ntrMflow.c

Last change on this file was 4597, checked in by nmedfort, 4 years ago

Upload of the CUDD library.

File size: 55.8 KB
Line 
1/**CFile***********************************************************************
2
3  FileName    [ntrMflow.c]
4
5  PackageName [ntr]
6
7  Synopsis    [Symbolic maxflow algorithm.]
8
9  Description [This file contains the functions that implement the
10  symbolic version of Dinits's maxflow algorithm described in the
11  ICCAD93 paper. The present implementation differs from the algorithm
12  described in the paper in that more than one matching techniques is
13  used. The technique of the paper is the one applied to
14  hourglass-type bilayers here.]
15
16  Author      [Fabio Somenzi, Gary Hachtel]
17
18  Copyright   [Copyright (c) 1995-2012, Regents of the University of Colorado
19
20  All rights reserved.
21
22  Redistribution and use in source and binary forms, with or without
23  modification, are permitted provided that the following conditions
24  are met:
25
26  Redistributions of source code must retain the above copyright
27  notice, this list of conditions and the following disclaimer.
28
29  Redistributions in binary form must reproduce the above copyright
30  notice, this list of conditions and the following disclaimer in the
31  documentation and/or other materials provided with the distribution.
32
33  Neither the name of the University of Colorado nor the names of its
34  contributors may be used to endorse or promote products derived from
35  this software without specific prior written permission.
36
37  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
38  "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
39  LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
40  FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
41  COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
42  INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
43  BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
44  LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
45  CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
46  LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
47  ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
48  POSSIBILITY OF SUCH DAMAGE.]
49
50******************************************************************************/
51
52#include "ntr.h"
53
54/*---------------------------------------------------------------------------*/
55/* Constant declarations                                                     */
56/*---------------------------------------------------------------------------*/
57
58#define MAXPHASE 1000
59#define MAXLAYER 1000
60#define MAXFPIT  100000
61#define MANY_TIMES 3.0
62
63#define PRUNE           /* If defined, enables pruning of E */
64
65/*---------------------------------------------------------------------------*/
66/* Stucture declarations                                                     */
67/*---------------------------------------------------------------------------*/
68
69/*---------------------------------------------------------------------------*/
70/* Type declarations                                                         */
71/*---------------------------------------------------------------------------*/
72
73typedef struct flowStatsStruct {
74    int pr;                     /* level of verbosity */
75    long start_time;            /* cpu time when the covering started */
76    int phases;                 /* number of phases */
77    int layers;                 /* number of layers */
78    int fpit;                   /* number of fixed point iterations */
79} flowStats;
80
81/*---------------------------------------------------------------------------*/
82/* Variable declarations                                                     */
83/*---------------------------------------------------------------------------*/
84
85#ifndef lint
86static char rcsid[] UTIL_UNUSED = "$Id: ntrMflow.c,v 1.8 2012/02/05 01:53:01 fabio Exp fabio $";
87#endif
88
89static DdNode *xcube, *ycube, *zcube;
90
91/*---------------------------------------------------------------------------*/
92/* Macro declarations                                                        */
93/*---------------------------------------------------------------------------*/
94
95/**AutomaticStart*************************************************************/
96
97/*---------------------------------------------------------------------------*/
98/* Static function prototypes                                                */
99/*---------------------------------------------------------------------------*/
100
101static void maximal_pull (DdManager *bdd, int l, DdNode *ty, DdNode **neW, DdNode **U, DdNode *E, DdNode **F, DdNode **x, DdNode **y, DdNode **z, int n, DdNode *pryz, DdNode *prxz, flowStats *stats);
102static void propagate_maximal_flow (DdManager *bdd, int m, DdNode **neW, DdNode **U, DdNode **x, DdNode **y, DdNode **z, int n, DdNode *pryz, DdNode *prxz, flowStats *stats);
103static void trellis (DdManager *bdd, int m, DdNode **neW, DdNode **U, DdNode **x, DdNode **y, DdNode **z, int n, DdNode *pryz, DdNode *prxz, flowStats *stats);
104static void rhombus (DdManager *bdd, int m, DdNode **neW, DdNode **U, DdNode **x, DdNode **y, DdNode **z, int n, DdNode *pryz, DdNode *prxz, flowStats *stats);
105static void hourglass (DdManager *bdd, int m, DdNode **neW, DdNode **U, DdNode **x, DdNode **y, DdNode **z, int n, DdNode *pryz, DdNode *prxz, flowStats *stats);
106static void maximal_push (DdManager *bdd, int l, DdNode **U, DdNode **F, DdNode **x, DdNode **y, DdNode **z, int n, DdNode *pryz, DdNode *prxz, flowStats *stats);
107static void trellisPush (DdManager *bdd, int m, DdNode **U, DdNode **x, DdNode **y, DdNode **z, int n, DdNode *pryz, DdNode *prxz, flowStats *stats);
108static void rhombusPush (DdManager *bdd, int m, DdNode **U, DdNode **x, DdNode **y, DdNode **z, int n, DdNode *pryz, DdNode *prxz, flowStats *stats);
109static void hourglassPush (DdManager *bdd, int m, DdNode **U, DdNode **x, DdNode **y, DdNode **z, int n, DdNode *pryz, DdNode *prxz, flowStats *stats);
110
111/**AutomaticEnd***************************************************************/
112
113/*---------------------------------------------------------------------------*/
114/* Definition of exported functions                                          */
115/*---------------------------------------------------------------------------*/
116
117
118/**Function********************************************************************
119
120  Synopsis    []
121
122  Description [This function implements Dinits's algorithm for (0-1)
123  max flow, using BDDs and a symbolic technique to trace multiple
124  edge-disjoint augmenting paths to complete a phase. The outer
125  forever loop is over phases, and the inner forever loop is to
126  propagate a (not yet) maximal flow of edge-disjoint augmenting paths
127  from one layer to the previous. The subprocedure call implements a
128  least fixed point iteration to compute a (not yet) maximal flow
129  update between layers. At the end of each phase (except the last
130  one) the flow is actually pushed from the source to the sink.
131
132  Data items:
133    <ul>
134    <li> sx(ty) BDD representations of s(t).
135    <li> x(y)   The variables encoding the from(to)-node u(v) of an
136                edge (u,v) in the given digraph.
137    <li> z      Another set of variables.
138    <li> E(x,y) The edge relation.
139    <li> F(x,y) BDD representation of the current flow, initialized to 0
140                for each arc, and updated by +1, -1, or 0 at the
141                end of each phase.
142    <li> Ms Mt  The maximum flow, that is, the cardinality of a minimum cut,
143                measured at the source and at the sink, respectively.
144                The two values should be identical.
145    <li> reached The set of nodes already visited in the BFS of the digraph.
146    <li> fos    fanout of the source node s.
147    <li> fit    fanin of the sink node t.
148    </ul>
149  ]
150
151  SideEffects [The flow realtion F and the cutset relation cut are returned
152  as side effects.]
153
154  SeeAlso     []
155
156******************************************************************************/
157double
158Ntr_maximum01Flow(
159  DdManager * bdd /* manager */,
160  DdNode * sx /* source node */,
161  DdNode * ty /* sink node */,
162  DdNode * E /* edge relation */,
163  DdNode ** F /* flow relation */,
164  DdNode ** cut /* cutset relation */,
165  DdNode ** x /* array of row variables */,
166  DdNode ** y /* array of column variables */,
167  DdNode ** z /* arrays of auxiliary variables */,
168  int  n /* number of variables in each array */,
169  int  pr /* verbosity level */)
170{
171    flowStats stats;
172    DdNode *one, *zero,
173#ifdef PRUNE
174             *EDC,          /* Edge don't care set */
175#endif
176             *reached,      /* states reached through useful edges */
177             *fos, *fit,    /* fanout of source, fanin of sink */
178             *rF, *rB, *tx,
179             *I, *P,
180             *w, *p, *q, *r,/* intemediate results */
181             *pryz, *prxz,  /* priority functions for disjoint path tracing */
182             **neW, **U;    /* arrays of BDDs for flow propagation */
183    int      i, j, l;
184    double   Ms, Mt;
185
186    /* Initialize debugging structure. */
187    stats.pr = pr;
188    stats.start_time = util_cpu_time();
189    stats.phases = 0;
190    stats.layers = 0;
191    stats.fpit = 0;
192
193    /* Allocate arrays for new (just reached vertex sets)
194    ** and U (useful edge sets).
195    */
196    U   = ALLOC(DdNode *, ((unsigned) MAXLAYER));
197    neW = ALLOC(DdNode *, ((unsigned) MAXLAYER));
198
199    one = Cudd_ReadOne(bdd);
200    zero = Cudd_Not(one);
201
202    /* Initialize xcube, ycube, and zcube (for abstractions). */
203    Cudd_Ref(xcube = Cudd_bddComputeCube(bdd,x,NULL,n));
204    Cudd_Ref(ycube = Cudd_bddComputeCube(bdd,y,NULL,n));
205    Cudd_Ref(zcube = Cudd_bddComputeCube(bdd,z,NULL,n));
206
207    /* Build the BDD for the priority functions. */
208    Cudd_Ref(pryz = Cudd_Dxygtdxz(bdd, n, x, y, z));
209    Cudd_Ref(prxz = Cudd_Dxygtdyz(bdd, n, x, y, z));
210    /* Now "randomize" by shuffling the x variables in pryz and the y
211    ** variables in prxz.
212    */
213    Cudd_Ref(p = Cudd_bddAdjPermuteX(bdd,pryz,x,n));
214    Cudd_RecursiveDeref(bdd,pryz);
215    pryz = p;
216    if(pr>2){(void) fprintf(stdout,"pryz");Cudd_PrintDebug(bdd,pryz,n*3,pr);}
217    Cudd_Ref(p = Cudd_bddAdjPermuteX(bdd,prxz,y,n));
218    Cudd_RecursiveDeref(bdd,prxz);
219    prxz = p;
220    if(pr>2){(void) fprintf(stdout,"prxz");Cudd_PrintDebug(bdd,prxz,n*3,pr);}
221
222#ifdef PRUNE
223    /* Build the edge don't care set and prune E. The edge don't care
224    ** set consists of the edges into the source(s), the edges out of the
225    ** sink(s), and the self-loops. These edges cannot contribute to the
226    ** maximum flow.
227    */
228    Cudd_Ref(p = Cudd_bddSwapVariables(bdd, sx, x, y, n));
229    Cudd_Ref(q = Cudd_bddSwapVariables(bdd, ty, x, y, n));
230    Cudd_Ref(r = Cudd_bddOr(bdd, p, q));
231    Cudd_RecursiveDeref(bdd,p);
232    Cudd_RecursiveDeref(bdd,q);
233    Cudd_Ref(p = Cudd_Xeqy(bdd, n, x, y));
234    Cudd_Ref(EDC = Cudd_bddOr(bdd, p, r));
235    Cudd_RecursiveDeref(bdd,p);
236    Cudd_RecursiveDeref(bdd,r);
237    if(pr>2){(void) fprintf(stdout,"EDC");Cudd_PrintDebug(bdd,EDC,n<<1,pr);}
238    Cudd_Ref(p = Cudd_bddAnd(bdd, E, Cudd_Not(EDC)));
239    Cudd_RecursiveDeref(bdd,EDC);
240    if(pr>0){(void) fprintf(stdout,"Given  E");Cudd_PrintDebug(bdd,E,n<<1,pr);}
241    E = p;
242    if(pr>0){(void) fprintf(stdout,"Pruned E");Cudd_PrintDebug(bdd,E,n<<1,pr);}
243#endif
244
245    /* Compute fanin of sink node t: it is an upper bound for the flow. */
246    Cudd_Ref(fit = Cudd_bddAnd(bdd, E, ty));
247    if (pr>2) {
248        /* Compute fanout of source node s. */
249        Cudd_Ref(fos = Cudd_bddAnd(bdd, E, sx));
250        (void) fprintf(stdout,"fos");Cudd_PrintDebug(bdd,fos,n<<1,pr);
251        Cudd_RecursiveDeref(bdd,fos);
252
253        (void) fprintf(stdout,"fit");Cudd_PrintDebug(bdd,fit,n<<1,pr);
254    }
255    /* t(x) is used to check for termination of forward traversal. */
256    Cudd_Ref(tx = Cudd_bddSwapVariables(bdd, ty, x, y, n));
257
258    /* \KW{Procedure}\ \ \PC{maximum\_flow}$(s,t,E(x,y)) */
259    Cudd_Ref(*F = zero);
260
261    for (i = 1; i < MAXPHASE; i++) {
262        stats.phases++;
263        if(pr>0){(void) fprintf(stdout,"## Starting Phase %d at time = %s\n",i,
264                util_print_time(util_cpu_time() - stats.start_time));}
265        /* new[0](x) = s(x);U^0(x,y)=E(x,y)\cdot s(x) \cdot \overline{F(x,y)};
266        ** reached=s; new[1](x)=\exists_xU^0(x,y);
267        */
268        Cudd_Ref(neW[0] = sx);
269        Cudd_Ref(p = Cudd_bddAnd(bdd, sx, Cudd_Not(*F)));
270        Cudd_Ref(U[0] = Cudd_bddAnd(bdd, p, E));
271        Cudd_RecursiveDeref(bdd,p);
272        Cudd_Ref(reached = sx);
273        Cudd_Ref(r = Cudd_bddExistAbstract(bdd, U[0], xcube));
274        Cudd_RecursiveDeref(bdd,U[0]);
275        Cudd_Ref(q = Cudd_bddSwapVariables(bdd, r, x, y, n));
276        Cudd_RecursiveDeref(bdd,r);
277        Cudd_Ref(neW[1] = Cudd_bddAnd(bdd, q, Cudd_Not(reached)));
278        Cudd_RecursiveDeref(bdd,q);
279        if(pr>0) {
280            (void) fprintf(stdout,"neW[1]");Cudd_PrintDebug(bdd,neW[1],n,pr);
281        }
282        for (l = 1; l < MAXLAYER; l++) {
283            if (neW[l] == zero) {       /* flow is maximum */
284                /* cut=reached(x) \cdot E(x,y) \cdot \overline{reached(y)} */
285                Cudd_Ref(p = Cudd_bddAnd(bdd, reached, E));
286                Cudd_Ref(q = Cudd_bddSwapVariables(bdd, reached, x, y, n));
287                Cudd_Ref(*cut = Cudd_bddAnd(bdd, p, Cudd_Not(q)));
288                Cudd_RecursiveDeref(bdd,p);
289                Cudd_RecursiveDeref(bdd,q);
290                Cudd_RecursiveDeref(bdd,reached);
291                for (j = 0; j <= l; j++)
292                    Cudd_RecursiveDeref(bdd,neW[j]);
293                goto endPhases;
294            }
295            /* As soon as we touch one sink node we stop traversal.
296            ** \KW{if} ($t\cdot new^{l} \neq 0$)
297            */
298            if (!Cudd_bddLeq(bdd, tx, Cudd_Not(neW[l]))) {
299                Cudd_RecursiveDeref(bdd,reached);
300                maximal_pull(bdd,l-1,ty,neW,U,E,F,x,y,z,n,pryz,prxz,&stats);
301                goto endLayers;
302            }
303            stats.layers++;
304            if(pr>2){(void) fprintf(stdout,"===== Layer %d =====\n",l);}
305            /* reached(x) = reached(x) + new^l(x) */
306            Cudd_Ref(w = Cudd_bddOr(bdd, reached, neW[l]));
307            Cudd_RecursiveDeref(bdd,reached);
308            reached = w;
309            /* I(y) = \exists_x((E(x,y) \cdot \overline{F(x,y)})
310            **        \cdot new^l(x))
311            */
312            Cudd_Ref(p = Cudd_bddAnd(bdd, E, Cudd_Not(*F)));
313            Cudd_Ref(I = Cudd_bddAndAbstract(bdd, p, neW[l], xcube));
314            if(pr>3){(void) fprintf(stdout,"I");Cudd_PrintDebug(bdd,I,n,pr);}
315            Cudd_RecursiveDeref(bdd,p);
316            /* rF(x)= I(x) \cdot \overline{reached(x)}) */
317            Cudd_Ref(p = Cudd_bddSwapVariables(bdd, I, x, y, n));
318            Cudd_RecursiveDeref(bdd,I);
319            Cudd_Ref(rF = Cudd_bddAnd(bdd, p, Cudd_Not(reached)));
320            Cudd_RecursiveDeref(bdd,p);
321            if(pr>2){(void) fprintf(stdout,"rF");Cudd_PrintDebug(bdd,rF,n,pr);}
322            /* P(x) = \exists_{y}(F(x,y) \cdot new^l(y)) */
323            Cudd_Ref(p = Cudd_bddSwapVariables(bdd, neW[l], x, y, n));
324            Cudd_Ref(P = Cudd_bddAndAbstract(bdd, *F, p, ycube));
325            Cudd_RecursiveDeref(bdd,p);
326            /* rB(x) = P(x) \cdot \overline{reached(x)}) */
327            Cudd_Ref(rB = Cudd_bddAnd(bdd, P, Cudd_Not(reached)));
328            Cudd_RecursiveDeref(bdd,P);
329            if(pr>2){(void) fprintf(stdout,"rB");Cudd_PrintDebug(bdd,rB,n,pr);}
330            /* new^{l+1}(x) = rF(x) + rB(x) */
331            Cudd_Ref(neW[l+1] = Cudd_bddOr(bdd, rF, rB));
332            Cudd_RecursiveDeref(bdd,rB);
333            Cudd_RecursiveDeref(bdd,rF);
334            if(pr>0) {
335                (void) fprintf(stdout,"new[%d]",l+1);
336                Cudd_PrintDebug(bdd,neW[l+1],n,pr);
337            }
338        } /* start next layer */
339        if (l==MAXLAYER) (void) fprintf(stderr,"ERROR! MAXLAYER exceeded.\n");
340        exit(3);
341endLayers:
342        maximal_push(bdd, l-1, U, F, x, y, z, n, pryz, prxz, &stats);
343        for (j = 0; j < l; j++) {
344            Cudd_RecursiveDeref(bdd,U[j]);
345            Cudd_RecursiveDeref(bdd,neW[j]);
346        }
347        Cudd_RecursiveDeref(bdd,neW[l]);
348        if (pr > 0) {
349            Cudd_Ref(p = Cudd_bddAnd(bdd, sx, *F));
350            Ms=Cudd_CountMinterm(bdd, p, n<<1);
351            (void) fprintf(stdout,"# Flow out of s: %g\n", Ms);
352            Cudd_RecursiveDeref(bdd,p);
353        }
354        if (Cudd_bddLeq(bdd, fit, *F)) {
355            Cudd_Ref(*cut = fit);
356            goto endPhases;
357        }
358    } /* start next phase */
359    if (i == MAXPHASE) (void) fprintf(stderr,"ERROR! MAXPHASE exceeded.\n");
360
361    /* Last phase is completed --- print flow results. */
362endPhases:
363    Cudd_RecursiveDeref(bdd,tx);
364
365    Cudd_Ref(q = Cudd_bddAnd(bdd, *F, sx));
366    Ms = Cudd_CountMinterm(bdd, q, n<<1);
367    Cudd_RecursiveDeref(bdd,q);
368
369    Cudd_Ref(q = Cudd_bddAnd(bdd, *F, ty));
370    Mt = Cudd_CountMinterm(bdd, q, n<<1);
371    Cudd_RecursiveDeref(bdd,q);
372
373    if (pr>1) (void) fprintf(stdout,"Mt= %g, Ms= %g \n", Mt, Ms);
374
375    if (pr>3){(void) fprintf(stdout,"pryz");Cudd_PrintDebug(bdd,pryz,n*3,pr);}
376    if (pr>3){(void) fprintf(stdout,"prxz");Cudd_PrintDebug(bdd,prxz,n*3,pr);}
377
378    if (pr>0) {
379        (void) fprintf(stdout,"#### Stats for maximum_flow ####\n");
380        (void) fprintf(stdout,"%d variables %d of which x[i]\n", Cudd_ReadSize(bdd), n);
381        (void) fprintf(stdout,"time       = %s\n",
382                util_print_time(util_cpu_time() - stats.start_time));
383        (void) fprintf(stdout,"phases     = %d\n", stats.phases);
384        (void) fprintf(stdout,"layers     = %d\n", stats.layers);
385        (void) fprintf(stdout,"FP iter.   = %d\n", stats.fpit);
386    }
387
388    Cudd_RecursiveDeref(bdd,fit);
389    Cudd_RecursiveDeref(bdd,pryz);
390    Cudd_RecursiveDeref(bdd,prxz);
391    Cudd_RecursiveDeref(bdd,xcube);
392    Cudd_RecursiveDeref(bdd,ycube);
393    Cudd_RecursiveDeref(bdd,zcube);
394#ifdef PRUNE
395    Cudd_RecursiveDeref(bdd,E);
396#endif
397
398    FREE(U);
399    FREE(neW);
400
401    return(Ms);
402
403} /* end of Ntr_maximum01Flow */
404
405
406/*---------------------------------------------------------------------------*/
407/* Definition of internal functions                                          */
408/*---------------------------------------------------------------------------*/
409
410/*---------------------------------------------------------------------------*/
411/* Definition of static functions                                            */
412/*---------------------------------------------------------------------------*/
413
414
415/**Function********************************************************************
416
417  Synopsis    [Selects set of edge-disjoint paths from layered network.]
418
419  Description [Selects set of edge-disjoint paths from layered
420  network.  maximal_pull is called when the BFS constructing the
421  layered graph reaches a sink. At this point the new sets of the
422  BFS have been formed, and we know every vertex in these sets is
423  reachable from the source vertex (or vertices) s. The new sets are
424  used to compute the set of useful edges exiting each layer to the
425  right. In each layer, propagate_maximal_flow is called to select a
426  maximal subset of these useful edges. This subset is then used to
427  prune new and U.]
428
429  SideEffects [None]
430
431  SeeAlso     [maximal_push]
432
433******************************************************************************/
434static void
435maximal_pull(
436  DdManager * bdd /* manager */,
437  int  l /* depth of layered network for current phase */,
438  DdNode * ty /* sink node */,
439  DdNode ** neW /* array of BFS layers */,
440  DdNode ** U /* array of useful edges */,
441  DdNode * E /* edge relation */,
442  DdNode ** F /* flow relation */,
443  DdNode ** x,
444  DdNode ** y,
445  DdNode ** z /* arrays of variables for rows and columns */,
446  int  n /* number of x variables */,
447  DdNode * pryz,
448  DdNode * prxz /* priority functions */,
449  flowStats * stats)
450{
451    DdNode *p, *q, *r,
452           *UF, *UB;
453    int    m,
454           pr;          /* Print control */
455
456    pr = stats->pr;
457
458    /* The useful edges of the last layer are all the empty edges into
459    ** the sink(s) from new[l].
460    ** U^{l}(x,y) = t(y)\cdot \overline{F(x,y)}\cdot E(x,y)\cdot new^{l}(x)
461    */
462    Cudd_Ref(p = Cudd_bddAnd(bdd, E, Cudd_Not(*F)));
463    Cudd_Ref(q = Cudd_bddAnd(bdd, neW[l], p));
464    Cudd_RecursiveDeref(bdd,p);
465    Cudd_Ref(U[l] = Cudd_bddAnd(bdd, ty, q));
466    Cudd_RecursiveDeref(bdd,q);
467    if(pr>1){(void) fprintf(stdout,"U[%d]",l);Cudd_PrintDebug(bdd,U[l],n<<1,pr);}
468    /* Eliminate from new[l] the states with no paths to the sink(s).
469    ** new^{l}(x)=\exists_y U^{l}(x,y)
470    */
471    Cudd_RecursiveDeref(bdd,neW[l]);
472    Cudd_Ref(neW[l] = Cudd_bddExistAbstract(bdd, U[l], ycube));
473
474    for (m = l; m >= 1; m--) {
475        /* Find usable backward edges from level m-1 to level m.
476        ** UB(x,y) = new^{m}(x) \cdot F(x,y) \cdot new^{m-1}(y)
477        */
478        Cudd_Ref(r = Cudd_bddSwapVariables(bdd, neW[m-1], x, y, n));
479        Cudd_Ref(q = Cudd_bddAnd(bdd, r, *F));
480        Cudd_RecursiveDeref(bdd,r);
481        Cudd_Ref(UB = Cudd_bddAnd(bdd, neW[m], q));
482        Cudd_RecursiveDeref(bdd,q);
483        if(pr>2){(void) fprintf(stdout,"UB");Cudd_PrintDebug(bdd,UB,n<<1,pr);}
484        /* Find usable forward edges.
485        ** UF(x,y) = new^{m}(y) \cdot \overline{F(x,y)} \cdot E(x,y)
486        ** \cdot new^{m-1}(x)
487        */
488        Cudd_Ref(p = Cudd_bddAnd(bdd, E, Cudd_Not(*F)));
489        Cudd_Ref(q = Cudd_bddAnd(bdd, neW[m-1], p));
490        Cudd_RecursiveDeref(bdd,p);
491        Cudd_Ref(r = Cudd_bddSwapVariables(bdd, neW[m], x, y, n));
492        Cudd_Ref(UF = Cudd_bddAnd(bdd, r, q));
493        Cudd_RecursiveDeref(bdd,q);
494        Cudd_RecursiveDeref(bdd,r);
495        if(pr>2){(void) fprintf(stdout,"UF");Cudd_PrintDebug(bdd,UF,n<<1,pr);}
496        /* U^{m-1}(x,y) = UB(y,x) + UF(x,y) */
497        Cudd_Ref(r = Cudd_bddSwapVariables(bdd, UB, x, y, n));
498        Cudd_RecursiveDeref(bdd,UB);
499        Cudd_Ref(U[m-1] = Cudd_bddOr(bdd, UF, r));
500        Cudd_RecursiveDeref(bdd,UF);
501        Cudd_RecursiveDeref(bdd,r);
502        if(pr>2){(void)fprintf(stdout,"U[%d]",m-1);
503                Cudd_PrintDebug(bdd,U[m-1],n<<1,pr);}
504        /* new[m-1](x) = \exists_{y}U^{m-1}(x,y) */
505        Cudd_RecursiveDeref(bdd,neW[m-1]);
506        Cudd_Ref(neW[m-1] = Cudd_bddExistAbstract(bdd, U[m-1], ycube));
507        /* Compute maximal disjoint interlayer edge set. */
508        propagate_maximal_flow(bdd, m, neW, U, x, y, z, n, pryz, prxz, stats);
509        if(pr>0) {
510            (void)fprintf(stdout,"U[%d]",m-1);
511            Cudd_PrintDebug(bdd,U[m-1],n<<1,pr);
512        }
513    } /* prune next layer */
514
515    return;
516
517} /* end of maximal_pull */
518
519
520/**Function********************************************************************
521
522  Synopsis    [Pulls flow though a layer.]
523
524  Description [Pulls flow though a layer. propagate_maximal_flow only
525  affects U[m-1] and new[m-1].  At the end of this function, the edges
526  in U[m] are guaranteed to drain all the flow supplied by the edges
527  in U[m-1]. This effect is obtained by pruning U[m-1]. After the
528  pruned U[m-1] is computed, new[m-1] is updated to keep track of what
529  nodes are still useful.
530
531  The pruning is performed without actually measuring the in-potential
532  and the out-potential of each node. Instead, pairs of nodes from U[m-1]
533  and U[m] are matched. To avoid counting, the procedure computes sets of
534  paths that connect layer m-1 to layer m+1 and are edge-disjoint.
535
536  Two paths from layer m-1 to layer m+1 are disjoint if they have distinct
537  end-point or if they have distinct middle points. What condition to
538  enforce depends on the "shape" of the layers.]
539
540  SideEffects [Changes U[m-1] and new[m-1]]
541
542  SeeAlso     [trellis rhombus hourglass]
543
544******************************************************************************/
545static void
546propagate_maximal_flow(
547  DdManager * bdd,
548  int  m /* center of current bilayer */,
549  DdNode ** neW /* array of reachable or useful nodes */,
550  DdNode ** U /* array of usable or useful edges */,
551  DdNode ** x,
552  DdNode ** y,
553  DdNode ** z /* arrays of variables for rows and columns */,
554  int  n /* number of x variables */,
555  DdNode * pryz,
556  DdNode * prxz /* priority functions */,
557  flowStats * stats)
558{
559    DdNode *rN;
560    double   mtl, mtc, mtr;     /* minterms for left, center, right levels */
561    int      pr;                /* print control */
562
563    pr = stats->pr;
564
565    mtl = Cudd_CountMinterm(bdd, neW[m-1], n);
566    mtc = Cudd_CountMinterm(bdd, neW[m], n);
567
568    /* rN(y) = \exists_x U^{m}(x,y) */
569    Cudd_Ref(rN = Cudd_bddExistAbstract(bdd, U[m], xcube));
570    mtr = Cudd_CountMinterm(bdd, rN, n);
571    Cudd_RecursiveDeref(bdd,rN);
572
573    if (pr > 0) {
574        (void) fprintf(stdout, "layer = %d mtl = %g  mtc = %g  mtr = %g",
575                       m, mtl, mtc, mtr);
576    }
577
578    if ((mtc > MANY_TIMES * mtl) || (mtc > MANY_TIMES * mtr)) {
579        if (pr>0) (void) fprintf(stdout, " R\n");
580        rhombus(bdd, m, neW, U, x, y, z, n, pryz, prxz, stats);
581    } else if (mtr > MANY_TIMES * mtc) {
582        if (pr>0) (void) fprintf(stdout, " H\n");
583        hourglass(bdd, m, neW, U, x, y, z, n, pryz, prxz, stats);
584    } else {
585        if (pr>0) (void) fprintf(stdout, " T\n");
586        trellis(bdd, m, neW, U, x, y, z, n, pryz, prxz, stats);
587    }
588    return;
589
590} /* end of propagate_maximal_flow */
591
592
593/**Function********************************************************************
594
595  Synopsis    [Selects edges from a trellis-type bilayer.]
596
597  Description [Selects edges from a trellis-type bilayer. Used to pull flow.
598  First a matching is found in the left layer. Then the paths are extended
599  (if possible) through the right layer. This process is repeated until a
600  maximal flow is found.]
601
602  SideEffects [None]
603
604  SeeAlso     [rhombus hourglass trellisPush]
605
606******************************************************************************/
607static void
608trellis(
609  DdManager * bdd,
610  int  m /* center level of current bilayer */,
611  DdNode ** neW /* array of node levels */,
612  DdNode ** U /* array of edge layers */,
613  DdNode ** x,
614  DdNode ** y,
615  DdNode ** z /* arrays of variables for rows and columns */,
616  int  n /* number of x variables */,
617  DdNode * pryz,
618  DdNode * prxz /* priority functions */,
619  flowStats * stats)
620{
621    DdNode *one, *zero,
622             *p, *q, *r,
623             *Uin,              /* edges to be matched from U[m-1] */
624             *Uout,             /* edges to be matched from U[m] */
625             *P,
626             *LU, *RU,          /* left-unique and right-unique sets of edges */
627             *D,
628             *Ml, *Mr;          /* nodes matched from left and right */
629    int      k,
630             pr;                /* print control */
631
632    pr = stats->pr;
633
634    one = Cudd_ReadOne(bdd);
635    zero = Cudd_Not(one);
636
637    /*Uout(x,y)=U^m(x,y)*/
638    Cudd_Ref(Uout = U[m]);
639    if(pr>3){(void)fprintf(stdout,"Uout");Cudd_PrintDebug(bdd,Uout,n<<1,pr);}
640    /*Uin(x,y)=U^{m-1}(x,y)*/
641    Cudd_Ref(Uin = U[m-1]);
642    if(pr>3){(void)fprintf(stdout,"Uin");Cudd_PrintDebug(bdd,Uin,n<<1,pr);}
643
644    for(k = 0; k < MAXFPIT; k++) {
645        stats->fpit++;
646        /*LU(x,y)=Uin(x,y)\cdot\overline{\exists_{z}(Uin(z,y)\cdot\Pi(x,z))}*/
647        Cudd_Ref(p = Cudd_bddSwapVariables(bdd, Uin, x, z, n));
648        Cudd_Ref(r = Cudd_bddAndAbstract(bdd, p, prxz, zcube));
649        Cudd_RecursiveDeref(bdd,p);
650        Cudd_Ref(LU = Cudd_bddAnd(bdd, Uin, Cudd_Not(r)));
651        Cudd_RecursiveDeref(bdd,r);
652        if(pr>3){(void)fprintf(stdout,"LU");Cudd_PrintDebug(bdd,LU,n<<1,pr);}
653        /*D(x,y)= LU(x,y)\cdot \overline{\exists_{z}(LU(x,z)\cdot \Pi(y,z))}*/
654        Cudd_Ref(p = Cudd_bddSwapVariables(bdd, LU, y, z, n));
655        Cudd_Ref(r = Cudd_bddAndAbstract(bdd, p, pryz, zcube));
656        Cudd_RecursiveDeref(bdd,p);
657        Cudd_Ref(D = Cudd_bddAnd(bdd, LU, Cudd_Not(r)));
658        Cudd_RecursiveDeref(bdd,r);
659        Cudd_RecursiveDeref(bdd,LU);
660        if(pr>3){(void)fprintf(stdout,"D");Cudd_PrintDebug(bdd,D,n<<1,pr);}
661        /*Ml(y)=\exists_{x}D(x,y)*/
662        Cudd_Ref(Ml = Cudd_bddExistAbstract(bdd, D, xcube));
663        if(pr>3){(void)fprintf(stdout,"Ml");Cudd_PrintDebug(bdd,Ml,n,pr);}
664        /*P(x,y)=Ml(x)\cdot Uout(x,y)*/
665        Cudd_Ref(p = Cudd_bddSwapVariables(bdd, Ml, x, y, n));
666        Cudd_Ref(P = Cudd_bddAnd(bdd, p, Uout));
667        Cudd_RecursiveDeref(bdd,p);
668        Cudd_RecursiveDeref(bdd,Ml);
669        if(pr>3){(void)fprintf(stdout,"P");Cudd_PrintDebug(bdd,P,n<<1,pr);}
670        /*RU(x,y)= P(x,y)\cdot \overline{\exists_{z}(P(x,z)\cdot \Pi(y,z))}*/
671        Cudd_Ref(p = Cudd_bddSwapVariables(bdd, P, y, z, n));
672        Cudd_Ref(r = Cudd_bddAndAbstract(bdd, p, pryz, zcube));
673        Cudd_RecursiveDeref(bdd,p);
674        Cudd_Ref(RU = Cudd_bddAnd(bdd, P, Cudd_Not(r)));
675        Cudd_RecursiveDeref(bdd,r);
676        Cudd_RecursiveDeref(bdd,P);
677        if(pr>3){(void)fprintf(stdout,"RU");Cudd_PrintDebug(bdd,RU,n<<1,pr);}
678        /*Mr(x)=\exists_{y}RU(x,y)*/
679        Cudd_Ref(Mr = Cudd_bddExistAbstract(bdd, RU, ycube));
680        if(pr>3){(void)fprintf(stdout,"Mr");Cudd_PrintDebug(bdd,Mr,n,pr);}
681        /*D(x,y)=D(x,y)\cdot Mr(y)*/
682        Cudd_Ref(p = Cudd_bddSwapVariables(bdd, Mr, x, y, n));
683        Cudd_RecursiveDeref(bdd,Mr);
684        Cudd_Ref(q = Cudd_bddAnd(bdd, D, p));
685        Cudd_RecursiveDeref(bdd,p);
686        Cudd_RecursiveDeref(bdd,D);
687        D = q;
688        if(pr>3){(void)fprintf(stdout,"D");Cudd_PrintDebug(bdd,D,n<<1,pr);}
689        /*Uin(x,y)=Uin(x,y)-D(x,y)*/
690        Cudd_Ref(p = Cudd_bddAnd(bdd, Uin, Cudd_Not(D)));
691        Cudd_RecursiveDeref(bdd,Uin);
692        Uin = p;
693        if(pr>3){(void)fprintf(stdout,"Uin");Cudd_PrintDebug(bdd,Uin,n<<1,pr);}
694        /*Uout(x,y)=Uout(x,y)-RU(x,y)*/
695        Cudd_Ref(p = Cudd_bddAnd(bdd, Uout, Cudd_Not(RU)));
696        Cudd_RecursiveDeref(bdd,Uout);
697        Cudd_RecursiveDeref(bdd,RU);
698        Uout = p;
699        if(pr>3){(void)fprintf(stdout,"Uout");Cudd_PrintDebug(bdd,Uout,n<<1,pr);}
700        /*\KW{if}(($D(x,y)=zero$)~~\KW{or}~~($Uin(x,y)=zero$)~~\KW{or}
701          ($Uout(x,y)=zero$))~~KW{break}*/
702        if ((D == zero)||(Uin == zero)||(Uout == zero)) {
703            Cudd_RecursiveDeref(bdd,D);
704            break;
705        }
706        Cudd_RecursiveDeref(bdd,D);
707
708    } /* Next least fixed point iteration with smaller Uin and Uout */
709    if (k == MAXFPIT) (void) fprintf(stderr,
710        "Trellis: WARNING! MAXFPIT (%d) exceeded processing Layer %d.\n",
711        MAXFPIT, m);
712
713    if (Uin != zero) {
714        /* $U^{m-1}(x,y) = U^{m-1}(x,y)-Uin(x,y)$*/
715        Cudd_Ref(p = Cudd_bddAnd(bdd, U[m-1], Cudd_Not(Uin)));
716        Cudd_RecursiveDeref(bdd,U[m-1]);
717        U[m-1] = p;
718        /* $new^{m-1}(x,y) = \esists_yU^{m-1}(x,y)$*/
719        Cudd_RecursiveDeref(bdd,neW[m-1]);
720        Cudd_Ref(neW[m-1] = Cudd_bddExistAbstract(bdd, U[m-1], ycube));
721    }
722    if(pr>2){(void)fprintf(stdout,"U[%d]",m-1); Cudd_PrintDebug(bdd,U[m-1],n<<1,pr);}
723    if(pr>2){(void)fprintf(stdout,"new[%d]",m-1);
724                Cudd_PrintDebug(bdd,neW[m-1],n,pr);}
725
726    Cudd_RecursiveDeref(bdd,Uin);
727    Cudd_RecursiveDeref(bdd,Uout);
728
729    return;
730
731} /* end of trellis */
732
733
734/**Function********************************************************************
735
736  Synopsis    [Selects edges from a rhombus-type bilayer.]
737
738  Description [Selects edges from a rhombus-type bilayer. Used to pull flow.
739  Makes the left layer left-unique and the right layer right-unique. Prunes
740  and repeats until convergence to a maximal flow. It makes sure that all
741  intermediate points of the two-arc paths are disjoint at each iteration.]
742
743  SideEffects [None]
744
745  SeeAlso     [trellis hourglass rhombusPush]
746
747******************************************************************************/
748static void
749rhombus(
750  DdManager * bdd,
751  int  m /* center of current bilayer */,
752  DdNode ** neW,
753  DdNode ** U /* arrays for flow propagation */,
754  DdNode ** x,
755  DdNode ** y,
756  DdNode ** z /* arrays of variables for rows and columns */,
757  int  n /* number of x variables */,
758  DdNode * pryz,
759  DdNode * prxz /* priority functions */,
760  flowStats * stats)
761{
762    DdNode *one, *zero,
763             *p, *q, *r,
764             *Uin,              /* edges to be matched from U[m-1] */
765             *Uout,             /* edges to be matched from U[m] */
766             *P,
767             *LU, *RU,          /* left-unique and right-unique sets of edges */
768             *Ml, *Mr;          /* nodes matched from left and right */
769    int      k,
770             pr;                /* print control */
771
772    pr = stats->pr;
773
774    one = Cudd_ReadOne(bdd);
775    zero = Cudd_Not(one);
776
777    /*Uout(x,y)=U^m(x,y)*/
778    Cudd_Ref(Uout = U[m]);
779    if(pr>3){(void)fprintf(stdout,"Uout");Cudd_PrintDebug(bdd,Uout,n<<1,pr);}
780
781    /*Uin(x,y)=U^{m-1}(x,y)*/
782    Cudd_Ref(Uin = U[m-1]);
783    if(pr>3){(void)fprintf(stdout,"Uin");Cudd_PrintDebug(bdd,Uin,n<<1,pr);}
784
785    for(k = 0; k < MAXFPIT; k++) {
786        stats->fpit++;
787        /*LU(x,y)=Uin(x,y)\cdot\overline{\exists_{z}(Uin(z,y)\cdot\Pi(x,z))}*/
788        Cudd_Ref(p = Cudd_bddSwapVariables(bdd, Uin, x, z, n));
789        Cudd_Ref(r = Cudd_bddAndAbstract(bdd, p, prxz, zcube));
790        Cudd_RecursiveDeref(bdd,p);
791        Cudd_Ref(LU = Cudd_bddAnd(bdd, Uin, Cudd_Not(r)));
792        Cudd_RecursiveDeref(bdd,r);
793        if(pr>3){(void)fprintf(stdout,"LU");Cudd_PrintDebug(bdd,LU,n<<1,pr);}
794        /*Ml(y)=\exists_{x}LU(x,y)*/
795        Cudd_Ref(Ml = Cudd_bddExistAbstract(bdd, LU, xcube));
796        if(pr>3){(void)fprintf(stdout,"Ml");Cudd_PrintDebug(bdd,Ml,n,pr);}
797        /*P(x,y)=Ml(x)\cdot Uout(x,y)*/
798        Cudd_Ref(p = Cudd_bddSwapVariables(bdd, Ml, x, y, n));
799        Cudd_Ref(P = Cudd_bddAnd(bdd, p, Uout));
800        Cudd_RecursiveDeref(bdd,p);
801        Cudd_RecursiveDeref(bdd,Ml);
802        if(pr>3){(void)fprintf(stdout,"P");Cudd_PrintDebug(bdd,P,n<<1,pr);}
803        /*RU(x,y)= P(x,y)\cdot \overline{\exists_{z}(P(x,z)\cdot \Pi(y,z))}*/
804        Cudd_Ref(p = Cudd_bddSwapVariables(bdd, P, y, z, n));
805        Cudd_Ref(r = Cudd_bddAndAbstract(bdd, p, pryz, zcube));
806        Cudd_RecursiveDeref(bdd,p);
807        Cudd_Ref(RU = Cudd_bddAnd(bdd, P, Cudd_Not(r)));
808        Cudd_RecursiveDeref(bdd,r);
809        Cudd_RecursiveDeref(bdd,P);
810        if(pr>3){(void)fprintf(stdout,"RU");Cudd_PrintDebug(bdd,RU,n<<1,pr);}
811        /*Mr(x)=\exists_{y}RU(x,y)*/
812        Cudd_Ref(Mr = Cudd_bddExistAbstract(bdd, RU, ycube));
813        if(pr>3){(void)fprintf(stdout,"Mr");Cudd_PrintDebug(bdd,Mr,n,pr);}
814        /*LU(x,y)=LU(x,y)\cdot Mr(y)*/
815        Cudd_Ref(p = Cudd_bddSwapVariables(bdd, Mr, x, y, n));
816        Cudd_RecursiveDeref(bdd,Mr);
817        Cudd_Ref(q = Cudd_bddAnd(bdd, LU, p));
818        Cudd_RecursiveDeref(bdd,p);
819        Cudd_RecursiveDeref(bdd,LU);
820        LU = q;
821        if(pr>3){(void)fprintf(stdout,"LU");Cudd_PrintDebug(bdd,LU,n<<1,pr);}
822        /*Uin(x,y)=Uin(x,y)-LU(x,y)*/
823        Cudd_Ref(p = Cudd_bddAnd(bdd, Uin, Cudd_Not(LU)));
824        Cudd_RecursiveDeref(bdd,Uin);
825        Uin = p;
826        if(pr>3){(void)fprintf(stdout,"Uin");Cudd_PrintDebug(bdd,Uin,n<<1,pr);}
827        /*Uout(x,y)=Uout(x,y)-RU(x,y)*/
828        Cudd_Ref(p = Cudd_bddAnd(bdd, Uout, Cudd_Not(RU)));
829        Cudd_RecursiveDeref(bdd,Uout);
830        Cudd_RecursiveDeref(bdd,RU);
831        Uout = p;
832        if(pr>3){(void)fprintf(stdout,"Uout");Cudd_PrintDebug(bdd,Uout,n<<1,pr);}
833        /*\KW{if}(($LU(x,y)=zero$)~~\KW{or}~~($Uin(x,y)=zero$)~~\KW{or}
834          ($Uout(x,y)=zero$))~~KW{break}*/
835        if((LU == zero)||(Uin == zero)||(Uout == zero)) {
836            Cudd_RecursiveDeref(bdd,LU);
837            break;
838        }
839        Cudd_RecursiveDeref(bdd,LU);
840
841    } /* Next least fixed point iteration with smaller Uin and Uout */
842    if (k == MAXFPIT) (void) fprintf(stderr,
843        "Rhombus: WARNING! MAXFPIT (%d) exceeded processing Layer %d.\n",
844        MAXFPIT, m);
845
846    if (Uin != zero) {
847        /* $U^{m-1}(x,y) = U^{m-1}(x,y)-Uin(x,y)$*/
848        Cudd_Ref(p = Cudd_bddAnd(bdd, U[m-1], Cudd_Not(Uin)));
849        Cudd_RecursiveDeref(bdd,U[m-1]);
850        U[m-1] = p;
851        /* $new^{m-1}(x,y) = \esists_yU^{m-1}(x,y)$*/
852        Cudd_RecursiveDeref(bdd,neW[m-1]);
853        Cudd_Ref(neW[m-1] = Cudd_bddExistAbstract(bdd, U[m-1], ycube));
854    }
855    if(pr>2){(void)fprintf(stdout,"U[%d]",m-1); Cudd_PrintDebug(bdd,U[m-1],n<<1,pr);}
856    if(pr>2){
857        (void)fprintf(stdout,"new[%d]",m-1);
858        Cudd_PrintDebug(bdd,neW[m-1],n,pr);
859    }
860    Cudd_RecursiveDeref(bdd,Uin);
861    Cudd_RecursiveDeref(bdd,Uout);
862
863    return;
864
865} /* end of rhombus */
866
867
868/**Function********************************************************************
869
870  Synopsis    [Selects edges from a hourglass-type bilayer.]
871
872  Description [Selects edges from a hourglass-type bilayer. Used to
873  pull flow.  Method described in paper. More general, but more
874  expensive than the others.]
875
876  SideEffects [None]
877
878  SeeAlso     [trellis rhombus hourglassPush]
879
880******************************************************************************/
881static void
882hourglass(
883  DdManager * bdd,
884  int  m /* center level of current bilayer */,
885  DdNode ** neW,
886  DdNode ** U /* arrays for flow propagation */,
887  DdNode ** x,
888  DdNode ** y,
889  DdNode ** z /* arrays of variables for rows and columns */,
890  int  n /* number of x variables */,
891  DdNode * pryz,
892  DdNode * prxz /* priority functions */,
893  flowStats * stats)
894{
895    DdNode *one, *zero,
896             *przy,
897             *p, *q, *r,
898             *Uin,              /* edges to be matched from U[m-1] */
899             *Uout,             /* edges to be matched from U[m] */
900             *P, *Q,
901             *PA, *D;
902    int      k,
903             pr;                /* print control */
904
905    pr = stats->pr;
906
907    one = Cudd_ReadOne(bdd);
908    zero = Cudd_Not(one);
909
910    /* Build priority function. */
911    Cudd_Ref(p = Cudd_bddSwapVariables(bdd, pryz, y, z, n));
912    Cudd_Ref(przy = Cudd_bddAdjPermuteX(bdd,p,x,n));
913    Cudd_RecursiveDeref(bdd,p);
914    if(pr>2){(void)fprintf(stdout,"przy");Cudd_PrintDebug(bdd,przy,n*3,pr);}
915
916    /*Uout(x,y)=U^m(x,y)*/
917    Cudd_Ref(Uout = U[m]);
918    if(pr>1){(void)fprintf(stdout,"Uout");Cudd_PrintDebug(bdd,Uout,n<<1,pr);}
919
920    /*Uin(x,y)=U^{m-1}(x,y)*/
921    Cudd_Ref(Uin = U[m-1]);
922    if(pr>1){(void)fprintf(stdout,"Uin");Cudd_PrintDebug(bdd,Uin,n<<1,pr);}
923
924    for(k = 0; k < MAXFPIT; k++) {
925        stats->fpit++;
926        /*P(x,y,z)=Uin(x,y)\cdot Uout(y,z)*/
927        Cudd_Ref(p = Cudd_bddSwapVariables(bdd, Uout, y, z, n));
928        if(pr>2){(void)fprintf(stdout,"p");Cudd_PrintDebug(bdd,p,n<<1,pr);}
929        Cudd_Ref(q = Cudd_bddSwapVariables(bdd, p, x, y, n));
930        Cudd_RecursiveDeref(bdd,p);
931        if(pr>2){(void)fprintf(stdout,"q");Cudd_PrintDebug(bdd,q,n<<1,pr);}
932        Cudd_Ref(P = Cudd_bddAnd(bdd, Uin, q));
933        Cudd_RecursiveDeref(bdd,q);
934        if(pr>1){(void)fprintf(stdout,"P");Cudd_PrintDebug(bdd,P,n*3,pr);}
935        /*PA(x,z)=\exists_yP(x,y,z)*/
936        Cudd_Ref(PA = Cudd_bddExistAbstract(bdd, P, ycube));
937        if(pr>2){(void)fprintf(stdout,"PA");Cudd_PrintDebug(bdd,PA,n<<1,pr);}
938        if(pr>3) {
939            Cudd_Ref(p = Cudd_bddExistAbstract(bdd, PA, xcube));
940            (void) fprintf(stdout,"p");Cudd_PrintDebug(bdd,p,n,pr);
941            Cudd_RecursiveDeref(bdd,p);
942            Cudd_Ref(p = Cudd_bddExistAbstract(bdd, PA, zcube));
943            (void) fprintf(stdout,"p");Cudd_PrintDebug(bdd,p,n,pr);
944            Cudd_RecursiveDeref(bdd,p);
945        }
946        /*Q(x,z)= PA(x,z)\cdot \overline{\exists_{y}(PA(x,y)\cdot \Pi(z,y))}*/
947        Cudd_Ref(p = Cudd_bddSwapVariables(bdd, PA, y, z, n));
948        Cudd_Ref(r = Cudd_bddAndAbstract(bdd, p, przy, ycube));
949        Cudd_RecursiveDeref(bdd,p);
950        Cudd_Ref(Q = Cudd_bddAnd(bdd, PA, Cudd_Not(r)));
951        Cudd_RecursiveDeref(bdd,r);
952        Cudd_RecursiveDeref(bdd,PA);
953        if(pr>2){(void)fprintf(stdout,"Q");Cudd_PrintDebug(bdd,Q,n<<1,pr);}
954        if(pr>3) {
955            Cudd_Ref(p = Cudd_bddExistAbstract(bdd, Q, xcube));
956            (void) fprintf(stdout,"p");Cudd_PrintDebug(bdd,p,n,pr);
957            Cudd_RecursiveDeref(bdd,p);
958        }
959        /*D(x,z)= Q(x,z)\cdot \overline{\exists_{y}(Q(y,z)\cdot \Pi(x,y))}*/
960        Cudd_Ref(p = Cudd_bddSwapVariables(bdd, Q, x, y, n));
961        Cudd_Ref(q = Cudd_bddSwapVariables(bdd, prxz, y, z, n));
962        Cudd_Ref(r = Cudd_bddAndAbstract(bdd, p, q, ycube));
963        Cudd_RecursiveDeref(bdd,p);
964        Cudd_RecursiveDeref(bdd,q);
965        Cudd_Ref(D = Cudd_bddAnd(bdd, Q, Cudd_Not(r)));
966        Cudd_RecursiveDeref(bdd,r);
967        Cudd_RecursiveDeref(bdd,Q);
968        if(pr>1){(void)fprintf(stdout,"D");Cudd_PrintDebug(bdd,D,n<<1,pr);}
969        /*P(x,y,z)=P(x,y,z)\cdot D(x,z)*/
970        Cudd_Ref(p = Cudd_bddAnd(bdd, P, D));
971        Cudd_RecursiveDeref(bdd,D);
972        Cudd_RecursiveDeref(bdd,P);
973        P = p;
974        if(pr>2){(void)fprintf(stdout,"P");Cudd_PrintDebug(bdd,P,n*3,pr);}
975        /*Uin(x,y)=Uin(x,y)-\exists_zP(x,y,z)*/
976        Cudd_Ref(p = Cudd_bddExistAbstract(bdd, P, zcube));
977        if(pr>3){(void)fprintf(stdout,"p");Cudd_PrintDebug(bdd,p,n<<1,pr);}
978        Cudd_Ref(q = Cudd_bddAnd(bdd, Uin, Cudd_Not(p)));
979        Cudd_RecursiveDeref(bdd,p);
980        Cudd_RecursiveDeref(bdd,Uin);
981        Uin = q;
982        if(pr>1){(void)fprintf(stdout,"Uin");Cudd_PrintDebug(bdd,Uin,n<<1,pr);}
983        /*Uout(x,y)=Uout(x,y)-\exists_zP(z,x,y)*/
984        Cudd_Ref(p = Cudd_bddSwapVariables(bdd, P, x, y, n));
985        if(pr>3){(void)fprintf(stdout,"p");Cudd_PrintDebug(bdd,p,n*3,pr);}
986        Cudd_Ref(r = Cudd_bddSwapVariables(bdd, p, y, z, n));
987        Cudd_RecursiveDeref(bdd,p);
988        if(pr>3){(void)fprintf(stdout,"r");Cudd_PrintDebug(bdd,r,n*3,pr);}
989        Cudd_Ref(p = Cudd_bddExistAbstract(bdd, r, zcube));
990        Cudd_RecursiveDeref(bdd,r);
991        if(pr>3){(void)fprintf(stdout,"p");Cudd_PrintDebug(bdd,p,n<<1,pr);}
992        Cudd_Ref(q = Cudd_bddAnd(bdd, Uout, Cudd_Not(p)));
993        Cudd_RecursiveDeref(bdd,p);
994        Cudd_RecursiveDeref(bdd,Uout);
995        Uout = q;
996        if(pr>1){(void)fprintf(stdout,"Uout");Cudd_PrintDebug(bdd,Uout,n<<1,pr);}
997        /*\KW{if}(($P(x,y,z)=zero$)~~\KW{or}~~($Uin(x,y)=zero$)~~\KW{or}
998          ($Uout(x,y)=zero$))~~KW{break}*/
999        if((P == zero)||(Uin == zero)||(Uout == zero)) {
1000            Cudd_RecursiveDeref(bdd,P);
1001            break;
1002        }
1003        Cudd_RecursiveDeref(bdd,P);
1004
1005    } /* Next least fixed point iteration with smaller P */
1006    if (k == MAXFPIT) (void) fprintf(stderr,
1007        "Hourglass: WARNING! MAXFPIT (%d) exceeded processing Layer %d.\n",
1008        MAXFPIT, m);
1009
1010    if (Uin != zero) {
1011        /* $U^{m-1}(x,y) = U^{m-1}(x,y)-Uin(x,y)$*/
1012        Cudd_Ref(p = Cudd_bddAnd(bdd, U[m-1], Cudd_Not(Uin)));
1013        Cudd_RecursiveDeref(bdd,U[m-1]);
1014        U[m-1] = p;
1015        /* $new^{m-1}(x,y) = \esists_yU^{m-1}(x,y)$*/
1016        Cudd_RecursiveDeref(bdd,neW[m-1]);
1017        Cudd_Ref(neW[m-1] = Cudd_bddExistAbstract(bdd, U[m-1], ycube));
1018    }
1019    if(pr>1){(void)fprintf(stdout,"U[%d]",m-1); Cudd_PrintDebug(bdd,U[m-1],n<<1,pr);}
1020    if(pr>1){(void)fprintf(stdout,"new[%d]",m-1);
1021             Cudd_PrintDebug(bdd,neW[m-1],n,pr);}
1022
1023    Cudd_RecursiveDeref(bdd,Uin);
1024    Cudd_RecursiveDeref(bdd,Uout);
1025    Cudd_RecursiveDeref(bdd,przy);
1026
1027    return;
1028
1029} /* end of hourglass */
1030
1031
1032/**Function********************************************************************
1033
1034  Synopsis    [Pushes flow through useful edges.]
1035
1036  Description [Pushes flow from the source(s) to the sink(s) through
1037  useful edges.]
1038
1039  SideEffects [None]
1040
1041  SeeAlso     []
1042
1043******************************************************************************/
1044static void
1045maximal_push(
1046  DdManager * bdd,
1047  int  l /* Depth of layered network for current phase */,
1048  DdNode ** U /* array of edge sets for flow propagation */,
1049  DdNode ** F /* edge and flow relations */,
1050  DdNode ** x,
1051  DdNode ** y,
1052  DdNode ** z /* arrays of variables for rows and columns */,
1053  int  n /* number of x variables */,
1054  DdNode * pryz,
1055  DdNode * prxz /* priority functions */,
1056  flowStats * stats)
1057{
1058    DdNode *p, *q, *r,
1059           *UT,
1060           *lN, *cN, *rN; /* left, center, right nodes of bilayer */
1061    double mtl, mtc, mtr;
1062    int    m,
1063           pr;    /* print control */
1064
1065    pr = stats->pr;
1066
1067    if (l == 0) {
1068        /* F(x,y) = F(x,y) + U^{0}(x,y) */
1069        Cudd_Ref(q = Cudd_bddOr(bdd, *F, U[0]));
1070        Cudd_RecursiveDeref(bdd,*F);
1071        *F = q;
1072        if(pr>3){(void) fprintf(stdout,"F");Cudd_PrintDebug(bdd,*F,n<<1,pr);}
1073        return;
1074    }
1075
1076    for (m = 1; m < l; m++) {
1077        /* lN(x) = \exists_y U^{m-1}(x,y) */
1078        Cudd_Ref(lN = Cudd_bddExistAbstract(bdd, U[m-1], ycube));
1079        mtl = Cudd_CountMinterm(bdd, lN, n);
1080        Cudd_RecursiveDeref(bdd,lN);
1081
1082        /* cN(y) = \exists_x U^{m-1}(x,y) */
1083        Cudd_Ref(cN = Cudd_bddExistAbstract(bdd, U[m-1], xcube));
1084        mtc = Cudd_CountMinterm(bdd, cN, n);
1085        Cudd_RecursiveDeref(bdd,cN);
1086
1087        /* rN(y) = \exists_x U^{m}(x,y) */
1088        Cudd_Ref(rN = Cudd_bddExistAbstract(bdd, U[m], xcube));
1089        mtr = Cudd_CountMinterm(bdd, rN, n);
1090        Cudd_RecursiveDeref(bdd,rN);
1091
1092        if (pr > 0) {
1093            (void) fprintf(stdout, "layer = %d mtl = %g  mtc = %g  mtr = %g",
1094                           m, mtl, mtc, mtr);
1095        }
1096        if ((mtc > MANY_TIMES * mtl) && !(mtr > MANY_TIMES * mtl)) {
1097            if (pr>0) (void) fprintf(stdout, " R\n");
1098            rhombusPush(bdd, m, U, x, y, z, n, pryz, prxz, stats);
1099        } else if (mtl > MANY_TIMES * mtc) {
1100            if (pr>0) (void) fprintf(stdout, " H\n");
1101            hourglassPush(bdd, m, U, x, y, z, n, pryz, prxz, stats);
1102        } else {
1103            if (pr>0) (void) fprintf(stdout, " T\n");
1104            trellisPush(bdd, m, U, x, y, z, n, pryz, prxz, stats);
1105        }
1106
1107        /* F(x,y) = F(x,y) + U^{m-1}(x,y) \cdot \overline{F(y,x)} */
1108        Cudd_Ref(p = Cudd_bddSwapVariables(bdd, *F, x, y, n));
1109        Cudd_Ref(q = Cudd_bddAnd(bdd, Cudd_Not(p), U[m-1]));
1110        Cudd_RecursiveDeref(bdd,p);
1111        Cudd_Ref(r = Cudd_bddOr(bdd, *F, q));
1112        Cudd_RecursiveDeref(bdd,q);
1113        Cudd_RecursiveDeref(bdd,*F);
1114        *F = r;
1115        if(pr>3){(void) fprintf(stdout,"F");Cudd_PrintDebug(bdd,*F,n<<1,pr);}
1116
1117        /* F(x,y) = F(x,y) - U^{m-1}(y,x) */
1118        Cudd_Ref(r = Cudd_bddSwapVariables(bdd, U[m-1], x, y, n));
1119        Cudd_Ref(q = Cudd_bddAnd(bdd, *F, Cudd_Not(r)));
1120        Cudd_RecursiveDeref(bdd,r);
1121        Cudd_RecursiveDeref(bdd,*F);
1122        *F = q;
1123        if(pr>3){(void) fprintf(stdout,"F");Cudd_PrintDebug(bdd,*F,n<<1,pr);}
1124
1125    } /* Push maximal flow to next layer */
1126
1127    /*F(x,y)=F(x,y)+U^{l-1}(x,y)\cdot\overline{F(y,x)}*/
1128    Cudd_Ref(p = Cudd_bddSwapVariables(bdd, *F, x, y, n));
1129    Cudd_Ref(q = Cudd_bddAnd(bdd, Cudd_Not(p), U[l-1]));
1130    Cudd_RecursiveDeref(bdd,p);
1131    Cudd_Ref(r = Cudd_bddOr(bdd, *F, q));
1132    Cudd_RecursiveDeref(bdd,q);
1133    Cudd_RecursiveDeref(bdd,*F);
1134    *F = r;
1135    if(pr>3){(void) fprintf(stdout,"F");Cudd_PrintDebug(bdd,*F,n<<1,pr);}
1136
1137    /*F(y,x)=F(y,x)-U^{l-1}(x,y)*/
1138    Cudd_Ref(r = Cudd_bddSwapVariables(bdd, U[l-1], x, y, n));
1139    Cudd_Ref(q = Cudd_bddAnd(bdd, *F, Cudd_Not(r)));
1140    Cudd_RecursiveDeref(bdd,r);
1141    Cudd_RecursiveDeref(bdd,*F);
1142    *F = q;
1143    if(pr>1){(void) fprintf(stdout,"F");Cudd_PrintDebug(bdd,*F,n<<1,pr);}
1144
1145    /*UT(x,y)=\exists_y(U^{l-1}(y,x))\cdot U^{l}(x,y)*/
1146    Cudd_Ref(p = Cudd_bddExistAbstract(bdd, U[l-1], xcube));
1147    if(pr>4){(void) fprintf(stdout,"p");Cudd_PrintDebug(bdd,p,n,pr);}
1148    Cudd_Ref(q = Cudd_bddSwapVariables(bdd, p, x, y, n));
1149    Cudd_RecursiveDeref(bdd,p);
1150    if(pr>4){(void) fprintf(stdout,"q");Cudd_PrintDebug(bdd,q,n,pr);}
1151    Cudd_Ref(UT = Cudd_bddAnd(bdd, U[l], q));
1152    Cudd_RecursiveDeref(bdd,q);
1153    if(pr>2){(void) fprintf(stdout,"UT");Cudd_PrintDebug(bdd,UT,n<<1,pr);}
1154
1155    /*F(x,y)=F(x,y)+UT(x,y)\cdot\overline{F(y,x)}*/
1156    Cudd_Ref(p = Cudd_bddSwapVariables(bdd, *F, x, y, n));
1157    Cudd_Ref(q = Cudd_bddAnd(bdd, Cudd_Not(p), UT));
1158    Cudd_RecursiveDeref(bdd,p);
1159    Cudd_Ref(r = Cudd_bddOr(bdd, *F, q));
1160    Cudd_RecursiveDeref(bdd,q);
1161    Cudd_RecursiveDeref(bdd,*F);
1162    *F = r;
1163    if(pr>3){(void) fprintf(stdout,"F");Cudd_PrintDebug(bdd,*F,n<<1,pr);}
1164
1165    /*F(y,x)=F(y,x)-UT(x,y)*/
1166    Cudd_Ref(r = Cudd_bddSwapVariables(bdd, UT, x, y, n));
1167    Cudd_RecursiveDeref(bdd,UT);
1168    Cudd_Ref(q = Cudd_bddAnd(bdd, *F, Cudd_Not(r)));
1169    Cudd_RecursiveDeref(bdd,r);
1170    Cudd_RecursiveDeref(bdd,*F);
1171    *F = q;
1172    if(pr>1){(void) fprintf(stdout,"F");Cudd_PrintDebug(bdd,*F,n<<1,pr);}
1173
1174    return;
1175
1176} /* end of maximal_push */
1177
1178
1179/**Function********************************************************************
1180
1181  Synopsis    []
1182
1183  Description []
1184
1185  SideEffects [None]
1186
1187  SeeAlso     []
1188
1189******************************************************************************/
1190static void
1191trellisPush(
1192  DdManager * bdd,
1193  int  m /* Current layer */,
1194  DdNode ** U /* Array of edge sets for flow propagation */,
1195  DdNode ** x,
1196  DdNode ** y,
1197  DdNode ** z /* Arrays of variables for rows and columns */,
1198  int  n /* Number of x variables */,
1199  DdNode * pryz,
1200  DdNode * prxz /* Priority functions */,
1201  flowStats * stats)
1202{
1203    DdNode *one, *zero,
1204             *p, *r,
1205             *Uin,              /* Edges to be matched from U[m-1] */
1206             *Uout,             /* Edges to be matched from U[m] */
1207             *RU, *LU,
1208             *P,
1209             *Ml;
1210
1211    int      i,
1212             pr;          /* Print control */
1213
1214    pr = stats->pr;
1215
1216    one = Cudd_ReadOne(bdd);
1217    zero = Cudd_Not(one);
1218
1219    /*Uout(x,y)=U^m(x,y)*/
1220    Cudd_Ref(Uout = U[m]);
1221    if(pr>3){(void)fprintf(stdout,"Uout");Cudd_PrintDebug(bdd,Uout,n<<1,pr);}
1222
1223    /*Uin(x,y)=U^{m-1}(x,y)*/
1224    Cudd_Ref(Uin = U[m-1]);
1225    if(pr>3){(void)fprintf(stdout,"Uin");Cudd_PrintDebug(bdd,Uin,n<<1,pr);}
1226
1227    for(i=0; i<MAXFPIT; i++) {
1228        stats->fpit++;
1229        /*LU(x,y)=Uin(x,y)\cdot\overline{\exists_{z}(Uin(z,y)\cdot\Pi(x,z))}*/
1230        Cudd_Ref(p = Cudd_bddSwapVariables(bdd, Uin, x, z, n));
1231        Cudd_Ref(r = Cudd_bddAndAbstract(bdd, p, prxz, zcube));
1232        Cudd_RecursiveDeref(bdd,p);
1233        Cudd_Ref(LU = Cudd_bddAnd(bdd, Uin, Cudd_Not(r)));
1234        Cudd_RecursiveDeref(bdd,r);
1235        if(pr>3){(void)fprintf(stdout,"LU");Cudd_PrintDebug(bdd,LU,n<<1,pr);}
1236
1237        /*Ml(y)=\exists_{x}LU(x,y)*/
1238        Cudd_Ref(Ml = Cudd_bddExistAbstract(bdd, LU, xcube));
1239        if(pr>3){(void)fprintf(stdout,"Ml");Cudd_PrintDebug(bdd,Ml,n,pr);}
1240
1241        /*P(x,y)=Ml(x)\cdot Uout(x,y)*/
1242        Cudd_Ref(p = Cudd_bddSwapVariables(bdd, Ml, x, y, n));
1243        Cudd_Ref(P = Cudd_bddAnd(bdd, p, Uout));
1244        Cudd_RecursiveDeref(bdd,p);
1245        Cudd_RecursiveDeref(bdd,Ml);
1246        if(pr>3){(void)fprintf(stdout,"P");Cudd_PrintDebug(bdd,P,n<<1,pr);}
1247
1248        /*RU(x,y)= P(x,y)\cdot \overline{\exists_{z}(P(x,z)\cdot \Pi(y,z))}*/
1249        Cudd_Ref(p = Cudd_bddSwapVariables(bdd, P, y, z, n));
1250        Cudd_Ref(r = Cudd_bddAndAbstract(bdd, p, pryz, zcube));
1251        Cudd_RecursiveDeref(bdd,p);
1252        Cudd_Ref(RU = Cudd_bddAnd(bdd, P, Cudd_Not(r)));
1253        Cudd_RecursiveDeref(bdd,r);
1254        Cudd_RecursiveDeref(bdd,P);
1255        if(pr>3){(void)fprintf(stdout,"RU");Cudd_PrintDebug(bdd,RU,n<<1,pr);}
1256
1257        /*Uin(x,y)=Uin(x,y)-LU(x,y)*/
1258        Cudd_Ref(p = Cudd_bddAnd(bdd, Uin, Cudd_Not(LU)));
1259        Cudd_RecursiveDeref(bdd,Uin);
1260        Uin = p;
1261        if(pr>3){(void)fprintf(stdout,"Uin");Cudd_PrintDebug(bdd,Uin,n<<1,pr);}
1262
1263        /*Uout(x,y)=Uout(x,y)-RU(x,y)*/
1264        Cudd_Ref(p = Cudd_bddAnd(bdd, Uout, Cudd_Not(RU)));
1265        Cudd_RecursiveDeref(bdd,Uout);
1266        Cudd_RecursiveDeref(bdd,RU);
1267        Uout = p;
1268        if(pr>3){(void)fprintf(stdout,"Uout");Cudd_PrintDebug(bdd,Uout,n<<1,pr);}
1269
1270        /*\KW{if}(($LU(x,y)=zero$)~~\KW{or}~~($Uin(x,y)=zero$))~~KW{break}*/
1271        if((LU == zero)||(Uin == zero)) {
1272            Cudd_RecursiveDeref(bdd,LU);
1273            break;
1274        }
1275
1276        Cudd_RecursiveDeref(bdd,LU);
1277
1278    } /* Next least fixed point iteration with smaller Uin and Uout */
1279    if (i == MAXFPIT) (void) fprintf(stderr,
1280        "TrellisPush: ERROR! MAXFPIT (%d) exceeded at layer %d.\n",
1281        MAXFPIT, m);
1282
1283    /* $U^{m}(x,y) = U^{m}(x,y)-Uout(x,y)$*/
1284    if (Uout != zero) {
1285        Cudd_Ref(p = Cudd_bddAnd(bdd, U[m], Cudd_Not(Uout)));
1286        Cudd_RecursiveDeref(bdd,U[m]);
1287        U[m] = p;
1288        if(pr>3){(void)fprintf(stdout,"U[%d]",m);
1289            Cudd_PrintDebug(bdd,U[m],n<<1,pr);}
1290    }
1291
1292    Cudd_RecursiveDeref(bdd,Uin);
1293    Cudd_RecursiveDeref(bdd,Uout);
1294
1295} /* end of trellisPush */
1296
1297
1298/**Function********************************************************************
1299
1300  Synopsis    []
1301
1302  Description []
1303
1304  SideEffects [None]
1305
1306  SeeAlso     []
1307
1308******************************************************************************/
1309static void
1310rhombusPush(
1311  DdManager * bdd,
1312  int  m /* Current layer */,
1313  DdNode ** U /* Array of edge sets for flow propagation */,
1314  DdNode ** x,
1315  DdNode ** y,
1316  DdNode ** z /* Arrays of variables for rows and columns */,
1317  int  n /* Number of x variables */,
1318  DdNode * pryz,
1319  DdNode * prxz /* Priority functions */,
1320  flowStats * stats)
1321{
1322    DdNode *one, *zero,
1323             *p, *r,
1324             *Uin,              /* Edges to be matched from U[m-1] */
1325             *Uout,             /* Edges to be matched from U[m] */
1326             *RU, *LU,
1327             *P,
1328             *Ml;
1329
1330    int      i,
1331             pr;          /* Print control */
1332
1333    pr = stats->pr;
1334
1335    one = Cudd_ReadOne(bdd);
1336    zero = Cudd_Not(one);
1337
1338    /*Uout(x,y)=U^m(x,y)*/
1339    Cudd_Ref(Uout = U[m]);
1340    if(pr>3){(void)fprintf(stdout,"Uout");Cudd_PrintDebug(bdd,Uout,n<<1,pr);}
1341
1342    /*Uin(x,y)=U^{m-1}(x,y)*/
1343    Cudd_Ref(Uin = U[m-1]);
1344    if(pr>3){(void)fprintf(stdout,"Uin");Cudd_PrintDebug(bdd,Uin,n<<1,pr);}
1345
1346    for(i = 0; i < MAXFPIT; i++) {
1347        stats->fpit++;
1348        /*RU(x,y)=Uin(x,y)\cdot\overline{\exists_{z}(Uin(x,z)\cdot\Pi(y,z))}*/
1349        Cudd_Ref(p = Cudd_bddSwapVariables(bdd, Uin, y, z, n));
1350        Cudd_Ref(r = Cudd_bddAndAbstract(bdd, p, pryz, zcube));
1351        Cudd_RecursiveDeref(bdd,p);
1352        Cudd_Ref(RU = Cudd_bddAnd(bdd, Uin, Cudd_Not(r)));
1353        Cudd_RecursiveDeref(bdd,r);
1354        if(pr>3){(void)fprintf(stdout,"RU");Cudd_PrintDebug(bdd,RU,n<<1,pr);}
1355
1356        /*Ml(y)=\exists_{x}RU(x,y)*/
1357        Cudd_Ref(Ml = Cudd_bddExistAbstract(bdd, RU, xcube));
1358        if(pr>3){(void)fprintf(stdout,"Ml");Cudd_PrintDebug(bdd,Ml,n,pr);}
1359
1360        /*P(x,y)=Ml(x)\cdot Uout(x,y)*/
1361        Cudd_Ref(p = Cudd_bddSwapVariables(bdd, Ml, x, y, n));
1362        Cudd_Ref(P = Cudd_bddAnd(bdd, p, Uout));
1363        Cudd_RecursiveDeref(bdd,p);
1364        Cudd_RecursiveDeref(bdd,Ml);
1365        if(pr>3){(void)fprintf(stdout,"P");Cudd_PrintDebug(bdd,P,n<<1,pr);}
1366
1367        /*LU(x,y)= P(x,y)\cdot \overline{\exists_{z}(P(z,y)\cdot \Pi(x,z))}*/
1368        Cudd_Ref(p = Cudd_bddSwapVariables(bdd, P, x, z, n));
1369        Cudd_Ref(r = Cudd_bddAndAbstract(bdd, p, prxz, zcube));
1370        Cudd_RecursiveDeref(bdd,p);
1371        Cudd_Ref(LU = Cudd_bddAnd(bdd, P, Cudd_Not(r)));
1372        Cudd_RecursiveDeref(bdd,r);
1373        Cudd_RecursiveDeref(bdd,P);
1374        if(pr>3){(void)fprintf(stdout,"LU");Cudd_PrintDebug(bdd,LU,n<<1,pr);}
1375
1376        /*Uin(x,y)=Uin(x,y)-RU(x,y)*/
1377        Cudd_Ref(p = Cudd_bddAnd(bdd, Uin, Cudd_Not(RU)));
1378        Cudd_RecursiveDeref(bdd,Uin);
1379        Uin = p;
1380        if(pr>3){(void)fprintf(stdout,"Uin");Cudd_PrintDebug(bdd,Uin,n<<1,pr);}
1381
1382        /*Uout(x,y)=Uout(x,y)-LU(x,y)*/
1383        Cudd_Ref(p = Cudd_bddAnd(bdd, Uout, Cudd_Not(LU)));
1384        Cudd_RecursiveDeref(bdd,Uout);
1385        Cudd_RecursiveDeref(bdd,LU);
1386        Uout = p;
1387        if(pr>3){(void)fprintf(stdout,"Uout");Cudd_PrintDebug(bdd,Uout,n<<1,pr);}
1388
1389        /*\KW{if}(($RU(x,y)=zero$)~~\KW{or}~~($Uin(x,y)=zero$))~~KW{break}*/
1390        if((RU == zero)||(Uin == zero)) {
1391            Cudd_RecursiveDeref(bdd,RU);
1392            break;
1393        }
1394
1395        Cudd_RecursiveDeref(bdd,RU);
1396
1397    } /* Next least fixed point iteration with smaller Uin and Uout */
1398    if (i == MAXFPIT) (void) fprintf(stderr,
1399        "RhombusPush: ERROR! MAXFPIT (%d) exceeded at layer %d.\n",
1400        MAXFPIT, m);
1401
1402    /* $U^{m}(x,y) = U^{m}(x,y)-Uout(x,y)$*/
1403    if (Uout != zero) {
1404        Cudd_Ref(p = Cudd_bddAnd(bdd, U[m], Cudd_Not(Uout)));
1405        Cudd_RecursiveDeref(bdd,U[m]);
1406        U[m] = p;
1407        if(pr>3){(void)fprintf(stdout,"U[%d]",m);
1408            Cudd_PrintDebug(bdd,U[m],n<<1,pr);}
1409    }
1410
1411    Cudd_RecursiveDeref(bdd,Uin);
1412    Cudd_RecursiveDeref(bdd,Uout);
1413
1414    return;
1415
1416} /* end of rhombusPush */
1417
1418
1419/**Function********************************************************************
1420
1421  Synopsis    []
1422
1423  Description []
1424
1425  SideEffects [None]
1426
1427  SeeAlso     []
1428
1429******************************************************************************/
1430static void
1431hourglassPush(
1432  DdManager * bdd,
1433  int  m /* Current layer */,
1434  DdNode ** U /* Array of edge sets for flow propagation */,
1435  DdNode ** x,
1436  DdNode ** y,
1437  DdNode ** z /* Arrays of variables for rows and columns */,
1438  int  n /* Number of x variables */,
1439  DdNode * pryz,
1440  DdNode * prxz /* Priority functions */,
1441  flowStats * stats)
1442{
1443    DdNode *one, *zero,
1444             *przy,
1445             *p, *q, *r,
1446             *Uin,              /* Edges to be matched from U[m-1] */
1447             *Uout,             /* Edges to be matched from U[m] */
1448             *P, *Q,
1449             *PA, *D;
1450
1451    int      i,
1452             pr;          /* Print control */
1453
1454    pr = stats->pr;
1455
1456    one = Cudd_ReadOne(bdd);
1457    zero = Cudd_Not(one);
1458
1459    /* Build priority function. */
1460    Cudd_Ref(p = Cudd_bddSwapVariables(bdd, pryz, y, z, n));
1461    Cudd_Ref(przy = Cudd_bddAdjPermuteX(bdd,p,x,n));
1462    Cudd_RecursiveDeref(bdd,p);
1463    if(pr>2){(void)fprintf(stdout,"przy");Cudd_PrintDebug(bdd,przy,n*3,pr);}
1464
1465    /*Uout(x,y)=U^m(x,y)*/
1466    Cudd_Ref(Uout = U[m]);
1467    if(pr>3){(void)fprintf(stdout,"Uout");Cudd_PrintDebug(bdd,Uout,n<<1,pr);}
1468
1469    /*Uin(x,y)=U^{m-1}(x,y)*/
1470    Cudd_Ref(Uin = U[m-1]);
1471    if(pr>3){(void)fprintf(stdout,"Uin");Cudd_PrintDebug(bdd,Uin,n<<1,pr);}
1472
1473    for(i = 0; i < MAXFPIT; i++) {
1474        stats->fpit++;
1475        /*P(x,y,z)=Uin(x,y)\cdot Uout(y,z)*/
1476        Cudd_Ref(p = Cudd_bddSwapVariables(bdd, Uout, y, z, n));
1477        if(pr>2){(void)fprintf(stdout,"p");Cudd_PrintDebug(bdd,p,n<<1,pr);}
1478        Cudd_Ref(q = Cudd_bddSwapVariables(bdd, p, x, y, n));
1479        Cudd_RecursiveDeref(bdd,p);
1480        if(pr>2){(void)fprintf(stdout,"q");Cudd_PrintDebug(bdd,q,n<<1,pr);}
1481        Cudd_Ref(P = Cudd_bddAnd(bdd, Uin, q));
1482        Cudd_RecursiveDeref(bdd,q);
1483        if(pr>1){(void)fprintf(stdout,"P");Cudd_PrintDebug(bdd,P,n*3,pr);}
1484
1485        /*PA(x,z)=\exists_yP(x,y,z)*/
1486        Cudd_Ref(PA = Cudd_bddExistAbstract(bdd, P, ycube));
1487        if(pr>2){(void)fprintf(stdout,"PA");Cudd_PrintDebug(bdd,PA,n<<1,pr);}
1488        if(pr>3) {
1489            Cudd_Ref(p = Cudd_bddExistAbstract(bdd, PA, xcube));
1490            (void) fprintf(stdout,"p");Cudd_PrintDebug(bdd,p,n,pr);
1491            Cudd_RecursiveDeref(bdd,p);
1492            Cudd_Ref(p = Cudd_bddExistAbstract(bdd, PA, zcube));
1493            (void) fprintf(stdout,"p");Cudd_PrintDebug(bdd,p,n,pr);
1494            Cudd_RecursiveDeref(bdd,p);
1495        }
1496
1497        /*Q(x,z)= PA(x,z)\cdot \overline{\exists_{y}(PA(x,y)\cdot \Pi(z,y))}*/
1498        Cudd_Ref(p = Cudd_bddSwapVariables(bdd, PA, y, z, n));
1499        Cudd_Ref(r = Cudd_bddAndAbstract(bdd, p, przy, ycube));
1500        Cudd_RecursiveDeref(bdd,p);
1501        Cudd_Ref(Q = Cudd_bddAnd(bdd, PA, Cudd_Not(r)));
1502        Cudd_RecursiveDeref(bdd,r);
1503        Cudd_RecursiveDeref(bdd,PA);
1504        if(pr>2){(void)fprintf(stdout,"Q");Cudd_PrintDebug(bdd,Q,n<<1,pr);}
1505        if(pr>3) {
1506            Cudd_Ref(p = Cudd_bddExistAbstract(bdd, Q, xcube));
1507            (void) fprintf(stdout,"p");Cudd_PrintDebug(bdd,p,n,pr);
1508            Cudd_RecursiveDeref(bdd,p);
1509        }
1510
1511        /*D(x,z)= Q(x,z)\cdot \overline{\exists_{y}(Q(y,z)\cdot \Pi(x,y))}*/
1512        Cudd_Ref(p = Cudd_bddSwapVariables(bdd, Q, x, y, n));
1513        Cudd_Ref(q = Cudd_bddSwapVariables(bdd, prxz, y, z, n));
1514        Cudd_Ref(r = Cudd_bddAndAbstract(bdd, p, q, ycube));
1515        Cudd_RecursiveDeref(bdd,p);
1516        Cudd_RecursiveDeref(bdd,q);
1517        Cudd_Ref(D = Cudd_bddAnd(bdd, Q, Cudd_Not(r)));
1518        Cudd_RecursiveDeref(bdd,r);
1519        Cudd_RecursiveDeref(bdd,Q);
1520        if(pr>1){(void)fprintf(stdout,"D");Cudd_PrintDebug(bdd,D,n<<1,pr);}
1521
1522        /*P(x,y,z)=P(x,y,z)\cdot D(x,z)*/
1523        Cudd_Ref(p = Cudd_bddAnd(bdd, P, D));
1524        Cudd_RecursiveDeref(bdd,D);
1525        Cudd_RecursiveDeref(bdd,P);
1526        P = p;
1527        if(pr>2){(void)fprintf(stdout,"P");Cudd_PrintDebug(bdd,P,n*3,pr);}
1528
1529        /*Uin(x,y)=Uin(x,y)-\exists_zP(x,y,z)*/
1530        Cudd_Ref(p = Cudd_bddExistAbstract(bdd, P, zcube));
1531        if(pr>3){(void)fprintf(stdout,"p");Cudd_PrintDebug(bdd,p,n<<1,pr);}
1532        Cudd_Ref(q = Cudd_bddAnd(bdd, Uin, Cudd_Not(p)));
1533        Cudd_RecursiveDeref(bdd,p);
1534        Cudd_RecursiveDeref(bdd,Uin);
1535        Uin = q;
1536        if(pr>1){(void)fprintf(stdout,"Uin");Cudd_PrintDebug(bdd,Uin,n<<1,pr);}
1537
1538        /*Uout(x,y)=Uout(x,y)-\exists_zP(z,x,y)*/
1539        Cudd_Ref(p = Cudd_bddSwapVariables(bdd, P, x, y, n));
1540        if(pr>3){(void)fprintf(stdout,"p");Cudd_PrintDebug(bdd,p,n*3,pr);}
1541        Cudd_Ref(r = Cudd_bddSwapVariables(bdd, p, y, z, n));
1542        Cudd_RecursiveDeref(bdd,p);
1543        if(pr>3){(void)fprintf(stdout,"r");Cudd_PrintDebug(bdd,r,n*3,pr);}
1544        Cudd_Ref(p = Cudd_bddExistAbstract(bdd, r, zcube));
1545        Cudd_RecursiveDeref(bdd,r);
1546        if(pr>3){(void)fprintf(stdout,"p");Cudd_PrintDebug(bdd,p,n<<1,pr);}
1547        Cudd_Ref(q = Cudd_bddAnd(bdd, Uout, Cudd_Not(p)));
1548        Cudd_RecursiveDeref(bdd,p);
1549        Cudd_RecursiveDeref(bdd,Uout);
1550        Uout = q;
1551        if(pr>1){(void)fprintf(stdout,"Uout");Cudd_PrintDebug(bdd,Uout,n<<1,pr);}
1552
1553        /*\KW{if}(($P(x,y,z)=zero$)~~\KW{or}~~($Uin(x,y)=zero$)~~\KW{or}
1554          ($Uout(x,y)=zero$))~~KW{break}*/
1555        if((P == zero)||(Uin == zero)||(Uout == zero)) {
1556            Cudd_RecursiveDeref(bdd,P);
1557            break;
1558        }
1559
1560        Cudd_RecursiveDeref(bdd,P);
1561
1562    } /* Next least fixed point iteration with smaller P */
1563    if (i == MAXFPIT) (void) fprintf(stderr,
1564        "HourglassPush: ERROR! MAXFPIT (%d) exceeded at layer %d.\n",
1565        MAXFPIT, m);
1566
1567    /* $U^{m}(x,y) = U^{m}(x,y)-Uout(x,y)$*/
1568    if (Uout != zero) {
1569        Cudd_Ref(p = Cudd_bddAnd(bdd, U[m], Cudd_Not(Uout)));
1570        Cudd_RecursiveDeref(bdd,U[m]);
1571        U[m] = p;
1572    }
1573    if(pr>1){(void)fprintf(stdout,"U[%d]",m); Cudd_PrintDebug(bdd,U[m],n<<1,pr);}
1574
1575    Cudd_RecursiveDeref(bdd,Uin);
1576    Cudd_RecursiveDeref(bdd,Uout);
1577    Cudd_RecursiveDeref(bdd,przy);
1578
1579    return;
1580
1581} /* end of hourglassPush */
Note: See TracBrowser for help on using the repository browser.