blkchol2.c

经济学专业代码

void spsuboutprod(const int *xjjc, double *xnz, const int mj, const int nj,
                  double *xk,const double xkk,const int mk, const int *relind)
{
  int i, j, jcol, bottomj;
  double xjk;

  ++xjjc;             /* now it points beyond bottom of columns */
  for(j = 0; j < nj; j++){
    jcol = mj - relind[j];       /* affected column */
    bottomj = xjjc[jcol];
    xjk = xk[j] / xkk;
    for(i = j; i < mk; i++)
      xnz[bottomj - relind[i]] -= xjk * xk[i];
    xk[j] = xjk;                     /* FINAL entry ljk */
  }
}

/* ************************************************************
  spadd  --  Let xj += xk, where the supernode "xj", has a sparsity
        structure. The relevant nonzeros of xj are accessed by a indirection,
	via "relind[:]".
  INPUT
     mj, nj  -  supernode "xj" has mj rows in its 1st column. In total, we
	        will update nj columns, corresponding to the 1st nj nonzero
                rows in xk.
     xjjc    - xjjc[0] is start of 1st column of xj (as index into xnz), etc.
     mk      -  length of xk.  WE ASSUME mk <= mj.
     relind  - (mj - relind[0:mk-1]) yields the locations in xj on which the
	       xk nonzeros will act.
     xk      -  mk * nk - nk*(nk-1)/2 matrix, with column lengths
	        mk, mk-1, mk-2,.. mk-(nj-1).  They'll be substracted from
                the entries in xj that are listed by relind.
  UPDATED
     xnz     -  On return, xj(relind,:) += xk
   ************************************************************ */
void spadd(const int *xjjc, double *xnz, const int mj, const int nj,
           const double *xk, const int mk, const int *relind)
{
  int i, j, jcol, bottomj,mkcol;

  ++xjjc;             /* now it points beyond bottom of columns */
  mkcol = mk;         /* mkcol = mk - j */
  for(j = 0; j < nj; j++){
    jcol = mj - relind[j];       /* affected column */
    bottomj = xjjc[jcol];
    for(i = j; i < mk; i++)
      xnz[bottomj - relind[i]] += xk[i];
    xk += (--mkcol);   /* xk(i:mk-1) is next column */
  }
}

/* ************************************************************      
   PROCEDURE precorrect  -  Apply corrections from affecting supernode
      (skipping subnodes with non-positive diagonal) on supernodal
      diagonal block in L-factor.
   INPUT
     ljc   - start of columns in lpr.
     d     - Length neqns vector. The diagonal in L*diag(d)*L'. Only
       d[firstk:nextk-1] will be used.
     irInv - For row-indices Jir of affected supernode, Jir[m-irInv[i]]  == i.
     nextj - Last subnode of affected supernode is nextj-1.
     firstk, nextk - subnodes of affecting supernode are firstk:nextk-1.
     Kir   - unfinished row indices of affecting supernode
     mk    - number of unfinished nonzeros in affecting supernode
     fwsiz - Allocated length of fwork.
   UPDATED
     lpr  - For each column k=firstk:nextk-1, and the affected columns j
       in node, DO  L(:,j) -= (ljk / lkk) * L(:,k),
       and store the definitive j-th row of L, viz. ljk /= lkk.
   WORKING ARRAYS
     relind - length mk integer array
     fwork  - length fwsiz vector, for storing -Xk * inv(LABK) * Xk'.
   RETURNS  ncolup, number of columns updated by snode k.
    if -1, then fwsiz is too small.
   ************************************************************ */
int precorrect(double *lpr, const int *ljc,const double *d, const int *irInv,
               const int nextj, const int *Kir, const int mk,
               const int firstk, const int nextk,
               int *relind, const int fwsiz, double *fwork)
{
  int i,j,k,ncolup,mj;
  double *xj;
/* ------------------------------------------------------------
   j = first subscript in k (== 1st affected column)
   i = last subscript in k
   ncolup = number of nz-rows in k corresponding to columns in node.
   mj = number of nonzeros in l(:,j), the 1st affected column
   ------------------------------------------------------------ */
  j = Kir[0];
  i = Kir[mk-1];
  if(i < nextj)
    ncolup = mk;
  else
    for(ncolup = 1; Kir[ncolup] < nextj; ncolup++);
  mj = ljc[j+1] - ljc[j];
/* ------------------------------------------------------------
   If nz-structure of k is a single block in structure of node,
   (i.e. irInv[Kir[0]] - irInv[Kir[mk-1]] == mk-1). The subnodes
   of "node" must then be consecutive and at the start.
   Thus, we use dense computations :
   ------------------------------------------------------------ */
  if(irInv[j] - irInv[i] < mk){
    xj = lpr + ljc[j];
    for(k = firstk; k < nextk; k++)
      if(d[k] > 0.0)                        /* Skip pivot when d[k] <= 0 */
        suboutprod(xj, mj, ncolup, lpr + ljc[k+1] - mk, d[k], mk);
  }
  else{
/* ------------------------------------------------------------
   Otherwise, the nz-indices of k are scattered within the structure of node.
   Let relind be the position of these nz's in node, COUNTED FROM THE BOTTOM.
   ------------------------------------------------------------*/
    for(i = 0; i < mk; i++)
      relind[i] = irInv[Kir[i]];
/* ------------------------------------------------------------
   If k is a single column, then perform update directly in lpr:
   ------------------------------------------------------------ */
    if(nextk - firstk == 1){
      if(d[firstk] > 0.0)                   /* Skip pivot when d[k] <= 0 */
        spsuboutprod(ljc+j,lpr,mj, ncolup, lpr + ljc[nextk]-mk,
                     d[firstk],mk, relind);
    }
    else{
/* ------------------------------------------------------------
   Multiple columns in affecting snode:
   1. compute the complete modification, and store it in fwork:
   fwork = -Xk * inv(LABK) * Xk'
   ------------------------------------------------------------ */
      if(fwsiz < mk * ncolup - ncolup*(ncolup-1)/2)
        return -1;
      for(k = firstk; k < nextk; k++)      /* find 1st positive diag */
        if(d[k] > 0.0)
          break;
      if(k < nextk){                       /* if any positive diag: */
        isminoutprod(fwork, ncolup, lpr + ljc[k+1] - mk, d[k], mk);
        for(++k; k < nextk; k++)           /* remaining cols */
          if(d[k] > 0.0)                   /* Skip pivot when d[k] <= 0 */
            suboutprod(fwork, mk, ncolup, lpr + ljc[k+1] - mk, d[k], mk);
/* ------------------------------------------------------------
   2. subtract fwork from the sparse columns of node, using relind.
   ------------------------------------------------------------ */
        spadd(ljc+j,lpr,mj, ncolup, fwork,mk, relind);
      } /* end exists positive diag */
    } /* end multiple affecting cols */
  } /* end of scattered case */
/* ------------------------------------------------------------
   RETURN number of columns updated, i.e. #subnodes in k that we finished.
   ------------------------------------------------------------ */
  return ncolup;
}

/* ************************************************************
   BLKLDL  --  Block-sparse L*D*L' Cholesky factorization.

   INPUT:
      neqns   - Order "m": L is neqns * neqns
      nsuper  - Number of supernodes (blocks).
      xsuper  - Length nsuper+1: first simple-node of each supernode
      snode   - Length neqns: snode(node) is the supernode containing "node".
      xlindx  - Length nsuper+1: Start of sparsity structure in lindx,
              for each supernode (all simple nodes in a supernode have the
              same nonzero-structure).
      lindx   - row indices, for each supernode.
      ljc     - Length neqns+1: start of the columns of L.
      perm    - Length neqns: reordering of pne->At columns in Cholesky.
      ub     - max(diag(x)) / maxu^2. No stability check for pivots > ub.
      maxu    - Force max(max(abs(L))) <= maxu (by adding low-rank diag).
      iwsiz, fwsiz - size of integer and floating-point working storage.
               See "WORKING ARRAYS" for required amount.
   UPDATED:
      Lpr     - On input, contains tril(X), on output, L is
	      such that   X = L*D*L'. For columns k where d[k]=0, L(:,k)
              contains the column updated upto pivot k-1.
      lb    - Length neqns. INPUT: cancelation threshold per pivot. Skip pivot
          if it drops below.
          OUTPUT: lb(skipIr) are values of low rank diag. matrix that is
          added before factorization.
   OUTPUT
      d      - length neqns vector, diagonal in L*diag(d)*L'.
      skipIr - length nskip (<= neqns) array. skipIr(1:nskip) lists the
        columns that have been skipped in the Cholesky. d[skipIr] = 0.
   WORKING ARRAYS:
      iwork  - Length iwsiz working array, used for
           link(nsuper), length(nsuper),
           irInv(neqns), relind(neqns),
           iwsiz = 2*m + 2 * nsuper
      fwork  - Length fwsiz. Used for fwork(L.tmpsiz) in precorrect.
           fwsiz = L.tmpsiz.
   ACKNOWLEDGMENT:
       Parts are inspired by F77-block Cholesky of Ng and Peyton (ORNL).
   RETURNS  nskip (<=neqns), number of skipped nodes. Length of skipIr.
     if -1 then not enough workspace (iwsiz, fwsiz) allocated.
   ************************************************************ */
int blkLDL(const int neqns, const int nsuper, const int *xsuper,
           const int *snode,  const int *xlindx, const int *lindx,
           double *lb,
           const int *ljc, double *lpr, double *d, const int *perm,
           const double ub, const double maxu, int *skipIr,
           int iwsiz, int *iwork, int fwsiz, double *fwork)
{
  const int *Jir;
  int *link, *length, *irInv, *relind, *ncolupLst;
  int node,nextj,i,j,nnzj,n,inz,  k,colk,mk,linkk, ncolup,snodei, nskip;
/* ------------------------------------------------------------
   Partition integer working array of size 2*(nsuper+neqns):
   iwork = [link(nsuper); length(nsuper); irInv(neqns); relind(neqns)].
   ------------------------------------------------------------ */
  if(iwsiz < 2 * (neqns + nsuper))
    return -1;
  link   = iwork;                    /* 2 times length nsuper: */
  length = link + nsuper;
  irInv  = length + nsuper;          /* 2 * length neqns: */
  relind    = irInv + neqns;
/* ------------------------------------------------------------
   ncolupLst(neqns) shares the same working array as irInv(neqns).
   Namely, at stage j=xsuper[node], irInv uses only entries >= j,
   whereas ncolupLst only applies to the "old" columns < j.
   ------------------------------------------------------------ */
  ncolupLst = irInv;
/* ------------------------------------------------------------
   Initialize: link = nsuper * ones(nsuper,1) (means END-OF-LIST)
   ------------------------------------------------------------ */
  for(node = 0; node < nsuper; node++)
    link[node] = nsuper;
/* ------------------------------------------------------------
   Initialize nskip = 0.
   ------------------------------------------------------------ */
  nskip = 0;
/* ------------------------------------------------------------
   For each supernode "node", start at subnode j = xsuper[node],
   having sparsity pattern Jir.
   ------------------------------------------------------------ */
  nextj = xsuper[0];
  for(node = 0; node < nsuper; node++){
    j = nextj;		                /* 1st col in node */
    nextj = xsuper[node+1];
    n = nextj - j;			/* length of node */
    Jir = lindx + xlindx[node];         /* row-indices for column j */
    nnzj = ljc[j+1] - ljc[j];           /* nnz( column j ) */
/* ------------------------------------------------------------
   Compute inverse of Jir, yielding position from the bottom:
   Jir[nnzj-irInv[i]]  == i
   This will be handy when adding a column with a sub-sparsity structure
   to column j.
   ------------------------------------------------------------ */
    getbwIrInv(irInv, Jir, nnzj);
/* ------------------------------------------------------------
   Apply corrections from relevant previous super-nodes;
   these snodes are
   node -> link[node] -> link[link[node]] -> ...
   ------------------------------------------------------------ */
    for(k = link[node]; k < nsuper; k = link[k]){
      if((ncolupLst[k] = precorrect(lpr,ljc,d,irInv, nextj,
                                    lindx + xlindx[k+1]-length[k],
                                    length[k],xsuper[k],xsuper[k+1],
                                    relind,fwsiz,fwork)) < 0)
        return -1;         /* fwsiz too small */
    }	
/* ------------------------------------------------------------
   DO DENSE CHOLESKY on the current supernode
   ------------------------------------------------------------ */
    cholonBlk(lpr + ljc[j],d+j, nnzj, n, j, ub, maxu, lb+j, skipIr,&nskip);
/* ------------------------------------------------------------
   insert each current affecting snode k into linked list of
   next supernode it will affect.
   ------------------------------------------------------------ */
    for(k = link[node]; k < nsuper; k = linkk){
      linkk = link[k];
      mk = (length[k] -= ncolupLst[k]);    /* unfinished nonzeros in k */
      if(mk){                              /* if not yet terminated: */
        i = lindx[xlindx[k+1]-mk];
        snodei = snode[i];
        link[k] = link[snodei];            /* prev. also affecting i */
        link[snodei] = k;                  /* next snode it'll affect */
      }
    }
/* ------------------------------------------------------------
   The same for current snode "node" itself:
   ------------------------------------------------------------ */
    if((length[node] = nnzj - n) > 0){
      i = Jir[n];                    /* 1st row outside snode */
      snodei = snode[i];
      link[node] = link[snodei];     /* prev. also affecting i */
      link[snodei] = node;
    }
    else
      length[node] = 0;              /* Supernode terminated */
  } /* node = 0:nsuper-1 */
/* ------------------------------------------------------------
   FINISHING: return the number of skipped pivots
   ------------------------------------------------------------ */
  return nskip;
}