/* % 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 #include #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 */ }