/*

Copyright (C) 1996, 1997, 1999, 2000, 2005, 2006, 2007 John W. Eaton
Copyright (C) 2008, 2009 Jaroslav Hajek
Copyright (C) 2008, 2009 VZLU Prague

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

*/

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

template <class MT>
static octave_value
maybe_set_triangular (const MT& m, MatrixType::matrix_type t = MatrixType::Upper)
{
  typedef typename MT::element_type T;
  octave_value retval;
  octave_idx_type r = m.rows (), c = m.columns ();
  if (r == c)
    {
      const T zero = T();
      octave_idx_type i = 0;
      for (;i != r && m(i,i) != zero; i++) ;
      if (i == r)
        retval = octave_value (m, MatrixType (t));
      else
        retval = m;
    }
  else
    retval = m;

  return retval;
}

// [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\
@ifnottex\n\
\n\
@example\n\
@code{min norm(A x - b)}\n\
@end example\n\
\n\
@end ifnottex\n\
for overdetermined systems of equations (i.e.,\n\
@iftex\n\
@tex\n\
$A$\n\
@end tex\n\
@end iftex\n\
@ifnottex\n\
@code{a}\n\
@end ifnottex\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\
@ifnottex\n\
@code{q * r = a} where @code{q} is an orthogonal matrix and @code{r} is\n\
upper triangular.\n\
@end ifnottex\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;

  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) = maybe_set_triangular (fact.R ());
			retval(0) = fact.Q ();
		      }
		      break;

		    default:
		      {
			FloatQRP fact (m, type);
                        if (type == QR::economy)
                          retval(2) = fact.Pvec ();
                        else
                          retval(2) = fact.P ();
			retval(1) = maybe_set_triangular (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) = maybe_set_triangular (fact.R ());
			retval(0) = fact.Q ();
		      }
		      break;

		    default:
		      {
			FloatComplexQRP fact (m, type);
                        if (type == QR::economy)
                          retval(2) = fact.Pvec ();
                        else
                          retval(2) = fact.P ();
			retval(1) = maybe_set_triangular (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) = maybe_set_triangular (fact.R ());
			retval(0) = fact.Q ();
		      }
		      break;

		    default:
		      {
			QRP fact (m, type);
                        if (type == QR::economy)
                          retval(2) = fact.Pvec ();
                        else
                          retval(2) = fact.P ();
			retval(1) = maybe_set_triangular (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) = maybe_set_triangular (fact.R ());
			retval(0) = fact.Q ();
		      }
		      break;

		    default:
		      {
			ComplexQRP fact (m, type);
                        if (type == QR::economy)
                          retval(2) = fact.Pvec ();
                        else
                          retval(2) = fact.P ();
			retval(1) = maybe_set_triangular (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));

*/

static
bool check_qr_dims (const octave_value& q, const octave_value& r,
                    bool allow_ecf = false)
{
  octave_idx_type m = q.rows (), k = r.rows (), n = r.columns ();
  return ((q.ndims () == 2 && r.ndims () == 2 && k == q.columns ())
            && (m == k || (allow_ecf && k == n && k < m)));
}

static 
bool check_index (const octave_value& i, bool vector_allowed = false)
{
  return ((i.is_real_type () || i.is_integer_type ()) 
          && (i.is_scalar_type () || vector_allowed));
}

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) or matrices with equal number of columns\n\
(rank-k update). Notice that the latter case is done as a sequence of rank-1 updates;\n\
thus, for k large enough, it will be both faster and more accurate to recompute\n\
the factorization from scratch.\n\
\n\
The QR factorization supplied may be either full\n\
(Q is square) or economized (R is square).\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 argu = args(2);
  octave_value argv = args(3);

  if (argq.is_numeric_type () && argr.is_numeric_type () 
      && argu.is_numeric_type () && argv.is_numeric_type ())
    {
      if (check_qr_dims (argq, argr, true))
        {
          if (argq.is_real_type () 
	      && argr.is_real_type () 
	      && argu.is_real_type () 
	      && argv.is_real_type ())
            {
	      // 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
    error ("qrupdate: expecting numeric arguments");

  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\
If @var{orient} is @code{\"col\"},\n\
@var{u} may be a matrix and @var{j} an index vector\n\
resulting in the QR@tie{}factorization of a matrix @var{B} such that\n\
@w{B(:,@var{j})} gives @var{u} and @w{B(:,@var{j}) = []} gives @var{A}.\n\
Notice that the latter case is done as a sequence of k insertions;\n\
thus, for k large enough, it will be both faster and more accurate to recompute\n\
the factorization from scratch.\n\
\n\
If @var{orient} is @code{\"col\"},\n\
the QR factorization supplied may be either full\n\
(Q is square) or economized (R is square).\n\
\n\
If @var{orient} is @code{\"row\"}, full factorization is needed.\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_numeric_type () && argr.is_numeric_type ()
      && argx.is_numeric_type ()
      && (nargin < 5 || args(4).is_string ()))
    {
      std::string orient = (nargin < 5) ? "col" : args(4).string_value ();

      bool col = orient == "col";

      if (col || orient == "row")
        if (check_qr_dims (argq, argr, col) 
            && (col || argx.rows () == 1))
          {
            if (check_index (argj, col))
              {
                MArray<octave_idx_type> j
		  = argj.octave_idx_type_vector_value ();

                if (argq.is_real_type () 
		    && argr.is_real_type () 
		    && argx.is_real_type ())
                  {
                    // 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 (col) 
			  fact.insert_col (x, j-1);
			else 
			  fact.insert_row (x.row (0), j(0)-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 (col) 
			  fact.insert_col (x, j-1);
			else 
			  fact.insert_row (x.row (0), j(0)-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 (col) 
			  fact.insert_col (x, j-1);
			else 
			  fact.insert_row (x.row (0), j(0)-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 (col) 
			  fact.insert_col (x, j-1);
			else 
			  fact.insert_row (x.row (0), j(0)-1);

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

              }
            else
              error ("qrinsert: invalid index");
          }
        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 @code{\"col\"},\n\
@var{j} may be an index vector\n\
resulting in the QR@tie{}factorization of a matrix @var{B} such that\n\
@w{A(:,@var{j}) = []} gives @var{B}.\n\
Notice that the latter case is done as a sequence of k deletions;\n\
thus, for k large enough, it will be both faster and more accurate to recompute\n\
the factorization from scratch.\n\
\n\
If @var{orient} is @code{\"col\"},\n\
the QR factorization supplied may be either full\n\
(Q is square) or economized (R is square).\n\
\n\
If @var{orient} is @code{\"row\"}, full factorization is needed.\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_numeric_type () && argr.is_numeric_type ()
      && (nargin < 4 || args(3).is_string ()))
    {
      std::string orient = (nargin < 4) ? "col" : args(3).string_value ();

      bool col = orient == "col";

      if (col || orient == "row")
        if (check_qr_dims (argq, argr, col))
          {
            if (check_index (argj, col))
              {
                MArray<octave_idx_type> j
		  = argj.octave_idx_type_vector_value ();

                if (argq.is_real_type ()
		    && argr.is_real_type ())
                  {
                    // 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 (col) 
                          fact.delete_col (j-1);
			else 
			  fact.delete_row (j(0)-1);

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

			QR fact (Q, R);

			if (col) 
                          fact.delete_col (j-1);
			else 
			  fact.delete_row (j(0)-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);

			if (col) 
                          fact.delete_col (j-1);
			else 
			  fact.delete_row (j(0)-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);

			if (col) 
                          fact.delete_col (j-1);
			else 
			  fact.delete_row (j(0)-1);

			retval(1) = fact.R ();
			retval(0) = fact.Q ();
		      }
                  }
              }
            else
              error ("qrdelete: invalid index");
          }
        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_numeric_type () && argr.is_numeric_type ())
    {
      if (check_qr_dims (argq, argr, true))
        {
          if (check_index (argi) && check_index (argj))
            {
              octave_idx_type i = argi.int_value ();
              octave_idx_type j = argj.int_value ();

              if (argq.is_real_type () 
                  && argr.is_real_type ())
                {
                  // 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: invalid index");
        }
      else
	error ("qrshift: dimensions mismatch");
    }
  else
    error ("qrshift: expecting numeric arguments");

  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: ***
*/