1 /* ========================================================================== */
2 /* === UMF_triplet ========================================================== */
3 /* ========================================================================== */
4
5 /* -------------------------------------------------------------------------- */
6 /* Copyright (c) 2005-2012 by Timothy A. Davis, http://www.suitesparse.com. */
7 /* All Rights Reserved. See ../Doc/License.txt for License. */
8 /* -------------------------------------------------------------------------- */
9
10 /*
11 Not user callable. Converts triplet input to column-oriented form.
12 Duplicate entries may exist (they are summed in the output). The columns
13 of the column-oriented form are in sorted order. The input is not modified.
14 Returns 1 if OK, 0 if an error occurred.
15
16 Compiled into four different routines for each version (di, dl, zi, zl),
17 for a total of 16 different routines.
18 */
19
20 #include "umf_internal.h"
21 #include "umf_triplet.h"
22
23 #ifdef DO_MAP
24 #ifdef DO_VALUES
UMF_triplet_map_x(Int n_row,Int n_col,Int nz,const Int Ti[],const Int Tj[],Int Ap[],Int Ai[],Int Rp[],Int Rj[],Int W[],Int RowCount[],const double Tx[],double Ax[],double Rx[],const double Tz[],double Az[],double Rz[],Int Map[],Int Map2[])25 GLOBAL Int UMF_triplet_map_x
26 #else
27 GLOBAL Int UMF_triplet_map_nox
28 #endif
29 #else
30 #ifdef DO_VALUES
31 GLOBAL Int UMF_triplet_nomap_x
32 #else
33 GLOBAL Int UMF_triplet_nomap_nox
34 #endif
35 #endif
36 (
37 Int n_row,
38 Int n_col,
39 Int nz,
40 const Int Ti [ ], /* size nz */
41 const Int Tj [ ], /* size nz */
42 Int Ap [ ], /* size n_col + 1 */
43 Int Ai [ ], /* size nz */
44 Int Rp [ ], /* size n_row + 1 */
45 Int Rj [ ], /* size nz */
46 Int W [ ], /* size max (n_row, n_col) */
47 Int RowCount [ ] /* size n_row */
48 #ifdef DO_VALUES
49 , const double Tx [ ] /* size nz */
50 , double Ax [ ] /* size nz */
51 , double Rx [ ] /* size nz */
52 #ifdef COMPLEX
53 , const double Tz [ ] /* size nz */
54 , double Az [ ] /* size nz */
55 , double Rz [ ] /* size nz */
56 #endif
57 #endif
58 #ifdef DO_MAP
59 , Int Map [ ] /* size nz */
60 , Int Map2 [ ] /* size nz */
61 #endif
62 )
63 {
64
65 /* ---------------------------------------------------------------------- */
66 /* local variables */
67 /* ---------------------------------------------------------------------- */
68
69 Int i, j, k, p, cp, p1, p2, pdest, pj ;
70 #ifdef DO_MAP
71 Int duplicates ;
72 #endif
73 #ifdef DO_VALUES
74 #ifdef COMPLEX
75 Int split = SPLIT (Tz) && SPLIT (Az) && SPLIT (Rz) ;
76 #endif
77 #endif
78
79 /* ---------------------------------------------------------------------- */
80 /* count the entries in each row (also counting duplicates) */
81 /* ---------------------------------------------------------------------- */
82
83 /* use W as workspace for row counts (including duplicates) */
84 for (i = 0 ; i < n_row ; i++)
85 {
86 W [i] = 0 ;
87 }
88
89 for (k = 0 ; k < nz ; k++)
90 {
91 i = Ti [k] ;
92 j = Tj [k] ;
93 if (i < 0 || i >= n_row || j < 0 || j >= n_col)
94 {
95 return (UMFPACK_ERROR_invalid_matrix) ;
96 }
97 W [i]++ ;
98 #ifndef NDEBUG
99 DEBUG1 ((ID " triplet: "ID" "ID" ", k, i, j)) ;
100 #ifdef DO_VALUES
101 {
102 Entry tt ;
103 ASSIGN (tt, Tx, Tz, k, split) ;
104 EDEBUG2 (tt) ;
105 DEBUG1 (("\n")) ;
106 }
107 #endif
108 #endif
109 }
110
111 /* ---------------------------------------------------------------------- */
112 /* compute the row pointers */
113 /* ---------------------------------------------------------------------- */
114
115 Rp [0] = 0 ;
116 for (i = 0 ; i < n_row ; i++)
117 {
118 Rp [i+1] = Rp [i] + W [i] ;
119 W [i] = Rp [i] ;
120 }
121
122 /* W is now equal to the row pointers */
123
124 /* ---------------------------------------------------------------------- */
125 /* construct the row form */
126 /* ---------------------------------------------------------------------- */
127
128 for (k = 0 ; k < nz ; k++)
129 {
130 p = W [Ti [k]]++ ;
131 #ifdef DO_MAP
132 Map [k] = p ;
133 #endif
134 Rj [p] = Tj [k] ;
135 #ifdef DO_VALUES
136 #ifdef COMPLEX
137 if (split)
138 {
139 Rx [p] = Tx [k] ;
140 Rz [p] = Tz [k] ;
141 }
142 else
143 {
144 Rx [2*p ] = Tx [2*k ] ;
145 Rx [2*p+1] = Tx [2*k+1] ;
146 }
147 #else
148 Rx [p] = Tx [k] ;
149 #endif
150 #endif
151 }
152
153 /* Rp stays the same, but W [i] is advanced to the start of row i+1 */
154
155 #ifndef NDEBUG
156 for (i = 0 ; i < n_row ; i++)
157 {
158 ASSERT (W [i] == Rp [i+1]) ;
159 }
160 #ifdef DO_MAP
161 for (k = 0 ; k < nz ; k++)
162 {
163 /* make sure that kth triplet is mapped correctly */
164 p = Map [k] ;
165 DEBUG1 (("First row map: Map ["ID"] = "ID"\n", k, p)) ;
166 i = Ti [k] ;
167 j = Tj [k] ;
168 ASSERT (j == Rj [p]) ;
169 ASSERT (Rp [i] <= p && p < Rp [i+1]) ;
170 }
171 #endif
172 #endif
173
174 /* ---------------------------------------------------------------------- */
175 /* sum up duplicates */
176 /* ---------------------------------------------------------------------- */
177
178 /* use W [j] to hold position in Ri/Rx/Rz of a_ij, for row i [ */
179
180 for (j = 0 ; j < n_col ; j++)
181 {
182 W [j] = EMPTY ;
183 }
184
185 #ifdef DO_MAP
186 duplicates = FALSE ;
187 #endif
188
189 for (i = 0 ; i < n_row ; i++)
190 {
191 p1 = Rp [i] ;
192 p2 = Rp [i+1] ;
193 pdest = p1 ;
194 /* At this point, W [j] < p1 holds true for all columns j, */
195 /* because Ri/Rx/Rz is stored in row oriented order. */
196 #ifndef NDEBUG
197 if (UMF_debug >= -2)
198 {
199 for (j = 0 ; j < n_col ; j++)
200 {
201 ASSERT (W [j] < p1) ;
202 }
203 }
204 #endif
205 for (p = p1 ; p < p2 ; p++)
206 {
207 j = Rj [p] ;
208 ASSERT (j >= 0 && j < n_col) ;
209 pj = W [j] ;
210 if (pj >= p1)
211 {
212 /* this column index, j, is already in row i, at position pj */
213 ASSERT (pj < p) ;
214 ASSERT (Rj [pj] == j) ;
215 #ifdef DO_MAP
216 Map2 [p] = pj ;
217 duplicates = TRUE ;
218 #endif
219 #ifdef DO_VALUES
220 /* sum the entry */
221 #ifdef COMPLEX
222 if (split)
223 {
224 Rx [pj] += Rx [p] ;
225 Rz [pj] += Rz [p] ;
226 }
227 else
228 {
229 Rx[2*pj ] += Rx[2*p ] ;
230 Rx[2*pj+1] += Rx[2*p+1] ;
231 }
232 #else
233 Rx [pj] += Rx [p] ;
234 #endif
235 #endif
236 }
237 else
238 {
239 /* keep the entry */
240 /* also keep track in W[j] of position of a_ij for case above */
241 W [j] = pdest ;
242 #ifdef DO_MAP
243 Map2 [p] = pdest ;
244 #endif
245 /* no need to move the entry if pdest is equal to p */
246 if (pdest != p)
247 {
248 Rj [pdest] = j ;
249 #ifdef DO_VALUES
250 #ifdef COMPLEX
251 if (split)
252 {
253 Rx [pdest] = Rx [p] ;
254 Rz [pdest] = Rz [p] ;
255 }
256 else
257 {
258 Rx [2*pdest ] = Rx [2*p ] ;
259 Rx [2*pdest+1] = Rx [2*p+1] ;
260 }
261 #else
262 Rx [pdest] = Rx [p] ;
263 #endif
264 #endif
265 }
266 pdest++ ;
267 }
268 }
269 RowCount [i] = pdest - p1 ;
270 }
271
272 /* done using W for position of a_ij ] */
273
274 /* ---------------------------------------------------------------------- */
275 /* merge Map and Map2 into a single Map */
276 /* ---------------------------------------------------------------------- */
277
278 #ifdef DO_MAP
279 if (duplicates)
280 {
281 for (k = 0 ; k < nz ; k++)
282 {
283 Map [k] = Map2 [Map [k]] ;
284 }
285 }
286 #ifndef NDEBUG
287 else
288 {
289 /* no duplicates, so no need to recompute Map */
290 for (k = 0 ; k < nz ; k++)
291 {
292 ASSERT (Map2 [k] == k) ;
293 }
294 }
295 for (k = 0 ; k < nz ; k++)
296 {
297 /* make sure that kth triplet is mapped correctly */
298 p = Map [k] ;
299 DEBUG1 (("Second row map: Map ["ID"] = "ID"\n", k, p)) ;
300 i = Ti [k] ;
301 j = Tj [k] ;
302 ASSERT (j == Rj [p]) ;
303 ASSERT (Rp [i] <= p && p < Rp [i+1]) ;
304 }
305 #endif
306 #endif
307
308 /* now the kth triplet maps to p = Map [k], and thus to Rj/Rx [p] */
309
310 /* ---------------------------------------------------------------------- */
311 /* count the entries in each column */
312 /* ---------------------------------------------------------------------- */
313
314 /* [ use W as work space for column counts of A */
315 for (j = 0 ; j < n_col ; j++)
316 {
317 W [j] = 0 ;
318 }
319
320 for (i = 0 ; i < n_row ; i++)
321 {
322 for (p = Rp [i] ; p < Rp [i] + RowCount [i] ; p++)
323 {
324 j = Rj [p] ;
325 ASSERT (j >= 0 && j < n_col) ;
326 W [j]++ ;
327 }
328 }
329
330 /* ---------------------------------------------------------------------- */
331 /* create the column pointers */
332 /* ---------------------------------------------------------------------- */
333
334 Ap [0] = 0 ;
335 for (j = 0 ; j < n_col ; j++)
336 {
337 Ap [j+1] = Ap [j] + W [j] ;
338 }
339 /* done using W as workspace for column counts of A ] */
340
341 for (j = 0 ; j < n_col ; j++)
342 {
343 W [j] = Ap [j] ;
344 }
345
346 /* ---------------------------------------------------------------------- */
347 /* construct the column form */
348 /* ---------------------------------------------------------------------- */
349
350 for (i = 0 ; i < n_row ; i++)
351 {
352 for (p = Rp [i] ; p < Rp [i] + RowCount [i] ; p++)
353 {
354 cp = W [Rj [p]]++ ;
355 #ifdef DO_MAP
356 Map2 [p] = cp ;
357 #endif
358 Ai [cp] = i ;
359 #ifdef DO_VALUES
360 #ifdef COMPLEX
361 if (split)
362 {
363 Ax [cp] = Rx [p] ;
364 Az [cp] = Rz [p] ;
365 }
366 else
367 {
368 Ax [2*cp ] = Rx [2*p ] ;
369 Ax [2*cp+1] = Rx [2*p+1] ;
370 }
371 #else
372 Ax [cp] = Rx [p] ;
373 #endif
374 #endif
375 }
376 }
377
378 /* ---------------------------------------------------------------------- */
379 /* merge Map and Map2 into a single Map */
380 /* ---------------------------------------------------------------------- */
381
382 #ifdef DO_MAP
383 for (k = 0 ; k < nz ; k++)
384 {
385 Map [k] = Map2 [Map [k]] ;
386 }
387 #endif
388
389 /* now the kth triplet maps to p = Map [k], and thus to Ai/Ax [p] */
390
391 #ifndef NDEBUG
392 for (j = 0 ; j < n_col ; j++)
393 {
394 ASSERT (W [j] == Ap [j+1]) ;
395 }
396
397 UMF_dump_col_matrix (
398 #ifdef DO_VALUES
399 Ax,
400 #ifdef COMPLEX
401 Az,
402 #endif
403 #else
404 (double *) NULL,
405 #ifdef COMPLEX
406 (double *) NULL,
407 #endif
408 #endif
409 Ai, Ap, n_row, n_col, nz) ;
410
411 #ifdef DO_MAP
412 for (k = 0 ; k < nz ; k++)
413 {
414 /* make sure that kth triplet is mapped correctly */
415 p = Map [k] ;
416 DEBUG1 (("Col map: Map ["ID"] = "ID"\t", k, p)) ;
417 i = Ti [k] ;
418 j = Tj [k] ;
419 ASSERT (i == Ai [p]) ;
420 DEBUG1 ((" i "ID" j "ID" Ap[j] "ID" p "ID" Ap[j+1] "ID"\n",
421 i, j, Ap [j], p, Ap [j+1])) ;
422 ASSERT (Ap [j] <= p && p < Ap [j+1]) ;
423 }
424 #endif
425 #endif
426
427 return (UMFPACK_OK) ;
428 }
429