Context Navigation

source: vis_dev/glu-2.1/src/cuBdd/cuddMatMult.c @ 10

Last change on this file since 10 was 8, checked in by cecile, 13 years ago
src glu
File size: 20.0 KB

Line
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, AB != BA
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 */

Note: See TracBrowser for help on using the repository browser.

Download in other formats: