vectril.c

经济学专业代码

/*
   y = vectril(x,K)
   For the PSD submatrices, we let Yk = tril(Xk,0) + triu(Xk,1)'
   Complex numbers are stored in SeDuMi's (real Xk; imag Xk) format.

% 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 NPARIN 2

/* ************************************************************
   PROCEDURE sptriujcT - Compute starting positions for storing
     nonzeros of rows triu(X(i,:),1). Note: X is nxn block stored
     as subvector in a long vector.
   INPUT
     xir   - integer subscript array of length xjc1.
     xjc0,xjc1 - xjc0 points to n x n sparse matrix X in xir. 
     first - subscript of X(0,0); vec(X) is stored as subvector in xir.
     n - order of X.
   OUTPUT
     triujc - length n-1 integer array; triujc[i] is starting index
       for storing triu(X(i+1,:),1), i=0:n-2.
       It starts at triujc[0] = nnz(tril(X)).
   RETURNS nnz(X). Note that
   nnz(X) = triujc[n-2] + nnz(triu(X(n-1,:),1)) <= xjc1 - xjc0.
   ************************************************************ */
int sptriujcT(int *triujc, const int *xir, const int xjc0, const int xjc1,
              const int first,const int n)
{
  int i,j,inz,jnz,jfirst,jlast;
/* ------------------------------------------------------------
   Observe that triu(X(n,:))=[] and hence only use triujc[0:n-2].
   ------------------------------------------------------------ */
  mxAssert(n > 0,"");
  memset(triujc, 0, (n-1) * sizeof(int));
  jnz = 0;        /* jnz = nnz(tril(X)) */
  jlast = 0;      /* index right after last activated column */
  for(inz = xjc0; inz < xjc1; inz++){
/* ------------------------------------------------------------
   Translate vec(x)-subscript into (i,j) subscript
   ------------------------------------------------------------ */
    if((i = xir[inz]) >= jlast){      /* move to new z-column */
      j = (i-first) / n;              /* col j = floor( (i-first)/n ) */
      if(j >= n)
        break;                        /* finished all columns of X */
      jfirst = first + n * j;
      jlast = jfirst + n;
    }
    mxAssert(i >= jfirst,"");
    i -= jfirst;
/* ------------------------------------------------------------
   x_{i,j} with i < j --> count as nonzero in row i of triu(x)
   ------------------------------------------------------------ */
    if(i < j)              /* Upper triangular nonzero in row i */
      ++triujc[i];    /* 0 <= i <= j-1 <= n-2 */
/* ------------------------------------------------------------
   x_{i,j} with i >= j -->  count as nonzero in tril(x).
   ------------------------------------------------------------ */
    else{
      ++jnz;
    }
  }
/* ------------------------------------------------------------
   Now jnz = nnz(tril(X)) and triujc[i] = nnz(triu(X(i+1,:),1)).
   We will now let triujc point to 1st nonzero position for storing
   jth row of triu(X,1) in row-wise format, as follows:
   [jnz, jnz + triujc[0], .., jnz + triujc[n-2]].
   ------------------------------------------------------------ */
  for(i = 0; i < n-1; i++){
    j = triujc[i];        /* nnz(triu(X(i,:),1)) */
    triujc[i] = jnz;      /* position to store triu(X(i,:),1) */
    jnz += j;            /* point to next available position */
  }
/* ------------------------------------------------------------
   Now jnz = triujc[n-1, NEW] + triujc[n-1,OLD] = nnz(X).
   ------------------------------------------------------------ */
  mxAssert(jnz == inz - xjc0,"");
  return jnz;
}

/* ************************************************************
   PROCEDURE sptrilandtriu - Compute y = [tril(X); triu(X,1)']
   INPUT
     xir   - integer subscript array of length xjc1.
     xpr   - vector of xjc1 nonzeros.
     xjc0,xjc1 - xjc0 points to n x n sparse matrix X in xir. 
     first - subscript of X(0,0); vec(X) is stored as subvector in xir.
     n - order of X.
   UPDATED
     triujc - length n-1 integer array. INPUT: triujc[i] is starting index
       for storing triu(X(i+1,:),1), i=0:n-2. It starts at triujc[0] =
       nnz(tril(X)). OUTPUT: triujc[i] points beyond last nonzero
       written in triu(X(i,:),1)'. Thus, triujc[n-2] = nnz(y).
   OUTPUT
     yir, ypr - length nnz(y) sparse vector, y = [tril(X); triu(X,1)'].
   RETURNS tilxnnz = nnz(tril(X)) if !skew, and nnz(tril(X,-1)) if skew.
   ************************************************************ */
int sptrilandtriu(int *yir, double *ypr, int *triujc, const int *xir,
                   const double *xpr, const int xjc0, const int xjc1,
                   const int first, const int n,const bool skew)
{
  int i,j,inz,jnz,knz,jfirst,jlast;
/* ------------------------------------------------------------
   Store [tril(X); triu(X,1)'] into y,
   using the column pointers of triu(X,1)' as computed in triujc.
   ------------------------------------------------------------ */
  jnz = 0;
  jlast = 0;      /* index right after last activated column */
  for(inz = xjc0; inz < xjc1; inz++){
/* ------------------------------------------------------------
   Translate vec(x)-subscript into (i,j) subscript
   ------------------------------------------------------------ */
    if((i = xir[inz]) >= jlast){      /* move to new z-column */
      j = (i-first) / n;              /* col j = floor( (i-first)/n ) */
      if(j >= n)
        break;                        /* finished all columns of X */
      jfirst = first + n * j;
      jlast = jfirst + n;
    }
    i -= jfirst;
/* ------------------------------------------------------------
   x_{i,j} with i < j --> store nonzero in row i of triu(x)
   ------------------------------------------------------------ */
    if(i < j){              /* Upper triangular nonzero in row i */
      knz = triujc[i];
      ++triujc[i];
      yir[knz] = first + n * i + j;   /* vec-index  after transposing */
      ypr[knz] = xpr[inz];
    }
/* ------------------------------------------------------------
   x_{i,j} with i >= j --> store nonzero of tril(x) or, for skew, tril(x,-1).
   ------------------------------------------------------------ */
    else if(!skew || i > j){
      yir[jnz] = jfirst + i;   /* vec-index */
      ypr[jnz] = xpr[inz];
      ++jnz;
    }
  }
/* ------------------------------------------------------------
   Return number of tril-nonzeros written.
   !skew ==> nnz(tril(X));  skew ==> nnz(tril(X,-1)).
   ------------------------------------------------------------ */
  return jnz;
}

/* ************************************************************
   PROCEDURE spadd - Let z = x + y
   INPUT
     xir, xpr, xnnz - sparse vector
     yir, ypr, ynnz - sparse vector
     iwsize - ynnz+2 + floor(log_2(ynnz+1))
   OUTPUT
     zir, zpr - length znnz arrays containing sparse output z = x + y.
   WORK
     iwork - length iwsize working array
   RETURNS znnz
   ************************************************************ */
int spadd2(int *zir, double *zpr, const int *xir, const double *xpr,
           const int xnnz, const int *yir, const double *ypr, const int ynnz,
           int iwsize, char *cfound, int *iwork)
{
  int inz,jnz,knz, i;
  int * yinx;
/* ------------------------------------------------------------
   Partition working array [yinx(ynnz+2), iwork(log_2(ynnz+1))].
   ------------------------------------------------------------ */
  yinx = iwork;
  iwork += ynnz + 2;
  iwsize -= ynnz + 2;
  intmbsearch(yinx, cfound, xir, xnnz, yir, ynnz, iwork, iwsize);
  jnz = yinx[1];
  memcpy(zir, xir, jnz * sizeof(int));
  memcpy(zpr, xpr, jnz * sizeof(double));
  for(i = 0; i < ynnz; i++){
    inz = yinx[i+1];
    if(cfound[i])
      zpr[jnz] = xpr[inz++] + ypr[i];
    else
      zpr[jnz] = ypr[i];
    zir[jnz++] = yir[i];
    knz = yinx[i+2]-inz;
    memcpy(zir + jnz, xir + inz, knz * sizeof(int));
    memcpy(zpr + jnz, xpr + inz, knz * sizeof(double));
    jnz += knz;
  }
  return jnz;
}
      
int spadd(int *zir, double *zpr, const int *xir, const double *xpr,
          const int xnnz, const int *yir, const double *ypr, const int ynnz,
          const int iwsize, char *cfound, int *iwork)
{
  if(xnnz < ynnz)
    return spadd2(zir,zpr, yir,ypr,ynnz, xir,xpr,xnnz, iwsize, cfound, iwork);
  else
    return spadd2(zir,zpr, xir,xpr,xnnz, yir,ypr,ynnz, iwsize, cfound, iwork);
}

/* ************************************************************
   PROCEDURE spsub - Let z = x - y
   INPUT
     xir, xpr, xnnz - sparse vector
     yir, ypr, ynnz - sparse vector
     iwsize - ynnz+2 + floor(log_2(ynnz+1))
   OUTPUT
     zir, zpr - length znnz arrays containing sparse output z = x - y.
   WORK
     iwork - length iwsize working array
   RETURNS znnz
   ************************************************************ */
int spsub(int *zir, double *zpr, const int *xir, const double *xpr,
          const int xnnz, const int *yir, const double *ypr, const int ynnz,
          int iwsize, char *cfound, int *iwork)
{
  int inz,jnz,knz, i;
  int *yinx;
/* ------------------------------------------------------------
   Partition working array [yinx(ynnz+2), iwork(log_2(ynnz+1))].
   ------------------------------------------------------------ */
  yinx = iwork;
  iwork += ynnz + 2;
  iwsize -= ynnz + 2;
  intmbsearch(yinx, cfound, xir, xnnz, yir, ynnz, iwork, iwsize);
  jnz = yinx[1];
  memcpy(zir, xir, jnz * sizeof(int));
  memcpy(zpr, xpr, jnz * sizeof(double));
  for(i = 0; i < ynnz; i++){
    inz = yinx[i+1];
    if(cfound[i])