aboutsummaryrefslogtreecommitdiffstats
path: root/test/monniaux/glpk-4.65/src/api/advbas.c
blob: 23067624617170526e66754f699b14ad87daa133 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
/* advbas.c (construct advanced initial LP basis) */

/***********************************************************************
*  This code is part of GLPK (GNU Linear Programming Kit).
*
*  Copyright (C) 2008-2016 Andrew Makhorin, Department for Applied
*  Informatics, Moscow Aviation Institute, Moscow, Russia. All rights
*  reserved. E-mail: <mao@gnu.org>.
*
*  GLPK is free software: you can redistribute it and/or modify it
*  under the terms of the GNU General Public License as published by
*  the Free Software Foundation, either version 3 of the License, or
*  (at your option) any later version.
*
*  GLPK is distributed in the hope that it will be useful, but WITHOUT
*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
*  License for more details.
*
*  You should have received a copy of the GNU General Public License
*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
***********************************************************************/

#include "env.h"
#include "prob.h"
#include "triang.h"

/***********************************************************************
*  NAME
*
*  glp_adv_basis - construct advanced initial LP basis
*
*  SYNOPSIS
*
*  void glp_adv_basis(glp_prob *P, int flags);
*
*  DESCRIPTION
*
*  The routine glp_adv_basis constructs an advanced initial LP basis
*  for the specified problem object.
*
*  The parameter flag is reserved for use in the future and should be
*  specified as zero.
*
*  NOTE
*
*  The routine glp_adv_basis should be called after the constraint
*  matrix has been scaled (if scaling is used). */

static int mat(void *info, int k, int ind[], double val[])
{     glp_prob *P = info;
      int m = P->m;
      int n = P->n;
      GLPROW **row = P->row;
      GLPCOL **col = P->col;
      GLPAIJ *aij;
      int i, j, len;
      if (k > 0)
      {  /* retrieve scaled row of constraint matrix */
         i = +k;
         xassert(1 <= i && i <= m);
         len = 0;
         if (row[i]->type == GLP_FX)
         {  for (aij = row[i]->ptr; aij != NULL; aij = aij->r_next)
            {  j = aij->col->j;
               if (col[j]->type != GLP_FX)
               {  len++;
                  ind[len] = j;
                  val[len] = aij->row->rii * aij->val * aij->col->sjj;
               }
            }
         }
      }
      else
      {  /* retrieve scaled column of constraint matrix */
         j = -k;
         xassert(1 <= j && j <= n);
         len = 0;
         if (col[j]->type != GLP_FX)
         {  for (aij = col[j]->ptr; aij != NULL; aij = aij->c_next)
            {  i = aij->row->i;
               if (row[i]->type == GLP_FX)
               {  len++;
                  ind[len] = i;
                  val[len] = aij->row->rii * aij->val * aij->col->sjj;
               }
            }
         }
      }
      return len;
}

void glp_adv_basis(glp_prob *P, int flags)
{     int i, j, k, m, n, min_mn, size, *rn, *cn;
      char *flag;
      if (flags != 0)
         xerror("glp_adv_basis: flags = %d; invalid flags\n", flags);
      m = P->m; /* number of rows */
      n = P->n; /* number of columns */
      if (m == 0 || n == 0)
      {  /* trivial case */
         glp_std_basis(P);
         goto done;
      }
      xprintf("Constructing initial basis...\n");
      /* allocate working arrays */
      min_mn = (m < n ? m : n);
      rn = talloc(1+min_mn, int);
      cn = talloc(1+min_mn, int);
      flag = talloc(1+m, char);
      /* make the basis empty */
      for (i = 1; i <= m; i++)
      {  flag[i] = 0;
         glp_set_row_stat(P, i, GLP_NS);
      }
      for (j = 1; j <= n; j++)
         glp_set_col_stat(P, j, GLP_NS);
      /* find maximal triangular part of the constraint matrix;
         to prevent including non-fixed rows and fixed columns in the
         triangular part, such rows and columns are temporarily made
         empty by the routine mat */
#if 1 /* FIXME: tolerance */
      size = triang(m, n, mat, P, 0.001, rn, cn);
#endif
      xassert(0 <= size && size <= min_mn);
      /* include in the basis non-fixed structural variables, whose
         columns constitute the triangular part */
      for (k = 1; k <= size; k++)
      {  i = rn[k];
         xassert(1 <= i && i <= m);
         flag[i] = 1;
         j = cn[k];
         xassert(1 <= j && j <= n);
         glp_set_col_stat(P, j, GLP_BS);
      }
      /* include in the basis appropriate auxiliary variables, whose
         unity columns preserve triangular form of the basis matrix */
      for (i = 1; i <= m; i++)
      {  if (flag[i] == 0)
         {  glp_set_row_stat(P, i, GLP_BS);
            if (P->row[i]->type != GLP_FX)
               size++;
         }
      }
      /* size of triangular part = (number of rows) - (number of basic
         fixed auxiliary variables) */
      xprintf("Size of triangular part is %d\n", size);
      /* deallocate working arrays */
      tfree(rn);
      tfree(cn);
      tfree(flag);
done: return;
}

/* eof */