[cvxopt] 21/64: Imported Upstream version 1.0
Andreas Tille
tille at debian.org
Wed Jul 20 11:23:50 UTC 2016
This is an automated email from the git hooks/post-receive script.
tille pushed a commit to branch master
in repository cvxopt.
commit 57e4df59e532c047fafde71dd540410ba530bff1
Author: Andreas Tille <tille at debian.org>
Date: Wed Jul 20 08:26:55 2016 +0200
Imported Upstream version 1.0
---
INSTALL | 4 +-
LICENSE | 2 +-
doc/base.tex | 18 +-
doc/base_sparse.tex | 16 +-
doc/coneprog.tex | 8 +-
doc/cvxopt.tex | 2 +-
doc/fftw.tex | 45 ++
doc/figures/floorplan.pdf | Bin 17887 -> 17885 bytes
doc/figures/normappr.pdf | Bin 21648 -> 21648 bytes
doc/figures/portfolio1.pdf | Bin 21164 -> 21164 bytes
doc/figures/portfolio2.pdf | Bin 22818 -> 22820 bytes
doc/intro.tex | 2 +-
doc/spsolvers.tex | 2 +
src/C/amd.c | 4 +-
src/C/base.c | 2 +-
src/C/blas.c | 2 +-
src/C/cholmod.c | 4 +-
src/C/cvxopt.h | 23 +-
src/C/dense.c | 44 +-
src/C/dsdp.c | 4 +-
src/C/fftw.c | 257 +++++++++++-
src/C/glpk.c | 335 ++++++++++++++-
src/C/gsl.c | 2 +-
src/C/lapack.c | 2 +-
src/C/misc.h | 2 +-
src/C/misc_solvers.c | 1004 ++++++++++++++++++++++++++++++++++++++++++++
src/C/sparse.c | 55 ++-
src/C/umfpack.c | 2 +-
src/python/__init__.py | 4 +-
src/python/coneprog.py | 54 ++-
src/python/cvxprog.py | 22 +-
src/python/info.py | 4 +-
src/python/misc.py | 615 +++++++++++++++------------
src/python/modeling.py | 2 +-
src/python/mosek.py | 180 +++++++-
src/python/printing.py | 3 +-
src/python/solvers.py | 2 +-
src/setup.py | 9 +-
38 files changed, 2359 insertions(+), 377 deletions(-)
diff --git a/INSTALL b/INSTALL
index f70676c..33b1616 100644
--- a/INSTALL
+++ b/INSTALL
@@ -1,6 +1,6 @@
-Installation instructions for CVXOPT Version 0.9.3.
+Installation instructions for CVXOPT Version 1.0.
-The package requires version 2.3 or newer of Python, and is built from
+The package requires version 2.5 or newer of Python, and is built from
source, so the header files and libraries for Python must be installed,
as well as the core binaries.
diff --git a/LICENSE b/LICENSE
index 9833a5e..f49560f 100644
--- a/LICENSE
+++ b/LICENSE
@@ -1,4 +1,4 @@
-CVXOPT version 0.9.3. Copyright (c) 2004-2008 J. Dahl and L. Vandenberghe.
+CVXOPT version 1.0. Copyright (c) 2004-2008 J. Dahl and L. Vandenberghe.
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
diff --git a/doc/base.tex b/doc/base.tex
index d602cce..6439b25 100644
--- a/doc/base.tex
+++ b/doc/base.tex
@@ -155,9 +155,9 @@ list (\ie, when the length of \var{x} is one, it can be replaced with
A \mtrx\ has the following attributes.
\begin{memberdesc}{size}
-A tuple with the dimensions of the matrix. This is a read-only
-attribute; operations that change the size of a matrix are not
-permitted.
+A tuple with the dimensions of the matrix. The size of the matrix
+can be changed by altering this attribute, as long as the number
+of elements in the matrix remains unchanged.
\end{memberdesc}
\begin{memberdesc}{typecode}
@@ -284,7 +284,7 @@ Note that Python rounds the result of an integer division towards minus
infinity.
The following in-place operations are also defined, but only if
-they do not change the type or the size of the matrix \var{A}:
+they do not change the type of the matrix \var{A}:
\begin{center}\begin{tabular}{l|l} \hline
In-place addition &
\code{\var{A}+=\var{B}}, \code{\var{A}+=\var{c}} \\ \hline
@@ -352,11 +352,11 @@ return new objects. Hence \samp{A = +B} is equivalent to
\end{verbatim}
\end{itemize}
-The restrictions on in-place operations follow the principle that once
-a matrix object is created, its size and type cannot be modified.
-The only matrix attributes that can be changed are the values of the
-elements. The values can be changed by in-place operations or by an
-indexed assignment, as explained in the next section.
+%The restrictions on in-place operations follow the principle that once
+%a matrix object is created, its size and type cannot be modified.
+%The only matrix attributes that can be changed are the values of the
+%elements. The values can be changed by in-place operations or by an
+%indexed assignment, as explained in the next section.
\section{Indexing and Slicing} \label{s-indexing}
diff --git a/doc/base_sparse.tex b/doc/base_sparse.tex
index 6cb65ee..254a4b1 100644
--- a/doc/base_sparse.tex
+++ b/doc/base_sparse.tex
@@ -228,7 +228,9 @@ in \code{V}. A read-only attribute.
\end{memberdesc}
\begin{memberdesc}{size}
-A tuple with the dimensions of the matrix. A read-only attribute.
+A tuple with the dimensions of the matrix. The size of the matrix
+can be changed by altering this attribute, as long as the number
+of elements in the matrix remains unchanged.
\end{memberdesc}
\begin{memberdesc}{CCS}
@@ -419,12 +421,12 @@ is interpreted as a scalar product \code{\var{A}*=\var{c}[0]}.)
As for dense operations, the in-place sparse operations do not return
a new matrix but modify the existing object \var{A}.
-The restrictions on in-place operations follow the principle that once
-a sparse matrix is created, its size and type cannot be modified.
-The only attributes that can be modified are the sparsity pattern and
-the numerical values of the nonzero elements.
-These attributes can be modified by in-place operations or by indexed
-assignments.
+%The restrictions on in-place operations follow the principle that once
+%a sparse matrix is created, its size and type cannot be modified.
+%The only attributes that can be modified are the sparsity pattern and
+%the numerical values of the nonzero elements.
+%These attributes can be modified by in-place operations or by indexed
+%assignments.
diff --git a/doc/coneprog.tex b/doc/coneprog.tex
index 1ee7bec..25b2477 100644
--- a/doc/coneprog.tex
+++ b/doc/coneprog.tex
@@ -32,7 +32,7 @@ and the algorithm parameters that control the cone programming solvers.
\begin{funcdesc}{conelp}{c, G, h\optional{, dims\optional{, A, b\optional{,
primalstart\optional{, dualstart\optional{, kktsolver}}}}}}
Solves a pair of primal and dual cone programs
-\BEQ \label{e-conelp}
+\[ %\label{e-conelp}
\begin{array}[t]{ll}
\mbox{minimize} & c^T x \\
\mbox{subject to} & G x + s = h \\ & Ax = b \\ & s \succeq 0
@@ -43,7 +43,7 @@ Solves a pair of primal and dual cone programs
\mbox{subject to} & G^T z + A^T y + c = 0 \\
& z \succeq 0.
\end{array}
-\EEQ
+\]
The primal variables are \tm{x} and the slack variable \tm{s}.
The dual variables are \tm{y} and \tm{z}. The inequalities are
interpreted as $s \in C$, $z\in C$, where \tm{C} is a cone defined as a
@@ -238,7 +238,7 @@ below.
\begin{funcdesc}{coneqp}{P, q\optional{, G, h\optional{,
dims\optional{, A, b\optional{, initvals\optional{, kktsolver}}}}}}
Solves a pair of primal and dual quadratic cone programs
-\BEQ \label{e-conelp}
+\[ %\label{e-coneqp}
\begin{array}[t]{ll}
\mbox{minimize} & (1/2) x^T Px + q^T x \\
\mbox{subject to} & G x + s = h \\ & Ax = b \\ & s \succeq 0
@@ -249,7 +249,7 @@ Solves a pair of primal and dual quadratic cone programs
(q+G^Tz+A^Ty) -h^T z - b^T y \\
\mbox{subject to} & q + G^T z + A^T y \in \Range(P) \\ & z \succeq 0.
\end{array}
-\EEQ
+\]
The primal variables are \tm{x} and the slack variable \tm{s}.
The dual variables are \tm{y} and \tm{z}. The inequalities are
interpreted as $s \in C$, $z\in C$, where \tm{C} is a cone defined as a
diff --git a/doc/cvxopt.tex b/doc/cvxopt.tex
index b265095..ab02ca9 100644
--- a/doc/cvxopt.tex
+++ b/doc/cvxopt.tex
@@ -82,7 +82,7 @@
\title{CVXOPT User's Guide}
\author{Joachim Dahl \& Lieven Vandenberghe}
-\date{Release 0.9.3 -- February 24, 2008}
+\date{Release 1.0 -- April 24, 2008}
\begin{document}
\Configure{crosslinks*}{next}{prev}{up}{}
diff --git a/doc/fftw.tex b/doc/fftw.tex
index 2f3ff6c..f5e0f9e 100644
--- a/doc/fftw.tex
+++ b/doc/fftw.tex
@@ -29,6 +29,21 @@ discrete Fourier transforms: if \var{X} has {\it n} rows,
\]
\end{funcdesc}
+The separable discrete two dimensional Fourier transform first
+computes the corresponding one dimensional tranform along the columns
+of the matrix, followed by the one dimensional transform along the
+rows of the matrix.
+
+\begin{funcdesc}{dft2}{X}
+Replaces a dense complex matrix with the two dimensional
+discrete Fourier transform.
+\end{funcdesc}
+
+\begin{funcdesc}{idft2}{X}
+Replaces a dense complex matrix with the inverse two dimensional
+discrete Fourier transform.
+\end{funcdesc}
+
\section{Discrete Cosine Transform}
\begin{funcdesc}{dct}{X\optional{, type=2}}
Replaces the columns of a dense real matrix with their discrete
@@ -60,6 +75,21 @@ Replaces the columns of a dense real matrix with the inverses
of the discrete cosine transforms defined above.
\end{funcdesc}
+The separable discrete two dimensional cosine transform first
+computes the corresponding one dimensional tranform along the columns
+of the matrix, followed by the one dimensional transform along the
+rows of the matrix.
+
+\begin{funcdesc}{dct2}{X\optional{, type=2}}
+Replaces a dense real matrix with the two dimensional
+discrete cosine transform.
+\end{funcdesc}
+
+\begin{funcdesc}{idct2}{X\optional{, type=2}}
+Replaces a dense real matrix with the inverse two dimensional
+discrete cosine transform.
+\end{funcdesc}
+
\section{Discrete Sine Transform}
\begin{funcdesc}{dst}{X\optional{, type=1}}
Replaces the columns of a dense real matrix with their discrete
@@ -89,3 +119,18 @@ These transforms are defined as follows
Replaces the columns of a dense real matrix with the inverses
of the discrete sine transforms defined above.
\end{funcdesc}
+
+The separable discrete two dimensional sine transform first
+computes the corresponding one dimensional tranform along the columns
+of the matrix, followed by the one dimensional transform along the
+rows of the matrix.
+
+\begin{funcdesc}{dst2}{X\optional{, type=1}}
+Replaces a dense real matrix with the two dimensional
+discrete sine transform.
+\end{funcdesc}
+
+\begin{funcdesc}{idst2}{X\optional{, type=1}}
+Replaces a dense real matrix with the inverse two dimensional
+discrete sine transform.
+\end{funcdesc}
diff --git a/doc/figures/floorplan.pdf b/doc/figures/floorplan.pdf
index 9852f21..c211716 100644
Binary files a/doc/figures/floorplan.pdf and b/doc/figures/floorplan.pdf differ
diff --git a/doc/figures/normappr.pdf b/doc/figures/normappr.pdf
index d1a9b8c..d726169 100644
Binary files a/doc/figures/normappr.pdf and b/doc/figures/normappr.pdf differ
diff --git a/doc/figures/portfolio1.pdf b/doc/figures/portfolio1.pdf
index efebecf..82c0312 100644
Binary files a/doc/figures/portfolio1.pdf and b/doc/figures/portfolio1.pdf differ
diff --git a/doc/figures/portfolio2.pdf b/doc/figures/portfolio2.pdf
index 0a603f3..79e2450 100644
Binary files a/doc/figures/portfolio2.pdf and b/doc/figures/portfolio2.pdf differ
diff --git a/doc/intro.tex b/doc/intro.tex
index 1209880..42b2d58 100644
--- a/doc/intro.tex
+++ b/doc/intro.tex
@@ -10,7 +10,7 @@ optimization applications straightforward by building on Python's
extensive standard library and on the strengths of Python as a
high-level programming language.
-Release 0.9.3 of CVXOPT includes routines for basic linear algebra
+Release 1.0 of CVXOPT includes routines for basic linear algebra
calculations, interfaces to efficient libraries for solving dense and
sparse linear equations, convex optimization solvers written in Python,
interfaces to a few other optimization libraries,
diff --git a/doc/spsolvers.tex b/doc/spsolvers.tex
index cada4df..41514e2 100644
--- a/doc/spsolvers.tex
+++ b/doc/spsolvers.tex
@@ -302,6 +302,8 @@ is accessed and the upper triangular part is ignored.
If \var{uplo} is \code{'U'}, only the upper triangular part of \var{A}
is accessed and the lower triangular part is ignored.
+If \var{p} is not provided, a reordering from the AMD library is used.
+
The symbolic factorization is returned as an opaque C object that
can be passed to \function{cholmod.numeric()}.
\end{funcdesc}
diff --git a/src/C/amd.c b/src/C/amd.c
index b7fddea..ccda063 100644
--- a/src/C/amd.c
+++ b/src/C/amd.c
@@ -1,7 +1,7 @@
/*
* Copyright 2004-2008 J. Dahl and L. Vandenberghe.
*
- * This file is part of CVXOPT version 0.9.3.
+ * This file is part of CVXOPT version 1.0.
*
* CVXOPT is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
@@ -65,7 +65,7 @@ static int get_param_idx(char *str, int *idx)
static int set_defaults(double *control)
{
- int_compat pos=0;
+ int_t pos=0;
int param_id;
PyObject *param, *key, *value;
char *keystr, err_str[100];
diff --git a/src/C/base.c b/src/C/base.c
index 25318cf..3ceeba3 100644
--- a/src/C/base.c
+++ b/src/C/base.c
@@ -1,7 +1,7 @@
/*
* Copyright 2004-2008 J. Dahl and L. Vandenberghe.
*
- * This file is part of CVXOPT version 0.9.3.
+ * This file is part of CVXOPT version 1.0.
*
* CVXOPT is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
diff --git a/src/C/blas.c b/src/C/blas.c
index 8798cc0..0a38fc3 100644
--- a/src/C/blas.c
+++ b/src/C/blas.c
@@ -1,7 +1,7 @@
/*
* Copyright 2004-2008 J. Dahl and L. Vandenberghe.
*
- * This file is part of CVXOPT version 0.9.3.
+ * This file is part of CVXOPT version 1.0.
*
* CVXOPT is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
diff --git a/src/C/cholmod.c b/src/C/cholmod.c
index 9b60615..92f57a0 100644
--- a/src/C/cholmod.c
+++ b/src/C/cholmod.c
@@ -1,7 +1,7 @@
/*
* Copyright 2004-2008 J. Dahl and L. Vandenberghe.
*
- * This file is part of CVXOPT version 0.9.3.
+ * This file is part of CVXOPT version 1.0.
*
* CVXOPT is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
@@ -69,7 +69,7 @@ static cholmod_common Common;
static int set_options(void)
{
- int_compat pos=0;
+ int_t pos=0;
PyObject *param, *key, *value;
char *keystr, err_str[100];
diff --git a/src/C/cvxopt.h b/src/C/cvxopt.h
index 7d85a07..61f6161 100644
--- a/src/C/cvxopt.h
+++ b/src/C/cvxopt.h
@@ -1,7 +1,7 @@
/*
* Copyright 2004-2008 J. Dahl and L. Vandenberghe.
*
- * This file is part of CVXOPT version 0.9.3.
+ * This file is part of CVXOPT version 1.0.
*
* CVXOPT is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
@@ -28,23 +28,26 @@
#define MAT_BUFZ(O) ((complex *)((matrix *)O)->buffer)
#endif
-
#ifndef __CVXOPT__
#define __CVXOPT__
+/*
+static void * malloc_aligned(size_t size) {
+ void *p = NULL;
+ if (!posix_memalign(&p, 16, size))
+ return p;
+ else
+ return NULL;
+}
+
+#define malloc(size) malloc_aligned(size)
+*/
+
#define INT 0
#define DOUBLE 1
#define COMPLEX 2
-/* compatibility between Python2.5 and lower versions */
-#if PY_VERSION_HEX < 0x02050000
-#define int_t Py_intptr_t
-#define lenfunc inquiry
-#define int_compat int
-#else
#define int_t Py_ssize_t
-#define int_compat Py_ssize_t
-#endif
typedef struct {
PyObject_HEAD
diff --git a/src/C/dense.c b/src/C/dense.c
index 49d9bb7..2ebf343 100644
--- a/src/C/dense.c
+++ b/src/C/dense.c
@@ -1,7 +1,7 @@
/*
* Copyright 2004-2008 J. Dahl and L. Vandenberghe.
*
- * This file is part of CVXOPT version 0.9.3.
+ * This file is part of CVXOPT version 1.0.
*
* CVXOPT is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
@@ -25,7 +25,6 @@
#include <complexobject.h>
-
/* NumPy array protocol */
typedef struct {
int version;
@@ -575,7 +574,7 @@ matrix * create_indexlist(int dim, PyObject *A)
}
/* slice */
else if (PySlice_Check(A)) {
- int_compat start, stop, step, lgt;
+ int_t start, stop, step, lgt;
if (PySlice_GetIndicesEx((PySliceObject*)A, dim,
&start, &stop, &step, &lgt) < 0) return NULL;
@@ -656,8 +655,8 @@ matrix_subscr(matrix* self, PyObject* args)
/* two slices, handled separately for speed */
if (PySlice_Check(argI) && PySlice_Check(argJ)) {
- int_compat rowstart, rowstop, rowstep, rowlgt;
- int_compat colstart, colstop, colstep, collgt;
+ int_t rowstart, rowstop, rowstep, rowlgt;
+ int_t colstart, colstop, colstep, collgt;
if ( (PySlice_GetIndicesEx((PySliceObject*)argI, MAT_NROWS(self),
&rowstart, &rowstop, &rowstep, &rowlgt) < 0) ||
@@ -793,8 +792,8 @@ matrix_ass_subscr(matrix* self, PyObject* args, PyObject* val)
if (PySlice_Check(argI) && PySlice_Check(argJ) &&
Matrix_Check(val) && MAT_ID(val) == MAT_ID(self)) {
- int_compat rowstart, rowstop, rowstep, rowlgt;
- int_compat colstart, colstop, colstep, collgt;
+ int_t rowstart, rowstop, rowstep, rowlgt;
+ int_t colstart, colstop, colstep, collgt;
if ( (PySlice_GetIndicesEx((PySliceObject*)argI, MAT_NROWS(self),
&rowstart, &rowstop, &rowstep, &rowlgt) < 0) ||
@@ -1067,6 +1066,34 @@ matrix_get_size(matrix *self, void *closure)
return t;
}
+static int
+matrix_set_size(matrix *self, PyObject *value, void *closure)
+{
+ if (value == NULL)
+ PY_ERR_INT(PyExc_TypeError, "size attribute cannot be deleted");
+
+ if (!PyTuple_Check(value) || PyTuple_Size(value) != 2)
+ PY_ERR_INT(PyExc_TypeError, "can only assign a 2-tuple to size");
+
+ if (!PyInt_Check(PyTuple_GET_ITEM(value, 0)) ||
+ !PyInt_Check(PyTuple_GET_ITEM(value, 1)))
+ PY_ERR_INT(PyExc_TypeError, "invalid size tuple");
+
+ int m = PyInt_AS_LONG(PyTuple_GET_ITEM(value, 0));
+ int n = PyInt_AS_LONG(PyTuple_GET_ITEM(value, 1));
+
+ if (m<0 || n<0)
+ PY_ERR_INT(PyExc_TypeError, "dimensions must be non-negative");
+
+ if (m*n != MAT_LGT(self))
+ PY_ERR_INT(PyExc_TypeError, "number of elements in matrix cannot change");
+
+ MAT_NROWS(self) = m;
+ MAT_NCOLS(self) = n;
+
+ return 0;
+}
+
static PyObject * matrix_get_typecode(matrix *self, void *closure)
{
return PyString_FromStringAndSize(TC_CHAR[self->id], 1);
@@ -1119,7 +1146,8 @@ static PyObject * matrix_get_H(matrix *self, void *closure)
}
static PyGetSetDef matrix_getsets[] = {
- {"size", (getter) matrix_get_size, NULL, "matrix dimensions"},
+ {"size", (getter) matrix_get_size, (setter) matrix_set_size,
+ "matrix dimensions"},
{"typecode", (getter) matrix_get_typecode, NULL, "typecode character"},
{"__array_struct__", (getter) matrix_array_struct, NULL,
"C object implementing the NumPy array protocol"},
diff --git a/src/C/dsdp.c b/src/C/dsdp.c
index 10f8200..fefa332 100644
--- a/src/C/dsdp.c
+++ b/src/C/dsdp.c
@@ -1,7 +1,7 @@
/*
* Copyright 2004-2008 J. Dahl and L. Vandenberghe.
*
- * This file is part of CVXOPT version 0.9.3.
+ * This file is part of CVXOPT version 1.0.
*
* CVXOPT is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
@@ -110,7 +110,7 @@ static PyObject* solvesdp(PyObject *self, PyObject *args,
*param, *key, *value;
int i, j, k, n, ml, l, mk, nnz, *lp_colptr=NULL, *lp_rowind=NULL,
incx, incy, lngth, maxm;
- int_compat pos=0;
+ int_t pos=0;
double *lp_values=NULL, *zlvals=NULL, *zk=NULL, r, beta=-1.0,
gamma=-1.0, tol;
dsdp_matrix **lmis=NULL;
diff --git a/src/C/fftw.c b/src/C/fftw.c
index 5080acf..de7eea9 100644
--- a/src/C/fftw.c
+++ b/src/C/fftw.c
@@ -1,7 +1,7 @@
/*
* Copyright 2004-2008 J. Dahl and L. Vandenberghe.
*
- * This file is part of CVXOPT version 0.9.3.
+ * This file is part of CVXOPT version 1.0.
*
* CVXOPT is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
@@ -34,7 +34,6 @@ static char doc_dft[] =
"ARGUMENTS\n"
"X A dense matrix of typecode 'z'.";
-
static PyObject *dft(PyObject *self, PyObject *args, PyObject *kwrds)
{
matrix *X;
@@ -59,6 +58,35 @@ static PyObject *dft(PyObject *self, PyObject *args, PyObject *kwrds)
return Py_BuildValue("");
}
+static char doc_dft2[] =
+ "2D-DFT of a matrix, X := dft2(X)\n\n"
+ "PURPOSE\n"
+ "Computes the 2D-DFT of a dense matrix X.\n\n"
+ "ARGUMENTS\n"
+ "X A dense matrix of typecode 'z'.";
+
+static PyObject *dft2(PyObject *self, PyObject *args, PyObject *kwrds)
+{
+ matrix *X;
+ char *kwlist[] = {"X", NULL};
+
+ if (!PyArg_ParseTupleAndKeywords(args, kwrds, "O", kwlist, &X))
+ return NULL;
+
+ if (!(Matrix_Check(X) && MAT_ID(X) == COMPLEX))
+ PY_ERR(PyExc_ValueError, "X must be a dense matrix with type 'z'");
+
+ int m = X->nrows, n = X->ncols;
+ if (m == 0 || n == 0) return Py_BuildValue("");
+
+ fftw_plan p = fftw_plan_dft_2d(n, m,
+ X->buffer, X->buffer, FFTW_FORWARD, FFTW_ESTIMATE);
+
+ fftw_execute(p);
+ fftw_destroy_plan(p);
+ return Py_BuildValue("");
+}
+
static char doc_idft[] =
"IDFT of a matrix, X := idft(X)\n\n"
"PURPOSE\n"
@@ -96,6 +124,41 @@ static PyObject *idft(PyObject *self, PyObject *args, PyObject *kwrds)
return Py_BuildValue("");
}
+static char doc_idft2[] =
+ "2D-IDFT of a matrix, X := idft2(X)\n\n"
+ "PURPOSE\n"
+ "Computes the inverse 2D-DFT of a dense matrix X.\n\n"
+ "ARGUMENTS\n"
+ "X A dense matrix of typecode 'z'.";
+
+static PyObject *idft2(PyObject *self, PyObject *args, PyObject *kwrds)
+{
+ matrix *X;
+ char *kwlist[] = {"X", NULL};
+
+ if (!PyArg_ParseTupleAndKeywords(args, kwrds, "O", kwlist, &X))
+ return NULL;
+
+ if (!(Matrix_Check(X) && MAT_ID(X) == COMPLEX))
+ PY_ERR(PyExc_ValueError, "X must be a dense matrix with type 'z'");
+
+ int m = X->nrows, n = X->ncols;
+ if (m == 0 || n == 0) return Py_BuildValue("");
+
+ fftw_plan p = fftw_plan_dft_2d(n, m,
+ X->buffer, X->buffer, FFTW_BACKWARD, FFTW_ESTIMATE);
+
+ fftw_execute(p);
+
+ number a;
+ a.z = 1.0/(m*n);
+ int mn = m*n, ix = 1;
+ zscal_(&mn, &a.z, MAT_BUFZ(X), &ix);
+
+ fftw_destroy_plan(p);
+ return Py_BuildValue("");
+}
+
static char doc_dct[] =
"DCT of a matrix.\n"
"X := dct(X, type=2)\n\n"
@@ -142,6 +205,52 @@ static PyObject *dct(PyObject *self, PyObject *args, PyObject *kwrds)
return Py_BuildValue("");
}
+static char doc_dct2[] =
+ "2D-DCT of a matrix.\n"
+ "X := dct2(X, type=2)\n\n"
+ "PURPOSE\n"
+ "Computes the 2D-DCT of a dense matrix X.\n\n"
+ "ARGUMENTS\n"
+ "X A dense matrix of typecode 'd'.\n\n"
+ "type integer from 1 to 4; chooses either DCT-I, DCT-II, \n"
+ " DCT-III or DCT-IV.";
+
+
+static PyObject *dct2(PyObject *self, PyObject *args, PyObject *kwrds)
+{
+ matrix *X;
+ int type = 2;
+ char *kwlist[] = {"X", "type", NULL};
+
+ if (!PyArg_ParseTupleAndKeywords(args, kwrds, "O|i", kwlist, &X, &type))
+ return NULL;
+
+ if (!(Matrix_Check(X) && MAT_ID(X) == DOUBLE))
+ PY_ERR(PyExc_ValueError, "X must be a dense matrix with type 'd'");
+
+ int m = X->nrows, n = X->ncols;
+ if (m == 0 || n == 0) return Py_BuildValue("");
+
+ fftw_r2r_kind kind;
+ switch(type) {
+ case 1:
+ kind = FFTW_REDFT00;
+ if (m <= 1) PY_ERR(PyExc_ValueError, "m must be greater than 1 for DCT-I");
+ break;
+ case 2: kind = FFTW_REDFT10; break;
+ case 3: kind = FFTW_REDFT01; break;
+ case 4: kind = FFTW_REDFT11; break;
+ default: PY_ERR(PyExc_ValueError, "type must be between 1 and 4");
+ }
+
+ fftw_plan p = fftw_plan_r2r_2d(n, m,
+ X->buffer, X->buffer, kind, kind, FFTW_ESTIMATE);
+
+ fftw_execute(p);
+ fftw_destroy_plan(p);
+ return Py_BuildValue("");
+}
+
static char doc_idct[] =
"IDCT of a matrix.\n"
"X := idct(X, type=2)\n\n"
@@ -192,6 +301,56 @@ static PyObject *idct(PyObject *self, PyObject *args, PyObject *kwrds)
return Py_BuildValue("");
}
+static char doc_idct2[] =
+ "2D-IDCT of a matrix.\n"
+ "X := idct2(X, type=2)\n\n"
+ "PURPOSE\n"
+ "Computes the 2D-IDCT of a dense matrix X.\n\n"
+ "ARGUMENTS\n"
+ "X A dense matrix of typecode 'd'.\n\n"
+ "type integer from 1 to 4; chooses the inverse transform for\n"
+ " either DCT-I, DCT-II, DCT-III or DCT-IV.";
+
+static PyObject *idct2(PyObject *self, PyObject *args, PyObject *kwrds)
+{
+ matrix *X;
+ int type = 2;
+ char *kwlist[] = {"X", "type", NULL};
+
+ if (!PyArg_ParseTupleAndKeywords(args, kwrds, "O|i", kwlist, &X, &type))
+ return NULL;
+
+ if (!(Matrix_Check(X) && MAT_ID(X) == DOUBLE))
+ PY_ERR(PyExc_ValueError, "X must be a dense matrix with type 'd'");
+
+ int m = X->nrows, n = X->ncols;
+ if (m == 0 || n == 0) return Py_BuildValue("");
+
+ fftw_r2r_kind kind;
+ switch(type) {
+ case 1:
+ kind = FFTW_REDFT00;
+ if (MIN(m,n) <= 1)
+ PY_ERR(PyExc_ValueError, "m and n must both greater than 1 for DCT2-I");
+ break;
+ case 2: kind = FFTW_REDFT01; break;
+ case 3: kind = FFTW_REDFT10; break;
+ case 4: kind = FFTW_REDFT11; break;
+ default: PY_ERR(PyExc_ValueError, "type must be between 1 and 4");
+ }
+
+ fftw_plan p = fftw_plan_r2r_2d(n, m,
+ X->buffer, X->buffer, kind, kind, FFTW_ESTIMATE);
+
+ double a = 1.0/(type == 1 ? MAX(1,2*(m-1))*MAX(1,2*(n-1)) : 4*m*n);
+ int mn = m*n, ix = 1;
+ dscal_(&mn, &a, MAT_BUFD(X), &ix);
+
+ fftw_execute(p);
+ fftw_destroy_plan(p);
+ return Py_BuildValue("");
+}
+
static char doc_dst[] =
"DST of a matrix.\n"
"X := dst(X, type=1)\n\n"
@@ -235,6 +394,48 @@ static PyObject *dst(PyObject *self, PyObject *args, PyObject *kwrds)
return Py_BuildValue("");
}
+static char doc_dst2[] =
+ "2D-DST of a matrix.\n"
+ "X := dst2(X, type=1)\n\n"
+ "PURPOSE\n"
+ "Computes the 2D-DST of a dense matrix X.\n\n"
+ "ARGUMENTS\n"
+ "X A dense matrix of typecode 'd'.\n\n"
+ "type integer from 1 to 4; chooses either DST-I, DST-II, \n"
+ " DST-III or DST-IV.";
+
+static PyObject *dst2(PyObject *self, PyObject *args, PyObject *kwrds)
+{
+ matrix *X;
+ int type = 1;
+ char *kwlist[] = {"X", "type", NULL};
+
+ if (!PyArg_ParseTupleAndKeywords(args, kwrds, "O|i", kwlist, &X, &type))
+ return NULL;
+
+ if (!(Matrix_Check(X) && MAT_ID(X) == DOUBLE))
+ PY_ERR(PyExc_ValueError, "X must be a dense matrix with type 'd'");
+
+ int m = X->nrows, n = X->ncols;
+ if (m == 0 || n == 0) return Py_BuildValue("");
+
+ fftw_r2r_kind kind;
+ switch(type) {
+ case 1: kind = FFTW_RODFT00; break;
+ case 2: kind = FFTW_RODFT10; break;
+ case 3: kind = FFTW_RODFT01; break;
+ case 4: kind = FFTW_RODFT11; break;
+ default: PY_ERR(PyExc_ValueError, "type must be between 1 and 4");
+ }
+
+ fftw_plan p = fftw_plan_r2r_2d(n, m,
+ X->buffer, X->buffer, kind, kind, FFTW_ESTIMATE);
+
+ fftw_execute(p);
+ fftw_destroy_plan(p);
+ return Py_BuildValue("");
+}
+
static char doc_idst[] =
"IDST of a matrix.\n"
"X := idst(X, type=1)\n\n"
@@ -283,13 +484,65 @@ static PyObject *idst(PyObject *self, PyObject *args, PyObject *kwrds)
return Py_BuildValue("");
}
+static char doc_idst2[] =
+ "2D-IDST of a matrix.\n"
+ "X := idst2(X, type=1)\n\n"
+ "PURPOSE\n"
+ "Computes the 2D-IDST of a dense matrix X.\n\n"
+ "ARGUMENTS\n"
+ "X A dense matrix of typecode 'd'.\n\n"
+ "type integer from 1 to 4; chooses the inverse transform for\n"
+ " either DST-I, DST-II, DST-III or DST-IV.";
+
+static PyObject *idst2(PyObject *self, PyObject *args, PyObject *kwrds)
+{
+ matrix *X;
+ int type = 1;
+ char *kwlist[] = {"X", "type", NULL};
+
+ if (!PyArg_ParseTupleAndKeywords(args, kwrds, "O|i", kwlist, &X, &type))
+ return NULL;
+
+ if (!(Matrix_Check(X) && MAT_ID(X) == DOUBLE))
+ PY_ERR(PyExc_ValueError, "X must be a dense matrix with type 'd'");
+
+ int m = X->nrows, n = X->ncols;
+ if (m == 0 || n == 0) return Py_BuildValue("");
+
+ fftw_r2r_kind kind;
+ switch(type) {
+ case 1: kind = FFTW_RODFT00; break;
+ case 2: kind = FFTW_RODFT01; break;
+ case 3: kind = FFTW_RODFT10; break;
+ case 4: kind = FFTW_RODFT11; break;
+ default: PY_ERR(PyExc_ValueError, "type must be between 1 and 4");
+ }
+ fftw_plan p = fftw_plan_r2r_2d(n, m,
+ X->buffer, X->buffer, kind, kind, FFTW_ESTIMATE);
+
+ fftw_execute(p);
+
+ double a = 1.0/(type == 1 ? MAX(1,2*(m+1))*MAX(1,2*(n+1)) : 4*m*n);
+ int mn = m*n, ix = 1;
+ dscal_(&mn, &a, MAT_BUFD(X), &ix);
+
+ fftw_destroy_plan(p);
+ return Py_BuildValue("");
+}
+
static PyMethodDef fftw_functions[] = {
{"dft", (PyCFunction) dft, METH_VARARGS|METH_KEYWORDS, doc_dft},
+ {"dft2", (PyCFunction) dft2, METH_VARARGS|METH_KEYWORDS, doc_dft2},
{"idft", (PyCFunction) idft, METH_VARARGS|METH_KEYWORDS, doc_idft},
+ {"idft2", (PyCFunction) idft2, METH_VARARGS|METH_KEYWORDS, doc_idft2},
{"dct", (PyCFunction) dct, METH_VARARGS|METH_KEYWORDS, doc_dct},
+ {"dct2", (PyCFunction) dct2, METH_VARARGS|METH_KEYWORDS, doc_dct2},
{"idct", (PyCFunction) idct, METH_VARARGS|METH_KEYWORDS, doc_idct},
+ {"idct2", (PyCFunction) idct2, METH_VARARGS|METH_KEYWORDS, doc_idct2},
{"dst", (PyCFunction) dst, METH_VARARGS|METH_KEYWORDS, doc_dst},
+ {"dst2", (PyCFunction) dst2, METH_VARARGS|METH_KEYWORDS, doc_dst2},
{"idst", (PyCFunction) idst, METH_VARARGS|METH_KEYWORDS, doc_idst},
+ {"idst2", (PyCFunction) idst2, METH_VARARGS|METH_KEYWORDS, doc_idst2},
{NULL} /* Sentinel */
};
diff --git a/src/C/glpk.c b/src/C/glpk.c
index b052c40..cd4bc9a 100644
--- a/src/C/glpk.c
+++ b/src/C/glpk.c
@@ -1,7 +1,7 @@
/*
- * Copyright 2004-2008 J. Dahl and L. Vandenberghe.
+ * Copyright 2004-2007 J. Dahl and L. Vandenberghe.
*
- * This file is part of CVXOPT version 0.9.3.
+ * This file is part of CVXOPT version 1.0.
*
* CVXOPT is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
@@ -18,16 +18,9 @@
*/
#include "cvxopt.h"
+#include "misc.h"
#include "glpk.h"
-#define err_dbl_mtrx(s) {PyErr_SetString(PyExc_TypeError, \
- s " must be a matrix with typecode 'd'"); \
- return NULL;}
-
-#define err_p_int(s) {PyErr_SetString(PyExc_ValueError, \
- s " must be a positive integer"); \
- return NULL;}
-
PyDoc_STRVAR(glpk__doc__,
"Interface to the simplex algorithm in GLPK.\n\n"
"The GLPK control parameters have the default values listed in \n"
@@ -96,15 +89,15 @@ static int get_param_idx(char *str, int *idx, char *type)
static char doc_simplex[] =
"Solves a linear program using GLPK.\n\n"
- "(status, x, z, y) = solvelp(c, G, h, A, b)\n"
- "(status, x, z) = solvelp(c, G, h)\n\n"
+ "(status, x, z, y) = lp(c, G, h, A, b)\n"
+ "(status, x, z) = lp(c, G, h)\n\n"
"PURPOSE\n"
- "(status, x, z, y) = solvelp(c, G, h, A, b) solves the pair\n"
+ "(status, x, z, y) = lp(c, G, h, A, b) solves the pair\n"
"of primal and dual LPs\n\n"
" minimize c'*x maximize -h'*z + b'*y\n"
" subject to G*x <= h subject to G'*z + A'*y + c = 0\n"
" A*x = b z >= 0.\n\n"
- "(status, x, z) = solvelp(c, G, h) solves the pair of primal\n"
+ "(status, x, z) = lp(c, G, h) solves the pair of primal\n"
"and dual LPs\n\n"
" minimize c'*x maximize -h'*z \n"
" subject to G*x <= h subject to G'*z + c = 0\n"
@@ -130,7 +123,7 @@ static PyObject *simplex(PyObject *self, PyObject *args,
PyObject *G, *A=NULL, *t=NULL, *param, *key, *value;
LPX *lp;
int m, n, p, i, j, k, nnz, nnzmax, *rn=NULL, *cn=NULL, param_id;
- int_compat pos=0;
+ int_t pos=0;
double *a=NULL, val;
char param_type, err_str[100], *keystr;
char *kwlist[] = {"c", "G", "h", "A", "b", NULL};
@@ -354,9 +347,319 @@ static PyObject *simplex(PyObject *self, PyObject *args,
}
+
+static char doc_integer[] =
+ "Solves a mixed integer linear program using GLPK.\n\n"
+ "(status, x) = ilp(c, G, h, A, b, I, B)\n\n"
+ "PURPOSE\n"
+ "Solves the mixed integer linear programming problem\n\n"
+ " minimize c'*x\n"
+ " subject to G*x <= h\n"
+ " A*x = b\n"
+ " x[I] are all integer\n"
+ " x[B] are all binary\n\n"
+ "ARGUMENTS\n"
+ "c nx1 dense 'd' matrix with n>=1\n\n"
+ "G mxn dense or sparse 'd' matrix with m>=1\n\n"
+ "h mx1 dense 'd' matrix\n\n"
+ "A pxn dense or sparse 'd' matrix with p>=0\n\n"
+ "b px1 dense 'd' matrix\n\n"
+ "I set with indices of integer variables\n\n"
+ "B set with indices of binary variables\n\n"
+ "status 'optimal', 'primal infeasible', 'dual infeasible', \n"
+ " 'invalid MIP formulation', 'maxiters exceeded', \n"
+ " 'time limit exceeded', 'unknown'\n\n"
+ "x an optimal solution if status is 'optimal';\n"
+ " None otherwise";
+
+static PyObject *integer(PyObject *self, PyObject *args,
+ PyObject *kwrds)
+{
+ matrix *c, *h, *b=NULL, *x=NULL;
+ PyObject *G, *A=NULL, *IntSet=NULL, *BinSet = NULL;
+ PyObject *t=NULL, *param, *key, *value;
+ LPX *lp;
+ int m, n, p, i, j, k, nnz, nnzmax, *rn=NULL, *cn=NULL, param_id;
+ int_t pos=0;
+ double *a=NULL, val;
+ char param_type, err_str[100], *keystr;
+ char *kwlist[] = {"c", "G", "h", "A", "b", "I", "B", NULL};
+
+ if (!PyArg_ParseTupleAndKeywords(args, kwrds, "OOO|OOOO", kwlist, &c,
+ &G, &h, &A, &b, &IntSet, &BinSet)) return NULL;
+
+ if ((Matrix_Check(G) && MAT_ID(G) != DOUBLE) ||
+ (SpMatrix_Check(G) && SP_ID(G) != DOUBLE) ||
+ (!Matrix_Check(G) && !SpMatrix_Check(G))){
+ PyErr_SetString(PyExc_TypeError, "G must be a 'd' matrix");
+ return NULL;
+ }
+ if ((m = Matrix_Check(G) ? MAT_NROWS(G) : SP_NROWS(G)) <= 0)
+ err_p_int("m");
+ if ((n = Matrix_Check(G) ? MAT_NCOLS(G) : SP_NCOLS(G)) <= 0)
+ err_p_int("n");
+
+ if (!Matrix_Check(h) || h->id != DOUBLE) err_dbl_mtrx("h");
+ if (h->nrows != m || h->ncols != 1){
+ PyErr_SetString(PyExc_ValueError, "incompatible dimensions");
+ return NULL;
+ }
+
+ if (A){
+ if ((Matrix_Check(A) && MAT_ID(A) != DOUBLE) ||
+ (SpMatrix_Check(A) && SP_ID(A) != DOUBLE) ||
+ (!Matrix_Check(A) && !SpMatrix_Check(A))){
+ PyErr_SetString(PyExc_ValueError, "A must be a dense "
+ "'d' matrix or a general sparse matrix");
+ return NULL;
+ }
+ if ((p = Matrix_Check(A) ? MAT_NROWS(A) : SP_NROWS(A)) < 0)
+ err_p_int("p");
+ if ((Matrix_Check(A) ? MAT_NCOLS(A) : SP_NCOLS(A)) != n){
+ PyErr_SetString(PyExc_ValueError, "incompatible "
+ "dimensions");
+ return NULL;
+ }
+ }
+ else p = 0;
+
+ if (b && (!Matrix_Check(b) || b->id != DOUBLE)) err_dbl_mtrx("b");
+ if ((b && (b->nrows != p || b->ncols != 1)) || (!b && p !=0 )){
+ PyErr_SetString(PyExc_ValueError, "incompatible dimensions");
+ return NULL;
+ }
+
+ if ((IntSet) && (!PyAnySet_Check(IntSet)))
+ PY_ERR_TYPE("invalid integer index set");
+
+ if ((BinSet) && (!PyAnySet_Check(BinSet)))
+ PY_ERR_TYPE("invalid binary index set");
+
+ lp = lpx_create_prob();
+ lpx_add_rows(lp, m+p);
+ lpx_add_cols(lp, n);
+
+ for (i=0; i<n; i++){
+ lpx_set_obj_coef(lp, i+1, MAT_BUFD(c)[i]);
+ lpx_set_col_bnds(lp, i+1, LPX_FR, 0.0, 0.0);
+ }
+ for (i=0; i<m; i++)
+ lpx_set_row_bnds(lp, i+1, LPX_UP, 0.0, MAT_BUFD(h)[i]);
+ for (i=0; i<p; i++)
+ lpx_set_row_bnds(lp, i+m+1, LPX_FX, MAT_BUFD(b)[i],
+ MAT_BUFD(b)[i]);
+
+ nnzmax = (SpMatrix_Check(G) ? SP_NNZ(G) : m*n ) +
+ ((A && SpMatrix_Check(A)) ? SP_NNZ(A) : p*n);
+ a = (double *) calloc(nnzmax+1, sizeof(double));
+ rn = (int *) calloc(nnzmax+1, sizeof(int));
+ cn = (int *) calloc(nnzmax+1, sizeof(int));
+ if (!a || !rn || !cn){
+ free(a); free(rn); free(cn); lpx_delete_prob(lp);
+ return PyErr_NoMemory();
+ }
+
+ nnz = 0;
+ if (SpMatrix_Check(G)) {
+ for (j=0; j<n; j++) for (k=SP_COL(G)[j]; k<SP_COL(G)[j+1]; k++)
+ if ((val = SP_VALD(G)[k]) != 0.0){
+ a[1+nnz] = val;
+ rn[1+nnz] = SP_ROW(G)[k]+1;
+ cn[1+nnz] = j+1;
+ nnz++;
+ }
+ }
+ else for (j=0; j<n; j++) for (i=0; i<m; i++)
+ if ((val = MAT_BUFD(G)[i+j*m]) != 0.0){
+ a[1+nnz] = val;
+ rn[1+nnz] = i+1;
+ cn[1+nnz] = j+1;
+ nnz++;
+ }
+
+ if (A && SpMatrix_Check(A)){
+ for (j=0; j<n; j++) for (k=SP_COL(A)[j]; k<SP_COL(A)[j+1]; k++)
+ if ((val = SP_VALD(A)[k]) != 0.0){
+ a[1+nnz] = val;
+ rn[1+nnz] = m+SP_ROW(A)[k]+1;
+ cn[1+nnz] = j+1;
+ nnz++;
+ }
+ }
+ else for (j=0; j<n; j++) for (i=0; i<p; i++)
+ if ((val = MAT_BUFD(A)[i+j*p]) != 0.0){
+ a[1+nnz] = val;
+ rn[1+nnz] = m+i+1;
+ cn[1+nnz] = j+1;
+ nnz++;
+ }
+
+ lpx_load_matrix(lp, nnz, rn, cn, a);
+ free(rn); free(cn); free(a);
+
+ if (!(t = PyTuple_New(2))) {
+ lpx_delete_prob(lp);
+ return PyErr_NoMemory();
+ }
+
+ if (!(param = PyObject_GetAttrString(glpk_module, "options"))
+ || !PyDict_Check(param)){
+ lpx_delete_prob(lp);
+ PyErr_SetString(PyExc_AttributeError,
+ "missing glpk.options dictionary");
+ return NULL;
+ }
+
+ while (PyDict_Next(param, &pos, &key, &value))
+ if ((keystr = PyString_AsString(key)) && get_param_idx(keystr,
+ ¶m_id, ¶m_type)){
+ if (param_type == 'i'){
+ if (!PyInt_Check(value)){
+ sprintf(err_str, "invalid value for integer "
+ "GLPK parameter: %-.20s", keystr);
+ PyErr_SetString(PyExc_ValueError, err_str);
+ lpx_delete_prob(lp);
+ Py_DECREF(param);
+ return NULL;
+ }
+ if (!strcmp("LPX_K_PRESOL", keystr) &&
+ PyInt_AS_LONG(value) != 1){
+ PyErr_Warn(PyExc_UserWarning, "ignoring value of "
+ "GLPK parameter 'LPX_K_PRESOL'");
+ }
+ else lpx_set_int_parm(lp, param_id,
+ PyInt_AS_LONG(value));
+ }
+ else {
+ if (!PyInt_Check(value) && !PyFloat_Check(value)){
+ sprintf(err_str, "invalid value for floating point "
+ "GLPK parameter: %-.20s", keystr);
+ PyErr_SetString(PyExc_ValueError, err_str);
+ lpx_delete_prob(lp);
+ Py_DECREF(param);
+ return NULL;
+ }
+ lpx_set_real_parm(lp, param_id,
+ PyFloat_AsDouble(value));
+ }
+ }
+ lpx_set_int_parm(lp, LPX_K_PRESOL, 1);
+ Py_DECREF(param);
+
+ if (IntSet) {
+ PyObject *iter = PySequence_Fast(IntSet, "Critical error: not sequence");
+
+ for (i=0; i<PySet_GET_SIZE(IntSet); i++) {
+
+ PyObject *tmp = PySequence_Fast_GET_ITEM(iter, i);
+ if (!PyInt_Check(tmp)) {
+ lpx_delete_prob(lp);
+ Py_DECREF(iter);
+ PY_ERR_TYPE("non-integer element in I");
+ }
+ int k = PyInt_AS_LONG(tmp);
+ if ((k < 0) || (k >= n)) {
+ lpx_delete_prob(lp);
+ Py_DECREF(iter);
+ PY_ERR(PyExc_IndexError, "index element out of range in I");
+ }
+ glp_set_col_kind(lp, k+1, GLP_IV);
+ }
+
+ Py_DECREF(iter);
+ }
+
+ if (BinSet) {
+ PyObject *iter = PySequence_Fast(BinSet, "Critical error: not sequence");
+
+ for (i=0; i<PySet_GET_SIZE(BinSet); i++) {
+
+ PyObject *tmp = PySequence_Fast_GET_ITEM(iter, i);
+ if (!PyInt_Check(tmp)) {
+ lpx_delete_prob(lp);
+ Py_DECREF(iter);
+ PY_ERR_TYPE("non-binary element in I");
+ }
+ int k = PyInt_AS_LONG(tmp);
+ if ((k < 0) || (k >= n)) {
+ lpx_delete_prob(lp);
+ Py_DECREF(iter);
+ PY_ERR(PyExc_IndexError, "index element out of range in B");
+ }
+ glp_set_col_kind(lp, k+1, GLP_BV);
+ }
+
+ Py_DECREF(iter);
+ }
+
+
+
+ switch (lpx_intopt(lp)){
+
+ case LPX_E_OK:
+
+ x = (matrix *) Matrix_New(n,1,DOUBLE);
+ if (!x) {
+ Py_XDECREF(t);
+ lpx_delete_prob(lp);
+ return PyErr_NoMemory();
+ }
+ PyTuple_SET_ITEM(t, 0, (PyObject *)
+ PyString_FromString("optimal"));
+
+ for (i=0; i<n; i++)
+ MAT_BUFD(x)[i] = lpx_mip_col_val(lp, i+1);
+ PyTuple_SET_ITEM(t, 1, (PyObject *) x);
+
+ lpx_delete_prob(lp);
+ return (PyObject *) t;
+
+ case LPX_E_FAULT:
+ PyTuple_SET_ITEM(t, 0, (PyObject *)
+ PyString_FromString("invalid MIP formulation"));
+
+ case LPX_E_NOPFS:
+ PyTuple_SET_ITEM(t, 0, (PyObject *)
+ PyString_FromString("primal infeasible"));
+
+ case LPX_E_NODFS:
+
+ PyTuple_SET_ITEM(t, 0, (PyObject *)
+ PyString_FromString("dual infeasible"));
+
+ case LPX_E_ITLIM:
+
+ PyTuple_SET_ITEM(t, 0, (PyObject *)
+ PyString_FromString("maxiters exceeded"));
+
+ case LPX_E_TMLIM:
+
+ PyTuple_SET_ITEM(t, 0, (PyObject *)
+ PyString_FromString("time limit exceeded"));
+
+ case LPX_E_SING:
+
+ PyTuple_SET_ITEM(t, 0, (PyObject *)
+ PyString_FromString("singular or ill-conditioned basis"));
+
+ default:
+
+ PyTuple_SET_ITEM(t, 0, (PyObject *)
+ PyString_FromString("unknown"));
+ }
+
+ lpx_delete_prob(lp);
+
+ PyTuple_SET_ITEM(t, 1, Py_BuildValue(""));
+ return (PyObject *) t;
+}
+
+
static PyMethodDef glpk_functions[] = {
- {"solvelp", (PyCFunction) simplex, METH_VARARGS|METH_KEYWORDS,
+ {"lp", (PyCFunction) simplex, METH_VARARGS|METH_KEYWORDS,
doc_simplex},
+ {"ilp", (PyCFunction) integer, METH_VARARGS|METH_KEYWORDS,
+ doc_integer},
{NULL} /* Sentinel */
};
diff --git a/src/C/gsl.c b/src/C/gsl.c
index c1e9edb..19e2c56 100644
--- a/src/C/gsl.c
+++ b/src/C/gsl.c
@@ -1,7 +1,7 @@
/*
* Copyright 2004-2008 J. Dahl and L. Vandenberghe.
*
- * This file is part of CVXOPT version 0.9.3.
+ * This file is part of CVXOPT version 1.0.
*
* CVXOPT is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
diff --git a/src/C/lapack.c b/src/C/lapack.c
index 2025695..e5bbe70 100644
--- a/src/C/lapack.c
+++ b/src/C/lapack.c
@@ -1,7 +1,7 @@
/*
* Copyright 2004-2008 J. Dahl and L. Vandenberghe.
*
- * This file is part of CVXOPT version 0.9.3.
+ * This file is part of CVXOPT version 1.0.
*
* CVXOPT is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
diff --git a/src/C/misc.h b/src/C/misc.h
index a0e7501..debb469 100644
--- a/src/C/misc.h
+++ b/src/C/misc.h
@@ -1,7 +1,7 @@
/*
* Copyright 2004-2008 J. Dahl and L. Vandenberghe.
*
- * This file is part of CVXOPT version 0.9.3.
+ * This file is part of CVXOPT version 1.0.
*
* CVXOPT is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
diff --git a/src/C/misc_solvers.c b/src/C/misc_solvers.c
new file mode 100644
index 0000000..14f7b44
--- /dev/null
+++ b/src/C/misc_solvers.c
@@ -0,0 +1,1004 @@
+/*
+ * Copyright 2004-2008 J. Dahl and L. Vandenberghe.
+ *
+ * This file is part of CVXOPT version 1.0.
+ *
+ * CVXOPT 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.
+ *
+ * CVXOPT 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 this program. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include "Python.h"
+#include "cvxopt.h"
+#include "misc.h"
+#include "math.h"
+#include "float.h"
+
+PyDoc_STRVAR(misc_solvers__doc__, "Miscellaneous functions used by the "
+ "CVXOPT solvers.");
+
+extern void dcopy_(int *n, double *x, int *incx, double *y, int *incy);
+extern double dnrm2_(int *n, double *x, int *incx);
+extern double ddot_(int *n, double *x, int *incx, double *y, int *incy);
+extern void dscal_(int *n, double *alpha, double *x, int *incx);
+extern void daxpy_(int *n, double *alpha, double *x, int *incx, double *y,
+ int *incy);
+extern void dtbmv_(char *uplo, char *trans, char *diag, int *n, int *k,
+ double *A, int *lda, double *x, int *incx);
+extern void dtbsv_(char *uplo, char *trans, char *diag, int *n, int *k,
+ double *A, int *lda, double *x, int *incx);
+extern void dgemv_(char* trans, int *m, int *n, double *alpha, double *A,
+ int *lda, double *x, int *incx, double *beta, double *y, int *incy);
+extern void dger_(int *m, int *n, double *alpha, double *x, int *incx,
+ double *y, int *incy, double *A, int *lda);
+extern void dtrmm_(char *side, char *uplo, char *transa, char *diag,
+ int *m, int *n, double *alpha, double *A, int *lda, double *B,
+ int *ldb);
+extern void dsyr2k_(char *uplo, char *trans, int *n, int *k, double *alpha,
+ double *A, int *lda, double *B, int *ldb, double *beta, double *C,
+ int *ldc);
+extern void dlacpy_(char *uplo, int *m, int *n, double *A, int *lda,
+ double *B, int *ldb);
+extern void dsyevr_(char *jobz, char *range, char *uplo, int *n, double *A,
+ int *ldA, double *vl, double *vu, int *il, int *iu, double *abstol,
+ int *m, double *W, double *Z, int *ldZ, int *isuppz, double *work,
+ int *lwork, int *iwork, int *liwork, int *info);
+extern void dsyevd_(char *jobz, char *uplo, int *n, double *A, int *ldA,
+ double *W, double *work, int *lwork, int *iwork, int *liwork,
+ int *info);
+
+
+static char doc_scale[] =
+ "Applies Nesterov-Todd scaling or its inverse.\n\n"
+ "scale(x, W, trans = 'N', inverse = 'N')\n\n"
+ "Computes\n\n"
+ " x := W*x (trans is 'N', inverse = 'N')\n"
+ " x := W^T*x (trans is 'T', inverse = 'N')\n"
+ " x := W^{-1}*x (trans is 'N', inverse = 'I')\n"
+ " x := W^{-T}*x (trans is 'T', inverse = 'I').\n\n"
+ "x is a dense 'd' matrix.\n\n"
+ "W is a dictionary with entries:\n\n"
+ "- W['dnl']: positive vector\n"
+ "- W['dnli']: componentwise inverse of W['dnl']\n"
+ "- W['d']: positive vector\n"
+ "- W['di']: componentwise inverse of W['d']\n"
+ "- W['v']: lists of 2nd order cone vectors with unit hyperbolic \n"
+ " norms\n"
+ "- W['beta']: list of positive numbers\n"
+ "- W['r']: list of square matrices\n"
+ "- W['rti']: list of square matrices. rti[k] is the inverse\n"
+ " transpose of r[k]. \n\n"
+ "The 'dnl' and 'dnli' entries are optional, and only present when \n"
+ "the function is called from the nonlinear solver.";
+
+static PyObject* scale(PyObject *self, PyObject *args, PyObject *kwrds)
+{
+ matrix *x, *d, *vk, *rk;
+ PyObject *W, *v, *beta, *r, *betak;
+ char trans = 'N', inverse = 'N';
+ int m, n, ind = 0, int0 = 0, int1 = 1, i, k, inc, len, ld, maxn, N;
+ double b, dbl0 = 0.0, dbl1 = 1.0, dblm1 = -1.0, dbl2 = 2.0, dbl5 = 0.5,
+ *wrk;
+ char *kwlist[] = {"x", "W", "trans", "inverse", NULL};
+
+ if (!PyArg_ParseTupleAndKeywords(args, kwrds, "OO|cc", kwlist,
+ &x, &W, &trans, &inverse)) return NULL;
+
+
+ /*
+ * Scaling for nonlinear component xk is xk := dnl .* xk; inverse is
+ * xk ./ dnl = dnli .* xk, where dnl = W['dnl'], dnli = W['dnli'].
+ */
+
+ if ((d = (inverse == 'N') ? (matrix *) PyDict_GetItemString(W, "dnl") :
+ (matrix *) PyDict_GetItemString(W, "dnli"))){
+ m = len(d);
+ for (i = 0; i < x->ncols; i++)
+ dtbmv_("L", "N", "N", &m, &int0, MAT_BUFD(d), &int1,
+ MAT_BUFD(x) + i*x->nrows, &int1);
+ ind += m;
+ }
+
+
+ /*
+ * Scaling for 'l' component xk is xk := d .* xk; inverse scaling is
+ * xk ./ d = di .* xk, where d = W['d'], di = W['di'].
+ */
+
+ if (!(d = (inverse == 'N') ? (matrix *) PyDict_GetItemString(W, "d") :
+ (matrix *) PyDict_GetItemString(W, "di"))){
+ PyErr_SetString(PyExc_KeyError, "missing item W['d'] or W['di']");
+ return NULL;
+ }
+ m = len(d);
+ for (i = 0; i < x->ncols; i++)
+ dtbmv_("L", "N", "N", &m, &int0, MAT_BUFD(d), &int1, MAT_BUFD(x)
+ + i*x->nrows + ind, &int1);
+ ind += m;
+
+
+ /*
+ * Scaling for 'q' component is
+ *
+ * xk := beta * (2*v*v' - J) * xk
+ * = beta * (2*v*(xk'*v)' - J*xk)
+ *
+ * where beta = W['beta'][k], v = W['v'][k], J = [1, 0; 0, -I].
+ *
+ * Inverse scaling is
+ *
+ * xk := 1/beta * (2*J*v*v'*J - J) * xk
+ * = 1/beta * (-J) * (2*v*((-J*xk)'*v)' + xk).
+ */
+
+ v = PyDict_GetItemString(W, "v");
+ beta = PyDict_GetItemString(W, "beta");
+ N = (int) PyList_Size(v);
+ if (!(wrk = (double *) calloc(x->ncols, sizeof(double))))
+ return PyErr_NoMemory();
+ for (k = 0; k < N; k++){
+ vk = (matrix *) PyList_GetItem(v, (Py_ssize_t) k);
+ m = vk->nrows;
+ if (inverse == 'I')
+ dscal_(&(x->ncols), &dblm1, MAT_BUFD(x) + ind, &(x->nrows));
+ ld = MAX(x->nrows, 1);
+ dgemv_("T", &m, &(x->ncols), &dbl1, MAT_BUFD(x) + ind, &ld,
+ MAT_BUFD(vk), &int1, &dbl0, wrk, &int1);
+ dscal_(&(x->ncols), &dblm1, MAT_BUFD(x) + ind, &(x->nrows));
+ dger_(&m, &(x->ncols), &dbl2, MAT_BUFD(vk), &int1, wrk, &int1,
+ MAT_BUFD(x) + ind, &ld);
+ if (inverse == 'I')
+ dscal_(&(x->ncols), &dblm1, MAT_BUFD(x) + ind, &(x->nrows));
+
+ betak = PyList_GetItem(beta, (Py_ssize_t) k);
+ b = PyFloat_AS_DOUBLE(betak);
+ if (inverse == 'I') b = 1.0 / b;
+ for (i = 0; i < x->ncols; i++)
+ dscal_(&m, &b, MAT_BUFD(x) + ind + i*x->nrows, &int1);
+ ind += m;
+ }
+ free(wrk);
+
+
+ /*
+ * Scaling for 's' component xk is
+ *
+ * xk := vec( r' * mat(xk) * r ) if trans = 'N'
+ * xk := vec( r * mat(xk) * r' ) if trans = 'T'.
+ *
+ * r is kth element of W['r'].
+ *
+ * Inverse scaling is
+ *
+ * xk := vec( rti * mat(xk) * rti' ) if trans = 'N'
+ * xk := vec( rti' * mat(xk) * rti ) if trans = 'T'.
+ *
+ * rti is kth element of W['rti'].
+ */
+
+ r = (inverse == 'N') ? PyDict_GetItemString(W, "r") :
+ PyDict_GetItemString(W, "rti");
+ N = (int) PyList_Size(r);
+ for (k = 0, maxn = 0; k < N; k++){
+ rk = (matrix *) PyList_GetItem(r, (Py_ssize_t) k);
+ maxn = MAX(maxn, rk->nrows);
+ }
+ if (!(wrk = (double *) calloc(maxn*maxn, sizeof(double))))
+ return PyErr_NoMemory();
+ for (k = 0; k < N; k++){
+ rk = (matrix *) PyList_GetItem(r, (Py_ssize_t) k);
+ n = rk->nrows;
+ for (i = 0; i < x->ncols; i++){
+
+ /* scale diagonal of rk by 0.5 */
+ inc = n + 1;
+ dscal_(&n, &dbl5, MAT_BUFD(x) + ind + i*x->nrows, &inc);
+
+ /* wrk = r*tril(x) if inverse is 'N' and trans is 'T' or
+ * inverse is 'I' and trans is 'N'
+ * wrk = tril(x)*r otherwise. */
+ len = n*n;
+ dcopy_(&len, MAT_BUFD(rk), &int1, wrk, &int1);
+ ld = MAX(1, n);
+ dtrmm_( (( inverse == 'N' && trans == 'T') || ( inverse == 'I'
+ && trans == 'N')) ? "R" : "L", "L", "N", "N", &n, &n,
+ &dbl1, MAT_BUFD(x) + ind + i*x->nrows, &ld, wrk, &ld);
+
+ /* x := (r*wrk' + wrk*r') if inverse is 'N' and trans is 'T'
+ * or inverse is 'I' and trans is 'N'
+ * x := (r'*wrk + wrk'*r) otherwise. */
+ dsyr2k_("L", ((inverse == 'N' && trans == 'T') ||
+ (inverse == 'I' && trans == 'N')) ? "N" : "T", &n, &n,
+ &dbl1, MAT_BUFD(rk), &ld, wrk, &ld, &dbl0, MAT_BUFD(x) +
+ ind + i*x->nrows, &ld);
+ }
+ ind += n*n;
+ }
+ free(wrk);
+
+ return Py_BuildValue("");
+}
+
+
+static char doc_scale2[] =
+ "Multiplication with square root of the Hessian.\n\n"
+ "scale2(lmbda, x, dims, mnl = 0, inverse = 'N')\n\n"
+ "Computes\n\n"
+ "Evaluates\n\n"
+ " x := H(lambda^{1/2}) * x (inverse is 'N')\n"
+ " x := H(lambda^{-1/2}) * x (inverse is 'I').\n\n"
+ "H is the Hessian of the logarithmic barrier.";
+
+static PyObject* scale2(PyObject *self, PyObject *args, PyObject *kwrds)
+{
+ matrix *lmbda, *x;
+ PyObject *dims, *O, *Ok;
+ char inverse = 'N';
+ double a, lx, x0, b, *c = NULL, *sql = NULL;
+ int m = 0, mk, i, j, len, int0 = 0, int1 = 1, maxn = 0, ind2;
+ char *kwlist[] = {"lmbda", "x", "dims", "mnl", "inverse", NULL};
+
+ if (!PyArg_ParseTupleAndKeywords(args, kwrds, "OOO|ic", kwlist, &lmbda,
+ &x, &dims, &m, &inverse)) return NULL;
+
+
+ /*
+ * For nonlinear and 'l' blocks:
+ *
+ * xk := xk ./ l (invers is 'N')
+ * xk := xk .* l (invers is 'I')
+ *
+ * where l is the first mnl + dims['l'] components of lmbda.
+ */
+
+ O = PyDict_GetItemString(dims, "l");
+ m += (int) PyInt_AsLong(O);
+ if (inverse == 'N')
+ dtbsv_("L", "N", "N", &m, &int0, MAT_BUFD(lmbda), &int1,
+ MAT_BUFD(x), &int1);
+ else
+ dtbmv_("L", "N", "N", &m, &int0, MAT_BUFD(lmbda), &int1,
+ MAT_BUFD(x), &int1);
+
+
+ /*
+ * For 'q' blocks, if inverse is 'N',
+ *
+ * xk := 1/a * [ l'*J*xk;
+ * xk[1:] - (xk[0] + l'*J*xk) / (l[0] + 1) * l[1:] ].
+ *
+ * If inverse is 'I',
+ *
+ * xk := a * [ l'*xk;
+ * xk[1:] + (xk[0] + l'*xk) / (l[0] + 1) * l[1:] ].
+ *
+ * a = sqrt(lambda_k' * J * lambda_k), l = lambda_k / a.
+ */
+
+ O = PyDict_GetItemString(dims, "q");
+ for (i = 0; i < (int) PyList_Size(O); i++){
+ Ok = PyList_GetItem(O, (Py_ssize_t) i);
+ mk = (int) PyInt_AsLong(Ok);
+ len = mk - 1;
+ a = dnrm2_(&len, MAT_BUFD(lmbda) + m + 1, &int1);
+ a = sqrt(MAT_BUFD(lmbda)[m] + a) * sqrt(MAT_BUFD(lmbda)[m] - a);
+ if (inverse == 'N')
+ lx = ( MAT_BUFD(lmbda)[m] * MAT_BUFD(x)[m] -
+ ddot_(&len, MAT_BUFD(lmbda) + m + 1, &int1, MAT_BUFD(x) + m
+ + 1, &int1) ) / a;
+ else
+ lx = ddot_(&mk, MAT_BUFD(lmbda) + m, &int1, MAT_BUFD(x) + m,
+ &int1) / a;
+ x0 = MAT_BUFD(x)[m];
+ MAT_BUFD(x)[m] = lx;
+ b = (x0 + lx) / (MAT_BUFD(lmbda)[m]/a + 1.0) / a;
+ if (inverse == 'N') b *= -1.0;
+ daxpy_(&len, &b, MAT_BUFD(lmbda) + m + 1, &int1,
+ MAT_BUFD(x) + m + 1, &int1);
+ if (inverse == 'N') a = 1.0 / a;
+ dscal_(&mk, &a, MAT_BUFD(x) + m, &int1);
+ m += mk;
+ }
+
+
+ /*
+ * For the 's' blocks, if inverse is 'N',
+ *
+ * xk := vec( diag(l)^{-1/2} * mat(xk) * diag(k)^{-1/2}).
+ *
+ * If inverse is 'I',
+ *
+ * xk := vec( diag(l)^{1/2} * mat(xk) * diag(k)^{1/2}).
+ *
+ * where l is kth block of lambda.
+ *
+ * We scale upper and lower triangular part of mat(xk) because the
+ * inverse operation will be applied to nonsymmetric matrices.
+ */
+
+ O = PyDict_GetItemString(dims, "s");
+ for (i = 0; i < (int) PyList_Size(O); i++){
+ Ok = PyList_GetItem(O, (Py_ssize_t) i);
+ maxn = MAX(maxn, (int) PyInt_AsLong(Ok));
+ }
+ if (!(c = (double *) calloc(maxn, sizeof(double))) ||
+ !(sql = (double *) calloc(maxn, sizeof(double)))){
+ free(c); free(sql);
+ return PyErr_NoMemory();
+ }
+ ind2 = m;
+ for (i = 0; i < (int) PyList_Size(O); i++){
+ Ok = PyList_GetItem(O, (Py_ssize_t) i);
+ mk = (int) PyInt_AsLong(Ok);
+ for (j = 0; j < mk; j++)
+ sql[j] = sqrt(MAT_BUFD(lmbda)[ind2 + j]);
+ for (j = 0; j < mk; j++){
+ dcopy_(&mk, sql, &int1, c, &int1);
+ b = sqrt(MAT_BUFD(lmbda)[ind2 + j]);
+ dscal_(&mk, &b, c, &int1);
+ if (inverse == 'N')
+ dtbsv_("L", "N", "N", &mk, &int0, c, &int1, MAT_BUFD(x) +
+ m + j*mk, &int1);
+ else
+ dtbmv_("L", "N", "N", &mk, &int0, c, &int1, MAT_BUFD(x) +
+ m + j*mk, &int1);
+ }
+ m += mk*mk;
+ ind2 += mk;
+ }
+ free(c); free(sql);
+
+ return Py_BuildValue("");
+}
+
+
+static char doc_pack[] =
+ "Copy x to y using packed storage.\n\n"
+ "pack(x, y, dims, mnl = 0, offsetx = 0, offsety = 0)\n\n"
+ "The vector x is an element of S, with the 's' components stored in\n"
+ "unpacked storage. On return, x is copied to y with the 's' \n"
+ "components stored in packed storage and the off-diagonal entries \n"
+ "scaled by sqrt(2).";
+
+static PyObject* pack(PyObject *self, PyObject *args, PyObject *kwrds)
+{
+ matrix *x, *y;
+ PyObject *O, *Ok, *dims;
+ double a;
+ int i, k, nlq = 0, ox = 0, oy = 0, np, iu, ip, int1 = 1, len, n;
+ char *kwlist[] = {"x", "y", "dims", "mnl", "offsetx", "offsety", NULL};
+
+ if (!PyArg_ParseTupleAndKeywords(args, kwrds, "OOO|iii", kwlist, &x,
+ &y, &dims, &nlq, &ox, &oy)) return NULL;
+
+ O = PyDict_GetItemString(dims, "l");
+ nlq += (int) PyInt_AsLong(O);
+
+ O = PyDict_GetItemString(dims, "q");
+ for (i = 0; i < (int) PyList_Size(O); i++){
+ Ok = PyList_GetItem(O, (Py_ssize_t) i);
+ nlq += (int) PyInt_AsLong(Ok);
+ }
+ dcopy_(&nlq, MAT_BUFD(x) + ox, &int1, MAT_BUFD(y) + oy, &int1);
+
+ O = PyDict_GetItemString(dims, "s");
+ for (i = 0, np = 0, iu = ox + nlq, ip = oy + nlq; i < (int)
+ PyList_Size(O); i++){
+ Ok = PyList_GetItem(O, (Py_ssize_t) i);
+ n = (int) PyInt_AsLong(Ok);
+ for (k = 0; k < n; k++){
+ len = n-k;
+ dcopy_(&len, MAT_BUFD(x) + iu + k*(n+1), &int1, MAT_BUFD(y) +
+ ip, &int1);
+ MAT_BUFD(y)[ip] /= sqrt(2.0);
+ ip += len;
+ }
+ np += n*(n+1)/2;
+ iu += n*n;
+ }
+
+ a = sqrt(2.0);
+ dscal_(&np, &a, MAT_BUFD(y) + oy + nlq, &int1);
+
+ return Py_BuildValue("");
+}
+
+
+static char doc_pack2[] =
+ "In-place version of pack().\n\n"
+ "pack2(x, y, dims, mnl = 0)\n\n"
+ "In-place version of pack(), which also accepts matrix arguments x.\n"
+ "The columns of x are elements of S, with the 's' components stored\n"
+ "in unpacked storage. On return, the 's' components are stored in\n"
+ "packed storage and the off-diagonal entries are scaled by sqrt(2).";
+
+static PyObject* pack2(PyObject *self, PyObject *args, PyObject *kwrds)
+{
+ matrix *x;
+ PyObject *O, *Ok, *dims;
+ double a = sqrt(2.0), *wrk;
+ int i, j, k, nlq = 0, iu, ip, len, n, maxn;
+ char *kwlist[] = {"x", "dims", "mnl", NULL};
+
+ if (!PyArg_ParseTupleAndKeywords(args, kwrds, "OO|i", kwlist, &x,
+ &dims, &nlq)) return NULL;
+
+ O = PyDict_GetItemString(dims, "l");
+ nlq += (int) PyInt_AsLong(O);
+
+ O = PyDict_GetItemString(dims, "q");
+ for (i = 0; i < (int) PyList_Size(O); i++){
+ Ok = PyList_GetItem(O, (Py_ssize_t) i);
+ nlq += (int) PyInt_AsLong(Ok);
+ }
+
+ O = PyDict_GetItemString(dims, "s");
+ for (i = 0, maxn = 0; i < (int) PyList_Size(O); i++){
+ Ok = PyList_GetItem(O, (Py_ssize_t) i);
+ maxn = MAX(maxn, (int) PyInt_AsLong(Ok));
+ }
+ if (!maxn) return Py_BuildValue("");
+ if (!(wrk = (double *) calloc(maxn * x->ncols, sizeof(double))))
+ return PyErr_NoMemory();
+
+ for (i = 0, iu = nlq, ip = nlq; i < (int) PyList_Size(O); i++){
+ Ok = PyList_GetItem(O, (Py_ssize_t) i);
+ n = (int) PyInt_AsLong(Ok);
+ for (k = 0; k < n; k++){
+ len = n-k;
+ dlacpy_(" ", &len, &(x->ncols), MAT_BUFD(x) + iu + k*(n+1),
+ &(x->nrows), wrk, &maxn);
+ for (j = 1; j < len; j++)
+ dscal_(&(x->ncols), &a, wrk + j, &maxn);
+ dlacpy_(" ", &len, &(x->ncols), wrk, &maxn, MAT_BUFD(x) + ip,
+ &(x->nrows));
+ ip += len;
+ }
+ iu += n*n;
+ }
+
+ free(wrk);
+ return Py_BuildValue("");
+}
+
+
+static char doc_unpack[] =
+ "Unpacks x into y.\n\n"
+ "unpack(x, y, dims, mnl = 0, offsetx = 0, offsety = 0)\n\n"
+ "The vector x is an element of S, with the 's' components stored in\n"
+ "unpacked storage and off-diagonal entries scaled by sqrt(2).\n"
+ "On return, x is copied to y with the 's' components stored in\n"
+ "unpacked storage.";
+
+static PyObject* unpack(PyObject *self, PyObject *args, PyObject *kwrds)
+{
+ matrix *x, *y;
+ PyObject *O, *Ok, *dims;
+ double a = 1.0 / sqrt(2.0);
+ int m = 0, ox = 0, oy = 0, int1 = 1, iu, ip, len, i, k, n;
+ char *kwlist[] = {"x", "y", "dims", "mnl", "offsetx", "offsety", NULL};
+
+ if (!PyArg_ParseTupleAndKeywords(args, kwrds, "OOO|iii", kwlist, &x,
+ &y, &dims, &m, &ox, &oy)) return NULL;
+
+ O = PyDict_GetItemString(dims, "l");
+ m += (int) PyInt_AsLong(O);
+
+ O = PyDict_GetItemString(dims, "q");
+ for (i = 0; i < (int) PyList_Size(O); i++){
+ Ok = PyList_GetItem(O, (Py_ssize_t) i);
+ m += (int) PyInt_AsLong(Ok);
+ }
+ dcopy_(&m, MAT_BUFD(x) + ox, &int1, MAT_BUFD(y) + oy, &int1);
+
+ O = PyDict_GetItemString(dims, "s");
+ for (i = 0, ip = ox + m, iu = oy + m; i < (int) PyList_Size(O); i++){
+ Ok = PyList_GetItem(O, (Py_ssize_t) i);
+ n = (int) PyInt_AsLong(Ok);
+ for (k = 0; k < n; k++){
+ len = n-k;
+ dcopy_(&len, MAT_BUFD(x) + ip, &int1, MAT_BUFD(y) + iu +
+ k*(n+1), &int1);
+ ip += len;
+ len -= 1;
+ dscal_(&len, &a, MAT_BUFD(y) + iu + k*(n+1) + 1, &int1);
+ }
+ iu += n*n;
+ }
+
+ return Py_BuildValue("");
+}
+
+
+static char doc_symm[] =
+ "Converts lower triangular matrix to symmetric.\n\n"
+ "symm(x, n, offset = 0)\n\n"
+ "Fills in the upper triangular part of the symmetric matrix stored\n"
+ "in x[offset : offset+n*n] using 'L' storage.";
+
+static PyObject* symm(PyObject *self, PyObject *args, PyObject *kwrds)
+{
+ matrix *x;
+ int n, ox = 0, k, len, int1 = 1;
+ char *kwlist[] = {"x", "n", "offset", NULL};
+
+ if (!PyArg_ParseTupleAndKeywords(args, kwrds, "Oi|i", kwlist, &x, &n,
+ &ox)) return NULL;
+
+ if (n > 1) for (k = 0; k < n; k++){
+ len = n-k-1;
+ dcopy_(&len, MAT_BUFD(x) + ox + k*(n+1) + 1, &int1, MAT_BUFD(x) +
+ ox + (k+1)*(n+1)-1, &n);
+ }
+ return Py_BuildValue("");
+}
+
+
+static char doc_sprod[] =
+ "The product x := (y o x).\n\n"
+ "sprod(x, y, dims, mnl = 0, diag = 'N')\n\n"
+ "If diag is 'D', the 's' part of y is diagonal and only the diagonal\n"
+ "is stored.";
+
+static PyObject* sprod(PyObject *self, PyObject *args, PyObject *kwrds)
+{
+ matrix *x, *y;
+ PyObject *dims, *O, *Ok;
+ int i, j, k, mk, len, maxn, ind = 0, ind2, int0 = 0, int1 = 1, ld;
+ double a, *A = NULL, dbl2 = 0.5, dbl0 = 0.0;
+ char diag = 'N';
+ char *kwlist[] = {"x", "y", "dims", "mnl", "diag", NULL};
+
+ if (!PyArg_ParseTupleAndKeywords(args, kwrds, "OOO|ic", kwlist, &x, &y,
+ &dims, &ind, &diag)) return NULL;
+
+
+ /*
+ * For nonlinear and 'l' blocks:
+ *
+ * yk o xk = yk .* xk
+ */
+
+ O = PyDict_GetItemString(dims, "l");
+ ind += (int) PyInt_AsLong(O);
+ dtbmv_("L", "N", "N", &ind, &int0, MAT_BUFD(y), &int1, MAT_BUFD(x),
+ &int1);
+
+
+ /*
+ * For 'q' blocks:
+ *
+ * [ l0 l1' ]
+ * yk o xk = [ ] * xk
+ * [ l1 l0*I ]
+ *
+ * where yk = (l0, l1).
+ */
+
+ O = PyDict_GetItemString(dims, "q");
+ for (i = 0; i < (int) PyList_Size(O); i++){
+ Ok = PyList_GetItem(O, (Py_ssize_t) i);
+ mk = (int) PyInt_AsLong(Ok);
+ a = ddot_(&mk, MAT_BUFD(y) + ind, &int1, MAT_BUFD(x) + ind, &int1);
+ len = mk - 1;
+ dscal_(&len, MAT_BUFD(y) + ind, MAT_BUFD(x) + ind + 1, &int1);
+ daxpy_(&len, MAT_BUFD(x) + ind, MAT_BUFD(y) + ind + 1, &int1,
+ MAT_BUFD(x) + ind + 1, &int1);
+ MAT_BUFD(x)[ind] = a;
+ ind += mk;
+ }
+
+
+ /*
+ * For the 's' blocks:
+ *
+ * yk o sk = .5 * ( Yk * mat(xk) + mat(xk) * Yk )
+ *
+ * where Yk = mat(yk) if diag is 'N' and Yk = diag(yk) if diag is 'D'.
+ */
+
+ O = PyDict_GetItemString(dims, "s");
+ for (i = 0, maxn = 0; i < (int) PyList_Size(O); i++){
+ Ok = PyList_GetItem(O, (Py_ssize_t) i);
+ maxn = MAX(maxn, (int) PyInt_AsLong(Ok));
+ }
+ if (diag == 'N'){
+ if (!(A = (double *) calloc(maxn * maxn, sizeof(double))))
+ return PyErr_NoMemory();
+ for (i = 0; i < (int) PyList_Size(O); ind += mk*mk, i++){
+ Ok = PyList_GetItem(O, (Py_ssize_t) i);
+ mk = (int) PyInt_AsLong(Ok);
+ len = mk*mk;
+ dcopy_(&len, MAT_BUFD(x) + ind, &int1, A, &int1);
+
+ if (mk > 1) for (k = 0; k < mk; k++){
+ len = mk - k - 1;
+ dcopy_(&len, A + k*(mk+1) + 1, &int1, A + (k+1)*(mk+1)-1,
+ &mk);
+ dcopy_(&len, MAT_BUFD(y) + ind + k*(mk+1) + 1, &int1,
+ MAT_BUFD(y) + ind + (k+1)*(mk+1)-1, &mk);
+ }
+
+ ld = MAX(1, mk);
+ dsyr2k_("L", "N", &mk, &mk, &dbl2, A, &ld, MAT_BUFD(y) + ind,
+ &ld, &dbl0, MAT_BUFD(x) + ind, &ld);
+ }
+ }
+ else {
+ if (!(A = (double *) calloc(maxn, sizeof(double))))
+ return PyErr_NoMemory();
+ for (i = 0, ind2 = ind; i < (int) PyList_Size(O); ind += mk*mk,
+ ind2 += mk, i++){
+ Ok = PyList_GetItem(O, (Py_ssize_t) i);
+ mk = (int) PyInt_AsLong(Ok);
+ for (k = 0; k < mk; k++){
+ len = mk - k;
+ dcopy_(&len, MAT_BUFD(y) + ind2 + k, &int1, A, &int1);
+ for (j = 0; j < len; j++) A[j] += MAT_BUFD(y)[ind2 + k];
+ dscal_(&len, &dbl2, A, &int1);
+ dtbmv_("L", "N", "N", &len, &int0, A, &int1, MAT_BUFD(x) +
+ ind + k * (mk+1), &int1);
+ }
+ }
+ }
+
+ free(A);
+ return Py_BuildValue("");
+}
+
+
+static char doc_sinv[] =
+ "The inverse of the product x := (y o x) when the 's' components of \n"
+ "y are diagonal.\n\n"
+ "sinv(x, y, dims, mnl = 0)";
+
+static PyObject* sinv(PyObject *self, PyObject *args, PyObject *kwrds)
+{
+ matrix *x, *y;
+ PyObject *dims, *O, *Ok;
+ int i, j, k, mk, len, maxn, ind = 0, ind2, int0 = 0, int1 = 1;
+ double a, c, d, alpha, *A = NULL, dbl2 = 0.5;
+ char *kwlist[] = {"x", "y", "dims", "mnl", NULL};
+
+ if (!PyArg_ParseTupleAndKeywords(args, kwrds, "OOO|i", kwlist, &x, &y,
+ &dims, &ind)) return NULL;
+
+
+ /*
+ * For nonlinear and 'l' blocks:
+ *
+ * yk o\ xk = yk .\ xk
+ */
+
+ O = PyDict_GetItemString(dims, "l");
+ ind += (int) PyInt_AsLong(O);
+ dtbsv_("L", "N", "N", &ind, &int0, MAT_BUFD(y), &int1, MAT_BUFD(x),
+ &int1);
+
+
+ /*
+ * For 'q' blocks:
+ *
+ * [ l0 -l1' ]
+ * yk o\ xk = 1/a^2 * [ ] * xk
+ * [ -l1 (a*I + l1*l1')/l0 ]
+ *
+ * where yk = (l0, l1) and a = l0^2 - l1'*l1.
+ */
+
+ O = PyDict_GetItemString(dims, "q");
+ for (i = 0; i < (int) PyList_Size(O); i++){
+ Ok = PyList_GetItem(O, (Py_ssize_t) i);
+ mk = (int) PyInt_AsLong(Ok);
+ len = mk - 1;
+ a = dnrm2_(&len, MAT_BUFD(y) + ind + 1, &int1);
+ a = (MAT_BUFD(y)[ind] + a) * (MAT_BUFD(y)[ind] - a);
+ c = MAT_BUFD(x)[ind];
+ d = ddot_(&len, MAT_BUFD(x) + ind + 1, &int1,
+ MAT_BUFD(y) + ind + 1, &int1);
+ MAT_BUFD(x)[ind] = c * MAT_BUFD(y)[ind] - d;
+ alpha = a / MAT_BUFD(y)[ind];
+ dscal_(&len, &alpha, MAT_BUFD(x) + ind + 1, &int1);
+ alpha = d / MAT_BUFD(y)[ind] - c;
+ daxpy_(&len, &alpha, MAT_BUFD(y) + ind + 1, &int1, MAT_BUFD(x) +
+ ind + 1, &int1);
+ alpha = 1.0 / a;
+ dscal_(&mk, &alpha, MAT_BUFD(x) + ind, &int1);
+ ind += mk;
+ }
+
+
+ /*
+ * For the 's' blocks:
+ *
+ * yk o\ sk = xk ./ gamma
+ *
+ * where gammaij = .5 * (yk_i + yk_j).
+ */
+
+ O = PyDict_GetItemString(dims, "s");
+ for (i = 0, maxn = 0; i < (int) PyList_Size(O); i++){
+ Ok = PyList_GetItem(O, (Py_ssize_t) i);
+ maxn = MAX(maxn, (int) PyInt_AsLong(Ok));
+ }
+ if (!(A = (double *) calloc(maxn, sizeof(double))))
+ return PyErr_NoMemory();
+ for (i = 0, ind2 = ind; i < (int) PyList_Size(O); ind += mk*mk,
+ ind2 += mk, i++){
+ Ok = PyList_GetItem(O, (Py_ssize_t) i);
+ mk = (int) PyInt_AsLong(Ok);
+ for (k = 0; k < mk; k++){
+ len = mk - k;
+ dcopy_(&len, MAT_BUFD(y) + ind2 + k, &int1, A, &int1);
+ for (j = 0; j < len; j++) A[j] += MAT_BUFD(y)[ind2 + k];
+ dscal_(&len, &dbl2, A, &int1);
+ dtbsv_("L", "N", "N", &len, &int0, A, &int1, MAT_BUFD(x) + ind
+ + k * (mk+1), &int1);
+ }
+ }
+
+ free(A);
+ return Py_BuildValue("");
+}
+
+
+
+static char doc_trisc[] =
+ "Sets the upper triangular part of the 's' components of x equal to\n"
+ "zero and scales the strictly lower triangular part\n\n"
+ "trisc(x, dims, offset = 0)";
+
+static PyObject* trisc(PyObject *self, PyObject *args, PyObject *kwrds)
+{
+ matrix *x;
+ double dbl0 = 0.0, dbl2 = 2.0;
+ int ox = 0, i, k, nk, len, int1 = 1;
+ PyObject *dims, *O, *Ok;
+ char *kwlist[] = {"x", "dims", "offset", NULL};
+
+ if (!PyArg_ParseTupleAndKeywords(args, kwrds, "OO|i", kwlist, &x,
+ &dims, &ox)) return NULL;
+
+ O = PyDict_GetItemString(dims, "l");
+ ox += (int) PyInt_AsLong(O);
+
+ O = PyDict_GetItemString(dims, "q");
+ for (i = 0; i < (int) PyList_Size(O); i++){
+ Ok = PyList_GetItem(O, (Py_ssize_t) i);
+ ox += (int) PyInt_AsLong(Ok);
+ }
+
+ O = PyDict_GetItemString(dims, "s");
+ for (k = 0; k < (int) PyList_Size(O); k++){
+ Ok = PyList_GetItem(O, (Py_ssize_t) k);
+ nk = (int) PyInt_AsLong(Ok);
+ for (i = 1; i < nk; i++){
+ len = nk - i;
+ dscal_(&len, &dbl0, MAT_BUFD(x) + ox + i*(nk+1) - 1, &nk);
+ dscal_(&len, &dbl2, MAT_BUFD(x) + ox + nk*(i-1) + i, &int1);
+ }
+ ox += nk*nk;
+ }
+
+ return Py_BuildValue("");
+}
+
+
+static char doc_triusc[] =
+ "Scales the strictly lower triangular part of the 's' components of\n"
+ "x by 0.5.\n\n"
+ "triusc(x, dims, offset = 0)";
+
+static PyObject* triusc(PyObject *self, PyObject *args, PyObject *kwrds)
+{
+ matrix *x;
+ double dbl5 = 0.5;
+ int ox = 0, i, k, nk, len, int1 = 1;
+ PyObject *dims, *O, *Ok;
+ char *kwlist[] = {"x", "dims", "offset", NULL};
+
+ if (!PyArg_ParseTupleAndKeywords(args, kwrds, "OO|i", kwlist, &x,
+ &dims, &ox)) return NULL;
+
+ O = PyDict_GetItemString(dims, "l");
+ ox += (int) PyInt_AsLong(O);
+
+ O = PyDict_GetItemString(dims, "q");
+ for (i = 0; i < (int) PyList_Size(O); i++){
+ Ok = PyList_GetItem(O, (Py_ssize_t) i);
+ ox += (int) PyInt_AsLong(Ok);
+ }
+
+ O = PyDict_GetItemString(dims, "s");
+ for (k = 0; k < (int) PyList_Size(O); k++){
+ Ok = PyList_GetItem(O, (Py_ssize_t) k);
+ nk = (int) PyInt_AsLong(Ok);
+ for (i = 1; i < nk; i++){
+ len = nk - i;
+ dscal_(&len, &dbl5, MAT_BUFD(x) + ox + nk*(i-1) + i, &int1);
+ }
+ ox += nk*nk;
+ }
+
+ return Py_BuildValue("");
+}
+
+
+static char doc_sdot[] =
+ "Inner product of two vectors in S.\n\n"
+ "sdot(x, y, dims, mnl= 0)";
+
+static PyObject* sdot(PyObject *self, PyObject *args, PyObject *kwrds)
+{
+ matrix *x, *y;
+ int m = 0, int1 = 1, i, k, nk, inc, len;
+ double a;
+ PyObject *dims, *O, *Ok;
+ char *kwlist[] = {"x", "y", "dims", "mnl", NULL};
+
+ if (!PyArg_ParseTupleAndKeywords(args, kwrds, "OOO|i", kwlist, &x, &y,
+ &dims, &m)) return NULL;
+
+ O = PyDict_GetItemString(dims, "l");
+ m += (int) PyInt_AsLong(O);
+
+ O = PyDict_GetItemString(dims, "q");
+ for (i = 0; i < (int) PyList_Size(O); i++){
+ Ok = PyList_GetItem(O, (Py_ssize_t) i);
+ m += (int) PyInt_AsLong(Ok);
+ }
+ a = ddot_(&m, MAT_BUFD(x), &int1, MAT_BUFD(y), &int1);
+
+ O = PyDict_GetItemString(dims, "s");
+ for (k = 0; k < (int) PyList_Size(O); k++){
+ Ok = PyList_GetItem(O, (Py_ssize_t) k);
+ nk = (int) PyInt_AsLong(Ok);
+ inc = nk+1;
+ a += ddot_(&nk, MAT_BUFD(x) + m, &inc, MAT_BUFD(y) + m, &inc);
+ for (i = 1; i < nk; i++){
+ len = nk - i;
+ a += 2.0 * ddot_(&len, MAT_BUFD(x) + m + i, &inc,
+ MAT_BUFD(y) + m + i, &inc);
+ }
+ m += nk*nk;
+ }
+
+ return Py_BuildValue("d", a);
+}
+
+
+static char doc_max_step[] =
+ "Returns min {t | x + t*e >= 0}\n\n."
+ "max_step(x, dims, mnl = 0, sigma = None)\n\n"
+ "e is defined as follows\n\n"
+ "- For the nonlinear and 'l' blocks: e is the vector of ones.\n"
+ "- For the 'q' blocks: e is the first unit vector.\n"
+ "- For the 's' blocks: e is the identity matrix.\n\n"
+ "When called with the argument sigma, also returns the eigenvalues\n"
+ "(in sigma) and the eigenvectors (in x) of the 's' components of x.\n";
+
+static PyObject* max_step(PyObject *self, PyObject *args, PyObject *kwrds)
+{
+ matrix *x, *sigma = NULL;
+ PyObject *dims, *O, *Ok;
+ int i, mk, len, maxn, ind = 0, ind2, int1 = 1, ld, Ns = 0, info, lwork,
+ *iwork = NULL, liwork, iwl, m;
+ double t = -FLT_MAX, dbl0 = 0.0, *work = NULL, wl, *Q = NULL,
+ *w = NULL;
+ char *kwlist[] = {"x", "dims", "mnl", "sigma", NULL};
+
+ if (!PyArg_ParseTupleAndKeywords(args, kwrds, "OO|iO", kwlist, &x,
+ &dims, &ind, &sigma)) return NULL;
+
+ O = PyDict_GetItemString(dims, "l");
+ ind += (int) PyInt_AsLong(O);
+ for (i = 0; i < ind; i++) t = MAX(t, -MAT_BUFD(x)[i]);
+
+ O = PyDict_GetItemString(dims, "q");
+ for (i = 0; i < (int) PyList_Size(O); i++){
+ Ok = PyList_GetItem(O, (Py_ssize_t) i);
+ mk = (int) PyInt_AsLong(Ok);
+ len = mk - 1;
+ t = MAX(t, dnrm2_(&len, MAT_BUFD(x) + ind + 1, &int1) -
+ MAT_BUFD(x)[ind]);
+ ind += mk;
+ }
+
+ O = PyDict_GetItemString(dims, "s");
+ Ns = (int) PyList_Size(O);
+ for (i = 0, maxn = 0; i < Ns; i++){
+ Ok = PyList_GetItem(O, (Py_ssize_t) i);
+ maxn = MAX(maxn, (int) PyInt_AsLong(Ok));
+ }
+ if (!maxn) return Py_BuildValue("d", (ind) ? t : 0.0);
+
+ lwork = -1;
+ liwork = -1;
+ ld = MAX(1, maxn);
+ if (sigma){
+ dsyevd_("V", "L", &maxn, NULL, &ld, NULL, &wl, &lwork, &iwl,
+ &liwork, &info);
+ }
+ else {
+ if (!(Q = (double *) calloc(maxn * maxn, sizeof(double))) ||
+ !(w = (double *) calloc(maxn, sizeof(double)))){
+ free(Q); free(w);
+ return PyErr_NoMemory();
+ }
+ dsyevr_("N", "I", "L", &maxn, NULL, &ld, &dbl0, &dbl0, &int1,
+ &int1, &dbl0, &maxn, NULL, NULL, &int1, NULL, &wl, &lwork,
+ &iwl, &liwork, &info);
+ }
+ lwork = (int) wl;
+ liwork = iwl;
+ if (!(work = (double *) calloc(lwork, sizeof(double))) ||
+ (!(iwork = (int *) calloc(liwork, sizeof(int))))){
+ free(Q); free(w); free(work); free(iwork);
+ return PyErr_NoMemory();
+ }
+ for (i = 0, ind2 = 0; i < Ns; i++){
+ Ok = PyList_GetItem(O, (Py_ssize_t) i);
+ mk = (int) PyInt_AsLong(Ok);
+ if (mk){
+ if (sigma){
+ dsyevd_("V", "L", &mk, MAT_BUFD(x) + ind, &mk,
+ MAT_BUFD(sigma) + ind2, work, &lwork, iwork, &liwork,
+ &info);
+ t = MAX(t, -MAT_BUFD(sigma)[ind2]);
+ }
+ else {
+ len = mk*mk;
+ dcopy_(&len, MAT_BUFD(x) + ind, &int1, Q, &int1);
+ ld = MAX(1, mk);
+ dsyevr_("N", "I", "L", &mk, Q, &mk, &dbl0, &dbl0, &int1,
+ &int1, &dbl0, &m, w, NULL, &int1, NULL, work, &lwork,
+ iwork, &liwork, &info);
+ t = MAX(t, -w[0]);
+ }
+ }
+ ind += mk*mk;
+ ind2 += mk;
+ }
+ free(work); free(iwork); free(Q); free(w);
+
+ return Py_BuildValue("d", (ind) ? t : 0.0);
+}
+
+static PyMethodDef misc_solvers__functions[] = {
+ {"scale", (PyCFunction) scale, METH_VARARGS|METH_KEYWORDS, doc_scale},
+ {"scale2", (PyCFunction) scale2, METH_VARARGS|METH_KEYWORDS,
+ doc_scale2},
+ {"pack", (PyCFunction) pack, METH_VARARGS|METH_KEYWORDS, doc_pack},
+ {"pack2", (PyCFunction) pack2, METH_VARARGS|METH_KEYWORDS, doc_pack2},
+ {"unpack", (PyCFunction) unpack, METH_VARARGS|METH_KEYWORDS,
+ doc_unpack},
+ {"symm", (PyCFunction) symm, METH_VARARGS|METH_KEYWORDS, doc_symm},
+ {"trisc", (PyCFunction) trisc, METH_VARARGS|METH_KEYWORDS, doc_trisc},
+ {"triusc", (PyCFunction) triusc, METH_VARARGS|METH_KEYWORDS,
+ doc_triusc},
+ {"sdot", (PyCFunction) sdot, METH_VARARGS|METH_KEYWORDS, doc_sdot},
+ {"sprod", (PyCFunction) sprod, METH_VARARGS|METH_KEYWORDS, doc_sprod},
+ {"sinv", (PyCFunction) sinv, METH_VARARGS|METH_KEYWORDS, doc_sinv},
+ {"max_step", (PyCFunction) max_step, METH_VARARGS|METH_KEYWORDS,
+ doc_max_step},
+ {NULL} /* Sentinel */
+};
+
+PyMODINIT_FUNC initmisc_solvers(void)
+{
+ PyObject *m;
+
+ m = Py_InitModule3("cvxopt.misc_solvers", misc_solvers__functions,
+ misc_solvers__doc__);
+
+ if (import_cvxopt() < 0) return;
+}
diff --git a/src/C/sparse.c b/src/C/sparse.c
index 10b00c0..89d2706 100644
--- a/src/C/sparse.c
+++ b/src/C/sparse.c
@@ -1,7 +1,7 @@
/*
* Copyright 2004-2008 J. Dahl and L. Vandenberghe.
*
- * This file is part of CVXOPT version 0.9.3.
+ * This file is part of CVXOPT version 1.0.
*
* CVXOPT is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
@@ -2490,8 +2490,43 @@ static int spmatrix_set_size(spmatrix *self, PyObject *value, void *closure)
{
if (!value) PY_ERR_INT(PyExc_TypeError,"size attribute cannot be deleted");
- PY_ERR_INT(PyExc_TypeError, "sparse reshape not implemented");
- return -1;
+ if (!PyTuple_Check(value) || PyTuple_Size(value) != 2)
+ PY_ERR_INT(PyExc_TypeError, "can only assign a 2-tuple to size");
+
+ if (!PyInt_Check(PyTuple_GET_ITEM(value, 0)) ||
+ !PyInt_Check(PyTuple_GET_ITEM(value, 1)))
+ PY_ERR_INT(PyExc_TypeError, "invalid size tuple");
+
+ int m = PyInt_AS_LONG(PyTuple_GET_ITEM(value, 0));
+ int n = PyInt_AS_LONG(PyTuple_GET_ITEM(value, 1));
+
+ if (m<0 || n<0)
+ PY_ERR_INT(PyExc_TypeError, "dimensions must be non-negative");
+
+ if (m*n != SP_NROWS(self)*SP_NCOLS(self))
+ PY_ERR_INT(PyExc_TypeError, "number of elements in matrix cannot change");
+
+ int_t *colptr = calloc((n+1),sizeof(int_t));
+ if (!colptr) PY_ERR_INT(PyExc_MemoryError, "insufficient memory");
+
+ int j, k, in, jn;
+ for (j=0; j<SP_NCOLS(self); j++) {
+ for (k=SP_COL(self)[j]; k<SP_COL(self)[j+1]; k++) {
+ jn = (SP_ROW(self)[k] + j*SP_NROWS(self)) / m;
+ in = (SP_ROW(self)[k] + j*SP_NROWS(self)) % m;
+ colptr[jn+1]++;
+ SP_ROW(self)[k] = in;
+ }
+ }
+
+ for (j=1; j<n+1; j++) colptr[j] += colptr[j-1];
+
+ free(SP_COL(self));
+ SP_COL(self) = colptr;
+ SP_NROWS(self) = m;
+ SP_NCOLS(self) = n;
+
+ return 0;
}
static PyObject * spmatrix_get_typecode(matrix *self, void *closure)
@@ -2924,7 +2959,7 @@ spmatrix_ass_subscr(spmatrix* self, PyObject* args, PyObject* value)
if (val_id > id)
PY_ERR_INT(PyExc_TypeError, "invalid type in assignment");
- /* assigment value is matrix or number ? */
+ /* assignment value is matrix or number ? */
if (PY_NUMBER(value)) {
if (convert_num[id](&val, value, 1, 0))
PY_ERR_INT(PyExc_TypeError, "invalid argument type");
@@ -2942,7 +2977,7 @@ spmatrix_ass_subscr(spmatrix* self, PyObject* args, PyObject* value)
/* single integer */
if (PyInt_Check(args)) {
if (itype != 'n')
- PY_ERR_INT(PyExc_IndexError, "incompatible sizes in assigment");
+ PY_ERR_INT(PyExc_IndexError, "incompatible sizes in assignment");
i = PyInt_AsLong(args);
if ( i<-SP_LGT(self) || i >= SP_LGT(self) )
@@ -2962,7 +2997,7 @@ spmatrix_ass_subscr(spmatrix* self, PyObject* args, PyObject* value)
/* integer matrix list */
if (PyList_Check(args) || Matrix_Check(args) || PySlice_Check(args)) {
-
+
if (!(Il = create_indexlist(SP_LGT(self), args))) {
if (decref_val) { Py_DECREF(value); }
return -1;
@@ -2974,12 +3009,12 @@ spmatrix_ass_subscr(spmatrix* self, PyObject* args, PyObject* value)
if (((itype == 'd') &&
((lgtI != MAT_LGT(value) || MAT_NCOLS(value) != 1))) ||
(((itype == 's') &&
- ((lgtI != SP_LGT(value)) || SP_NROWS(value) != 1)))) {
+ ((lgtI != SP_LGT(value)) || SP_NCOLS(value) != 1)))) {
if (!Matrix_Check(args)) { Py_DECREF(Il); }
if (decref_val) { Py_DECREF(value); }
- PY_ERR_INT(PyExc_TypeError, "incompatible sizes in assigment");
+ PY_ERR_INT(PyExc_TypeError, "incompatible sizes in assignment");
}
-
+
/* ass. argument is dense matrix or number */
if (itype == 'd' || itype == 'n') {
@@ -3222,7 +3257,7 @@ spmatrix_ass_subscr(spmatrix* self, PyObject* args, PyObject* value)
if (!Matrix_Check(argI)) { Py_DECREF(Il); }
if (!Matrix_Check(argJ)) { Py_DECREF(Jl); }
if (decref_val) { Py_DECREF(value); }
- PY_ERR_INT(PyExc_TypeError, "incompatible size of assigment");
+ PY_ERR_INT(PyExc_TypeError, "incompatible size of assignment");
}
/* ass. argument is dense matrix or number */
diff --git a/src/C/umfpack.c b/src/C/umfpack.c
index 4edc535..ccea2bc 100644
--- a/src/C/umfpack.c
+++ b/src/C/umfpack.c
@@ -1,7 +1,7 @@
/*
* Copyright 2004-2008 J. Dahl and L. Vandenberghe.
*
- * This file is part of CVXOPT version 0.9.3.
+ * This file is part of CVXOPT version 1.0.
*
* CVXOPT is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
diff --git a/src/python/__init__.py b/src/python/__init__.py
index dbd9b75..48f8e8a 100644
--- a/src/python/__init__.py
+++ b/src/python/__init__.py
@@ -34,7 +34,7 @@ def normal(nrows, ncols=1, mean=0.0, std=1.0):
from cvxopt.base import matrix
from random import gauss
return matrix([gauss(mean, std) for k in xrange(nrows*ncols)],
- (nrows,ncols) )
+ (nrows,ncols), 'd' )
return gsl.normal(nrows, ncols, mean, std)
@@ -65,7 +65,7 @@ def uniform(nrows, ncols=1, a=0, b=1):
from cvxopt.base import matrix
from random import uniform
return matrix([uniform(a, b) for k in xrange(nrows*ncols)],
- (nrows,ncols) )
+ (nrows,ncols), 'd' )
return gsl.uniform(nrows, ncols, a, b)
diff --git a/src/python/coneprog.py b/src/python/coneprog.py
index 9d10f9f..575f2ef 100644
--- a/src/python/coneprog.py
+++ b/src/python/coneprog.py
@@ -4,7 +4,7 @@ Solver for linear, second-order cone and semidefinite programming.
# Copyright 2004-2008 J. Dahl and L. Vandenberghe.
#
-# This file is part of CVXOPT version 0.9.3.
+# This file is part of CVXOPT version 1.0.
#
# CVXOPT is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
@@ -278,6 +278,9 @@ def conelp(c, G, h, dims = None, A = None, b = None, primalstart = None,
EXPON = 3
STEP = 0.99
+ try: DEBUG = options['debug']
+ except KeyError: DEBUG = False
+
try: MAXITERS = options['maxiters']
except KeyError: MAXITERS = 100
else:
@@ -431,7 +434,7 @@ def conelp(c, G, h, dims = None, A = None, b = None, primalstart = None,
elif kktsolver == 'qr':
factor = misc.kkt_qr(G, dims, A)
elif kktsolver == 'chol':
- factor = kkt_chol(G, dims, A)
+ factor = misc.kkt_chol(G, dims, A)
else:
factor = misc.kkt_chol2(G, dims, A)
def kktsolver(W):
@@ -656,6 +659,7 @@ def conelp(c, G, h, dims = None, A = None, b = None, primalstart = None,
z[ind : ind+m*m : m+1] += a
ind += m**2
+
tau, kappa = 1.0, 1.0
rx, hrx = xnewcopy(c), xnewcopy(c)
@@ -916,7 +920,7 @@ def conelp(c, G, h, dims = None, A = None, b = None, primalstart = None,
# but applies iterative refinement.
if iters == 0:
- if refinement:
+ if refinement or DEBUG:
wx, wy = xnewcopy(c), ynewcopy(b)
wz, ws = matrix(0.0, (cdim, 1)), matrix(0.0, (cdim, 1))
wtau, wkappa = matrix(0.0), matrix(0.0)
@@ -926,7 +930,7 @@ def conelp(c, G, h, dims = None, A = None, b = None, primalstart = None,
wtau2, wkappa2 = matrix(0.0), matrix(0.0)
def f6(x, y, z, tau, s, kappa):
- if refinement:
+ if refinement or DEBUG:
xcopy(x, wx)
ycopy(y, wy)
blas.copy(z, wz)
@@ -950,6 +954,17 @@ def conelp(c, G, h, dims = None, A = None, b = None, primalstart = None,
tau[0] += wtau2[0]
blas.axpy(ws2, s)
kappa[0] += wkappa2[0]
+ if DEBUG:
+ res(x, y, z, tau, s, kappa, wx, wy, wz, wtau, ws, wkappa,
+ W, dg, lmbda)
+ print "KKT residuals"
+ print " 'x': %e" %math.sqrt(xdot(wx, wx))
+ print " 'y': %e" %math.sqrt(ydot(wy, wy))
+ print " 'z': %e" %misc.snrm2(wz, dims)
+ print " 'tau': %e" %abs(wtau[0])
+ print " 's': %e" %misc.snrm2(ws, dims)
+ print " 'kappa': %e" %abs(wkappa[0])
+
mu = blas.nrm2(lmbda)**2 / (1 + cdim_diag)
sigma = 0.0
@@ -1011,6 +1026,7 @@ def conelp(c, G, h, dims = None, A = None, b = None, primalstart = None,
misc.sprod(ws3, dz, dims)
wkappa3 = dtau[0] * dkappa[0]
+
# Maximum step to boundary.
#
# If i is 1, also compute eigenvalue decomposition of the 's'
@@ -1396,6 +1412,13 @@ def coneqp(P, q, G = None, h = None, dims = None, A = None, b = None,
STEP = 0.99
EXPON = 3
+ try: DEBUG = options['debug']
+ except KeyError: DEBUG = False
+
+ # Use Mehrotra correction or not.
+ try: correction = options['use_correction']
+ except KeyError: correction = True
+
try: MAXITERS = options['maxiters']
except KeyError: MAXITERS = 100
else:
@@ -1687,7 +1710,6 @@ def coneqp(P, q, G = None, h = None, dims = None, A = None, b = None,
sigs = matrix(0.0, (sum(dims['s']), 1))
sigz = matrix(0.0, (sum(dims['s']), 1))
-
if show_progress:
print "% 10s% 12s% 10s% 8s% 7s" %("pcost", "dcost", "gap", "pres",
"dres")
@@ -1816,7 +1838,7 @@ def coneqp(P, q, G = None, h = None, dims = None, A = None, b = None,
# iterative refinement.
if iters == 0:
- if refinement:
+ if refinement or DEBUG:
wx, wy = xnewcopy(q), ynewcopy(b)
wz, ws = matrix(0.0, (cdim,1)), matrix(0.0, (cdim,1))
if refinement:
@@ -1824,7 +1846,7 @@ def coneqp(P, q, G = None, h = None, dims = None, A = None, b = None,
wz2, ws2 = matrix(0.0, (cdim,1)), matrix(0.0, (cdim,1))
def f4(x, y, z, s):
- if refinement:
+ if refinement or DEBUG:
xcopy(x, wx)
ycopy(y, wy)
blas.copy(z, wz)
@@ -1841,6 +1863,14 @@ def coneqp(P, q, G = None, h = None, dims = None, A = None, b = None,
yaxpy(wy2, y)
blas.axpy(wz2, z)
blas.axpy(ws2, s)
+ if DEBUG:
+ res(x, y, z, s, wx, wy, wz, ws, W, lmbda)
+ print "KKT residuals:"
+ print " 'x': %e" %math.sqrt(xdot(wx, wx))
+ print " 'y': %e" %math.sqrt(ydot(wy, wy))
+ print " 'z': %e" %misc.snrm2(wz, dims)
+ print " 's': %e" %misc.snrm2(ws, dims)
+
mu = gap / (dims['l'] + len(dims['q']) + sum(dims['s']))
sigma, eta = 0.0, 0.0
@@ -1862,7 +1892,7 @@ def coneqp(P, q, G = None, h = None, dims = None, A = None, b = None,
# = -lmbdasq - dsa o dza + sigma * mu * e (if i is 1),
# where ds, dz are solution for i is 0.
blas.scal(0.0, ds)
- if i == 1:
+ if correction and i == 1:
blas.axpy(ws3, ds, alpha = -1.0)
blas.axpy(lmbdasq, ds, n = dims['l'] + sum(dims['q']),
alpha = -1.0)
@@ -1894,7 +1924,7 @@ def coneqp(P, q, G = None, h = None, dims = None, A = None, b = None,
raise ArithmeticError, "singular KKT matrix"
# Save ds o dz for Mehrotra correction
- if i == 0:
+ if correction and i == 0:
blas.copy(ds, ws3)
misc.sprod(ws3, dz, dims)
@@ -2134,7 +2164,7 @@ def lp(c, G, h, A = None, b = None, solver = None, primalstart = None,
except ImportError: raise ValueError, "invalid option "\
"(solver='glpk'): cvxopt.glpk is not installed"
glpk.options = options
- status, x, z, y = glpk.solvelp(c,G,h,A,b)
+ status, x, z, y = glpk.lp(c,G,h,A,b)
if status == 'optimal':
s = matrix(h)
base.gemv(G, x, s, alpha=-1.0, beta=1.0)
@@ -2153,7 +2183,7 @@ def lp(c, G, h, A = None, b = None, solver = None, primalstart = None,
mosek.options = options['MOSEK']
else:
mosek.options = {}
- solsta, x, z, y = mosek.solvelp(c, G, h, A, b)
+ solsta, x, z, y = mosek.lp(c, G, h, A, b)
if solsta is msk.solsta.optimal:
s = matrix(h)
base.gemv(G, x, s, alpha=-1.0, beta=1.0)
@@ -2784,7 +2814,7 @@ def qp(P, q, G = None, h = None, A = None, b = None, solver = None,
mosek.options = options['MOSEK']
else:
mosek.options = {}
- solsta, x, z, y = mosek.solveqp(P, q, G, h, A, b)
+ solsta, x, z, y = mosek.qp(P, q, G, h, A, b)
m = G.size[0]
if solsta == pymosek.solsta.optimal:
diff --git a/src/python/cvxprog.py b/src/python/cvxprog.py
index 857973c..ee2e5c1 100644
--- a/src/python/cvxprog.py
+++ b/src/python/cvxprog.py
@@ -8,7 +8,7 @@ to the quadratic programming solver from MOSEK.
# Copyright 2004-2008 J. Dahl and L. Vandenberghe.
#
-# This file is part of CVXOPT version 0.9.3.
+# This file is part of CVXOPT version 1.0.
#
# CVXOPT is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
@@ -318,6 +318,9 @@ def cpl(c, F, G = None, h = None, dims = None, A = None, b = None,
EXPON = 3
MAX_RELAXED_ITERS = 8
+ try: DEBUG = options['debug']
+ except KeyError: DEBUG = False
+
try: MAXITERS = options['maxiters']
except KeyError: MAXITERS = 100
else:
@@ -548,11 +551,10 @@ def cpl(c, F, G = None, h = None, dims = None, A = None, b = None,
print "% 10s% 12s% 10s% 8s% 7s" %("pcost", "dcost", "gap", "pres",
"dres")
-
relaxed_iters = 0
for iters in xrange(MAXITERS):
- if refinement:
+ if refinement or DEBUG:
# We need H to compute residuals of KKT equations.
f, Df, H = F(x, z[:mnl])
else:
@@ -578,7 +580,7 @@ def cpl(c, F, G = None, h = None, dims = None, A = None, b = None,
"a user-provided kktsolver"
fDf = Df
- if refinement:
+ if refinement or DEBUG:
if type(H) is matrix or type(H) is spmatrix:
if customx: raise ValueError, "use of non-vector type "\
"for x requires function valued H"
@@ -811,7 +813,7 @@ def cpl(c, F, G = None, h = None, dims = None, A = None, b = None,
# iterative refinement.
if iters == 0:
- if refinement:
+ if refinement or DEBUG:
wx, wy = xnewcopy(c), ynewcopy(b)
wz = matrix(0.0, (mnl + cdim, 1))
ws = matrix(0.0, (mnl + cdim, 1))
@@ -821,7 +823,7 @@ def cpl(c, F, G = None, h = None, dims = None, A = None, b = None,
ws2 = matrix(0.0, (mnl + cdim, 1))
def f4(x, y, z, s):
- if refinement:
+ if refinement or DEBUG:
xcopy(x, wx)
ycopy(y, wy)
blas.copy(z, wz)
@@ -838,6 +840,14 @@ def cpl(c, F, G = None, h = None, dims = None, A = None, b = None,
yaxpy(wy2, y)
blas.axpy(wz2, z)
blas.axpy(ws2, s)
+ if DEBUG:
+ res(x, y, z, s, wx, wy, wz, ws)
+ print "KKT residuals:"
+ print " 'x': %e" %math.sqrt(xdot(wx, wx))
+ print " 'y': %e" %math.sqrt(ydot(wy, wy))
+ print " 'z': %e" %misc.snrm2(wz, dims, mnl)
+ print " 's': %e" %misc.snrm2(ws, dims, mnl)
+
sigma, eta = 0.0, 0.0
for i in [0, 1]:
diff --git a/src/python/info.py b/src/python/info.py
index a5f54ba..ed4e178 100644
--- a/src/python/info.py
+++ b/src/python/info.py
@@ -1,8 +1,8 @@
-version = '0.9.3'
+version = '1.0'
def license(): print(
"""
-CVXOPT version 0.9.3. Copyright (c) 2004-2008 J. Dahl and L. Vandenberghe.
+CVXOPT version 1.0. Copyright (c) 2004-2008 J. Dahl and L. Vandenberghe.
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
diff --git a/src/python/misc.py b/src/python/misc.py
index 8672260..1b7378d 100644
--- a/src/python/misc.py
+++ b/src/python/misc.py
@@ -1,6 +1,6 @@
# Copyright 2004-2008 J. Dahl and L. Vandenberghe.
#
-# This file is part of CVXOPT version 0.9.3.
+# This file is part of CVXOPT version 1.0.
#
# CVXOPT is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
@@ -16,36 +16,43 @@
# along with this program. If not, see <http://www.gnu.org/licenses/>.
import math
-from cvxopt import base, blas, lapack, cholmod
+from cvxopt import base, blas, lapack, cholmod, misc_solvers
from cvxopt.base import matrix, spmatrix
__all__ = []
+use_C = True
-def scale(x, W, trans = 'N', inverse = 'N'):
-
- # Computes
- #
- # x := W*x (trans is 'N', inverse = 'N')
- # x := W^T*x (trans is 'T', inverse = 'N')
- # x := W^{-1}*x (trans is 'N', inverse = 'I')
- # x := W^{-T}*x (trans is 'T', inverse = 'I').
- #
- # x is a dense 'd' matrix.
- #
- # W is a dictionary with entries:
- #
- # - W['dnl']: positive vector
- # - W['dnli']: componentwise inverse of W['dnl']
- # - W['d']: positive vector
- # - W['di']: componentwise inverse of W['d']
- # - W['v']: lists of 2nd order cone vectors with unit hyperbolic norms
- # - W['beta']: list of positive numbers
- # - W['r']: list of square matrices
- # - W['rti']: list of square matrices. rti[k] is the inverse transpose
- # of r[k].
- #
- # The 'dnl' and 'dnli' entries are optional, and only present when the
- # function is called from the nonlinear solver.
+if use_C:
+ scale = misc_solvers.scale
+else:
+ def scale(x, W, trans = 'N', inverse = 'N'):
+ """
+ Applies Nesterov-Todd scaling or its inverse.
+
+ Computes
+
+ x := W*x (trans is 'N', inverse = 'N')
+ x := W^T*x (trans is 'T', inverse = 'N')
+ x := W^{-1}*x (trans is 'N', inverse = 'I')
+ x := W^{-T}*x (trans is 'T', inverse = 'I').
+
+ x is a dense 'd' matrix.
+
+ W is a dictionary with entries:
+
+ - W['dnl']: positive vector
+ - W['dnli']: componentwise inverse of W['dnl']
+ - W['d']: positive vector
+ - W['di']: componentwise inverse of W['d']
+ - W['v']: lists of 2nd order cone vectors with unit hyperbolic norms
+ - W['beta']: list of positive numbers
+ - W['r']: list of square matrices
+ - W['rti']: list of square matrices. rti[k] is the inverse transpose
+ of r[k].
+
+ The 'dnl' and 'dnli' entries are optional, and only present when the
+ function is called from the nonlinear solver.
+ """
ind = 0
@@ -155,12 +162,18 @@ def scale(x, W, trans = 'N', inverse = 'N'):
ind += n**2
-def scale2(lmbda, x, dims, mnl = 0, inverse = 'N'):
+if use_C:
+ scale2 = misc_solvers.scale2
+else:
+ def scale2(lmbda, x, dims, mnl = 0, inverse = 'N'):
+ """
+ Evaluates
- # x := H(lambda^{1/2}) * x (inverse is 'N')
- # x := H(lambda^{-1/2}) * x (inverse is 'I')
- #
- # H is the Hessian of the logarithmic barrier.
+ x := H(lambda^{1/2}) * x (inverse is 'N')
+ x := H(lambda^{-1/2}) * x (inverse is 'I').
+
+ H is the Hessian of the logarithmic barrier.
+ """
# For the nonlinear and 'l' blocks,
@@ -233,11 +246,13 @@ def scale2(lmbda, x, dims, mnl = 0, inverse = 'N'):
def compute_scaling(s, z, lmbda, dims, mnl = None):
+ """
+ Returns the Nesterov-Todd scaling W at points s and z, and stores the
+ scaled variable in lmbda.
+
+ W * z = W^{-T} * s = lmbda.
- # Returns the Nesterov-Todd scaling W at points s and z, and stores
- # the scaled variable in lmbda.
- #
- # W * z = W^{-T} * s = lmbda.
+ """
W = {}
@@ -403,19 +418,20 @@ def compute_scaling(s, z, lmbda, dims, mnl = None):
def update_scaling(W, lmbda, s, z):
-
- # Updates the Nesterov-Todd scaling matrix W and the scaled variable
- # lmbda so that on exit
- #
- # W * zt = W^{-T} * st = lmbda.
- #
- # On entry, the nonlinear, 'l' and 'q' components of the arguments s
- # and z contain W^{-T}*st and W*zt, i.e, the new iterates in the
- # current scaling.
- #
- # The 's' components contain the factors Ls, Lz in a factorization of
- # the new iterates in the current scaling, W^{-T}*st = Ls*Ls',
- # W*zt = Lz*Lz'.
+ """
+ Updates the Nesterov-Todd scaling matrix W and the scaled variable
+ lmbda so that on exit
+
+ W * zt = W^{-T} * st = lmbda.
+
+ On entry, the nonlinear, 'l' and 'q' components of the arguments s
+ and z contain W^{-T}*st and W*zt, i.e, the new iterates in the current
+ scaling.
+
+ The 's' components contain the factors Ls, Lz in a factorization of
+ the new iterates in the current scaling, W^{-T}*st = Ls*Ls',
+ W*zt = Lz*Lz'.
+ """
# Nonlinear and 'l' blocks
@@ -616,15 +632,20 @@ def update_scaling(W, lmbda, s, z):
ind3 += m
-def pack(x, y, dims, mnl = 0, offsetx = 0, offsety = 0):
-
- # The vector x is an element of S, with the 's' components stored
- # in unpacked storage. On return, x is copied to y with the 's'
- # components stored in packed storage and the off-diagonal entries
- # scaled by sqrt(2).
+if use_C:
+ pack = misc_solvers.pack
+else:
+ def pack(x, y, dims, mnl = 0, offsetx = 0, offsety = 0):
+ """
+ Copy x to y using packed storage.
+
+ The vector x is an element of S, with the 's' components stored in
+ unpacked storage. On return, x is copied to y with the 's' components
+ stored in packed storage and the off-diagonal entries scaled by
+ sqrt(2).
+ """
nlq = mnl + dims['l'] + sum(dims['q'])
- np = sum([ n*(n+1)/2 for n in dims['s'] ])
blas.copy(x, y, n = nlq, offsetx = offsetx, offsety = offsety)
iu, ip = offsetx + nlq, offsety + nlq
for n in dims['s']:
@@ -633,18 +654,46 @@ def pack(x, y, dims, mnl = 0, offsetx = 0, offsety = 0):
y[ip] /= math.sqrt(2)
ip += n-k
iu += n**2
+ np = sum([ n*(n+1)/2 for n in dims['s'] ])
blas.scal(math.sqrt(2.0), y, n = np, offset = offsety+nlq)
-def unpack(x, y, dims, mnl = 0, offsetx = 0, offsety = 0):
+if use_C:
+ pack2 = misc_solvers.pack2
+else:
+ def pack2(x, dims, mnl = 0):
+ """
+ In-place version of pack(), which also accepts matrix arguments x.
+ The columns of x are elements of S, with the 's' components stored in
+ unpacked storage. On return, the 's' components are stored in packed
+ storage and the off-diagonal entries are scaled by sqrt(2).
+ """
+
+ if not dims['s']: return
+ iu = mnl + dims['l'] + sum(dims['q'])
+ ip = iu
+ for n in dims['s']:
+ for k in xrange(n):
+ x[ip, :] = x[iu + (n+1)*k, :]
+ x[ip + 1 : ip+n-k, :] = x[iu + (n+1)*k + 1: iu + n*(k+1), :] \
+ * math.sqrt(2.0)
+ ip += n - k
+ iu += n**2
+ np = sum([ n*(n+1)/2 for n in dims['s'] ])
+
- # The vector x is an element of S, with the 's' components stored
- # in unpacked storage and off-diagonal entries scaled by sqrt(2).
- # On return, x is copied to y with the 's' components stored in
- # unpacked storage and off-diagonal entries scaled by sqrt(2).
+if use_C:
+ unpack = misc_solvers.unpack
+else:
+ def unpack(x, y, dims, mnl = 0, offsetx = 0, offsety = 0):
+ """
+ The vector x is an element of S, with the 's' components stored
+ in unpacked storage and off-diagonal entries scaled by sqrt(2).
+ On return, x is copied to y with the 's' components stored in
+ unpacked storage.
+ """
nlq = mnl + dims['l'] + sum(dims['q'])
- nu = sum([ n**2 for n in dims['s'] ])
blas.copy(x, y, n = nlq, offsetx = offsetx, offsety = offsety)
iu, ip = offsety+nlq, offsetx+nlq
for n in dims['s']:
@@ -653,12 +702,17 @@ def unpack(x, y, dims, mnl = 0, offsetx = 0, offsety = 0):
y[iu+k*(n+1)] *= math.sqrt(2)
ip += n-k
iu += n**2
+ nu = sum([ n**2 for n in dims['s'] ])
blas.scal(1.0/math.sqrt(2.0), y, n = nu, offset = offsety+nlq)
-def sdot(x, y, dims, mnl = 0):
-
- # Returns the inner product of two vectors in S.
+if use_C:
+ sdot = misc_solvers.sdot
+else:
+ def sdot(x, y, dims, mnl = 0):
+ """
+ Inner product of two vectors in S.
+ """
ind = mnl + dims['l'] + sum(dims['q'])
a = blas.dot(x, y, n = ind)
@@ -699,85 +753,116 @@ def sdot2(x, y):
def snrm2(x, dims, mnl = 0):
-
- # Returns the norm of a vector in S
+ """
+ Returns the norm of a vector in S
+ """
return math.sqrt(sdot(x, x, dims, mnl))
-def sgemv(A, x, y, dims, trans = 'N', alpha = 1.0, beta = 0.0, m = None,
- n = None, offsetA = 0, offsetx = 0, offsety = 0):
+if use_C:
+ trisc = misc_solvers.trisc
+else:
+ def trisc(x, dims, offset = 0):
+ """
+ Sets upper triangular part of the 's' components of x equal to zero
+ and scales the strictly lower triangular part by 2.0.
+ """
+
+ m = dims['l'] + sum(dims['q']) + sum([ k**2 for k in dims['s'] ])
+ ind = offset + dims['l'] + sum(dims['q'])
+ for mk in dims['s']:
+ for j in xrange(1, mk):
+ blas.scal(0.0, x, n = mk-j, inc = mk, offset =
+ ind + j*(mk + 1) - 1)
+ blas.scal(2.0, x, offset = ind + mk*(j-1) + j, n = mk-j)
+ ind += mk**2
- # A is a matrix or spmatrix of size (N, n) where
- #
- # N = dims['l'] + sum(dims['q']) + sum( k**2 for k in dims['s'] ).
- #
- # If trans is 'N':
- #
- # y := alpha*A*x + beta * y (trans = 'N').
- #
- # x is a vector of length n. y is a vector of length N.
- #
- # If trans is 'T':
- #
- # y := alpha*A'*x + beta * y (trans = 'T').
- #
- # x is a vector of length N. y is a vector of length n.
- #
- # The 's' components in S are stored in unpacked 'L' storage.
- if m is None: m = A.size[0]
- if n is None: n = A.size[1]
+if use_C:
+ triusc = misc_solvers.triusc
+else:
+ def triusc(x, dims, offset = 0):
+ """
+ Scales the strictly lower triangular part of the 's' components of x
+ by 0.5.
+ """
+
+ m = dims['l'] + sum(dims['q']) + sum([ k**2 for k in dims['s'] ])
+ ind = offset + dims['l'] + sum(dims['q'])
+ for mk in dims['s']:
+ for j in xrange(1, mk):
+ blas.scal(0.5, x, offset = ind + mk*(j-1) + j, n = mk-j)
+ ind += mk**2
- if trans == 'T' and alpha:
- ind = offsetx + dims['l'] + sum(dims['q'])
- for mk in dims['s']:
- # Set upper triangular part of x to zero and scale strict
- # lower triangular part by 2.
- for j in xrange(1, mk):
- blas.scal(0.0, x, n = mk-j, inc = mk, offset =
- ind + j*(mk + 1) - 1)
- blas.scal(2.0, x, offset = ind + mk*(j-1) + j, n = mk-j)
- ind += mk**2
+def sgemv(A, x, y, dims, trans = 'N', alpha = 1.0, beta = 0.0, n = None,
+ offsetA = 0, offsetx = 0, offsety = 0):
+ """
+ Matrix-vector multiplication.
+
+ A is a matrix or spmatrix of size (m, n) where
+
+ N = dims['l'] + sum(dims['q']) + sum( k**2 for k in dims['s'] )
+
+ representing a method from R^n to S.
+
+ If trans is 'N':
+
+ y := alpha*A*x + beta * y (trans = 'N').
+
+ x is a vector of length n. y is a vector of length N.
+
+ If trans is 'T':
+
+ y := alpha*A'*x + beta * y (trans = 'T').
+
+ x is a vector of length N. y is a vector of length n.
+
+ The 's' components in S are stored in unpacked 'L' storage.
+ """
+
+ m = dims['l'] + sum(dims['q']) + sum([ k**2 for k in dims['s'] ])
+ if n is None: n = A.size[1]
+ if trans == 'T' and alpha: trisc(x, dims, offsetx)
base.gemv(A, x, y, trans = trans, alpha = alpha, beta = beta, m = m,
n = n, offsetA = offsetA, offsetx = offsetx, offsety = offsety)
-
- if trans == 'T' and alpha:
- ind = offsetx + dims['l'] + sum(dims['q'])
- for mk in dims['s']:
- # Scale strict lower triangular part of x by 0.5.
- for j in xrange(1, mk):
- blas.scal(0.5, x, offset = ind + mk*(j-1) + j, n = mk-j)
- ind += mk**2
+ if trans == 'T' and alpha: triusc(x, dims, offsetx)
def jdot(x, y, n = None, offsetx = 0, offsety = 0):
-
- # Returns x' * J * y, where J = [1, 0; 0, -I].
+ """
+ Returns x' * J * y, where J = [1, 0; 0, -I].
+ """
if n is None:
- if len(x) != len(y): raise ValueError, "x and y must have the "\
+ if len(x) != len(y): raise ValueError, "x and y must have the "\
"same length"
- n = len(x)
+ n = len(x)
return x[offsetx] * y[offsety] - blas.dot(x, y, n = n-1,
offsetx = offsetx + 1, offsety = offsety + 1)
def jnrm2(x, n = None, offset = 0):
-
- # Returns sqrt(x' * J * x) where J = [1, 0; 0, -I], for a vector
- # x in a second order cone.
+ """
+ Returns sqrt(x' * J * x) where J = [1, 0; 0, -I], for a vector
+ x in a second order cone.
+ """
if n is None: n = len(x)
a = blas.nrm2(x, n = n-1, offset = offset+1)
return math.sqrt(x[offset] - a) * math.sqrt(x[offset] + a)
-def symm(x, n, offset = 0):
-
- # Fills in the upper triangular part of the symmetric matrix stored in
- # x[offset : offset+n*n] using 'L' storage.
+if use_C:
+ symm = misc_solvers.symm
+else:
+ def symm(x, n, offset = 0):
+ """
+ Converts lower triangular matrix to symmetric.
+ Fills in the upper triangular part of the symmetric matrix stored in
+ x[offset : offset+n*n] using 'L' storage.
+ """
if n <= 1: return
for i in xrange(n-1):
@@ -785,10 +870,14 @@ def symm(x, n, offset = 0):
offset + (i+1)*(n+1) - 1, incy = n, n = n-i-1)
-def sprod(x, y, dims, mnl = 0, diag = 'N'):
-
- # The product x := (y o x). If diag is 'D', the 's' part of y is
- # diagonal and only the diagonal is stored.
+if use_C:
+ sprod = misc_solvers.sprod
+else:
+ def sprod(x, y, dims, mnl = 0, diag = 'N'):
+ """
+ The product x := (y o x). If diag is 'D', the 's' part of y is
+ diagonal and only the diagonal is stored.
+ """
# For the nonlinear and 'l' blocks:
@@ -800,9 +889,9 @@ def sprod(x, y, dims, mnl = 0, diag = 'N'):
# For 'q' blocks:
#
- # [ lo l1' ]
+ # [ l0 l1' ]
# yk o xk = [ ] * xk
- # [ lo l0*I ]
+ # [ l1 l0*I ]
#
# where yk = (l0, l1).
@@ -852,9 +941,10 @@ def sprod(x, y, dims, mnl = 0, diag = 'N'):
def ssqr(x, y, dims, mnl = 0):
-
- # The product x := y o y. The 's' components of y are diagonal and
- # only the diagonals of x and y are stored.
+ """
+ The product x := y o y. The 's' components of y are diagonal and
+ only the diagonals of x and y are stored.
+ """
blas.copy(y, x)
blas.tbmv(y, x, n = mnl + dims['l'], k = 0, ldA = 1)
@@ -867,10 +957,14 @@ def ssqr(x, y, dims, mnl = 0):
offsetx = ind)
-def sinv(x, y, dims, mnl = 0):
-
- # The inverse product x := (y o\ x), when the 's' components of y are
- # diagonal.
+if use_C:
+ sinv = misc_solvers.sinv
+else:
+ def sinv(x, y, dims, mnl = 0):
+ """
+ The inverse product x := (y o\ x), when the 's' components of y are
+ diagonal.
+ """
# For the nonlinear and 'l' blocks:
#
@@ -916,16 +1010,20 @@ def sinv(x, y, dims, mnl = 0):
ind2 += m
-def max_step(x, dims, mnl = 0, sigma = None):
-
- # Returns min {t | x + t*e >= 0}, where e is defined as follows
- #
- # - For the nonlinear and 'l' blocks: e is the vector of ones.
- # - For the 'q' blocks: e is the first unit vector.
- # - For the 's' blocks: e is the identity matrix.
- #
- # When called with the argument sigma, also returns the eigenvalues
- # (in sigma) and the eigenvectors (in x) of the 's' components of x.
+if use_C:
+ max_step = misc_solvers.max_step
+else:
+ def max_step(x, dims, mnl = 0, sigma = None):
+ """
+ Returns min {t | x + t*e >= 0}, where e is defined as follows
+
+ - For the nonlinear and 'l' blocks: e is the vector of ones.
+ - For the 'q' blocks: e is the first unit vector.
+ - For the 's' blocks: e is the identity matrix.
+
+ When called with the argument sigma, also returns the eigenvalues
+ (in sigma) and the eigenvectors (in x) of the 's' components of x.
+ """
t = []
ind = mnl + dims['l']
@@ -953,25 +1051,26 @@ def max_step(x, dims, mnl = 0, sigma = None):
def kkt_ldl(G, dims, A, mnl = 0):
-
- # Solution of KKT equations by a dense LDL factorization of the
- # 3 x 3 system.
- #
- # Returns a function that (1) computes the LDL factorization of
- #
- # [ H A' GG'*W^{-1} ]
- # [ A 0 0 ],
- # [ W^{-T}*GG 0 -I ]
- #
- # given H, Df, W, where GG = [Df; G], and (2) returns a function for
- # solving
- #
- # [ H A' GG' ] [ ux ] [ bx ]
- # [ A 0 0 ] * [ uy ] = [ by ].
- # [ GG 0 -W'*W ] [ uz ] [ bz ]
- #
- # H is n x n, A is p x n, Df is mnl x n, G is N x n where
- # N = dims['l'] + sum(dims['q']) + sum( k**2 for k in dims['s'] ).
+ """
+ Solution of KKT equations by a dense LDL factorization of the
+ 3 x 3 system.
+
+ Returns a function that (1) computes the LDL factorization of
+
+ [ H A' GG'*W^{-1} ]
+ [ A 0 0 ],
+ [ W^{-T}*GG 0 -I ]
+
+ given H, Df, W, where GG = [Df; G], and (2) returns a function for
+ solving
+
+ [ H A' GG' ] [ ux ] [ bx ]
+ [ A 0 0 ] * [ uy ] = [ by ].
+ [ GG 0 -W'*W ] [ uz ] [ bz ]
+
+ H is n x n, A is p x n, Df is mnl x n, G is N x n where
+ N = dims['l'] + sum(dims['q']) + sum( k**2 for k in dims['s'] ).
+ """
p, n = A.size
ldK = n + p + mnl + dims['l'] + sum(dims['q']) + sum([ k*(k+1)/2 for k
@@ -1022,25 +1121,26 @@ def kkt_ldl(G, dims, A, mnl = 0):
def kkt_ldl2(G, dims, A, mnl = 0):
-
- # Solution of KKT equations by a dense LDL factorization of the 2 x 2
- # system.
- #
- # Returns a function that (1) computes the LDL factorization of
- #
- # [ H + GG' * W^{-1} * W^{-T} * GG A' ]
- # [ ]
- # [ A 0 ]
- #
- # given H, Df, W, where GG = [Df; G], and (2) returns a function for
- # solving
- #
- # [ H A' GG' ] [ ux ] [ bx ]
- # [ A 0 0 ] * [ uy ] = [ by ].
- # [ GG 0 -W'*W ] [ uz ] [ bz ]
- #
- # H is n x n, A is p x n, Df is mnl x n, G is N x n where
- # N = dims['l'] + sum(dims['q']) + sum( k**2 for k in dims['s'] ).
+ """
+ Solution of KKT equations by a dense LDL factorization of the 2 x 2
+ system.
+
+ Returns a function that (1) computes the LDL factorization of
+
+ [ H + GG' * W^{-1} * W^{-T} * GG A' ]
+ [ ]
+ [ A 0 ]
+
+ given H, Df, W, where GG = [Df; G], and (2) returns a function for
+ solving
+
+ [ H A' GG' ] [ ux ] [ bx ]
+ [ A 0 0 ] * [ uy ] = [ by ].
+ [ GG 0 -W'*W ] [ uz ] [ bz ]
+
+ H is n x n, A is p x n, Df is mnl x n, G is N x n where
+ N = dims['l'] + sum(dims['q']) + sum( k**2 for k in dims['s'] ).
+ """
p, n = A.size
ldK = n + p
@@ -1106,28 +1206,29 @@ def kkt_ldl2(G, dims, A, mnl = 0):
def kkt_chol(G, dims, A, mnl = 0):
-
- # Solution of KKT equations by reduction to a 2 x 2 system, a QR
- # factorization to eliminate the equality constraints, and a dense
- # Cholesky factorization of order n-p.
- #
- # Computes the QR factorization
- #
- # A' = [Q1, Q2] * [R; 0]
- #
- # and returns a function that (1) computes the Cholesky factorization
- #
- # Q_2^T * (H + GG^T * W^{-1} * W^{-T} * GG) * Q2 = L * L^T,
- #
- # given H, Df, W, where GG = [Df; G], and (2) returns a function for
- # solving
- #
- # [ H A' GG' ] [ ux ] [ bx ]
- # [ A 0 0 ] * [ uy ] = [ by ].
- # [ GG 0 -W'*W ] [ uz ] [ bz ]
- #
- # H is n x n, A is p x n, Df is mnl x n, G is N x n where
- # N = dims['l'] + sum(dims['q']) + sum( k**2 for k in dims['s'] ).
+ """
+ Solution of KKT equations by reduction to a 2 x 2 system, a QR
+ factorization to eliminate the equality constraints, and a dense
+ Cholesky factorization of order n-p.
+
+ Computes the QR factorization
+
+ A' = [Q1, Q2] * [R; 0]
+
+ and returns a function that (1) computes the Cholesky factorization
+
+ Q_2^T * (H + GG^T * W^{-1} * W^{-T} * GG) * Q2 = L * L^T,
+
+ given H, Df, W, where GG = [Df; G], and (2) returns a function for
+ solving
+
+ [ H A' GG' ] [ ux ] [ bx ]
+ [ A 0 0 ] * [ uy ] = [ by ].
+ [ GG 0 -W'*W ] [ uz ] [ bz ]
+
+ H is n x n, A is p x n, Df is mnl x n, G is N x n where
+ N = dims['l'] + sum(dims['q']) + sum( k**2 for k in dims['s'] ).
+ """
p, n = A.size
cdim = mnl + dims['l'] + sum(dims['q']) + sum([ k**2 for k in
@@ -1144,7 +1245,6 @@ def kkt_chol(G, dims, A, mnl = 0):
lapack.geqrf(QA, tauA)
Gs = matrix(0.0, (cdim, n))
- g = matrix(0.0, (cdim, 1))
K = matrix(0.0, (n,n))
bzp = matrix(0.0, (cdim_pckd, 1))
yy = matrix(0.0, (p,1))
@@ -1164,13 +1264,11 @@ def kkt_chol(G, dims, A, mnl = 0):
Gs[:mnl, :] = Df
Gs[mnl:, :] = G
scale(Gs, W, trans = 'T', inverse = 'I')
- for k in xrange(n):
- g[:] = Gs[:, k]
- pack(g, Gs, dims, mnl, offsety = k*Gs.size[0])
+ pack2(Gs, dims, mnl)
# K = [Q1, Q2]' * (H + Gs' * Gs) * [Q1, Q2].
- if H is not None: K[:,:] = H
- blas.syrk(Gs, K, beta = 1.0, k = cdim_pckd, trans = 'T')
+ blas.syrk(Gs, K, k = cdim_pckd, trans = 'T')
+ if H is not None: K[:,:] += H
symm(K, n)
lapack.ormqr(QA, tauA, K, side = 'L', trans = 'T')
lapack.ormqr(QA, tauA, K, side = 'R')
@@ -1247,32 +1345,33 @@ def kkt_chol(G, dims, A, mnl = 0):
def kkt_chol2(G, dims, A, mnl = 0):
-
- # Solution of KKT equations by reduction to a 2 x 2 system, a sparse
- # or dense Cholesky factorization of order n to eliminate the 1,1
- # block, and a sparse or dense Cholesky factorization of order p.
- # Implemented only for problems with no second-order or semidefinite
- # cone constraints.
- #
- # Returns a function that (1) computes Cholesky factorizations of
- # the matrices
- #
- # S = H + GG' * W^{-1} * W^{-T} * GG,
- # K = A * S^{-1} *A'
- #
- # or (if K is singular in the first call to the function), the matrices
- #
- # S = H + GG' * W^{-1} * W^{-T} * GG + A' * A,
- # K = A * S^{-1} * A',
- #
- # given H, Df, W, where GG = [Df; G], and (2) returns a function for
- # solving
- #
- # [ H A' GG' ] [ ux ] [ bx ]
- # [ A 0 0 ] * [ uy ] = [ by ].
- # [ GG 0 -W'*W ] [ uz ] [ bz ]
- #
- # H is n x n, A is p x n, Df is mnl x n, G is dims['l'] x n.
+ """
+ Solution of KKT equations by reduction to a 2 x 2 system, a sparse
+ or dense Cholesky factorization of order n to eliminate the 1,1
+ block, and a sparse or dense Cholesky factorization of order p.
+ Implemented only for problems with no second-order or semidefinite
+ cone constraints.
+
+ Returns a function that (1) computes Cholesky factorizations of
+ the matrices
+
+ S = H + GG' * W^{-1} * W^{-T} * GG,
+ K = A * S^{-1} *A'
+
+ or (if K is singular in the first call to the function), the matrices
+
+ S = H + GG' * W^{-1} * W^{-T} * GG + A' * A,
+ K = A * S^{-1} * A',
+
+ given H, Df, W, where GG = [Df; G], and (2) returns a function for
+ solving
+
+ [ H A' GG' ] [ ux ] [ bx ]
+ [ A 0 0 ] * [ uy ] = [ by ].
+ [ GG 0 -W'*W ] [ uz ] [ bz ]
+
+ H is n x n, A is p x n, Df is mnl x n, G is dims['l'] x n.
+ """
if dims['q'] or dims['s']:
raise ValueError, "kktsolver option 'kkt_chol2' is implemented "\
@@ -1464,28 +1563,29 @@ def kkt_chol2(G, dims, A, mnl = 0):
def kkt_qr(G, dims, A):
-
- # Solution of KKT equations with zero 1,1 block, by eliminating the
- # equality constraints via a QR factorization, and solving the
- # reduced KKT system by another QR factorization.
- #
- # Computes the QR factorization
- #
- # A' = [Q1, Q2] * [R1; 0]
- #
- # and returns a function that (1) computes the QR factorization
- #
- # W^{-T} * G * Q2 = Q3 * R3
- #
- # (with columns of W^{-T}*G in packed storage), and (2) returns a
- # function for # solving
- #
- # [ 0 A' G' ] [ ux ] [ bx ]
- # [ A 0 0 ] * [ uy ] = [ by ].
- # [ G 0 -W'*W ] [ uz ] [ bz ]
- #
- # A is p x n and G is N x n where N = dims['l'] + sum(dims['q']) +
- # sum( k**2 for k in dims['s'] ).
+ """
+ Solution of KKT equations with zero 1,1 block, by eliminating the
+ equality constraints via a QR factorization, and solving the
+ reduced KKT system by another QR factorization.
+
+ Computes the QR factorization
+
+ A' = [Q1, Q2] * [R1; 0]
+
+ and returns a function that (1) computes the QR factorization
+
+ W^{-T} * G * Q2 = Q3 * R3
+
+ (with columns of W^{-T}*G in packed storage), and (2) returns a
+ function for # solving
+
+ [ 0 A' G' ] [ ux ] [ bx ]
+ [ A 0 0 ] * [ uy ] = [ by ].
+ [ G 0 -W'*W ] [ uz ] [ bz ]
+
+ A is p x n and G is N x n where N = dims['l'] + sum(dims['q']) +
+ sum( k**2 for k in dims['s'] ).
+ """
p, n = A.size
cdim = dims['l'] + sum(dims['q']) + sum([ k**2 for k in dims['s'] ])
@@ -1502,7 +1602,6 @@ def kkt_qr(G, dims, A):
Gs = matrix(0.0, (cdim, n))
tauG = matrix(0.0, (n-p,1))
- g = matrix(0.0, (cdim, 1))
u = matrix(0.0, (cdim_pckd, 1))
vv = matrix(0.0, (n,1))
w = matrix(0.0, (cdim_pckd, 1))
@@ -1512,9 +1611,7 @@ def kkt_qr(G, dims, A):
# Gs = W^{-T}*G, in packed storage.
Gs[:,:] = G
scale(Gs, W, trans = 'T', inverse = 'I')
- for k in xrange(n):
- g[:] = Gs[:, k]
- pack(g, Gs, dims, offsety = k*Gs.size[0])
+ pack2(Gs, dims)
# Gs := [ Gs1, Gs2 ]
# = Gs * [ Q1, Q2 ]
@@ -1522,7 +1619,7 @@ def kkt_qr(G, dims, A):
# QR factorization Gs2 := [ Q3, Q4 ] * [ R3; 0 ]
lapack.geqrf(Gs, tauG, n = n-p, m = cdim_pckd, offsetA =
- Gs.size[0]*p)
+ Gs.size[0]*p)
def solve(x, y, z):
diff --git a/src/python/modeling.py b/src/python/modeling.py
index 788b36f..aa5dc88 100644
--- a/src/python/modeling.py
+++ b/src/python/modeling.py
@@ -7,7 +7,7 @@ piecewise-linear objective and constraint functions.
# Copyright 2004-2008 J. Dahl and L. Vandenberghe.
#
-# This file is part of CVXOPT version 0.9.3.
+# This file is part of CVXOPT version 1.0.
#
# CVXOPT is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
diff --git a/src/python/mosek.py b/src/python/mosek.py
index 7fdc4ca..7592b7e 100644
--- a/src/python/mosek.py
+++ b/src/python/mosek.py
@@ -4,7 +4,7 @@ CVXOPT interface for MOSEK 5.0
# Copyright 2004-2008 J. Dahl and L. Vandenberghe.
#
-# This file is part of CVXOPT version 0.9.3.
+# This file is part of CVXOPT version 1.0.
#
# CVXOPT is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
@@ -32,7 +32,7 @@ inf = 0.0 # numeric value doesn't matter
options = {}
-def solvelp(c, G, h, A=None, b=None):
+def lp(c, G, h, A=None, b=None):
"""
Solves a pair of primal and dual LPs
@@ -43,7 +43,7 @@ def solvelp(c, G, h, A=None, b=None):
using MOSEK 5.0.
- (solsta, x, y, z) = solvelp(c, G, h, A=None, b=None).
+ (solsta, x, y, z) = lp(c, G, h, A=None, b=None).
Input arguments
@@ -185,7 +185,6 @@ def solvelp(c, G, h, A=None, b=None):
else:
return (solsta[0], x, z, y)
-
def conelp(c, G, h, dims = None):
"""
Solves a pair of primal and dual SOCPs
@@ -244,7 +243,7 @@ def conelp(c, G, h, dims = None):
"""
if dims is None:
- (solsta, x, y, z) = solvelp(c, G, h)
+ (solsta, x, y, z) = lp(c, G, h)
return (solsta, x, z, None)
try:
@@ -538,7 +537,7 @@ def socp(c, Gl = None, hl = None, Gq = None, hq = None):
return (solsta[0], x, zl, zq)
-def solveqp(P, q, G=None, h=None, A=None, b=None):
+def qp(P, q, G=None, h=None, A=None, b=None):
"""
Solves a quadratic program
@@ -548,7 +547,7 @@ def solveqp(P, q, G=None, h=None, A=None, b=None):
using MOSEK 5.0.
- (solsta, x, z, y) = solveqp(P, q, G=None, h=None, A=None, b=None)
+ (solsta, x, z, y) = qp(P, q, G=None, h=None, A=None, b=None)
Return values
@@ -702,5 +701,172 @@ def solveqp(P, q, G=None, h=None, A=None, b=None):
return (solsta[0], x, z, y)
+def ilp(c, G, h, A=None, b=None, I=None):
+ """
+ Solves the mixed integer LP
+
+ minimize c'*x
+ subject to G*x + s = h
+ A*x = b
+ s >= 0
+ xi integer, forall i in I
+
+ using MOSEK 5.0.
+
+ (solsta, x) = ilp(c, G, h, A=None, b=None, I=None).
+
+ Input arguments
+
+ G is m x n, h is m x 1, A is p x n, b is p x 1. G and A must be
+ dense or sparse 'd' matrices. h and b are dense 'd' matrices
+ with one column. The default values for A and b are empty
+ matrices with zero rows.
+
+ I is a Python set with indices of integer elements of x. By
+ default all elements in x are constrained to be integer, i.e.,
+ the default value of I is I = set(range(n))
+
+ Dual variables are not returned for MOSEK.
+
+ Return values
+
+ solsta the solution status.
+ If solsta is solsta.integer_optimal,
+ then x contains the solution.
+ If solsta is solsta.unknown,
+ then x is None
+
+ Other return values for solsta include:
+ mosek.solsta.near_integer_optimal
+ in which case the x value may not be well-defined,
+ c.f., the MOSEK documentation.
+
+ x the solution
+
+ Options are passed to MOSEK solvers via the mosek.options
+ dictionary, e.g., the following turns off output from
+ the MOSEK solvers
+
+ >>> mosek.options = {iparam.log:0}
+
+ see chapter 14 of the MOSEK Python API manual.
+ """
+ if type(c) is not matrix or c.typecode != 'd' or c.size[1] != 1:
+ raise TypeError, "'c' must be a dense column matrix"
+ n = c.size[0]
+ if n < 1: raise ValueError, "number of variables must be at least 1"
+
+ if (type(G) is not matrix and type(G) is not spmatrix) or \
+ G.typecode != 'd' or G.size[1] != n:
+ raise TypeError, "'G' must be a dense or sparse 'd' matrix "\
+ "with %d columns" %n
+ m = G.size[0]
+ if m is 0: raise ValueError, "m cannot be 0"
+
+ if type(h) is not matrix or h.typecode != 'd' or h.size != (m,1):
+ raise TypeError, "'h' must be a 'd' matrix of size (%d,1)" %m
+
+ if A is None: A = spmatrix([], [], [], (0,n), 'd')
+ if (type(A) is not matrix and type(A) is not spmatrix) or \
+ A.typecode != 'd' or A.size[1] != n:
+ raise TypeError, "'A' must be a dense or sparse 'd' matrix "\
+ "with %d columns" %n
+ p = A.size[0]
+ if b is None: b = matrix(0.0, (0,1))
+ if type(b) is not matrix or b.typecode != 'd' or b.size != (p,1):
+ raise TypeError, "'b' must be a dense matrix of size (%d,1)" %p
+
+ c = array(c)
+
+ if I is None: I = set(range(n))
+
+ if type(I) is not set:
+ raise TypeError, "invalid argument for integer index set"
+
+ for i in I:
+ if type(i) is not int:
+ raise TypeError, "invalid integer index set I"
+
+ if len(I) > 0 and min(I) < 0: raise IndexError, \
+ "negative element in integer index set I"
+ if len(I) > 0 and max(I) > n-1: raise IndexError, \
+ "maximum element in in integer index set I is larger than n-1"
+
+ bkc = m*[ msk.boundkey.up ] + p*[ msk.boundkey.fx ]
+ blc = m*[ -inf ] + [ bi for bi in b ]
+ buc = matrix([h, b])
+
+ bkx = n*[msk.boundkey.fr]
+ blx = n*[ -inf ]
+ bux = n*[ +inf ]
+
+ colptr, asub, acof = sparse([G,A]).CCS
+ aptrb, aptre = colptr[:-1], colptr[1:]
+
+ mskhandle = msk.mosek ()
+ env = mskhandle.Env ()
+ env.set_Stream (msk.streamtype.log, streamprinter)
+
+ env.init ()
+ task = env.Task(0,0)
+ task.set_Stream (msk.streamtype.log, streamprinter)
+
+ # set MOSEK options
+ for (param, val) in options.items():
+ if str(param)[:6] == "iparam":
+ task.putintparam(param, val)
+ elif str(param)[:6] == "dparam":
+ task.putdouparam(param, val)
+ elif str(param)[:6] == "sparam":
+ task.putstrparam(param, val)
+ else:
+ raise ValueError, "invalid MOSEK parameter: "+str(param)
+
+ task.inputdata (m+p, # number of constraints
+ n, # number of variables
+ array(c), # linear objective coefficients
+ 0.0, # objective fixed value
+ array(aptrb),
+ array(aptre),
+ array(asub),
+ array(acof),
+ array(bkc),
+ array(blc),
+ array(buc),
+ array(bkx),
+ array(blx),
+ array(bux))
+
+ task.putobjsense(msk.objsense.minimize)
+
+ # Define integer variables
+ if len(I) > 0:
+ task.putvartypelist(array(I),
+ array(len(I)*[ msk.variabletype.type_int ]))
+
+ task.putintparam (msk.iparam.mio_mode, msk.miomode.satisfied)
+
+ task.optimize()
+
+ task.solutionsummary (msk.streamtype.msg);
+
+ prosta, solsta = list(), list()
+ if len(I) > 0:
+ task.getsolutionstatus(msk.soltype.itg, prosta, solsta)
+ else:
+ task.getsolutionstatus(msk.soltype.bas, prosta, solsta)
+
+ x = zeros(n, Float)
+ if len(I) > 0:
+ task.getsolutionslice(msk.soltype.itg, msk.solitem.xx, 0, n, x)
+ else:
+ task.getsolutionslice(msk.soltype.bas, msk.solitem.xx, 0, n, x)
+ x = matrix(x)
+
+ if (solsta[0] is msk.solsta.unknown):
+ return (solsta[0], None)
+ else:
+ return (solsta[0], x)
+
diff --git a/src/python/printing.py b/src/python/printing.py
index 3e7ab68..f5f00a0 100644
--- a/src/python/printing.py
+++ b/src/python/printing.py
@@ -1,6 +1,6 @@
# Copyright 2004-2008 J. Dahl and L. Vandenberghe.
#
-# This file is part of CVXOPT version 0.9.3.
+# This file is part of CVXOPT version 1.0.
#
# CVXOPT is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
@@ -99,7 +99,6 @@ def spmatrix_str_default(X):
if height < 0: height = maxint
if width*height is 0: return ""
- if len(X) is 0: return ""
rlist = range(0,min(m,height))
clist = range(0,min(n,width))
diff --git a/src/python/solvers.py b/src/python/solvers.py
index b647376..757984f 100644
--- a/src/python/solvers.py
+++ b/src/python/solvers.py
@@ -15,7 +15,7 @@ options: dictionary with customizable algorithm parameters.
# Copyright 2004-2008 J. Dahl and L. Vandenberghe.
#
-# This file is part of CVXOPT version 0.9.3.
+# This file is part of CVXOPT version 1.0.
#
# CVXOPT is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
diff --git a/src/setup.py b/src/setup.py
index 67c3483..d6e2e24 100644
--- a/src/setup.py
+++ b/src/setup.py
@@ -148,11 +148,16 @@ amd = Extension('amd',
sources = [ 'C/amd.c' ] + [ 'C/SuiteSparse/AMD/Source/' + s for s in
listdir('C/SuiteSparse/AMD/Source') if s[-2:] == '.c' ])
-extmods += [base, blas, lapack, umfpack, cholmod, amd]
+misc_solvers = Extension('misc_solvers', libraries = ['lapack', 'blas'],
+ library_dirs = [ ATLAS_LIB_DIR ],
+ define_macros = MACROS,
+ sources = ['C/misc_solvers.c'] )
+
+extmods += [base, blas, lapack, umfpack, cholmod, amd, misc_solvers]
setup (name = 'cvxopt',
description = 'Convex optimization package',
- version = '0.9.3',
+ version = '1.0',
long_description = '''
CVXOPT is a free software package for convex optimization based on the
Python programming language. It can be used with the interactive Python
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/cvxopt.git
More information about the debian-science-commits
mailing list