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
156
157
158
159
160
161
162
163
164
|
// This file is part of libigl, a simple c++ geometry processing library.
//
// Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
//
// This Source Code Form is subject to the terms of the Mozilla Public License
// v. 2.0. If a copy of the MPL was not distributed with this file, You can
// obtain one at http://mozilla.org/MPL/2.0/.
#include "mosek_linprog.h"
#include "../mosek/mosek_guarded.h"
#include "../harwell_boeing.h"
#include <limits>
#include <cmath>
#include <vector>
IGL_INLINE bool igl::mosek::mosek_linprog(
const Eigen::VectorXd & c,
const Eigen::SparseMatrix<double> & A,
const Eigen::VectorXd & lc,
const Eigen::VectorXd & uc,
const Eigen::VectorXd & lx,
const Eigen::VectorXd & ux,
Eigen::VectorXd & x)
{
// variables for mosek task, env and result code
MSKenv_t env;
// Create the MOSEK environment
mosek_guarded(MSK_makeenv(&env,NULL));
// initialize mosek environment
#if MSK_VERSION_MAJOR <= 7
mosek_guarded(MSK_initenv(env));
#endif
const bool ret = mosek_linprog(c,A,lc,uc,lx,ux,env,x);
MSK_deleteenv(&env);
return ret;
}
IGL_INLINE bool igl::mosek::mosek_linprog(
const Eigen::VectorXd & c,
const Eigen::SparseMatrix<double> & A,
const Eigen::VectorXd & lc,
const Eigen::VectorXd & uc,
const Eigen::VectorXd & lx,
const Eigen::VectorXd & ux,
const MSKenv_t & env,
Eigen::VectorXd & x)
{
// following http://docs.mosek.com/7.1/capi/Linear_optimization.html
using namespace std;
// number of constraints
const int m = A.rows();
// number of variables
const int n = A.cols();
vector<double> vAv;
vector<int> vAri,vAcp;
int nr;
harwell_boeing(A,nr,vAv,vAri,vAcp);
MSKtask_t task;
// Create the optimization task
mosek_guarded(MSK_maketask(env,m,n,&task));
// no threads
mosek_guarded(MSK_putintparam(task,MSK_IPAR_NUM_THREADS,1));
if(m>0)
{
// Append 'm' empty constraints, the constrainst will initially have no
// bounds
mosek_guarded(MSK_appendcons(task,m));
}
mosek_guarded(MSK_appendvars(task,n));
const auto & key = [](const double lxj, const double uxj) ->
MSKboundkeye
{
MSKboundkeye k = MSK_BK_FR;
if(isfinite(lxj) && isfinite(uxj))
{
if(lxj == uxj)
{
k = MSK_BK_FX;
}else{
k = MSK_BK_RA;
}
}else if(isfinite(lxj))
{
k = MSK_BK_LO;
}else if(isfinite(uxj))
{
k = MSK_BK_UP;
}
return k;
};
// loop over variables
for(int j = 0;j<n;j++)
{
if(c.size() > 0)
{
// Set linear term c_j in the objective
mosek_guarded(MSK_putcj(task,j,c(j)));
}
// Set constant bounds on variable j
const double lxj = lx.size()>0?lx[j]:-numeric_limits<double>::infinity();
const double uxj = ux.size()>0?ux[j]: numeric_limits<double>::infinity();
mosek_guarded(MSK_putvarbound(task,j,key(lxj,uxj),lxj,uxj));
if(m>0)
{
// Input column j of A
mosek_guarded(
MSK_putacol(
task,
j,
vAcp[j+1]-vAcp[j],
&vAri[vAcp[j]],
&vAv[vAcp[j]])
);
}
}
// loop over constraints
for(int i = 0;i<m;i++)
{
// Set constraint bounds for row i
const double lci = lc.size()>0?lc[i]:-numeric_limits<double>::infinity();
const double uci = uc.size()>0?uc[i]: numeric_limits<double>::infinity();
mosek_guarded(MSK_putconbound(task,i,key(lci,uci),lci,uci));
}
// Now the optimizer has been prepared
MSKrescodee trmcode;
// run the optimizer
mosek_guarded(MSK_optimizetrm(task,&trmcode));
// Get status
MSKsolstae solsta;
MSK_getsolsta (task,MSK_SOL_ITR,&solsta);
bool success = false;
switch(solsta)
{
case MSK_SOL_STA_OPTIMAL:
case MSK_SOL_STA_NEAR_OPTIMAL:
x.resize(n);
/* Request the basic solution. */
MSK_getxx(task,MSK_SOL_BAS,x.data());
success = true;
break;
case MSK_SOL_STA_DUAL_INFEAS_CER:
case MSK_SOL_STA_PRIM_INFEAS_CER:
case MSK_SOL_STA_NEAR_DUAL_INFEAS_CER:
case MSK_SOL_STA_NEAR_PRIM_INFEAS_CER:
//printf("Primal or dual infeasibility certificate found.\n");
break;
case MSK_SOL_STA_UNKNOWN:
//printf("The status of the solution could not be determined.\n");
break;
default:
//printf("Other solution status.");
break;
}
MSK_deletetask(&task);
return success;
}
|
