Logo Search packages:      
Sourcecode: semidef-oct version File versions  Download package

sp.cc

// semidef package, (c) 1994 Lieven Vandenberghe and Stephen Boyd.
//
// Implemented Sept. 1998 for Octave by A. S. Hodel a.s.hodel@eng.auburn.edu

//#define DEBUG
//#define DEBUG_INPUT_PARAMETERS

#include <octave/config.h>
#include <iostream.h>
#include <math.h>
#include <string.h>
#include <strstream.h>
#include <utils.h>

#include <octave/defun-dld.h>
#include <octave/error.h>
#include <octave/f77-fcn.h>
#include <octave/gripes.h>
#include <octave/help.h>
#include <octave/oct-obj.h>
#include <octave/oct-map.h>
#include <octave/ov.h>
#include <octave/pager.h>
#include <octave/symtab.h>
#include <octave/variables.h>

#if defined(DEBUG) || defined(DEBUG_INPUT_PARAMETERS)
#include <octave/pr-output.h>
#endif

#include "ocst.h"
#include "sp.h"
#ifdef userusage
#include <sys/resource.h>
#else
#include <sys/times.h>
#endif

DEFUN_DLD(sp, args, nargout,
"\n\
 [x,Z,upper_limit,infostr,time] = ...\n\
          sp(F,blck_szs,c,x0,Z0,nu,abstol,reltol,tv,maxiters);\n\
\n\
 Solves semidefinite program\n\
\n\
  minimize    c'*x\n\
  subject to  F_0 + x_1*F_1 + ... + x_m*F_m  >= 0\n\
\n\
 and its dual\n\
\n\
  maximize    -Tr F_0 Z\n\
  subject to  Tr F_i Z = c_i, i=1,...,m\n\
              Z >= 0\n\
\n\
 Convergence criterion:\n\
 (1) maxiters is exceeded\n\
 (2) duality gap is less than abstol\n\
 (3) primal and dual objective are both positive and\n\
     duality gap is less than (reltol * dual objective)\n\
     or primal and dual objective are both negative and\n\
     duality gap is less than (reltol * minus the primal objective)\n\
 (4) reltol is negative and primal objective is less than tv\n\
 (5) reltol is negative and dual objective is greater than tv\n\
\n\
\n\
 Input arguments:\n\
 - F:        matrix with m+1 columns,\n\
             F = [ F_0^1(:)  F_1^1(:) ... F_m^1(:) ]\n\
                 [ F_0^2(:)  F_1^2(:) ... F_m^2(:) ]\n\
                 [ ...        ...     ...    ...   ]\n\
                 [ F_0^L(:)  F_1^L(:) ... F_m^L(:) ]\n\
             F_i^j is the jth diagonal block of F_i\n\
 - blck_szs: L-vector with dimensions of the diagonal blocks\n\
             the jth block of F (F_i^j has dimension blck_szs(j) )\n\
 - c:        m-vector, primal objective\n\
 - x0:       m-vector, must be strictly primal feasible\n\
 - Z0:       [ Z0^1(:); ... ; Z0^L(:) ],\n\
             must be strictly dual feasible\n\
 - nu:       >= 1.0\n\
 - abstol:   absolute tolerance\n\
 - reltol:   negative tolerance\n\
 - tv:       target value\n\
 - maxiters: max. no of iterations (maxiters >= 0)\n\
\n\
 Output arguments:\n\
 - x, Z:     last primal and dual iterates\n\
 - upper_limit[2]:    upper_limit[0] is c'*x;  upper_limit[1] is -Tr F_0*Z\n\
 - infostr:  'maxiters exceeded', 'target reached',\n\
             'unachievable target', 'relative accuracy reached',\n\
             'absolute accuracy reached', or 'error occurred in sp'\n\
 - time:     [ user time (sec.), system time (sec.), no of iters ]\n\
")
{
  octave_value_list retval;
  retval(0) = (double) 0;
  retval(1) = (double) 0;
  retval(2) = (double) 0;
  retval(3) = (double) 0;
  retval(4) = (double) 0;

  int nargin = args.length ();
  if( nargin != 10 || nargout != 5)
  {
    error("sp: %d input arguments (need 10); %d output argument (need 5)\n",
      nargin, nargout);
    print_usage("sp:");
    return retval;
  }

  // pull arguments out of passed structure.
  int errflg=0;
  Matrix FF = get_matrix(args(0),"sp: ",1,errflg);
  if(errflg) return retval;
  
  #ifdef DEBUG_INPUT_PARAMETERS
  cout << "sp: FF = " << endl;
  octave_print_internal(cout,FF,0);
  cout << endl;
  #endif

  ColumnVector dblksiz = get_vector(args(1),"sp: ",2,errflg);
  if(errflg) return retval;

  #ifdef DEBUG_INPUT_PARAMETERS
  cout << "sp: dblksiz = " << endl;
  octave_print_internal(cout, (Matrix) dblksiz,0);
  cout << endl;
  #endif

  ColumnVector cc = get_vector(args(2),"sp: ",3,errflg); // primal objective
  if(errflg) return retval;

  #ifdef DEBUG_INPUT_PARAMETERS
  cout << "sp: cc = " << endl;
  octave_print_internal(cout, (Matrix) cc,0);
  cout << endl;
  #endif

  ColumnVector x0 = get_vector(args(3),"sp: ",4,errflg); // initial vector
  if(errflg) return retval;
  
  #ifdef DEBUG_INPUT_PARAMETERS
  cout << "sp: x0 = " << endl;
  octave_print_internal(cout,(Matrix) x0, 0);
  cout << endl;
  #endif

  ColumnVector Z0 = get_vector(args(4),"sp: ",5,errflg); // dual feasible
  if(errflg) return retval;

  #ifdef DEBUG_INPUT_PARAMETERS
  cout << "sp: Z0 = " << endl;
  octave_print_internal(cout, (Matrix) Z0,0);
  cout << endl;
  #endif

  // extract nu (scalar >= 1.0)
  double nu = get_positive_real_scalar(args(5),"sp: ",6, errflg); 
  if(errflg) return retval;
  if (nu < 1.0)
  {
    error("sp: argument 6 (nu=%e) must be >= 1.0",nu);
    return retval;
  }

  #ifdef DEBUG_INPUT_PARAMETERS
  cout << "sp: nu = " << nu << endl;
  #endif

  double abstol = get_positive_real_scalar(args(6),"sp: ", 7, errflg);
  if(errflg) return retval;

  #ifdef DEBUG_INPUT_PARAMETERS
  cout << "sp: abstol = " << abstol << endl;
  #endif

  // extract reltol (negative, scalar)
  double reltol = get_real_scalar(args(7), "sp: ", 8, errflg);
  if(errflg) return retval;
  #ifdef DEBUG_INPUT_PARAMETERS
  cout << "sp: reltol = " << reltol << endl;
  #endif

  double tv = get_real_scalar(args(8), "sp: ", 9,errflg); // target value
  if(errflg) return retval;

  #ifdef DEBUG_INPUT_PARAMETERS
  cout << "sp: tv = " << tv <<endl;
  #endif

  int maxiters = (int) get_positive_real_scalar(args(9),"sp: ",10,errflg);
  if(errflg) return retval;
  #ifdef DEBUG_INPUT_PARAMETERS
  cout << "sp: maxiters = " << maxiters << endl;
  #endif

  
  //  ********************************************************
  // Process argument dimensions
  int ii;
  MArray<int> blksiz( dblksiz.length() );
  for(ii = 0 ; ii < dblksiz.length() ; ii++)
     blksiz(ii) = (int) dblksiz(ii);
  
  // Get matrix argument dimensions: FF = Fn x Fm
  int Fn = FF.rows();
  int Fm = FF.columns();
  int Fm1 = Fm-1;

  if(Fn == 0 || Fm == 0 || Fm1 == 0)
  {
    error("FF(%d,%d): bad dimensions",Fn,Fm);
    return retval;
  }

  // verify dimensions of blksiz
  for( ii = 0 ; ii < blksiz.length() ; ii++)
  {
    if( blksiz(ii) < 0 )
    {
      error("blksiz(%d) = %e < 0", ii+1, blksiz(ii+1));
      return retval;
    }
  }
  
  // compute various dimensions
  int nn,         // matrices F_i are nn x nn
      packed_size,      // number of entries of F_i in packed storage
      unpacked_size,    // number of entries of F_i in unpacked storage
      max_n ;           // max block size

  nn = packed_size = unpacked_size = max_n = 0;
  for (ii=0;  ii<blksiz.length();  ii++)
  {
    nn += blksiz(ii);
    packed_size += (blksiz(ii) * (blksiz(ii)+1))/2;
    unpacked_size += SQR(blksiz(ii));
    max_n = MAX(max_n, blksiz(ii)); 
  }

  // verify dimension of FF (unpacked storage);
  if (Fn  != unpacked_size) 
  {
    error("F(%d,%d) does not match blck_szs(expect %d rows).\n", 
      Fn, Fm1, unpacked_size);
    return retval;
  }
  
  // copy Z0 (5th input argument) to ZZ_packed (for 2nd output arg) in 
  // packed storage
  if(Z0.length() != unpacked_size)
  {
    error("sp: length(Z0)=%d, unpacked size should be %d",
      Z0.length(), unpacked_size);
    return retval;
  }

  // pack FF into FF_packed; loop over matrices F_j, store only 
  //  lower triangular portion
  // Also pack Z0 into ZZ_packed
  // of each block in FF_packed
  Matrix FF_packed(packed_size,Fm);
  ColumnVector ZZ_packed(unpacked_size);  // big just in case sp needs it
  int kk;
  #ifdef DEBUG
  for(ii = 0 ; ii < packed_size ; ii++)
    for(kk = 0 ; kk < Fm ; kk++)
      FF_packed(ii,kk) = -100*ii - kk - 101;
  cout << "initial FF_Packed:" << endl;
  octave_print_internal(cout,FF_packed,0);
  cout << endl;
  #endif

  // iterate over matrices
  int pidx = 0;   // index of next element to put into FF_packed
  int uidx = 0;   // index of next column to copy from unpacked FF

  // iterate over blocks within this matrix
  for( int blknum = 0 ; blknum < blksiz.length() ; blknum++)
  {
    // iterate over columns in this block
    for(int colnum = 0 ; colnum < blksiz(blknum) ; colnum++)
    {
      // copy the lower triangular portion of the column
      for(kk=0 ; kk < (blksiz(blknum) - colnum) ; kk++)
      {
        // iterate over matrices F_i
        for( int matnum = 0; matnum < FF.columns() ; matnum++)
          FF_packed(pidx+kk,matnum) = FF(uidx+kk+colnum,matnum);

        // also pack Z0 into ZZ_packed
        ZZ_packed(pidx+kk) = Z0(uidx+kk+colnum);
        #ifdef DEBUG
        cout << "block " << blknum << ", column " << colnum << "entry "
            << kk << ", FF_packed=" << endl;
        octave_print_internal(cout,FF_packed,0);
        cout << endl;
        #endif
      }

      // advance indices to the next column
      pidx += (blksiz(blknum) - colnum);
      uidx += blksiz(blknum);      // go to next column
    }
  }

  if(cc.length() != Fm1)
  {
    error("sp: length(cc)=%d, columns(FF)-1=%d",cc.length(), Fm1);
    return retval;
  }

  if(x0.length() != Fm1)
  {
    error("sp: length(x0)=%d, columns(FF)-1=%d",x0.length(), Fm1);
    return retval;
  }

  // upper_limit (third output arg)
  ColumnVector upper_limit(2);

  //
  // allocate work space
  //
  int lwork = (Fm1+2)*packed_size + unpacked_size + 2*nn + 
         MAX( 3*Fm, MAX( Fm+packed_size*NB, 3*max_n + max_n*(max_n+1) )); 
  ColumnVector work(lwork);
  MArray<int> iwork(Fm);

  //
  // call sp
  //
  // initialize utime, stime
  #ifdef userusage
  struct rusage stats;
  int info2 = getrusage(RUSAGE_SELF,&stats);
  double utime = - (double) stats.ru_utime.tv_sec
         - ((double) stats.ru_utime.tv_usec)/1e6;
  double stime = - (double) stats.ru_stime.tv_sec
         - ((double) stats.ru_stime.tv_usec)/1e6;
  #else
  struct tms stats;
  int info2 = times(&stats);
  double utime = - (double) stats.tms_utime / CLK_TCK;
  double stime = - (double) stats.tms_stime / CLK_TCK;
  #endif

  #ifdef DEBUG
  cout << "sp(cc): utime = " << utime << ", stime = " << stime << endl;
  #endif

  // sp
  ColumnVector xx = x0;
  int info;
  #ifdef DEBUG
  cout << "calling sp:" << endl;
  cout << "Fm1=" << Fm1 << endl
      << "# blocks = " << blksiz.length() << endl
      << "Packed FF =" << endl;
  octave_print_internal(cout,FF_packed,0);
  cout <<  endl
      << "blocksizes=";
  for(ii = 0 ; ii < blksiz.length() ; ii++)
    cout << " " << blksiz(ii) ;
  cout << endl;
  #endif
  int info3 = sp(Fm1, blksiz.length(), FF_packed.fortran_vec(),
      blksiz.fortran_vec(), cc.fortran_vec(), 
      xx.fortran_vec(), ZZ_packed.fortran_vec (), upper_limit.fortran_vec (), 
      nu, abstol, reltol, tv, &maxiters, work.fortran_vec(), lwork, 
      iwork.fortran_vec(), &info);
  #ifdef DEBUG
  cout << "back from sp: info=" << info << ", info3=" << info3 << endl;
  #endif

  // update utime, stime
  #ifdef userusage
  info2 = getrusage(RUSAGE_SELF,&stats);
  utime += (double) stats.ru_utime.tv_sec
         + ((double) stats.ru_utime.tv_usec)/1e6;
  stime += (double) stats.ru_stime.tv_sec
         + ((double) stats.ru_stime.tv_usec)/1e6;
  #else
  info2 = times(&stats);
  utime += (double) stats.tms_utime/CLK_TCK;
  stime += (double) stats.tms_stime/CLK_TCK;
  #endif
  #ifdef DEBUG
  cout << "sp(cc): utime = " << utime << ", stime = " << stime << endl;
  #endif

  //
  // unpack ZZ_packed into ZZ
  ColumnVector ZZ(unpacked_size);
  pidx = 0;   // index of next element to put into FF_packed
  uidx = 0;   // index of next column to copy from unpacked FF

  // iterate over blocks within this matrix
  for( int blknum = 0 ; blknum < blksiz.length() ; blknum++)
  {
    // iterate over columns in this block
    for(int colnum = 0 ; colnum < blksiz(blknum) ; colnum++)
    {
      // copy the lower triangular portion of the column
      for(kk=0 ; kk < (blksiz(blknum) - colnum) ; kk++)
        ZZ(uidx+kk+colnum) = ZZ_packed(pidx+kk);

      // copy the upper triangular portion of the column
      for(kk=1 ; kk < colnum ; kk++)
        ZZ(uidx-kk+colnum) = ZZ_packed(pidx+kk);

      // advance indices to the next column
      pidx += (blksiz(blknum) - colnum);
      uidx += blksiz(blknum);        // go to next column
    }
  }

  // information string output
  ostrstream infostr;
  switch (info) 
  {
  case 1:
    infostr << "maxiters exceeded";
    break;
  case 2:  infostr << "absolute tolerance reached";
           break;
  case 3:
    infostr << "relative tolerance reached";
    break;
  case 4:
    infostr << "target value reached";
    break;
  case 5:
     infostr << "target value unachievable";
     break;
  default:
    #ifdef DEBUG
    cout << "info switch: info = " << info << endl;
    #endif
    infostr << "error occurred in sp"; 
  }
  infostr << ends;

  // time
  ColumnVector time_data(3);
  time_data(0) = utime;
  time_data(1) = stime;
  time_data(2) = (double) maxiters;
  if (info3) error("Error in sp: %s", infostr.str());

  retval(0) = xx;
  #ifdef DEBUG
  cout << "retval(0) = xx=" << xx << endl;
  #endif

  retval(1) = ZZ;
  #ifdef DEBUG
  cout << "retval(1) = ZZ = " << ZZ << endl;
  #endif

  retval(2) = upper_limit;
  #ifdef DEBUG
  cout << "retval(2) = " << upper_limit << endl;
  #endif

  retval(3) = infostr.str();
  #ifdef DEBUG
  cout << "retval(3) = " << infostr.str() << endl;
  #endif

  retval(4) = time_data;
  #ifdef DEBUG
  cout << "retval(4) = " << time_data << endl;
  #endif

  #ifdef DEBUG
  cout << "sp(cc): returning" << endl;
  #endif
  return retval;
}

/*
;;; Local Variables: ***
;;; mode: C++ ***
;;; End: ***
*/

Generated by  Doxygen 1.6.0   Back to index