/*
   x = psdframeit(lab,frms,K);

   Computes x = FRM*lab.
   The frame "frms" is a product-form Householder reflection. 

    This file is part of SeDuMi 1.03BETA
    Copyright (C) 1999 Jos F. Sturm
    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., 675 Mass Ave, Cambridge, MA 02139, USA.

*/
#include <math.h>
#include <string.h>
#include "mex.h"
#include "blksdp.h"
#include "reflect.h"

#define X_OUT plhs[0]
#define NPAROUT 1

#define LAB_IN prhs[0]
#define FRMS_IN prhs[1]
#define K_IN prhs[2]
#define NPARIN 3

/* ************************************************************
   PROCEDURE psdframeit
   INPUT
     frms - lenud+hlen eigenvector basis/ orthogonal matrix Q,
       in Householder form.
     lab - sdplen vector of eigenvalues
     sdpNL sdpN-array: order of each PSD block, i.e. K.s.
     rsdpN number of real symmetric blocks: 0<= rsdpN <= sdpN.
     sdpN number of PSD blocks, i.e. length(K.s).
   UPDATED
     x - lenud vector. On input zeros. On output, x = Q*LAB*Q'.
   WORKING ARRAY:
     fwork length max(rmaxn,2*hmaxn).
   ************************************************************ */
void psdframeit(double *x, const double *frms, const double *lab,
                const int *sdpNL,const int rsdpN,const int sdpN,
                double *fwork)
{
  int k,nk,nksqr, i,inz;
  const double *beta;
  for(k = 0; k < rsdpN; k++){                /* real symmetric */
    nk = sdpNL[k];
    nksqr = SQR(nk);
    for(i = 0, inz = 0; i < nk; i++, inz += nk+1)  /* Let X = diag(lab) */
      x[inz] = lab[i];
    qtxq(x, frms + nksqr - nk, frms, nk, fwork);
    tril2sym(x,nk);
    x += nksqr; frms += nksqr;
    lab += nk;
  }
/* ------------------------------------------------------------
   For complex Hermitian, we have X = Qb' * diag(lab) * Qb,
   where Qb = Q_1 * ... * Q_{n-1} * diag(q(:,n)),
   with q(:,n) a complex sign-vector.
   Format: [RE q(:,1:n-1), RE q(:,n), IM q(:,1:n-1), IM q(:,n), beta].
   ------------------------------------------------------------ */
  for(; k < sdpN; k++){                    /* complex Hermitian */
    nk = sdpNL[k];
    nksqr = SQR(nk);
    for(i = 0, inz = 0; i < nk; i++, inz += nk+1)  /* Let X = diag(lab) */
      x[inz] = lab[i];
    beta = frms + 2*nksqr;                   /* beta = frms(:,2*n+1) */
    prpiqtxq(x,x+nksqr, beta, frms,frms+nksqr, nk, fwork);
    tril2herm(x,x+nksqr, nk);
    nksqr += nksqr;
    x += nksqr; frms += nksqr + nk;          /* skip also beta. */
    lab += nk;
  }
}

/* ============================================================
   MAIN: MEXFUNCTION
   ============================================================ */
/* ************************************************************
   PROCEDURE mexFunction - Entry for Matlab
   ************************************************************ */
void mexFunction(const int nlhs, mxArray *plhs[],
  const int nrhs, const mxArray *prhs[])
{
  int i,lendiag, lenfull, lenud,qsize;
  double *x, *fwork;
  const double *lab,*frms;
  int *sdpNL;
  coneK cK;
/* ------------------------------------------------------------
   Check for proper number of arguments 
   ------------------------------------------------------------ */
  if(nrhs < NPARIN)
    mexErrMsgTxt("psdframeit requires more input arguments.");
  if(nlhs > NPAROUT)
    mexErrMsgTxt("psdframeit generates less output arguments.");
/* ------------------------------------------------------------
   Disassemble cone K structure
   ------------------------------------------------------------ */
  conepars(K_IN, &cK);
/* ------------------------------------------------------------
   Get statistics of cone K structure
   ------------------------------------------------------------ */
  lenud = cK.rDim + cK.hDim;
  qsize = lenud + cK.hLen;
  lenfull = cK.lpN +  cK.qDim + lenud;
  lendiag = cK.lpN + 2 * cK.lorN + cK.rLen + cK.hLen;
/* ------------------------------------------------------------
   Get inputs lab,frms
   ------------------------------------------------------------ */
  lab = mxGetPr(LAB_IN);
  if(mxGetM(LAB_IN) * mxGetN(LAB_IN) != cK.rLen + cK.hLen)
    if(mxGetM(LAB_IN) * mxGetN(LAB_IN) != lendiag)
      mexErrMsgTxt("lab size mismatch");
    else
      lab += cK.lpN + 2 * cK.lorN;
  if(mxGetM(FRMS_IN) * mxGetN(FRMS_IN) != qsize)
    mexErrMsgTxt("frms size mismatch");
  frms = mxGetPr(FRMS_IN);
/* ------------------------------------------------------------
   Allocate output x
   ------------------------------------------------------------ */
  X_OUT =  mxCreateDoubleMatrix(lenud, 1, mxREAL);
  x = mxGetPr(X_OUT);
/* ------------------------------------------------------------
   Allocate working array fwork(max(rmaxn,2*hmaxn))
   integer working array sdpNL(sdpN).
   ------------------------------------------------------------ */
  fwork = (double *) mxCalloc(MAX(1,MAX(cK.rMaxn,2*cK.hMaxn)),sizeof(double));
  sdpNL = (int *) mxCalloc(MAX(1,cK.sdpN), sizeof(int));
/* ------------------------------------------------------------
   double to integer
   ------------------------------------------------------------ */
  for(i = 0; i < cK.sdpN; i++)
    sdpNL[i] = cK.sdpNL[i];
/* ------------------------------------------------------------
   PSD: X = Qb' * diag(lab) * Qb,          Qb = Q_1*Q_2*..*Q_{n-1}
   where Q_k = I-ck*ck'/betak is an elementary Householder reflection.
   Format: frms = [c1, .., c_{n-1}, beta].
   VERY INEFFICIENT !!!!
   ------------------------------------------------------------ */
  psdframeit(x, frms,lab,sdpNL,cK.rsdpN,cK.sdpN,fwork);
/* ------------------------------------------------------------
   Release working array
   ------------------------------------------------------------ */
  mxFree(sdpNL);
  mxFree(fwork);
}
