1 /* prrngs.c (print sensitivity analysis report) */
2 
3 /***********************************************************************
4 *  This code is part of GLPK (GNU Linear Programming Kit).
5 *  Copyright (C) 2009-2016 Free Software Foundation, Inc.
6 *  Written by Andrew Makhorin <mao@gnu.org>.
7 *
8 *  GLPK is free software: you can redistribute it and/or modify it
9 *  under the terms of the GNU General Public License as published by
10 *  the Free Software Foundation, either version 3 of the License, or
11 *  (at your option) any later version.
12 *
13 *  GLPK is distributed in the hope that it will be useful, but WITHOUT
14 *  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 *  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
16 *  License for more details.
17 *
18 *  You should have received a copy of the GNU General Public License
19 *  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
20 ***********************************************************************/
21 
22 #include "env.h"
23 #include "prob.h"
24 
25 #define xfprintf glp_format
26 
format(char buf[13+1],double x)27 static char *format(char buf[13+1], double x)
28 {     /* format floating-point number in MPS/360-like style */
29       if (x == -DBL_MAX)
30          strcpy(buf, "         -Inf");
31       else if (x == +DBL_MAX)
32          strcpy(buf, "         +Inf");
33       else if (fabs(x) <= 999999.99998)
34       {  sprintf(buf, "%13.5f", x);
35 #if 1
36          if (strcmp(buf, "      0.00000") == 0 ||
37              strcmp(buf, "     -0.00000") == 0)
38             strcpy(buf, "       .     ");
39          else if (memcmp(buf, "      0.", 8) == 0)
40             memcpy(buf, "       .", 8);
41          else if (memcmp(buf, "     -0.", 8) == 0)
42             memcpy(buf, "      -.", 8);
43 #endif
44       }
45       else
46          sprintf(buf, "%13.6g", x);
47       return buf;
48 }
49 
glp_print_ranges(glp_prob * P,int len,const int list[],int flags,const char * fname)50 int glp_print_ranges(glp_prob *P, int len, const int list[],
51       int flags, const char *fname)
52 {     /* print sensitivity analysis report */
53       glp_file *fp = NULL;
54       GLPROW *row;
55       GLPCOL *col;
56       int m, n, pass, k, t, numb, type, stat, var1, var2, count, page,
57          ret;
58       double lb, ub, slack, coef, prim, dual, value1, value2, coef1,
59          coef2, obj1, obj2;
60       const char *name, *limit;
61       char buf[13+1];
62       /* sanity checks */
63 #if 0 /* 04/IV-2016 */
64       if (P == NULL || P->magic != GLP_PROB_MAGIC)
65          xerror("glp_print_ranges: P = %p; invalid problem object\n",
66             P);
67 #endif
68       m = P->m, n = P->n;
69       if (len < 0)
70          xerror("glp_print_ranges: len = %d; invalid list length\n",
71             len);
72       if (len > 0)
73       {  if (list == NULL)
74             xerror("glp_print_ranges: list = %p: invalid parameter\n",
75                list);
76          for (t = 1; t <= len; t++)
77          {  k = list[t];
78             if (!(1 <= k && k <= m+n))
79                xerror("glp_print_ranges: list[%d] = %d; row/column numb"
80                   "er out of range\n", t, k);
81          }
82       }
83       if (flags != 0)
84          xerror("glp_print_ranges: flags = %d; invalid parameter\n",
85             flags);
86       if (fname == NULL)
87          xerror("glp_print_ranges: fname = %p; invalid parameter\n",
88             fname);
89       if (glp_get_status(P) != GLP_OPT)
90       {  xprintf("glp_print_ranges: optimal basic solution required\n");
91          ret = 1;
92          goto done;
93       }
94       if (!glp_bf_exists(P))
95       {  xprintf("glp_print_ranges: basis factorization required\n");
96          ret = 2;
97          goto done;
98       }
99       /* start reporting */
100       xprintf("Write sensitivity analysis report to '%s'...\n", fname);
101       fp = glp_open(fname, "w");
102       if (fp == NULL)
103       {  xprintf("Unable to create '%s' - %s\n", fname, get_err_msg());
104          ret = 3;
105          goto done;
106       }
107       page = count = 0;
108       for (pass = 1; pass <= 2; pass++)
109       for (t = 1; t <= (len == 0 ? m+n : len); t++)
110       {  if (t == 1) count = 0;
111          k = (len == 0 ? t : list[t]);
112          if (pass == 1 && k > m || pass == 2 && k <= m)
113             continue;
114          if (count == 0)
115          {  xfprintf(fp, "GLPK %-4s - SENSITIVITY ANALYSIS REPORT%73sPa"
116                "ge%4d\n", glp_version(), "", ++page);
117             xfprintf(fp, "\n");
118             xfprintf(fp, "%-12s%s\n", "Problem:",
119                P->name == NULL ? "" : P->name);
120             xfprintf(fp, "%-12s%s%s%.10g (%s)\n", "Objective:",
121                P->obj == NULL ? "" : P->obj,
122                P->obj == NULL ? "" : " = ", P->obj_val,
123                P->dir == GLP_MIN ? "MINimum" :
124                P->dir == GLP_MAX ? "MAXimum" : "???");
125             xfprintf(fp, "\n");
126             xfprintf(fp, "%6s %-12s %2s %13s %13s %13s  %13s %13s %13s "
127                "%s\n", "No.", pass == 1 ? "Row name" : "Column name",
128                "St", "Activity", pass == 1 ? "Slack" : "Obj coef",
129                "Lower bound", "Activity", "Obj coef", "Obj value at",
130                "Limiting");
131             xfprintf(fp, "%6s %-12s %2s %13s %13s %13s  %13s %13s %13s "
132                "%s\n", "", "", "", "", "Marginal", "Upper bound",
133                "range", "range", "break point", "variable");
134             xfprintf(fp, "------ ------------ -- ------------- --------"
135                "----- -------------  ------------- ------------- ------"
136                "------- ------------\n");
137          }
138          if (pass == 1)
139          {  numb = k;
140             xassert(1 <= numb && numb <= m);
141             row = P->row[numb];
142             name = row->name;
143             type = row->type;
144             lb = glp_get_row_lb(P, numb);
145             ub = glp_get_row_ub(P, numb);
146             coef = 0.0;
147             stat = row->stat;
148             prim = row->prim;
149             if (type == GLP_FR)
150                slack = - prim;
151             else if (type == GLP_LO)
152                slack = lb - prim;
153             else if (type == GLP_UP || type == GLP_DB || type == GLP_FX)
154                slack = ub - prim;
155             dual = row->dual;
156          }
157          else
158          {  numb = k - m;
159             xassert(1 <= numb && numb <= n);
160             col = P->col[numb];
161             name = col->name;
162             lb = glp_get_col_lb(P, numb);
163             ub = glp_get_col_ub(P, numb);
164             coef = col->coef;
165             stat = col->stat;
166             prim = col->prim;
167             slack = 0.0;
168             dual = col->dual;
169          }
170          if (stat != GLP_BS)
171          {  glp_analyze_bound(P, k, &value1, &var1, &value2, &var2);
172             if (stat == GLP_NF)
173                coef1 = coef2 = coef;
174             else if (stat == GLP_NS)
175                coef1 = -DBL_MAX, coef2 = +DBL_MAX;
176             else if (stat == GLP_NL && P->dir == GLP_MIN ||
177                      stat == GLP_NU && P->dir == GLP_MAX)
178                coef1 = coef - dual, coef2 = +DBL_MAX;
179             else
180                coef1 = -DBL_MAX, coef2 = coef - dual;
181             if (value1 == -DBL_MAX)
182             {  if (dual < -1e-9)
183                   obj1 = +DBL_MAX;
184                else if (dual > +1e-9)
185                   obj1 = -DBL_MAX;
186                else
187                   obj1 = P->obj_val;
188             }
189             else
190                obj1 = P->obj_val + dual * (value1 - prim);
191             if (value2 == +DBL_MAX)
192             {  if (dual < -1e-9)
193                   obj2 = -DBL_MAX;
194                else if (dual > +1e-9)
195                   obj2 = +DBL_MAX;
196                else
197                   obj2 = P->obj_val;
198             }
199             else
200                obj2 = P->obj_val + dual * (value2 - prim);
201          }
202          else
203          {  glp_analyze_coef(P, k, &coef1, &var1, &value1, &coef2,
204                &var2, &value2);
205             if (coef1 == -DBL_MAX)
206             {  if (prim < -1e-9)
207                   obj1 = +DBL_MAX;
208                else if (prim > +1e-9)
209                   obj1 = -DBL_MAX;
210                else
211                   obj1 = P->obj_val;
212             }
213             else
214                obj1 = P->obj_val + (coef1 - coef) * prim;
215             if (coef2 == +DBL_MAX)
216             {  if (prim < -1e-9)
217                   obj2 = -DBL_MAX;
218                else if (prim > +1e-9)
219                   obj2 = +DBL_MAX;
220                else
221                   obj2 = P->obj_val;
222             }
223             else
224                obj2 = P->obj_val + (coef2 - coef) * prim;
225          }
226          /*** first line ***/
227          /* row/column number */
228          xfprintf(fp, "%6d", numb);
229          /* row/column name */
230          xfprintf(fp, " %-12.12s", name == NULL ? "" : name);
231          if (name != NULL && strlen(name) > 12)
232             xfprintf(fp, "%s\n%6s %12s", name+12, "", "");
233          /* row/column status */
234          xfprintf(fp, " %2s",
235             stat == GLP_BS ? "BS" : stat == GLP_NL ? "NL" :
236             stat == GLP_NU ? "NU" : stat == GLP_NF ? "NF" :
237             stat == GLP_NS ? "NS" : "??");
238          /* row/column activity */
239          xfprintf(fp, " %s", format(buf, prim));
240          /* row slack, column objective coefficient */
241          xfprintf(fp, " %s", format(buf, k <= m ? slack : coef));
242          /* row/column lower bound */
243          xfprintf(fp, " %s", format(buf, lb));
244          /* row/column activity range */
245          xfprintf(fp, "  %s", format(buf, value1));
246          /* row/column objective coefficient range */
247          xfprintf(fp, " %s", format(buf, coef1));
248          /* objective value at break point */
249          xfprintf(fp, " %s", format(buf, obj1));
250          /* limiting variable name */
251          if (var1 != 0)
252          {  if (var1 <= m)
253                limit = glp_get_row_name(P, var1);
254             else
255                limit = glp_get_col_name(P, var1 - m);
256             if (limit != NULL)
257                xfprintf(fp, " %s", limit);
258          }
259          xfprintf(fp, "\n");
260          /*** second line ***/
261          xfprintf(fp, "%6s %-12s %2s %13s", "", "", "", "");
262          /* row/column reduced cost */
263          xfprintf(fp, " %s", format(buf, dual));
264          /* row/column upper bound */
265          xfprintf(fp, " %s", format(buf, ub));
266          /* row/column activity range */
267          xfprintf(fp, "  %s", format(buf, value2));
268          /* row/column objective coefficient range */
269          xfprintf(fp, " %s", format(buf, coef2));
270          /* objective value at break point */
271          xfprintf(fp, " %s", format(buf, obj2));
272          /* limiting variable name */
273          if (var2 != 0)
274          {  if (var2 <= m)
275                limit = glp_get_row_name(P, var2);
276             else
277                limit = glp_get_col_name(P, var2 - m);
278             if (limit != NULL)
279                xfprintf(fp, " %s", limit);
280          }
281          xfprintf(fp, "\n");
282          xfprintf(fp, "\n");
283          /* print 10 items per page */
284          count = (count + 1) % 10;
285       }
286       xfprintf(fp, "End of report\n");
287 #if 0 /* FIXME */
288       xfflush(fp);
289 #endif
290       if (glp_ioerr(fp))
291       {  xprintf("Write error on '%s' - %s\n", fname, get_err_msg());
292          ret = 4;
293          goto done;
294       }
295       ret = 0;
296 done: if (fp != NULL) glp_close(fp);
297       return ret;
298 }
299 
300 /* eof */
301