spscale.c

经济学专业代码

  xd = fwork + SQR(n);
/* ------------------------------------------------------------
   Let dxcols={j: X(:,j)!= all-0}, m = |dxcols|,
   fwork = D*X(:,dxcols).
   ------------------------------------------------------------ */
  m = realdmulx(fwork, dxcols, d, xpr, xir, pxjc0, xjc1, first, n);
/* ------------------------------------------------------------
   Let xd = fwork' = X(:,dxcols)'*D.
   ------------------------------------------------------------ */
  realtransp(xd, fwork, n,m);
/* ------------------------------------------------------------
   Let fwork = D(dxcols,:).     Both xd and fwork are m x n.
   ------------------------------------------------------------ */
  for(j = 0, inz = 0; j < n; j++, d += n)
    for(i = 0; i < m; i++)
      fwork[inz++] = d[dxcols[i]];        /* fwork(i,j) = d(dxcols(i),j) */
/* ------------------------------------------------------------
   For all target subscripts (i,j), let
   zij = (D*sym(X)*D)_ij = [ DXD_ij + DXD_ji ] /2.
   Note that DXD_ij = xd(:,i)' * fwork(:,j).   (m mults)
   ------------------------------------------------------------ */
  jlast = 0;      /* index right after last activated column */
  for(inz = 0; inz < znnz; inz++){
    if((i = zir[inz]) >= jlast){      /* move to new z-column */
      j = (i-first) / n;              /* col j = floor( (i-first)/n ) */
      if(j >= n)
        break;                        /* finished all columns of Z */
      jfirst = first + n * j;
      jlast = jfirst + n;
      icol = m * j;
      dj = fwork + icol;
      xdj = xd + icol;
      if(i - jfirst == j)             /* zjj = xdj' * dj */
        z[i] = realdot(xdj,dj, m);
      else{              /* zij = (xdj'*di + xdi'*dj) / 2 */
        icol = (i - jfirst) * m;
        z[i] = (realdot(xdj,fwork + icol, m) + realdot(xd + icol, dj, m)) / 2;
      }
    }
    else{               /* zij = (xdj'*di + xdi'*dj) / 2 */
      icol = (i - jfirst) * m;
      z[i] = (realdot(xdj,fwork + icol, m) + realdot(xd + icol, dj, m)) / 2;
    }
  }
}

/* ************************************************************
   PROCEDURE spcpxdxd - Computes Z = D*herm(X)*D with D = D' and
        herm(X) = (X+X')/2. X sparse, and Z has pre-chosen sparsity
        structure zir (typically subset of tril-entries).
      Z,D, and X are complex matrices, stored as [vec(RE); vec(IM)].
   INPUT
     zir, zjc0, zjc1 - Target sparsity structure. We compute z(i) only for
       i in zir(zjc0:zjc1-1) with i < first + 2*n^2.
       zjc0 such that zir(zjc0) >= first.
     d      - full 2*(n x n) matrix D.
     xpr, xir, xjc1 - sparse vector, xjc1 nonzeros.
     first  - Subscript of X(1,1) and Z(1,1), i.e.
       x(first:first+n^2-1) = vec(real(X)), x(first+n^2:end)=vec(imag(X)).
     n      - order of D, X.
   UPDATED
     *pxjc0 - On input, points to first nonzero with index in range
       first:first+2*n^2-1. On output, points just BEYOND.
   OUTPUT
     Z - lenfull vector. Only z(first:first+2*n^2-1) is affected, and from this
       only the indices listed in zir. (Other entries unaffected.)
   WORK
     dxcols - length n integer array.
     fwork  - 4 * n^2 vector.
   RETURNS zjc0_NEW
   ************************************************************ */
void spcpxdxd(double *z, const int *zir, const int znnz,
              const double *d,
              const double *xpr, const int *xir, int *pxjc0, const int xjc1,
              const int first, const int n, double *fwork, int *dxcols)
{
  int inz, i, icol, j, jfirst, jlast, m, nsqr, imgfirst;
  double *dx, *dxpi, *fworkpi;
  double zi;
  const double *dj, *djpi, *djx, *djxpi;
  nsqr = SQR(n);
/* ------------------------------------------------------------
   Partition 4*n^2 WORKING array fwork into [fwork(2*n^2), dxRows(2*n^2)].
   ------------------------------------------------------------ */
  fworkpi = fwork + nsqr;
  dx = fworkpi + nsqr;
  dxpi = dx + nsqr;
/* ------------------------------------------------------------
   Let dxcols={j: X(:,j)!= all-0}, m = |dxcols|,
   fwork = RE(D*X(:,dxcols)), fwork+n^2 = IM D*X(:,dxcols).
   ------------------------------------------------------------ */
  m = cpxdmulx(fwork, dxcols, d, xpr, xir, pxjc0, xjc1, first, n);
/* ------------------------------------------------------------
   Let dxRows = D*X(:,dxcols) in row-stacking (instead of usual col-stack).
   ------------------------------------------------------------ */
  realtransp(dx, fwork, n,m);
  realtransp(dxpi, fworkpi, n,m);
/* ------------------------------------------------------------
   Let fwork = D(dxcols,:) in (usual) column stacking; m x n.
   ------------------------------------------------------------ */
  for(j = 0, inz = 0; j < n; j++, d += n)
    for(i = 0; i < m; i++)
      fwork[inz++] = d[dxcols[i]];        /* fwork(i,j) = RE d(dxcols(i),j) */
  for(j = 0, inz = 0; j < n; j++, d += n)
    for(i = 0; i < m; i++)
      fworkpi[inz++] = d[dxcols[i]];     /* fworkpi(i,j) = IM d(dxcols(i),j) */
/* ------------------------------------------------------------
   For all target subscripts (i,j) in RE Z, let
   zij = RE (D*herm(X)*D)_ij = RE [ DXD_ij + DXD_ji ] /2.
   Note that DXD_ij = dx(i,:) * fwork(:,j).   (m mults)
   Since dx is in row-stacking, we can simply browse thru memory.
   ------------------------------------------------------------ */
  jlast = 0;      /* index right after last activated column */
  for(inz = 0; inz < znnz; inz++){
    if((i = zir[inz]) >= jlast){      /* move to new z-column */
      j = (i-first) / n;              /* col j = floor( (i-first)/n ) */
      if(j >= n)
        break;                        /* finished all columns of RE(Z) */
      jfirst = first + n * j;
      jlast = jfirst + n;
      icol = m * j;
      dj = fwork + icol;
      djpi = dj + nsqr;
      djx = dx + icol;
      djxpi = djx + nsqr;
      if(i - jfirst == j)             /* zjj = djx * dj - IMdjx*IMdj * */
        z[i] = realdot(djx,dj, m) - realdot(djxpi,djpi, m);
      else{              /* zij = RE (dix*dj + conj(djx*di)) / 2 */
        icol = (i - jfirst) * m;
        zi = realdot(djx,fwork + icol, m) - realdot(djxpi,fworkpi+icol,m);
        zi += realdot(dx + icol, dj, m) - realdot(dxpi+icol, djpi, m);
        z[i] = zi / 2;
      }
    }
    else{               /* zij = RE (dix*dj + conj(djx*di)) / 2 */
      icol = (i - jfirst) * m;
      zi = realdot(djx,fwork + icol, m) - realdot(djxpi,fworkpi+icol,m);
      zi += realdot(dx + icol, dj, m) - realdot(dxpi+icol, djpi, m);
      z[i] = zi / 2;
    }
  }
/* ------------------------------------------------------------
   Remains: target subscripts (i,j) in IM Z.
   zij = IM (D*herm(X)*D)_ij = IM [ DXD_ij + DXD_ji ] /2.
   Note that DXD_ij = dx(i,:) * fwork(:,j).   (m mults)
   Since dx is in row-stacking, we can simply browse thru memory.
   ------------------------------------------------------------ */
  imgfirst = first + nsqr;
  for(; inz < znnz; inz++){
    if((i = zir[inz]) >= jlast){      /* move to new z-column */
      j = (i-imgfirst) / n;              /* col j = floor( (i-first)/n ) */
      if(j >= n)
        break;                        /* finished all columns of IM(Z) */
      jfirst = imgfirst + n * j;
      jlast = jfirst + n;
      icol = m * j;
      dj = fwork + icol;
      djpi = dj + nsqr;
      djx = dx + icol;
      djxpi = djx + nsqr;
      mxAssert(i - jfirst != j,"");          /* Herm => diag(imag(Z))=0 */
      icol = (i - jfirst) * m;     /* zij = IM (dix*dj + conj(djx*di)) / 2 */
      zi = realdot(dx + icol, djpi, m) + realdot(dxpi+icol, dj, m);
      zi -= realdot(djx,fworkpi + icol, m) + realdot(djxpi,fwork+icol,m);
      z[i] = zi / 2;
    }
    else{               /* zij = IM (dix*dj + conj(djx*di)) / 2 */
      icol = (i - jfirst) * m;
      zi = realdot(dx + icol, djpi, m) + realdot(dxpi+icol, dj, m);
      zi -= realdot(djx,fworkpi + icol, m) + realdot(djxpi,fwork+icol,m);
      z[i] = zi / 2;
    }
  }
}

/* ************************************************************
   PROCEDURE spsqrscale - Computes z = D(d^2)x for PSD-part.
   INPUT
     zir, zjc1 - Target sparsity structure. We compute z(i) only for
       i in zir(0:zjc1-1).
     d      - lenud vector containing full n x n matrices Dk (=Ud'*Ud).
     xpr, xir, xjc1 - sparse vector, xjc1 nonzeros.
     xjc0 - Points to first nonzero with index in range
       blkstart[0]:blkstart[psdN].
     blkstart  - length psdN+1 array such that 
       x(blkstart[k]:blkstart[k+1]) = vec(Xk) for all PSD blocks k.
     xblk - length psdDim array, with k = xblk(i-blkstart[0]) iff
       blkstart[k] <= i < blkstart[k+1], k=0:nblk-1.
       psdDim:=blkstart[end]-blkstart[0].
     psdNL     - length psdN array, is K.s.
     rpsdN     - number of real symmetric PSD blocks.
   OUTPUT
     z - lenfull vector. Only indices listed in zir are affected.
       Other entries are unaffected (not even set to 0).
     blks - length blksnnz <= length(K.s), lists PSD-blocks where
       z has nonzeros.
   WORK
     fwork - 2 * (rMaxn^2 + 2*hMaxn^2) vector.  (<= 4*max(K.s)^2)
     iwork - length n=MAX(rMaxn,hMaxn) integer array. (n=max(K.s))
   RETURNS blksnnz, number of entries written into blks.
   ************************************************************ */
int spsqrscale(double *z, int *blks, const int *zjc, const int *zir,
               const int *znnz, const double *d,
               const int *xir, const double *xpr, int xjc0, const int xjc1,
               const int *blkstart, const int *xblk, const int *psdNL,
               const int rpsdN, double *fwork, int *iwork)
{
  int blksnnz, k;

  blksnnz = 0;
  d -= blkstart[0];        /* Now, d + blkstart[k] gives Dk */
  while(xjc0 < xjc1){
/* ------------------------------------------------------------
   For each nonzero block xblk[.] in x, Let zk = vec(tril(Dk * Xk * Dk)).
   ------------------------------------------------------------ */
    if((k = xblk[xir[xjc0]]) >= rpsdN)
      break;                                /* 1st nz Hermitian PSD block */
    sprealdxd(z, zir+zjc[k], znnz[k], d + blkstart[k], xpr, xir, &xjc0, xjc1,
              blkstart[k], psdNL[k], fwork, iwork);
    blks[blksnnz++] = k;
  }
/* ------------------------------------------------------------
   Hermitian PSD blocks
   ------------------------------------------------------------ */
  while(xjc0 < xjc1){
    k = xblk[xir[xjc0]];
    spcpxdxd(z, zir+zjc[k], znnz[k], d + blkstart[k], xpr, xir, &xjc0, xjc1,
             blkstart[k], psdNL[k], fwork, iwork);
    blks[blksnnz++] = k;
  }
  return blksnnz;
}

#ifdef OLDSEDUMI
/* =========================  P S D :  ========================= 
   D(d)x = vec(U' * X * U), only at zir positions.
   Typically, zir is subset of tril (or of triu).
   ============================================================ */

/* ************************************************************
   PROCEDURE sprealutxu - Computes Z(perm,perm) = U' * X * U.
        Z has pre-chosen sparsity structure zir
	(typically subset of tril-entries [or of triu-entries]).
   INPUT
     zir, zjc0, zjc1 - Target sparsity structure. We compute z(i) only for
       i in zir(zjc0:zjc1-1) with i < first + n^2.
       zjc0 such that zir(zjc0) >= first.
     u      - full n x n matrix, triu(u) = U.
     invperm - length n array, invperm(perm)=0:n-1. Used to translate Z-index
       to u- and x-index.
     x - n x n full symmetric matrix containing triu(X).
     first  - Subscript of X(1,1) and Z(1,1), i.e.
       x(first:first+n^2-1) = vec(X).
     n      - order of U, X.
   OUTPUT
     Z - lenfull vector. Only z(first:first+n^2-1) is affected, and from this
       only the indices listed in zir. (Other entries unaffected.)