1 /* ========================================================================== */
2 /* ssmult_transpose */
3 /* ========================================================================== */
4
5 /* C = A' or A.' where the input matrix A may have unsorted columns. The output
6 C is always returned with sorted columns.
7 */
8
9 #include "ssmult.h"
10
ssdump(const mxArray * A)11 void ssdump (const mxArray *A)
12 {
13 Int j, p, m, n, *Ap, *Ai, anz ;
14 double *Ax, *Az ;
15
16 m = mxGetM (A) ;
17 n = mxGetN (A) ;
18 Ap = (Int *) mxGetJc (A) ;
19 Ai = (Int *) mxGetIr (A) ;
20 Ax = mxGetPr (A) ;
21 Az = mxGetPi (A) ;
22 anz = Ap [n] ;
23
24 printf ("%d by %d\n", m, n) ;
25 for (j = 0 ; j < n ; j++)
26 {
27 printf ("column %d:\n", j) ;
28 for (p = Ap [j] ; p < Ap [j+1] ; p++)
29 {
30 printf ("A(%d,%d) = %g", Ai [p], j, Ax [p]) ;
31 if (Az != NULL) printf (" %g ", Az [p]) ;
32 printf ("\n") ;
33 }
34 }
35 printf ("\n") ;
36 }
37
38
ssmult_transpose(const mxArray * A,int conj)39 mxArray *ssmult_transpose /* returns C = A' or A.' */
40 (
41 const mxArray *A,
42 int conj /* compute A' if true, compute A.' if false */
43 )
44 {
45 Int *Cp, *Ci, *Ap, *Ai, *W ;
46 double *Cx, *Cz, *Ax, *Az ;
47 mxArray *C ;
48 Int p, pend, q, i, j, n, m, anz, cnz ;
49 int C_is_complex ;
50
51 /* ---------------------------------------------------------------------- */
52 /* get inputs */
53 /* ---------------------------------------------------------------------- */
54
55 m = mxGetM (A) ;
56 n = mxGetN (A) ;
57 Ap = (Int *) mxGetJc (A) ;
58 Ai = (Int *) mxGetIr (A) ;
59 Ax = mxGetPr (A) ;
60 Az = mxGetPi (A) ;
61 anz = Ap [n] ;
62 C_is_complex = mxIsComplex (A) ;
63
64 /* ---------------------------------------------------------------------- */
65 /* allocate C but do not initialize it */
66 /* ---------------------------------------------------------------------- */
67
68 cnz = MAX (anz, 1) ;
69 C = mxCreateSparse (0, 0, 0, C_is_complex ? mxCOMPLEX : mxREAL) ;
70 mxFree (mxGetJc (C)) ;
71 mxFree (mxGetIr (C)) ;
72 mxFree (mxGetPr (C)) ;
73 mxFree (mxGetPi (C)) ;
74 Cp = mxMalloc ((m+1) * sizeof (Int)) ;
75 Ci = mxMalloc (cnz * sizeof (Int)) ;
76 Cx = mxMalloc (cnz * sizeof (double)) ;
77 Cz = C_is_complex ? mxMalloc (cnz * sizeof (double)) : NULL ;
78 mxSetJc (C, (mwIndex *) Cp) ;
79 mxSetIr (C, (mwIndex *) Ci) ;
80 mxSetPr (C, Cx) ;
81 mxSetPi (C, Cz) ;
82 mxSetNzmax (C, cnz) ;
83 mxSetM (C, n) ;
84 mxSetN (C, m) ;
85
86 /* ---------------------------------------------------------------------- */
87 /* allocate workspace */
88 /* ---------------------------------------------------------------------- */
89
90 W = mxCalloc (m, sizeof (Int)) ;
91
92 /* ---------------------------------------------------------------------- */
93 /* compute row counts */
94 /* ---------------------------------------------------------------------- */
95
96 for (p = 0 ; p < anz ; p++)
97 {
98 W [Ai [p]]++ ;
99 }
100
101 /* ---------------------------------------------------------------------- */
102 /* compute column pointers of C and copy back into W */
103 /* ---------------------------------------------------------------------- */
104
105 for (p = 0, i = 0 ; i < m ; i++)
106 {
107 Cp [i] = p ;
108 p += W [i] ;
109 W [i] = Cp [i] ;
110 }
111 Cp [m] = p ;
112
113 /* ---------------------------------------------------------------------- */
114 /* C = A' */
115 /* ---------------------------------------------------------------------- */
116
117 p = 0 ;
118 if (!C_is_complex)
119 {
120 /* C = A' (real case) */
121 for (j = 0 ; j < n ; j++)
122 {
123 pend = Ap [j+1] ;
124 for ( ; p < pend ; p++)
125 {
126 q = W [Ai [p]]++ ; /* find position for C(j,i) */
127 Ci [q] = j ; /* place A(i,j) as entry C(j,i) */
128 Cx [q] = Ax [p] ;
129 }
130 }
131 }
132 else if (conj)
133 {
134 /* C = A' (complex conjugate) */
135 for (j = 0 ; j < n ; j++)
136 {
137 pend = Ap [j+1] ;
138 for ( ; p < pend ; p++)
139 {
140 q = W [Ai [p]]++ ; /* find position for C(j,i) */
141 Ci [q] = j ; /* place A(i,j) as entry C(j,i) */
142 Cx [q] = Ax [p] ;
143 Cz [q] = -Az [p] ;
144 }
145 }
146 }
147 else
148 {
149 /* C = A.' (complex case) */
150 for (j = 0 ; j < n ; j++)
151 {
152 pend = Ap [j+1] ;
153 for ( ; p < pend ; p++)
154 {
155 q = W [Ai [p]]++ ; /* find position for C(j,i) */
156 Ci [q] = j ; /* place A(i,j) as entry C(j,i) */
157 Cx [q] = Ax [p] ;
158 Cz [q] = Az [p] ;
159 }
160 }
161 }
162
163 /* ---------------------------------------------------------------------- */
164 /* free workspace and return result */
165 /* ---------------------------------------------------------------------- */
166
167 mxFree (W) ;
168 return (C) ;
169 }
170