Page MenuHomec4science
Files F121995616

bwblkslv2.c
No OneTemporary
Actions

File Metadata

Created
Tue, Jul 15, 06:05

bwblkslv2.c
View Options

/*
% 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 "mex.h"
#include "blkchol.h"
/* ============================================================
BACKWARD SOLVE:
============================================================ */
/* ************************************************************
PROCEDURE bwsolve -- Solve y from L'*y = b, where
L is lower-triangular.
INPUT
Ljc, Lir, Lpr - sparse lower triangular matrix
xsuper - starting column in L for each (dense) supernode.
nsuper - number of super nodes
UPDATED
y - full xsuper[nsuper]-vector, yOUTPUT = L' \ yINPUT.
************************************************************ */
void bwsolve(double *y, const mwIndex *Ljc, const mwIndex *Lir, const double *Lpr,
const mwIndex *xsuper, const mwIndex nsuper)
{
mwIndex jsup,j,inz,k;
double yj, ljj;
/* ------------------------------------------------------------
For each supernode jsup:
------------------------------------------------------------ */
j = xsuper[nsuper]; /* column after current snode (j=m)*/
for(jsup = nsuper; jsup > 0; jsup--){
for(k = 0; k < xsuper[jsup] - xsuper[jsup-1]; k++){
/* ------------------------------------------------------------
The equation L(:,j)'*yNEW = yOLD(j), yields
y(j) -= L(j+1:m,j)'*y.
------------------------------------------------------------ */
inz = Ljc[j-1];
inz++; /* jump over diagonal entry */
yj = realdot(Lpr+inz, y+j, k);
for(inz += k; inz < Ljc[j]; inz++)
yj += Lpr[inz] * y[Lir[inz]];
y[--j] -= yj;
}
}
}
/* ************************************************************
PROCEDURE partbwsolve -- Solve y from L(0:m-1,0:m-1)'*y = b, where
L is lower-triangular.
INPUT
L=(ljc,lpr,xlindx,lindx) - sparse lower triangular matrix
xsuper - starting column in L for each (dense) supernode.
nsuper - number of super nodes (for NW-subblock)
m - order of L'- (sub) block, xsuper[nsuper-1] < m <= xsuper[nsuper]
To solve with the complete L, choose m = xsuper[nsuper].
UPDATED
y - full m-vector, yOUT = L(0:m-1,0:m-1)'\yIN.
************************************************************ */
void partbwsolve(double *y, const mwIndex *ljc, const double *lpr,
const mwIndex *xlindx, const mwIndex *lindx, const mwIndex *xsuper,
const mwIndex nsuper, const mwIndex m)
{
mwIndex jsup,j,inz,k,i, ixfirst,ixnz;
double yj;
/* ------------------------------------------------------------
For each supernode jsup:
Let ixfirst point to the 1st row-subscript below current supernode.
------------------------------------------------------------ */
j = m; /* column/row after current entry */
for(jsup = nsuper; jsup > 0; jsup--){
ixfirst = xlindx[jsup-1] + (xsuper[jsup] - xsuper[jsup-1]);
/* ------------------------------------------------------------
Case 1: L has sparse nonzeros below current supernode, but not
beyond m:
------------------------------------------------------------ */
if(ixfirst < xlindx[jsup])
if(lindx[xlindx[jsup] - 1] < m)
for(k = 0; j > xsuper[jsup-1]; k++){
/* ------------------------------------------------------------
The equation L(:,j)'*y=b(j), yields
y(j) = b(j)-L(j+1:m,j)'*y.
------------------------------------------------------------ */
inz = ljc[j-1] + 1; /* jump over diagonal entry */
yj = realdot(lpr+inz, y+j, k);
for(inz += k, ixnz = ixfirst; inz < ljc[j]; inz++, ixnz++)
yj += lpr[inz] * y[lindx[ixnz]];
y[--j] -= yj;
}
else
/* ------------------------------------------------------------
2nd CASE: L(:,j) has nonzeros beyond the requested order m,
so use while i < m.
------------------------------------------------------------ */
for(k = 0; j > xsuper[jsup-1]; k++){
inz = ljc[j-1] + 1; /* jump over diagonal entry */
yj = realdot(lpr+inz, y+j, k);
inz += k; ixnz = ixfirst;
for(i=lindx[ixnz]; i < m; i = lindx[++ixnz])
yj += lpr[inz++] * y[i];
y[--j] -= yj;
}
else
/* ------------------------------------------------------------
3rd CASE: L(:,j) has only nonzeros in snode jsup: purely dense.
------------------------------------------------------------ */
for(k = 1, --j; j > xsuper[jsup-1]; k++){
yj = realdot(lpr + ljc[j-1] + 1, y+j, k);
y[--j] -= yj;
}
}
}

Event Timeline

c4science · Help