new file mode 100644
--- /dev/null
+++ b/scripts/sparse/pcr.m
@@ -0,0 +1,437 @@
+## Copyright (C) 2004 Piotr Krzyzanowski <piotr.krzyzanowski@mimuw.edu.pl>
+##
+## 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 2, 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, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{x} =} pcr (@var{A}, @var{b}, @var{tol}, @var{maxit}, @var{M}, @var{x0}, @dots{})
+## @deftypefnx {Function File} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}] =} pcr (@dots{})
+## 
+## Solves the linear system of equations @code{@var{A} * @var{x} =
+## @var{b}} by means of the  Preconditioned Conjugate Residuals iterative
+## method. The input arguments are
+##
+## @itemize
+## @item
+## @var{A} can be either a square (preferably sparse) matrix or a
+## function handle, inline function or string containing the name
+## of a function which computes @code{@var{A} * @var{x}}. In principle
+## @var{A} should be symmetric and non-singular; if @code{pcr}
+## finds @var{A} to be numerically singular, you will get a warning
+## message and the @var{flag} output parameter will be set.
+## 
+## @item
+## @var{b} is the right hand side vector.
+## 
+## @item
+## @var{tol} is the required relative tolerance for the residual error,
+## @code{@var{b} - @var{A} * @var{x}}. The iteration stops if @code{norm
+## (@var{b} - @var{A} * @var{x}) <= @var{tol} * norm (@var{b} - @var{A} *
+## @var{x0})}. If @var{tol} is empty or is omitted, the function sets
+## @code{@var{tol} = 1e-6} by default.
+## 
+## @item
+## @var{maxit} is the maximum allowable number of iterations; if
+## @code{[]} is supplied for @code{maxit}, or @code{pcr} has less
+## arguments, a default value equal to 20 is used.
+##
+## @item
+## @var{M} is the (left) preconditioning matrix, so that the iteration is
+## (theoretically) equivalent to solving by @code{pcr} @code{@var{P} *
+## @var{x} = @var{M} \ @var{b}}, with @code{@var{P} = @var{M} \ @var{A}}.
+## Note that a proper choice of the preconditioner may dramatically
+## improve the overall performance of the method. Instead of matrix
+## @var{M}, the user may pass a function which returns the results of 
+## applying the inverse of @var{M} to a vector (usually this is the
+## preferred way of using the preconditioner). If @code{[]} is supplied
+## for @var{M}, or @var{M} is omitted, no preconditioning is applied.
+## 
+## @item
+## @var{x0} is the initial guess. If @var{x0} is empty or omitted, the 
+## function sets @var{x0} to a zero vector by default.
+## @end itemize
+## 
+## The arguments which follow @var{x0} are treated as parameters, and
+## passed in a proper way to any of the functions (@var{A} or @var{M})
+## which are passed to @code{pcr}. See the examples below for further
+## details. The output arguments are
+##
+## @itemize
+## @item
+## @var{x} is the computed approximation to the solution of
+## @code{@var{A} * @var{x} = @var{b}}.
+## 
+## @item
+## @var{flag} reports on the convergence. @code{@var{flag} = 0} means
+## the solution converged and the tolerance criterion given by @var{tol}
+## is satisfied. @code{@var{flag} = 1} means that the @var{maxit} limit
+## for the iteration count was reached. @code{@var{flag} = 3} reports t
+## @code{pcr} breakdown, see [1] for details.
+## 
+## @item
+## @var{relres} is the ratio of the final residual to its initial value,
+## measured in the Euclidean norm.
+## 
+## @item
+## @var{iter} is the actual number of iterations performed.
+##
+## @item 
+## @var{resvec} describes the convergence history of the method,
+## so that @code{@var{resvec} (i)} contains the Euclidean norms of the 
+## residualafter the (@var{i}-1)-th iteration, @code{@var{i} =
+## 1,2,...@var{iter}+1}.
+## @end itemize
+## 
+## Let us consider a trivial problem with a diagonal matrix (we exploit the
+## sparsity of A) 
+## 
+## @example
+## @group
+## 	N = 10; 
+## 	A = diag([1:N]); A = sparse(A);  
+## 	b = rand(N,1);
+## @end group
+## @end example
+## 
+## @sc{Example 1:} Simplest use of @code{pcr}
+## 
+## @example
+##   x = pcr(A, b)
+## @end example
+## 
+## @sc{Example 2:} @code{pcr} with a function which computes
+## @code{@var{A} * @var{x}}.
+##
+## @example
+## @group
+##   function y = applyA(x) 
+##     y = [1:10]'.*x; 
+##   endfunction
+## 
+##   x = pcr('applyA',b)
+## @end group
+## @end example
+## 
+## @sc{Example 3:}  Preconditioned iteration, with full diagnostics. The
+## preconditioner (quite strange, because even the original matrix
+## @var{A} is trivial) is defined as a function
+## 
+## @example
+## @group
+##   function y = applyM(x)		
+##     K = floor(length(x)-2); 
+##     y = x; 
+##     y(1:K) = x(1:K)./[1:K]';	
+##   endfunction
+## 
+##   [x, flag, relres, iter, resvec] = pcr(A,b,[],[],'applyM')
+##   semilogy([1:iter+1], resvec);
+## @end group
+## @end example
+##
+## @sc{Example 4:} Finally, a preconditioner which depends on a
+## parameter @var{k}.
+## 
+## @example
+## @group
+##   function y = applyM(x, varargin)
+##     K = varargin@{1@}; 
+##     y = x; y(1:K) = x(1:K)./[1:K]';	 
+##   endfunction
+## 
+##   [x, flag, relres, iter, resvec] = pcr(A,b,[],[],'applyM',[],3)
+## @end group
+## @end example
+## 
+## @sc{References}
+## 
+## 	[1] W. Hackbusch, "Iterative Solution of Large Sparse Systems of
+##  	Equations", section 9.5.4; Springer, 1994
+##
+## @seealso{sparse, pcg}
+## @end deftypefn
+
+## REVISION HISTORY
+##
+## 2004-08-14, Piotr Krzyzanowski <piotr.krzyzanowski@mimuw.edu.pl>
+## 
+## 	Added 4 demos and 4 tests
+##  
+## 2004-08-01, Piotr Krzyzanowski <piotr.krzyzanowski@mimuw.edu.pl>
+## 
+## 	First version of pcr code
+
+function [x, flag, relres, iter, resvec] = pcr(A,b,tol,maxit,M,x0,varargin)
+
+  breakdown = false;
+
+  if ((nargin < 6) || isempty(x0))
+    x = zeros(size(b));
+  else
+    x = x0;
+  endif
+
+  if (nargin < 5)
+    M = [];
+  endif
+
+  if ((nargin < 4) || isempty(maxit))
+    maxit = 20;
+  endif
+
+  maxit = maxit + 2;
+
+  if ((nargin < 3) || isempty(tol))
+    tol = 1e-6;
+  endif
+
+  if (nargin < 2)
+    print_usage();
+  endif
+
+  ##  init
+  if (isnumeric(A))		# is A a matrix?
+    r = b-A*x;
+  else				# then A should be a function!
+    r = b-feval(A,x,varargin{:});
+  endif
+
+  if (isnumeric(M))		# is M a matrix?
+    if isempty(M)		# if M is empty, use no precond
+      p = r;
+    else			# otherwise, apply the precond
+      p = M\r;
+    endif
+  else				# then M should be a function!
+    p = feval(M,r,varargin{:});
+  endif
+
+  iter = 2;
+
+  b_bot_old = 1;
+  q_old = p_old = s_old = zeros(size(x));
+
+  if (isnumeric(A))		# is A a matrix?
+    q = A*p;
+  else				# then A should be a function!
+    q = feval(A,p,varargin{:});
+  endif
+	
+  resvec(1) = abs(norm(r)); 
+
+  ## iteration
+  while ((resvec(iter-1) > tol*resvec(1)) && (iter < maxit))
+	
+    if (isnumeric(M))		# is M a matrix?
+      if isempty(M)		# if M is empty, use no precond
+	s = q;
+      else			# otherwise, apply the precond
+	s = M\q;
+      endif
+    else			# then M should be a function!
+      s = feval(M,q,varargin{:});
+    endif
+    b_top = r'*s;
+    b_bot = q'*s;
+	
+    if (b_bot == 0.0)
+      breakdown = true;
+      break;
+    endif
+    lambda = b_top/b_bot;
+	
+    x = x + lambda*p;
+    r = r - lambda*q;
+	
+    if (isnumeric(A))		# is A a matrix?
+      t = A*s;
+    else			# then A should be a function!
+      t = feval(A,s,varargin{:});
+    endif
+	
+    alpha0 = (t'*s)/b_bot;
+    alpha1 = (t'*s_old)/b_bot_old;
+	
+    p_temp = p; q_temp = q;
+    p = s - alpha0*p - alpha1*p_old;
+    q = t - alpha0*q - alpha1*q_old;
+	
+    s_old = s;
+    p_old = p_temp;
+    q_old = q_temp;
+    b_bot_old = b_bot;
+	
+	
+    resvec(iter) = abs(norm(r));
+    iter = iter + 1;
+  endwhile
+
+  flag = 0;
+  relres = resvec(iter-1)./resvec(1);
+  iter = iter - 2;
+  if (iter >= (maxit-2))
+    flag = 1;
+    if (nargout < 2)
+      warning("PCR: maximum number of iterations (%d) reached\n", iter);
+      warning("The initial residual norm was reduced %g times.\n", 1.0/relres);
+    endif
+  else
+    if ((nargout < 2) && (~breakdown))
+      fprintf(stderr, "PCR: converged in %d iterations. \n", iter);
+      fprintf(stderr, "The initial residual norm was reduced %g times.\n",...
+	      1.0/relres);
+    endif
+  endif
+  if (breakdown)
+    flag = 3;
+    if (nargout < 2)
+      warning("PCR: breakdown occured.\n");
+      warning("System matrix singular or preconditioner indefinite?\n");
+    endif
+  endif
+
+endfunction
+
+%!demo
+%!
+%!	# Simplest usage of PCR (see also 'help pcr')
+%!
+%!	N = 20; 
+%!	A = diag(linspace(-3.1,3,N)); b = rand(N,1); y = A\b; #y is the true solution
+%!  	x = pcr(A,b);
+%!	printf('The solution relative error is %g\n', norm(x-y)/norm(y));
+%!
+%!	# You shouldn't be afraid if PCR issues some warning messages in this
+%!	# example: watch out in the second example, why it takes N iterations 
+%!	# of PCR to converge to (a very accurate, by the way) solution
+%!demo
+%!
+%!	# Full output from PCR
+%!	# We use this output to plot the convergence history  
+%!
+%!	N = 20; 
+%!	A = diag(linspace(-3.1,30,N)); b = rand(N,1); X = A\b; #X is the true solution
+%!  	[x, flag, relres, iter, resvec] = pcr(A,b);
+%!	printf('The solution relative error is %g\n', norm(x-X)/norm(X));
+%!	title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||/||b||)');
+%!	semilogy([0:iter],resvec/resvec(1),'o-g;relative residual;');
+%!demo
+%!
+%!	# Full output from PCR
+%!	# We use indefinite matrix based on the Hilbert matrix, with one
+%!	# strongly negative eigenvalue
+%!	# Hilbert matrix is extremely ill conditioned, so is ours, 
+%!	# and that's why PCR WILL have problems
+%!
+%!	N = 10; 
+%!	A = hilb(N); A(1,1)=-A(1,1); b = rand(N,1); X = A\b; #X is the true solution
+%!	printf('Condition number of A is   %g\n', cond(A));
+%!  	[x, flag, relres, iter, resvec] = pcr(A,b,[],200);
+%!	if (flag == 3)
+%!	  printf('PCR breakdown. System matrix is [close to] singular\n');
+%!	end
+%!	title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||)');
+%!	semilogy([0:iter],resvec,'o-g;absolute residual;');
+%!demo
+%!
+%!	# Full output from PCR
+%!	# We use an indefinite matrix based on the 1-D Laplacian matrix for A, 
+%!	# and here we have cond(A) = O(N^2)
+%!	# That's the reason we need some preconditioner; here we take
+%!	# a very simple and not powerful Jacobi preconditioner, 
+%!	# which is the diagonal of A
+%!
+%!	# Note that we use here indefinite preconditioners!
+%!
+%!	N = 100; 
+%!	A = zeros(N,N);
+%!	for i=1:N-1 # form 1-D Laplacian matrix
+%!		A(i:i+1,i:i+1) = [2 -1; -1 2];
+%!	endfor
+%!	A = [A, zeros(size(A)); zeros(size(A)), -A];
+%!	b = rand(2*N,1); X = A\b; #X is the true solution
+%!	maxit = 80;
+%!	printf('System condition number is %g\n',cond(A));
+%!	# No preconditioner: the convergence is very slow!
+%!
+%!  	[x, flag, relres, iter, resvec] = pcr(A,b,[],maxit);
+%!	title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||)');
+%!	semilogy([0:iter],resvec,'o-g;NO preconditioning: absolute residual;');
+%!
+%!	pause(1);
+%!	# Test Jacobi preconditioner: it will not help much!!!
+%!
+%!	M = diag(diag(A)); # Jacobi preconditioner
+%!  	[x, flag, relres, iter, resvec] = pcr(A,b,[],maxit,M);
+%!	hold on;
+%!	semilogy([0:iter],resvec,'o-r;JACOBI preconditioner: absolute residual;');
+%!
+%!	pause(1);
+%!	# Test nonoverlapping block Jacobi preconditioner: this one should give
+%!	# some convergence speedup!
+%!
+%!	M = zeros(N,N);k=4;
+%!	for i=1:k:N # get k x k diagonal blocks of A
+%!		M(i:i+k-1,i:i+k-1) = A(i:i+k-1,i:i+k-1);
+%!	endfor
+%!	M = [M, zeros(size(M)); zeros(size(M)), -M];
+%!  	[x, flag, relres, iter, resvec] = pcr(A,b,[],maxit,M);
+%!	semilogy([0:iter],resvec,'o-b;BLOCK JACOBI preconditioner: absolute residual;');
+%!	hold off;
+%!test
+%!
+%!	#solve small indefinite diagonal system
+%!
+%!	N = 10; 
+%!	A = diag(linspace(-10.1,10,N)); b = ones(N,1); X = A\b; #X is the true solution
+%!  	[x, flag] = pcr(A,b,[],N+1);
+%!	assert(norm(x-X)/norm(X)<1e-10);
+%!	assert(flag,0);
+%!
+%!test
+%!
+%!	#solve tridiagonal system, do not converge in default 20 iterations
+%!	#should perform max allowable default number of iterations
+%!
+%!	N = 100; 
+%!	A = zeros(N,N);
+%!	for i=1:N-1 # form 1-D Laplacian matrix
+%!		A(i:i+1,i:i+1) = [2 -1; -1 2];
+%!	endfor
+%!	b = ones(N,1); X = A\b; #X is the true solution
+%!  	[x, flag, relres, iter, resvec] = pcr(A,b,1e-12);
+%!	assert(flag,1);
+%!	assert(relres>0.6);
+%!	assert(iter,20);
+%!
+%!test
+%!
+%!	#solve tridiagonal system with 'prefect' preconditioner
+%!	#converges in one iteration
+%!
+%!	N = 100; 
+%!	A = zeros(N,N);
+%!	for i=1:N-1 # form 1-D Laplacian matrix
+%!		A(i:i+1,i:i+1) = [2 -1; -1 2];
+%!	endfor
+%!	b = ones(N,1); X = A\b; #X is the true solution
+%!  	[x, flag, relres, iter] = pcr(A,b,[],[],A,b);
+%!	assert(norm(x-X)/norm(X)<1e-6);
+%!	assert(relres<1e-6);
+%!	assert(flag,0);
+%!	assert(iter,1); #should converge in one iteration
+%!