Page MenuHomec4science

adendotd.c
No OneTemporary

File Metadata

Created
Sat, Jul 12, 21:33

adendotd.c

/*
% Ad = Adendotd(dense, d, sparAd, Ablk, blkstart)
% ADENDOTD Computes d[k]'*Aj[k] for Lorentz blocks that are to be factored
% by dpr1fact.
%
% SEE ALSO sedumi
% ********** INTERNAL FUNCTION OF SEDUMI **********
function Ad = Adendotd(dense, d, sparAd, Ablk, 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 "mex.h"
#include "blksdp.h"
#define AD_OUT plhs[0]
#define NPAROUT 1
#define DENSE_IN prhs[0]
#define D_IN prhs[1]
#define ADOTD_IN prhs[2]
#define ABLK_IN prhs[3]
#define BLKSTART_IN prhs[4]
#define NPARIN 5
/* ************************************************************
PROCEDURE adendotd
INPUT
aden - dense.A(:,dense.l+1:end) sparse m x (nq+nden) matrix,
with aden.jc[0] possibly nonzero.
adotd - sparse m x nq array, has ai[k]'*d[k] for k in q, where
the ai's are the sparse part of the A-matrix. We still need to
add contribution from dense part, resulting in Ad.
d1 - length |K.q| vector. We will use d1(q) entries.
d2 - length firstQ+(sym(K.q)-|K.q|) vector. We use entries d2(dencols),
where dencols >=firstQ.
q - length nq array: dense lorentz blocks
dencols - length nden array: dense lorentz norm-bound columns. These
are global subscripts, at or beyond firstQ.
blkend - length nq array, listing 1-beyond-last subscript of Lorentz
norm bound blocks listed in q.
nq - number of dense lorentz blocks
nden - number of dense lorentz norm-bound columns
fwork - length m vector.
UPDATED
ad - sparse m x nq. ad.ir and ad.jc are INPUTS, ad.pr is OUTPUT.
On output, has (ai[k]+Adeni[k])'*d[k] for k in q.
************************************************************ */
void adendotd(jcir ad,jcir adotd,jcir aden,const double *d1,const double *d2,
const mwIndex *q,const mwIndex *dencols,
const mwIndex *blkend,const mwIndex nq,const mwIndex nden, double *fwork)
{
mwIndex inz, i,j,k;
const mwIndex *aden2jc;
double dj;
/* ------------------------------------------------------------
Initialize (Lorentz norm-bound part):
1) aden2jc(0:nden) points to dense columns
2) j is next dense column to handle, inz point to next nonzero
------------------------------------------------------------ */
j = 0;
aden2jc = aden.jc + nq; /* jump over Lorentz trace columns*/
inz = aden2jc[j];
for(k = 0; k < nq; k++){
/* ------------------------------------------------------------
Set fwork = all-0;
------------------------------------------------------------ */
for(i = ad.jc[k]; i < ad.jc[k+1]; i++) /* fwork = all-0 */
fwork[ad.ir[i]] = 0.0;
/* ------------------------------------------------------------
Let fwork = adotd(:,k) (Contribution from SPARSE part of A)
------------------------------------------------------------ */
for(i = adotd.jc[k]; i < adotd.jc[k+1]; i++)
fwork[adotd.ir[i]] = adotd.pr[i];
/* ------------------------------------------------------------
Let fwork += d1(q(k)) * aden(:,k) (Contribution Lorentz-trace)
------------------------------------------------------------ */
dj = d1[q[k]];
for(i = aden.jc[k]; i < aden.jc[k+1]; i++)
fwork[aden.ir[i]] += dj * aden.pr[i];
/* ------------------------------------------------------------
Add contribution of dense Lorentz-norm-bound columns, i.e.
let fwork += sum_j{d2(dencols[j]) * Aden(:,j) | dencols[j]<blkend[k]}
------------------------------------------------------------ */
for(; j < nden; j++){
if((i = dencols[j]) >= blkend[k])
break; /* Break if beyond block k */
dj = d2[i];
for(; inz < aden2jc[j+1]; inz++)
fwork[aden.ir[inz]] += dj * aden.pr[inz];
}
/* ------------------------------------------------------------
Store ad(:,k) = fwork
------------------------------------------------------------ */
for(i = ad.jc[k]; i < ad.jc[k+1]; i++)
ad.pr[i] = fwork[ad.ir[i]];
}
}
/* ============================================================
MAIN: MEXFUNCTION
============================================================ */
/* ************************************************************
PROCEDURE mexFunction - Entry for Matlab
************************************************************ */
void mexFunction(int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[])
{
const mxArray *MY_FIELD;
mwIndex i,j,firstQ, m,nden, nl, nq, lorN;
mwIndex *q, *dencols, *blkend;
const double *d1, *d2, *qPr, *dencolsPr, *blkstartPr;
double *fwork;
jcir ad, aden,adotd;
/* ------------------------------------------------------------
Check for proper number of arguments
------------------------------------------------------------ */
mxAssert(nrhs >= NPARIN, "adendotd requires more input arguments");
mxAssert(nlhs <= NPAROUT, "adendotd produces less output arguments");
/* ------------------------------------------------------------
DISASSEMBLE dense structure: dense.{cols,l,q,A}
------------------------------------------------------------ */
mxAssert(mxIsStruct(DENSE_IN),"dense should be a structure.");
MY_FIELD = mxGetField(DENSE_IN,(mwIndex)0,"l"); /* dense.l */
mxAssert( MY_FIELD != NULL, "Missing field dense.l.");
nl = (mwIndex) mxGetScalar(MY_FIELD); /* double to mwIndex */
MY_FIELD = mxGetField(DENSE_IN,(mwIndex)0,"q"); /* dense.q */
mxAssert( MY_FIELD != NULL, "Missing field dense.q.");
nq = mxGetM(MY_FIELD) * mxGetN(MY_FIELD);
qPr = mxGetPr(MY_FIELD);
MY_FIELD = mxGetField(DENSE_IN,(mwIndex)0,"cols"); /* dense.cols */
mxAssert( MY_FIELD != NULL, "Missing field dense.cols.");
nden = mxGetM(MY_FIELD) * mxGetN(MY_FIELD) - nl - nq;
mxAssert(nden >= 0, "dense.q size mismatch.");
dencolsPr = mxGetPr(MY_FIELD) + nl + nq; /* Skip LP and Q-tr*/
MY_FIELD = mxGetField(DENSE_IN,(mwIndex)0,"A"); /* dense.A */
mxAssert( MY_FIELD != NULL, "Missing field dense.A.");
mxAssert(mxIsSparse(MY_FIELD), "dense.A must be sparse");
m = mxGetM(MY_FIELD);
mxAssert(mxGetN(MY_FIELD) - nl == nq + nden, "dense.A size mismatch");
aden.jc = mxGetJc(MY_FIELD) + nl; /* Skip LP part */
aden.ir = mxGetIr(MY_FIELD);
aden.pr = mxGetPr(MY_FIELD);
/* ------------------------------------------------------------
DISASSEMBLE d structure: d.{q1,q2}
------------------------------------------------------------ */
mxAssert(mxIsStruct(D_IN), "d should be a structure.");
MY_FIELD = mxGetField(D_IN,(mwIndex)0,"q1"); /* d.q1 */
mxAssert( MY_FIELD != NULL, "Missing field d.q1.");
lorN = mxGetM(MY_FIELD) * mxGetN(MY_FIELD);
d1 = mxGetPr(MY_FIELD);
MY_FIELD = mxGetField(D_IN,(mwIndex)0,"q2"); /* d.q2 */
mxAssert( MY_FIELD != NULL, "Missing field d.q2.");
d2 = mxGetPr(MY_FIELD);
/* ------------------------------------------------------------
Get inputs adotd (contains Ad from sparse A in dense.qs blocks),
blkstart (partitions d2 into Lorentz norm-bound blocks)
------------------------------------------------------------ */
mxAssert(mxIsSparse(ADOTD_IN), "sparAD must be sparse"); /* adotd */
mxAssert((m == mxGetM(ADOTD_IN) || nq <= 0) && nq == mxGetN(ADOTD_IN), "Size mismatch sparAD");
adotd.jc = mxGetJc(ADOTD_IN);
adotd.ir = mxGetIr(ADOTD_IN);
adotd.pr = mxGetPr(ADOTD_IN);
blkstartPr = mxGetPr(BLKSTART_IN); /* blkstart */
mxAssert(lorN +1 == mxGetM(BLKSTART_IN) * mxGetN(BLKSTART_IN), "blkstart size mismatch");
/* ------------------------------------------------------------
Create working arrays q(nq), dencols(nden), fwork(m),
blkend(nq)
------------------------------------------------------------ */
q = (mwIndex *) mxCalloc(MAX(1,nq), sizeof(mwIndex));
dencols = (mwIndex *) mxCalloc(MAX(1,nden), sizeof(mwIndex));
blkend = (mwIndex *) mxCalloc(MAX(1,nq), sizeof(mwIndex));
fwork = (double *) mxCalloc(MAX(m,1), sizeof(double));
/* ------------------------------------------------------------
Convert to integer C-style; dencols, q, blkstart(q+1)
------------------------------------------------------------ */
for(i = 0; i < nden; i++){
j = (mwIndex) dencolsPr[i];
mxAssert(j>0,"");
dencols[i] = --j;
}
for(i = 0; i < nq; i++){
j = (mwIndex) qPr[i];
mxAssert(j>0,"");
q[i] = --j;
}
/* ------------------------------------------------------------
Let firstQ point to subscript of 1st Lorentz norm-bound variable
------------------------------------------------------------ */
firstQ = (mwIndex) blkstartPr[0]; /* double to mwIndex */
mxAssert(firstQ>0,"");
--firstQ; /* Fortran to C */
for(i = 0; i < nq; i++){
j = (mwIndex) blkstartPr[q[i] + 1]; /* F-double to C-mwIndex */
mxAssert(j>0,"");
blkend[i] = --j;
}
/* ------------------------------------------------------------
Create output: Ad = Ablk
------------------------------------------------------------ */
AD_OUT = mxDuplicateArray(ABLK_IN); /* Ad = Ablk */
ad.jc = mxGetJc(AD_OUT);
ad.ir = mxGetIr(AD_OUT);
ad.pr = mxGetPr(AD_OUT);
/* ------------------------------------------------------------
The real job is done here:
------------------------------------------------------------ */
adendotd(ad,adotd,aden,d1,d2 - firstQ,q,dencols,blkend,nq,nden, fwork);
/* ------------------------------------------------------------
Release working arrays
------------------------------------------------------------ */
mxFree(fwork);
mxFree(dencols);
mxFree(q);
}

Event Timeline