/*

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, and qrdelete 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 "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_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;
}

/*

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
              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
              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;
}
/*
%!test
%! 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 ];
%!
%! [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
%! A = [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 ];
%!
%! u = [0.20351 + 0.05401i;
%!      0.13141 + 0.43708i;
%!      0.29808 + 0.08789i;
%!      0.69821 + 0.38844i;
%!      0.74871 + 0.25821i ];
%!
%! v = [0.85839 + 0.29468i;
%!      0.20820 + 0.93090i;
%!      0.86184 + 0.34689i ];
%!
%! [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)
*/

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.int_value ();

            if (j >= 1 && j <= (row ? n : m)+1)
              {
                if (argq.is_real_matrix () 
		    && argr.is_real_matrix () 
		    && argx.is_real_matrix ())
                  {
                    // real case
                    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);
                    else 
                      fact.insert_col (x, j);

                    retval(1) = fact.R ();
                    retval(0) = fact.Q ();
                  }
                else
                  {
                    // complex case
                    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);
                    else 
                      fact.insert_col (x, j);

                    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
%! 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 ];
%!
%! x = [0.85082;  
%!      0.76426;  
%!      0.42883;  
%!      0.53010;  
%!      0.80683 ];
%!
%! [Q,R] = qr(A);
%! [Q,R] = qrinsert(Q,R,3,x);
%! 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) x A(:,3)]),Inf) < norm(A)*1e1*eps)
%! 
%!test
%! A = [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 ];
%!
%! x = [0.20351 + 0.05401i;
%!      0.13141 + 0.43708i;
%!      0.29808 + 0.08789i;
%!      0.69821 + 0.38844i;
%!      0.74871 + 0.25821i ];
%!
%! [Q,R] = qr(A);
%! [Q,R] = qrinsert(Q,R,3,x);
%! 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) x A(:,3)]),Inf) < norm(A)*1e1*eps)
%!
%!test
%! 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 ];
%!
%! 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
%! A = [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 ];
%!
%! x = [0.20351 + 0.05401i  0.13141 + 0.43708i  0.29808 + 0.08789i ];
%!
%! [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)
*/

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
                    Matrix Q = argq.matrix_value ();
                    Matrix R = argr.matrix_value ();

                    QR fact (Q, R);

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

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

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

                    ComplexQR fact (Q, R);

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

                        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
%! A = [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(A);
%! [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 - [A(:,1:2) A(:,4)]),Inf) < norm(A)*1e1*eps)
%! 
%!test
%! A = [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(A);
%! [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 - [A(:,1:2) A(:,4)]),Inf) < norm(A)*1e1*eps)
%!
%!test
%! A = [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(A);
%! [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 - [A(1:2,:);A(4:5,:)]),Inf) < norm(A)*1e1*eps)
%! 
%!test
%! A = [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(A);
%! [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 - [A(1:2,:);A(4:5,:)]),Inf) < norm(A)*1e1*eps)
*/

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