Mercurial > hg > octave-lojdl
changeset 8338:a35bf360b919
Add the cgs and treelayout functions
author | Radek Salac <salac.r@gmail.com> |
---|---|
date | Fri, 21 Nov 2008 15:03:03 +0100 |
parents | e02242c54c49 |
children | 18c4ded8612a |
files | doc/interpreter/contributors.in scripts/ChangeLog scripts/sparse/Makefile.in scripts/sparse/cgs.m scripts/sparse/treelayout.m |
diffstat | 5 files changed, 376 insertions(+), 2 deletions(-) [+] |
line wrap: on
line diff
--- a/doc/interpreter/contributors.in +++ b/doc/interpreter/contributors.in @@ -170,6 +170,7 @@ Olli Saarela Toni Saarela Juhani Saastamoinen +Radek Salac Ben Sapp Alois Schloegl Michel D. Schmid
--- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -1,3 +1,8 @@ +2008-11-21 Radek Salac <salac.r@gmail.com> + + * sparse/cgs.m, sparse/treelayout.m: New functions. + * sparse/Makefile.in (SOURCES): Add them here. + 2008-11-14 David Bateman <dbateman@free.fr> * plot/__go_draw_axes__.m (do_tics_1): Support the minorick properties
--- a/scripts/sparse/Makefile.in +++ b/scripts/sparse/Makefile.in @@ -32,10 +32,10 @@ INSTALL_PROGRAM = @INSTALL_PROGRAM@ INSTALL_DATA = @INSTALL_DATA@ -SOURCES = colperm.m etreeplot.m gplot.m nonzeros.m normest.m \ +SOURCES = cgs.m colperm.m etreeplot.m gplot.m nonzeros.m normest.m \ pcg.m pcr.m spalloc.m spaugment.m spconvert.m spdiags.m speye.m \ spfun.m sphcat.m spones.m sprand.m sprandn.m sprandsym.m spstats.m \ - spvcat.m spy.m treeplot.m + spvcat.m spy.m treelayout.m treeplot.m DISTFILES = $(addprefix $(srcdir)/, Makefile.in $(SOURCES))
new file mode 100644 --- /dev/null +++ b/scripts/sparse/cgs.m @@ -0,0 +1,160 @@ +## Copyright (C) 2008 Radek Salac +## +## This file is part of Octave. +## +## Octave 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 3 of the License, or (at +## your option) any later version. +## +## Octave 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 Octave; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {} cgs (@var{A}, @var{b}) +## @deftypefnx {Function File} {} cgs (@var{A}, @var{b}, @var{tol}, @var{maxit}, @var{M1}, @var{M2}, @var{x0}) +## This procedure attempts to solve a system of linear equations A*x = b for x. +## The @var{A} must be square, symmetric and positive definite real matrix N*N. +## The @var{b} must be a one column vector with a length of N. +## The @var{tol} specifies the tolerance of the method, default value is 1e-6. +## The @var{maxit} specifies the maximum number of iteration, default value is MIN(20,N). +## The @var{M1} specifies a preconditioner, can also be a function handler which returns M\X. +## The @var{M2} combined with @var{M1} defines preconditioner as preconditioner=M1*M2. +## The @var{x0} is initial guess, default value is zeros(N,1). +## +## @end deftypefn + +function [x, flag, relres, iter, resvec] = cgs (A, b, tol, maxit, M1, M2, x0) + + if (nargin < 2 || nargin > 7 || nargout > 5) + print_usage (); + elseif (!isnumeric (A) || rows (A) != columns (A)) + error ("cgs: first argument must be a n-by-n matrix"); + elseif (!isvector (b)) + error ("cgs: b must be a vector"); + elseif (rows (A) != rows (b)) + error ("cgs: first and second argument must have the same number of rows"); + elseif (nargin > 2 && !isscalar (tol)) + error ("cgs: tol must be a scalar"); + elseif (nargin > 3 && !isscalar (maxit)) + error ("cgs: maxit must be a scalar"); + elseif (nargin > 4 && ismatrix (M1) && (rows (M1) != rows (A) || columns (M1) != columns (A))) + error ("cgs: M1 must have the same number of rows and columns as A"); + elseif (nargin > 5 && (!ismatrix (M2) || rows (M2) != rows (A) || columns (M2) != columns (A))) + error ("cgs: M2 must have the same number of rows and columns as A"); + elseif (nargin > 6 && !isvector (x0)) + error ("cgs: x0 must be a vector"); + elseif (nargin > 6 && rows (x0) != rows (b)) + error ("cgs: xO must have the same number of rows as b"); + endif + + ## default toleration + if (nargin < 3) + tol = 1e-6; + endif + + ## default maximum number of iteration + if (nargin < 4) + maxit = min (rows (b),20); + endif + + + ## left preconditioner + precon = []; + if (nargin == 5) + precon = M1; + elseif (nargin > 5) + if (isparse(M1) && issparse(M2)) + precon = @(x) M1 * (M2 * x); + else + precon = M1 * M2; + endif + endif + + ## precon can by also function + if (nargin > 4 && isnumeric (precon)) + + ## we can compute inverse preconditioner and use quicker algorithm + if (det (precon) != 0) + precon=inv (precon); + else + error ("cgs: preconditioner is ill conditioned"); + endif + + ## we must make test if preconditioner isn't ill conditioned + if (isinf (cond (precon))); + error ("cgs: preconditioner is ill conditioned"); + endif + endif + + ## specifies initial estimate x0 + if (nargin < 7) + x = zeros (rows (b), 1); + else + x = x0; + endif + + relres = b - A * x; + ## vector of the residual norms for each iteration + resvec = [norm(relres)]; + ro = 0; + norm_b = norm (b); + ## default behaviour we don't reach tolerance tol within maxit iterations + flag = 1; + for iter = 1 : maxit + + if (nargin > 4 && isnumeric (precon)) + ## we have computed inverse matrix so we can use quick algorithm + z = precon * relres; + elseif (nargin > 4) + ## our preconditioner is a function + z = feval (precon, relres); + else + ## we don't use preconditioning + z = relres; + endif + + ## cache + ro_old = ro; + ro = relres' * z; + if (iter == 1) + p = z; + else + beta = ro / ro_old; + p = z + beta * p; + endif + q = A * p; #cache + alpha = ro / (p' * q); + x = x + alpha * p; + + relres = relres - alpha * q; + resvec = [resvec; norm(relres)]; + + relres_distance = resvec (end) / norm_b; + if (relres_distance <= tol) + ## we reach tolerance tol within maxit iterations + flag = 0; + break; + elseif (resvec (end) == resvec (end - 1) ) + ## the method stagnates + flag = 3; + break; + endif + endfor; + + relres = relres_distance; +endfunction + + + +%!demo +%! % Solve system of A*x=b +%! A=[5 -1 3;-1 2 -2;3 -2 3] +%! b=[7;-1;4] +%! [a,b,c,d,e]=cgs(A,b)
new file mode 100644 --- /dev/null +++ b/scripts/sparse/treelayout.m @@ -0,0 +1,208 @@ +## Copyright (C) 2008 Ivana Varekova & Radek Salac +## +## This file is part of Octave. +## +## Octave 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 3 of the License, or (at +## your option) any later version. +## +## Octave 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 Octave; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {} treelayout (@var{Tree}) +## @deftypefnx {Function File} {} treelayout (@var{Tree}, @var{Permutation}) +## treelayout lays out a tree or a forest. The first argument @var{Tree} is a vector of +## predecessors, optional parameter @var{Permutation} is an optional postorder permutation. +## The complexity of the algorithm is O(n) in +## terms of time and memory requirements. +## @seealso{etreeplot, gplot,treeplot} +## @end deftypefn + +function [XCoordinate, YCoordinate, Height, s] = treelayout (Tree, Permutation) + if (nargin < 1 || nargin > 2 || nargout > 4) + print_usage (); + elseif (! isvector (Tree) || rows (Tree) != 1 || ! isnumeric (Tree) + || any (Tree > length (Tree)) || any (Tree < 0) ) + error ("treelayout: the first input argument must be a vector of predecessors"); + else + ## make it a row vector + Tree = Tree(:)'; + + ## the count of nodes of the graph + NodNumber = length (Tree); + ## the number of children + ChildNumber = zeros (1, NodNumber + 1); + + + ## checking vector of predecessors + for i = 1 : NodNumber + if (Tree (i) < i) + ## this part of graph was checked before + continue; + endif + + ## Try to find cicle in this part of graph + ## we use modified Floyd's cycle-finding algorithm + tortoise = Tree (i); + hare = Tree (tortoise); + + while (tortoise != hare) + ## we end after find a cicle or when we reach a checked part of graph + + if (hare < i) + ## this part of graph was checked before + break + endif + + tortoise = Tree (tortoise); + ## hare will move faster than tortoise so in cicle hare + ## must reach tortoise + hare = Tree (Tree (hare)); + + endwhile + + if (tortoise == hare) + ## if hare reach tortoise we find cicle + error ("treelayout: vector of predecessors has bad format"); + endif + + endfor + ## vector of predecessors has right format + + for i = 1:NodNumber + ## VecOfChild is helping vector which is used to speed up the + ## choose of descendant nodes + + ChildNumber (Tree (i) + 1) = ChildNumber (Tree (i) + 1) + 1; + endfor + + Pos = 1; + for i = 1 : NodNumber + 1 + Start (i) = Pos; + Help (i) = Pos; + Pos += ChildNumber (i); + Stop (i) = Pos; + endfor + + if (nargin == 1) + for i = 1 : NodNumber + VecOfChild (Help (Tree (i) + 1)) = i; + Help (Tree (i) + 1) = Help (Tree (i) + 1) + 1; + endfor + else + VecOfChild = Permutation; + endif + + + ## the number of "parent" (actual) node (it's descendants will be + ## browse in the next iteration) + ParNumber = 0; + + ## the x-coordinate of the left most descendant of "parent node" + ## this value is increased in each leaf + LeftMost = 0; + + ## the level of "parent" node (root level is NodNumber) + Level = NodNumber; + + ## NodNumber - Max is the height of this graph + Max = NodNumber; + + ## main stack - each item consists of two numbers - the number of + ## node and the number it's of parent node on the top of stack + ## there is "parent node" + St = [-1, 0]; + + #number of vertices s in the top-level separator + s = 0; + # flag which says if we are in top level separator + topLevel = 1; + ## the top of the stack + while (ParNumber != -1) + if (Start(ParNumber + 1) < Stop(ParNumber + 1)) + idx = VecOfChild (Start (ParNumber + 1) : Stop (ParNumber + 1) - 1); + else + idx = zeros (1, 0); + endif + + ## add to idx the vector of parent descendants + St = [St ; [idx', ones(fliplr(size(idx))) * ParNumber]]; + + # we are in top level separator when we have one children + ## and the flag is 1 + if (columns(idx) == 1 && topLevel ==1 ) + s += 1; + else + # we arent in top level separator now + topLevel = 0; + endif + ## if there is not any descendant of "parent node": + if (St(end,2) != ParNumber) + LeftMost = LeftMost + 1; + XCoordinateR(ParNumber) = LeftMost; + Max = min (Max, Level); + if ((length(St) > 1) && (find((shift(St,1)-St) == 0) >1) + && St(end,2) != St(end-1,2)) + ## return to the nearest branching the position to return + ## position is the position on the stack, where should be + ## started further search (there are two nodes which has the + ## same parent node) + + Position = (find ((shift (St(:, 2), 1) - St(:, 2)) == 0))(end)+1; + ParNumberVec = St(Position : end, 2); + + ## the vector of removed nodes (the content of stack form + ## position to end) + + Level = Level + length(ParNumberVec); + + ## the level have to be decreased + + XCoordinateR(ParNumberVec) = LeftMost; + St(Position:end, :) = []; + endif + + ## remove the next node from "searched branch" + + St(end, :) = []; + ## choose new "parent node" + ParNumber = St(end, 1); + ## if there is another branch start to search it + if (ParNumber != -1) + YCoordinate(ParNumber) = Level; + XCoordinateL(ParNumber) = LeftMost + 1; + endif + else + + ## there were descendants of "parent nod" choose the last of + ## them and go on through it + Level--; + ParNumber = St(end, 1); + YCoordinate(ParNumber) = Level; + XCoordinateL(ParNumber) = LeftMost+1; + endif + endwhile + + ## calculate the x coordinates (the known values are the position + ## of most left and most right descendants) + XCoordinate = (XCoordinateL + XCoordinateR) / 2; + + Height = NodNumber - Max - 1; + endif +endfunction + +%!demo +%! % Compute a simple tree layout +%! [x,y,h,s]=treelayout([0 1 2 2]) + +%!demo +%! % Compute a simple tree layout with defined postorder permutation +%! [x,y,h,s]=treelayout([0 1 2 2],[1 2 3 4])