spscale.c

经济学专业代码

/*
% 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 <string.h>
#include "mex.h"
#include "blksdp.h"

/* ************************************************************
   PROCEDURE realtransp - Let Y = X', where X is m x n full.
   You may also view this as going from column- towards row-stacking.
   ************************************************************ */
void realtransp(double *y, const double *x, const int m, const int n)
{
  int i, j, ji;
/* ------------------------------------------------------------
   Y(n x m), X(m x n).
   For each column y:= yj, j=0:m-1
   Let yj[i] = xji, i = 0:n-1.
   ------------------------------------------------------------ */
  for(j = 0; j < m; j++, y += n)
    for(i = 0, ji = j; i < n; i++, ji += m)
      y[i] = x[ji];            /* yij = xji */
}

/* =========================  P S D :  ========================= 
   D(d^2)x = vec(D * X * D), only at zir positions, with X sparse.
   ============================================================ */

/* ************************************************************
   PROCEDURE realDmulX - Computes Y = D*X, where X is sparse.
   INPUT
     d  - full n x n matrix D.
     xpr, xir, xjc1 - sparse vector, xjc1 nonzeros.
     first - Subscript of X(1,1), i.e. x(first:first+n^2-1) = vec(X).
     n - order of D, X.
   UPDATED
     *pxjc0 - On input, points to first nonzero with index in range
       first:first+n^2-1. On output, points to first nonzero AFTER this range.
   OUTPUT
     y  - Full n x knz output matrix. Y = D*X(:,ycols).
     ycols - length knz <= n integer array, listing all columns where X
       and hence Y is nonzero.
   RETURNS knz = length(ycols), number of nonzero columns.
   ************************************************************ */
int realdmulx(double *y, int *ycols, const double *d,
              const double *xpr, const int *xir, int *pxjc0, const int xjc1,
              const int first, const int n)
{
  int knz, jfirst, jlast, inz, i, j;
  knz = 0;        /* length(ycols) */
  jlast = 0;      /* index right after last activated column */
  y -= n;         /* point to dummy column, which will be skipped */
/* ------------------------------------------------------------
   For every nonzero xij, Let Y += D*(xij * ei*ej'), i.e.
   yj += xij * di. Let ycols list the columns where yj is nonzero.
   Store only those columns in y. (ycols is increasing)
   ------------------------------------------------------------ */
  for(inz = *pxjc0; inz < xjc1; inz++)
    if((i = xir[inz]) >= jlast){      /* we need to create new y-column */
      j = (i-first) / n;              /* col j = floor( (i-first)/n ) */
      if(j >= n)
        break;                        /* finished all columns of X */
      ycols[knz++] = j;               /* open new column */
      y += n;
      jfirst = first + n * j;
      jlast = jfirst + n;
      i = i - jfirst;                  /* i = rem( (i-first)/n ) */
      scalarmul(y, xpr[inz], d + i * n, n);               /* yj = xij * di */
    }
    else
      addscalarmul(y, xpr[inz], d + (i-jfirst) * n, n);   /* yj += xij * di */
/* ------------------------------------------------------------
   RETURN xjc0NEW = inz, pointer to next entry in sparse x.
   knz = length(ycols), the number of nonzero columns in y.
   ------------------------------------------------------------ */
  *pxjc0 = inz;
  return knz;
}

/* ************************************************************
   PROCEDURE cpxDmulX - Computes Y = D*X, where X is sparse.
      Y,D, and X are complex matrices, stored as [vec(RE); vec(IM)].
   INPUT
     d  - full 2*(n x n) matrix D.
     xpr, xir, xjc1 - sparse vector, xjc1 nonzeros (some RE, some IM).
     first - Subscript of RE X(1,1), i.e. x(first:first+n^2-1) = vec(RE(X)),
        x(first+n^2:first+2n^2-1) = vec(IM(X)).
     n - order of D, X.
   UPDATED
     *pxjc0 - On input, points to first nonzero with index in range
       first:first+2*n^2-1. On output, first nonzero AFTER this range.
   OUTPUT
     y  - Full 2*(n x n) output matrix. y(1:n*knz) = RE D*X(:,ycols),
       y(n^2+(1:n*knz)) = IM D*X(:,ycols).
     ycols - length knz <= n integer array, listing all columns where X
       and hence Y is nonzero.
   RETURNS knz = length(ycols), number of nonzero columns.
   ************************************************************ */
int cpxdmulx(double *y, int *ycols, const double *d,
             const double *xpr, const int *xir, int *pxjc0, const int xjc1,
             const int first, const int n)
{
  int jfirst, jlast, inz, jnz, knz, i, j, icol, imgfirst, nsqr, ncols;
  double *ypr, *ypi;
  const double *dpi;
  char found;
  knz = 0;        /* length(ycols) */
  jlast = 0;      /* index right after last activated column */
  nsqr = SQR(n);
  dpi = d + nsqr;
/* ------------------------------------------------------------
   For every REAL nonzero xij, Let Y += D*(xij * ei*ej'), i.e.
   RE(yj) += xij * di, and IM(yj) += xij *IM(di).
   Let ycols list the columns where yj is nonzero due to REAL X.
   Store only those columns in y. (ycols is increasing - for now)
   ------------------------------------------------------------ */
  ypr = y-n;         /* point to dummy column, which will be skipped */
  ypi = ypr + nsqr;
  for(inz = *pxjc0; inz < xjc1; inz++)
    if((i = xir[inz]) >= jlast){      /* we need to create new y-column */
      j = (i-first) / n;              /* col j = floor( (i-first)/n ) */
      if(j >= n)
        break;                        /* finished all real columns of X */
      ycols[knz++] = j;               /* open new column */
      ypr += n;
      ypi += n;
      jfirst = first + n * j;
      jlast = jfirst + n;
      icol = (i - jfirst) * n;               /* i = rem( (i-first)/n ) */
      scalarmul(ypr, xpr[inz], d + icol, n);               /* yj = xij * di */
      scalarmul(ypi, xpr[inz], dpi + icol, n);
    }
    else{
      icol = (i-jfirst) * n;
      addscalarmul(ypr, xpr[inz], d + icol, n); /* yj += xij * di */
      addscalarmul(ypi, xpr[inz], dpi + icol, n);
    }
/* ------------------------------------------------------------
   Finished with RE(X), yielding ncols:=knz columns in Y. Use
   jnz to browse through these existing cols while processing IMG(X).
   ------------------------------------------------------------ */
  ncols = knz;            /* #existing cols, due to RE(X) */
  jnz = 0;                /* pointer into ycols(1:ncols) */
/* ------------------------------------------------------------
   For every IMAG nonzero xij, Let Y += iD*(xij * ei*ej'), where
   iD := sqrt(-1)*D. Thus:
   RE(yj) -= xij * IM(di), and IM(yj) += xij *RE(di).
   New nonzero cols of y are listed in ycols(ncols+1:knz).
   This means that ycols consists of 2 disjoint increasing lists.
   ------------------------------------------------------------ */
  imgfirst = first + nsqr;
  for(; inz < xjc1; inz++)
    if((i = xir[inz]) >= jlast){      /* we need to create new y-column */
      j = (i-imgfirst) / n;           /* col j = floor( (i-first)/n ) */
      if(j >= n)
        break;                        /* finished all columns of X */
      jfirst = imgfirst + n * j;
      jlast = jfirst + n;
      icol = (i - jfirst) * n;               /* i = rem( (i-first)/n ) */
/* ------------------------------------------------------------
   1st time to use column j while processing IM(X). See whether
   yj already in ycols (=>found=1). Otherwise, create new yj.
   ------------------------------------------------------------ */
      for(found = 0; jnz < ncols; jnz++)
        if(ycols[jnz] >= j){
          found = (ycols[jnz] == j);
          break;                        /* ycols(1:ncols) is increasing */
        }
      if(!found){
        ypr = y + n * knz;              /* IM(X) results in new yj */
        ypi = ypr + nsqr;
        ycols[knz++] = j;               /* open new column */
        scalarmul(ypr, -xpr[inz], dpi + icol, n);     /* yj = xij * (i*di) */
        scalarmul(ypi, xpr[inz], d + icol, n);
      }
      else{
        ypr = y + n * jnz;             /* yj already exists */
        ypi = ypr + nsqr;
        addscalarmul(ypr, -xpr[inz], dpi + icol, n); /* yj += xij * (i*di) */
        addscalarmul(ypi, xpr[inz], d + icol, n);
      }
    }
    else{
      icol = (i-jfirst) * n;
      addscalarmul(ypr, -xpr[inz], dpi + icol, n); /* yj += xij * (i*di) */
      addscalarmul(ypi, xpr[inz], d + icol, n);
    }
/* ------------------------------------------------------------
   RETURN xjc0NEW = inz, pointer to next entry in sparse x.
   knz = length(ycols), the number of nonzero columns in y.
   ------------------------------------------------------------ */
  *pxjc0 = inz;
  return knz;
}

/* ************************************************************
   PROCEDURE sprealdxd - Computes Z = D*sym(X)*D with D = D' and
        sym(X) = (X+X')/2. X sparse, and Z has pre-chosen sparsity
        structure zir (typically subset of tril-entries).
   INPUT
     zir, zjc0, zjc1 - Target sparsity structure. We compute z(i) only for
       i in zir(zjc0:zjc1-1) with i < first + n^2.
       zjc0 such that zir(zjc0) >= first.
     d      - full n x n matrix D.
     xpr, xir, xjc1 - sparse vector, xjc1 nonzeros.
     first  - Subscript of X(1,1) and Z(1,1), i.e.
       x(first:first+n^2-1) = vec(X).
     n      - order of D, X.
   UPDATED
     *pxjc0 - On input, points to first nonzero with index in range
       first:first+n^2-1. On output, points to first nonzero AFTER this range.
   OUTPUT
     Z - lenfull vector. Only z(first:first+n^2-1) is affected, and from this
       only the indices listed in zir. (Remaining entries unaffected.)
   WORK
     dxcols - length n integer array.
     fwork  - 2 * n^2 vector.
   ************************************************************ */
void sprealdxd(double *z, const int *zir, const int znnz,
               const double *d,
               const double *xpr, const int *xir, int *pxjc0, const int xjc1,
               const int first, const int n, double *fwork, int *dxcols)
{
  int inz, i, icol, j, jfirst, jlast, m;
  double *xd;
  const double *dj, *xdj;
/* ------------------------------------------------------------
   Partition 2*n^2 WORKING array fwork into [fwork(n^2), xd(n^2)].
   ------------------------------------------------------------ */