qreshape.c

经济学专业代码

        k = ipos[inz+1] - 1;
        iwork[knz] = k;
        ipos[knz] = inz;               /* start of block k in yir */
        cfound[knz++] = 0;
      }
    ipos[knz] = ynnz;
    blknnz = knz;
  }
/* ------------------------------------------------------------
   Update all subscripts of norm-bound entries: += lorN-(k+1).
   ------------------------------------------------------------ */
  for(knz = 0; knz < blknnz; knz++){
    k = iwork[knz];
    if(++k < lorN){
      inz = ipos[knz] + cfound[knz];
      intadd(yir + inz, lorN-k, ipos[knz+1] - inz);
    }
  }
/* ------------------------------------------------------------
   Condense to a list of nonzero tr-entries
   ------------------------------------------------------------ */
  for(knz = 0; knz < blknnz; knz++)
    if(!cfound[knz])
      break;
  for(k = knz+1; k < blknnz; k++)
    if(cfound[k]){
      iwork[knz] = iwork[k];    /* block number (= new subscript) */
      ipos[knz] = ipos[k];      /* (old) position in yir */
      knz++;
    }
/* ------------------------------------------------------------
   Let blknnz be the number of nonzero trace entries
   ------------------------------------------------------------ */
  blknnz = knz;
  if(blknnz > 0){
/* ------------------------------------------------------------
   temporarily store nonzeros of tr-entries in fwork, then
   move other downzeros downwards to fill-up the gaps
   ------------------------------------------------------------ */
    for(knz = 0; knz < blknnz; knz++)
      fwork[knz] = ypr[ipos[knz]];
    for(knz = 1; knz < blknnz; knz++){
      inz = ipos[blknnz - knz - 1] + 1;     /* just beyond tr-entry */
      k = ipos[blknnz - knz] - inz;         /* amount before next tr-entry */
      if(k > 0){
        memmove(yir+inz+knz, yir+inz, k * sizeof(int));
        memmove(ypr+inz+knz, ypr+inz, k * sizeof(double));
      }
    }
    if(ipos[0] > 0){                   /* nonzeros before 1st nz-tr entry */
      memmove(yir+blknnz, yir, ipos[0] * sizeof(int));
      memmove(ypr+blknnz, ypr, ipos[0] * sizeof(double));
    }
/* ------------------------------------------------------------
   Re-insert the (saved) nonzero tr-entries at the start
   ------------------------------------------------------------ */
    memcpy(ypr, fwork, blknnz * sizeof(double));
    memcpy(yir, iwork, blknnz * sizeof(int));
    intadd(yir, blks[0], blknnz);
  }
}

/* ************************************************************
   PROCEDURE inverse of qreshape0, Transforms from [x1(lorN); x2-blocks]
     to pure block-wise [x[1], x[2],.. x[lorN]].
   INPUT
     blks - lorN array; blks[k]:blk[k+1]-1 are the subscripts of Lorentz
        block k=0:lorN-2.
     lorN - number of Lorentz blocks
   UPDATED
     y - entries in y(blk[0]:blks[lorN-1]) are affected (reshuffled)
   WORK
     fwork - length lorN-1 vector.
   ************************************************************ */
void qreshape1(double *y, const int *blks, const int lorN, double *fwork)
{
  int i,j,k;
  if(lorN <=1)
    return;             /* Nothing to do if only 1 block */
/* ------------------------------------------------------------
   Save y(2:lorN), which are the trace ("x1")-entries of block 2:lorN.
   ------------------------------------------------------------ */
  memcpy(fwork, y+blks[0] + 1, (lorN-1) * sizeof(double));
/* ------------------------------------------------------------
   For each block k = 0:lorN-2, shift y2[k] (lorN-1)-k positions upwards.
   ------------------------------------------------------------ */
  j = lorN;
  for(k = 1; k < lorN; k++){
    --j;                          /* lorN-1 >= j = lorN-k >= 1 */
    i = blks[k-1]+1;              /* new start of block y2[k-1] */
    memmove(y+i,y+i+j, (blks[k]-i) * sizeof(double));
  }
/* ------------------------------------------------------------
   Let y(iTr) = x(1:lorN), where iTr = cumsum([1 lorNL]) = blks.
   ------------------------------------------------------------ */
  for(k = 1; k < lorN; k++)
    y[blks[k]] = fwork[k-1];       /* y[blks[0]] is left unaffected */
}

/* ************************************************************
   PROCEDURE inverse of spqreshape0, Transforms from [x1(lorN); x2-blocks]
     to pure block-wise [x[1], x[2],.. x[lorN]].
   INPUT
     ynnz - nnz(y)
     blks - lorN+1 array; blks[0]:blks[1]-1 are the subscripts of the lorN
        trace values "y1[0:lorN-1]". Blks[k+1]:blk[k+2]-1 are the subscripts
        of Lorentz block y2[k], k = 0:lorN-2.
     lorN - number of Lorentz blocks
     iwsize - length of iwork, should be 2*(1+lorN).
   UPDATED
     yir, ypr - entries with subscripts in y(blk[0]:blks[lorN]) are
       affected (reshuffled)
   WORK
     cfound - length lorN character working array.
     iwork - length iwsize = 2 + lorN + MAX(floor(log(lorN+1)/log(2)), lorN-1),
        which is at most 2*(1+lorN).
     fwork - length lorN-1 vector.
   ************************************************************ */
/* still to be implemented */

/* ============================================================
   MAIN: MEXFUNCTION
   ============================================================ */
/* ************************************************************
   PROCEDURE mexFunction - Entry for Matlab
   ************************************************************ */
void mexFunction(const int nlhs, mxArray *plhs[],
  const int nrhs, const mxArray *prhs[])
{
 int i,j, ifirst, N,m, iwsize;
 coneK cK;
 int *iwork, *blks;
 double *fwork;
 char *cwork;
 bool flag, bsparse;
 jcir y;

/* ------------------------------------------------------------
   Check for proper number of arguments 
   ------------------------------------------------------------ */
 mxAssert(nrhs >= NPARIN, "qreshape requires more input arguments.");
 mxAssert(nlhs <= NPAROUT, "qreshape generates 1 output argument.");
/* ------------------------------------------------------------
   Disassemble cone K structure
   ------------------------------------------------------------ */
 conepars(K_IN, &cK);
/* ------------------------------------------------------------
   Get inputs x, flag
   ------------------------------------------------------------ */
 N = mxGetM(X_IN);                    /* dimemsion */
 m = mxGetN(X_IN);                    /* number of cols */
 ifirst = 0;                          /* 1st Lorentz index */
 if(mxGetM(X_IN) != cK.qDim){
   mxAssert(mxGetM(X_IN) >= cK.lpN + cK.qDim, "x size mismatch");
   ifirst = cK.lpN;
 }
/* ------------------------------------------------------------
   Allocate output y(N,m)
   ------------------------------------------------------------ */
 Y_OUT = mxDuplicateArray(X_IN);                   /* Y = X_IN */
 if((bsparse = mxIsSparse(Y_OUT))){
   y.jc = mxGetJc(Y_OUT);
   y.ir = mxGetIr(Y_OUT);
 }
 y.pr = mxGetPr(Y_OUT);
 flag = mxGetScalar(FLAG_IN);
/* ------------------------------------------------------------
   Allocate working array iwork(2*(1+lorN)), fwork(lorN), blks(lorN),
   and cwork(lorN).
   ------------------------------------------------------------ */
 iwsize = 2 + 2 * cK.lorN;
 iwork = (int *) mxCalloc(iwsize, sizeof(int));
 fwork = (double *) mxCalloc(MAX(cK.lorN,1), sizeof(double));
 blks = (int *) mxCalloc(cK.lorN+1, sizeof(int));
 for(i = 1; i <= cK.lorN; i++)
   blks[i] = cK.lorNL[i-1];           /* float to int */
 cwork = (char *) mxCalloc(MAX(1,cK.lorN), sizeof(char));
/* ------------------------------------------------------------
   Let iwork(0:lorN) = cumsum([ifirst, K.q(1:end)])
   ------------------------------------------------------------ */
 j = ifirst;
 blks[0] = j;
 for(i = 1; i <= cK.lorN; i++){
   j += blks[i];
   blks[i] = j;
 }
/* ------------------------------------------------------------
   The real job:
   ------------------------------------------------------------ */
 if(bsparse){
   mxAssert(flag != 1, "qreshape(x,1,K) cannot handle sparse x.");
#ifdef DO_ALL
   if(flag == 1)
     spqreshape1(y.pr, blks, cK.lorN, iwsize, iwork,fwork,cwork);
   else
#endif
     for(j = 0; j < m; j++)
       spqreshape0(y.ir+y.jc[j],y.pr+y.jc[j],y.jc[j+1]-y.jc[j],
                   blks, cK.lorN, iwsize, cwork, iwork,fwork);
 }
 else
   if(flag == 1)
     for(j = 0; j < m; j++)
       qreshape1(y.pr+j*N, blks, cK.lorN, fwork);
   else
     for(j = 0; j < m; j++)
       qreshape0(y.pr+j*N, blks, cK.lorN, fwork);
/* ------------------------------------------------------------
   RELEASE WORKING ARRAY
   ------------------------------------------------------------ */
 mxFree(iwork);
 mxFree(fwork);
 mxFree(cwork);
 mxFree(blks);
}