findblks.c

/* ************************************************************
%                         Ablk = findblks(At,Ablkjc,blk0,blk1,blkstart)
% FINDBLKS  Find nonzero blocks
%   in A, with subscripts per column bounded bij Ablkjc([blk0,blk1]),
%   block partitioned by blkstart.
%   If blk0 < 1 (blk1 > size(Ablkjc,2)) then start (stop) searching at column
%   start (end) of A.
%
%  SEE ALSO partitA.
% ******************** INTERNAL FUNCTION OF SEDUMI ********************

function Ablk = findblks(At,Ablkjc,blk0,blk1,blkstart) --  Find nonzero blocks

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

#define ABLK_OUT plhs[0]
#define NPAROUT 1

#define AT_IN prhs[0]
#define ABLKJC_IN prhs[1]
#define BLK0_IN prhs[2]
#define BLK1_IN prhs[3]
#define BLKSTART_IN prhs[4]
#define NPARIN 5

/* ************************************************************
   PROCEDURE findblks
   INPUT
     Ajc1, Ajc2, Air, m  - sparse N x m matrix, we only consider nonzeros
       in Air[Ajc1[k]:Ajc2[k]-1], k=0:m-1.
     blkstart, nblk - length nblk integer array of subscripts.
     blkstartm1 - length nblk array, blkstartm1[k]=blkstart[k]-1
     iwsize - length of iwork, iwsize = nblk+2+floor(log(1+nblk)/log(2)).
   OUTPUT
     Ablkjc, Ablkir - sparse nblk x m matrix, less than sum(Ajc2-Ajc1)
       nonzeros.
   WORK
     cfound - length nblk char work array
     iwork  - length iwsize = nblk+2+floor(log(1+nblk)/log(2)) work array.
   ************************************************************ */
void findblks(mwIndex *Ablkir, mwIndex *Ablkjc, const mwIndex *Ajc1,const mwIndex *Ajc2,
              const mwIndex *Air, const mwIndex *blkstart, const mwIndex *blkstartm1,
              const mwIndex m,const mwIndex nblk,
              mwIndex iwsize, bool *cfound, mwIndex *iwork)
{
  mwIndex i,j,inz,ajnnz;
  mwIndex *ipos;
/* ------------------------------------------------------------
   Partition working array into ipos(nblk+2), iwork.
   ------------------------------------------------------------ */
  ipos = iwork;
  iwork += nblk+2;
  iwsize -= nblk+2;
  inz = 0;
  for(j = 0; j < m; j++){
    Ablkjc[j] = inz;
/* ------------------------------------------------------------
   If A(:,j) has more nonzeros than blkstart, we search blkstart
   ------------------------------------------------------------ */
    if((ajnnz = Ajc2[j]-Ajc1[j]) > nblk){
      intmbsearch(ipos, cfound, Air+Ajc1[j], ajnnz,
                  blkstart, nblk, iwork, iwsize);
      for(i = 0; i < nblk; i++)
        if(ipos[i+2] > ipos[i+1])
          Ablkir[inz++] = i;
    }
    else{
/* ------------------------------------------------------------
   If A(:,j) has less nonzeros than blkstart, we search those nonzeros
   within blkstartm1. The position of the nonzero is then the block number+1.
   ------------------------------------------------------------ */
      intmbsearch(ipos, cfound, blkstartm1,nblk, Air+Ajc1[j],ajnnz,
                  iwork, iwsize);
      for(i = 0; i < ajnnz; i++)
        if(ipos[i+1] > ipos[i])         /* New block number ? */
          Ablkir[inz++] = ipos[i+1] - 1;   /* ipos is block number + 1 */
    }
  }
/* ------------------------------------------------------------
   Close last column of Ablk
   ------------------------------------------------------------ */
  Ablkjc[m] = inz;
}

/* ************************************************************
   PROCEDURE mexFunction - Entry for Matlab
   ************************************************************ */
void mexFunction(int nlhs, mxArray *plhs[],
                int nrhs, const mxArray *prhs[])
{
  jcir At, Ablk;
  mwIndex i,j, nblk,m, blknnz, njc, iwsize, blk0,blk1;
  mwIndex *iwork, *Ajc, *blkstart;
  const double *blkstartPr, *AjcPr;
  bool *cwork;
  bool isblk0negative;
/* ------------------------------------------------------------
   Check for proper number of arguments
   ------------------------------------------------------------ */
  mxAssert(nrhs >= NPARIN, "findblks requires more input arguments.");
  mxAssert(nlhs <= NPAROUT, "findblks produces less output arguments.");
/* --------------------------------------------------
   GET inputs At, blkstart, Ablkjc, blk0, blk1
   -------------------------------------------------- */
  mxAssert(mxIsSparse(AT_IN), "At must be a sparse matrix.");
  At.jc = mxGetJc(AT_IN);
  At.ir = mxGetIr(AT_IN);
  m = mxGetN(AT_IN);
  nblk = mxGetM(BLKSTART_IN) * mxGetN(BLKSTART_IN) - 1;
  blkstartPr = mxGetPr(BLKSTART_IN);
  AjcPr = mxGetPr(ABLKJC_IN);
  mxAssert(m == mxGetM(ABLKJC_IN), "Ablkjc size mismatch.");
  njc = mxGetN(ABLKJC_IN);
  blk0 = (mwIndex) mxGetScalar(BLK0_IN);           /* double to mwIndex */
  isblk0negative=0;
  mxAssert(blk0>0,"");
  if(blk0>0)
    --blk0;  /* Fortran to C */
  else
    isblk0negative=1;
  if(mxGetM(BLK1_IN) * mxGetN(BLK1_IN) != 1)
    blk1 = njc;                           /*default to end */
  else{
    blk1 = (mwIndex) mxGetScalar(BLK1_IN);   /* double to mwIndex (thus inf not allowed) */
    mxAssert(blk1>0,"");
    --blk1;                                /* Fortran to C */
  }
/* ------------------------------------------------------------
   Allocate working array iwork(nblk+2+log_2(1+nblk)),
   blkstart(2*nblk), Ajc(2*m)
   char cwork(nblk)
   ------------------------------------------------------------ */
  iwsize = nblk + 2 + (mwIndex) floor(log(1.0+nblk)/log(2.0));
  iwork = (mwIndex *) mxCalloc(iwsize, sizeof(mwIndex));
  blkstart = (mwIndex *) mxCalloc(MAX(2*nblk,1), sizeof(mwIndex));
  Ajc = (mwIndex *) mxCalloc(MAX(2*m,1), sizeof(mwIndex));
  cwork = (bool *) mxCalloc(MAX(nblk,1), sizeof(bool));
/* ------------------------------------------------------------
   Translate blkstart from Fortran-double to C-mwIndex
   ------------------------------------------------------------ */
  for(i = 0; i < nblk; i++){                         /* to integers */
    j = (mwIndex) blkstartPr[i];
    mxAssert(j>0,"");
    blkstart[i] = --j;
    mxAssert(j>0,"");    
    blkstart[nblk+i] = --j;          /* blkstart minus 1 */
  }
/* ------------------------------------------------------------
   Convert Ajc from double to mwIndex:
   ------------------------------------------------------------ */
  mxAssert(blk0 < njc, "Ablkjc size mismatches blk0.");
  if(isblk0negative)
    memcpy(Ajc,At.jc,m*sizeof(mwIndex));          /* default: start of column */
  else
    for(i = 0; i < m; i++){                         /* to integers */
      Ajc[i] = (mwIndex) AjcPr[m*blk0 + i];
    }
  mxAssert(blk1 >= 0, "blk1 must be positive.");
  if(blk1 >= njc)
    memcpy(Ajc+m,At.jc+1,m*sizeof(mwIndex));      /* default: end of column */
  else
    for(i = 0; i < m; i++){                         /* to integers */
      Ajc[m+i] = (mwIndex) AjcPr[blk1*m + i];
    }
/* ------------------------------------------------------------
   Ablk = sparse(nblk,m,blknnz);
   ------------------------------------------------------------ */
  blknnz = 0;
  for(i = 0; i < m; i++)
    blknnz += Ajc[m+i]-Ajc[i];              /* upper bound on nnz blocks */
  blknnz = MAX(blknnz,1);
  ABLK_OUT = mxCreateSparse(nblk,m, blknnz,mxREAL);
  Ablk.jc = mxGetJc(ABLK_OUT);
  Ablk.ir = mxGetIr(ABLK_OUT);
/* ------------------------------------------------------------
   The real job:
   ------------------------------------------------------------ */
  findblks(Ablk.ir,Ablk.jc, Ajc,Ajc+m,At.ir, blkstart,blkstart+nblk,
           m,nblk, iwsize,cwork,iwork);
/* ------------------------------------------------------------
   REALLOC (shrink) Ablk to Ablk.jc[m] nonzeros.
   ------------------------------------------------------------ */
  mxAssert(Ablk.jc[m] <= blknnz,"");
  blknnz = MAX(Ablk.jc[m],1);
  if((Ablk.ir = (mwIndex *) mxRealloc(Ablk.ir, blknnz * sizeof(mwIndex))) == NULL)
    mexErrMsgTxt("Memory allocation error");
  if((Ablk.pr = (double *) mxRealloc(mxGetPr(ABLK_OUT), blknnz*sizeof(double)))
     == NULL)
    mexErrMsgTxt("Memory allocation error");
  mxSetPr(ABLK_OUT,Ablk.pr);
  mxSetIr(ABLK_OUT,Ablk.ir);
  mxSetNzmax(ABLK_OUT,blknnz);
  for(i = 0; i < blknnz; i++)
    Ablk.pr[i] = 1.0;
/* ------------------------------------------------------------
   Release working arrays
   ------------------------------------------------------------ */
  mxFree(cwork);
  mxFree(iwork);
  mxFree(Ajc);
  mxFree(blkstart);
}