← Back to branch summary

~ubuntu-branches/ubuntu/raring/python-scipy/raring-proposed

« back to all changes in this revision

Viewing changes to Lib/sandbox/pysparse/umfpack/umf_triplet.c

  • Committer: Bazaar Package Importer
  • Author(s): Matthias Klose
  • Date: 2007-01-07 14:12:12 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20070107141212-mm0ebkh5b37hcpzn
http://bugs.debian.org/402783
* Remove build dependency on python-numpy-dev.
* python-scipy: Depend on python-numpy instead of python-numpy-dev.
* Package builds on other archs than i386. Closes: #402783.

Show diffs side-by-side

added added

removed removed

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