1 /**CFile***********************************************************************
3 FileName [cuddMatMult.c]
7 Synopsis [Matrix multiplication functions.]
9 Description [External procedures included in this module:
11 <li> Cudd_addMatrixMultiply()
12 <li> Cudd_addTimesPlus()
13 <li> Cudd_addTriangle()
14 <li> Cudd_addOuterSum()
16 Static procedures included in this module:
19 <li> addTriangleRecur()
20 <li> cuddAddOuterSumRecur()
23 Author [Fabio Somenzi]
25 Copyright [Copyright (c) 1995-2004, Regents of the University of Colorado
29 Redistribution and use in source and binary forms, with or without
30 modification, are permitted provided that the following conditions
33 Redistributions of source code must retain the above copyright
34 notice, this list of conditions and the following disclaimer.
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.
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.
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.]
57 ******************************************************************************/
63 /*---------------------------------------------------------------------------*/
64 /* Constant declarations */
65 /*---------------------------------------------------------------------------*/
68 /*---------------------------------------------------------------------------*/
69 /* Stucture declarations */
70 /*---------------------------------------------------------------------------*/
73 /*---------------------------------------------------------------------------*/
74 /* Type declarations */
75 /*---------------------------------------------------------------------------*/
78 /*---------------------------------------------------------------------------*/
79 /* Variable declarations */
80 /*---------------------------------------------------------------------------*/
83 static char rcsid
[] DD_UNUSED
= "$Id: cuddMatMult.c,v 1.17 2004/08/13 18:04:50 fabio Exp $";
86 /*---------------------------------------------------------------------------*/
87 /* Macro declarations */
88 /*---------------------------------------------------------------------------*/
91 /**AutomaticStart*************************************************************/
93 /*---------------------------------------------------------------------------*/
94 /* Static function prototypes */
95 /*---------------------------------------------------------------------------*/
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
);
101 /**AutomaticEnd***************************************************************/
104 /*---------------------------------------------------------------------------*/
105 /* Definition of exported functions */
106 /*---------------------------------------------------------------------------*/
108 /**Function********************************************************************
110 Synopsis [Calculates the product of two matrices represented as
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.]
124 SeeAlso [Cudd_addTimesPlus Cudd_addTriangle Cudd_bddAndAbstract]
126 ******************************************************************************/
128 Cudd_addMatrixMultiply(
138 /* Array vars says what variables are "summation" variables. */
140 vars
= ALLOC(int,nvars
);
142 dd
->errorCode
= CUDD_MEMORY_OUT
;
145 for (i
= 0; i
< nvars
; i
++) {
148 for (i
= 0; i
< nz
; i
++) {
149 vars
[z
[i
]->index
] = 1;
154 res
= addMMRecur(dd
,A
,B
,-1,vars
);
155 } while (dd
->reordered
== 1);
159 } /* end of Cudd_addMatrixMultiply */
162 /**Function********************************************************************
164 Synopsis [Calculates the product of two matrices represented as
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.]
177 SeeAlso [Cudd_addMatrixMultiply]
179 ******************************************************************************/
188 DdNode
*w
, *cube
, *tmp
, *res
;
190 tmp
= Cudd_addApply(dd
,Cudd_addTimes
,A
,B
);
191 if (tmp
== NULL
) return(NULL
);
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
));
197 Cudd_RecursiveDeref(dd
,tmp
);
201 Cudd_RecursiveDeref(dd
,cube
);
204 res
= Cudd_addExistAbstract(dd
,tmp
,cube
);
206 Cudd_RecursiveDeref(dd
,tmp
);
207 Cudd_RecursiveDeref(dd
,cube
);
211 Cudd_RecursiveDeref(dd
,cube
);
212 Cudd_RecursiveDeref(dd
,tmp
);
216 } /* end of Cudd_addTimesPlus */
219 /**Function********************************************************************
221 Synopsis [Performs the triangulation step for the shortest path
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. ]
235 SeeAlso [Cudd_addMatrixMultiply Cudd_bddAndAbstract]
237 ******************************************************************************/
250 vars
= ALLOC(int, nvars
);
252 dd
->errorCode
= CUDD_MEMORY_OUT
;
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
);
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
);
274 } /* end of Cudd_addTriangle */
277 /**Function********************************************************************
279 Synopsis [Takes the minimum of a matrix and the outer sum of two vectors.]
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.]
290 ******************************************************************************/
302 res
= cuddAddOuterSumRecur(dd
, M
, r
, c
);
303 } while (dd
->reordered
== 1);
306 } /* end of Cudd_addOuterSum */
309 /*---------------------------------------------------------------------------*/
310 /* Definition of internal functions */
311 /*---------------------------------------------------------------------------*/
314 /*---------------------------------------------------------------------------*/
315 /* Definition of static functions */
316 /*---------------------------------------------------------------------------*/
318 /**Function********************************************************************
320 Synopsis [Performs the recursive step of Cudd_addMatrixMultiply.]
322 Description [Performs the recursive step of Cudd_addMatrixMultiply.
323 Returns a pointer to the result if successful; NULL otherwise.]
327 ******************************************************************************/
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 */
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
;
356 if (A
== zero
|| B
== zero
) {
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.
366 value
= cuddV(A
) * cuddV(B
);
367 for (i
= 0; i
< dd
->size
; i
++) {
369 if (dd
->perm
[i
] > topP
) {
370 value
*= (CUDD_VALUE_TYPE
) 2;
374 res
= cuddUniqueConst(dd
, value
);
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.
389 topA
= cuddI(dd
,A
->index
); topB
= cuddI(dd
,B
->index
);
390 topV
= ddMin(topA
,topB
);
392 cacheOp
= (DD_CTFP
) addMMRecur
;
393 res
= cuddCacheLookup2(dd
,cacheOp
,A
,B
);
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
401 if (res
== zero
) return(res
);
403 for (i
= 0; i
< dd
->size
; i
++) {
405 if (dd
->perm
[i
] > topP
&& (unsigned) dd
->perm
[i
] < topV
) {
412 add_scale
= cuddUniqueConst(dd
,(CUDD_VALUE_TYPE
)scale
);
413 if (add_scale
== NULL
) {
414 Cudd_RecursiveDeref(dd
, res
);
418 scaled
= cuddAddApplyRecur(dd
,Cudd_addTimes
,res
,add_scale
);
419 if (scaled
== NULL
) {
420 Cudd_RecursiveDeref(dd
, add_scale
);
421 Cudd_RecursiveDeref(dd
, res
);
425 Cudd_RecursiveDeref(dd
, add_scale
);
426 Cudd_RecursiveDeref(dd
, res
);
433 /* compute the cofactors */
447 t
= addMMRecur(dd
, At
, Bt
, (int)topV
, vars
);
448 if (t
== NULL
) return(NULL
);
450 e
= addMMRecur(dd
, Ae
, Be
, (int)topV
, vars
);
452 Cudd_RecursiveDeref(dd
, t
);
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.
463 res
= (t
== e
) ? t
: cuddUniqueInter(dd
,index
,t
,e
);
465 Cudd_RecursiveDeref(dd
, t
);
466 Cudd_RecursiveDeref(dd
, e
);
473 /* we have simultaneously split on the columns of A and
474 ** the rows of B. The two subresults must be added.
476 res
= cuddAddApplyRecur(dd
,Cudd_addPlus
,t
,e
);
478 Cudd_RecursiveDeref(dd
, t
);
479 Cudd_RecursiveDeref(dd
, e
);
483 Cudd_RecursiveDeref(dd
, t
);
484 Cudd_RecursiveDeref(dd
, e
);
487 cuddCacheInsert2(dd
,cacheOp
,A
,B
,res
);
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.
497 for (i
= 0; i
< dd
->size
; i
++) {
499 if (dd
->perm
[i
] > topP
&& (unsigned) dd
->perm
[i
] < topV
) {
505 add_scale
= cuddUniqueConst(dd
,(CUDD_VALUE_TYPE
)scale
);
506 if (add_scale
== NULL
) {
507 Cudd_RecursiveDeref(dd
, res
);
511 scaled
= cuddAddApplyRecur(dd
,Cudd_addTimes
,res
,add_scale
);
512 if (scaled
== NULL
) {
513 Cudd_RecursiveDeref(dd
, res
);
514 Cudd_RecursiveDeref(dd
, add_scale
);
518 Cudd_RecursiveDeref(dd
, add_scale
);
519 Cudd_RecursiveDeref(dd
, res
);
526 } /* end of addMMRecur */
529 /**Function********************************************************************
531 Synopsis [Performs the recursive step of Cudd_addTriangle.]
533 Description [Performs the recursive step of Cudd_addTriangle. Returns
534 a pointer to the result if successful; NULL otherwise.]
538 ******************************************************************************/
547 DdNode
*fv
, *fvn
, *gv
, *gvn
, *t
, *e
, *res
;
548 CUDD_VALUE_TYPE value
;
549 int top
, topf
, topg
, index
;
552 if (f
== DD_PLUS_INFINITY(dd
) || g
== DD_PLUS_INFINITY(dd
)) {
553 return(DD_PLUS_INFINITY(dd
));
556 if (cuddIsConstant(f
) && cuddIsConstant(g
)) {
557 value
= cuddV(f
) + cuddV(g
);
558 res
= cuddUniqueConst(dd
, value
);
567 if (f
->ref
!= 1 || g
->ref
!= 1) {
568 res
= cuddCacheLookup(dd
, DD_ADD_TRIANGLE_TAG
, f
, g
, cube
);
574 topf
= cuddI(dd
,f
->index
); topg
= cuddI(dd
,g
->index
);
575 top
= ddMin(topf
,topg
);
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
;}
580 t
= addTriangleRecur(dd
, fv
, gv
, vars
, cube
);
581 if (t
== NULL
) return(NULL
);
583 e
= addTriangleRecur(dd
, fvn
, gvn
, vars
, cube
);
585 Cudd_RecursiveDeref(dd
, t
);
590 index
= dd
->invperm
[top
];
591 if (vars
[index
] < 0) {
592 res
= (t
== e
) ? t
: cuddUniqueInter(dd
,index
,t
,e
);
594 Cudd_RecursiveDeref(dd
, t
);
595 Cudd_RecursiveDeref(dd
, e
);
601 res
= cuddAddApplyRecur(dd
,Cudd_addMinimum
,t
,e
);
603 Cudd_RecursiveDeref(dd
, t
);
604 Cudd_RecursiveDeref(dd
, e
);
608 Cudd_RecursiveDeref(dd
, t
);
609 Cudd_RecursiveDeref(dd
, e
);
613 if (f
->ref
!= 1 || g
->ref
!= 1) {
614 cuddCacheInsert(dd
, DD_ADD_TRIANGLE_TAG
, f
, g
, cube
, res
);
619 } /* end of addTriangleRecur */
622 /**Function********************************************************************
624 Synopsis [Performs the recursive step of Cudd_addOuterSum.]
626 Description [Performs the recursive step of Cudd_addOuterSum.
627 Returns a pointer to the result if successful; NULL otherwise.]
633 ******************************************************************************/
635 cuddAddOuterSumRecur(
641 DdNode
*P
, *R
, *Mt
, *Me
, *rt
, *re
, *ct
, *ce
, *Rt
, *Re
;
642 int topM
, topc
, topr
;
646 /* Check special cases. */
647 if (r
== DD_PLUS_INFINITY(dd
) || c
== DD_PLUS_INFINITY(dd
)) return(M
);
649 if (cuddIsConstant(c
) && cuddIsConstant(r
)) {
650 R
= cuddUniqueConst(dd
,Cudd_V(c
)+Cudd_V(r
));
652 if (cuddIsConstant(M
)) {
653 if (cuddV(R
) <= cuddV(M
)) {
657 Cudd_RecursiveDeref(dd
,R
);
661 P
= Cudd_addApply(dd
,Cudd_addMinimum
,R
,M
);
663 Cudd_RecursiveDeref(dd
,R
);
669 /* Check the cache. */
670 R
= cuddCacheLookup(dd
,DD_ADD_OUT_SUM_TAG
,M
,r
,c
);
671 if (R
!= NULL
) return(R
);
673 topM
= cuddI(dd
,M
->index
); topr
= cuddI(dd
,r
->index
);
674 topc
= cuddI(dd
,c
->index
);
675 v
= ddMin(topM
,ddMin(topr
,topc
));
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
; }
682 /* Recursively solve. */
683 Rt
= cuddAddOuterSumRecur(dd
,Mt
,rt
,ct
);
684 if (Rt
== NULL
) return(NULL
);
686 Re
= cuddAddOuterSumRecur(dd
,Me
,re
,ce
);
688 Cudd_RecursiveDeref(dd
, Rt
);
692 index
= dd
->invperm
[v
];
693 R
= (Rt
== Re
) ? Rt
: cuddUniqueInter(dd
,index
,Rt
,Re
);
695 Cudd_RecursiveDeref(dd
, Rt
);
696 Cudd_RecursiveDeref(dd
, Re
);
702 /* Store the result in the cache. */
703 cuddCacheInsert(dd
,DD_ADD_OUT_SUM_TAG
,M
,r
,c
,R
);
707 } /* end of cuddAddOuterSumRecur */