Page MenuHomec4science
Files F121914377

partitA.c
No OneTemporary
Actions

File Metadata

Created
Mon, Jul 14, 19:32

partitA.c
View Options

/* ************************************************************
% Ablkjc = partitA(At,blkstart)
% PARTITA Partition columns of A according to the subscripts listed in blkstart.
%
% SEE ALSO sedumi
% ******************** INTERNAL FUNCTION OF SEDUMI ********************
function Ablkjc = partitA(At,blkstart) -- Partition columns of A
according to the subscripts listed in blkstart
% 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> /* floor and log */
#include "mex.h"
#include "blksdp.h"
#define ABLKJC_OUT plhs[0]
#define NPAROUT 1
#define AT_IN prhs[0]
#define BLKSTART_IN prhs[1]
#define NPARIN 2
void intadd(mwIndex *x, const mwIndex y, const mwIndex n)
{
mwIndex i;
for(i = 0; i < n; i++)
x[i] += y;
}
/* ************************************************************
PROCEDURE partitA
INPUT
Ajc, Air, m - sparse N x m matrix
blkstart, nblk - length nblk integer array of subscripts.
iwsize - length of iwork, iwsize = floor(log(1+nblk)/log(2)).
OUTPUT
Ablkjc - length (nblk+1)*m array. Rows 1+(1:nblk) list 1st nonzero
with subscript at or beyond blkstart.
WORK
cfound - length nblk char work array
iwork - length iwsize = floor(log(1+nblk)/log(2)) work array.
************************************************************ */
void partitA(mwIndex *Ablkjc, const mwIndex *Ajc,const mwIndex *Air,
const mwIndex *blkstart, const mwIndex m,const mwIndex nblk,
const mwIndex iwsize, bool *cfound, mwIndex *iwork)
{
mwIndex j, L;
L = nblk+2;
for(j = 0; j < m; j++)
intmbsearch(Ablkjc + j*L, cfound, Air+Ajc[j], Ajc[j+1]-Ajc[j],
blkstart, nblk, iwork, iwsize);
for(j = 0; j < m; j++)
intadd(Ablkjc + j*L, Ajc[j],L);
}
/* ************************************************************
PROCEDURE mexFunction - Entry for Matlab
************************************************************ */
void mexFunction(int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[])
{
jcir At;
mwIndex i,j, nblk,m, L, iwsize;
mwIndex *iwork, *Ablkjc, *blkstart;
const mwIndex *rowj;
double *AblkjcPr;
const double *blkstartPr;
bool *cwork;
/* ------------------------------------------------------------
Check for proper number of arguments
------------------------------------------------------------ */
mxAssert(nrhs >= NPARIN, "partitA requires more input arguments.");
mxAssert(nlhs <= NPAROUT, "partitA produces less output arguments.");
/* --------------------------------------------------
GET inputs At, blkstart
-------------------------------------------------- */
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);
blkstartPr = mxGetPr(BLKSTART_IN);
/* ------------------------------------------------------------
Allocate working array Ablkjc((nblk+2) * m), iwork(log_2(1+nblk)),
blkstart(nblk)
------------------------------------------------------------ */
iwsize = (mwIndex) floor(log(1.0+nblk)/log(2.0));
iwork = (mwIndex *) mxCalloc(MAX(iwsize,1), sizeof(mwIndex));
Ablkjc = (mwIndex *) mxCalloc(MAX((nblk+2)*m,1), sizeof(mwIndex));
blkstart = (mwIndex *) mxCalloc(MAX(nblk,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;
}
/* ------------------------------------------------------------
The real job:
------------------------------------------------------------ */
partitA(Ablkjc, At.jc,At.ir, blkstart, m,nblk, iwsize,cwork,iwork);
/* ------------------------------------------------------------
Create output Ablkjc m x nblk.
------------------------------------------------------------ */
ABLKJC_OUT = mxCreateDoubleMatrix(m, nblk, mxREAL);
AblkjcPr = mxGetPr(ABLKJC_OUT);
rowj = Ablkjc;
L = nblk+2;
for(j = 0; j < nblk; j++){
++rowj;
for(i = 0; i < m; i++)
AblkjcPr[i] = (double) rowj[i*L]; /* convert mwIndex to double */
AblkjcPr += m;
}
/* ------------------------------------------------------------
Release working arrays
------------------------------------------------------------ */
mxFree(cwork);
mxFree(iwork);
mxFree(Ablkjc);
mxFree(blkstart);
}

Event Timeline

c4science · Help