blkchol2.c

经济学专业代码

/*

% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA


*/
#include <math.h>
#include <string.h>
#include "blksdp.h"
#include "mex.h"

/* ------------------------------------------------------------
   PROTOTYPES:
   ------------------------------------------------------------ */
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);

/* ************************************************************
   TIME-CRITICAL PROCEDURE -- isscalarmul(x,alpha,n)
   Computes x *= alpha using LEVEL 8 loop-unrolling.
   ************************************************************ */
void isscalarmul(double *x, const double alpha, const int n)
{
  int i;
  
  for(i=0; i< n-7; ){          /* LEVEL 8 */
    x[i++] *= alpha;
    x[i++] *= alpha;
    x[i++] *= alpha;
    x[i++] *= alpha;
    x[i++] *= alpha;
    x[i++] *= alpha;
    x[i++] *= alpha;
    x[i++] *= alpha;
  }
/* ------------------------------------------------------------
   Now, i in {n-7, n-6, ..., n}. Do the last n-i elements.
   ------------------------------------------------------------ */
  if(i < n-3){                           /* LEVEL 4 */
    x[i++] *= alpha;
    x[i++] *= alpha;
    x[i++] *= alpha;
    x[i++] *= alpha;
  }
  if(i < n-1){                           /* LEVEL 2 */
    x[i++] *= alpha;
    x[i++] *= alpha;
  }
  if(i < n)                              /* LEVEL 1 */
    x[i] *= alpha;
}

/* ************************************************************
   PROCEDURE maxabs - computes inf-norm
   INPUT
     x - vector of length n
     n - length of x.
   RETURNS y = norm(x,inf).
   ************************************************************ */
double maxabs(const double *x,const int n)
{
  int i;
  double y,xi;
  y = 0.0;
  for(i = 0; i < n; i++)
    if((xi = fabs(x[i])) > y)
      y = xi;
  return y;
}

/* ************************************************************
   PROCEDURE cholonBlk - CHOLESKY on a dense diagonal block.
            Also updates nonzeros below this diagonal block -
            they need merely be divided by the scalar diagonals
            "lkk" afterwards.
   INPUT
     m      - number of rows (length of the first column).
     ncols  - number of columns in the supernode.(n <= m)
     lb     - Length ncols. Skip k-th pivot if drops below lb[k].
     ub     - max(diag(x)) / maxu^2. No stability check for pivots > ub.
     maxu   - Max. acceptable |lik|/lkk when lkk suffers cancelation.
     first - global column number of column 0. This is used only to insert
        the global column numbers into skipIr.
   UPDATED
     x  - On input, contains the columns of the supernode to
          be factored. On output, contains the factored columns of
          the supernode.
     skipIr - Lists skipped pivots with their global column number
        in 0:neqns-1. Active range is first:first+ncols-1.
        Skipped if d(k) suffers cancelation and max(abs(L(:,k)) > maxu.
     *pnskip - nnz in skip; *pnskip <= order N of sparse matrix.
   OUTPUT
     d - Length ncols. Diagonal in L*diag(d)*L' with diag(L)=all-1.
   ************************************************************ */
void cholonBlk(double *x, double *d, int m, const int ncols, const int first,
               const double ub, const double maxu, double *lb,
               int *skipIr, int *pnskip)
{
  int inz,i,k,n,coltail, nskip;
  double xkk, xik, ubk;
  double *xi;
/* ------------------------------------------------------------
   Initialize:
   ------------------------------------------------------------ */
  n = ncols;
  nskip = *pnskip;
  inz = 0;
  coltail = m - ncols;
  for(k = 0; k < ncols; k++, --m, --n){
/* -------------------------------------------------------
   Let xkk = L(k,k), ubk = max(|xik|) / maxu.
   ------------------------------------------------------- */
    xkk = x[inz];
    if(xkk > lb[k]){ /* now xkk > 0 */
      if(xkk < ub){
        ubk = maxabs(x+inz+1,m-1) / maxu;
        if(xkk < ubk){
/* ------------------------------------------------------------
   If we need to add on diagonal, store this in (skipIr, lb(k)).
   ------------------------------------------------------------ */
          skipIr[nskip++] = first + k;
	  lb[k] = ubk - xkk;           /* amount added on diagonal */
	  xkk = ubk;
	}
      }
/* --------------------------------------------------------------
   Set dk = xkk, lkk = 1 (for LDL').
   -------------------------------------------------------------- */
      d[k] = xkk;                   /* now d[k] > 0 MEANS NO-SKIPPING */
      x[inz] = 1.0;
      xi = x + inz + m;                 /* point to next column */
      ++inz;
/* --------------------------------------------------------------
   REGULAR JOB: correct remaining n-k cols with col k.
   x(k+1:m,k+1:n) -= x(k+1:m,k) * x(k+1:n,k)' / xkk
   x(k+1:n,k) /= xkk,
   -------------------------------------------------------------- */
      for(i = 1; i < n; i++){
        xik = x[inz] / xkk;
        subscalarmul(xi, xik, x+inz, m-i);
        x[inz++] = xik;
        xi += m-i;
      }
      inz += coltail;                 /* Let inz point to next column */
    }
/* ------------------------------------------------------------
   If skipping is enabled and this pivot is too small:
   1) don't touch L(k:end,k): allows pivot delaying if desired.
   2) List first+k in skipIr. Set dk = 0 (MEANS SKIPPING).
   -------------------------------------------------------------- */
    else{
      skipIr[nskip++] = first + k;
      d[k] = 0.0;                 /* tag "0": means column skipped in LDL'.*/
      inz += m;                   /* Don't touch nor use L(k:end,k) */
    }
  } /* k=0:ncols-1 */
/* ------------------------------------------------------------
   Return updated number of added or skipped pivots.
   ------------------------------------------------------------ */
  *pnskip = nskip;
}

/* ************************************************************
  getbwIrInv  --  Inverse of the subscript function: given a subscript,
      irInv yields the position, counted FROM THE BOTTOM of the sparse column.

  INPUT PARAMETERS -
     nnz    - LENGTH OF THE FIRST COLUMN OF THE SUPERNODE,
              INCLUDING THE DIAGONAL ENTRY.
     Lir    - Lir[0:nnz-1] ARE THE ROW INDICES OF THE NONZEROS
              OF THE FIRST COLUMN OF THE SUPERNODE.
  OUTPUT PARAMETERS - 
     irInv - On return, irInv[Lir[0:nnz-1]] = nnz:-1:1, so that
		           Lir[nnz-irInv[i]]  == i
             The position of subscript "xij" is thus
			   xjc[j+1] - irInv[i].
   ************************************************************ */
void getbwIrInv(int *irInv, const int *Lir, const int nnz)
{
  int inz,bwinz;

  bwinz = nnz;
  for(inz = 0; inz < nnz; inz++, bwinz--)
    irInv[Lir[inz]] = bwinz;               /* bwinz = nnz:-1:1 */
}

/* ************************************************************
  suboutprod  --  Computes update from a single previous column "xk" on
		a supernode "xj", using dense computations.
  INPUT
     mj, nj  -  supernode "xj" is mj x nj.  More precisely, the column
                lengths are {mj, mj-1, ..., mj-(nj-1)}.
     xkk     -  scalar, the 1st nj entries in xk are divided by this number.
     mk      -  length of xk.  WE ASSUME mk <= mj.  Only 1st mk rows in xj
                are updated.
  UPDATED
     xj  -  On return, xj -= xk*xk(0:nj-1)'/xkk
     xk  -  On return, xk(0:nj-1) /= xkk
   ************************************************************ */
void suboutprod(double *xj, int mj, const int nj, double *xk,
                const double xkk, int mk)
{
  int j;
  double xjk;

  for(j = 0; j < nj; j++){
    xjk = xk[0] / xkk;
    subscalarmul(xj, xjk, xk, mk);   /* xj -= xjk * xk */
    xk[0] = xjk;                     /* FINAL entry ljk */
    xj += mj;                    /* point to next column which is 1 shorter */
    --mj; --mk; ++xk;
  }
}

/* ************************************************************
  isminoutprod  --  Computes update from a column "xk" and stores it in "xj",
	       using dense computations. If "xkk<=0", then let xj = 0.
  INPUT
     mk, nj  -  output "xj" is mk x nj - nj*(nj-1)/2. Its column lengths are
	        {mk, mk-1, ..., mk-(nj-1)}.
     xkk     -  scalar, the 1st nj entries in xk are divided by this number.
  OUTPUT
     xj      -  On return, xj = -xk*xk(0:nj-1)'/xkk       (NOTE THE MINUS !)
                BUT: if xkk <= 0, then xj = zeros(nj*(2m-nj+1)/2,1).
  UPDATED
     xk      -  On return, xk(0:nj-1) /= xkk if xkk > 0, otherwise unchanged.
   ************************************************************ */
void isminoutprod(double *xj, const int nj, double *xk, const double xkk,
                  int mk)
{
  int j;
  double xjk;

  if(xkk > 0.0)   /* if not phase 2 node */
    for(j = 0; j < nj; j++){
      xjk = xk[0] / xkk;
      memcpy(xj,xk,mk * sizeof(double));
      isscalarmul(xj, -xjk, mk);          /* xj = -xjk * xk */
      xk[0] = xjk;                     /* FINAL entry ljk */
      xj += mk;                /* point to next column which is 1 shorter */
      --mk; ++xk;
    }
  else  /* initialize to all-0 if phase-2 node */
    fzeros(xj,(nj * (mk + mk-nj + 1))/2);
}

/* ************************************************************
   spsuboutprod  --  Computes update from a single previous column "xk" on
		     a supernode "xj", with a different sparsity structure.
                     The relevant nonzeros of xj are accessed by a single
                     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 nonzeros
                in xk.
     xjjc    - xjjc[0] is start of 1st column of xj (as index into xnz), etc.
     xkk     -  scalar, the 1st nj entries in xk are divided by this number.
     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.
  UPDATED
     xnz  -  On return, xj(relind,:) -= xk*xk(0:nj-1)'/xkk
     xk   -  On return, xk(0:nj-1) /= xkk
   ************************************************************ */