1 /* ========================================================================== */
2 /* === MatrixOps/cholmod_drop =============================================== */
3 /* ========================================================================== */
4
5 /* -----------------------------------------------------------------------------
6 * CHOLMOD/MatrixOps Module. Copyright (C) 2005-2006, Timothy A. Davis
7 * http://www.suitesparse.com
8 * -------------------------------------------------------------------------- */
9
10 /* Drop small entries from A, and entries in the ignored part of A if A
11 * is symmetric. None of the matrix operations drop small numerical entries
12 * from a matrix, except for this one. NaN's and Inf's are kept.
13 *
14 * workspace: none
15 *
16 * Supports pattern and real matrices, complex and zomplex not supported.
17 */
18
19 #ifndef NGPL
20 #ifndef NMATRIXOPS
21
22 #include "cholmod_internal.h"
23 #include "cholmod_matrixops.h"
24
25
26 /* ========================================================================== */
27 /* === cholmod_drop ========================================================= */
28 /* ========================================================================== */
29
CHOLMOD(drop)30 int CHOLMOD(drop)
31 (
32 /* ---- input ---- */
33 double tol, /* keep entries with absolute value > tol */
34 /* ---- in/out --- */
35 cholmod_sparse *A, /* matrix to drop entries from */
36 /* --------------- */
37 cholmod_common *Common
38 )
39 {
40 double aij ;
41 double *Ax ;
42 Int *Ap, *Ai, *Anz ;
43 Int packed, i, j, nrow, ncol, p, pend, nz, values ;
44
45 /* ---------------------------------------------------------------------- */
46 /* check inputs */
47 /* ---------------------------------------------------------------------- */
48
49 RETURN_IF_NULL_COMMON (FALSE) ;
50 RETURN_IF_NULL (A, FALSE) ;
51 RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_REAL, FALSE) ;
52 Common->status = CHOLMOD_OK ;
53 ASSERT (CHOLMOD(dump_sparse) (A, "A predrop", Common) >= 0) ;
54
55 /* ---------------------------------------------------------------------- */
56 /* get inputs */
57 /* ---------------------------------------------------------------------- */
58
59 Ap = A->p ;
60 Ai = A->i ;
61 Ax = A->x ;
62 Anz = A->nz ;
63 packed = A->packed ;
64 ncol = A->ncol ;
65 nrow = A->nrow ;
66 values = (A->xtype != CHOLMOD_PATTERN) ;
67 nz = 0 ;
68
69 if (values)
70 {
71
72 /* ------------------------------------------------------------------ */
73 /* drop small numerical entries from A, and entries in ignored part */
74 /* ------------------------------------------------------------------ */
75
76 if (A->stype > 0)
77 {
78
79 /* -------------------------------------------------------------- */
80 /* A is symmetric, with just upper triangular part stored */
81 /* -------------------------------------------------------------- */
82
83 for (j = 0 ; j < ncol ; j++)
84 {
85 p = Ap [j] ;
86 pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
87 Ap [j] = nz ;
88 for ( ; p < pend ; p++)
89 {
90 i = Ai [p] ;
91 aij = Ax [p] ;
92 if (i <= j && (fabs (aij) > tol || IS_NAN (aij)))
93 {
94 Ai [nz] = i ;
95 Ax [nz] = aij ;
96 nz++ ;
97 }
98 }
99 }
100
101 }
102 else if (A->stype < 0)
103 {
104
105 /* -------------------------------------------------------------- */
106 /* A is symmetric, with just lower triangular part stored */
107 /* -------------------------------------------------------------- */
108
109 for (j = 0 ; j < ncol ; j++)
110 {
111 p = Ap [j] ;
112 pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
113 Ap [j] = nz ;
114 for ( ; p < pend ; p++)
115 {
116 i = Ai [p] ;
117 aij = Ax [p] ;
118 if (i >= j && (fabs (aij) > tol || IS_NAN (aij)))
119 {
120 Ai [nz] = i ;
121 Ax [nz] = aij ;
122 nz++ ;
123 }
124 }
125 }
126 }
127 else
128 {
129
130 /* -------------------------------------------------------------- */
131 /* both parts of A present, just drop small entries */
132 /* -------------------------------------------------------------- */
133
134 for (j = 0 ; j < ncol ; j++)
135 {
136 p = Ap [j] ;
137 pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
138 Ap [j] = nz ;
139 for ( ; p < pend ; p++)
140 {
141 i = Ai [p] ;
142 aij = Ax [p] ;
143 if (fabs (aij) > tol || IS_NAN (aij))
144 {
145 Ai [nz] = i ;
146 Ax [nz] = aij ;
147 nz++ ;
148 }
149 }
150 }
151 }
152 Ap [ncol] = nz ;
153
154 /* reduce A->i and A->x in size */
155 ASSERT (MAX (1,nz) <= A->nzmax) ;
156 CHOLMOD(reallocate_sparse) (nz, A, Common) ;
157 ASSERT (Common->status >= CHOLMOD_OK) ;
158
159 }
160 else
161 {
162
163 /* ------------------------------------------------------------------ */
164 /* consider only the pattern of A */
165 /* ------------------------------------------------------------------ */
166
167 /* Note that cholmod_band_inplace calls cholmod_reallocate_sparse */
168 if (A->stype > 0)
169 {
170 CHOLMOD(band_inplace) (0, ncol, 0, A, Common) ;
171 }
172 else if (A->stype < 0)
173 {
174 CHOLMOD(band_inplace) (-nrow, 0, 0, A, Common) ;
175 }
176 }
177
178 ASSERT (CHOLMOD(dump_sparse) (A, "A dropped", Common) >= 0) ;
179 return (TRUE) ;
180 }
181 #endif
182 #endif
183