distr/cudd/cuddMatMult.c

   1 /**CFile***********************************************************************
   2
   3   FileName    [cuddMatMult.c]
   4
   5   PackageName [cudd]
   6
   7   Synopsis    [Matrix multiplication functions.]
   8
   9   Description [External procedures included in this module:
  10                 <ul>
  11                 <li> Cudd_addMatrixMultiply()
  12                 <li> Cudd_addTimesPlus()
  13                 <li> Cudd_addTriangle()
  14                 <li> Cudd_addOuterSum()
  15                 </ul>
  16         Static procedures included in this module:
  17                 <ul>
  18                 <li> addMMRecur()
  19                 <li> addTriangleRecur()
  20                 <li> cuddAddOuterSumRecur()
  21                 </ul>]
  22
  23   Author      [Fabio Somenzi]
  24
  25   Copyright   [Copyright (c) 1995-2004, Regents of the University of Colorado
  26
  27   All rights reserved.
  28
  29   Redistribution and use in source and binary forms, with or without
  30   modification, are permitted provided that the following conditions
  31   are met:
  32
  33   Redistributions of source code must retain the above copyright
  34   notice, this list of conditions and the following disclaimer.
  35
  36   Redistributions in binary form must reproduce the above copyright
  37   notice, this list of conditions and the following disclaimer in the
  38   documentation and/or other materials provided with the distribution.
  39
  40   Neither the name of the University of Colorado nor the names of its
  41   contributors may be used to endorse or promote products derived from
  42   this software without specific prior written permission.
  43
  44   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
  45   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
  46   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
  47   FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
  48   COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
  49   INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
  50   BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
  51   LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
  52   CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  53   LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
  54   ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  55   POSSIBILITY OF SUCH DAMAGE.]
  56
  57 ******************************************************************************/
  58
  59 #include "util.h"
  60 #include "cuddInt.h"
  61
  62
  63 /*---------------------------------------------------------------------------*/
  64 /* Constant declarations                                                     */
  65 /*---------------------------------------------------------------------------*/
  66
  67
  68 /*---------------------------------------------------------------------------*/
  69 /* Stucture declarations                                                     */
  70 /*---------------------------------------------------------------------------*/
  71
  72
  73 /*---------------------------------------------------------------------------*/
  74 /* Type declarations                                                         */
  75 /*---------------------------------------------------------------------------*/
  76
  77
  78 /*---------------------------------------------------------------------------*/
  79 /* Variable declarations                                                     */
  80 /*---------------------------------------------------------------------------*/
  81
  82 #ifndef lint
  83 static char rcsid[] DD_UNUSED = "$Id: cuddMatMult.c,v 1.17 2004/08/13 18:04:50 fabio Exp $";
  84 #endif
  85
  86 /*---------------------------------------------------------------------------*/
  87 /* Macro declarations                                                        */
  88 /*---------------------------------------------------------------------------*/
  89
  90
  91 /**AutomaticStart*************************************************************/
  92
  93 /*---------------------------------------------------------------------------*/
  94 /* Static function prototypes                                                */
  95 /*---------------------------------------------------------------------------*/
  96
  97 static DdNode * addMMRecur (DdManager *dd, DdNode *A, DdNode *B, int topP, int *vars);
  98 static DdNode * addTriangleRecur (DdManager *dd, DdNode *f, DdNode *g, int *vars, DdNode *cube);
  99 static DdNode * cuddAddOuterSumRecur (DdManager *dd, DdNode *M, DdNode *r, DdNode *c);
 100
 101 /**AutomaticEnd***************************************************************/
 102
 103
 104 /*---------------------------------------------------------------------------*/
 105 /* Definition of exported functions                                          */
 106 /*---------------------------------------------------------------------------*/
 107
 108 /**Function********************************************************************
 109
 110   Synopsis [Calculates the product of two matrices represented as
 111   ADDs.]
 112
 113   Description [Calculates the product of two matrices, A and B,
 114   represented as ADDs. This procedure implements the quasiring multiplication
 115   algorithm.  A is assumed to depend on variables x (rows) and z
 116   (columns).  B is assumed to depend on variables z (rows) and y
 117   (columns).  The product of A and B then depends on x (rows) and y
 118   (columns).  Only the z variables have to be explicitly identified;
 119   they are the "summation" variables.  Returns a pointer to the
 120   result if successful; NULL otherwise.]
 121
 122   SideEffects [None]
 123
 124   SeeAlso     [Cudd_addTimesPlus Cudd_addTriangle Cudd_bddAndAbstract]
 125
 126 ******************************************************************************/
 127 DdNode *
 128 Cudd_addMatrixMultiply(
 129   DdManager * dd,
 130   DdNode * A,
 131   DdNode * B,
 132   DdNode ** z,
 133   int  nz)
 134 {
 135     int i, nvars, *vars;
 136     DdNode *res;
 137
 138     /* Array vars says what variables are "summation" variables. */
 139     nvars = dd->size;
 140     vars = ALLOC(int,nvars);
 141     if (vars == NULL) {
 142         dd->errorCode = CUDD_MEMORY_OUT;
 143         return(NULL);
 144     }
 145     for (i = 0; i < nvars; i++) {
 146         vars[i] = 0;
 147     }
 148     for (i = 0; i < nz; i++) {
 149         vars[z[i]->index] = 1;
 150     }
 151
 152     do {
 153         dd->reordered = 0;
 154         res = addMMRecur(dd,A,B,-1,vars);
 155     } while (dd->reordered == 1);
 156     FREE(vars);
 157     return(res);
 158
 159 } /* end of Cudd_addMatrixMultiply */
 160
 161
 162 /**Function********************************************************************
 163
 164   Synopsis    [Calculates the product of two matrices represented as
 165   ADDs.]
 166
 167   Description [Calculates the product of two matrices, A and B,
 168   represented as ADDs, using the CMU matrix by matrix multiplication
 169   procedure by Clarke et al..  Matrix A has x's as row variables and z's
 170   as column variables, while matrix B has z's as row variables and y's
 171   as column variables. Returns the pointer to the result if successful;
 172   NULL otherwise. The resulting matrix has x's as row variables and y's
 173   as column variables.]
 174
 175   SideEffects [None]
 176
 177   SeeAlso     [Cudd_addMatrixMultiply]
 178
 179 ******************************************************************************/
 180 DdNode *
 181 Cudd_addTimesPlus(
 182   DdManager * dd,
 183   DdNode * A,
 184   DdNode * B,
 185   DdNode ** z,
 186   int  nz)
 187 {
 188     DdNode *w, *cube, *tmp, *res;
 189     int i;
 190     tmp = Cudd_addApply(dd,Cudd_addTimes,A,B);
 191     if (tmp == NULL) return(NULL);
 192     Cudd_Ref(tmp);
 193     Cudd_Ref(cube = DD_ONE(dd));
 194     for (i = nz-1; i >= 0; i--) {
 195          w = Cudd_addIte(dd,z[i],cube,DD_ZERO(dd));
 196          if (w == NULL) {
 197             Cudd_RecursiveDeref(dd,tmp);
 198             return(NULL);
 199          }
 200          Cudd_Ref(w);
 201          Cudd_RecursiveDeref(dd,cube);
 202          cube = w;
 203     }
 204     res = Cudd_addExistAbstract(dd,tmp,cube);
 205     if (res == NULL) {
 206         Cudd_RecursiveDeref(dd,tmp);
 207         Cudd_RecursiveDeref(dd,cube);
 208         return(NULL);
 209     }
 210     Cudd_Ref(res);
 211     Cudd_RecursiveDeref(dd,cube);
 212     Cudd_RecursiveDeref(dd,tmp);
 213     Cudd_Deref(res);
 214     return(res);
 215
 216 } /* end of Cudd_addTimesPlus */
 217
 218
 219 /**Function********************************************************************
 220
 221   Synopsis    [Performs the triangulation step for the shortest path
 222   computation.]
 223
 224   Description [Implements the semiring multiplication algorithm used in
 225   the triangulation step for the shortest path computation.  f
 226   is assumed to depend on variables x (rows) and z (columns).  g is
 227   assumed to depend on variables z (rows) and y (columns).  The product
 228   of f and g then depends on x (rows) and y (columns).  Only the z
 229   variables have to be explicitly identified; they are the
 230   "abstraction" variables.  Returns a pointer to the result if
 231   successful; NULL otherwise. ]
 232
 233   SideEffects [None]
 234
 235   SeeAlso     [Cudd_addMatrixMultiply Cudd_bddAndAbstract]
 236
 237 ******************************************************************************/
 238 DdNode *
 239 Cudd_addTriangle(
 240   DdManager * dd,
 241   DdNode * f,
 242   DdNode * g,
 243   DdNode ** z,
 244   int  nz)
 245 {
 246     int    i, nvars, *vars;
 247     DdNode *res, *cube;
 248
 249     nvars = dd->size;
 250     vars = ALLOC(int, nvars);
 251     if (vars == NULL) {
 252         dd->errorCode = CUDD_MEMORY_OUT;
 253         return(NULL);
 254     }
 255     for (i = 0; i < nvars; i++) vars[i] = -1;
 256     for (i = 0; i < nz; i++) vars[z[i]->index] = i;
 257     cube = Cudd_addComputeCube(dd, z, NULL, nz);
 258     if (cube == NULL) {
 259         FREE(vars);
 260         return(NULL);
 261     }
 262     cuddRef(cube);
 263
 264     do {
 265         dd->reordered = 0;
 266         res = addTriangleRecur(dd, f, g, vars, cube);
 267     } while (dd->reordered == 1);
 268     if (res != NULL) cuddRef(res);
 269     Cudd_RecursiveDeref(dd,cube);
 270     if (res != NULL) cuddDeref(res);
 271     FREE(vars);
 272     return(res);
 273
 274 } /* end of Cudd_addTriangle */
 275
 276
 277 /**Function********************************************************************
 278
 279   Synopsis    [Takes the minimum of a matrix and the outer sum of two vectors.]
 280
 281   Description [Takes the pointwise minimum of a matrix and the outer
 282   sum of two vectors.  This procedure is used in the Floyd-Warshall
 283   all-pair shortest path algorithm.  Returns a pointer to the result if
 284   successful; NULL otherwise.]
 285
 286   SideEffects [None]
 287
 288   SeeAlso     []
 289
 290 ******************************************************************************/
 291 DdNode *
 292 Cudd_addOuterSum(
 293   DdManager *dd,
 294   DdNode *M,
 295   DdNode *r,
 296   DdNode *c)
 297 {
 298     DdNode *res;
 299
 300     do {
 301         dd->reordered = 0;
 302         res = cuddAddOuterSumRecur(dd, M, r, c);
 303     } while (dd->reordered == 1);
 304     return(res);
 305
 306 } /* end of Cudd_addOuterSum */
 307
 308
 309 /*---------------------------------------------------------------------------*/
 310 /* Definition of internal functions                                          */
 311 /*---------------------------------------------------------------------------*/
 312
 313
 314 /*---------------------------------------------------------------------------*/
 315 /* Definition of static functions                                            */
 316 /*---------------------------------------------------------------------------*/
 317
 318 /**Function********************************************************************
 319
 320   Synopsis    [Performs the recursive step of Cudd_addMatrixMultiply.]
 321
 322   Description [Performs the recursive step of Cudd_addMatrixMultiply.
 323   Returns a pointer to the result if successful; NULL otherwise.]
 324
 325   SideEffects [None]
 326
 327 ******************************************************************************/
 328 static DdNode *
 329 addMMRecur(
 330   DdManager * dd,
 331   DdNode * A,
 332   DdNode * B,
 333   int  topP,
 334   int * vars)
 335 {
 336     DdNode *zero,
 337            *At,         /* positive cofactor of first operand */
 338            *Ae,         /* negative cofactor of first operand */
 339            *Bt,         /* positive cofactor of second operand */
 340            *Be,         /* negative cofactor of second operand */
 341            *t,          /* positive cofactor of result */
 342            *e,          /* negative cofactor of result */
 343            *scaled,     /* scaled result */
 344            *add_scale,  /* ADD representing the scaling factor */
 345            *res;
 346     int i;              /* loop index */
 347     double scale;       /* scaling factor */
 348     int index;          /* index of the top variable */
 349     CUDD_VALUE_TYPE value;
 350     unsigned int topA, topB, topV;
 351     DD_CTFP cacheOp;
 352
 353     statLine(dd);
 354     zero = DD_ZERO(dd);
 355
 356     if (A == zero || B == zero) {
 357         return(zero);
 358     }
 359
 360     if (cuddIsConstant(A) && cuddIsConstant(B)) {
 361         /* Compute the scaling factor. It is 2^k, where k is the
 362         ** number of summation variables below the current variable.
 363         ** Indeed, these constants represent blocks of 2^k identical
 364         ** constant values in both A and B.
 365         */
 366         value = cuddV(A) * cuddV(B);
 367         for (i = 0; i < dd->size; i++) {
 368             if (vars[i]) {
 369                 if (dd->perm[i] > topP) {
 370                     value *= (CUDD_VALUE_TYPE) 2;
 371                 }
 372             }
 373         }
 374         res = cuddUniqueConst(dd, value);
 375         return(res);
 376     }
 377
 378     /* Standardize to increase cache efficiency. Clearly, A*B != B*A
 379     ** in matrix multiplication. However, which matrix is which is
 380     ** determined by the variables appearing in the ADDs and not by
 381     ** which one is passed as first argument.
 382     */
 383     if (A > B) {
 384         DdNode *tmp = A;
 385         A = B;
 386         B = tmp;
 387     }
 388
 389     topA = cuddI(dd,A->index); topB = cuddI(dd,B->index);
 390     topV = ddMin(topA,topB);
 391
 392     cacheOp = (DD_CTFP) addMMRecur;
 393     res = cuddCacheLookup2(dd,cacheOp,A,B);
 394     if (res != NULL) {
 395         /* If the result is 0, there is no need to normalize.
 396         ** Otherwise we count the number of z variables between
 397         ** the current depth and the top of the ADDs. These are
 398         ** the missing variables that determine the size of the
 399         ** constant blocks.
 400         */
 401         if (res == zero) return(res);
 402         scale = 1.0;
 403         for (i = 0; i < dd->size; i++) {
 404             if (vars[i]) {
 405                 if (dd->perm[i] > topP && (unsigned) dd->perm[i] < topV) {
 406                     scale *= 2;
 407                 }
 408             }
 409         }
 410         if (scale > 1.0) {
 411             cuddRef(res);
 412             add_scale = cuddUniqueConst(dd,(CUDD_VALUE_TYPE)scale);
 413             if (add_scale == NULL) {
 414                 Cudd_RecursiveDeref(dd, res);
 415                 return(NULL);
 416             }
 417             cuddRef(add_scale);
 418             scaled = cuddAddApplyRecur(dd,Cudd_addTimes,res,add_scale);
 419             if (scaled == NULL) {
 420                 Cudd_RecursiveDeref(dd, add_scale);
 421                 Cudd_RecursiveDeref(dd, res);
 422                 return(NULL);
 423             }
 424             cuddRef(scaled);
 425             Cudd_RecursiveDeref(dd, add_scale);
 426             Cudd_RecursiveDeref(dd, res);
 427             res = scaled;
 428             cuddDeref(res);
 429         }
 430         return(res);
 431     }
 432
 433     /* compute the cofactors */
 434     if (topV == topA) {
 435         At = cuddT(A);
 436         Ae = cuddE(A);
 437     } else {
 438         At = Ae = A;
 439     }
 440     if (topV == topB) {
 441         Bt = cuddT(B);
 442         Be = cuddE(B);
 443     } else {
 444         Bt = Be = B;
 445     }
 446
 447     t = addMMRecur(dd, At, Bt, (int)topV, vars);
 448     if (t == NULL) return(NULL);
 449     cuddRef(t);
 450     e = addMMRecur(dd, Ae, Be, (int)topV, vars);
 451     if (e == NULL) {
 452         Cudd_RecursiveDeref(dd, t);
 453         return(NULL);
 454     }
 455     cuddRef(e);
 456
 457     index = dd->invperm[topV];
 458     if (vars[index] == 0) {
 459         /* We have split on either the rows of A or the columns
 460         ** of B. We just need to connect the two subresults,
 461         ** which correspond to two submatrices of the result.
 462         */
 463         res = (t == e) ? t : cuddUniqueInter(dd,index,t,e);
 464         if (res == NULL) {
 465             Cudd_RecursiveDeref(dd, t);
 466             Cudd_RecursiveDeref(dd, e);
 467             return(NULL);
 468         }
 469         cuddRef(res);
 470         cuddDeref(t);
 471         cuddDeref(e);
 472     } else {
 473         /* we have simultaneously split on the columns of A and
 474         ** the rows of B. The two subresults must be added.
 475         */
 476         res = cuddAddApplyRecur(dd,Cudd_addPlus,t,e);
 477         if (res == NULL) {
 478             Cudd_RecursiveDeref(dd, t);
 479             Cudd_RecursiveDeref(dd, e);
 480             return(NULL);
 481         }
 482         cuddRef(res);
 483         Cudd_RecursiveDeref(dd, t);
 484         Cudd_RecursiveDeref(dd, e);
 485     }
 486
 487     cuddCacheInsert2(dd,cacheOp,A,B,res);
 488
 489     /* We have computed (and stored in the computed table) a minimal
 490     ** result; that is, a result that assumes no summation variables
 491     ** between the current depth of the recursion and its top
 492     ** variable. We now take into account the z variables by properly
 493     ** scaling the result.
 494     */
 495     if (res != zero) {
 496         scale = 1.0;
 497         for (i = 0; i < dd->size; i++) {
 498             if (vars[i]) {
 499                 if (dd->perm[i] > topP && (unsigned) dd->perm[i] < topV) {
 500                     scale *= 2;
 501                 }
 502             }
 503         }
 504         if (scale > 1.0) {
 505             add_scale = cuddUniqueConst(dd,(CUDD_VALUE_TYPE)scale);
 506             if (add_scale == NULL) {
 507                 Cudd_RecursiveDeref(dd, res);
 508                 return(NULL);
 509             }
 510             cuddRef(add_scale);
 511             scaled = cuddAddApplyRecur(dd,Cudd_addTimes,res,add_scale);
 512             if (scaled == NULL) {
 513                 Cudd_RecursiveDeref(dd, res);
 514                 Cudd_RecursiveDeref(dd, add_scale);
 515                 return(NULL);
 516             }
 517             cuddRef(scaled);
 518             Cudd_RecursiveDeref(dd, add_scale);
 519             Cudd_RecursiveDeref(dd, res);
 520             res = scaled;
 521         }
 522     }
 523     cuddDeref(res);
 524     return(res);
 525
 526 } /* end of addMMRecur */
 527
 528
 529 /**Function********************************************************************
 530
 531   Synopsis    [Performs the recursive step of Cudd_addTriangle.]
 532
 533   Description [Performs the recursive step of Cudd_addTriangle. Returns
 534   a pointer to the result if successful; NULL otherwise.]
 535
 536   SideEffects [None]
 537
 538 ******************************************************************************/
 539 static DdNode *
 540 addTriangleRecur(
 541   DdManager * dd,
 542   DdNode * f,
 543   DdNode * g,
 544   int * vars,
 545   DdNode *cube)
 546 {
 547     DdNode *fv, *fvn, *gv, *gvn, *t, *e, *res;
 548     CUDD_VALUE_TYPE value;
 549     int top, topf, topg, index;
 550
 551     statLine(dd);
 552     if (f == DD_PLUS_INFINITY(dd) || g == DD_PLUS_INFINITY(dd)) {
 553         return(DD_PLUS_INFINITY(dd));
 554     }
 555
 556     if (cuddIsConstant(f) && cuddIsConstant(g)) {
 557         value = cuddV(f) + cuddV(g);
 558         res = cuddUniqueConst(dd, value);
 559         return(res);
 560     }
 561     if (f < g) {
 562         DdNode *tmp = f;
 563         f = g;
 564         g = tmp;
 565     }
 566
 567     if (f->ref != 1 || g->ref != 1) {
 568         res = cuddCacheLookup(dd, DD_ADD_TRIANGLE_TAG, f, g, cube);
 569         if (res != NULL) {
 570             return(res);
 571         }
 572     }
 573
 574     topf = cuddI(dd,f->index); topg = cuddI(dd,g->index);
 575     top = ddMin(topf,topg);
 576
 577     if (top == topf) {fv = cuddT(f); fvn = cuddE(f);} else {fv = fvn = f;}
 578     if (top == topg) {gv = cuddT(g); gvn = cuddE(g);} else {gv = gvn = g;}
 579
 580     t = addTriangleRecur(dd, fv, gv, vars, cube);
 581     if (t == NULL) return(NULL);
 582     cuddRef(t);
 583     e = addTriangleRecur(dd, fvn, gvn, vars, cube);
 584     if (e == NULL) {
 585         Cudd_RecursiveDeref(dd, t);
 586         return(NULL);
 587     }
 588     cuddRef(e);
 589
 590     index = dd->invperm[top];
 591     if (vars[index] < 0) {
 592         res = (t == e) ? t : cuddUniqueInter(dd,index,t,e);
 593         if (res == NULL) {
 594             Cudd_RecursiveDeref(dd, t);
 595             Cudd_RecursiveDeref(dd, e);
 596             return(NULL);
 597         }
 598         cuddDeref(t);
 599         cuddDeref(e);
 600     } else {
 601         res = cuddAddApplyRecur(dd,Cudd_addMinimum,t,e);
 602         if (res == NULL) {
 603             Cudd_RecursiveDeref(dd, t);
 604             Cudd_RecursiveDeref(dd, e);
 605             return(NULL);
 606         }
 607         cuddRef(res);
 608         Cudd_RecursiveDeref(dd, t);
 609         Cudd_RecursiveDeref(dd, e);
 610         cuddDeref(res);
 611     }
 612
 613     if (f->ref != 1 || g->ref != 1) {
 614         cuddCacheInsert(dd, DD_ADD_TRIANGLE_TAG, f, g, cube, res);
 615     }
 616
 617     return(res);
 618
 619 } /* end of addTriangleRecur */
 620
 621
 622 /**Function********************************************************************
 623
 624   Synopsis    [Performs the recursive step of Cudd_addOuterSum.]
 625
 626   Description [Performs the recursive step of Cudd_addOuterSum.
 627   Returns a pointer to the result if successful; NULL otherwise.]
 628
 629   SideEffects [None]
 630
 631   SeeAlso     []
 632
 633 ******************************************************************************/
 634 static DdNode *
 635 cuddAddOuterSumRecur(
 636   DdManager *dd,
 637   DdNode *M,
 638   DdNode *r,
 639   DdNode *c)
 640 {
 641     DdNode *P, *R, *Mt, *Me, *rt, *re, *ct, *ce, *Rt, *Re;
 642     int topM, topc, topr;
 643     int v, index;
 644
 645     statLine(dd);
 646     /* Check special cases. */
 647     if (r == DD_PLUS_INFINITY(dd) || c == DD_PLUS_INFINITY(dd)) return(M);
 648
 649     if (cuddIsConstant(c) && cuddIsConstant(r)) {
 650         R = cuddUniqueConst(dd,Cudd_V(c)+Cudd_V(r));
 651         cuddRef(R);
 652         if (cuddIsConstant(M)) {
 653             if (cuddV(R) <= cuddV(M)) {
 654                 cuddDeref(R);
 655                 return(R);
 656             } else {
 657                 Cudd_RecursiveDeref(dd,R);
 658                 return(M);
 659             }
 660         } else {
 661             P = Cudd_addApply(dd,Cudd_addMinimum,R,M);
 662             cuddRef(P);
 663             Cudd_RecursiveDeref(dd,R);
 664             cuddDeref(P);
 665             return(P);
 666         }
 667     }
 668
 669     /* Check the cache. */
 670     R = cuddCacheLookup(dd,DD_ADD_OUT_SUM_TAG,M,r,c);
 671     if (R != NULL) return(R);
 672
 673     topM = cuddI(dd,M->index); topr = cuddI(dd,r->index);
 674     topc = cuddI(dd,c->index);
 675     v = ddMin(topM,ddMin(topr,topc));
 676
 677     /* Compute cofactors. */
 678     if (topM == v) { Mt = cuddT(M); Me = cuddE(M); } else { Mt = Me = M; }
 679     if (topr == v) { rt = cuddT(r); re = cuddE(r); } else { rt = re = r; }
 680     if (topc == v) { ct = cuddT(c); ce = cuddE(c); } else { ct = ce = c; }
 681
 682     /* Recursively solve. */
 683     Rt = cuddAddOuterSumRecur(dd,Mt,rt,ct);
 684     if (Rt == NULL) return(NULL);
 685     cuddRef(Rt);
 686     Re = cuddAddOuterSumRecur(dd,Me,re,ce);
 687     if (Re == NULL) {
 688         Cudd_RecursiveDeref(dd, Rt);
 689         return(NULL);
 690     }
 691     cuddRef(Re);
 692     index = dd->invperm[v];
 693     R = (Rt == Re) ? Rt : cuddUniqueInter(dd,index,Rt,Re);
 694     if (R == NULL) {
 695         Cudd_RecursiveDeref(dd, Rt);
 696         Cudd_RecursiveDeref(dd, Re);
 697         return(NULL);
 698     }
 699     cuddDeref(Rt);
 700     cuddDeref(Re);
 701
 702     /* Store the result in the cache. */
 703     cuddCacheInsert(dd,DD_ADD_OUT_SUM_TAG,M,r,c,R);
 704
 705     return(R);
 706
 707 } /* end of cuddAddOuterSumRecur */