git.oblomov.eu Git - ohcount/blob - test/src_dir/octave1.m

   1 ## Copyright (C) 2006, Regents of the University of California  -*- mode: octave; -*-
   2 ##
   3 ## This program is free software distributed under the "modified" or
   4 ## 3-clause BSD license appended to this file.
   5
   6 function varargout = toledolu(LU)
   7   ## -*- texinfo -*-
   8   ## @deftypefn{Function File} {[@var{L}, @var{U}, @var{P}]} = toledolu(@var{A})
   9   ## @deftypefnx{Function File} {[@var{L}, @var{U}]} = toledolu(@var{A})
  10   ## @deftypefnx{Function File} {@var{LUP}} = toledolu(@var{A})
  11   ##
  12   ## Note: This returns a vector of indices for @var{P} and not a permutation
  13   ## matrix.
  14   ##
  15   ## Factors @var{P}*@var{A}=@var{L}*@var{U} by Sivan Toledo's recursive factorization algorithm
  16   ## with partial pivoting.  While not as fast as the built-in LU, this
  17   ## is significantly faster than the standard, unblocked algorithm
  18   ## while remaining relatively easy to modify.
  19   ##
  20   ## See the help for lu for details about the other calling forms.
  21   ##
  22   ## For the algorithm, see
  23   ## @itemize
  24   ## @item
  25   ## Toledo, Sivan. "Locality of reference in LU decomposition with
  26   ## partial pivoting," SIAM J. of Matrix Analysis and Applications,
  27   ## v18, n4, 1997. DOI: 10.1137/S0895479896297744
  28   ## @end itemize
  29   ##
  30   ## @seealso{lu}
  31   ##
  32   ## @end deftypefn
  33
  34   ## Author: Jason Riedy <ejr@cs.berkeley.edu>
  35   ## Keywords: linear-algebra, LU, factorization
  36   ## Version: 0.2
  37
  38   ## This version isn't *quite* the same as Toledo's algorithm.  I use a
  39   ## doubling approach rather than using recursion.  So non-power-of-2
  40   ## columns likely will be slightly different, but that shouldn't
  41   ## affect the 'optimality' by more than a small constant factor.
  42
  43   ## Also, I don't handle ncol > nrow optimally.  The code factors the
  44   ## first nrow columns and then updates the remaining ncol-nrow columns
  45   ## with L.
  46
  47   ## Might be worth eating the memory cost and tracking L separately.
  48   ## The eye(n)+tril(LU,-1) could be expensive.
  49
  50   switch (nargout)
  51     case 0
  52       return;
  53     case {1,2,3}
  54     otherwise
  55       usage ("[L,U,P] = lu(A), [P\\L, U] = lu(A), or (P\\L-I+U) = lu(A)");
  56   endswitch
  57
  58   [nrow, ncol] = size(LU);
  59   nstep = min(nrow, ncol);
  60
  61   Pswap = zeros(nstep, 1);
  62
  63   for j=1:nstep,
  64     [pval, pind] = max(abs(LU(j:nrow, j)));
  65     pind = pind + j - 1;
  66     Pswap(j) = pind;
  67
  68     kahead = bitand(j, 1+bitcmp(j)); # last 1 bit in j
  69     kstart = j+1-kahead;
  70     kcols = min(kahead, nstep-j);
  71
  72     inds = kstart : j;
  73     Ucol = j+1 : j+kcols;
  74     Lrow = j+1 : nrow;
  75
  76     ## permute just this column
  77     if (pind != j)
  78       tmp = LU(pind, j);
  79       LU(pind, j) = LU(j,j);
  80       LU(j,j) = tmp;
  81     endif
  82     ## apply pending permutations to L
  83     n_to_piv = 1;
  84     ipivstart = j;
  85     jpivstart = j - n_to_piv;
  86     while (n_to_piv < kahead)
  87       pivcols = jpivstart : jpivstart+n_to_piv-1;
  88       for ipiv = ipivstart:j,
  89         pind = Pswap(ipiv);
  90         if (pind != ipiv)
  91           tmp = LU(pind, pivcols);
  92           LU(pind, pivcols) = LU(ipiv, pivcols);
  93           LU(ipiv, pivcols) = tmp;
  94         endif
  95       endfor
  96       ipivstart -= n_to_piv;
  97       n_to_piv *= 2;
  98       jpivstart -= n_to_piv;
  99     endwhile
 100
 101     if (LU(j,j) != 0.0 && !isnan(LU(j,j))),
 102       LU(j+1:nrow,j) /= LU(j,j);
 103     endif
 104
 105     if 0 == kcols, break; endif
 106
 107     ## permute U to match perm already applied to L
 108     for k = inds,
 109       tmp = LU(Pswap(k), Ucol);
 110       LU(Pswap(k), Ucol) = LU(k, Ucol);
 111       LU(k, Ucol) = tmp;
 112     endfor
 113
 114     LU(inds, Ucol) = (eye(kahead) + tril(LU(inds, inds),-1)) \ LU(inds, Ucol);
 115     LU(Lrow, Ucol) -= LU(Lrow, inds) * LU(inds, Ucol);
 116   endfor
 117
 118   ## handle pivot permutations in L out from the last step
 119   npived = bitand(nstep, 1+bitcmp(nstep));
 120   j = nstep-npived;
 121   while (j > 0)
 122     n_to_piv = bitand(j, 1+bitcmp(j));
 123
 124     pivcols = j-n_to_piv+1 : j;
 125     for ipiv = j+1:nstep,
 126       pind = Pswap(ipiv);
 127       if (pind != ipiv)
 128         tmp = LU(pind, pivcols);
 129         LU(pind, pivcols) = LU(ipiv, pivcols);
 130         LU(ipiv, pivcols) = tmp;
 131       endif
 132     endfor
 133
 134     j -= n_to_piv;
 135   endwhile
 136
 137   if (nrow < ncol),
 138     Ucol = nrow+1 : ncol;
 139     inds = 1:nrow;
 140     for k = inds,
 141       tmp = LU(Pswap(k), Ucol);
 142       LU(Pswap(k), Ucol) = LU(k, Ucol);
 143       LU(k, Ucol) = tmp;
 144     endfor
 145     LU(inds, Ucol) = (eye(nrow) + tril(LU(inds, inds),-1)) \ LU(inds, Ucol);
 146   endif
 147
 148   if (nargout == 1)
 149     varargout{1} = LU;
 150     return;
 151   endif
 152
 153   if nrow == ncol,
 154     L = eye(nrow) + tril(LU, -1);
 155     varargout{2} = triu(LU);
 156   elseif nrow < ncol,
 157     L = eye(nrow) + tril(LU, -1)(:,1:nrow);
 158     varargout{2} = triu(LU);
 159   else # nrow > ncol
 160     L = tril(LU, -1);
 161     for k=1:ncol,
 162       L(k,k) = 1;
 163     endfor
 164     varargout{2} = triu(LU)(1:ncol,:);
 165   endif
 166
 167   if (nargout == 2)
 168     for j = 1:nstep,
 169       pind = Pswap(j);
 170       tmp = L(pind,:);
 171       L(pind,:) = L(j,:);
 172       L(j,:) = tmp;
 173     endfor
 174   else # nargout == 3
 175     P = 1:nrow;
 176     for j = 1:nstep,
 177       tmp = P(j);
 178       P(j) = P(Pswap(j));
 179       P(Pswap(j)) = tmp;
 180     endfor
 181     varargout{3} = P;
 182   endif
 183   varargout{1} = L;
 184
 185 endfunction
 186
 187 %!test
 188 %!  M = 15;
 189 %!  N = 15;
 190 %!  A = rand(M,N);
 191 %!  [L,U,P] = toledolu(A);
 192 %!  assert(norm(L*U-A(P,:),inf), 0, M**2*N*eps)
 193
 194 %!test
 195 %!  M = 16;
 196 %!  N = 16;
 197 %!  A = rand(M,N);
 198 %!  [L,U,P] = toledolu(A);
 199 %!  assert(norm(L*U-A(P,:),inf), 0, M**2*N*eps)
 200
 201 %!test
 202 %!  M = 17;
 203 %!  N = 17;
 204 %!  A = rand(M,N);
 205 %!  [L,U,P] = toledolu(A);
 206 %!  assert(norm(L*U-A(P,:),inf), 0, M**2*N*eps)
 207
 208 %!test
 209 %!  M = 8;
 210 %!  N = 17;
 211 %!  A = rand(M,N);
 212 %!  [L,U,P] = toledolu(A);
 213 %!  assert(norm(L*U-A(P,:),inf), 0, M**2*N*eps)
 214
 215 %!test
 216 %!  M = 8;
 217 %!  N = 15;
 218 %!  A = rand(M,N);
 219 %!  [L,U,P] = toledolu(A);
 220 %!  assert(norm(L*U-A(P,:),inf), 0, M**2*N*eps)
 221
 222 %!test
 223 %!  M = 7;
 224 %!  N = 17;
 225 %!  A = rand(M,N);
 226 %!  [L,U,P] = toledolu(A);
 227 %!  assert(norm(L*U-A(P,:),inf), 0, M**2*N*eps)
 228
 229 %!test
 230 %!  M = 7;
 231 %!  N = 15;
 232 %!  A = rand(M,N);
 233 %!  [L,U,P] = toledolu(A);
 234 %!  assert(norm(L*U-A(P,:),inf), 0, M**2*N*eps)
 235
 236 %!test
 237 %!  M = 17;
 238 %!  N = 8;
 239 %!  A = rand(M,N);
 240 %!  [L,U,P] = toledolu(A);
 241 %!  assert(norm(L*U-A(P,:),inf), 0, M**2*N*eps)
 242
 243 %!test
 244 %!  M = 15;
 245 %!  N = 8;
 246 %!  A = rand(M,N);
 247 %!  [L,U,P] = toledolu(A);
 248 %!  assert(norm(L*U-A(P,:),inf), 0, M**2*N*eps)
 249
 250 %!test
 251 %!  M = 17;
 252 %!  N = 7;
 253 %!  A = rand(M,N);
 254 %!  [L,U,P] = toledolu(A);
 255 %!  assert(norm(L*U-A(P,:),inf), 0, M**2*N*eps)
 256
 257 %!test
 258 %!  M = 15;
 259 %!  N = 7;
 260 %!  A = rand(M,N);
 261 %!  [L,U,P] = toledolu(A);
 262 %!  assert(norm(L*U-A(P,:),inf), 0, M**2*N*eps)
 263
 264 %!test
 265 %!  M = 31;
 266 %!  N = 17;
 267 %!  A = rand(M,N);
 268 %!  [L,U,P] = toledolu(A);
 269 %!  assert(norm(L*U-A(P,:),inf), 0, M**2*N*eps)
 270
 271 %!test
 272 %!  M = 11;
 273 %!  N = 29;
 274 %!  A = rand(M,N);
 275 %!  [L,U,P] = toledolu(A);
 276 %!  assert(norm(L*U-A(P,:),inf), 0, M**2*N*eps)
 277
 278 ## Copyright (c) 2006, Regents of the University of California
 279 ## All rights reserved.
 280 ## Redistribution and use in source and binary forms, with or without
 281 ## modification, are permitted provided that the following conditions are met:
 282 ##
 283 ##     * Redistributions of source code must retain the above copyright
 284 ##       notice, this list of conditions and the following disclaimer.
 285 ##     * Redistributions in binary form must reproduce the above copyright
 286 ##       notice, this list of conditions and the following disclaimer in the
 287 ##       documentation and/or other materials provided with the distribution.
 288 ##     * Neither the name of the University of California, Berkeley nor the
 289 ##       names of its contributors may be used to endorse or promote products
 290 ##       derived from this software without specific prior written permission.
 291 ##
 292 ## THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND ANY
 293 ## EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 294 ## WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 295 ## DISCLAIMED. IN NO EVENT SHALL THE REGENTS AND CONTRIBUTORS BE LIABLE FOR ANY
 296 ## DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 297 ## (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 298 ## LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 299 ## ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 300 ## (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 301 ## SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.