makereal.c

经济学专业代码

/* ************************************************************
%                                                y = makereal(x,K,cpx)
% MAKEREAL  Converts matrix in MATLAB-complex format to internal
%  SeDuMi format.
%
% SEE ALSO sedumi.
% ******************** INTERNAL FUNCTION OF SEDUMI ********************

% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% 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
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program 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, write to the Free Software
% Foundation, Inc.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA

   ************************************************************ */

#include <math.h>
#include <string.h>
#include "mex.h"
#include "blksdp.h"

#define Y_OUT plhs[0]
#define NPAROUT 1

#define X_IN prhs[0]
#define K_IN prhs[1]
#define CPX_IN prhs[2]
#define NPARIN 3

/* ************************************************************
   PROCEDURE intscalaradd - Let y = d + x.
   ************************************************************ */
void intscalaradd(int *y, const int d, const int *x, const int n)
{
  int i;
  for(i = 0; i < n; i++)
    y[i] = x[i] + d;
}

/* ************************************************************
   PROCEDURE writenz
   ************************************************************ */
void writenz(int *yir, double *ypr, int *pjnz, const int *xir,
             const double *xpr, const int jc0, const int jc1, const int d)
{
  int i,jnz;
  double yi;
  jnz = *pjnz;
/* ------------------------------------------------------------
   Write nonzeros from x into y that are REALLY nonzero; also
   translate subscript by d.
   ------------------------------------------------------------ */
  for(i = jc0; i < jc1; i++)
    if( (yi = xpr[i]) != 0.0){          /* Really nonzero? */
      ypr[jnz] = yi;
      yir[jnz] = xir[i] + d;
      jnz++;
    }
  *pjnz = jnz;
}

/* ************************************************************
   PROCEDURE fmakereal
   INPUT
     xir,xpi - sparse imaginary vector with xnnz "nonzeros" in LP/Lorentz.
        If xpi == NULL, then all imaginary data is zero.
     xnnz - length of xir,xpi, before PSD part.
     cpxf - length nf integer array, listing free imaginary vars.
     nf - length of cpxf
     iwsize Length of iwork, should be 2+nf + floor(log_2(1+nf)).
   OUTPUT
     yir, ypr - sparse real output vector, ynnz nonzeros.
   WORK ARRAYS
     cfound - length MAXN := MAX(nf,nx,2*ns) character working array.
     iwork - lengt iwsize working array. Needs
      iwsize >= 2+nf + floor(log_2(1+nf)).
   RETURNS ynnz (<=2*xnnz), number of free imag nonzeros in y.
   ************************************************************ */
int fmakereal(int *yir, double *ypr, const int *xir, const double *xpi,
	     const int xnnz, const int *cpxf, const int nf,
	     char *cfound, int *iwork, int iwsize)
{
  int i,jnz;
  int *ipos;
  double yj;

  if(xpi == (double *) NULL)
    return 0;                  /* No imaginary nonzeros */
/* ------------------------------------------------------------
   Partition WORKING ARRAY:
   ------------------------------------------------------------ */
  ipos = iwork;
  iwork += 2 + nf;
  iwsize -= 2 + nf;
/* ------------------------------------------------------------
   Locate cpx.f in xir(1:xnnz)
   ------------------------------------------------------------ */
  if(intmbsearch(ipos, cfound, xir, xnnz, cpxf, nf, iwork, iwsize) != 0)
    mexErrMsgTxt("Out of working space");
/* ------------------------------------------------------------
   Write y = sparse(imag(x(cpx.f))), the free imaginary components
   ------------------------------------------------------------ */
  jnz = 0;
  for(i = 0; i < nf; i++)
    if(cfound[i] != 0){
      if( (yj = xpi[ipos[i+1]]) != 0.0){
	ypr[jnz] = yj;
	yir[jnz] = i;
	jnz++;
      }
    }
/* ------------------------------------------------------------
   RETURN number of (free imag) nonzeros in y
   ------------------------------------------------------------ */
  return jnz;
}


/* ************************************************************
   PROCEDURE xmakereal
   INPUT
     xir,xpr,xpi - sparse vector with xnnz nonzeros in LP/Lor.
     xnnz - length of xir,xpr,xpi; counting only LP/Lor nonzeros.
     idelta - Subscript adjustment for x(1), i.e. nf.
     cpxx - length nx integer array, listing Lorentz constrained
       imaginary vars.
     nx - length of cpxx.
     iwsize Length of iwork, should be 2+nx + floor(log_2(1+nx)).
   OUTPUT
     yir, ypr - sparse real output vector, ynnz nonzeros.
   WORK ARRAYS
     cfound - length nx character working array.
     iwork - lengt iwsize working array. Needs
      iwsize >= 2+nx + floor(log_2(1+nx)).
   RETURNS ynnz (<=2*xnnz), number of nonzeros in y.
   ************************************************************ */
int xmakereal(int *yir, double *ypr,
	      const int *xir, const double *xpr, const double *xpi,
	      const int xnnz, const int idelta, const int *cpxx, const int nx,
	      char *cfound, int *iwork, int iwsize)
{
  int i,j,k,inz,jnz;
  int *ipos;
  double yj;
/* ------------------------------------------------------------
   Partition WORKING ARRAY:
   ------------------------------------------------------------ */
  ipos = iwork;
  iwork += 2 + nx;
  iwsize -= 2 + nx;
/* ------------------------------------------------------------
   Locate cpx.x in xir(1:xnnz)
   ------------------------------------------------------------ */
  if(intmbsearch(ipos, cfound, xir, xnnz, cpxx, nx, iwork, iwsize) != 0)
    mexErrMsgTxt("Out of working space");
/* ------------------------------------------------------------
   Write y(nf:nf+dimflqr+nx) = merge( real(x), imag(x(cpx.xcomplex)) )
   ------------------------------------------------------------ */
  memcpy(ypr, xpr, ipos[1] * sizeof(double));
  if(idelta != 0)
    intscalaradd(yir, idelta, xir, ipos[1]);
  else
    memcpy(yir, xir, ipos[1] * sizeof(int));
  jnz = ipos[1];
  for(i = 0; i < nx; ){
    if(cfound[i++]){
      inz = ipos[i];
      j = xir[inz] + idelta + i;         /* new index of imag part */
      if((yj = xpr[inz]) != 0.0){
	ypr[jnz] = yj;
	yir[jnz] = j-1;              /* index of real part */
	jnz++;
      }
      if(xpi != (double *) NULL)
        if((yj = xpi[inz]) != 0.0){
          ypr[jnz] = yj;
          yir[jnz] = j;                /* imag part */
          jnz++;
        }
      inz++;   /* point to first nonzero in x beyond cpx.x(i) */
    }
    else
      inz = ipos[i];    /*  point to first nonzero in x beyond cpx.x(i) */
/* ------------------------------------------------------------
   Let y(nf+i+(cpx.x(i)+1:cpx.x(i+1)-1)) = real(x(cpx.x(i)+1:cpx.x(i+1)-1))
   If i==nx, then this is simply the remainder.
   ------------------------------------------------------------ */
    k = ipos[i+1] - inz;                  /* number of nonzeros to copy */
    memcpy(ypr + jnz, xpr + inz, k * sizeof(double));
    intscalaradd(yir+jnz, idelta + i, xir+inz, k);     /* adjust subscript */
    jnz += k;
  }
/* ------------------------------------------------------------
   RETURN number of nonzeros written into y
   ------------------------------------------------------------ */
  return jnz;
}

/* ************************************************************
   PROCEDURE smakereal
   INPUT
     xir,xpr,xpi - sparse vector with xnnz nonzeros, in PSD ONLY!.
     xnnz - length of xir,xpr,xpi (PSD only).
     idelta - subscript adjustment of 1st PSD entry
     cpxsi - length twons increasing integer array. The old subscripts
       of the kth Hermitian block in x are cpxsi[2*k]:cpxsi[2*k+1]-1.
     twons - twons/2 is number of Hermitian PSD blocks.
     lenfull - length of full(x)-vector, 1 beyond last possible subscript.
     iwsize Length of iwork, should be 2*(1+ns) + floor(log_2(1+2*ns)).
   OUTPUT
     yir, ypr - sparse real output vector, ynnz nonzeros.
   WORK ARRAYS
     cfound - length MAXN := MAX(nf,nx,2*ns) character working array.
     iwork - lengt iwsize working array. Needs
      iwsize >= 2+2*ns + floor(log_2(1+2*ns)).
   RETURNS ynnz (<=2*xnnz), number of nonzeros in y.
   ************************************************************ */
int smakereal(int *yir, double *ypr,
	      const int *xir, const double *xpr, const double *xpi,
	      const int xnnz, const int idelta,
	      const int *cpxsi, const int twons, const int lenfull,
	      char *cfound, int *iwork, int iwsize)
{
  int i,j,k,inz,jnz;
  int *ipos;
/* ------------------------------------------------------------
   Partition WORKING ARRAY:
   ------------------------------------------------------------ */
  ipos = iwork;         /* length 2*(1+ns) array */
  iwork += 2 + twons;
  iwsize -= 2 + twons;
/* ------------------------------------------------------------
   Locate cpxsi in x(jcs:end); these mark the start+ends of Hermitian blocks,
   and hence the complement are real blocks.
   ------------------------------------------------------------ */
  if(intmbsearch(ipos, cfound, xir, xnnz, cpxsi, twons, iwork, iwsize) != 0)
    mexErrMsgTxt("Out of working space");