--- a/doc/interpreter/contributors.in	Wed Nov 19 16:55:47 2008 +0100
+++ b/doc/interpreter/contributors.in	Fri Nov 21 15:03:03 2008 +0100
@@ -170,6 +170,7 @@
 Olli Saarela
 Toni Saarela
 Juhani Saastamoinen
+Radek Salac
 Ben Sapp
 Alois Schloegl
 Michel D. Schmid

--- a/scripts/ChangeLog	Wed Nov 19 16:55:47 2008 +0100
+++ b/scripts/ChangeLog	Fri Nov 21 15:03:03 2008 +0100
@@ -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	Wed Nov 19 16:55:47 2008 +0100
+++ b/scripts/sparse/Makefile.in	Fri Nov 21 15:03:03 2008 +0100
@@ -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))
 

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/sparse/cgs.m	Fri Nov 21 15:03:03 2008 +0100
@@ -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)

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/sparse/treelayout.m	Fri Nov 21 15:03:03 2008 +0100
@@ -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])