symbfwblk.c

matlab中uwb波形优化算法经常会使用的工具包:SeDuMi_1_1R3.

/*
 x = symbfwblk(L, b)

% 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"

#define X_OUT  plhs[0]
#define NPAROUT 1

#define L_IN      prhs[0]
#define B_IN      prhs[1]
#define NPARIN 2

/* ************************************************************
   PROCEDURE snodeCompress  - Compressed subscripts based on
                    supernodal partition.
   INPUT
     ljc, lir - uncompressed nz-structure of L
     xsuper - supernodal partition (length nsuper+1).
     nsuper - number of supernodes
   OUTPUT
     xlindx - length nsuper+1, the columns in lindx.
     lindx  - compressed subscript array: L(:,xsuper). Should be allocated
       to ljc[m], so that there are certainly enough entries.
     snode  - length m = xsuper[nsuper+1]. Maps subnode to its supernode.
   ************************************************************ */
void snodeCompress(int *xlindx,int *lindx,int *snode,
                   const int *ljc,const int *lir,const int *xsuper,
                   const int nsuper)
{
  int j, jsup, ix, collen, jcol;
/* ------------------------------------------------------------
   SNODE: map each column to the supernode containing it
   ------------------------------------------------------------ */
  j = xsuper[0];
  for(jsup = 0; jsup < nsuper; jsup++){
    while(j < xsuper[jsup + 1])
      snode[j++] = jsup;
  }
/* ------------------------------------------------------------
   COMPRESS SUBSCRIPTS:
    Let (xlindx,lindx) = ljc(xsuper(:)), i.e store only once
    for each snode, instead of once per column.
   ------------------------------------------------------------ */
  for(ix = 0, jsup = 0; jsup < nsuper; jsup++){
    xlindx[jsup] = ix;
    jcol = xsuper[jsup];
    collen = ljc[jcol+1] - ljc[jcol];
    memcpy(lindx + ix, lir + ljc[jcol], collen * sizeof(int));
    ix += collen;
  }
  xlindx[nsuper] = ix;
}

/* ************************************************************
   PROCEDURE getnzfwlj - find nonzero supernodes in L\e_{xsuper[jsup]}.
   INPUT
     jsup  - starting supernode to process from.
     snode, xsuper - supernodal partition.
         xsuper(nsuper+1): xsuper(j):xsuper(j+1)-1 is jth supernode
         snode(m): j=snode(i) means xsuper(j) <= i < xsuper(j+1).
     xlindx,lindx  - compressed subscript array.
         xlindx(nsuper+1): lindx(xlindx(j):xlindx(j+1)-1) are the subscripts
         for supernode j.
   UPDATED
     processed - Sets processed[jsup] = 1 if jsup is in L\e_{xsuper[jsup]}.
     snodefrom - Lists first relevant subnode of each supernode 
       where processed[jsup]=1.
   REMARK - caller has to set processed[jsup]=1; getnzfwlj only does
     this for the child supernodes.
   ************************************************************ */
void getnzfwlj(int *snodefrom, char *processed, int jsup,
              const int *snode, const int *xsuper,
              const int *xlindx, const int *lindx)
{
  int i,j;

  j = xsuper[jsup+1] - xsuper[jsup];
  while(xlindx[jsup] + j < xlindx[jsup + 1]){
    i = lindx[xlindx[jsup] + j];
    jsup = snode[i];         /* next affected snode */
/* ------------------------------------------------------------
   If jsup has already been processed, then we can stop here, after
   making sure that i >= snodefrom[jsup].
   ------------------------------------------------------------ */
    if(processed[jsup]){
      if(i < snodefrom[jsup])
        snodefrom[jsup] = i;
      break;             /* STOP */
    }
/* ------------------------------------------------------------
   Otherwise, we process and link through to next affected supernode
   ------------------------------------------------------------ */
    processed[jsup] = 1;
    snodefrom[jsup] = i;
    j = xsuper[jsup+1] - xsuper[jsup];
  }
}

/* ************************************************************
   PROCEDURE getnzsuper - Compute sparsity structure of L\b(perm), by
       determining nonzero-supernodes, and starting subnodes within
       them (each supernode is a dense diag block in L).
   INPUT
     bir, bnnz     - bir(bnnz) lists the row-indices of vector b.
     invperm       - int(m) Is s.t. perm[invperm[i]] = i.
     snode, xsuper - supernodal partition.
         xsuper(nsuper+1): xsuper(j):xsuper(j+1)-1 is jth supernode
         snode(m): j=snode(i) means xsuper(j) <= i < xsuper(j+1).
     xlindx,lindx  - compressed subscript array.
         xlindx(nsuper+1): lindx(xlindx(j):xlindx(j+1)-1) are the subscripts
         for supernode j.
   UPDATED
     processed - char(nsuper) array. On input all-0, on output
       processed[jsup] = 1 iff jsup is in nz structure of L\b.
   OUTPUT
     snodefrom - Length nsuper array. Lists first relevant subnode of
       each supernode where processed[jsup]=1.
   ************************************************************ */
void getnzsuper(int *snodefrom, char *processed,
                const int *bir, const int bnnz,
                const int *invperm, const int *snode, const int *xsuper,
                const int *xlindx, const int *lindx)
{
  int inz,i,jsup;
/* ------------------------------------------------------------
   We'll process each supernode jsup = snode[ bir[ inz ] ], to
   find all supernodes in x, L*x = b, and the first relevant
   subnode of each supernode, snodefrom[jsup].
   ------------------------------------------------------------ */
  if(bnnz <= 0)
    return;
  for(inz = 0; inz < bnnz; inz++){
    i = invperm[bir[inz]];            /* We're interested in b(perm) */
    jsup = snode[i];
/* ------------------------------------------------------------
   If jsup has not yet been processed, then find supernodes involved
   in solving L*x = e_i.
   ------------------------------------------------------------ */
    if(!processed[jsup]){
      snodefrom[jsup] = i;
      getnzfwlj(snodefrom,processed,jsup, snode,xsuper,xlindx,lindx);
      processed[jsup] = 1;
    }
/* ------------------------------------------------------------
   Otherwise, we only need to make sure that i >= snodefrom[jsup].
   ------------------------------------------------------------ */
    else if(i < snodefrom[jsup])
      snodefrom[jsup] = i;
  }
}    

/* ************************************************************
   PROCEDURE symbfwmat - Computes nz-structur of x = L\b(perm,:).
   INPUT
     bjc, bir - nz-structure of m x n RHS-matrix b.
     invperm       - int(m) Is s.t. perm[invperm[i]] = i.