SOCP.cpp

/*
This file is part of the BIAS library (Basic ImageAlgorithmS).

Copyright (C) 2003-2009    (see file CONTACT for details)
  Multimediale Systeme der Informationsverarbeitung
  Institut fuer Informatik
  Christian-Albrechts-Universitaet Kiel


BIAS is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version.

BIAS 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 Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License
along with BIAS; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

#include "SOCP.hh"

extern "C" {
#include "socp/socp.h"
}

using namespace BIAS;

SOCP::SOCP()
{
  MaxIters_ = 100;
  AbsTolerance_ = 1e-6;
  RelTolerance_ = 1e-4;
  TargetValue_ = 0.0;
  Nu_ = 10.0;
}

SOCP::~SOCP()
{
}

int SOCP::Compute(int L, int m, const std::vector<int> &N,
                  const Vector<double> &f,
                  const std::vector< Matrix<double> > &A,
                  const std::vector< Vector<double> > &b,
                  const std::vector< Vector<double> > &c,
                  const std::vector<double> &d,
                  Vector<double> &x, std::vector< Vector<double> > &z,
                  std::vector<double> &w)
{
#ifdef BIAS_DEBUG
  // Check sizes of input matrices and vectors
  BIASASSERT((int)f.Size() == m);
  BIASASSERT((int)x.Size() == m);
  BIASASSERT((int)N.size() == L);
  BIASASSERT((int)A.size() == L);
  BIASASSERT((int)b.size() == L);
  BIASASSERT((int)c.size() == L);
  BIASASSERT((int)d.size() == L);
  BIASASSERT((int)z.size() == L);
  BIASASSERT((int)w.size() == L);
  for (int i = 0; i < L; i++) {
    BIASASSERT((int)A[i].GetRows() == N[i] && (int)A[i].GetCols() == m);
    BIASASSERT((int)b[i].Size() == N[i]);
    BIASASSERT((int)c[i].Size() == m);
    BIASASSERT((int)z[i].Size() == N[i]);
  }
#endif

  int n = 0;
  int *in_N = new int[L];
  for (int i = 0; i < L; i++)
  {
    n += N[i];
    in_N[i] = N[i] + 1;
  }
  double *in_f = new double[m];
  double *in_Ac = new double[(n+L)*m];
  double *in_bd = new double[n+L];
  double *inout_x = new double[m];
  double *inout_zw = new double[n+L];
  double *ptr = in_Ac;
  for (int i = 0; i < L; i++) {
    for (int col = 0; col < m; col++) {
      for (int row = 0; row < N[i]; row++)
        *(ptr++) = A[i][row][col];
      *(ptr++) = c[i][col];
    }
  }
  ptr = in_bd;
  for (int i = 0; i < L; i++) {
    for (int row = 0; row < N[i]; row++)
      *(ptr++) = b[i][row];
    *(ptr++) = d[i];
  }
  for (int j = 0; j < m; j++) {
    in_f[j] = f[j];
    inout_x[j] = x[j];
  }
  ptr = inout_zw;
  for (int i = 0; i < L; i++) {
    for (int row = 0; row < N[i]; row++)
      *(ptr++) = z[i][row];
    *(ptr++) = w[i];
  }
  
  // Compute workspace size
  int dblSize, intSize, mhist, nhist;
  socp_getwork(L, in_N, m, MaxIters_, 0, &mhist, &nhist, &dblSize, &intSize);
  double* dblWorkspace = new double[dblSize];
  int* intWorkspace = new int[intSize];
  
  int iter = MaxIters_, info = 0;
  int res = socp(L, in_N, m, in_f, in_Ac, in_bd, inout_x, inout_zw,
                 AbsTolerance_, RelTolerance_, TargetValue_, &iter, Nu_,
                 &info, 0, (double*)0, dblWorkspace, intWorkspace);

  // Write results back
  for (int j = 0; j < m; j++)
    x[j] = inout_x[j];
  ptr = inout_zw;
  for (int i = 0; i < L; i++) {
    for (int row = 0; row < N[i]; row++)
      z[i][row] = *(ptr++);
    w[i] = *(ptr++);
  }
  
  // Release memory
  delete[] in_N;
  delete[] in_f;
  delete[] in_Ac;
  delete[] in_bd;
  delete[] inout_x;
  delete[] inout_zw;
  delete[] dblWorkspace;
  delete[] intWorkspace;
  
  return (res != 0) ? -1 : 0;
}

int SOCP::Compute(int m, int n, const Vector<double> &f,
                  const Matrix<double> &A, const Vector<double> &b,
                  const Vector<double> &c, double d,
                  Vector<double> &x, Vector<double> &z, double &w)
{
#ifdef BIAS_DEBUG
  // Check sizes of input matrix and vectors
  BIASASSERT((int)f.Size() == m);
  BIASASSERT((int)x.Size() == m);
  BIASASSERT((int)A.GetRows() == n && (int)A.GetCols() == m);
  BIASASSERT((int)b.Size() == n);
  BIASASSERT((int)c.Size() == m);
  BIASASSERT((int)z.Size() == n)
#endif

  const int L = 1;
  int in_N = n + 1;
  double *in_f = new double[m];
  double *in_Ac = new double[(n+L)*m];
  double *in_bd = new double[n+L];
  double *inout_x = new double[m];
  double *inout_zw = new double[n+L];
  double *ptr = in_Ac;
  for (int col = 0; col < m; col++) {
    for (int row = 0; row < n; row++)
      *(ptr++) = A[row][col];
    *(ptr++) = c[col];
  }
  ptr = in_bd;
  for (int row = 0; row < n; row++)
    *(ptr++) = b[row];
  *(ptr++) = d;
  for (int j = 0; j < m; j++) {
    in_f[j] = f[j];
    inout_x[j] = x[j];
  }
  ptr = inout_zw;
  for (int row = 0; row < n; row++)
    *(ptr++) = z[row];
  *(ptr++) = w;

  // Compute workspace size
  int dblSize, intSize, mhist, nhist;
  socp_getwork(L, &in_N, m, MaxIters_, 0, &mhist, &nhist, &dblSize, &intSize);
  double* dblWorkspace = new double[dblSize];
  int* intWorkspace = new int[intSize];

  int iter = MaxIters_, info = 0;
  int res = socp(L, &in_N, m, in_f, in_Ac, in_bd, inout_x, inout_zw,
                 AbsTolerance_, RelTolerance_, TargetValue_, &iter, Nu_,
                 &info, 0, (double*)0, dblWorkspace, intWorkspace);

  // Write results back
  for (int j = 0; j < m; j++)
    x[j] = inout_x[j];
  ptr = inout_zw;
  for (int row = 0; row < n; row++)
    z[row] = *(ptr++);
  w = *(ptr++);

  // Release memory
  delete[] in_f;
  delete[] in_Ac;
  delete[] in_bd;
  delete[] inout_x;
  delete[] inout_zw;
  delete[] dblWorkspace;
  delete[] intWorkspace;

  return (res != 0) ? -1 : 0;
}