Line data Source code
1 : #include "cs.h"
2 : #include "scs.h"
3 : /* -----------------------------------------------------------------
4 : * NB: this is a subset of the routines in the CSPARSE package by
5 : * Tim Davis et. al., for the full package please visit
6 : * http://www.cise.ufl.edu/research/sparse/CSparse/
7 : * ----------------------------------------------------------------- */
8 :
9 : /* wrapper for malloc */
10 : static void *scs_cs_malloc(scs_int n, scs_int size) {
11 3290 : return (scs_malloc(n * size));
12 : }
13 :
14 : /* wrapper for calloc */
15 : static void *scs_cs_calloc(scs_int n, scs_int size) {
16 1645 : return (scs_calloc(n, size));
17 : }
18 :
19 : /* wrapper for free */
20 : static void *scs_cs_free(void *p) {
21 6580 : if (p)
22 5922 : scs_free(p); /* free p if it is not already SCS_NULL */
23 : return (SCS_NULL); /* return SCS_NULL to simplify the use of cs_free */
24 : }
25 :
26 : /* C = compressed-column form of a triplet matrix T */
27 329 : scs_cs *scs_cs_compress(const scs_cs *T) {
28 : scs_int m, n, nz, p, k, *Cp, *Ci, *w, *Ti, *Tj;
29 : scs_float *Cx, *Tx;
30 : scs_cs *C;
31 329 : m = T->m;
32 329 : n = T->n;
33 329 : Ti = T->i;
34 329 : Tj = T->p;
35 329 : Tx = T->x;
36 329 : nz = T->nz;
37 329 : C = scs_cs_spalloc(m, n, nz, Tx != SCS_NULL, 0); /* allocate result */
38 329 : w = scs_cs_calloc(n, sizeof (scs_int)); /* get workspace */
39 329 : if (C == SCS_NULL || w == SCS_NULL)
40 0 : return (scs_cs_done(C, w, SCS_NULL, 0)); /* out of memory */
41 329 : Cp = C->p;
42 329 : Ci = C->i;
43 329 : Cx = C->x;
44 155313 : for (k = 0; k < nz; k++)
45 154984 : w[Tj[k]]++; /* column counts */
46 329 : scs_cs_cumsum(Cp, w, n); /* column pointers */
47 155313 : for (k = 0; k < nz; k++) {
48 154984 : Ci[p = w[Tj[k]]++] = Ti[k]; /* A(i,j) is the pth entry in C */
49 154984 : if (Cx)
50 154984 : Cx[p] = Tx[k];
51 : }
52 329 : return (scs_cs_done(C, w, SCS_NULL, 1)); /* success; free w and return C */
53 : }
54 :
55 658 : scs_cs *scs_cs_done(scs_cs *C, void *w, void *x, scs_int ok) {
56 : scs_cs_free(w); /* free workspace */
57 : scs_cs_free(x);
58 658 : return (ok ? C : scs_cs_spfree(C)); /* return result if OK, else free it */
59 : }
60 :
61 987 : scs_cs *scs_cs_spalloc(scs_int m, scs_int n, scs_int nzmax, scs_int values,
62 : scs_int triplet) {
63 987 : scs_cs *A = scs_cs_calloc(1, sizeof (scs_cs)); /* allocate the cs struct */
64 987 : if (A == SCS_NULL) return (SCS_NULL); /* out of memory */
65 987 : A->m = m; /* define dimensions and nzmax */
66 987 : A->n = n;
67 987 : A->nzmax = nzmax = MAX(nzmax, 1);
68 987 : A->nz = triplet ? 0 : -1; /* allocate triplet or comp.col */
69 1974 : A->p = scs_cs_malloc(triplet ? nzmax : n + 1, sizeof (scs_int));
70 987 : A->i = scs_cs_malloc(nzmax, sizeof (scs_int));
71 1974 : A->x = values ? scs_cs_malloc(nzmax, sizeof (scs_float)) : SCS_NULL;
72 987 : return ((!A->p || !A->i || (values && !A->x)) ? scs_cs_spfree(A) : A);
73 : }
74 :
75 1316 : scs_cs *scs_cs_spfree(scs_cs *A) {
76 1316 : if (A == SCS_NULL) return (SCS_NULL); /* do nothing if A already SCS_NULL */
77 1316 : scs_cs_free(A->p);
78 1316 : scs_cs_free(A->i);
79 1316 : scs_cs_free(A->x);
80 : return ((scs_cs *) scs_cs_free(A)); /* free the cs struct and return SCS_NULL */
81 : }
82 :
83 658 : scs_float scs_cs_cumsum(scs_int *p, scs_int *c, scs_int n) {
84 658 : scs_int i, nz = 0;
85 658 : scs_float nz2 = 0;
86 658 : if (p == SCS_NULL || c == SCS_NULL) return (-1); /* check inputs */
87 127984 : for (i = 0; i < n; i++) {
88 127984 : p[i] = nz;
89 127984 : nz += c[i];
90 127984 : nz2 += c[i]; /* also in scs_float to avoid scs_int overflow */
91 127984 : c[i] = p[i]; /* also copy p[0..n-1] back into c[0..n-1]*/
92 : }
93 658 : p[n] = nz;
94 658 : return (nz2); /* return sum (c [0..n-1]) */
95 : }
96 :
97 329 : scs_int *scs_cs_pinv(scs_int const *p, scs_int n) {
98 : scs_int k, *pinv;
99 329 : if (p == SCS_NULL) return (SCS_NULL); /* p = SCS_NULL denotes identity */
100 329 : pinv = scs_cs_malloc(n, sizeof (scs_int)); /* allocate result */
101 329 : if (pinv == SCS_NULL)
102 : return (SCS_NULL); /* out of memory */
103 63992 : for (k = 0; k < n; k++)
104 63992 : pinv[p[k]] = k; /* invert the permutation */
105 : return (pinv); /* return result */
106 : }
107 :
108 329 : scs_cs *scs_cs_symperm(const scs_cs *A, const scs_int *pinv, scs_int values) {
109 : scs_int i, j, p, q, i2, j2, n, *Ap, *Ai, *Cp, *Ci, *w;
110 : scs_float *Cx, *Ax;
111 : scs_cs *C;
112 329 : n = A->n;
113 329 : Ap = A->p;
114 329 : Ai = A->i;
115 329 : Ax = A->x;
116 329 : C = scs_cs_spalloc(n, n, Ap[n], values && (Ax != SCS_NULL),
117 : 0); /* alloc result*/
118 329 : w = scs_cs_calloc(n, sizeof (scs_int)); /* get workspace */
119 329 : if (C == SCS_NULL || w == SCS_NULL)
120 0 : return (scs_cs_done(C, w, SCS_NULL, 0)); /* out of memory */
121 329 : Cp = C->p;
122 329 : Ci = C->i;
123 329 : Cx = C->x;
124 64321 : for (j = 0; j < n; j++) /* count entries in each column of C */ {
125 63992 : j2 = pinv ? pinv[j] : j; /* column j of A is column j2 of C */
126 218976 : for (p = Ap[j]; p < Ap[j + 1]; p++) {
127 154984 : i = Ai[p];
128 154984 : if (i > j)
129 0 : continue; /* skip lower triangular part of A */
130 154984 : i2 = pinv ? pinv[i] : i; /* row i of A is row i2 of C */
131 154984 : w[MAX(i2, j2)]++; /* column count of C */
132 : }
133 : }
134 329 : scs_cs_cumsum(Cp, w, n); /* compute column pointers of C */
135 64321 : for (j = 0; j < n; j++) {
136 63992 : j2 = pinv ? pinv[j] : j; /* column j of A is column j2 of C */
137 218976 : for (p = Ap[j]; p < Ap[j + 1]; p++) {
138 154984 : i = Ai[p];
139 154984 : if (i > j)
140 0 : continue; /* skip lower triangular part of A*/
141 154984 : i2 = pinv ? pinv[i] : i; /* row i of A is row i2 of C */
142 154984 : Ci[q = w[MAX(i2, j2)]++] = MIN(i2, j2);
143 154984 : if (Cx)
144 154984 : Cx[q] = Ax[p];
145 : }
146 : }
147 329 : return (scs_cs_done(C, w, SCS_NULL, 1)); /* success; free workspace, return C */
148 : }
|
