1/* $Id: CoinPresolveZeros.cpp 1373 2011-01-03 23:57:44Z lou $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others. All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#include <stdio.h>
7#include <math.h>
8
9#include "CoinHelperFunctions.hpp"
10#include "CoinPresolveMatrix.hpp"
11#include "CoinPresolveZeros.hpp"
12
13#if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
14#include "CoinPresolvePsdebug.hpp"
15#endif
16
17namespace { // begin unnamed file-local namespace
18
19// Sees how many zeros there are
20static int count_col_zeros (int ncheckcols, int * checkcols,
21 const CoinBigIndex *mcstrt, double *colels,// int *hrow,
22 int *hincol)
23{
24 int nactions = 0;
25 int i;
26
27 for (i=0; i<ncheckcols; i++) {
28 int col = checkcols[i];
29 CoinBigIndex kcs = mcstrt[col];
30 CoinBigIndex kce = mcstrt[col] + hincol[col];
31 CoinBigIndex k;
32
33 for (k=kcs; k<kce; ++k) {
34 if (fabs(colels[k]) < ZTOLDP) {
35 nactions++;
36 }
37 }
38 }
39 return (nactions);
40}
41
42// Sees how many zeros there are
43static int count_col_zeros2 (int ncheckcols, int * checkcols,
44 const CoinBigIndex *mcstrt, double *colels,// int *hrow,
45 int *hincol)
46{
47 int nactions = 0;
48
49 for (int col=0; col<ncheckcols; col++) {
50 CoinBigIndex kcs = mcstrt[col];
51 CoinBigIndex kce = mcstrt[col] + hincol[col];
52 CoinBigIndex k;
53
54 for (k=kcs; k<kce; ++k) {
55 if (fabs(colels[k]) < ZTOLDP) {
56 checkcols[nactions++]=col;
57 }
58 }
59 }
60 return (nactions);
61}
62
63// searches the cols in checkcols for zero entries.
64// creates a dropped_zero entry for each one; doesn't check for out-of-memory.
65// returns number of zeros found.
66
67static int drop_col_zeros (int ncheckcols, int *checkcols,
68 const CoinBigIndex *mcstrt, double *colels, int *hrow,
69 int *hincol, presolvehlink *clink,
70 dropped_zero *actions)
71{
72 typedef dropped_zero action;
73 int nactions = 0;
74 int i;
75
76 for (i=0; i<ncheckcols; i++) {
77 int col = checkcols[i];
78 CoinBigIndex kcs = mcstrt[col];
79 CoinBigIndex kce = mcstrt[col] + hincol[col];
80 CoinBigIndex k;
81
82 for (k=kcs; k<kce; ++k) {
83 if (fabs(colels[k]) < ZTOLDP) {
84 actions[nactions].col = col;
85 actions[nactions].row = hrow[k];
86
87# if PRESOLVE_DEBUG
88 if (nactions == 0)
89 printf("ZEROS: ");
90 else
91 if (nactions%10 == 0)
92 printf("\n") ;
93 printf("(%d,%d) ", hrow[k], col);
94# endif
95
96 nactions++;
97
98 colels[k] = colels[kce-1];
99 hrow[k] = hrow[kce-1];
100 kce--;
101 hincol[col]--;
102
103 --k; // redo this position
104 }
105 }
106 if (hincol[col] == 0)
107 PRESOLVE_REMOVE_LINK(clink,col) ;
108 }
109
110# if PRESOLVE_DEBUG
111 if (nactions)
112 printf("\n");
113# endif
114
115 return (nactions);
116}
117
118// very similar to col, but without the buffer and reads zeros
119
120static void drop_row_zeros(int nzeros, const dropped_zero *zeros,
121 const CoinBigIndex *mrstrt, double *rowels, int *hcol,
122 int *hinrow, presolvehlink *rlink)
123{
124 int i;
125 for (i=0; i<nzeros; i++) {
126 int row = zeros[i].row;
127 CoinBigIndex krs = mrstrt[row];
128 CoinBigIndex kre = mrstrt[row] + hinrow[row];
129 CoinBigIndex k;
130
131 for (k=krs; k<kre; k++) {
132 if (fabs(rowels[k]) < ZTOLDP) {
133 rowels[k] = rowels[kre-1];
134 hcol[k] = hcol[kre-1];
135 kre--;
136 hinrow[row]--;
137
138 --k; // redo this position
139 }
140 }
141 if (hinrow[row] == 0)
142 PRESOLVE_REMOVE_LINK(rlink,row) ;
143 }
144}
145
146} // end unnamed file-local namespace
147
148const CoinPresolveAction
149 *drop_zero_coefficients_action::presolve (CoinPresolveMatrix *prob,
150 int *checkcols,
151 int ncheckcols,
152 const CoinPresolveAction *next)
153{
154 double *colels = prob->colels_;
155 int *hrow = prob->hrow_;
156 CoinBigIndex *mcstrt = prob->mcstrt_;
157 int *hincol = prob->hincol_;
158 presolvehlink *clink = prob->clink_ ;
159 presolvehlink *rlink = prob->rlink_ ;
160
161 // int i;
162 int nzeros;
163 if (ncheckcols==prob->ncols_) {
164 // can do faster
165 nzeros = count_col_zeros2(ncheckcols,checkcols,
166 mcstrt,colels,/*hrow,*/hincol);
167 } else {
168 nzeros = count_col_zeros(ncheckcols,checkcols,
169 mcstrt,colels,/*hrow,*/hincol);
170 }
171 if (nzeros == 0) {
172 return (next);
173 } else {
174 dropped_zero * zeros = new dropped_zero[nzeros];
175
176 nzeros=drop_col_zeros((ncheckcols==prob->ncols_) ? nzeros : ncheckcols,
177 checkcols,
178 mcstrt,colels,hrow,hincol,clink,
179 zeros);
180 double *rowels = prob->rowels_;
181 int *hcol = prob->hcol_;
182 CoinBigIndex *mrstrt = prob->mrstrt_;
183 int *hinrow = prob->hinrow_;
184 // int nrows = prob->nrows_;
185
186# if PRESOLVE_SUMMARY
187 printf("NZEROS: %d\n", nzeros);
188# endif
189
190 // make the row rep consistent
191 drop_row_zeros(nzeros,zeros,mrstrt,rowels,hcol,hinrow,rlink) ;
192
193 dropped_zero *zeros1 = new dropped_zero[nzeros];
194 CoinMemcpyN(zeros, nzeros, zeros1);
195
196 delete [] zeros;
197 return (new drop_zero_coefficients_action(nzeros, zeros1, next));
198 }
199}
200
201
202const CoinPresolveAction *drop_zero_coefficients(CoinPresolveMatrix *prob,
203 const CoinPresolveAction *next)
204{
205 int ncheck = prob->ncols_;
206 int *checkcols = new int[ncheck];
207
208 if (prob->anyProhibited()) {
209 // some prohibited
210 ncheck=0;
211 for (int i=0; i<prob->ncols_; i++)
212 if (!prob->colProhibited(i))
213 checkcols[ncheck++] = i;
214 }
215
216 const CoinPresolveAction *retval = drop_zero_coefficients_action::presolve(prob,
217 checkcols, ncheck,
218 next);
219 delete[]checkcols;
220 return (retval);
221}
222
223void drop_zero_coefficients_action::postsolve(CoinPostsolveMatrix *prob) const
224{
225 const int nzeros = nzeros_;
226 const dropped_zero *const zeros = zeros_;
227
228 double *colels = prob->colels_;
229 int *hrow = prob->hrow_;
230 CoinBigIndex *mcstrt = prob->mcstrt_;
231 int *hincol = prob->hincol_;
232 int *link = prob->link_;
233 CoinBigIndex &free_list = prob->free_list_;
234
235 for (const dropped_zero *z = &zeros[nzeros-1]; zeros<=z; z--) {
236 int irow = z->row;
237 int jcol = z->col;
238
239 {
240 CoinBigIndex k = free_list;
241 assert(k >= 0 && k < prob->bulk0_) ;
242 free_list = link[free_list];
243 hrow[k] = irow;
244 colels[k] = 0.0;
245 link[k] = mcstrt[jcol];
246 mcstrt[jcol] = k;
247 }
248
249 hincol[jcol]++;
250 }
251
252# if PRESOLVE_CONSISTENCY
253 presolve_check_free_list(prob) ;
254# endif
255
256}
257