R! wrote:

> Well thanks for your opinion, I really appreciate it.

> If you think that you can improve my code in any way please do it,but if you

> don't have time I will try it myself. But it's your idea so I think that you

> have the right to implement it :-).

>

> Thanks for all.

>

> Radek Salac.

>

>

>

>

Ok, I committed the attached changeset credited to you.

D.

--

David Bateman

[hidden email]
Motorola Labs - Paris +33 1 69 35 48 04 (Ph)

Parc Les Algorithmes, Commune de St Aubin +33 6 72 01 06 33 (Mob)

91193 Gif-Sur-Yvette FRANCE +33 1 69 35 77 01 (Fax)

The information contained in this communication has been classified as:

[x] General Business Information

[ ] Motorola Internal Use Only

[ ] Motorola Confidential Proprietary

# HG changeset patch

# User Radek Salac <

[hidden email]>

# Date 1227276183 -3600

# Node ID 5d8461010f118881bd040c02e8bbebd6ba81ee77

# Parent 67d5e0fd3c7dffbc4fba9ada0a0399c423a14a25

Add the cgs and treelayout functions

diff --git a/doc/interpreter/contributors.in b/doc/interpreter/contributors.in

--- 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

diff --git a/scripts/ChangeLog b/scripts/ChangeLog

--- a/scripts/ChangeLog

+++ b/scripts/ChangeLog

@@ -1,3 +1,8 @@

+2008-11-21 Radek Salac <

[hidden email]>

+

+ * sparse/cgs.m, sparse/treelayout.m: New functions.

+ * sparse/Makefile.in (SOURCES): Add them here.

+

2008-11-14 David Bateman <

[hidden email]>

* plot/__go_draw_axes__.m (do_tics_1): Support the minorick properties

diff --git a/scripts/sparse/Makefile.in b/scripts/sparse/Makefile.in

--- 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))

diff --git a/scripts/sparse/cgs.m b/scripts/sparse/cgs.m

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)

diff --git a/scripts/sparse/treelayout.m b/scripts/sparse/treelayout.m

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])