SCM

SCM Repository

[matrix] Annotation of /pkg/src/triplet_to_col.c
ViewVC logotype

Annotation of /pkg/src/triplet_to_col.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 47 - (view) (download) (as text)

1 : bates 47 /* Based on UMF_triplet.c from UMFPACK which carries the following copyright */
2 :    
3 :     /* -------------------------------------------------------------------------- */
4 :     /* UMFPACK Version 4.1 (Apr. 30, 2003), Copyright (c) 2003 by Timothy A. */
5 :     /* Davis. All Rights Reserved. See ../README for License. */
6 :     /* email: davis@cise.ufl.edu CISE Department, Univ. of Florida. */
7 :     /* web: http://www.cise.ufl.edu/research/sparse/umfpack */
8 :     /* -------------------------------------------------------------------------- */
9 :    
10 :     #include <R_ext/RS.h>
11 :    
12 :     void triplet_to_col
13 :     (
14 :     int n_row,
15 :     int n_col,
16 :     int nz,
17 :     const int Ti [ ], /* size nz */
18 :     const int Tj [ ], /* size nz */
19 :     const double Tx [ ], /* size nz */
20 :     int Ap [ ], /* size n_col + 1 */
21 :     int Ai [ ], /* size nz */
22 :     double Ax [ ] /* size nz */
23 :     )
24 :     {
25 :     int i, j, k, p, cp, p1, p2, pdest, pj;
26 :     int *Rp = Calloc((n_row + 1), int),
27 :     *Rj = Calloc(nz, int),
28 :     *W = Calloc((n_row > n_col) ? n_row : n_col, int),
29 :     *RowCount = Calloc(n_row, int);
30 :     double *Rx = Calloc(nz, double);
31 :    
32 :     /* count the entries in each row (including duplicates) */
33 :     /* use W as workspace for row counts (including duplicates) */
34 :     memset(W, 0, sizeof(int) * n_row);
35 :     for (k = 0 ; k < nz ; k++) {
36 :     i = Ti [k] ;
37 :     j = Tj [k] ;
38 :     if (i < 0 || i >= n_row || j < 0 || j >= n_col)
39 :     error("entry %d in input has row %d and column %d", k, i, j);
40 :     W [i]++ ;
41 :     }
42 :     /* compute the row pointers */
43 :     Rp [0] = 0 ;
44 :     for (i = 0 ; i < n_row ; i++) {
45 :     Rp [i+1] = Rp [i] + W [i] ;
46 :     W [i] = Rp [i] ;
47 :     }
48 :     /* W is now equal to the row pointers */
49 :    
50 :     /* ---------------------------------------------------------------------- */
51 :     /* construct the row form */
52 :     /* ---------------------------------------------------------------------- */
53 :     for (k = 0 ; k < nz ; k++) {
54 :     p = W [Ti [k]]++ ;
55 :     Rj [p] = Tj [k] ;
56 :     Rx [p] = Tx [k] ;
57 :     }
58 :     /* Rp stays the same, but W [i] is advanced to the start of row i+1 */
59 :     /* ---------------------------------------------------------------------- */
60 :     /* sum up duplicates */
61 :     /* ---------------------------------------------------------------------- */
62 :     /* use W [j] to hold position in Ri/Rx/Rz of a_ij, for row i [ */
63 :     for (j = 0 ; j < n_col ; j++) {
64 :     W [j] = -1;
65 :     }
66 :     for (i = 0 ; i < n_row ; i++) {
67 :     p1 = Rp [i] ;
68 :     p2 = Rp [i+1] ;
69 :     pdest = p1 ;
70 :     /* At this point, W [j] < p1 holds true for all columns j, */
71 :     /* because Ri/Rx/Rz is stored in row oriented order. */
72 :     for (p = p1; p < p2; p++) {
73 :     j = Rj [p] ;
74 :     pj = W [j] ;
75 :     if (pj >= p1) {
76 :     /* this column index, j, is already in row i, at position pj */
77 :     /* sum the entry */
78 :     Rx [pj] += Rx [p] ;
79 :     } else {
80 :     /* keep the entry */
81 :     /* also keep track in W[j] of position of a_ij for case above */
82 :     W [j] = pdest ;
83 :     /* no need to move the entry if pdest is equal to p */
84 :     if (pdest != p)
85 :     {
86 :     Rj [pdest] = j ;
87 :     Rx [pdest] = Rx [p] ;
88 :     }
89 :     pdest++ ;
90 :     }
91 :     }
92 :     RowCount [i] = pdest - p1 ;
93 :     }
94 :     /* done using W for position of a_ij ] */
95 :     /* ---------------------------------------------------------------------- */
96 :     /* count the entries in each column */
97 :     /* ---------------------------------------------------------------------- */
98 :     /* [ use W as work space for column counts of A */
99 :     memset(W, 0, sizeof(int) * n_col);
100 :     for (i = 0 ; i < n_row ; i++) {
101 :     for (p = Rp [i] ; p < Rp [i] + RowCount [i] ; p++) {
102 :     j = Rj [p] ;
103 :     W [j]++ ;
104 :     }
105 :     }
106 :     /* ---------------------------------------------------------------------- */
107 :     /* create the column pointers */
108 :     /* ---------------------------------------------------------------------- */
109 :     Ap [0] = 0 ;
110 :     for (j = 0 ; j < n_col ; j++) {
111 :     Ap [j+1] = Ap [j] + W [j] ;
112 :     }
113 :     /* done using W as workspace for column counts of A ] */
114 :     for (j = 0 ; j < n_col ; j++) {
115 :     W [j] = Ap [j] ;
116 :     }
117 :     /* ---------------------------------------------------------------------- */
118 :     /* construct the column form */
119 :     /* ---------------------------------------------------------------------- */
120 :    
121 :     for (i = 0 ; i < n_row ; i++) {
122 :     for (p = Rp [i] ; p < Rp [i] + RowCount [i] ; p++) {
123 :     cp = W [Rj [p]]++ ;
124 :     Ai [cp] = i ;
125 :     Ax [cp] = Rx [p] ;
126 :     }
127 :     }
128 :     Free(Rp); Free(Rj); Free(W); Free(RowCount); Free(Rx);
129 :     }

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge