fwdpr1.c

/*
%                                                         y = fwdpr1(Lden, b)
% FWDPR1  Solves "PROD_k L(pk,betak) * y = b", where
%    where L(p,beta) = eye(n) + tril(p*beta',-1).
%
% SEE ALSO sedumi, dpr1fact, bwdpr1
% **********  INTERNAL FUNCTION OF SEDUMI **********
function y = fwdpr1(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 = fwdpr1(Lden, b) */
#define Y_OUT plhs[0]
#define NPAROUT 1

#define LDEN_IN prhs[0]
#define B_IN prhs[1]
#define NPARIN 2

/* ************************************************************
   PROCEDURE fwprodform - 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.
     beta, betajc - Column beta(betajc[j]:betajc[j+1]-1) provides the
       beta vector in the j-th L-factor, L(p,beta)=eye(m) + tril(p*beta',-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 fwprodform(double *y, const mwIndex *xsuper, const mwIndex *perm,
                const double *p, const double *beta, const mwIndex *betajc,
                const char *ordered, const mwIndex n)
{
  mwIndex k,nk, mk;
/* ------------------------------------------------------------
   Forward solve L(pk,betak) * yNEXT = yPREV   for k=0,1,..,n-1, resp.
   ------------------------------------------------------------ */
  for(k = 1; k <= n; k++){
    mk = xsuper[k];
    nk = betajc[k] - betajc[k-1];
    if(ordered[k-1]){
      fwipr1o(y, perm, p, beta, mk, nk);
      perm += mk;
    }
    else
      fwipr1(y, p, beta, mk, nk);
    beta += nk;
    p += mk;
  }
}

/* ============================================================
   MAIN: MEXFUNCTION
   ============================================================ */
/* ************************************************************
   PROCEDURE mexFunction - Entry for Matlab
       y = fwdpr1(Lden, b)
   ************************************************************ */
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;
 const double *beta, *betajcPr, *orderedPr, *pivpermPr, *p;
 mwIndex *betajc, *pivperm;
 double *y, *fwork;
 jcir dz;
/* ------------------------------------------------------------
   Check for proper number of arguments 
   ------------------------------------------------------------ */
  mxAssert(nrhs >= NPARIN, "fwdpr1 requires more input arguments.");
  mxAssert(nlhs <= NPAROUT, "fwdpr1 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");
  mxAssert( MY_FIELD!= NULL, "Missing field Lden.betajc.");    /* 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");
    mxAssert( MY_FIELD != NULL, "Missing field Lden.p.");         /* p */
    p = mxGetPr(MY_FIELD);
    MY_FIELD = mxGetField(LDEN_IN,(mwIndex)0,"dopiv");
    mxAssert( MY_FIELD != NULL, "Missing field Lden.dopiv.");  /* 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: fwork = PROD_L\b(perm,j), y(perm,j) = fwork.
   ------------------------------------------------------------ */
    for(j = 0; j < n; j++, y += m){
      for(i = 0; i < dznnz; i++)            /* fwork = y(dzir) */
        fwork[i] = y[dz.ir[i]];
      fwprodform(fwork, dz.jc, pivperm, p, beta, betajc, ordered, nden);
      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 */
}