COIN-OR::LEMON - Graph Library

source: glpk-cmake/src/glpscl.c @ 1:c445c931472f

Last change on this file since 1:c445c931472f was 1:c445c931472f, checked in by Alpar Juttner <alpar@…>, 10 years ago

Import glpk-4.45

  • Generated files and doc/notes are removed
File size: 15.7 KB
Line 
1/* glpscl.c */
2
3/***********************************************************************
4*  This code is part of GLPK (GNU Linear Programming Kit).
5*
6*  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
7*  2009, 2010 Andrew Makhorin, Department for Applied Informatics,
8*  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
9*  E-mail: <mao@gnu.org>.
10*
11*  GLPK is free software: you can redistribute it and/or modify it
12*  under the terms of the GNU General Public License as published by
13*  the Free Software Foundation, either version 3 of the License, or
14*  (at your option) any later version.
15*
16*  GLPK is distributed in the hope that it will be useful, but WITHOUT
17*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
19*  License for more details.
20*
21*  You should have received a copy of the GNU General Public License
22*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
23***********************************************************************/
24
25#include "glpapi.h"
26
27/***********************************************************************
28*  min_row_aij - determine minimal |a[i,j]| in i-th row
29*
30*  This routine returns minimal magnitude of (non-zero) constraint
31*  coefficients in i-th row of the constraint matrix.
32*
33*  If the parameter scaled is zero, the original constraint matrix A is
34*  assumed. Otherwise, the scaled constraint matrix R*A*S is assumed.
35*
36*  If i-th row of the matrix is empty, the routine returns 1. */
37
38static double min_row_aij(glp_prob *lp, int i, int scaled)
39{     GLPAIJ *aij;
40      double min_aij, temp;
41      xassert(1 <= i && i <= lp->m);
42      min_aij = 1.0;
43      for (aij = lp->row[i]->ptr; aij != NULL; aij = aij->r_next)
44      {  temp = fabs(aij->val);
45         if (scaled) temp *= (aij->row->rii * aij->col->sjj);
46         if (aij->r_prev == NULL || min_aij > temp)
47            min_aij = temp;
48      }
49      return min_aij;
50}
51
52/***********************************************************************
53*  max_row_aij - determine maximal |a[i,j]| in i-th row
54*
55*  This routine returns maximal magnitude of (non-zero) constraint
56*  coefficients in i-th row of the constraint matrix.
57*
58*  If the parameter scaled is zero, the original constraint matrix A is
59*  assumed. Otherwise, the scaled constraint matrix R*A*S is assumed.
60*
61*  If i-th row of the matrix is empty, the routine returns 1. */
62
63static double max_row_aij(glp_prob *lp, int i, int scaled)
64{     GLPAIJ *aij;
65      double max_aij, temp;
66      xassert(1 <= i && i <= lp->m);
67      max_aij = 1.0;
68      for (aij = lp->row[i]->ptr; aij != NULL; aij = aij->r_next)
69      {  temp = fabs(aij->val);
70         if (scaled) temp *= (aij->row->rii * aij->col->sjj);
71         if (aij->r_prev == NULL || max_aij < temp)
72            max_aij = temp;
73      }
74      return max_aij;
75}
76
77/***********************************************************************
78*  min_col_aij - determine minimal |a[i,j]| in j-th column
79*
80*  This routine returns minimal magnitude of (non-zero) constraint
81*  coefficients in j-th column of the constraint matrix.
82*
83*  If the parameter scaled is zero, the original constraint matrix A is
84*  assumed. Otherwise, the scaled constraint matrix R*A*S is assumed.
85*
86*  If j-th column of the matrix is empty, the routine returns 1. */
87
88static double min_col_aij(glp_prob *lp, int j, int scaled)
89{     GLPAIJ *aij;
90      double min_aij, temp;
91      xassert(1 <= j && j <= lp->n);
92      min_aij = 1.0;
93      for (aij = lp->col[j]->ptr; aij != NULL; aij = aij->c_next)
94      {  temp = fabs(aij->val);
95         if (scaled) temp *= (aij->row->rii * aij->col->sjj);
96         if (aij->c_prev == NULL || min_aij > temp)
97            min_aij = temp;
98      }
99      return min_aij;
100}
101
102/***********************************************************************
103*  max_col_aij - determine maximal |a[i,j]| in j-th column
104*
105*  This routine returns maximal magnitude of (non-zero) constraint
106*  coefficients in j-th column of the constraint matrix.
107*
108*  If the parameter scaled is zero, the original constraint matrix A is
109*  assumed. Otherwise, the scaled constraint matrix R*A*S is assumed.
110*
111*  If j-th column of the matrix is empty, the routine returns 1. */
112
113static double max_col_aij(glp_prob *lp, int j, int scaled)
114{     GLPAIJ *aij;
115      double max_aij, temp;
116      xassert(1 <= j && j <= lp->n);
117      max_aij = 1.0;
118      for (aij = lp->col[j]->ptr; aij != NULL; aij = aij->c_next)
119      {  temp = fabs(aij->val);
120         if (scaled) temp *= (aij->row->rii * aij->col->sjj);
121         if (aij->c_prev == NULL || max_aij < temp)
122            max_aij = temp;
123      }
124      return max_aij;
125}
126
127/***********************************************************************
128*  min_mat_aij - determine minimal |a[i,j]| in constraint matrix
129*
130*  This routine returns minimal magnitude of (non-zero) constraint
131*  coefficients in the constraint matrix.
132*
133*  If the parameter scaled is zero, the original constraint matrix A is
134*  assumed. Otherwise, the scaled constraint matrix R*A*S is assumed.
135*
136*  If the matrix is empty, the routine returns 1. */
137
138static double min_mat_aij(glp_prob *lp, int scaled)
139{     int i;
140      double min_aij, temp;
141      min_aij = 1.0;
142      for (i = 1; i <= lp->m; i++)
143      {  temp = min_row_aij(lp, i, scaled);
144         if (i == 1 || min_aij > temp)
145            min_aij = temp;
146      }
147      return min_aij;
148}
149
150/***********************************************************************
151*  max_mat_aij - determine maximal |a[i,j]| in constraint matrix
152*
153*  This routine returns maximal magnitude of (non-zero) constraint
154*  coefficients in the constraint matrix.
155*
156*  If the parameter scaled is zero, the original constraint matrix A is
157*  assumed. Otherwise, the scaled constraint matrix R*A*S is assumed.
158*
159*  If the matrix is empty, the routine returns 1. */
160
161static double max_mat_aij(glp_prob *lp, int scaled)
162{     int i;
163      double max_aij, temp;
164      max_aij = 1.0;
165      for (i = 1; i <= lp->m; i++)
166      {  temp = max_row_aij(lp, i, scaled);
167         if (i == 1 || max_aij < temp)
168            max_aij = temp;
169      }
170      return max_aij;
171}
172
173/***********************************************************************
174*  eq_scaling - perform equilibration scaling
175*
176*  This routine performs equilibration scaling of rows and columns of
177*  the constraint matrix.
178*
179*  If the parameter flag is zero, the routine scales rows at first and
180*  then columns. Otherwise, the routine scales columns and then rows.
181*
182*  Rows are scaled as follows:
183*
184*                         n
185*     a'[i,j] = a[i,j] / max |a[i,j]|,  i = 1,...,m.
186*                        j=1
187*
188*  This makes the infinity (maximum) norm of each row of the matrix
189*  equal to 1.
190*
191*  Columns are scaled as follows:
192*
193*                         n
194*     a'[i,j] = a[i,j] / max |a[i,j]|,  j = 1,...,n.
195*                        i=1
196*
197*  This makes the infinity (maximum) norm of each column of the matrix
198*  equal to 1. */
199
200static void eq_scaling(glp_prob *lp, int flag)
201{     int i, j, pass;
202      double temp;
203      xassert(flag == 0 || flag == 1);
204      for (pass = 0; pass <= 1; pass++)
205      {  if (pass == flag)
206         {  /* scale rows */
207            for (i = 1; i <= lp->m; i++)
208            {  temp = max_row_aij(lp, i, 1);
209               glp_set_rii(lp, i, glp_get_rii(lp, i) / temp);
210            }
211         }
212         else
213         {  /* scale columns */
214            for (j = 1; j <= lp->n; j++)
215            {  temp = max_col_aij(lp, j, 1);
216               glp_set_sjj(lp, j, glp_get_sjj(lp, j) / temp);
217            }
218         }
219      }
220      return;
221}
222
223/***********************************************************************
224*  gm_scaling - perform geometric mean scaling
225*
226*  This routine performs geometric mean scaling of rows and columns of
227*  the constraint matrix.
228*
229*  If the parameter flag is zero, the routine scales rows at first and
230*  then columns. Otherwise, the routine scales columns and then rows.
231*
232*  Rows are scaled as follows:
233*
234*     a'[i,j] = a[i,j] / sqrt(alfa[i] * beta[i]),  i = 1,...,m,
235*
236*  where:
237*                n                        n
238*     alfa[i] = min |a[i,j]|,  beta[i] = max |a[i,j]|.
239*               j=1                      j=1
240*
241*  This allows decreasing the ratio beta[i] / alfa[i] for each row of
242*  the matrix.
243*
244*  Columns are scaled as follows:
245*
246*     a'[i,j] = a[i,j] / sqrt(alfa[j] * beta[j]),  j = 1,...,n,
247*
248*  where:
249*                m                        m
250*     alfa[j] = min |a[i,j]|,  beta[j] = max |a[i,j]|.
251*               i=1                      i=1
252*
253*  This allows decreasing the ratio beta[j] / alfa[j] for each column
254*  of the matrix. */
255
256static void gm_scaling(glp_prob *lp, int flag)
257{     int i, j, pass;
258      double temp;
259      xassert(flag == 0 || flag == 1);
260      for (pass = 0; pass <= 1; pass++)
261      {  if (pass == flag)
262         {  /* scale rows */
263            for (i = 1; i <= lp->m; i++)
264            {  temp = min_row_aij(lp, i, 1) * max_row_aij(lp, i, 1);
265               glp_set_rii(lp, i, glp_get_rii(lp, i) / sqrt(temp));
266            }
267         }
268         else
269         {  /* scale columns */
270            for (j = 1; j <= lp->n; j++)
271            {  temp = min_col_aij(lp, j, 1) * max_col_aij(lp, j, 1);
272               glp_set_sjj(lp, j, glp_get_sjj(lp, j) / sqrt(temp));
273            }
274         }
275      }
276      return;
277}
278
279/***********************************************************************
280*  max_row_ratio - determine worst scaling "quality" for rows
281*
282*  This routine returns the worst scaling "quality" for rows of the
283*  currently scaled constraint matrix:
284*
285*              m
286*     ratio = max ratio[i],
287*             i=1
288*  where:
289*                 n              n
290*     ratio[i] = max |a[i,j]| / min |a[i,j]|,  1 <= i <= m,
291*                j=1            j=1
292*
293*  is the scaling "quality" of i-th row. */
294
295static double max_row_ratio(glp_prob *lp)
296{     int i;
297      double ratio, temp;
298      ratio = 1.0;
299      for (i = 1; i <= lp->m; i++)
300      {  temp = max_row_aij(lp, i, 1) / min_row_aij(lp, i, 1);
301         if (i == 1 || ratio < temp) ratio = temp;
302      }
303      return ratio;
304}
305
306/***********************************************************************
307*  max_col_ratio - determine worst scaling "quality" for columns
308*
309*  This routine returns the worst scaling "quality" for columns of the
310*  currently scaled constraint matrix:
311*
312*              n
313*     ratio = max ratio[j],
314*             j=1
315*  where:
316*                 m              m
317*     ratio[j] = max |a[i,j]| / min |a[i,j]|,  1 <= j <= n,
318*                i=1            i=1
319*
320*  is the scaling "quality" of j-th column. */
321
322static double max_col_ratio(glp_prob *lp)
323{     int j;
324      double ratio, temp;
325      ratio = 1.0;
326      for (j = 1; j <= lp->n; j++)
327      {  temp = max_col_aij(lp, j, 1) / min_col_aij(lp, j, 1);
328         if (j == 1 || ratio < temp) ratio = temp;
329      }
330      return ratio;
331}
332
333/***********************************************************************
334*  gm_iterate - perform iterative geometric mean scaling
335*
336*  This routine performs iterative geometric mean scaling of rows and
337*  columns of the constraint matrix.
338*
339*  The parameter it_max specifies the maximal number of iterations.
340*  Recommended value of it_max is 15.
341*
342*  The parameter tau specifies a minimal improvement of the scaling
343*  "quality" on each iteration, 0 < tau < 1. It means than the scaling
344*  process continues while the following condition is satisfied:
345*
346*     ratio[k] <= tau * ratio[k-1],
347*
348*  where ratio = max |a[i,j]| / min |a[i,j]| is the scaling "quality"
349*  to be minimized, k is the iteration number. Recommended value of tau
350*  is 0.90. */
351
352static void gm_iterate(glp_prob *lp, int it_max, double tau)
353{     int k, flag;
354      double ratio = 0.0, r_old;
355      /* if the scaling "quality" for rows is better than for columns,
356         the rows are scaled first; otherwise, the columns are scaled
357         first */
358      flag = (max_row_ratio(lp) > max_col_ratio(lp));
359      for (k = 1; k <= it_max; k++)
360      {  /* save the scaling "quality" from previous iteration */
361         r_old = ratio;
362         /* determine the current scaling "quality" */
363         ratio = max_mat_aij(lp, 1) / min_mat_aij(lp, 1);
364#if 0
365         xprintf("k = %d; ratio = %g\n", k, ratio);
366#endif
367         /* if improvement is not enough, terminate scaling */
368         if (k > 1 && ratio > tau * r_old) break;
369         /* otherwise, perform another iteration */
370         gm_scaling(lp, flag);
371      }
372      return;
373}
374
375/***********************************************************************
376*  NAME
377*
378*  scale_prob - scale problem data
379*
380*  SYNOPSIS
381*
382*  #include "glpscl.h"
383*  void scale_prob(glp_prob *lp, int flags);
384*
385*  DESCRIPTION
386*
387*  The routine scale_prob performs automatic scaling of problem data
388*  for the specified problem object. */
389
390static void scale_prob(glp_prob *lp, int flags)
391{     static const char *fmt =
392         "%s: min|aij| = %10.3e  max|aij| = %10.3e  ratio = %10.3e\n";
393      double min_aij, max_aij, ratio;
394      xprintf("Scaling...\n");
395      /* cancel the current scaling effect */
396      glp_unscale_prob(lp);
397      /* report original scaling "quality" */
398      min_aij = min_mat_aij(lp, 1);
399      max_aij = max_mat_aij(lp, 1);
400      ratio = max_aij / min_aij;
401      xprintf(fmt, " A", min_aij, max_aij, ratio);
402      /* check if the problem is well scaled */
403      if (min_aij >= 0.10 && max_aij <= 10.0)
404      {  xprintf("Problem data seem to be well scaled\n");
405         /* skip scaling, if required */
406         if (flags & GLP_SF_SKIP) goto done;
407      }
408      /* perform iterative geometric mean scaling, if required */
409      if (flags & GLP_SF_GM)
410      {  gm_iterate(lp, 15, 0.90);
411         min_aij = min_mat_aij(lp, 1);
412         max_aij = max_mat_aij(lp, 1);
413         ratio = max_aij / min_aij;
414         xprintf(fmt, "GM", min_aij, max_aij, ratio);
415      }
416      /* perform equilibration scaling, if required */
417      if (flags & GLP_SF_EQ)
418      {  eq_scaling(lp, max_row_ratio(lp) > max_col_ratio(lp));
419         min_aij = min_mat_aij(lp, 1);
420         max_aij = max_mat_aij(lp, 1);
421         ratio = max_aij / min_aij;
422         xprintf(fmt, "EQ", min_aij, max_aij, ratio);
423      }
424      /* round scale factors to nearest power of two, if required */
425      if (flags & GLP_SF_2N)
426      {  int i, j;
427         for (i = 1; i <= lp->m; i++)
428            glp_set_rii(lp, i, round2n(glp_get_rii(lp, i)));
429         for (j = 1; j <= lp->n; j++)
430            glp_set_sjj(lp, j, round2n(glp_get_sjj(lp, j)));
431         min_aij = min_mat_aij(lp, 1);
432         max_aij = max_mat_aij(lp, 1);
433         ratio = max_aij / min_aij;
434         xprintf(fmt, "2N", min_aij, max_aij, ratio);
435      }
436done: return;
437}
438
439/***********************************************************************
440*  NAME
441*
442*  glp_scale_prob - scale problem data
443*
444*  SYNOPSIS
445*
446*  void glp_scale_prob(glp_prob *lp, int flags);
447*
448*  DESCRIPTION
449*
450*  The routine glp_scale_prob performs automatic scaling of problem
451*  data for the specified problem object.
452*
453*  The parameter flags specifies scaling options used by the routine.
454*  Options can be combined with the bitwise OR operator and may be the
455*  following:
456*
457*  GLP_SF_GM      perform geometric mean scaling;
458*  GLP_SF_EQ      perform equilibration scaling;
459*  GLP_SF_2N      round scale factors to nearest power of two;
460*  GLP_SF_SKIP    skip scaling, if the problem is well scaled.
461*
462*  The parameter flags may be specified as GLP_SF_AUTO, in which case
463*  the routine chooses scaling options automatically. */
464
465void glp_scale_prob(glp_prob *lp, int flags)
466{     if (flags & ~(GLP_SF_GM | GLP_SF_EQ | GLP_SF_2N | GLP_SF_SKIP |
467                    GLP_SF_AUTO))
468         xerror("glp_scale_prob: flags = 0x%02X; invalid scaling option"
469            "s\n", flags);
470      if (flags & GLP_SF_AUTO)
471         flags = (GLP_SF_GM | GLP_SF_EQ | GLP_SF_SKIP);
472      scale_prob(lp, flags);
473      return;
474}
475
476/* eof */
Note: See TracBrowser for help on using the repository browser.