spscale.c

经济学专业代码

   WORK
     fwork  - 2 * n^2 vector.
   RETURNS zjc0_NEW
   ************************************************************ */
int sprealutxu(double *z, const int *zir, const int zjc0, const int zjc1,
               const double *u, const int *invperm, const double *x,
               const int first, const int n, double *fwork)
{
  int inz, i, irow, j, icol, jcol, jfirst, jlast, m;
  double *xu;
  const double *uj, *xuj;
/* ------------------------------------------------------------
   Partition 2*n^2 WORKING array fwork into [fwork(n^2), xu(n^2)].
   ------------------------------------------------------------ */
  xu = fwork + SQR(n);
/* ------------------------------------------------------------
   Let fwork = X in full symmetric storage
   ------------------------------------------------------------ */
  memcpy(fwork, x, SQR(n) * sizeof(double));
  triu2sym(fwork,n);
/* ------------------------------------------------------------
   Compute xu = triu(fwork'*U):   1^2 + 2^2 + ... + (n-1)^2 + n^2 mults.
   For each j, let uj point to u(0,j); m is
   remaining length of column j, i.e. m=j.      (SEE realutxu.)
   ------------------------------------------------------------ */
  inz = 0; jcol = 0;
  for(j = 0, m = n; j < n; j++, jcol += n){
    for(i = 0, icol = 0; i <= j; i++, icol += n)
      xu[inz++] = realdot(fwork+icol,u+jcol,j+1);  /* x(:,i)'*u(0:j,j) */
    inz += --m;                                    /* skip tril */
  }
/* ------------------------------------------------------------
   For all target subscripts (i,j), let
   zij = (U'*XU)_ij [ = (U'*XU)_ji ]
   ------------------------------------------------------------ */
  jlast = 0;                  /* index right after last activated column */
  for(inz = zjc0; inz < zjc1; 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;
      j = invperm[j];                 /* change to U- and X-index */
      jcol = n * j;
      uj = u + jcol;
      xuj = xu + jcol;
    }
    irow = invperm[i - jfirst];       /* U- and X-index */
    if(irow <= j)            /* zij = u(0:i,i)' * xu(0:i,j) IF i<=j */
      z[i] = realdot(u + irow * n, xuj, irow + 1);
    else{                    /* zij = u(0:j,j)' * xu(0:j,i) IF i>j */
      z[i] = realdot(uj, xu + irow * n, j + 1);
    }
  }
/* ------------------------------------------------------------
   RETURN inz, next position in zir (after range first:first+n^2-1).
   ------------------------------------------------------------ */
  return inz;
}

/* ************************************************************
   PROCEDURE spcpxutxu - Computes Z(perm,perm) = U' * X * U.
        Z has pre-chosen sparsity structure zir
	(typically subset of tril-entries [or of triu-entries]).
        U,X and Z 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.
     u      - full 2*(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 - 2*(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(RE X), x(first+n^2:end) = vec(IM(X))..
     n      - order of U, X.
   OUTPUT
     Z - lenfull vector. Only z(first:first+2n^2-1) is affected, and from this
       only the indices listed in zir. (Other entries unaffected.)
   WORK
     fwork  - 4 * n^2 vector.
   RETURNS zjc0_NEW
   ************************************************************ */
int spcpxutxu(double *z, const int *zir, const int zjc0, const int zjc1,
              const double *u, const int *invperm, const double *x,
              const int first, const int n, double *fwork)
{
  int inz, i, irow, j, icol, jcol, jfirst, jlast, m, nsqr, imgfirst;
  double *xu, *xupi, *fworkpi;
  const double *uj, *ujpi, *xuj, *xujpi, *upi;
/* ------------------------------------------------------------
   Partition 4*n^2 WORKING array fwork into [fwork(2*n^2), xu(2*n^2)].
   ------------------------------------------------------------ */
  nsqr = SQR(n);
  upi = u + nsqr;
  fworkpi = fwork + nsqr;
  xu = fworkpi + nsqr;
  xupi = xu + nsqr;
/* ------------------------------------------------------------
   Let fwork = X in full hermitian storage
   ------------------------------------------------------------ */
  memcpy(fwork, x, 2*nsqr * sizeof(double));
  triu2herm(fwork,fworkpi,n);
/* ------------------------------------------------------------
   Compute xu = triu(fwork'*U)   ( =triu(XU), since X is Hermitian.)   
   For each j, let uj point to u(0,j); m is
   remaining length of column j, i.e. m=j.      (SEE prpiutxu.)
   ------------------------------------------------------------ */
  inz = 0; jcol = 0;
  for(j = 0, m = n; j < n; j++, jcol += n){
    for(i = 0, icol = 0; i <= j; i++, icol += n){
/* xu(i,j) = x(:,i)'*u(0:j,j), assume IM u(j,j) == 0 */
      xu[inz] = realdot(fwork+icol,u+jcol,j+1) +
        realdot(fworkpi+icol,upi+jcol,j);
      xupi[inz] = realdot(fwork+icol,upi+jcol,j) -
        realdot(fworkpi+icol,u+jcol,j+1);
      inz++;
    }
    inz += --m;                                /* skip tril */
  }
/* ------------------------------------------------------------
   For all target subscripts (i,j) in RE Z, let
   zij = RE (U'*XU)_ij [ = RE (U'*XU)_ji ]
   ------------------------------------------------------------ */
  jlast = 0;                  /* index right after last activated column */
  for(inz = zjc0; inz < zjc1; 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;
      j = invperm[j];                 /* change to U- and X-index */
      jcol = n * j;
      uj = u + jcol;
      ujpi = uj + nsqr;
      xuj = xu + jcol;
      xujpi = xuj + nsqr;
    }
    irow = invperm[i - jfirst];       /* U- and X-index */
    icol = irow * n;
    if(irow <= j)            /* zij = RE u(0:i,i)' * xu(0:i,j) IF i<=j */
      z[i] = realdot(u+icol, xuj, irow+1) + realdot(upi+icol, xujpi, irow);
    else{                    /* zij = RE u(0:j,j)' * xu(0:j,i) IF i>j */
      z[i] = realdot(uj, xu+icol, j+1) + realdot(ujpi, xupi+icol, j);
    }
  }
/* ------------------------------------------------------------
   Remains target subscripts (i,j) in IM Z. Let
   zij = IM (U'*XU)_ij [ = -IM (U'*XU)_ji ]
   ------------------------------------------------------------ */
  imgfirst = first + nsqr;
  for(; inz < zjc1; 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;
      j = invperm[j];                 /* change to U- and X-index */
      jcol = n * j;
      uj = u + jcol;
      ujpi = uj + nsqr;
      xuj = xu + jcol;
      xujpi = xuj + nsqr;
    }
    irow = invperm[i - jfirst];       /* U- and X-index */
    icol = irow * n;
    if(irow <= j)            /* zij = IM u(0:i,i)' * xu(0:i,j) IF i<=j */
      z[i] = realdot(u+icol, xujpi, irow+1) - realdot(upi+icol, xuj, irow);
    else                     /* zij = -IM u(0:j,j)' * xu(0:j,i) IF i>j */
      z[i] = realdot(ujpi, xu+icol, j) - realdot(uj, xupi+icol, j+1);
  }
/* ------------------------------------------------------------
   RETURN inz, next position in zir (after range first:first+2*n^2-1).
   ------------------------------------------------------------ */
  return inz;
}

/* ************************************************************
   PROCEDURE blkpsdscale - Computes z = D(d)x = U'*X*U for PSD-part.
     Typically, X = D(d)aj = vec(triu(U*Aj*U')).
   INPUT
     zir, zjc1 - Target sparsity structure. We compute z(i) only for
       i in zir(0:zjc1-1).
     u      - lenud vector containing full n x n matrices Uk
     invperm - lenpsd array, containing inverse perm for Uk for each psd-block
     x      - vector of full blocks, viz. those listed in xblk.
     xblk, blkjc0, blkjc1 - the psd-blocks stored in x.
     blkstart  - length psdN+1 array such that
       z(blkstart[k]:blkstart[k+1]) = vec(Zk) for all PSD blocks k.
     psdNL     - length psdN array, is K.s.
     cumpsdNL  - cumsum([0, psdNL]), i.e. start of block k in invperm.
     rpsdN     - number of real symmetric PSD blocks.
   OUTPUT (PRE-INITIALIZED)
     z - lenfull vector. Only indices listed in zir are affected.
       Other entries are unaffected (not even set to 0).
   WORK
     fwork  - 2 * (rMaxn^2 + 2*hMaxn^2) vector.  (<= 4*max(K.s)^2)
   RETURNS zjc0_NEW (should be zjc1)
   ************************************************************ */
int blkpsdscale(double *z, const int *zir, const int zjc1,
		const double *u, const int *invperm, const double *x,
		const int *xblk, const int blkjc0, const int blkjc1,
		const int *blkstart, const int *psdNL, const int *cumpsdNL,
		const int rpsdN, double *fwork)
{
  int inz, knz, k, nk, lastzsub;

  inz = 0;                 /* points into zir */
  if(inz >= zjc1)
    return inz;            /* return if no target subscripts */
  lastzsub = zir[zjc1-1];
  u -= blkstart[0];        /* Now, u + blkstart[k] gives Uk */
  for(knz = blkjc0; knz < blkjc1; knz++){
/* ------------------------------------------------------------
   For each nonzero block xblk[k] in x, Let zk = vec(tril(Dk * Xk * Dk)).
   ------------------------------------------------------------ */
    k = xblk[knz];
    if(k >= rpsdN)
      break;                                /* 1st nz Hermitian PSD block */
    if(lastzsub < blkstart[k])
      break;                                /* already beyond last target i*/
    nk = psdNL[k];
    while(zir[inz] < blkstart[k])           /* move to block k, viz. */
      z[zir[inz++]] = 0.0;                  /* all-0 where x is all-0 */
    inz = sprealutxu(z, zir, inz, zjc1, u + blkstart[k], invperm+cumpsdNL[k],
		     x, blkstart[k], nk, fwork);
    x += SQR(nk);                           /* point to next nonzero block */
  }
/* ------------------------------------------------------------
   Hermitian PSD blocks
   ------------------------------------------------------------ */
  for(; knz < blkjc1; knz++){
    k = xblk[knz];
    if(lastzsub < blkstart[k])
      break;                                /* already beyond last target i*/
    nk = psdNL[k];
    while(zir[inz] < blkstart[k])           /* move to block k, viz. */
      z[zir[inz++]] = 0.0;                  /* all-0 where x is all-0 */
    inz = spcpxutxu(z, zir, inz, zjc1, u + blkstart[k], invperm+cumpsdNL[k],
                    x, blkstart[k], nk, fwork);
    x += 2*SQR(nk);                         /* point to next nonzero block */
  }
/* ------------------------------------------------------------
   No remaining x blocks ==> remaining z-entries all-0.
   ------------------------------------------------------------ */
  for(; inz < zjc1; inz++)
    z[zir[inz]] = 0.0;                  /* all-0 where x is all-0 */
/* ------------------------------------------------------------
   RETURN next zir-position (should be zjc1)
   ------------------------------------------------------------ */
  return inz;       /* arrived at end of z-vector */
}
#endif