/***********************************************************************
* mexquadprod.c : C mex file 
*
*   x = mexquadprod(blk,A,t); 
* 
*   Input: blk   = [n1, n2,   nk]
*          A     = [A1 ... Ak] = mx(n1+n2+...+nk) sparse matrix
*          t     = [t1' ...tk']' = (n1+n2+...+nk) dim. vector
*
*   Output: [A1*t1 A2*t2  ... Ak*tk]    
* 
* SDPT3: version 3.0
* Copyright (c) 1997 by
* K.C. Toh, M.J. Todd, R.H. Tutuncu
* Last Modified: 2 Feb 01   
***********************************************************************/

#include <math.h>
#include <mex.h>

/**********************************************************
* compute [A1*t1 , A2*t2, ... , An*tn]
**********************************************************/
void compute(double *A, int isspA, int *irA, int *jcA, double *t, 
           double *B, int *irB, int *jcB, int numblk, int *cumblk, 
           double *At, int m) 

{  int idx, i, j, jn, l, jstart, jend;
   int  k, kstart, kend, r, jm;
   double tmp; 
   
   idx=0;
   if (isspA) { 
      for (l=0; l<numblk; ++l) {  
         jstart = cumblk[l]; jend = cumblk[l+1];
         for (j=jstart; j<jend; ++j){
            kstart = jcA[j];  kend = jcA[j+1];
            tmp = t[j];
            for (k=kstart; k<kend; ++k) {
               r = irA[k];
	       At[r] += tmp*A[k]; }  
         }
         for (k=0; k<m; ++k) {
            tmp = At[k];
	    if (tmp != 0) {
               B[idx] = tmp; 
               irB[idx] = k;
               At[k] = 0;  
               idx++; }
         }
         jcB[l+1] = idx;    
      }
   }
   else { 
      for (l=0; l<numblk; ++l) {  
         jstart = cumblk[l]; jend = cumblk[l+1];
         for (j=jstart; j<jend; ++j){
	    jm = j*m; 
            tmp = t[j];
            for (k=0; k<m; ++k) {
	       At[k] += A[k+jm]*tmp; }  
         }
         for (k=0; k<m; ++k) {
            tmp = At[k];
	    if (tmp != 0) {
               B[idx] = tmp; 
               irB[idx] = k;
               At[k] = 0;  
               idx++; }
         }
         jcB[l+1] = idx;    
      }
   }
   return;
}
/**********************************************************
* 
***********************************************************/
void mexFunction(
      int nlhs,   mxArray  *plhs[], 
      int nrhs,   const mxArray  *prhs[] )

{    mxArray  *blk_cell_pr;
     double   *A,  *B, *t,  *blksize, *At, *ttmp; 
     int      *irA, *jcA, *irB, *jcB, *irt, *jct;
     int      *cumblk;
     int      mblk, NZmax, isspA, isspt;

     int   subs[2];
     int   nsubs=2; 
     int   index, m, n, l, numblk, cols, k, r;

/* CHECK FOR PROPER NUMBER OF ARGUMENTS */

   if (nrhs < 3){
      mexErrMsgTxt("mexquadprod: requires 3 input arguments."); }
   if (nlhs > 1){ 
      mexErrMsgTxt("mexquadprod: requires 1 output argument."); }

/* CHECK THE DIMENSIONS */

    mblk = mxGetM(prhs[0]); 
    if (mblk > 1) { 
       mexErrMsgTxt("mexquadprod: blk can have only 1 row."); }   
    if (!mxIsCell(prhs[0])) {
        mexErrMsgTxt("mexqops: blk must be a cell array."); }
    subs[0]=0; subs[1]=1; 
    index = mxCalcSingleSubscript(prhs[0],nsubs,subs); 
    blk_cell_pr = mxGetCell(prhs[0],index);
    numblk = mxGetN(blk_cell_pr); 
    blksize = mxGetPr(blk_cell_pr);
    cumblk = mxCalloc(numblk+1,sizeof(int)); 
    for (l=0; l<numblk; l++) { 
        cols = (int)blksize[l]; 
        cumblk[l+1] = cumblk[l] + cols;  
    } 
    /***** assign pointers *****/
    m = mxGetM(prhs[1]); 
    n = mxGetN(prhs[1]); 
    A = mxGetPr(prhs[1]); 
    isspA = mxIsSparse(prhs[1]); 
    if (isspA) {
       irA = mxGetIr(prhs[1]); 
       jcA = mxGetJc(prhs[1]); 
       NZmax = mxGetNzmax(prhs[1]); }
    else {
       NZmax = m*n; }
    isspt = mxIsSparse(prhs[2]); 
    if (isspt) { 
       irt = mxGetIr(prhs[2]); 
       jct = mxGetJc(prhs[2]); 
       ttmp = mxGetPr(prhs[2]); 
       t = mxCalloc(n,sizeof(double)); 
       for (k=jct[0]; k<jct[1]; ++k) { 
	  r = irt[k]; t[r] = ttmp[k]; }
    }
    else { 
       t = mxGetPr(prhs[2]); }
    if (mxGetM(prhs[2]) != n) {
       mexErrMsgTxt("mexquadprod: dim not compatible."); }
    /***** create return argument *****/
     plhs[0] = mxCreateSparse(m,numblk,NZmax,mxREAL); 
     B = mxGetPr(plhs[0]);
     irB = mxGetIr(plhs[0]); 
     jcB = mxGetJc(plhs[0]); 
     jcB[0] = 0; 
     At = mxCalloc(m,sizeof(double));
     /***** Do the computations in a subroutine *****/
     compute(A,isspA,irA,jcA,t,B,irB,jcB,numblk,cumblk,At,m); 
     mxFree(At); mxFree(cumblk);  
     return;
 }
/**********************************************************/

