1 /* ========================================================================== */
2 /* === ssmultsym ============================================================ */
3 /* ========================================================================== */
4
5 /* s = ssmultsym (A,B) computes nnz(A*B), and the flops and memory required to
6 * compute it, where A and B are both sparse. Either A or B, or both, can be
7 * complex. Memory usage includes C itself, and workspace. If C is m-by-n,
8 * ssmultsym requires only 4*m bytes for 32-bit MATLAB, 8*m for 64-bit MATLAB.
9 *
10 * Copyright 2007-2009, Timothy A. Davis, http://www.suitesparse.com
11 */
12
13 #include "ssmult.h"
14
15 /* -------------------------------------------------------------------------- */
16 /* ssmultsym mexFunction */
17 /* -------------------------------------------------------------------------- */
18
mexFunction(int nargout,mxArray * pargout[],int nargin,const mxArray * pargin[])19 void mexFunction
20 (
21 int nargout,
22 mxArray *pargout [ ],
23 int nargin,
24 const mxArray *pargin [ ]
25 )
26 {
27 double cnz, flops, mem, e, multadds ;
28 Int *Ap, *Ai, *Bp, *Bi, *Flag ;
29 Int Anrow, Ancol, Bnrow, Bncol, i, j, k, pb, pa, pbend, paend, mark,
30 A_is_complex, B_is_complex, C_is_complex, pbstart, cjnz ;
31 static const char *snames [ ] =
32 {
33 "nz", /* # of nonzeros in C=A*B */
34 "flops", /* flop count required to compute C=A*B */
35 "memory" /* memory requirement in bytes */
36 } ;
37
38 /* ---------------------------------------------------------------------- */
39 /* get inputs and workspace */
40 /* ---------------------------------------------------------------------- */
41
42 if (nargin != 2 || nargout > 1)
43 {
44 mexErrMsgTxt ("Usage: s = ssmultsym (A,B)") ;
45 }
46
47 Ap = (Int *) mxGetJc (pargin [0]) ;
48 Ai = (Int *) mxGetIr (pargin [0]) ;
49 Anrow = mxGetM (pargin [0]) ;
50 Ancol = mxGetN (pargin [0]) ;
51 A_is_complex = mxIsComplex (pargin [0]) ;
52
53 Bp = (Int *) mxGetJc (pargin [1]) ;
54 Bi = (Int *) mxGetIr (pargin [1]) ;
55 Bnrow = mxGetM (pargin [1]) ;
56 Bncol = mxGetN (pargin [1]) ;
57 B_is_complex = mxIsComplex (pargin [1]) ;
58
59 if (Ancol != Bnrow || !mxIsSparse (pargin [0]) || !mxIsSparse (pargin [1]))
60 {
61 mexErrMsgTxt ("wrong dimensions, or A and B not sparse") ;
62 }
63
64 Flag = mxCalloc (Anrow, sizeof (Int)) ; /* workspace */
65
66 /* ---------------------------------------------------------------------- */
67 /* compute # of nonzeros in result, flop count, and memory */
68 /* ---------------------------------------------------------------------- */
69
70 pb = 0 ;
71 cnz = 0 ;
72 multadds = 0 ;
73 for (j = 0 ; j < Bncol ; j++) /* compute C(:,j) */
74 {
75 mark = j+1 ;
76 pbstart = Bp [j] ;
77 pbend = Bp [j+1] ;
78 cjnz = 0 ;
79 for ( ; pb < pbend ; pb++)
80 {
81 k = Bi [pb] ; /* nonzero entry B(k,j) */
82 pa = Ap [k] ;
83 paend = Ap [k+1] ;
84 multadds += (paend - pa) ; /* one mult-add per entry in A(:,k) */
85 if (cjnz == Anrow)
86 {
87 /* C(:,j) is already completely dense; no need to scan A.
88 * Continue scanning B(:,j) to compute flop count. */
89 continue ;
90 }
91 for ( ; pa < paend ; pa++)
92 {
93 i = Ai [pa] ; /* nonzero entry A(i,k) */
94 if (Flag [i] != mark)
95 {
96 /* C(i,j) is a new nonzero */
97 Flag [i] = mark ; /* mark i as appearing in C(:,j) */
98 cjnz++ ;
99 }
100 }
101 }
102 cnz += cjnz ;
103 }
104
105 C_is_complex = A_is_complex || B_is_complex ;
106 e = (C_is_complex ? 2 : 1) ;
107
108 mem =
109 Anrow * sizeof (Int) /* Flag */
110 + e * Anrow * sizeof (double) /* W */
111 + (Bncol+1) * sizeof (Int) /* Cp */
112 + cnz * sizeof (Int) /* Ci */
113 + e * cnz * sizeof (double) ; /* Cx and Cx */
114
115 if (A_is_complex)
116 {
117 if (B_is_complex)
118 {
119 /* all of C, A, and B are complex */
120 flops = 8 * multadds - 2 * cnz ;
121 }
122 else
123 {
124 /* C and A are complex, B is real */
125 flops = 4 * multadds - 2 * cnz ;
126 }
127 }
128 else
129 {
130 if (B_is_complex)
131 {
132 /* C and B are complex, A is real */
133 flops = 4 * multadds - 2 * cnz ;
134 }
135 else
136 {
137 /* all of C, A, and B are real */
138 flops = 2 * multadds - cnz ;
139 }
140 }
141
142 /* ---------------------------------------------------------------------- */
143 /* return result */
144 /* ---------------------------------------------------------------------- */
145
146 pargout [0] = mxCreateStructMatrix (1, 1, 3, snames) ;
147 mxSetFieldByNumber (pargout [0], 0, 0, mxCreateDoubleScalar (cnz)) ;
148 mxSetFieldByNumber (pargout [0], 0, 1, mxCreateDoubleScalar (flops)) ;
149 mxSetFieldByNumber (pargout [0], 0, 2, mxCreateDoubleScalar (mem)) ;
150 }
151