/*

Copyright (C) 1996, 1997, 1999, 2000, 2005, 2006, 2007 John W. Eaton

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

*/

// The qrupdate, qrinsert, qrdelete and qrshift functions were written by
// Jaroslav Hajek <highegg@gmail.com>, Copyright (C) 2008  VZLU
// Prague, a.s., Czech Republic.

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include "CmplxQR.h"
#include "CmplxQRP.h"
#include "dbleQR.h"
#include "dbleQRP.h"
#include "fCmplxQR.h"
#include "fCmplxQRP.h"
#include "floatQR.h"
#include "floatQRP.h"
#include "SparseQR.h"
#include "SparseCmplxQR.h"


#include "defun-dld.h"
#include "error.h"
#include "gripes.h"
#include "oct-obj.h"
#include "utils.h"

// [Q, R] = qr (X):      form Q unitary and R upper triangular such
//                        that Q * R = X
//
// [Q, R] = qr (X, 0):    form the economy decomposition such that if X is
//                        m by n then only the first n columns of Q are
//                        computed.
//
// [Q, R, P] = qr (X):    form QRP factorization of X where
//                        P is a permutation matrix such that
//                        A * P = Q * R
//
// [Q, R, P] = qr (X, 0): form the economy decomposition with
//                        permutation vector P such that Q * R = X (:, P)
//
// qr (X) alone returns the output of the LAPACK routine dgeqrf, such
// that R = triu (qr (X))

DEFUN_DLD (qr, args, nargout,
  "-*- texinfo -*-\n\
@deftypefn {Loadable Function} {[@var{q}, @var{r}, @var{p}] =} qr (@var{a})\n\
@deftypefnx {Loadable Function} {[@var{q}, @var{r}, @var{p}] =} qr (@var{a}, '0')\n\
@cindex QR factorization\n\
Compute the QR factorization of @var{a}, using standard @sc{Lapack}\n\
subroutines.  For example, given the matrix @code{a = [1, 2; 3, 4]},\n\
\n\
@example\n\
[q, r] = qr (a)\n\
@end example\n\
\n\
@noindent\n\
returns\n\
\n\
@example\n\
q =\n\
\n\
  -0.31623  -0.94868\n\
  -0.94868   0.31623\n\
\n\
r =\n\
\n\
  -3.16228  -4.42719\n\
   0.00000  -0.63246\n\
@end example\n\
\n\
The @code{qr} factorization has applications in the solution of least\n\
squares problems\n\
@iftex\n\
@tex\n\
$$\n\
\\min_x \\left\\Vert A x - b \\right\\Vert_2\n\
$$\n\
@end tex\n\
@end iftex\n\
@ifinfo\n\
\n\
@example\n\
@code{min norm(A x - b)}\n\
@end example\n\
\n\
@end ifinfo\n\
for overdetermined systems of equations (i.e.,\n\
@iftex\n\
@tex\n\
$A$\n\
@end tex\n\
@end iftex\n\
@ifinfo\n\
@code{a}\n\
@end ifinfo\n\
 is a tall, thin matrix).  The QR factorization is\n\
@iftex\n\
@tex\n\
$QR = A$ where $Q$ is an orthogonal matrix and $R$ is upper triangular.\n\
@end tex\n\
@end iftex\n\
@ifinfo\n\
@code{q * r = a} where @code{q} is an orthogonal matrix and @code{r} is\n\
upper triangular.\n\
@end ifinfo\n\
\n\
If given a second argument of '0', @code{qr} returns an economy-sized\n\
QR factorization, omitting zero rows of @var{R} and the corresponding\n\
columns of @var{Q}.\n\
\n\
If the matrix @var{a} is full, the permuted QR factorization\n\
@code{[@var{q}, @var{r}, @var{p}] = qr (@var{a})} forms the QR factorization\n\
such that the diagonal entries of @code{r} are decreasing in magnitude\n\
order.  For example,given the matrix @code{a = [1, 2; 3, 4]},\n\
\n\
@example\n\
[q, r, p] = qr(a)\n\
@end example\n\
\n\
@noindent\n\
returns\n\
\n\
@example\n\
q = \n\
\n\
  -0.44721  -0.89443\n\
  -0.89443   0.44721\n\
\n\
r =\n\
\n\
  -4.47214  -3.13050\n\
   0.00000   0.44721\n\
\n\
p =\n\
\n\
   0  1\n\
   1  0\n\
@end example\n\
\n\
The permuted @code{qr} factorization @code{[q, r, p] = qr (a)}\n\
factorization allows the construction of an orthogonal basis of\n\
@code{span (a)}.\n\
\n\
If the matrix @var{a} is sparse, then compute the sparse QR factorization\n\
of @var{a}, using @sc{CSparse}. As the matrix @var{Q} is in general a full\n\
matrix, this function returns the @var{Q}-less factorization @var{r} of\n\
@var{a}, such that @code{@var{r} = chol (@var{a}' * @var{a})}.\n\
\n\
If the final argument is the scalar @code{0} and the number of rows is\n\
larger than the number of columns, then an economy factorization is\n\
returned. That is @var{r} will have only @code{size (@var{a},1)} rows.\n\
\n\
If an additional matrix @var{b} is supplied, then @code{qr} returns\n\
@var{c}, where @code{@var{c} = @var{q}' * @var{b}}. This allows the\n\
least squares approximation of @code{@var{a} \\ @var{b}} to be calculated\n\
as\n\
\n\
@example\n\
[@var{c},@var{r}] = spqr (@var{a},@var{b})\n\
@var{x} = @var{r} \\ @var{c}\n\
@end example\n\
@end deftypefn")
{
  octave_value_list retval;

  int nargin = args.length ();

  if (nargin < 1 || nargin > (args(0).is_sparse_type() ? 3 : 2))
    {
      print_usage ();
      return retval;
    }

  octave_value arg = args(0);

  int arg_is_empty = empty_arg ("qr", arg.rows (), arg.columns ());

  if (arg_is_empty < 0)
    return retval;
  else if (arg_is_empty > 0)
    return octave_value_list (3, Matrix ());

  if (arg.is_sparse_type ())
    {
      bool economy = false;
      bool is_cmplx = false;
      int have_b = 0;

      if (arg.is_complex_type ())
	is_cmplx = true;
      if (nargin > 1)
	{
	  have_b = 1;
	  if (args(nargin-1).is_scalar_type ())
	    {
	      int val = args(nargin-1).int_value ();
	      if (val == 0)
		{
		  economy = true;
		  have_b = (nargin > 2 ? 2 : 0);
		}
	    }
	  if (have_b > 0 && args(have_b).is_complex_type ())
	    is_cmplx = true;
	}
	
      if (!error_state)
	{
	  if (have_b && nargout < 2)
	    error ("qr: incorrect number of output arguments");
	  else if (is_cmplx)
	    {
	      SparseComplexQR q (arg.sparse_complex_matrix_value ());
	      if (!error_state)
		{
		  if (have_b > 0)
		    {
		      retval(1) = q.R (economy);
		      retval(0) = q.C (args(have_b).complex_matrix_value ());
		      if (arg.rows() < arg.columns())
			warning ("qr: non minimum norm solution for under-determined problem");
		    }
		  else if (nargout > 1)
		    {
		      retval(1) = q.R (economy);
		      retval(0) = q.Q ();
		    }
		  else
		    retval(0) = q.R (economy);
		}
	    }
	  else
	    {
	      SparseQR q (arg.sparse_matrix_value ());
	      if (!error_state)
		{
		  if (have_b > 0)
		    {
		      retval(1) = q.R (economy);
		      retval(0) = q.C (args(have_b).matrix_value ());
		      if (args(0).rows() < args(0).columns())
			warning ("qr: non minimum norm solution for under-determined problem");
		    }
		  else if (nargout > 1)
		    {
		      retval(1) = q.R (economy);
		      retval(0) = q.Q ();
		    }
		  else
		    retval(0) = q.R (economy);
		}
	    }
	}
    }
  else
    {
      QR::type type = (nargout == 0 || nargout == 1) ? QR::raw
	: (nargin == 2 ? QR::economy : QR::std);

      if (arg.is_single_type ())
	{
	  if (arg.is_real_type ())
	    {
	      FloatMatrix m = arg.float_matrix_value ();

	      if (! error_state)
		{
		  switch (nargout)
		    {
		    case 0:
		    case 1:
		      {
			FloatQR fact (m, type);
			retval(0) = fact.R ();
		      }
		      break;

		    case 2:
		      {
			FloatQR fact (m, type);
			retval(1) = fact.R ();
			retval(0) = fact.Q ();
		      }
		      break;

		    default:
		      {
			FloatQRP fact (m, type);
			retval(2) = fact.P ();
			retval(1) = fact.R ();
			retval(0) = fact.Q ();
		      }
		      break;
		    }
		}
	    }
	  else if (arg.is_complex_type ())
	    {
	      FloatComplexMatrix m = arg.float_complex_matrix_value ();

	      if (! error_state)
		{
		  switch (nargout)
		    {
		    case 0:
		    case 1:
		      {
			FloatComplexQR fact (m, type);
			retval(0) = fact.R ();
		      }
		      break;

		    case 2:
		      {
			FloatComplexQR fact (m, type);
			retval(1) = fact.R ();
			retval(0) = fact.Q ();
		      }
		      break;

		    default:
		      {
			FloatComplexQRP fact (m, type);
			retval(2) = fact.P ();
			retval(1) = fact.R ();
			retval(0) = fact.Q ();
		      }
		      break;
		    }
		}
	    }
	}
      else
	{
	  if (arg.is_real_type ())
	    {
	      Matrix m = arg.matrix_value ();

	      if (! error_state)
		{
		  switch (nargout)
		    {
		    case 0:
		    case 1:
		      {
			QR fact (m, type);
			retval(0) = fact.R ();
		      }
		      break;

		    case 2:
		      {
			QR fact (m, type);
			retval(1) = fact.R ();
			retval(0) = fact.Q ();
		      }
		      break;

		    default:
		      {
			QRP fact (m, type);
			retval(2) = fact.P ();
			retval(1) = fact.R ();
			retval(0) = fact.Q ();
		      }
		      break;
		    }
		}
	    }
	  else if (arg.is_complex_type ())
	    {
	      ComplexMatrix m = arg.complex_matrix_value ();

	      if (! error_state)
		{
		  switch (nargout)
		    {
		    case 0:
		    case 1:
		      {
			ComplexQR fact (m, type);
			retval(0) = fact.R ();
		      }
		      break;

		    case 2:
		      {
			ComplexQR fact (m, type);
			retval(1) = fact.R ();
			retval(0) = fact.Q ();
		      }
		      break;

		    default:
		      {
			ComplexQRP fact (m, type);
			retval(2) = fact.P ();
			retval(1) = fact.R ();
			retval(0) = fact.Q ();
		      }
		      break;
		    }
		}
	    }
	  else
	    gripe_wrong_type_arg ("qr", arg);
	}
    }

  return retval;
}

/*

%!test
%! a = [0, 2, 1; 2, 1, 2];
%! 
%! [q, r] = qr (a);
%! 
%! [qe, re] = qr (a, 0);
%! 
%! assert (q * r, a, sqrt (eps));
%! assert (qe * re, a, sqrt (eps));

%!test
%! a = [0, 2, 1; 2, 1, 2];
%! 
%! [q, r, p] = qr (a);  # not giving right dimensions. FIXME
%! 
%! [qe, re, pe] = qr (a, 0);
%! 
%! assert (q * r, a * p, sqrt (eps));
%! assert (qe * re, a(:, pe), sqrt (eps));

%!test
%! a = [0, 2; 2, 1; 1, 2];
%! 
%! [q, r] = qr (a);
%! 
%! [qe, re] = qr (a, 0);
%! 
%! assert (q * r, a, sqrt (eps));
%! assert (qe * re, a, sqrt (eps));

%!test
%! a = [0, 2; 2, 1; 1, 2];
%! 
%! [q, r, p] = qr (a);
%! 
%! [qe, re, pe] = qr (a, 0);
%! 
%! assert (q * r, a * p, sqrt (eps));
%! assert (qe * re, a(:, pe), sqrt (eps));

%!error <Invalid call to qr.*> qr ();
%!error <Invalid call to qr.*> qr ([1, 2; 3, 4], 0, 2);

%!function retval = testqr (q, r, a, p)
%!  tol = 100*eps (class(q));
%!  retval = 0;
%!  if (nargin == 3)
%!    n1 = norm (q*r-a);
%!    n2 = norm (q'*q-eye(columns(q)));
%!    retval = (n1 < tol && n2 < tol);
%!  else
%!    n1 = norm (q'*q-eye(columns(q)));
%!    retval = (n1 < tol);
%!    if (isvector (p))
%!      n2 = norm (q*r-a(:,p));
%!      retval = (retval && n2 < tol);
%!    else
%!      n2 = norm (q*r - a*p);
%!      retval = (retval && n2 < tol);
%!    endif
%!  endif
%!test
%! 
%! t = ones (24, 1);
%! j = 1;
%! 
%! if false # eliminate big matrix tests
%!   a = rand(5000,20);
%!   [q,r]=qr(a,0); t(j++) = testqr(q,r,a);
%!   [q,r]=qr(a',0); t(j++) = testqr(q,r,a');
%!   [q,r,p]=qr(a,0); t(j++) = testqr(q,r,a,p);
%!   [q,r,p]=qr(a',0); t(j++) = testqr(q,r,a',p);
%! 
%!   a = a+1i*eps;
%!   [q,r]=qr(a,0); t(j++) = testqr(q,r,a);
%!   [q,r]=qr(a',0); t(j++) = testqr(q,r,a');
%!   [q,r,p]=qr(a,0); t(j++) = testqr(q,r,a,p);
%!   [q,r,p]=qr(a',0); t(j++) = testqr(q,r,a',p);
%! endif
%! 
%! a = [ ones(1,15); sqrt(eps)*eye(15) ];
%! [q,r]=qr(a); t(j++) = testqr(q,r,a);
%! [q,r]=qr(a'); t(j++) = testqr(q,r,a');
%! [q,r,p]=qr(a); t(j++) = testqr(q,r,a,p);
%! [q,r,p]=qr(a'); t(j++) = testqr(q,r,a',p);
%! 
%! a = a+1i*eps;
%! [q,r]=qr(a); t(j++) = testqr(q,r,a);
%! [q,r]=qr(a'); t(j++) = testqr(q,r,a');
%! [q,r,p]=qr(a); t(j++) = testqr(q,r,a,p);
%! [q,r,p]=qr(a'); t(j++) = testqr(q,r,a',p);
%! 
%! a = [ ones(1,15); sqrt(eps)*eye(15) ];
%! [q,r]=qr(a,0); t(j++) = testqr(q,r,a);
%! [q,r]=qr(a',0); t(j++) = testqr(q,r,a');
%! [q,r,p]=qr(a,0); t(j++) = testqr(q,r,a,p);
%! [q,r,p]=qr(a',0); t(j++) = testqr(q,r,a',p);
%! 
%! a = a+1i*eps;
%! [q,r]=qr(a,0); t(j++) = testqr(q,r,a);
%! [q,r]=qr(a',0); t(j++) = testqr(q,r,a');
%! [q,r,p]=qr(a,0); t(j++) = testqr(q,r,a,p);
%! [q,r,p]=qr(a',0); t(j++) = testqr(q,r,a',p);
%! 
%! a = [
%! 611   196  -192   407    -8   -52   -49    29
%! 196   899   113  -192   -71   -43    -8   -44
%! -192   113   899   196    61    49     8    52
%! 407  -192   196   611     8    44    59   -23
%! -8   -71    61     8   411  -599   208   208
%! -52   -43    49    44  -599   411   208   208
%! -49    -8     8    59   208   208    99  -911
%! 29   -44    52   -23   208   208  -911    99
%! ];
%! [q,r] = qr(a);
%! 
%! assert(all (t) && norm(q*r-a) < 5000*eps);

%!test
%! a = single ([0, 2, 1; 2, 1, 2]);
%! 
%! [q, r] = qr (a);
%! 
%! [qe, re] = qr (a, 0);
%! 
%! assert (q * r, a, sqrt (eps ('single')));
%! assert (qe * re, a, sqrt (eps ('single')));

%!test
%! a = single([0, 2, 1; 2, 1, 2]);
%! 
%! [q, r, p] = qr (a);  # not giving right dimensions. FIXME
%! 
%! [qe, re, pe] = qr (a, 0);
%! 
%! assert (q * r, a * p, sqrt (eps('single')));
%! assert (qe * re, a(:, pe), sqrt (eps('single')));

%!test
%! a = single([0, 2; 2, 1; 1, 2]);
%! 
%! [q, r] = qr (a);
%! 
%! [qe, re] = qr (a, 0);
%! 
%! assert (q * r, a, sqrt (eps('single')));
%! assert (qe * re, a, sqrt (eps('single')));

%!test
%! a = single([0, 2; 2, 1; 1, 2]);
%! 
%! [q, r, p] = qr (a);
%! 
%! [qe, re, pe] = qr (a, 0);
%! 
%! assert (q * r, a * p, sqrt (eps('single')));
%! assert (qe * re, a(:, pe), sqrt (eps('single')));

%!error <Invalid call to qr.*> qr ();
%!error <Invalid call to qr.*> qr ([1, 2; 3, 4], 0, 2);

%!test
%! 
%! t = ones (24, 1);
%! j = 1;
%! 
%! if false # eliminate big matrix tests
%!   a = rand(5000,20);
%!   [q,r]=qr(a,0); t(j++) = testqr(q,r,a);
%!   [q,r]=qr(a',0); t(j++) = testqr(q,r,a');
%!   [q,r,p]=qr(a,0); t(j++) = testqr(q,r,a,p);
%!   [q,r,p]=qr(a',0); t(j++) = testqr(q,r,a',p);
%! 
%!   a = a+1i*eps('single');
%!   [q,r]=qr(a,0); t(j++) = testqr(q,r,a);
%!   [q,r]=qr(a',0); t(j++) = testqr(q,r,a');
%!   [q,r,p]=qr(a,0); t(j++) = testqr(q,r,a,p);
%!   [q,r,p]=qr(a',0); t(j++) = testqr(q,r,a',p);
%! endif
%! 
%! a = [ ones(1,15); sqrt(eps('single'))*eye(15) ];
%! [q,r]=qr(a); t(j++) = testqr(q,r,a);
%! [q,r]=qr(a'); t(j++) = testqr(q,r,a');
%! [q,r,p]=qr(a); t(j++) = testqr(q,r,a,p);
%! [q,r,p]=qr(a'); t(j++) = testqr(q,r,a',p);
%! 
%! a = a+1i*eps('single');
%! [q,r]=qr(a); t(j++) = testqr(q,r,a);
%! [q,r]=qr(a'); t(j++) = testqr(q,r,a');
%! [q,r,p]=qr(a); t(j++) = testqr(q,r,a,p);
%! [q,r,p]=qr(a'); t(j++) = testqr(q,r,a',p);
%! 
%! a = [ ones(1,15); sqrt(eps('single'))*eye(15) ];
%! [q,r]=qr(a,0); t(j++) = testqr(q,r,a);
%! [q,r]=qr(a',0); t(j++) = testqr(q,r,a');
%! [q,r,p]=qr(a,0); t(j++) = testqr(q,r,a,p);
%! [q,r,p]=qr(a',0); t(j++) = testqr(q,r,a',p);
%! 
%! a = a+1i*eps('single');
%! [q,r]=qr(a,0); t(j++) = testqr(q,r,a);
%! [q,r]=qr(a',0); t(j++) = testqr(q,r,a');
%! [q,r,p]=qr(a,0); t(j++) = testqr(q,r,a,p);
%! [q,r,p]=qr(a',0); t(j++) = testqr(q,r,a',p);
%! 
%! a = [
%! 611   196  -192   407    -8   -52   -49    29
%! 196   899   113  -192   -71   -43    -8   -44
%! -192   113   899   196    61    49     8    52
%! 407  -192   196   611     8    44    59   -23
%! -8   -71    61     8   411  -599   208   208
%! -52   -43    49    44  -599   411   208   208
%! -49    -8     8    59   208   208    99  -911
%! 29   -44    52   -23   208   208  -911    99
%! ];
%! [q,r] = qr(a);
%! 
%! assert(all (t) && norm(q*r-a) < 5000*eps('single'));

%% The deactivated tests below can't be tested till rectangular back-subs is
%% implemented for sparse matrices.

%!testif HAVE_CXSPARSE
%! n = 20; d= 0.2;
%! a = sprandn(n,n,d)+speye(n,n);
%! r = qr(a);
%! assert(r'*r,a'*a,1e-10)

%!testif HAVE_CXSPARSE
%! n = 20; d= 0.2;
%! a = sprandn(n,n,d)+speye(n,n);
%! q = symamd(a);
%! a = a(q,q);
%! r = qr(a);
%! assert(r'*r,a'*a,1e-10)

%!testif HAVE_CXSPARSE
%! n = 20; d= 0.2;
%! a = sprandn(n,n,d)+speye(n,n);
%! [c,r] = qr(a,ones(n,1));
%! assert (r\c,full(a)\ones(n,1),10e-10)

%!testif HAVE_CXSPARSE
%! n = 20; d= 0.2;
%! a = sprandn(n,n,d)+speye(n,n);
%! b = randn(n,2);
%! [c,r] = qr(a,b);
%! assert (r\c,full(a)\b,10e-10)

%% Test under-determined systems!!
%!#testif HAVE_CXSPARSE
%! n = 20; d= 0.2;
%! a = sprandn(n,n+1,d)+speye(n,n+1);
%! b = randn(n,2);
%! [c,r] = qr(a,b);
%! assert (r\c,full(a)\b,10e-10)

%!testif HAVE_CXSPARSE
%! n = 20; d= 0.2;
%! a = 1i*sprandn(n,n,d)+speye(n,n);
%! r = qr(a);
%! assert(r'*r,a'*a,1e-10)

%!testif HAVE_CXSPARSE
%! n = 20; d= 0.2;
%! a = 1i*sprandn(n,n,d)+speye(n,n);
%! q = symamd(a);
%! a = a(q,q);
%! r = qr(a);
%! assert(r'*r,a'*a,1e-10)

%!testif HAVE_CXSPARSE
%! n = 20; d= 0.2;
%! a = 1i*sprandn(n,n,d)+speye(n,n);
%! [c,r] = qr(a,ones(n,1));
%! assert (r\c,full(a)\ones(n,1),10e-10)

%!testif HAVE_CXSPARSE
%! n = 20; d= 0.2;
%! a = 1i*sprandn(n,n,d)+speye(n,n);
%! b = randn(n,2);
%! [c,r] = qr(a,b);
%! assert (r\c,full(a)\b,10e-10)

%% Test under-determined systems!!
%!#testif HAVE_CXSPARSE
%! n = 20; d= 0.2;
%! a = 1i*sprandn(n,n+1,d)+speye(n,n+1);
%! b = randn(n,2);
%! [c,r] = qr(a,b);
%! assert (r\c,full(a)\b,10e-10)

%!error qr(sprandn(10,10,0.2),ones(10,1));

*/

DEFUN_DLD (qrupdate, args, ,
  "-*- texinfo -*-\n\
@deftypefn {Loadable Function} {[@var{Q1}, @var{R1}] =} qrupdate (@var{Q}, @var{R}, @var{u}, @var{v})\n\
Given a QR@tie{}factorization of a real or complex matrix\n\
@w{@var{A} = @var{Q}*@var{R}}, @var{Q}@tie{}unitary and\n\
@var{R}@tie{}upper trapezoidal, return the QR@tie{}factorization\n\
of @w{@var{A} + @var{u}*@var{v}'}, where @var{u} and @var{v} are\n\
column vectors (rank-1 update).\n\
\n\
If the matrix @var{Q} is not square, the matrix @var{A} is updated by\n\
Q*Q'*u*v' instead of u*v'.\n\
@seealso{qr, qrinsert, qrdelete}\n\
@end deftypefn")
{
  octave_idx_type nargin = args.length ();
  octave_value_list retval;

  if (nargin != 4)
    {
      print_usage ();
      return retval;
    }

  octave_value argq = args(0);
  octave_value argr = args(1);
  octave_value argu = args(2);
  octave_value argv = args(3);

  if (argq.is_matrix_type () && argr.is_matrix_type () 
      && argu.is_matrix_type () && argv.is_matrix_type ())
    {
      octave_idx_type m = argq.rows ();
      octave_idx_type n = argr.columns ();
      octave_idx_type k = argq.columns ();

      if (argr.rows () == k
          && argu.rows () == m && argu.columns () == 1
          && argv.rows () == n && argv.columns () == 1)
        {
          if (argq.is_real_matrix () 
	      && argr.is_real_matrix () 
	      && argu.is_real_matrix () 
	      && argv.is_real_matrix ())
            {
	      // all real case
	      if (argq.is_single_type () 
		  || argr.is_single_type () 
		  || argu.is_single_type () 
		  || argv.is_single_type ())
		{
		  FloatMatrix Q = argq.float_matrix_value ();
		  FloatMatrix R = argr.float_matrix_value ();
		  FloatMatrix u = argu.float_matrix_value ();
		  FloatMatrix v = argv.float_matrix_value ();

		  FloatQR fact (Q, R);
		  fact.update (u, v);

		  retval(1) = fact.R ();
		  retval(0) = fact.Q ();
		}
	      else
		{
		  Matrix Q = argq.matrix_value ();
		  Matrix R = argr.matrix_value ();
		  Matrix u = argu.matrix_value ();
		  Matrix v = argv.matrix_value ();

		  QR fact (Q, R);
		  fact.update (u, v);

		  retval(1) = fact.R ();
		  retval(0) = fact.Q ();
		}
            }
          else
            {
              // complex case
	      if (argq.is_single_type () 
		  || argr.is_single_type () 
		  || argu.is_single_type () 
		  || argv.is_single_type ())
		{
		  FloatComplexMatrix Q = argq.float_complex_matrix_value ();
		  FloatComplexMatrix R = argr.float_complex_matrix_value ();
		  FloatComplexMatrix u = argu.float_complex_matrix_value ();
		  FloatComplexMatrix v = argv.float_complex_matrix_value ();

		  FloatComplexQR fact (Q, R);
		  fact.update (u, v);
              
		  retval(1) = fact.R ();
		  retval(0) = fact.Q ();
		}
	      else
		{
		  ComplexMatrix Q = argq.complex_matrix_value ();
		  ComplexMatrix R = argr.complex_matrix_value ();
		  ComplexMatrix u = argu.complex_matrix_value ();
		  ComplexMatrix v = argv.complex_matrix_value ();

		  ComplexQR fact (Q, R);
		  fact.update (u, v);
              
		  retval(1) = fact.R ();
		  retval(0) = fact.Q ();
		}
            }
        }
      else
	error ("qrupdate: dimensions mismatch");
    }
  else
    print_usage ();

  return retval;
}
/*
%!shared A, u, v, Ac, uc, vc
%! A = [0.091364  0.613038  0.999083;
%!      0.594638  0.425302  0.603537;
%!      0.383594  0.291238  0.085574;
%!      0.265712  0.268003  0.238409;
%!      0.669966  0.743851  0.445057 ];
%!
%! u = [0.85082;  
%!      0.76426;  
%!      0.42883;  
%!      0.53010;  
%!      0.80683 ];
%!
%! v = [0.98810;
%!      0.24295;
%!      0.43167 ];
%!
%! Ac = [0.620405 + 0.956953i  0.480013 + 0.048806i  0.402627 + 0.338171i;
%!      0.589077 + 0.658457i  0.013205 + 0.279323i  0.229284 + 0.721929i;
%!      0.092758 + 0.345687i  0.928679 + 0.241052i  0.764536 + 0.832406i;
%!      0.912098 + 0.721024i  0.049018 + 0.269452i  0.730029 + 0.796517i;
%!      0.112849 + 0.603871i  0.486352 + 0.142337i  0.355646 + 0.151496i ];
%!
%! uc = [0.20351 + 0.05401i;
%!      0.13141 + 0.43708i;
%!      0.29808 + 0.08789i;
%!      0.69821 + 0.38844i;
%!      0.74871 + 0.25821i ];
%!
%! vc = [0.85839 + 0.29468i;
%!      0.20820 + 0.93090i;
%!      0.86184 + 0.34689i ];
%!

%!test
%! [Q,R] = qr(A);
%! [Q,R] = qrupdate(Q,R,u,v);
%! assert(norm(vec(Q'*Q - eye(5)),Inf) < 1e1*eps)
%! assert(norm(vec(triu(R)-R),Inf) == 0)
%! assert(norm(vec(Q*R - A - u*v'),Inf) < norm(A)*1e1*eps)
%! 
%!test
%! [Q,R] = qr(Ac);
%! [Q,R] = qrupdate(Q,R,uc,vc);
%! assert(norm(vec(Q'*Q - eye(5)),Inf) < 1e1*eps)
%! assert(norm(vec(triu(R)-R),Inf) == 0)
%! assert(norm(vec(Q*R - Ac - uc*vc'),Inf) < norm(Ac)*1e1*eps)

%!test
%! [Q,R] = qr(single(A));
%! [Q,R] = qrupdate(Q,R,single(u),single(v));
%! assert(norm(vec(Q'*Q - eye(5,'single')),Inf) < 1e1*eps('single'))
%! assert(norm(vec(triu(R)-R),Inf) == 0)
%! assert(norm(vec(Q*R - single(A) - single(u)*single(v)'),Inf) < norm(single(A))*1e1*eps('single'))
%! 
%!test
%! [Q,R] = qr(single(Ac));
%! [Q,R] = qrupdate(Q,R,single(uc),single(vc));
%! assert(norm(vec(Q'*Q - eye(5,'single')),Inf) < 1e1*eps('single'))
%! assert(norm(vec(triu(R)-R),Inf) == 0)
%! assert(norm(vec(Q*R - single(Ac) - single(uc)*single(vc)'),Inf) < norm(single(Ac))*1e1*eps('single'))
*/

DEFUN_DLD (qrinsert, args, ,
  "-*- texinfo -*-\n\
@deftypefn {Loadable Function} {[@var{Q1}, @var{R1}] =} qrinsert (@var{Q}, @var{R}, @var{j}, @var{x}, @var{orient})\n\
Given a QR@tie{}factorization of a real or complex matrix\n\
@w{@var{A} = @var{Q}*@var{R}}, @var{Q}@tie{}unitary and\n\
@var{R}@tie{}upper trapezoidal, return the QR@tie{}factorization of\n\
@w{[A(:,1:j-1) x A(:,j:n)]}, where @var{u} is a column vector to be\n\
inserted into @var{A} (if @var{orient} is @code{\"col\"}), or the\n\
QR@tie{}factorization of @w{[A(1:j-1,:);x;A(:,j:n)]}, where @var{x}\n\
is a row vector to be inserted into @var{A} (if @var{orient} is\n\
@code{\"row\"}).\n\
\n\
The default value of @var{orient} is @code{\"col\"}.\n\
\n\
If @var{orient} is @code{\"col\"} and the matrix @var{Q} is not square,\n\
then what gets inserted is the projection of @var{u} onto the space\n\
spanned by columns of @var{Q}, i.e. Q*Q'*u.\n\
\n\
If @var{orient} is @code{\"row\"}, @var{Q} must be square.\n\
@seealso{qr, qrupdate, qrdelete}\n\
@end deftypefn")
{
  octave_idx_type nargin = args.length ();
  octave_value_list retval;

  if (nargin < 4 || nargin > 5)
    {
      print_usage ();
      return retval;
    }
  
  octave_value argq = args(0);
  octave_value argr = args(1);
  octave_value argj = args(2);
  octave_value argx = args(3);
      
  if (argq.is_matrix_type () && argr.is_matrix_type ()
      && argj.is_scalar_type () && argx.is_matrix_type ()
      && (nargin < 5 || args(4).is_string ()))
    {
      octave_idx_type m = argq.rows ();
      octave_idx_type n = argr.columns ();
      octave_idx_type k = argq.columns ();

      std::string orient = (nargin < 5) ? "col" : args(4).string_value ();

      bool row = orient == "row";

      if (row || orient == "col")
        if (argr.rows () == k 
            && (! row || m == k)
            && argx.rows () == (row ? 1 : m)
            && argx.columns () == (row ? n : 1))
          {
            octave_idx_type j = argj.idx_type_value ();

            if (j >= 1 && j <= (row ? n : m)+1)
              {
                if (argq.is_real_matrix () 
		    && argr.is_real_matrix () 
		    && argx.is_real_matrix ())
                  {
                    // real case
		    if (argq.is_single_type () 
			|| argr.is_single_type () 
			|| argx.is_single_type ())
		      {
			FloatMatrix Q = argq.float_matrix_value ();
			FloatMatrix R = argr.float_matrix_value ();
			FloatMatrix x = argx.float_matrix_value ();

			FloatQR fact (Q, R);

			if (row) 
			  fact.insert_row (x, j-1);
			else 
			  fact.insert_col (x, j-1);

			retval(1) = fact.R ();
			retval(0) = fact.Q ();

		      }
		    else
		      {
			Matrix Q = argq.matrix_value ();
			Matrix R = argr.matrix_value ();
			Matrix x = argx.matrix_value ();

			QR fact (Q, R);

			if (row) 
			  fact.insert_row (x, j-1);
			else 
			  fact.insert_col (x, j-1);

			retval(1) = fact.R ();
			retval(0) = fact.Q ();

		      }
                  }
                else
                  {
                    // complex case
		    if (argq.is_single_type () 
			|| argr.is_single_type () 
			|| argx.is_single_type ())
		      {
			FloatComplexMatrix Q = argq.float_complex_matrix_value ();
			FloatComplexMatrix R = argr.float_complex_matrix_value ();
			FloatComplexMatrix x = argx.float_complex_matrix_value ();

			FloatComplexQR fact (Q, R);

			if (row) 
			  fact.insert_row (x, j-1);
			else 
			  fact.insert_col (x, j-1);

			retval(1) = fact.R ();
			retval(0) = fact.Q ();
		      }
		    else
		      {
			ComplexMatrix Q = argq.complex_matrix_value ();
			ComplexMatrix R = argr.complex_matrix_value ();
			ComplexMatrix x = argx.complex_matrix_value ();

			ComplexQR fact (Q, R);

			if (row) 
			  fact.insert_row (x, j-1);
			else 
			  fact.insert_col (x, j-1);

			retval(1) = fact.R ();
			retval(0) = fact.Q ();
		      }
                  }

              }
            else
              error ("qrinsert: index j out of range");
          }
        else
          error ("qrinsert: dimension mismatch");

      else
        error ("qrinsert: orient must be \"col\" or \"row\"");
    }
  else
    print_usage ();

  return retval;
}

/*
%!test
%! [Q,R] = qr(A);
%! [Q,R] = qrinsert(Q,R,3,u);
%! assert(norm(vec(Q'*Q - eye(5)),Inf) < 1e1*eps)
%! assert(norm(vec(triu(R)-R),Inf) == 0)
%! assert(norm(vec(Q*R - [A(:,1:2) u A(:,3)]),Inf) < norm(A)*1e1*eps)
%!test
%! [Q,R] = qr(Ac);
%! [Q,R] = qrinsert(Q,R,3,uc);
%! assert(norm(vec(Q'*Q - eye(5)),Inf) < 1e1*eps)
%! assert(norm(vec(triu(R)-R),Inf) == 0)
%! assert(norm(vec(Q*R - [Ac(:,1:2) uc Ac(:,3)]),Inf) < norm(Ac)*1e1*eps)
%!test
%! x = [0.85082  0.76426  0.42883 ];
%!
%! [Q,R] = qr(A);
%! [Q,R] = qrinsert(Q,R,3,x,'row');
%! assert(norm(vec(Q'*Q - eye(6)),Inf) < 1e1*eps)
%! assert(norm(vec(triu(R)-R),Inf) == 0)
%! assert(norm(vec(Q*R - [A(1:2,:);x;A(3:5,:)]),Inf) < norm(A)*1e1*eps)
%!test
%! x = [0.20351 + 0.05401i  0.13141 + 0.43708i  0.29808 + 0.08789i ];
%!
%! [Q,R] = qr(Ac);
%! [Q,R] = qrinsert(Q,R,3,x,'row');
%! assert(norm(vec(Q'*Q - eye(6)),Inf) < 1e1*eps)
%! assert(norm(vec(triu(R)-R),Inf) == 0)
%! assert(norm(vec(Q*R - [Ac(1:2,:);x;Ac(3:5,:)]),Inf) < norm(Ac)*1e1*eps)

%!test
%! [Q,R] = qr(single(A));
%! [Q,R] = qrinsert(Q,R,3,single(u));
%! assert(norm(vec(Q'*Q - eye(5,'single')),Inf) < 1e1*eps('single'))
%! assert(norm(vec(triu(R)-R),Inf) == 0)
%! assert(norm(vec(Q*R - single([A(:,1:2) u A(:,3)])),Inf) < norm(single(A))*1e1*eps('single'))
%!test
%! [Q,R] = qr(single(Ac));
%! [Q,R] = qrinsert(Q,R,3,single(uc));
%! assert(norm(vec(Q'*Q - eye(5,'single')),Inf) < 1e1*eps('single'))
%! assert(norm(vec(triu(R)-R),Inf) == 0)
%! assert(norm(vec(Q*R - single([Ac(:,1:2) uc Ac(:,3)])),Inf) < norm(single(Ac))*1e1*eps('single'))
%!test
%! x = single([0.85082  0.76426  0.42883 ]);
%!
%! [Q,R] = qr(single(A));
%! [Q,R] = qrinsert(Q,R,3,x,'row');
%! assert(norm(vec(Q'*Q - eye(6,'single')),Inf) < 1e1*eps('single'))
%! assert(norm(vec(triu(R)-R),Inf) == 0)
%! assert(norm(vec(Q*R - single([A(1:2,:);x;A(3:5,:)])),Inf) < norm(single(A))*1e1*eps('single'))
%!test
%! x = single([0.20351 + 0.05401i  0.13141 + 0.43708i  0.29808 + 0.08789i ]);
%!
%! [Q,R] = qr(single(Ac));
%! [Q,R] = qrinsert(Q,R,3,x,'row');
%! assert(norm(vec(Q'*Q - eye(6,'single')),Inf) < 1e1*eps('single'))
%! assert(norm(vec(triu(R)-R),Inf) == 0)
%! assert(norm(vec(Q*R - single([Ac(1:2,:);x;Ac(3:5,:)])),Inf) < norm(single(Ac))*1e1*eps('single'))
*/

DEFUN_DLD (qrdelete, args, ,
  "-*- texinfo -*-\n\
@deftypefn {Loadable Function} {[@var{Q1}, @var{R1}] =} qrdelete (@var{Q}, @var{R}, @var{j}, @var{orient})\n\
Given a QR@tie{}factorization of a real or complex matrix\n\
@w{@var{A} = @var{Q}*@var{R}}, @var{Q}@tie{}unitary and\n\
@var{R}@tie{}upper trapezoidal, return the QR@tie{}factorization of\n\
@w{[A(:,1:j-1) A(:,j+1:n)]}, i.e. @var{A} with one column deleted\n\
(if @var{orient} is \"col\"), or the QR@tie{}factorization of\n\
@w{[A(1:j-1,:);A(:,j+1:n)]}, i.e. @var{A} with one row deleted (if\n\
@var{orient} is \"row\").\n\
\n\
The default value of @var{orient} is \"col\".\n\
\n\
If @var{orient} is \"col\", the matrix @var{Q} is not required to\n\
be square.\n\
\n\
For @sc{Matlab} compatibility, if @var{Q} is nonsquare on input, the\n\
updated factorization is always stripped to the economical form, i.e.\n\
@code{columns (Q) == rows (R) <= columns (R)}.\n\
\n\
To get the less intelligent but more natural behaviour when @var{Q}\n\
retains it shape and @var{R} loses one column, set @var{orient} to\n\
\"col+\" instead.\n\
\n\
If @var{orient} is \"row\", @var{Q} must be square.\n\
@seealso{qr, qrinsert, qrupdate}\n\
@end deftypefn")
{
  octave_idx_type nargin = args.length ();
  octave_value_list retval;

  if (nargin < 3 || nargin > 4)
    {
      print_usage ();
      return retval;
    }

  octave_value argq = args(0);
  octave_value argr = args(1);
  octave_value argj = args(2);

  if (argq.is_matrix_type () && argr.is_matrix_type () && argj.is_scalar_type ()
      && (nargin < 4 || args(3).is_string ()))
    {
      octave_idx_type m = argq.rows ();
      octave_idx_type k = argq.columns ();
      octave_idx_type n = argr.columns ();

      std::string orient = (nargin < 4) ? "col" : args(3).string_value ();

      bool row = orient == "row";
      bool colp = orient == "col+";

      if (row || colp || orient == "col")
        if (argr.rows () == k
            && (! row || m == k))
          {
            octave_idx_type j = argj.scalar_value ();
            if (j >= 1 && j <= (row ? n : m))
              {
                if (argq.is_real_matrix ()
		    && argr.is_real_matrix ())
                  {
                    // real case
		    if (argq.is_single_type ()
			|| argr.is_single_type ())
		      {
			FloatMatrix Q = argq.float_matrix_value ();
			FloatMatrix R = argr.float_matrix_value ();

			FloatQR fact (Q, R);

			if (row) 
			  fact.delete_row (j-1);
			else 
			  {
			    fact.delete_col (j-1);

			    if (! colp && k < m)
			      fact.economize ();
			  }

			retval(1) = fact.R ();
			retval(0) = fact.Q ();
		      }
		    else
		      {
			Matrix Q = argq.matrix_value ();
			Matrix R = argr.matrix_value ();

			QR fact (Q, R);

			if (row) 
			  fact.delete_row (j-1);
			else 
			  {
			    fact.delete_col (j-1);

			    if (! colp && k < m)
			      fact.economize ();
			  }

			retval(1) = fact.R ();
			retval(0) = fact.Q ();
		      }
                  }
                else
                  {
                    // complex case
		    if (argq.is_single_type ()
			|| argr.is_single_type ())
		      {
			FloatComplexMatrix Q = argq.float_complex_matrix_value ();
			FloatComplexMatrix R = argr.float_complex_matrix_value ();

			FloatComplexQR fact (Q, R);

			if (row) 
			  fact.delete_row (j-1);
			else 
			  {
			    fact.delete_col (j-1);

			    if (! colp && k < m)
			      fact.economize ();
			  }

			retval(1) = fact.R ();
			retval(0) = fact.Q ();
		      }
		    else
		      {
			ComplexMatrix Q = argq.complex_matrix_value ();
			ComplexMatrix R = argr.complex_matrix_value ();

			ComplexQR fact (Q, R);

			if (row) 
			  fact.delete_row (j-1);
			else 
			  {
			    fact.delete_col (j-1);

			    if (! colp && k < m)
			      fact.economize ();
			  }

			retval(1) = fact.R ();
			retval(0) = fact.Q ();
		      }
                  }
              }
            else
              error ("qrdelete: index j out of range");
          }
        else
          error ("qrdelete: dimension mismatch");

      else
        error ("qrdelete: orient must be \"col\" or \"row\"");
    }
  else
    print_usage ();

  return retval;
}
 
/*
%!test
%! AA = [0.091364  0.613038  0.027504  0.999083;
%!       0.594638  0.425302  0.562834  0.603537;
%!       0.383594  0.291238  0.742073  0.085574;
%!       0.265712  0.268003  0.783553  0.238409;
%!       0.669966  0.743851  0.457255  0.445057 ];
%!
%! [Q,R] = qr(AA);
%! [Q,R] = qrdelete(Q,R,3);
%! assert(norm(vec(Q'*Q - eye(5)),Inf) < 1e1*eps)
%! assert(norm(vec(triu(R)-R),Inf) == 0)
%! assert(norm(vec(Q*R - [AA(:,1:2) AA(:,4)]),Inf) < norm(AA)*1e1*eps)
%! 
%!test
%! AA = [0.364554 + 0.993117i  0.669818 + 0.510234i  0.426568 + 0.041337i  0.847051 + 0.233291i;
%!       0.049600 + 0.242783i  0.448946 + 0.484022i  0.141155 + 0.074420i  0.446746 + 0.392706i;
%!       0.581922 + 0.657416i  0.581460 + 0.030016i  0.219909 + 0.447288i  0.201144 + 0.069132i;
%!       0.694986 + 0.000571i  0.682327 + 0.841712i  0.807537 + 0.166086i  0.192767 + 0.358098i;
%!       0.945002 + 0.066788i  0.350492 + 0.642638i  0.579629 + 0.048102i  0.600170 + 0.636938i ] * I;
%!
%! [Q,R] = qr(AA);
%! [Q,R] = qrdelete(Q,R,3);
%! assert(norm(vec(Q'*Q - eye(5)),Inf) < 1e1*eps)
%! assert(norm(vec(triu(R)-R),Inf) == 0)
%! assert(norm(vec(Q*R - [AA(:,1:2) AA(:,4)]),Inf) < norm(AA)*1e1*eps)
%!
%!test
%! AA = [0.091364  0.613038  0.027504  0.999083;
%!       0.594638  0.425302  0.562834  0.603537;
%!       0.383594  0.291238  0.742073  0.085574;
%!       0.265712  0.268003  0.783553  0.238409;
%!       0.669966  0.743851  0.457255  0.445057 ];
%!
%! [Q,R] = qr(AA);
%! [Q,R] = qrdelete(Q,R,3,'row');
%! assert(norm(vec(Q'*Q - eye(4)),Inf) < 1e1*eps)
%! assert(norm(vec(triu(R)-R),Inf) == 0)
%! assert(norm(vec(Q*R - [AA(1:2,:);AA(4:5,:)]),Inf) < norm(AA)*1e1*eps)
%! 
%!test
%! AA = [0.364554 + 0.993117i  0.669818 + 0.510234i  0.426568 + 0.041337i  0.847051 + 0.233291i;
%!       0.049600 + 0.242783i  0.448946 + 0.484022i  0.141155 + 0.074420i  0.446746 + 0.392706i;
%!       0.581922 + 0.657416i  0.581460 + 0.030016i  0.219909 + 0.447288i  0.201144 + 0.069132i;
%!       0.694986 + 0.000571i  0.682327 + 0.841712i  0.807537 + 0.166086i  0.192767 + 0.358098i;
%!       0.945002 + 0.066788i  0.350492 + 0.642638i  0.579629 + 0.048102i  0.600170 + 0.636938i ] * I;
%!
%! [Q,R] = qr(AA);
%! [Q,R] = qrdelete(Q,R,3,'row');
%! assert(norm(vec(Q'*Q - eye(4)),Inf) < 1e1*eps)
%! assert(norm(vec(triu(R)-R),Inf) == 0)
%! assert(norm(vec(Q*R - [AA(1:2,:);AA(4:5,:)]),Inf) < norm(AA)*1e1*eps)

%!test
%! AA = single([0.091364  0.613038  0.027504  0.999083;
%!              0.594638  0.425302  0.562834  0.603537;
%!              0.383594  0.291238  0.742073  0.085574;
%!              0.265712  0.268003  0.783553  0.238409;
%!              0.669966  0.743851  0.457255  0.445057 ]);
%!
%! [Q,R] = qr(AA);
%! [Q,R] = qrdelete(Q,R,3);
%! assert(norm(vec(Q'*Q - eye(5,'single')),Inf) < 1e1*eps('single'))
%! assert(norm(vec(triu(R)-R),Inf) == 0)
%! assert(norm(vec(Q*R - [AA(:,1:2) AA(:,4)]),Inf) < norm(AA)*1e1*eps('single'))
%! 
%!test
%! AA = single([0.364554 + 0.993117i  0.669818 + 0.510234i  0.426568 + 0.041337i  0.847051 + 0.233291i;
%!              0.049600 + 0.242783i  0.448946 + 0.484022i  0.141155 + 0.074420i  0.446746 + 0.392706i;
%!              0.581922 + 0.657416i  0.581460 + 0.030016i  0.219909 + 0.447288i  0.201144 + 0.069132i;
%!              0.694986 + 0.000571i  0.682327 + 0.841712i  0.807537 + 0.166086i  0.192767 + 0.358098i;
%!              0.945002 + 0.066788i  0.350492 + 0.642638i  0.579629 + 0.048102i  0.600170 + 0.636938i ]) * I;
%!
%! [Q,R] = qr(AA);
%! [Q,R] = qrdelete(Q,R,3);
%! assert(norm(vec(Q'*Q - eye(5,'single')),Inf) < 1e1*eps('single'))
%! assert(norm(vec(triu(R)-R),Inf) == 0)
%! assert(norm(vec(Q*R - [AA(:,1:2) AA(:,4)]),Inf) < norm(AA)*1e1*eps('single'))
%!
%!test
%! AA = single([0.091364  0.613038  0.027504  0.999083;
%!              0.594638  0.425302  0.562834  0.603537;
%!              0.383594  0.291238  0.742073  0.085574;
%!              0.265712  0.268003  0.783553  0.238409;
%!              0.669966  0.743851  0.457255  0.445057 ]);
%!
%! [Q,R] = qr(AA);
%! [Q,R] = qrdelete(Q,R,3,'row');
%! assert(norm(vec(Q'*Q - eye(4,'single')),Inf) < 1e1*eps('single'))
%! assert(norm(vec(triu(R)-R),Inf) == 0)
%! assert(norm(vec(Q*R - [AA(1:2,:);AA(4:5,:)]),Inf) < norm(AA)*1e1*eps('single'))
%! 
%!test
%! AA = single([0.364554 + 0.993117i  0.669818 + 0.510234i  0.426568 + 0.041337i  0.847051 + 0.233291i;
%!              0.049600 + 0.242783i  0.448946 + 0.484022i  0.141155 + 0.074420i  0.446746 + 0.392706i;
%!              0.581922 + 0.657416i  0.581460 + 0.030016i  0.219909 + 0.447288i  0.201144 + 0.069132i;
%!              0.694986 + 0.000571i  0.682327 + 0.841712i  0.807537 + 0.166086i  0.192767 + 0.358098i;
%!              0.945002 + 0.066788i  0.350492 + 0.642638i  0.579629 + 0.048102i  0.600170 + 0.636938i ]) * I;
%!
%! [Q,R] = qr(AA);
%! [Q,R] = qrdelete(Q,R,3,'row');
%! assert(norm(vec(Q'*Q - eye(4,'single')),Inf) < 1e1*eps('single'))
%! assert(norm(vec(triu(R)-R),Inf) == 0)
%! assert(norm(vec(Q*R - [AA(1:2,:);AA(4:5,:)]),Inf) < norm(AA)*1e1*eps('single'))
*/

DEFUN_DLD (qrshift, args, ,
  "-*- texinfo -*-\n\
@deftypefn {Loadable Function} {[@var{Q1}, @var{R1}] =} qrshift (@var{Q}, @var{R}, @var{i}, @var{j})\n\
Given a QR@tie{}factorization of a real or complex matrix\n\
@w{@var{A} = @var{Q}*@var{R}}, @var{Q}@tie{}unitary and\n\
@var{R}@tie{}upper trapezoidal, return the QR@tie{}factorization\n\
of @w{@var{A}(:,p)}, where @w{p} is the permutation @*\n\
@code{p = [1:i-1, shift(i:j, 1), j+1:n]} if @w{@var{i} < @var{j}} @*\n\
 or @*\n\
@code{p = [1:j-1, shift(j:i,-1), i+1:n]} if @w{@var{j} < @var{i}}. @*\n\
\n\
@seealso{qr, qrinsert, qrdelete}\n\
@end deftypefn")
{
  octave_idx_type nargin = args.length ();
  octave_value_list retval;

  if (nargin != 4)
    {
      print_usage ();
      return retval;
    }

  octave_value argq = args(0);
  octave_value argr = args(1);
  octave_value argi = args(2);
  octave_value argj = args(3);

  if (argq.is_matrix_type () && argr.is_matrix_type () 
      && argi.is_real_scalar () && argj.is_real_scalar ())
    {
      octave_idx_type n = argr.columns ();
      octave_idx_type k = argq.columns ();

      if (argr.rows () == k)
        {
          octave_idx_type i = argi.scalar_value ();
          octave_idx_type j = argj.scalar_value ();
          if (i > 1 && i <= n && j > 1 && j <= n)
            {
              if (argq.is_real_matrix () 
                  && argr.is_real_matrix ())
                {
                  // all real case
		  if (argq.is_single_type () 
		      && argr.is_single_type ())
		    {
		      FloatMatrix Q = argq.float_matrix_value ();
		      FloatMatrix R = argr.float_matrix_value ();

		      FloatQR fact (Q, R);
		      fact.shift_cols (i-1, j-1);

		      retval(1) = fact.R ();
		      retval(0) = fact.Q ();
		    }
		  else
		    {
		      Matrix Q = argq.matrix_value ();
		      Matrix R = argr.matrix_value ();

		      QR fact (Q, R);
		      fact.shift_cols (i-1, j-1);

		      retval(1) = fact.R ();
		      retval(0) = fact.Q ();
		    }
                }
              else
                {
                  // complex case
		  if (argq.is_single_type () 
		      && argr.is_single_type ())
		    {
		      FloatComplexMatrix Q = argq.float_complex_matrix_value ();
		      FloatComplexMatrix R = argr.float_complex_matrix_value ();

		      FloatComplexQR fact (Q, R);
		      fact.shift_cols (i-1, j-1);
                  
		      retval(1) = fact.R ();
		      retval(0) = fact.Q ();
		    }
		  else
		    {
		      ComplexMatrix Q = argq.complex_matrix_value ();
		      ComplexMatrix R = argr.complex_matrix_value ();

		      ComplexQR fact (Q, R);
		      fact.shift_cols (i-1, j-1);
                  
		      retval(1) = fact.R ();
		      retval(0) = fact.Q ();
		    }
                }
            }
          else
            error ("qrshift: index out of range");
        }
      else
	error ("qrshift: dimensions mismatch");
    }
  else
    print_usage ();

  return retval;
}
/*
%!test
%! AA = A.';
%! i = 2; j = 4; p = [1:i-1, shift(i:j,-1), j+1:5];
%!
%! [Q,R] = qr(AA);
%! [Q,R] = qrshift(Q,R,i,j);
%! assert(norm(vec(Q'*Q - eye(3)),Inf) < 1e1*eps)
%! assert(norm(vec(triu(R)-R),Inf) == 0)
%! assert(norm(vec(Q*R - AA(:,p)),Inf) < norm(AA)*1e1*eps)
%! 
%! j = 2; i = 4; p = [1:j-1, shift(j:i,+1), i+1:5];
%!
%! [Q,R] = qr(AA);
%! [Q,R] = qrshift(Q,R,i,j);
%! assert(norm(vec(Q'*Q - eye(3)),Inf) < 1e1*eps)
%! assert(norm(vec(triu(R)-R),Inf) == 0)
%! assert(norm(vec(Q*R - AA(:,p)),Inf) < norm(AA)*1e1*eps)
%! 
%!test
%! AA = Ac.';
%! i = 2; j = 4; p = [1:i-1, shift(i:j,-1), j+1:5];
%!
%! [Q,R] = qr(AA);
%! [Q,R] = qrshift(Q,R,i,j);
%! assert(norm(vec(Q'*Q - eye(3)),Inf) < 1e1*eps)
%! assert(norm(vec(triu(R)-R),Inf) == 0)
%! assert(norm(vec(Q*R - AA(:,p)),Inf) < norm(AA)*1e1*eps)
%! 
%! j = 2; i = 4; p = [1:j-1, shift(j:i,+1), i+1:5];
%!
%! [Q,R] = qr(AA);
%! [Q,R] = qrshift(Q,R,i,j);
%! assert(norm(vec(Q'*Q - eye(3)),Inf) < 1e1*eps)
%! assert(norm(vec(triu(R)-R),Inf) == 0)
%! assert(norm(vec(Q*R - AA(:,p)),Inf) < norm(AA)*1e1*eps)


%!test
%! AA = single (A).';
%! i = 2; j = 4; p = [1:i-1, shift(i:j,-1), j+1:5];
%!
%! [Q,R] = qr(AA);
%! [Q,R] = qrshift(Q,R,i,j);
%! assert(norm(vec(Q'*Q - eye(3,'single')),Inf) < 1e1*eps('single'))
%! assert(norm(vec(triu(R)-R),Inf) == 0)
%! assert(norm(vec(Q*R - AA(:,p)),Inf) < norm(AA)*1e1*eps('single'))
%! 
%! j = 2; i = 4; p = [1:j-1, shift(j:i,+1), i+1:5];
%!
%! [Q,R] = qr(AA);
%! [Q,R] = qrshift(Q,R,i,j);
%! assert(norm(vec(Q'*Q - eye(3,'single')),Inf) < 1e1*eps('single'))
%! assert(norm(vec(triu(R)-R),Inf) == 0)
%! assert(norm(vec(Q*R - AA(:,p)),Inf) < norm(AA)*1e1*eps('single'))
%! 
%!test
%! AA = single(Ac).';
%! i = 2; j = 4; p = [1:i-1, shift(i:j,-1), j+1:5];
%!
%! [Q,R] = qr(AA);
%! [Q,R] = qrshift(Q,R,i,j);
%! assert(norm(vec(Q'*Q - eye(3,'single')),Inf) < 1e1*eps('single'))
%! assert(norm(vec(triu(R)-R),Inf) == 0)
%! assert(norm(vec(Q*R - AA(:,p)),Inf) < norm(AA)*1e1*eps('single'))
%! 
%! j = 2; i = 4; p = [1:j-1, shift(j:i,+1), i+1:5];
%!
%! [Q,R] = qr(AA);
%! [Q,R] = qrshift(Q,R,i,j);
%! assert(norm(vec(Q'*Q - eye(3,'single')),Inf) < 1e1*eps('single'))
%! assert(norm(vec(triu(R)-R),Inf) == 0)
%! assert(norm(vec(Q*R - AA(:,p)),Inf) < norm(AA)*1e1*eps('single'))
*/

/*
;;; Local Variables: ***
;;; mode: C++ ***
;;; End: ***
*/