Page MenuHomec4science
Files F122103116

bwdpr1.c
No OneTemporary
Actions

File Metadata

Created
Tue, Jul 15, 19:18

bwdpr1.c
View Options

/*
% y = bwdpr1(Lden, b)
% BWDPR1 Solves "PROD_k L(pk,betak)' * y = b", where
% L(p,beta) = eye(n) + tril(p*beta',-1).
%
% SEE ALSO sedumi, dpr1fact, fwdpr1
% ********** INTERNAL FUNCTION OF SEDUMI **********
function y = bwdpr1(Lden, 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 <math.h>
#include <string.h>
#include "mex.h"
#include "blksdp.h"
/* y = bwdpr1fact(Lden,b) */
#define Y_OUT plhs[0]
#define NPAROUT 1
#define LDEN_IN prhs[0]
#define B_IN prhs[1]
#define NPARIN 2
/* ------------------------------------------------------------
PROCEDURE bwipr1 - I(dentity) P(lus) R(ank)1 forward solve.
INPUT:
p - length m floats
beta - length n floats
m, n - order of p and beta, resp. (n <= m)
UPDATED:
y - Length m. On input, contains the rhs. On output, the solution to
L(p,beta)'*yNEW = yOLD. This updates only y(0:n-1).
------------------------------------------------------------ */
void bwipr1(double *y, const double *p, const double *beta,
const mwIndex m, const mwIndex n)
{
mwIndex i;
double yi,t;
if(n < 1) /* If L = I, y remains the same */
return;
/* ------------------------------------------------------------
Let t = p(n:m-1)' * y
------------------------------------------------------------ */
t = realdot(p+n, y+n, m-n);
/* ------------------------------------------------------------
Solve yi for i = n-1:-1:0, from
(eye(m)+triu(beta*p',1)) * yNEW = yOLD,
i.e. yNEW(i) + betai*t = yOLD(i), with t := p(i+1:n-1)'*y.
------------------------------------------------------------ */
for(i = n; i > 0; i--){
yi = (y[i-1] -= t * beta[i-1]);
t += p[i-1] * yi;
}
}
/* ------------------------------------------------------------
PROCEDURE bwipr1o - I(dentity) P(lus) R(ank)1 forward solve, O(rdered).
INPUT:
perm - length m permutation on p and y.
p - length m floats
beta - length n floats
m, n - order of p and beta, resp. (n <= m)
UPDATED:
y - Length m. On input, contains the rhs. On output, the solution to
L(p,beta)'*yNEW = yOLD. This updates only y(0:n-1).
------------------------------------------------------------ */
void bwipr1o(double *y, const mwIndex *perm, const double *p, const double *beta,
const mwIndex m, const mwIndex n)
{
mwIndex i, permi;
double yi,t;
if(n < 1) /* If L = I, y remains the same */
return;
/* ------------------------------------------------------------
Let t = p(perm(n:m-1))' * y(perm(n:m-1))
------------------------------------------------------------ */
for(t = 0.0, i = m-1; i >= n; i--){
permi = perm[i];
t += p[permi] * y[permi];
}
/* ------------------------------------------------------------
Solve yi for i = n-1:-1:0, from
(eye(m)+triu(beta*p',1)) * yNEW = yOLD,
i.e. yNEW(i) + betai*t = yOLD(i), with t := p(i+1:n-1)'*y.
------------------------------------------------------------ */
for(i=n; i > 0; i--){
permi = perm[i-1];
yi = (y[permi] -= t * beta[i-1]);
t += p[permi] * yi;
}
}
/* ************************************************************
PROCEDURE bwprodform - Solves (PROD_j L(pj,betaj))' * yNEW = yOLD.
INPUT
p - nonzeros of sparse m x n matrix P. Has xsuper(j+1) nonzeros in
column j.
xsuper - xsuper(j+1) is number of nonzeros in p(:,j).
perm - lists pivot order for columns where ordered(j)==1.
ordered - ordered[j]==1 iff p(:,j) and beta(L:,j) have been reordered;
the original row numbers are in perm(:,j).
n - number of columns in p, beta.
UPDATED
y - length m vector. On input, the rhs. On output the solution to
(PROD_j L(pj,betaj))' * yNEW = yOLD.
************************************************************ */
void bwprodform(double *y, const mwIndex *xsuper, const mwIndex *perm,
const double *p, const double *beta, const mwIndex *betajc,
const char *ordered, const mwIndex n, mwIndex pnnz,
mwIndex permnnz)
{
mwIndex k,nk, mk;
/* ------------------------------------------------------------
Backward solve L(pk,betak) * yNEXT = yPREV for k=n-1:-1:0.
------------------------------------------------------------ */
for(k = n; k > 0; k--){
mk = xsuper[k];
nk = betajc[k] - betajc[k-1];
pnnz -= mk;
if(ordered[k-1]){
permnnz -= mk;
bwipr1o(y, perm+permnnz, p+pnnz, beta+betajc[k-1], mk, nk);
}
else
bwipr1(y, p+pnnz, beta+betajc[k-1], mk, nk);
}
mxAssert(pnnz == 0,"");
mxAssert(permnnz == 0 || permnnz == 1,"");
}
/* ============================================================
MAIN: MEXFUNCTION
============================================================ */
/* ************************************************************
PROCEDURE mexFunction - Entry for Matlab
************************************************************ */
void mexFunction(const int nlhs, mxArray *plhs[],
const int nrhs, const mxArray *prhs[])
{
const mxArray *MY_FIELD;
char *ordered;
mwIndex m,n,nden,dznnz, i,j, permnnz, pnnz;
const double *beta, *betajcPr, *orderedPr, *pivpermPr, *p;
mwIndex *betajc, *pivperm;
double *y, *fwork;
jcir dz;
/* ------------------------------------------------------------
Check for proper number of arguments
------------------------------------------------------------ */
mxAssert(nrhs >= NPARIN, "bwdpr1 requires more input arguments.");
mxAssert(nlhs <= NPAROUT, "bwdpr1 generates less output arguments.");
/* ------------------------------------------------------------
Get input b
------------------------------------------------------------ */
mxAssert(!mxIsSparse(B_IN), "b should be full");
m = mxGetM(B_IN);
n = mxGetN(B_IN);
/* ------------------------------------------------------------
Create output y as a copy of b
------------------------------------------------------------ */
Y_OUT = mxDuplicateArray(B_IN); /* Y_OUT = B_IN */
y = mxGetPr(Y_OUT);
/* ------------------------------------------------------------
Disassemble dense-update structure Lden (p,xsuper,beta,betajc,rowperm)
NOTE: if there are no dense columns, then simply let y = b, and return.
------------------------------------------------------------ */
mxAssert(mxIsStruct(LDEN_IN), "Parameter `Lden' should be a structure.");
MY_FIELD = mxGetField(LDEN_IN,(mwIndex)0,"betajc"); /* betajc */
mxAssert( MY_FIELD != NULL, "Missing field Lden.betajc.");
nden = mxGetM(MY_FIELD) * mxGetN(MY_FIELD) - 1;
/* If no dense columns: get out immediately ! */
if(nden == 0)
return;
else{
betajcPr = mxGetPr(MY_FIELD);
MY_FIELD = mxGetField(LDEN_IN,(mwIndex)0,"p"); /* p */
mxAssert( MY_FIELD != NULL, "Missing field Lden.p.");
p = mxGetPr(MY_FIELD);
pnnz = mxGetM(MY_FIELD) * mxGetN(MY_FIELD);
MY_FIELD = mxGetField(LDEN_IN,(mwIndex)0,"dopiv"); /* dopiv */
mxAssert( MY_FIELD != NULL, "Missing field Lden.dopiv.");
mxAssert(mxGetM(MY_FIELD) * mxGetN(MY_FIELD) == nden, "Size mismatch Lden.dopiv.");
orderedPr = mxGetPr(MY_FIELD);
MY_FIELD = mxGetField(LDEN_IN,(mwIndex)0,"pivperm"); /* pivperm */
mxAssert( MY_FIELD != NULL, "Missing field Lden.pivperm.");
pivpermPr = mxGetPr(MY_FIELD);
permnnz = mxGetM(MY_FIELD) * mxGetN(MY_FIELD);
MY_FIELD = mxGetField(LDEN_IN,(mwIndex)0,"beta"); /* beta */
mxAssert( MY_FIELD != NULL, "Missing field Lden.beta.");
beta = mxGetPr(MY_FIELD);
MY_FIELD = mxGetField(LDEN_IN,(mwIndex)0,"dz"); /* dz */
mxAssert( MY_FIELD != NULL, "Missing field Lden.dz.");
mxAssert(mxGetM(MY_FIELD) == m && mxGetN(MY_FIELD) == nden, "Lden.dz size mismatch.");
mxAssert(mxIsSparse(MY_FIELD), "Lden.dz must be sparse.");
dz.jc = mxGetJc(MY_FIELD);
dz.ir = mxGetIr(MY_FIELD); /* (rowperm) */
dznnz = dz.jc[nden];
mxAssert(dznnz <= m, "Lden.dz size mismatch.");
/* ------------------------------------------------------------
Allocate working arrays mwIndex: betajc(nden+1), pivperm(permnnz),
char: ordered(nden)
double: fwork(dznnz)
------------------------------------------------------------ */
betajc = (mwIndex *) mxCalloc(nden + 1,sizeof(mwIndex));
pivperm = (mwIndex *) mxCalloc(MAX(permnnz,1),sizeof(mwIndex));
ordered = (char *) mxCalloc(nden,sizeof(char)); /* nden > 0 */
fwork = (double *) mxCalloc(MAX(dznnz,1), sizeof(double));
/* ------------------------------------------------------------
Convert betajcPr, ordered, pivperm to mwIndex
------------------------------------------------------------ */
for(i = 0; i <= nden; i++){
j = betajcPr[i];
betajc[i] = --j;
}
for(i = 0; i < nden; i++){
ordered[i] = orderedPr[i];
}
for(i = 0; i < permnnz; i++){
pivperm[i] = pivpermPr[i];
}
MY_FIELD = mxGetField(LDEN_IN,(mwIndex)0,"beta");
mxAssert(mxGetM(MY_FIELD) * mxGetN(MY_FIELD) == betajc[nden], "Size mismatch Lden.beta.");
/* ------------------------------------------------------------
The actual job is done here: y(perm,:) = PROD_L'\b(perm,:).
------------------------------------------------------------ */
for(j = 0; j < n; j++, y += m){
for(i = 0; i < dznnz; i++) /* fwork = y(dzir) */
fwork[i] = y[dz.ir[i]];
bwprodform(fwork, dz.jc, pivperm, p, beta, betajc, ordered, nden,
pnnz, permnnz);
for(i = 0; i < dznnz; i++) /* y(dzir) = fwork */
y[dz.ir[i]] = fwork[i];
}
/* ------------------------------------------------------------
Release working arrays
------------------------------------------------------------ */
mxFree(fwork);
mxFree(ordered);
mxFree(pivperm);
mxFree(betajc);
} /* if dense columns */
}

Event Timeline

c4science · Help