# HG changeset patch
# User John W. Eaton <jwe@octave.org>
# Date 1455695841 18000
# Node ID eb1524b07fe367bc26305df42e30ec4ad5d851a2
# Parent  48d82f74243ec89b7d7bf24d862c16650e3a2a44
better use of templates for qr classes

* liboctave/numeric/qr.h, liboctave/numeric/qr.cc: New files for qr
classes generated from CmplxQR.cc, CmplxQR.h, base-qr.cc, base-qr.h,
dbleQR.cc, dbleQR.h, fCmplxQR.cc, fCmplxQR.h, floatQR.cc, and
floatQR.h with classes converted to templates.

* liboctave/numeric/module.mk: Update.

* qz.cc, qr.cc, CmplxQRP.cc, CmplxQRP.h, dbleQRP.cc, dbleQRP.h,
fCmplxQRP.cc fCmplxQRP.h, floatQRP.cc, floatQRP.h, mx-defs.h,
mx-ext.h: Use new classes.

diff -r 48d82f74243e -r eb1524b07fe3 libinterp/corefcn/qz.cc
--- a/libinterp/corefcn/qz.cc	Tue Feb 16 23:34:44 2016 -0500
+++ b/libinterp/corefcn/qz.cc	Wed Feb 17 02:57:21 2016 -0500
@@ -37,11 +37,9 @@
 #include <iostream>
 #include <iomanip>
 
-#include "CmplxQRP.h"
-#include "CmplxQR.h"
-#include "dbleQR.h"
 #include "f77-fcn.h"
 #include "lo-math.h"
+#include "qr.h"
 #include "quit.h"
 
 #include "defun.h"
@@ -616,7 +614,7 @@
       // Complex case.
 
       // The QR decomposition of cbb.
-      ComplexQR cbqr (cbb);
+      qr<ComplexMatrix> cbqr (cbb);
       // The R matrix of QR decomposition for cbb.
       cbb = cbqr.R ();
       // (Q*)caa for following work.
@@ -681,7 +679,7 @@
 #endif
 
       // Compute the QR factorization of bb.
-      QR bqr (bb);
+      qr<Matrix> bqr (bb);
 
 #ifdef DEBUG
       std::cout << "qz: qr (bb) done; now peforming qz decomposition"
diff -r 48d82f74243e -r eb1524b07fe3 libinterp/dldfcn/qr.cc
--- a/libinterp/dldfcn/qr.cc	Tue Feb 16 23:34:44 2016 -0500
+++ b/libinterp/dldfcn/qr.cc	Wed Feb 17 02:57:21 2016 -0500
@@ -26,14 +26,11 @@
 #  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 "qr.h"
 #include "sparse-qr.h"
 
 
@@ -45,7 +42,7 @@
 
 template <typename MT>
 static octave_value
-get_qr_r (const base_qr<MT>& fact)
+get_qr_r (const qr<MT>& fact)
 {
   MT R = fact.R ();
   if (R.is_square () && fact.regular ())
@@ -54,6 +51,15 @@
     return R;
 }
 
+template <typename T>
+static typename qr<T>::type
+qr_type (int nargin, int nargout)
+{
+  return ((nargout == 0 || nargout == 1)
+          ? qr<T>::raw
+          : (nargin == 2) ? qr<T>::economy : qr<T>::std);
+}
+
 // [Q, R] = qr (X):      form Q unitary and R upper triangular such
 //                        that Q * R = X
 //
@@ -270,14 +276,13 @@
     }
   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 ())
             {
+              qr<FloatMatrix>::type type
+                = qr_type<FloatMatrix> (nargin, nargout);
+
               FloatMatrix m = arg.float_matrix_value ();
 
               switch (nargout)
@@ -285,14 +290,14 @@
                 case 0:
                 case 1:
                   {
-                    FloatQR fact (m, type);
+                    qr<FloatMatrix> fact (m, type);
                     retval = ovl (fact.R ());
                   }
                   break;
 
                 case 2:
                   {
-                    FloatQR fact (m, type);
+                    qr<FloatMatrix> fact (m, type);
                     retval = ovl (fact.Q (), get_qr_r (fact));
                   }
                   break;
@@ -301,7 +306,7 @@
                   {
                     FloatQRP fact (m, type);
 
-                    if (type == QR::economy)
+                    if (type == qr<FloatMatrix>::economy)
                       retval = ovl (fact.Q (), get_qr_r (fact), fact.Pvec ());
                     else
                       retval = ovl (fact.Q (), get_qr_r (fact), fact.P ());
@@ -311,6 +316,9 @@
             }
           else if (arg.is_complex_type ())
             {
+              qr<FloatComplexMatrix>::type type
+                = qr_type<FloatComplexMatrix> (nargin, nargout);
+
               FloatComplexMatrix m = arg.float_complex_matrix_value ();
 
               switch (nargout)
@@ -318,14 +326,14 @@
                 case 0:
                 case 1:
                   {
-                    FloatComplexQR fact (m, type);
+                    qr<FloatComplexMatrix> fact (m, type);
                     retval = ovl (fact.R ());
                   }
                   break;
 
                 case 2:
                   {
-                    FloatComplexQR fact (m, type);
+                    qr<FloatComplexMatrix> fact (m, type);
                     retval = ovl (fact.Q (), get_qr_r (fact));
                   }
                   break;
@@ -333,7 +341,7 @@
                 default:
                   {
                     FloatComplexQRP fact (m, type);
-                    if (type == QR::economy)
+                    if (type == qr<FloatComplexMatrix>::economy)
                       retval = ovl (fact.Q (), get_qr_r (fact), fact.Pvec ());
                     else
                       retval = ovl (fact.Q (), get_qr_r (fact), fact.P ());
@@ -346,6 +354,8 @@
         {
           if (arg.is_real_type ())
             {
+              qr<Matrix>::type type = qr_type<Matrix> (nargin, nargout);
+
               Matrix m = arg.matrix_value ();
 
               switch (nargout)
@@ -353,14 +363,14 @@
                 case 0:
                 case 1:
                   {
-                    QR fact (m, type);
+                    qr<Matrix> fact (m, type);
                     retval = ovl (fact.R ());
                   }
                   break;
 
                 case 2:
                   {
-                    QR fact (m, type);
+                    qr<Matrix> fact (m, type);
                     retval = ovl (fact.Q (), get_qr_r (fact));
                   }
                   break;
@@ -368,7 +378,7 @@
                 default:
                   {
                     QRP fact (m, type);
-                    if (type == QR::economy)
+                    if (type == qr<Matrix>::economy)
                       retval = ovl (fact.Q (), get_qr_r (fact), fact.Pvec ());
                     else
                       retval = ovl (fact.Q (), get_qr_r (fact), fact.P ());
@@ -378,6 +388,9 @@
             }
           else if (arg.is_complex_type ())
             {
+              qr<ComplexMatrix>::type type
+                = qr_type<ComplexMatrix> (nargin, nargout);
+
               ComplexMatrix m = arg.complex_matrix_value ();
 
               switch (nargout)
@@ -385,14 +398,14 @@
                 case 0:
                 case 1:
                   {
-                    ComplexQR fact (m, type);
+                    qr<ComplexMatrix> fact (m, type);
                     retval = ovl (fact.R ());
                   }
                   break;
 
                 case 2:
                   {
-                    ComplexQR fact (m, type);
+                    qr<ComplexMatrix> fact (m, type);
                     retval = ovl (fact.Q (), get_qr_r (fact));
                   }
                   break;
@@ -400,7 +413,7 @@
                 default:
                   {
                     ComplexQRP fact (m, type);
-                    if (type == QR::economy)
+                    if (type == qr<ComplexMatrix>::economy)
                       retval = ovl (fact.Q (), get_qr_r (fact), fact.Pvec ());
                     else
                       retval = ovl (fact.Q (), get_qr_r (fact), fact.P ());
@@ -763,7 +776,7 @@
           FloatMatrix u = argu.float_matrix_value ();
           FloatMatrix v = argv.float_matrix_value ();
 
-          FloatQR fact (Q, R);
+          qr<FloatMatrix> fact (Q, R);
           fact.update (u, v);
 
           retval = ovl (fact.Q (), get_qr_r (fact));
@@ -775,7 +788,7 @@
           Matrix u = argu.matrix_value ();
           Matrix v = argv.matrix_value ();
 
-          QR fact (Q, R);
+          qr<Matrix> fact (Q, R);
           fact.update (u, v);
 
           retval = ovl (fact.Q (), get_qr_r (fact));
@@ -792,7 +805,7 @@
           FloatComplexMatrix u = argu.float_complex_matrix_value ();
           FloatComplexMatrix v = argv.float_complex_matrix_value ();
 
-          FloatComplexQR fact (Q, R);
+          qr<FloatComplexMatrix> fact (Q, R);
           fact.update (u, v);
 
           retval = ovl (fact.Q (), get_qr_r (fact));
@@ -804,7 +817,7 @@
           ComplexMatrix u = argu.complex_matrix_value ();
           ComplexMatrix v = argv.complex_matrix_value ();
 
-          ComplexQR fact (Q, R);
+          qr<ComplexMatrix> fact (Q, R);
           fact.update (u, v);
 
           retval = ovl (fact.Q (), get_qr_r (fact));
@@ -947,7 +960,7 @@
           FloatMatrix R = argr.float_matrix_value ();
           FloatMatrix x = argx.float_matrix_value ();
 
-          FloatQR fact (Q, R);
+          qr<FloatMatrix> fact (Q, R);
 
           if (col)
             fact.insert_col (x, j-one);
@@ -962,7 +975,7 @@
           Matrix R = argr.matrix_value ();
           Matrix x = argx.matrix_value ();
 
-          QR fact (Q, R);
+          qr<Matrix> fact (Q, R);
 
           if (col)
             fact.insert_col (x, j-one);
@@ -985,7 +998,7 @@
           FloatComplexMatrix x =
             argx.float_complex_matrix_value ();
 
-          FloatComplexQR fact (Q, R);
+          qr<FloatComplexMatrix> fact (Q, R);
 
           if (col)
             fact.insert_col (x, j-one);
@@ -1000,7 +1013,7 @@
           ComplexMatrix R = argr.complex_matrix_value ();
           ComplexMatrix x = argx.complex_matrix_value ();
 
-          ComplexQR fact (Q, R);
+          qr<ComplexMatrix> fact (Q, R);
 
           if (col)
             fact.insert_col (x, j-one);
@@ -1138,7 +1151,7 @@
           FloatMatrix Q = argq.float_matrix_value ();
           FloatMatrix R = argr.float_matrix_value ();
 
-          FloatQR fact (Q, R);
+          qr<FloatMatrix> fact (Q, R);
 
           if (col)
             fact.delete_col (j-one);
@@ -1152,7 +1165,7 @@
           Matrix Q = argq.matrix_value ();
           Matrix R = argr.matrix_value ();
 
-          QR fact (Q, R);
+          qr<Matrix> fact (Q, R);
 
           if (col)
             fact.delete_col (j-one);
@@ -1172,7 +1185,7 @@
           FloatComplexMatrix R =
             argr.float_complex_matrix_value ();
 
-          FloatComplexQR fact (Q, R);
+          qr<FloatComplexMatrix> fact (Q, R);
 
           if (col)
             fact.delete_col (j-one);
@@ -1186,7 +1199,7 @@
           ComplexMatrix Q = argq.complex_matrix_value ();
           ComplexMatrix R = argr.complex_matrix_value ();
 
-          ComplexQR fact (Q, R);
+          qr<ComplexMatrix> fact (Q, R);
 
           if (col)
             fact.delete_col (j-one);
@@ -1364,7 +1377,7 @@
           FloatMatrix Q = argq.float_matrix_value ();
           FloatMatrix R = argr.float_matrix_value ();
 
-          FloatQR fact (Q, R);
+          qr<FloatMatrix> fact (Q, R);
           fact.shift_cols (i-1, j-1);
 
           retval = ovl (fact.Q (), get_qr_r (fact));
@@ -1374,7 +1387,7 @@
           Matrix Q = argq.matrix_value ();
           Matrix R = argr.matrix_value ();
 
-          QR fact (Q, R);
+          qr<Matrix> fact (Q, R);
           fact.shift_cols (i-1, j-1);
 
           retval = ovl (fact.Q (), get_qr_r (fact));
@@ -1389,7 +1402,7 @@
           FloatComplexMatrix Q = argq.float_complex_matrix_value ();
           FloatComplexMatrix R = argr.float_complex_matrix_value ();
 
-          FloatComplexQR fact (Q, R);
+          qr<FloatComplexMatrix> fact (Q, R);
           fact.shift_cols (i-1, j-1);
 
           retval = ovl (fact.Q (), get_qr_r (fact));
@@ -1399,7 +1412,7 @@
           ComplexMatrix Q = argq.complex_matrix_value ();
           ComplexMatrix R = argr.complex_matrix_value ();
 
-          ComplexQR fact (Q, R);
+          qr<ComplexMatrix> fact (Q, R);
           fact.shift_cols (i-1, j-1);
 
           retval = ovl (fact.Q (), get_qr_r (fact));
diff -r 48d82f74243e -r eb1524b07fe3 liboctave/numeric/CmplxQR.cc
--- a/liboctave/numeric/CmplxQR.cc	Tue Feb 16 23:34:44 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,684 +0,0 @@
-/*
-
-Copyright (C) 1994-2015 John W. Eaton
-Copyright (C) 2008-2009 Jaroslav Hajek
-Copyright (C) 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 "f77-fcn.h"
-#include "lo-error.h"
-#include "Range.h"
-#include "idx-vector.h"
-#include "oct-locbuf.h"
-
-#include "base-qr.cc"
-
-template class base_qr<ComplexMatrix>;
-
-extern "C"
-{
-  F77_RET_T
-  F77_FUNC (zgeqrf, ZGEQRF) (const octave_idx_type&, const octave_idx_type&,
-                             Complex*, const octave_idx_type&, Complex*,
-                             Complex*, const octave_idx_type&,
-                             octave_idx_type&);
-
-  F77_RET_T
-  F77_FUNC (zungqr, ZUNGQR) (const octave_idx_type&, const octave_idx_type&,
-                             const octave_idx_type&, Complex*,
-                             const octave_idx_type&, Complex*, Complex*,
-                             const octave_idx_type&, octave_idx_type&);
-
-#ifdef HAVE_QRUPDATE
-
-  F77_RET_T
-  F77_FUNC (zqr1up, ZQR1UP) (const octave_idx_type&, const octave_idx_type&,
-                             const octave_idx_type&, Complex*,
-                             const octave_idx_type&, Complex*,
-                             const octave_idx_type&, Complex*,
-                             Complex*, Complex*, double*);
-
-  F77_RET_T
-  F77_FUNC (zqrinc, ZQRINC) (const octave_idx_type&, const octave_idx_type&,
-                             const octave_idx_type&, Complex*,
-                             const octave_idx_type&, Complex*,
-                             const octave_idx_type&, const octave_idx_type&,
-                             const Complex*, double*);
-
-  F77_RET_T
-  F77_FUNC (zqrdec, ZQRDEC) (const octave_idx_type&, const octave_idx_type&,
-                             const octave_idx_type&, Complex*,
-                             const octave_idx_type&, Complex*,
-                             const octave_idx_type&, const octave_idx_type&,
-                             double*);
-
-  F77_RET_T
-  F77_FUNC (zqrinr, ZQRINR) (const octave_idx_type&, const octave_idx_type&,
-                             Complex*, const octave_idx_type&, Complex*,
-                             const octave_idx_type&, const octave_idx_type&,
-                             const Complex*, double*);
-
-  F77_RET_T
-  F77_FUNC (zqrder, ZQRDER) (const octave_idx_type&, const octave_idx_type&,
-                             Complex*, const octave_idx_type&, Complex*,
-                             const octave_idx_type&, const octave_idx_type&,
-                             Complex*, double*);
-
-  F77_RET_T
-  F77_FUNC (zqrshc, ZQRSHC) (const octave_idx_type&, const octave_idx_type&,
-                             const octave_idx_type&, Complex*,
-                             const octave_idx_type&, Complex*,
-                             const octave_idx_type&, const octave_idx_type&,
-                             const octave_idx_type&, Complex*, double*);
-
-#endif
-}
-
-ComplexQR::ComplexQR (const ComplexMatrix& a, qr_type_t qr_type)
-{
-  init (a, qr_type);
-}
-
-void
-ComplexQR::init (const ComplexMatrix& a, qr_type_t qr_type)
-{
-  octave_idx_type m = a.rows ();
-  octave_idx_type n = a.cols ();
-
-  octave_idx_type min_mn = m < n ? m : n;
-  OCTAVE_LOCAL_BUFFER (Complex, tau, min_mn);
-
-  octave_idx_type info = 0;
-
-  ComplexMatrix afact = a;
-  if (m > n && qr_type == qr_type_std)
-    afact.resize (m, m);
-
-  if (m > 0)
-    {
-      // workspace query.
-      Complex clwork;
-      F77_XFCN (zgeqrf, ZGEQRF, (m, n, afact.fortran_vec (), m, tau,
-                                 &clwork, -1, info));
-
-      // allocate buffer and do the job.
-      octave_idx_type lwork = clwork.real ();
-      lwork = std::max (lwork, static_cast<octave_idx_type> (1));
-      OCTAVE_LOCAL_BUFFER (Complex, work, lwork);
-      F77_XFCN (zgeqrf, ZGEQRF, (m, n, afact.fortran_vec (), m, tau,
-                                 work, lwork, info));
-    }
-
-  form (n, afact, tau, qr_type);
-}
-
-void ComplexQR::form (octave_idx_type n, ComplexMatrix& afact,
-                      Complex *tau, qr_type_t qr_type)
-{
-  octave_idx_type m = afact.rows ();
-  octave_idx_type min_mn = std::min (m, n);
-  octave_idx_type info;
-
-  if (qr_type == qr_type_raw)
-    {
-      for (octave_idx_type j = 0; j < min_mn; j++)
-        {
-          octave_idx_type limit = j < min_mn - 1 ? j : min_mn - 1;
-          for (octave_idx_type i = limit + 1; i < m; i++)
-            afact.elem (i, j) *= tau[j];
-        }
-
-      r = afact;
-    }
-  else
-    {
-      // Attempt to minimize copying.
-      if (m >= n)
-        {
-          // afact will become q.
-          q = afact;
-          octave_idx_type k = qr_type == qr_type_economy ? n : m;
-          r = ComplexMatrix (k, n);
-          for (octave_idx_type j = 0; j < n; j++)
-            {
-              octave_idx_type i = 0;
-              for (; i <= j; i++)
-                r.xelem (i, j) = afact.xelem (i, j);
-              for (; i < k; i++)
-                r.xelem (i, j) = 0;
-            }
-          afact = ComplexMatrix (); // optimize memory
-        }
-      else
-        {
-          // afact will become r.
-          q = ComplexMatrix (m, m);
-          for (octave_idx_type j = 0; j < m; j++)
-            for (octave_idx_type i = j + 1; i < m; i++)
-              {
-                q.xelem (i, j) = afact.xelem (i, j);
-                afact.xelem (i, j) = 0;
-              }
-          r = afact;
-        }
-
-
-      if (m > 0)
-        {
-          octave_idx_type k = q.columns ();
-          // workspace query.
-          Complex clwork;
-          F77_XFCN (zungqr, ZUNGQR, (m, k, min_mn, q.fortran_vec (), m, tau,
-                                     &clwork, -1, info));
-
-          // allocate buffer and do the job.
-          octave_idx_type lwork = clwork.real ();
-          lwork = std::max (lwork, static_cast<octave_idx_type> (1));
-          OCTAVE_LOCAL_BUFFER (Complex, work, lwork);
-          F77_XFCN (zungqr, ZUNGQR, (m, k, min_mn, q.fortran_vec (), m, tau,
-                                     work, lwork, info));
-        }
-    }
-}
-
-#ifdef HAVE_QRUPDATE
-
-void
-ComplexQR::update (const ComplexColumnVector& u, const ComplexColumnVector& v)
-{
-  octave_idx_type m = q.rows ();
-  octave_idx_type n = r.columns ();
-  octave_idx_type k = q.columns ();
-
-  if (u.numel () != m || v.numel () != n)
-    (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
-
-  ComplexColumnVector utmp = u;
-  ComplexColumnVector vtmp = v;
-  OCTAVE_LOCAL_BUFFER (Complex, w, k);
-  OCTAVE_LOCAL_BUFFER (double, rw, k);
-  F77_XFCN (zqr1up, ZQR1UP, (m, n, k, q.fortran_vec (),
-                             m, r.fortran_vec (), k,
-                             utmp.fortran_vec (), vtmp.fortran_vec (),
-                             w, rw));
-}
-
-void
-ComplexQR::update (const ComplexMatrix& u, const ComplexMatrix& v)
-{
-  octave_idx_type m = q.rows ();
-  octave_idx_type n = r.columns ();
-  octave_idx_type k = q.columns ();
-
-  if (u.rows () != m || v.rows () != n || u.cols () != v.cols ())
-    (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
-
-  OCTAVE_LOCAL_BUFFER (Complex, w, k);
-  OCTAVE_LOCAL_BUFFER (double, rw, k);
-  for (volatile octave_idx_type i = 0; i < u.cols (); i++)
-    {
-      ComplexColumnVector utmp = u.column (i);
-      ComplexColumnVector vtmp = v.column (i);
-      F77_XFCN (zqr1up, ZQR1UP, (m, n, k, q.fortran_vec (),
-                                 m, r.fortran_vec (), k,
-                                 utmp.fortran_vec (), vtmp.fortran_vec (),
-                                 w, rw));
-    }
-}
-
-void
-ComplexQR::insert_col (const ComplexColumnVector& u, octave_idx_type j)
-{
-  octave_idx_type m = q.rows ();
-  octave_idx_type n = r.columns ();
-  octave_idx_type k = q.columns ();
-
-  if (u.numel () != m)
-    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
-  if (j < 0 || j > n)
-    (*current_liboctave_error_handler) ("qrinsert: index out of range");
-
-  if (k < m)
-    {
-      q.resize (m, k+1);
-      r.resize (k+1, n+1);
-    }
-  else
-    {
-      r.resize (k, n+1);
-    }
-
-  ComplexColumnVector utmp = u;
-  OCTAVE_LOCAL_BUFFER (double, rw, k);
-  F77_XFCN (zqrinc, ZQRINC, (m, n, k, q.fortran_vec (), q.rows (),
-                             r.fortran_vec (), r.rows (), j + 1,
-                             utmp.data (), rw));
-}
-
-void
-ComplexQR::insert_col (const ComplexMatrix& u, const Array<octave_idx_type>& j)
-{
-  octave_idx_type m = q.rows ();
-  octave_idx_type n = r.columns ();
-  octave_idx_type k = q.columns ();
-
-  Array<octave_idx_type> jsi;
-  Array<octave_idx_type> js = j.sort (jsi, 0, ASCENDING);
-  octave_idx_type nj = js.numel ();
-  bool dups = false;
-  for (octave_idx_type i = 0; i < nj - 1; i++)
-    dups = dups && js(i) == js(i+1);
-
-  if (dups)
-    (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
-  if (u.numel () != m || u.columns () != nj)
-    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
-  if (nj > 0 && (js(0) < 0 || js(nj-1) > n))
-    (*current_liboctave_error_handler) ("qrinsert: index out of range");
-
-  if (nj > 0)
-    {
-      octave_idx_type kmax = std::min (k + nj, m);
-      if (k < m)
-        {
-          q.resize (m, kmax);
-          r.resize (kmax, n + nj);
-        }
-      else
-        {
-          r.resize (k, n + nj);
-        }
-
-      OCTAVE_LOCAL_BUFFER (double, rw, kmax);
-      for (volatile octave_idx_type i = 0; i < js.numel (); i++)
-        {
-          octave_idx_type ii = i;
-          ComplexColumnVector utmp = u.column (jsi(i));
-          F77_XFCN (zqrinc, ZQRINC, (m, n + ii, std::min (kmax, k + ii),
-                                     q.fortran_vec (), q.rows (),
-                                     r.fortran_vec (), r.rows (), js(ii) + 1,
-                                     utmp.data (), rw));
-        }
-    }
-}
-
-void
-ComplexQR::delete_col (octave_idx_type j)
-{
-  octave_idx_type m = q.rows ();
-  octave_idx_type k = r.rows ();
-  octave_idx_type n = r.columns ();
-
-  if (j < 0 || j > n-1)
-    (*current_liboctave_error_handler) ("qrdelete: index out of range");
-
-  OCTAVE_LOCAL_BUFFER (double, rw, k);
-  F77_XFCN (zqrdec, ZQRDEC, (m, n, k, q.fortran_vec (), q.rows (),
-                             r.fortran_vec (), r.rows (), j + 1, rw));
-
-  if (k < m)
-    {
-      q.resize (m, k-1);
-      r.resize (k-1, n-1);
-    }
-  else
-    {
-      r.resize (k, n-1);
-    }
-}
-
-void
-ComplexQR::delete_col (const Array<octave_idx_type>& j)
-{
-  octave_idx_type m = q.rows ();
-  octave_idx_type n = r.columns ();
-  octave_idx_type k = q.columns ();
-
-  Array<octave_idx_type> jsi;
-  Array<octave_idx_type> js = j.sort (jsi, 0, DESCENDING);
-  octave_idx_type nj = js.numel ();
-  bool dups = false;
-  for (octave_idx_type i = 0; i < nj - 1; i++)
-    dups = dups && js(i) == js(i+1);
-
-  if (dups)
-    (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
-  if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0))
-    (*current_liboctave_error_handler) ("qrinsert: index out of range");
-
-  if (nj > 0)
-    {
-      OCTAVE_LOCAL_BUFFER (double, rw, k);
-      for (volatile octave_idx_type i = 0; i < js.numel (); i++)
-        {
-          octave_idx_type ii = i;
-          F77_XFCN (zqrdec, ZQRDEC, (m, n - ii, k == m ? k : k - ii,
-                                     q.fortran_vec (), q.rows (),
-                                     r.fortran_vec (), r.rows (),
-                                     js(ii) + 1, rw));
-        }
-      if (k < m)
-        {
-          q.resize (m, k - nj);
-          r.resize (k - nj, n - nj);
-        }
-      else
-        {
-          r.resize (k, n - nj);
-        }
-
-    }
-}
-
-void
-ComplexQR::insert_row (const ComplexRowVector& u, octave_idx_type j)
-{
-  octave_idx_type m = r.rows ();
-  octave_idx_type n = r.columns ();
-  octave_idx_type k = std::min (m, n);
-
-  if (! q.is_square () || u.numel () != n)
-    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
-  if (j < 0 || j > m)
-    (*current_liboctave_error_handler) ("qrinsert: index out of range");
-
-  q.resize (m + 1, m + 1);
-  r.resize (m + 1, n);
-  ComplexRowVector utmp = u;
-  OCTAVE_LOCAL_BUFFER (double, rw, k);
-  F77_XFCN (zqrinr, ZQRINR, (m, n, q.fortran_vec (), q.rows (),
-                             r.fortran_vec (), r.rows (),
-                             j + 1, utmp.fortran_vec (), rw));
-
-}
-
-void
-ComplexQR::delete_row (octave_idx_type j)
-{
-  octave_idx_type m = r.rows ();
-  octave_idx_type n = r.columns ();
-
-  if (! q.is_square ())
-    (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch");
-  if (j < 0 || j > m-1)
-    (*current_liboctave_error_handler) ("qrdelete: index out of range");
-
-  OCTAVE_LOCAL_BUFFER (Complex, w, m);
-  OCTAVE_LOCAL_BUFFER (double, rw, m);
-  F77_XFCN (zqrder, ZQRDER, (m, n, q.fortran_vec (), q.rows (),
-                             r.fortran_vec (), r.rows (), j + 1,
-                             w, rw));
-
-  q.resize (m - 1, m - 1);
-  r.resize (m - 1, n);
-}
-
-void
-ComplexQR::shift_cols (octave_idx_type i, octave_idx_type j)
-{
-  octave_idx_type m = q.rows ();
-  octave_idx_type k = r.rows ();
-  octave_idx_type n = r.columns ();
-
-  if (i < 0 || i > n-1 || j < 0 || j > n-1)
-    (*current_liboctave_error_handler) ("qrshift: index out of range");
-
-  OCTAVE_LOCAL_BUFFER (Complex, w, k);
-  OCTAVE_LOCAL_BUFFER (double, rw, k);
-  F77_XFCN (zqrshc, ZQRSHC, (m, n, k,
-                             q.fortran_vec (), q.rows (),
-                             r.fortran_vec (), r.rows (),
-                             i + 1, j + 1, w, rw));
-}
-
-#else
-
-// Replacement update methods.
-
-void
-ComplexQR::update (const ComplexColumnVector& u, const ComplexColumnVector& v)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type m = q.rows ();
-  octave_idx_type n = r.columns ();
-
-  if (u.numel () != m || v.numel () != n)
-    (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
-
-  init (q*r + ComplexMatrix (u) * ComplexMatrix (v).hermitian (),
-        get_type ());
-}
-
-void
-ComplexQR::update (const ComplexMatrix& u, const ComplexMatrix& v)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type m = q.rows ();
-  octave_idx_type n = r.columns ();
-
-  if (u.rows () != m || v.rows () != n || u.cols () != v.cols ())
-    (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
-
-  init (q*r + u * v.hermitian (), get_type ());
-}
-
-static
-ComplexMatrix insert_col (const ComplexMatrix& a, octave_idx_type i,
-                          const ComplexColumnVector& x)
-{
-  ComplexMatrix retval (a.rows (), a.columns () + 1);
-  retval.assign (idx_vector::colon, idx_vector (0, i),
-                 a.index (idx_vector::colon, idx_vector (0, i)));
-  retval.assign (idx_vector::colon, idx_vector (i), x);
-  retval.assign (idx_vector::colon, idx_vector (i+1, retval.columns ()),
-                 a.index (idx_vector::colon, idx_vector (i, a.columns ())));
-  return retval;
-}
-
-static
-ComplexMatrix insert_row (const ComplexMatrix& a, octave_idx_type i,
-                          const ComplexRowVector& x)
-{
-  ComplexMatrix retval (a.rows () + 1, a.columns ());
-  retval.assign (idx_vector (0, i), idx_vector::colon,
-                 a.index (idx_vector (0, i), idx_vector::colon));
-  retval.assign (idx_vector (i), idx_vector::colon, x);
-  retval.assign (idx_vector (i+1, retval.rows ()), idx_vector::colon,
-                 a.index (idx_vector (i, a.rows ()), idx_vector::colon));
-  return retval;
-}
-
-static
-ComplexMatrix delete_col (const ComplexMatrix& a, octave_idx_type i)
-{
-  ComplexMatrix retval = a;
-  retval.delete_elements (1, idx_vector (i));
-  return retval;
-}
-
-static
-ComplexMatrix delete_row (const ComplexMatrix& a, octave_idx_type i)
-{
-  ComplexMatrix retval = a;
-  retval.delete_elements (0, idx_vector (i));
-  return retval;
-}
-
-static
-ComplexMatrix shift_cols (const ComplexMatrix& a,
-                          octave_idx_type i, octave_idx_type j)
-{
-  octave_idx_type n = a.columns ();
-  Array<octave_idx_type> p (dim_vector (n, 1));
-  for (octave_idx_type k = 0; k < n; k++) p(k) = k;
-  if (i < j)
-    {
-      for (octave_idx_type k = i; k < j; k++) p(k) = k+1;
-      p(j) = i;
-    }
-  else if (j < i)
-    {
-      p(j) = i;
-      for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1;
-    }
-
-  return a.index (idx_vector::colon, idx_vector (p));
-}
-
-void
-ComplexQR::insert_col (const ComplexColumnVector& u, octave_idx_type j)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type m = q.rows ();
-  octave_idx_type n = r.columns ();
-
-  if (u.numel () != m)
-    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
-  if (j < 0 || j > n)
-    (*current_liboctave_error_handler) ("qrinsert: index out of range");
-
-  init (::insert_col (q*r, j, u), get_type ());
-}
-
-void
-ComplexQR::insert_col (const ComplexMatrix& u, const Array<octave_idx_type>& j)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type m = q.rows ();
-  octave_idx_type n = r.columns ();
-
-  Array<octave_idx_type> jsi;
-  Array<octave_idx_type> js = j.sort (jsi, 0, ASCENDING);
-  octave_idx_type nj = js.numel ();
-  bool dups = false;
-  for (octave_idx_type i = 0; i < nj - 1; i++)
-    dups = dups && js(i) == js(i+1);
-
-  if (dups)
-    (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
-  if (u.numel () != m || u.columns () != nj)
-    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
-  if (nj > 0 && (js(0) < 0 || js(nj-1) > n))
-    (*current_liboctave_error_handler) ("qrinsert: index out of range");
-
-  if (nj > 0)
-    {
-      ComplexMatrix a = q*r;
-      for (octave_idx_type i = 0; i < js.numel (); i++)
-        a = ::insert_col (a, js(i), u.column (i));
-      init (a, get_type ());
-    }
-}
-
-void
-ComplexQR::delete_col (octave_idx_type j)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type n = r.columns ();
-
-  if (j < 0 || j > n-1)
-    (*current_liboctave_error_handler) ("qrdelete: index out of range");
-
-  init (::delete_col (q*r, j), get_type ());
-}
-
-void
-ComplexQR::delete_col (const Array<octave_idx_type>& j)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type n = r.columns ();
-
-  Array<octave_idx_type> jsi;
-  Array<octave_idx_type> js = j.sort (jsi, 0, DESCENDING);
-  octave_idx_type nj = js.numel ();
-  bool dups = false;
-  for (octave_idx_type i = 0; i < nj - 1; i++)
-    dups = dups && js(i) == js(i+1);
-
-  if (dups)
-    (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
-  if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0))
-    (*current_liboctave_error_handler) ("qrinsert: index out of range");
-
-  if (nj > 0)
-    {
-      ComplexMatrix a = q*r;
-      for (octave_idx_type i = 0; i < js.numel (); i++)
-        a = ::delete_col (a, js(i));
-      init (a, get_type ());
-    }
-}
-
-void
-ComplexQR::insert_row (const ComplexRowVector& u, octave_idx_type j)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type m = r.rows ();
-  octave_idx_type n = r.columns ();
-
-  if (! q.is_square () || u.numel () != n)
-    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
-  if (j < 0 || j > m)
-    (*current_liboctave_error_handler) ("qrinsert: index out of range");
-
-  init (::insert_row (q*r, j, u), get_type ());
-}
-
-void
-ComplexQR::delete_row (octave_idx_type j)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type m = r.rows ();
-
-  if (! q.is_square ())
-    (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch");
-  if (j < 0 || j > m-1)
-    (*current_liboctave_error_handler) ("qrdelete: index out of range");
-
-  init (::delete_row (q*r, j), get_type ());
-}
-
-void
-ComplexQR::shift_cols (octave_idx_type i, octave_idx_type j)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type n = r.columns ();
-
-  if (i < 0 || i > n-1 || j < 0 || j > n-1)
-    (*current_liboctave_error_handler) ("qrshift: index out of range");
-
-  init (::shift_cols (q*r, i, j), get_type ());
-}
-
-#endif
diff -r 48d82f74243e -r eb1524b07fe3 liboctave/numeric/CmplxQR.h
--- a/liboctave/numeric/CmplxQR.h	Tue Feb 16 23:34:44 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,78 +0,0 @@
-/*
-
-Copyright (C) 1994-2015 John W. Eaton
-Copyright (C) 2008-2009 Jaroslav Hajek
-
-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/>.
-
-*/
-
-#if ! defined (octave_CmplxQR_h)
-#define octave_CmplxQR_h 1
-
-#include "octave-config.h"
-
-#include <iosfwd>
-
-#include "CMatrix.h"
-#include "CColVector.h"
-#include "CRowVector.h"
-#include "base-qr.h"
-
-
-class
-OCTAVE_API
-ComplexQR : public base_qr<ComplexMatrix>
-{
-public:
-
-  ComplexQR (void) : base_qr<ComplexMatrix> () { }
-
-  ComplexQR (const ComplexMatrix&, qr_type_t = qr_type_std);
-
-  ComplexQR (const ComplexMatrix& qx, const ComplexMatrix& rx)
-    : base_qr<ComplexMatrix> (qx, rx) { }
-
-  ComplexQR (const ComplexQR& a) : base_qr<ComplexMatrix> (a) { }
-
-  void init (const ComplexMatrix&, qr_type_t = qr_type_std);
-
-  void update (const ComplexColumnVector& u, const ComplexColumnVector& v);
-
-  void update (const ComplexMatrix& u, const ComplexMatrix& v);
-
-  void insert_col (const ComplexColumnVector& u, octave_idx_type j);
-
-  void insert_col (const ComplexMatrix& u, const Array<octave_idx_type>& j);
-
-  void delete_col (octave_idx_type j);
-
-  void delete_col (const Array<octave_idx_type>& j);
-
-  void insert_row (const ComplexRowVector& u, octave_idx_type j);
-
-  void delete_row (octave_idx_type j);
-
-  void shift_cols (octave_idx_type i, octave_idx_type j);
-
-protected:
-
-  void form (octave_idx_type n, ComplexMatrix& afact,
-             Complex *tau, qr_type_t qr_type);
-};
-
-#endif
diff -r 48d82f74243e -r eb1524b07fe3 liboctave/numeric/CmplxQRP.cc
--- a/liboctave/numeric/CmplxQRP.cc	Tue Feb 16 23:34:44 2016 -0500
+++ b/liboctave/numeric/CmplxQRP.cc	Wed Feb 17 02:57:21 2016 -0500
@@ -42,18 +42,18 @@
                              octave_idx_type&);
 }
 
-// It would be best to share some of this code with ComplexQR class...
+// It would be best to share some of this code with qr<ComplexMatrix> class...
 
-ComplexQRP::ComplexQRP (const ComplexMatrix& a, qr_type_t qr_type)
-  : ComplexQR (), p ()
+ComplexQRP::ComplexQRP (const ComplexMatrix& a, type qr_type)
+  : qr<ComplexMatrix> (), p ()
 {
   init (a, qr_type);
 }
 
 void
-ComplexQRP::init (const ComplexMatrix& a, qr_type_t qr_type)
+ComplexQRP::init (const ComplexMatrix& a, type qr_type)
 {
-  assert (qr_type != qr_type_raw);
+  assert (qr_type != qr<ComplexMatrix>::raw);
 
   octave_idx_type m = a.rows ();
   octave_idx_type n = a.cols ();
@@ -64,7 +64,7 @@
   octave_idx_type info = 0;
 
   ComplexMatrix afact = a;
-  if (m > n && qr_type == qr_type_std)
+  if (m > n && qr_type == qr<ComplexMatrix>::std)
     afact.resize (m, m);
 
   MArray<octave_idx_type> jpvt (dim_vector (n, 1), 0);
diff -r 48d82f74243e -r eb1524b07fe3 liboctave/numeric/CmplxQRP.h
--- a/liboctave/numeric/CmplxQRP.h	Tue Feb 16 23:34:44 2016 -0500
+++ b/liboctave/numeric/CmplxQRP.h	Wed Feb 17 02:57:21 2016 -0500
@@ -27,27 +27,31 @@
 
 #include <iosfwd>
 
-#include "CmplxQR.h"
 #include "PermMatrix.h"
-#include "dColVector.h"
+#include "CColVector.h"
+#include "CMatrix.h"
+#include "CRowVector.h"
+#include "qr.h"
 
 class
 OCTAVE_API
-ComplexQRP : public ComplexQR
+ComplexQRP : public qr<ComplexMatrix>
 {
 public:
 
-  ComplexQRP (void) : ComplexQR (), p () { }
+  typedef qr<ComplexMatrix>::type type;
+
+  ComplexQRP (void) : qr<ComplexMatrix> (), p () { }
 
-  ComplexQRP (const ComplexMatrix&, qr_type_t = qr_type_std);
+  ComplexQRP (const ComplexMatrix&, type = qr<ComplexMatrix>::std);
 
-  ComplexQRP (const ComplexQRP& a) : ComplexQR (a), p (a.p) { }
+  ComplexQRP (const ComplexQRP& a) : qr<ComplexMatrix> (a), p (a.p) { }
 
   ComplexQRP& operator = (const ComplexQRP& a)
   {
     if (this != &a)
       {
-        ComplexQR::operator = (a);
+        qr<ComplexMatrix>::operator = (a);
         p = a.p;
       }
     return *this;
@@ -55,7 +59,7 @@
 
   ~ComplexQRP (void) { }
 
-  void init (const ComplexMatrix&, qr_type_t = qr_type_std);
+  void init (const ComplexMatrix&, type = qr<ComplexMatrix>::std);
 
   PermMatrix P (void) const { return p; }
 
diff -r 48d82f74243e -r eb1524b07fe3 liboctave/numeric/base-qr.cc
--- a/liboctave/numeric/base-qr.cc	Tue Feb 16 23:34:44 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,82 +0,0 @@
-/*
-
-Copyright (C) 2009-2015 Jaroslav Hajek
-
-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 "base-qr.h"
-
-template <typename qr_type>
-base_qr<qr_type>::base_qr (const qr_type& q_arg, const qr_type& r_arg)
-  : q (q_arg), r (r_arg)
-{
-  octave_idx_type q_nr = q.rows ();
-  octave_idx_type q_nc = q.columns ();
-  octave_idx_type r_nr = r.rows ();
-  octave_idx_type r_nc = r.columns ();
-
-  if (! (q_nc == r_nr && (q_nr == q_nc || (q_nr > q_nc && r_nr == r_nc))))
-    {
-      q = qr_type ();
-      r = qr_type ();
-
-      (*current_liboctave_error_handler) ("QR dimensions mismatch");
-    }
-}
-
-template <typename qr_type>
-qr_type_t
-base_qr<qr_type>::get_type (void) const
-{
-  qr_type_t retval;
-
-  if (! q.is_empty () && q.is_square ())
-    retval = qr_type_std;
-  else if (q.rows () > q.columns () && r.is_square ())
-    retval = qr_type_economy;
-  else
-    retval = qr_type_raw;
-
-  return retval;
-}
-
-template <typename qr_type>
-bool
-base_qr<qr_type>::regular (void) const
-{
-  bool retval = true;
-
-  octave_idx_type k = std::min (r.rows (), r.columns ());
-
-  for (octave_idx_type i = 0; i < k; i++)
-    {
-      if (r(i, i) == qr_elt_type ())
-        {
-          retval = false;
-          break;
-        }
-    }
-
-  return retval;
-}
-
diff -r 48d82f74243e -r eb1524b07fe3 liboctave/numeric/base-qr.h
--- a/liboctave/numeric/base-qr.h	Tue Feb 16 23:34:44 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,82 +0,0 @@
-/*
-
-Copyright (C) 2009-2015 Jaroslav Hajek
-
-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/>.
-
-*/
-
-#if ! defined (octave_base_qr_h)
-#define octave_base_qr_h 1
-
-#include "octave-config.h"
-
-#include "MArray.h"
-#include "dColVector.h"
-#include "PermMatrix.h"
-
-enum qr_type_t
-{
-  qr_type_std,
-  qr_type_raw,
-  qr_type_economy
-};
-
-template <typename qr_type>
-class
-base_qr
-{
-public:
-
-  typedef typename qr_type::element_type qr_elt_type;
-
-  base_qr (void) : q (), r () { }
-
-  base_qr (const qr_type& q, const qr_type& r);
-
-  base_qr (const base_qr& a) : q (a.q), r (a.r) { }
-
-  base_qr& operator = (const base_qr& a)
-  {
-    if (this != &a)
-      {
-        q = a.q;
-        r = a.r;
-      }
-
-    return *this;
-  }
-
-  virtual ~base_qr (void) { }
-
-  qr_type Q (void) const { return q; }
-
-  qr_type R (void) const { return r; }
-
-  qr_type_t get_type (void) const;
-
-  bool regular (void) const;
-
-protected:
-
-  qr_type q;
-  qr_type r;
-};
-
-extern void warn_qrupdate_once (void);
-
-#endif
diff -r 48d82f74243e -r eb1524b07fe3 liboctave/numeric/dbleQR.cc
--- a/liboctave/numeric/dbleQR.cc	Tue Feb 16 23:34:44 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,695 +0,0 @@
-/*
-
-Copyright (C) 1994-2015 John W. Eaton
-Copyright (C) 2008-2009 Jaroslav Hajek
-Copyright (C) 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 "dbleQR.h"
-#include "f77-fcn.h"
-#include "lo-error.h"
-#include "Range.h"
-#include "idx-vector.h"
-#include "oct-locbuf.h"
-
-#include "base-qr.cc"
-
-template class base_qr<Matrix>;
-
-extern "C"
-{
-  F77_RET_T
-  F77_FUNC (dgeqrf, DGEQRF) (const octave_idx_type&, const octave_idx_type&,
-                             double*, const octave_idx_type&, double*,
-                             double*, const octave_idx_type&,
-                             octave_idx_type&);
-
-  F77_RET_T
-  F77_FUNC (dorgqr, DORGQR) (const octave_idx_type&, const octave_idx_type&,
-                             const octave_idx_type&, double*,
-                             const octave_idx_type&, double*, double*,
-                             const octave_idx_type&, octave_idx_type&);
-
-#ifdef HAVE_QRUPDATE
-
-  F77_RET_T
-  F77_FUNC (dqr1up, DQR1UP) (const octave_idx_type&, const octave_idx_type&,
-                             const octave_idx_type&, double*,
-                             const octave_idx_type&, double*,
-                             const octave_idx_type&, double*, double*, double*);
-
-  F77_RET_T
-  F77_FUNC (dqrinc, DQRINC) (const octave_idx_type&, const octave_idx_type&,
-                             const octave_idx_type&, double*,
-                             const octave_idx_type&, double*,
-                             const octave_idx_type&, const octave_idx_type&,
-                             const double*, double*);
-
-  F77_RET_T
-  F77_FUNC (dqrdec, DQRDEC) (const octave_idx_type&, const octave_idx_type&,
-                             const octave_idx_type&, double*,
-                             const octave_idx_type&, double*,
-                             const octave_idx_type&, const octave_idx_type&,
-                             double*);
-
-  F77_RET_T
-  F77_FUNC (dqrinr, DQRINR) (const octave_idx_type&, const octave_idx_type&,
-                             double*, const octave_idx_type&, double*,
-                             const octave_idx_type&, const octave_idx_type&,
-                             const double*, double*);
-
-  F77_RET_T
-  F77_FUNC (dqrder, DQRDER) (const octave_idx_type&, const octave_idx_type&,
-                             double*, const octave_idx_type&, double*,
-                             const octave_idx_type&, const octave_idx_type&,
-                             double*);
-
-  F77_RET_T
-  F77_FUNC (dqrshc, DQRSHC) (const octave_idx_type&, const octave_idx_type&,
-                             const octave_idx_type&, double*,
-                             const octave_idx_type&, double*,
-                             const octave_idx_type&, const octave_idx_type&,
-                             const octave_idx_type&, double*);
-
-#endif
-}
-
-const QR::type QR::raw, QR::std, QR::economy;
-
-QR::QR (const Matrix& a, qr_type_t qr_type)
-{
-  init (a, qr_type);
-}
-
-void
-QR::init (const Matrix& a, qr_type_t qr_type)
-{
-  octave_idx_type m = a.rows ();
-  octave_idx_type n = a.cols ();
-
-  octave_idx_type min_mn = m < n ? m : n;
-  OCTAVE_LOCAL_BUFFER (double, tau, min_mn);
-
-  octave_idx_type info = 0;
-
-  Matrix afact = a;
-  if (m > n && qr_type == qr_type_std)
-    afact.resize (m, m);
-
-  if (m > 0)
-    {
-      // workspace query.
-      double rlwork;
-      F77_XFCN (dgeqrf, DGEQRF, (m, n, afact.fortran_vec (), m, tau,
-                                 &rlwork, -1, info));
-
-      // allocate buffer and do the job.
-      octave_idx_type lwork = rlwork;
-      lwork = std::max (lwork, static_cast<octave_idx_type> (1));
-      OCTAVE_LOCAL_BUFFER (double, work, lwork);
-      F77_XFCN (dgeqrf, DGEQRF, (m, n, afact.fortran_vec (), m, tau,
-                                 work, lwork, info));
-    }
-
-  form (n, afact, tau, qr_type);
-}
-
-void QR::form (octave_idx_type n, Matrix& afact,
-               double *tau, qr_type_t qr_type)
-{
-  octave_idx_type m = afact.rows ();
-  octave_idx_type min_mn = std::min (m, n);
-  octave_idx_type info;
-
-  if (qr_type == qr_type_raw)
-    {
-      for (octave_idx_type j = 0; j < min_mn; j++)
-        {
-          octave_idx_type limit = j < min_mn - 1 ? j : min_mn - 1;
-          for (octave_idx_type i = limit + 1; i < m; i++)
-            afact.elem (i, j) *= tau[j];
-        }
-
-      r = afact;
-    }
-  else
-    {
-      // Attempt to minimize copying.
-      if (m >= n)
-        {
-          // afact will become q.
-          q = afact;
-          octave_idx_type k = qr_type == qr_type_economy ? n : m;
-          r = Matrix (k, n);
-          for (octave_idx_type j = 0; j < n; j++)
-            {
-              octave_idx_type i = 0;
-              for (; i <= j; i++)
-                r.xelem (i, j) = afact.xelem (i, j);
-              for (; i < k; i++)
-                r.xelem (i, j) = 0;
-            }
-          afact = Matrix (); // optimize memory
-        }
-      else
-        {
-          // afact will become r.
-          q = Matrix (m, m);
-          for (octave_idx_type j = 0; j < m; j++)
-            for (octave_idx_type i = j + 1; i < m; i++)
-              {
-                q.xelem (i, j) = afact.xelem (i, j);
-                afact.xelem (i, j) = 0;
-              }
-          r = afact;
-        }
-
-
-      if (m > 0)
-        {
-          octave_idx_type k = q.columns ();
-          // workspace query.
-          double rlwork;
-          F77_XFCN (dorgqr, DORGQR, (m, k, min_mn, q.fortran_vec (), m, tau,
-                                     &rlwork, -1, info));
-
-          // allocate buffer and do the job.
-          octave_idx_type lwork = rlwork;
-          lwork = std::max (lwork, static_cast<octave_idx_type> (1));
-          OCTAVE_LOCAL_BUFFER (double, work, lwork);
-          F77_XFCN (dorgqr, DORGQR, (m, k, min_mn, q.fortran_vec (), m, tau,
-                                     work, lwork, info));
-        }
-    }
-}
-
-#ifdef HAVE_QRUPDATE
-
-void
-QR::update (const ColumnVector& u, const ColumnVector& v)
-{
-  octave_idx_type m = q.rows ();
-  octave_idx_type n = r.columns ();
-  octave_idx_type k = q.columns ();
-
-  if (u.numel () != m || v.numel () != n)
-    (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
-
-  ColumnVector utmp = u;
-  ColumnVector vtmp = v;
-  OCTAVE_LOCAL_BUFFER (double, w, 2*k);
-  F77_XFCN (dqr1up, DQR1UP, (m, n, k, q.fortran_vec (),
-                             m, r.fortran_vec (), k,
-                             utmp.fortran_vec (), vtmp.fortran_vec (), w));
-}
-
-void
-QR::update (const Matrix& u, const Matrix& v)
-{
-  octave_idx_type m = q.rows ();
-  octave_idx_type n = r.columns ();
-  octave_idx_type k = q.columns ();
-
-  if (u.rows () != m || v.rows () != n || u.cols () != v.cols ())
-    (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
-
-  OCTAVE_LOCAL_BUFFER (double, w, 2*k);
-  for (volatile octave_idx_type i = 0; i < u.cols (); i++)
-    {
-      ColumnVector utmp = u.column (i);
-      ColumnVector vtmp = v.column (i);
-      F77_XFCN (dqr1up, DQR1UP, (m, n, k, q.fortran_vec (),
-                                 m, r.fortran_vec (), k,
-                                 utmp.fortran_vec (), vtmp.fortran_vec (),
-                                 w));
-    }
-}
-
-void
-QR::insert_col (const ColumnVector& u, octave_idx_type j)
-{
-  octave_idx_type m = q.rows ();
-  octave_idx_type n = r.columns ();
-  octave_idx_type k = q.columns ();
-
-  if (u.numel () != m)
-    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
-  if (j < 0 || j > n)
-    (*current_liboctave_error_handler) ("qrinsert: index out of range");
-
-  if (k < m)
-    {
-      q.resize (m, k+1);
-      r.resize (k+1, n+1);
-    }
-  else
-    {
-      r.resize (k, n+1);
-    }
-
-  ColumnVector utmp = u;
-  OCTAVE_LOCAL_BUFFER (double, w, k);
-  F77_XFCN (dqrinc, DQRINC, (m, n, k, q.fortran_vec (), q.rows (),
-                             r.fortran_vec (), r.rows (), j + 1,
-                             utmp.data (), w));
-}
-
-void
-QR::insert_col (const Matrix& u, const Array<octave_idx_type>& j)
-{
-  octave_idx_type m = q.rows ();
-  octave_idx_type n = r.columns ();
-  octave_idx_type k = q.columns ();
-
-  Array<octave_idx_type> jsi;
-  Array<octave_idx_type> js = j.sort (jsi, 0, ASCENDING);
-  octave_idx_type nj = js.numel ();
-  bool dups = false;
-  for (octave_idx_type i = 0; i < nj - 1; i++)
-    dups = dups && js(i) == js(i+1);
-
-  if (dups)
-    (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
-  if (u.numel () != m || u.columns () != nj)
-    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
-  if (nj > 0 && (js(0) < 0 || js(nj-1) > n))
-    (*current_liboctave_error_handler) ("qrinsert: index out of range");
-
-  if (nj > 0)
-    {
-      octave_idx_type kmax = std::min (k + nj, m);
-      if (k < m)
-        {
-          q.resize (m, kmax);
-          r.resize (kmax, n + nj);
-        }
-      else
-        {
-          r.resize (k, n + nj);
-        }
-
-      OCTAVE_LOCAL_BUFFER (double, w, kmax);
-      for (volatile octave_idx_type i = 0; i < js.numel (); i++)
-        {
-          octave_idx_type ii = i;
-          ColumnVector utmp = u.column (jsi(i));
-          F77_XFCN (dqrinc, DQRINC, (m, n + ii, std::min (kmax, k + ii),
-                                     q.fortran_vec (), q.rows (),
-                                     r.fortran_vec (), r.rows (), js(ii) + 1,
-                                     utmp.data (), w));
-        }
-    }
-}
-
-void
-QR::delete_col (octave_idx_type j)
-{
-  octave_idx_type m = q.rows ();
-  octave_idx_type k = r.rows ();
-  octave_idx_type n = r.columns ();
-
-  if (j < 0 || j > n-1)
-    (*current_liboctave_error_handler) ("qrdelete: index out of range");
-
-  OCTAVE_LOCAL_BUFFER (double, w, k);
-  F77_XFCN (dqrdec, DQRDEC, (m, n, k, q.fortran_vec (), q.rows (),
-                             r.fortran_vec (), r.rows (), j + 1, w));
-
-  if (k < m)
-    {
-      q.resize (m, k-1);
-      r.resize (k-1, n-1);
-    }
-  else
-    {
-      r.resize (k, n-1);
-    }
-}
-
-void
-QR::delete_col (const Array<octave_idx_type>& j)
-{
-  octave_idx_type m = q.rows ();
-  octave_idx_type n = r.columns ();
-  octave_idx_type k = q.columns ();
-
-  Array<octave_idx_type> jsi;
-  Array<octave_idx_type> js = j.sort (jsi, 0, DESCENDING);
-  octave_idx_type nj = js.numel ();
-  bool dups = false;
-  for (octave_idx_type i = 0; i < nj - 1; i++)
-    dups = dups && js(i) == js(i+1);
-
-  if (dups)
-    (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
-  if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0))
-    (*current_liboctave_error_handler) ("qrinsert: index out of range");
-
-  if (nj > 0)
-    {
-      OCTAVE_LOCAL_BUFFER (double, w, k);
-      for (volatile octave_idx_type i = 0; i < js.numel (); i++)
-        {
-          octave_idx_type ii = i;
-          F77_XFCN (dqrdec, DQRDEC, (m, n - ii, k == m ? k : k - ii,
-                                     q.fortran_vec (), q.rows (),
-                                     r.fortran_vec (), r.rows (),
-                                     js(ii) + 1, w));
-        }
-      if (k < m)
-        {
-          q.resize (m, k - nj);
-          r.resize (k - nj, n - nj);
-        }
-      else
-        {
-          r.resize (k, n - nj);
-        }
-
-    }
-}
-
-void
-QR::insert_row (const RowVector& u, octave_idx_type j)
-{
-  octave_idx_type m = r.rows ();
-  octave_idx_type n = r.columns ();
-  octave_idx_type k = std::min (m, n);
-
-  if (! q.is_square () || u.numel () != n)
-    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
-  if (j < 0 || j > m)
-    (*current_liboctave_error_handler) ("qrinsert: index out of range");
-
-  q.resize (m + 1, m + 1);
-  r.resize (m + 1, n);
-  RowVector utmp = u;
-  OCTAVE_LOCAL_BUFFER (double, w, k);
-  F77_XFCN (dqrinr, DQRINR, (m, n, q.fortran_vec (), q.rows (),
-                             r.fortran_vec (), r.rows (),
-                             j + 1, utmp.fortran_vec (), w));
-
-}
-
-void
-QR::delete_row (octave_idx_type j)
-{
-  octave_idx_type m = r.rows ();
-  octave_idx_type n = r.columns ();
-
-  if (! q.is_square ())
-    (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch");
-  if (j < 0 || j > m-1)
-    (*current_liboctave_error_handler) ("qrdelete: index out of range");
-
-  OCTAVE_LOCAL_BUFFER (double, w, 2*m);
-  F77_XFCN (dqrder, DQRDER, (m, n, q.fortran_vec (), q.rows (),
-                             r.fortran_vec (), r.rows (), j + 1,
-                             w));
-
-  q.resize (m - 1, m - 1);
-  r.resize (m - 1, n);
-}
-
-void
-QR::shift_cols (octave_idx_type i, octave_idx_type j)
-{
-  octave_idx_type m = q.rows ();
-  octave_idx_type k = r.rows ();
-  octave_idx_type n = r.columns ();
-
-  if (i < 0 || i > n-1 || j < 0 || j > n-1)
-    (*current_liboctave_error_handler) ("qrshift: index out of range");
-
-  OCTAVE_LOCAL_BUFFER (double, w, 2*k);
-  F77_XFCN (dqrshc, DQRSHC, (m, n, k,
-                             q.fortran_vec (), q.rows (),
-                             r.fortran_vec (), r.rows (),
-                             i + 1, j + 1, w));
-}
-
-#else
-
-// Replacement update methods.
-
-void
-QR::update (const ColumnVector& u, const ColumnVector& v)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type m = q.rows ();
-  octave_idx_type n = r.columns ();
-
-  if (u.numel () != m || v.numel () != n)
-    (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
-
-  init (q*r + Matrix (u) * Matrix (v).transpose (), get_type ());
-}
-
-void
-QR::update (const Matrix& u, const Matrix& v)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type m = q.rows ();
-  octave_idx_type n = r.columns ();
-
-  if (u.rows () != m || v.rows () != n || u.cols () != v.cols ())
-    (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
-
-  init (q*r + u * v.transpose (), get_type ());
-}
-
-static
-Matrix insert_col (const Matrix& a, octave_idx_type i,
-                   const ColumnVector& x)
-{
-  Matrix retval (a.rows (), a.columns () + 1);
-  retval.assign (idx_vector::colon, idx_vector (0, i),
-                 a.index (idx_vector::colon, idx_vector (0, i)));
-  retval.assign (idx_vector::colon, idx_vector (i), x);
-  retval.assign (idx_vector::colon, idx_vector (i+1, retval.columns ()),
-                 a.index (idx_vector::colon, idx_vector (i, a.columns ())));
-  return retval;
-}
-
-static
-Matrix insert_row (const Matrix& a, octave_idx_type i,
-                   const RowVector& x)
-{
-  Matrix retval (a.rows () + 1, a.columns ());
-  retval.assign (idx_vector (0, i), idx_vector::colon,
-                 a.index (idx_vector (0, i), idx_vector::colon));
-  retval.assign (idx_vector (i), idx_vector::colon, x);
-  retval.assign (idx_vector (i+1, retval.rows ()), idx_vector::colon,
-                 a.index (idx_vector (i, a.rows ()), idx_vector::colon));
-  return retval;
-}
-
-static
-Matrix delete_col (const Matrix& a, octave_idx_type i)
-{
-  Matrix retval = a;
-  retval.delete_elements (1, idx_vector (i));
-  return retval;
-}
-
-static
-Matrix delete_row (const Matrix& a, octave_idx_type i)
-{
-  Matrix retval = a;
-  retval.delete_elements (0, idx_vector (i));
-  return retval;
-}
-
-static
-Matrix shift_cols (const Matrix& a,
-                   octave_idx_type i, octave_idx_type j)
-{
-  octave_idx_type n = a.columns ();
-  Array<octave_idx_type> p (dim_vector (n, 1));
-  for (octave_idx_type k = 0; k < n; k++) p(k) = k;
-  if (i < j)
-    {
-      for (octave_idx_type k = i; k < j; k++) p(k) = k+1;
-      p(j) = i;
-    }
-  else if (j < i)
-    {
-      p(j) = i;
-      for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1;
-    }
-
-  return a.index (idx_vector::colon, idx_vector (p));
-}
-
-void
-QR::insert_col (const ColumnVector& u, octave_idx_type j)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type m = q.rows ();
-  octave_idx_type n = r.columns ();
-
-  if (u.numel () != m)
-    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
-  if (j < 0 || j > n)
-    (*current_liboctave_error_handler) ("qrinsert: index out of range");
-
-  init (::insert_col (q*r, j, u), get_type ());
-}
-
-void
-QR::insert_col (const Matrix& u, const Array<octave_idx_type>& j)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type m = q.rows ();
-  octave_idx_type n = r.columns ();
-
-  Array<octave_idx_type> jsi;
-  Array<octave_idx_type> js = j.sort (jsi, 0, ASCENDING);
-  octave_idx_type nj = js.numel ();
-  bool dups = false;
-  for (octave_idx_type i = 0; i < nj - 1; i++)
-    dups = dups && js(i) == js(i+1);
-
-  if (dups)
-    (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
-  if (u.numel () != m || u.columns () != nj)
-    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
-  if (nj > 0 && (js(0) < 0 || js(nj-1) > n))
-    (*current_liboctave_error_handler) ("qrinsert: index out of range");
-
-  if (nj > 0)
-    {
-      Matrix a = q*r;
-      for (octave_idx_type i = 0; i < js.numel (); i++)
-        a = ::insert_col (a, js(i), u.column (i));
-      init (a, get_type ());
-    }
-}
-
-void
-QR::delete_col (octave_idx_type j)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type n = r.columns ();
-
-  if (j < 0 || j > n-1)
-    (*current_liboctave_error_handler) ("qrdelete: index out of range");
-
-  init (::delete_col (q*r, j), get_type ());
-}
-
-void
-QR::delete_col (const Array<octave_idx_type>& j)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type n = r.columns ();
-
-  Array<octave_idx_type> jsi;
-  Array<octave_idx_type> js = j.sort (jsi, 0, DESCENDING);
-  octave_idx_type nj = js.numel ();
-  bool dups = false;
-  for (octave_idx_type i = 0; i < nj - 1; i++)
-    dups = dups && js(i) == js(i+1);
-
-  if (dups)
-    (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
-  if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0))
-    (*current_liboctave_error_handler) ("qrinsert: index out of range");
-
-  if (nj > 0)
-    {
-      Matrix a = q*r;
-      for (octave_idx_type i = 0; i < js.numel (); i++)
-        a = ::delete_col (a, js(i));
-      init (a, get_type ());
-    }
-}
-
-void
-QR::insert_row (const RowVector& u, octave_idx_type j)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type m = r.rows ();
-  octave_idx_type n = r.columns ();
-
-  if (! q.is_square () || u.numel () != n)
-    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
-  if (j < 0 || j > m)
-    (*current_liboctave_error_handler) ("qrinsert: index out of range");
-
-  init (::insert_row (q*r, j, u), get_type ());
-}
-
-void
-QR::delete_row (octave_idx_type j)
-{
-  octave_idx_type m = r.rows ();
-
-  if (! q.is_square ())
-    (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch");
-  if (j < 0 || j > m-1)
-    (*current_liboctave_error_handler) ("qrdelete: index out of range");
-
-  init (::delete_row (q*r, j), get_type ());
-}
-
-void
-QR::shift_cols (octave_idx_type i, octave_idx_type j)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type n = r.columns ();
-
-  if (i < 0 || i > n-1 || j < 0 || j > n-1)
-    (*current_liboctave_error_handler) ("qrshift: index out of range");
-
-  init (::shift_cols (q*r, i, j), get_type ());
-}
-
-#endif
-
-void
-warn_qrupdate_once (void)
-{
-  static bool warned = false;
-
-  if (! warned)
-    {
-      (*current_liboctave_warning_with_id_handler)
-        ("Octave:missing-dependency",
-         "In this version of Octave, QR & Cholesky updating routines "
-         "simply update the matrix and recalculate factorizations. "
-         "To use fast algorithms, link Octave with the qrupdate library. "
-         "See <http://sourceforge.net/projects/qrupdate>.");
-
-      warned = true;
-    }
-}
diff -r 48d82f74243e -r eb1524b07fe3 liboctave/numeric/dbleQR.h
--- a/liboctave/numeric/dbleQR.h	Tue Feb 16 23:34:44 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,84 +0,0 @@
-/*
-
-Copyright (C) 1994-2015 John W. Eaton
-Copyright (C) 2008-2009 Jaroslav Hajek
-
-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/>.
-
-*/
-
-#if ! defined (octave_dbleQR_h)
-#define octave_dbleQR_h 1
-
-#include "octave-config.h"
-
-#include <iosfwd>
-
-#include "dMatrix.h"
-#include "dColVector.h"
-#include "dRowVector.h"
-#include "base-qr.h"
-
-class
-OCTAVE_API
-QR : public base_qr<Matrix>
-{
-public:
-
-  // Import them here to allow the QR:: prefix.
-  typedef qr_type_t type;
-
-  static const type std = qr_type_std;
-  static const type raw = qr_type_raw;
-  static const type economy = qr_type_economy;
-
-  QR (void) : base_qr<Matrix> () { }
-
-  QR (const Matrix&, qr_type_t = qr_type_std);
-
-  QR (const Matrix& qx, const Matrix& rx)
-    : base_qr<Matrix> (qx, rx) { }
-
-  QR (const QR& a) : base_qr<Matrix> (a) { }
-
-  void init (const Matrix&, qr_type_t);
-
-  void update (const ColumnVector& u, const ColumnVector& v);
-
-  void update (const Matrix& u, const Matrix& v);
-
-  void insert_col (const ColumnVector& u, octave_idx_type j);
-
-  void insert_col (const Matrix& u, const Array<octave_idx_type>& j);
-
-  void delete_col (octave_idx_type j);
-
-  void delete_col (const Array<octave_idx_type>& j);
-
-  void insert_row (const RowVector& u, octave_idx_type j);
-
-  void delete_row (octave_idx_type j);
-
-  void shift_cols (octave_idx_type i, octave_idx_type j);
-
-protected:
-
-  void form (octave_idx_type n, Matrix& afact,
-             double *tau, qr_type_t qr_type);
-};
-
-#endif
diff -r 48d82f74243e -r eb1524b07fe3 liboctave/numeric/dbleQRP.cc
--- a/liboctave/numeric/dbleQRP.cc	Tue Feb 16 23:34:44 2016 -0500
+++ b/liboctave/numeric/dbleQRP.cc	Wed Feb 17 02:57:21 2016 -0500
@@ -41,18 +41,18 @@
                              const octave_idx_type&, octave_idx_type&);
 }
 
-// It would be best to share some of this code with QR class...
+// It would be best to share some of this code with qr<Matrix> class...
 
-QRP::QRP (const Matrix& a, qr_type_t qr_type)
-  : QR (), p ()
+QRP::QRP (const Matrix& a, type qr_type)
+  : qr<Matrix> (), p ()
 {
   init (a, qr_type);
 }
 
 void
-QRP::init (const Matrix& a, qr_type_t qr_type)
+QRP::init (const Matrix& a, type qr_type)
 {
-  assert (qr_type != qr_type_raw);
+  assert (qr_type != qr<Matrix>::raw);
 
   octave_idx_type m = a.rows ();
   octave_idx_type n = a.cols ();
@@ -63,7 +63,7 @@
   octave_idx_type info = 0;
 
   Matrix afact = a;
-  if (m > n && qr_type == qr_type_std)
+  if (m > n && qr_type == qr<Matrix>::std)
     afact.resize (m, m);
 
   MArray<octave_idx_type> jpvt (dim_vector (n, 1), 0);
diff -r 48d82f74243e -r eb1524b07fe3 liboctave/numeric/dbleQRP.h
--- a/liboctave/numeric/dbleQRP.h	Tue Feb 16 23:34:44 2016 -0500
+++ b/liboctave/numeric/dbleQRP.h	Wed Feb 17 02:57:21 2016 -0500
@@ -27,27 +27,31 @@
 
 #include <iosfwd>
 
-#include "dbleQR.h"
 #include "PermMatrix.h"
 #include "dColVector.h"
+#include "dMatrix.h"
+#include "dRowVector.h"
+#include "qr.h"
 
 class
 OCTAVE_API
-QRP : public QR
+QRP : public qr<Matrix>
 {
 public:
 
-  QRP (void) : QR (), p () { }
+  typedef qr<Matrix>::type type;
+
+  QRP (void) : qr<Matrix> (), p () { }
 
-  QRP (const Matrix&, qr_type_t = qr_type_std);
+  QRP (const Matrix&, type = qr<Matrix>::std);
 
-  QRP (const QRP& a) : QR (a), p (a.p) { }
+  QRP (const QRP& a) : qr<Matrix> (a), p (a.p) { }
 
   QRP& operator = (const QRP& a)
   {
     if (this != &a)
       {
-        QR::operator = (a);
+        qr<Matrix>::operator = (a);
         p = a.p;
       }
 
@@ -56,7 +60,7 @@
 
   ~QRP (void) { }
 
-  void init (const Matrix&, qr_type_t = qr_type_std);
+  void init (const Matrix&, type = qr<Matrix>::std);
 
   PermMatrix P (void) const { return p; }
 
diff -r 48d82f74243e -r eb1524b07fe3 liboctave/numeric/fCmplxQR.cc
--- a/liboctave/numeric/fCmplxQR.cc	Tue Feb 16 23:34:44 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,694 +0,0 @@
-/*
-
-Copyright (C) 1994-2015 John W. Eaton
-Copyright (C) 2008-2009 Jaroslav Hajek
-Copyright (C) 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 "fCmplxQR.h"
-#include "f77-fcn.h"
-#include "lo-error.h"
-#include "Range.h"
-#include "idx-vector.h"
-#include "oct-locbuf.h"
-
-#include "base-qr.cc"
-
-template class base_qr<FloatComplexMatrix>;
-
-extern "C"
-{
-  F77_RET_T
-  F77_FUNC (cgeqrf, CGEQRF) (const octave_idx_type&, const octave_idx_type&,
-                             FloatComplex*, const octave_idx_type&,
-                             FloatComplex*, FloatComplex*,
-                             const octave_idx_type&, octave_idx_type&);
-
-  F77_RET_T
-  F77_FUNC (cungqr, CUNGQR) (const octave_idx_type&, const octave_idx_type&,
-                             const octave_idx_type&, FloatComplex*,
-                             const octave_idx_type&, FloatComplex*,
-                             FloatComplex*, const octave_idx_type&,
-                             octave_idx_type&);
-
-#ifdef HAVE_QRUPDATE
-
-  F77_RET_T
-  F77_FUNC (cqr1up, CQR1UP) (const octave_idx_type&, const octave_idx_type&,
-                             const octave_idx_type&, FloatComplex*,
-                             const octave_idx_type&, FloatComplex*,
-                             const octave_idx_type&, FloatComplex*,
-                             FloatComplex*, FloatComplex*, float*);
-
-  F77_RET_T
-  F77_FUNC (cqrinc, CQRINC) (const octave_idx_type&, const octave_idx_type&,
-                             const octave_idx_type&, FloatComplex*,
-                             const octave_idx_type&, FloatComplex*,
-                             const octave_idx_type&,const octave_idx_type&,
-                             const FloatComplex*, float*);
-
-  F77_RET_T
-  F77_FUNC (cqrdec, CQRDEC) (const octave_idx_type&, const octave_idx_type&,
-                             const octave_idx_type&, FloatComplex*,
-                             const octave_idx_type&, FloatComplex*,
-                             const octave_idx_type&, const octave_idx_type&,
-                             float*);
-
-  F77_RET_T
-  F77_FUNC (cqrinr, CQRINR) (const octave_idx_type&, const octave_idx_type&,
-                             FloatComplex*, const octave_idx_type&,
-                             FloatComplex*, const octave_idx_type&,
-                             const octave_idx_type&, const FloatComplex*,
-                             float*);
-
-  F77_RET_T
-  F77_FUNC (cqrder, CQRDER) (const octave_idx_type&, const octave_idx_type&,
-                             FloatComplex*, const octave_idx_type&,
-                             FloatComplex*, const octave_idx_type&,
-                             const octave_idx_type&, FloatComplex*, float*);
-
-  F77_RET_T
-  F77_FUNC (cqrshc, CQRSHC) (const octave_idx_type&, const octave_idx_type&,
-                             const octave_idx_type&, FloatComplex*,
-                             const octave_idx_type&, FloatComplex*,
-                             const octave_idx_type&, const octave_idx_type&,
-                             const octave_idx_type&, FloatComplex*,
-                             float*);
-
-#endif
-}
-
-FloatComplexQR::FloatComplexQR (const FloatComplexMatrix& a, qr_type_t qr_type)
-{
-  init (a, qr_type);
-}
-
-void
-FloatComplexQR::init (const FloatComplexMatrix& a, qr_type_t qr_type)
-{
-  octave_idx_type m = a.rows ();
-  octave_idx_type n = a.cols ();
-
-  octave_idx_type min_mn = m < n ? m : n;
-  OCTAVE_LOCAL_BUFFER (FloatComplex, tau, min_mn);
-
-  octave_idx_type info = 0;
-
-  FloatComplexMatrix afact = a;
-  if (m > n && qr_type == qr_type_std)
-    afact.resize (m, m);
-
-  if (m > 0)
-    {
-      // workspace query.
-      FloatComplex clwork;
-      F77_XFCN (cgeqrf, CGEQRF, (m, n, afact.fortran_vec (), m, tau,
-                                 &clwork, -1, info));
-
-      // allocate buffer and do the job.
-      octave_idx_type lwork = clwork.real ();
-      lwork = std::max (lwork, static_cast<octave_idx_type> (1));
-      OCTAVE_LOCAL_BUFFER (FloatComplex, work, lwork);
-      F77_XFCN (cgeqrf, CGEQRF, (m, n, afact.fortran_vec (), m, tau,
-                                 work, lwork, info));
-    }
-
-  form (n, afact, tau, qr_type);
-}
-
-void FloatComplexQR::form (octave_idx_type n, FloatComplexMatrix& afact,
-                           FloatComplex *tau, qr_type_t qr_type)
-{
-  octave_idx_type m = afact.rows ();
-  octave_idx_type min_mn = std::min (m, n);
-  octave_idx_type info;
-
-  if (qr_type == qr_type_raw)
-    {
-      for (octave_idx_type j = 0; j < min_mn; j++)
-        {
-          octave_idx_type limit = j < min_mn - 1 ? j : min_mn - 1;
-          for (octave_idx_type i = limit + 1; i < m; i++)
-            afact.elem (i, j) *= tau[j];
-        }
-
-      r = afact;
-    }
-  else
-    {
-      // Attempt to minimize copying.
-      if (m >= n)
-        {
-          // afact will become q.
-          q = afact;
-          octave_idx_type k = qr_type == qr_type_economy ? n : m;
-          r = FloatComplexMatrix (k, n);
-          for (octave_idx_type j = 0; j < n; j++)
-            {
-              octave_idx_type i = 0;
-              for (; i <= j; i++)
-                r.xelem (i, j) = afact.xelem (i, j);
-              for (; i < k; i++)
-                r.xelem (i, j) = 0;
-            }
-          afact = FloatComplexMatrix (); // optimize memory
-        }
-      else
-        {
-          // afact will become r.
-          q = FloatComplexMatrix (m, m);
-          for (octave_idx_type j = 0; j < m; j++)
-            for (octave_idx_type i = j + 1; i < m; i++)
-              {
-                q.xelem (i, j) = afact.xelem (i, j);
-                afact.xelem (i, j) = 0;
-              }
-          r = afact;
-        }
-
-
-      if (m > 0)
-        {
-          octave_idx_type k = q.columns ();
-          // workspace query.
-          FloatComplex clwork;
-          F77_XFCN (cungqr, CUNGQR, (m, k, min_mn, q.fortran_vec (), m, tau,
-                                     &clwork, -1, info));
-
-          // allocate buffer and do the job.
-          octave_idx_type lwork = clwork.real ();
-          lwork = std::max (lwork, static_cast<octave_idx_type> (1));
-          OCTAVE_LOCAL_BUFFER (FloatComplex, work, lwork);
-          F77_XFCN (cungqr, CUNGQR, (m, k, min_mn, q.fortran_vec (), m, tau,
-                                     work, lwork, info));
-        }
-    }
-}
-
-#ifdef HAVE_QRUPDATE
-
-void
-FloatComplexQR::update (const FloatComplexColumnVector& u,
-                        const FloatComplexColumnVector& v)
-{
-  octave_idx_type m = q.rows ();
-  octave_idx_type n = r.columns ();
-  octave_idx_type k = q.columns ();
-
-  if (u.numel () != m || v.numel () != n)
-    (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
-
-  FloatComplexColumnVector utmp = u;
-  FloatComplexColumnVector vtmp = v;
-  OCTAVE_LOCAL_BUFFER (FloatComplex, w, k);
-  OCTAVE_LOCAL_BUFFER (float, rw, k);
-  F77_XFCN (cqr1up, CQR1UP, (m, n, k, q.fortran_vec (),
-                             m, r.fortran_vec (), k,
-                             utmp.fortran_vec (), vtmp.fortran_vec (),
-                             w, rw));
-}
-
-void
-FloatComplexQR::update (const FloatComplexMatrix& u,
-                        const FloatComplexMatrix& v)
-{
-  octave_idx_type m = q.rows ();
-  octave_idx_type n = r.columns ();
-  octave_idx_type k = q.columns ();
-
-  if (u.rows () != m || v.rows () != n || u.cols () != v.cols ())
-    (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
-
-  OCTAVE_LOCAL_BUFFER (FloatComplex, w, k);
-  OCTAVE_LOCAL_BUFFER (float, rw, k);
-  for (volatile octave_idx_type i = 0; i < u.cols (); i++)
-    {
-      FloatComplexColumnVector utmp = u.column (i);
-      FloatComplexColumnVector vtmp = v.column (i);
-      F77_XFCN (cqr1up, CQR1UP, (m, n, k, q.fortran_vec (),
-                                 m, r.fortran_vec (), k,
-                                 utmp.fortran_vec (), vtmp.fortran_vec (),
-                                 w, rw));
-    }
-}
-
-void
-FloatComplexQR::insert_col (const FloatComplexColumnVector& u,
-                            octave_idx_type j)
-{
-  octave_idx_type m = q.rows ();
-  octave_idx_type n = r.columns ();
-  octave_idx_type k = q.columns ();
-
-  if (u.numel () != m)
-    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
-  if (j < 0 || j > n)
-    (*current_liboctave_error_handler) ("qrinsert: index out of range");
-
-  if (k < m)
-    {
-      q.resize (m, k+1);
-      r.resize (k+1, n+1);
-    }
-  else
-    {
-      r.resize (k, n+1);
-    }
-
-  FloatComplexColumnVector utmp = u;
-  OCTAVE_LOCAL_BUFFER (float, rw, k);
-  F77_XFCN (cqrinc, CQRINC, (m, n, k, q.fortran_vec (), q.rows (),
-                             r.fortran_vec (), r.rows (), j + 1,
-                             utmp.data (), rw));
-}
-
-void
-FloatComplexQR::insert_col (const FloatComplexMatrix& u,
-                            const Array<octave_idx_type>& j)
-{
-  octave_idx_type m = q.rows ();
-  octave_idx_type n = r.columns ();
-  octave_idx_type k = q.columns ();
-
-  Array<octave_idx_type> jsi;
-  Array<octave_idx_type> js = j.sort (jsi, 0, ASCENDING);
-  octave_idx_type nj = js.numel ();
-  bool dups = false;
-  for (octave_idx_type i = 0; i < nj - 1; i++)
-    dups = dups && js(i) == js(i+1);
-
-  if (dups)
-    (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
-  if (u.numel () != m || u.columns () != nj)
-    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
-  if (nj > 0 && (js(0) < 0 || js(nj-1) > n))
-    (*current_liboctave_error_handler) ("qrinsert: index out of range");
-
-  if (nj > 0)
-    {
-      octave_idx_type kmax = std::min (k + nj, m);
-      if (k < m)
-        {
-          q.resize (m, kmax);
-          r.resize (kmax, n + nj);
-        }
-      else
-        {
-          r.resize (k, n + nj);
-        }
-
-      OCTAVE_LOCAL_BUFFER (float, rw, kmax);
-      for (volatile octave_idx_type i = 0; i < js.numel (); i++)
-        {
-          octave_idx_type ii = i;
-          F77_XFCN (cqrinc, CQRINC, (m, n + ii, std::min (kmax, k + ii),
-                                     q.fortran_vec (), q.rows (),
-                                     r.fortran_vec (), r.rows (), js(ii) + 1,
-                                     u.column (jsi(i)).data (), rw));
-        }
-    }
-}
-
-void
-FloatComplexQR::delete_col (octave_idx_type j)
-{
-  octave_idx_type m = q.rows ();
-  octave_idx_type k = r.rows ();
-  octave_idx_type n = r.columns ();
-
-  if (j < 0 || j > n-1)
-    (*current_liboctave_error_handler) ("qrdelete: index out of range");
-
-  OCTAVE_LOCAL_BUFFER (float, rw, k);
-  F77_XFCN (cqrdec, CQRDEC, (m, n, k, q.fortran_vec (), q.rows (),
-                             r.fortran_vec (), r.rows (), j + 1, rw));
-
-  if (k < m)
-    {
-      q.resize (m, k-1);
-      r.resize (k-1, n-1);
-    }
-  else
-    {
-      r.resize (k, n-1);
-    }
-}
-
-void
-FloatComplexQR::delete_col (const Array<octave_idx_type>& j)
-{
-  octave_idx_type m = q.rows ();
-  octave_idx_type n = r.columns ();
-  octave_idx_type k = q.columns ();
-
-  Array<octave_idx_type> jsi;
-  Array<octave_idx_type> js = j.sort (jsi, 0, DESCENDING);
-  octave_idx_type nj = js.numel ();
-  bool dups = false;
-  for (octave_idx_type i = 0; i < nj - 1; i++)
-    dups = dups && js(i) == js(i+1);
-
-  if (dups)
-    (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
-  if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0))
-    (*current_liboctave_error_handler) ("qrinsert: index out of range");
-
-  if (nj > 0)
-    {
-      OCTAVE_LOCAL_BUFFER (float, rw, k);
-      for (volatile octave_idx_type i = 0; i < js.numel (); i++)
-        {
-          octave_idx_type ii = i;
-          F77_XFCN (cqrdec, CQRDEC, (m, n - ii, k == m ? k : k - ii,
-                                     q.fortran_vec (), q.rows (),
-                                     r.fortran_vec (), r.rows (),
-                                     js(ii) + 1, rw));
-        }
-      if (k < m)
-        {
-          q.resize (m, k - nj);
-          r.resize (k - nj, n - nj);
-        }
-      else
-        {
-          r.resize (k, n - nj);
-        }
-
-    }
-}
-
-void
-FloatComplexQR::insert_row (const FloatComplexRowVector& u, octave_idx_type j)
-{
-  octave_idx_type m = r.rows ();
-  octave_idx_type n = r.columns ();
-  octave_idx_type k = std::min (m, n);
-
-  if (! q.is_square () || u.numel () != n)
-    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
-  if (j < 0 || j > m)
-    (*current_liboctave_error_handler) ("qrinsert: index out of range");
-
-  q.resize (m + 1, m + 1);
-  r.resize (m + 1, n);
-  FloatComplexRowVector utmp = u;
-  OCTAVE_LOCAL_BUFFER (float, rw, k);
-  F77_XFCN (cqrinr, CQRINR, (m, n, q.fortran_vec (), q.rows (),
-                             r.fortran_vec (), r.rows (),
-                             j + 1, utmp.fortran_vec (), rw));
-
-}
-
-void
-FloatComplexQR::delete_row (octave_idx_type j)
-{
-  octave_idx_type m = r.rows ();
-  octave_idx_type n = r.columns ();
-
-  if (! q.is_square ())
-    (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch");
-  if (j < 0 || j > m-1)
-    (*current_liboctave_error_handler) ("qrdelete: index out of range");
-
-  OCTAVE_LOCAL_BUFFER (FloatComplex, w, m);
-  OCTAVE_LOCAL_BUFFER (float, rw, m);
-  F77_XFCN (cqrder, CQRDER, (m, n, q.fortran_vec (), q.rows (),
-                             r.fortran_vec (), r.rows (), j + 1,
-                             w, rw));
-
-  q.resize (m - 1, m - 1);
-  r.resize (m - 1, n);
-}
-
-void
-FloatComplexQR::shift_cols (octave_idx_type i, octave_idx_type j)
-{
-  octave_idx_type m = q.rows ();
-  octave_idx_type k = r.rows ();
-  octave_idx_type n = r.columns ();
-
-  if (i < 0 || i > n-1 || j < 0 || j > n-1)
-    (*current_liboctave_error_handler) ("qrshift: index out of range");
-
-  OCTAVE_LOCAL_BUFFER (FloatComplex, w, k);
-  OCTAVE_LOCAL_BUFFER (float, rw, k);
-  F77_XFCN (cqrshc, CQRSHC, (m, n, k,
-                             q.fortran_vec (), q.rows (),
-                             r.fortran_vec (), r.rows (),
-                             i + 1, j + 1, w, rw));
-}
-
-#else
-
-// Replacement update methods.
-
-void
-FloatComplexQR::update (const FloatComplexColumnVector& u,
-                        const FloatComplexColumnVector& v)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type m = q.rows ();
-  octave_idx_type n = r.columns ();
-
-  if (u.numel () != m || v.numel () != n)
-    (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
-
-  init (q*r + FloatComplexMatrix (u) * FloatComplexMatrix (v).hermitian (),
-        get_type ());
-}
-
-void
-FloatComplexQR::update (const FloatComplexMatrix& u,
-                        const FloatComplexMatrix& v)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type m = q.rows ();
-  octave_idx_type n = r.columns ();
-
-  if (u.rows () != m || v.rows () != n || u.cols () != v.cols ())
-    (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
-
-  init (q*r + u * v.hermitian (), get_type ());
-}
-
-static
-FloatComplexMatrix insert_col (const FloatComplexMatrix& a, octave_idx_type i,
-                               const FloatComplexColumnVector& x)
-{
-  FloatComplexMatrix retval (a.rows (), a.columns () + 1);
-  retval.assign (idx_vector::colon, idx_vector (0, i),
-                 a.index (idx_vector::colon, idx_vector (0, i)));
-  retval.assign (idx_vector::colon, idx_vector (i), x);
-  retval.assign (idx_vector::colon, idx_vector (i+1, retval.columns ()),
-                 a.index (idx_vector::colon, idx_vector (i, a.columns ())));
-  return retval;
-}
-
-static
-FloatComplexMatrix insert_row (const FloatComplexMatrix& a, octave_idx_type i,
-                               const FloatComplexRowVector& x)
-{
-  FloatComplexMatrix retval (a.rows () + 1, a.columns ());
-  retval.assign (idx_vector (0, i), idx_vector::colon,
-                 a.index (idx_vector (0, i), idx_vector::colon));
-  retval.assign (idx_vector (i), idx_vector::colon, x);
-  retval.assign (idx_vector (i+1, retval.rows ()), idx_vector::colon,
-                 a.index (idx_vector (i, a.rows ()), idx_vector::colon));
-  return retval;
-}
-
-static
-FloatComplexMatrix delete_col (const FloatComplexMatrix& a, octave_idx_type i)
-{
-  FloatComplexMatrix retval = a;
-  retval.delete_elements (1, idx_vector (i));
-  return retval;
-}
-
-static
-FloatComplexMatrix delete_row (const FloatComplexMatrix& a, octave_idx_type i)
-{
-  FloatComplexMatrix retval = a;
-  retval.delete_elements (0, idx_vector (i));
-  return retval;
-}
-
-static
-FloatComplexMatrix shift_cols (const FloatComplexMatrix& a,
-                               octave_idx_type i, octave_idx_type j)
-{
-  octave_idx_type n = a.columns ();
-  Array<octave_idx_type> p (dim_vector (n, 1));
-  for (octave_idx_type k = 0; k < n; k++) p(k) = k;
-  if (i < j)
-    {
-      for (octave_idx_type k = i; k < j; k++) p(k) = k+1;
-      p(j) = i;
-    }
-  else if (j < i)
-    {
-      p(j) = i;
-      for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1;
-    }
-
-  return a.index (idx_vector::colon, idx_vector (p));
-}
-
-void
-FloatComplexQR::insert_col (const FloatComplexColumnVector& u,
-                            octave_idx_type j)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type m = q.rows ();
-  octave_idx_type n = r.columns ();
-
-  if (u.numel () != m)
-    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
-  if (j < 0 || j > n)
-    (*current_liboctave_error_handler) ("qrinsert: index out of range");
-
-  init (::insert_col (q*r, j, u), get_type ());
-}
-
-void
-FloatComplexQR::insert_col (const FloatComplexMatrix& u,
-                            const Array<octave_idx_type>& j)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type m = q.rows ();
-  octave_idx_type n = r.columns ();
-
-  Array<octave_idx_type> jsi;
-  Array<octave_idx_type> js = j.sort (jsi, 0, ASCENDING);
-  octave_idx_type nj = js.numel ();
-  bool dups = false;
-  for (octave_idx_type i = 0; i < nj - 1; i++)
-    dups = dups && js(i) == js(i+1);
-
-  if (dups)
-    (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
-  if (u.numel () != m || u.columns () != nj)
-    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
-  if (nj > 0 && (js(0) < 0 || js(nj-1) > n))
-    (*current_liboctave_error_handler) ("qrinsert: index out of range");
-
-  if (nj > 0)
-    {
-      FloatComplexMatrix a = q*r;
-      for (octave_idx_type i = 0; i < js.numel (); i++)
-        a = ::insert_col (a, js(i), u.column (i));
-      init (a, get_type ());
-    }
-}
-
-void
-FloatComplexQR::delete_col (octave_idx_type j)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type n = r.columns ();
-
-  if (j < 0 || j > n-1)
-    (*current_liboctave_error_handler) ("qrdelete: index out of range");
-
-  init (::delete_col (q*r, j), get_type ());
-}
-
-void
-FloatComplexQR::delete_col (const Array<octave_idx_type>& j)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type n = r.columns ();
-
-  Array<octave_idx_type> jsi;
-  Array<octave_idx_type> js = j.sort (jsi, 0, DESCENDING);
-  octave_idx_type nj = js.numel ();
-  bool dups = false;
-  for (octave_idx_type i = 0; i < nj - 1; i++)
-    dups = dups && js(i) == js(i+1);
-
-  if (dups)
-    (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
-  if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0))
-    (*current_liboctave_error_handler) ("qrinsert: index out of range");
-
-  if (nj > 0)
-    {
-      FloatComplexMatrix a = q*r;
-      for (octave_idx_type i = 0; i < js.numel (); i++)
-        a = ::delete_col (a, js(i));
-      init (a, get_type ());
-    }
-}
-
-void
-FloatComplexQR::insert_row (const FloatComplexRowVector& u, octave_idx_type j)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type m = r.rows ();
-  octave_idx_type n = r.columns ();
-
-  if (! q.is_square () || u.numel () != n)
-    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
-  if (j < 0 || j > m)
-    (*current_liboctave_error_handler) ("qrinsert: index out of range");
-
-  init (::insert_row (q*r, j, u), get_type ());
-}
-
-void
-FloatComplexQR::delete_row (octave_idx_type j)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type m = r.rows ();
-
-  if (! q.is_square ())
-    (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch");
-  if (j < 0 || j > m-1)
-    (*current_liboctave_error_handler) ("qrdelete: index out of range");
-
-  init (::delete_row (q*r, j), get_type ());
-}
-
-void
-FloatComplexQR::shift_cols (octave_idx_type i, octave_idx_type j)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type n = r.columns ();
-
-  if (i < 0 || i > n-1 || j < 0 || j > n-1)
-    (*current_liboctave_error_handler) ("qrshift: index out of range");
-
-  init (::shift_cols (q*r, i, j), get_type ());
-}
-
-#endif
diff -r 48d82f74243e -r eb1524b07fe3 liboctave/numeric/fCmplxQR.h
--- a/liboctave/numeric/fCmplxQR.h	Tue Feb 16 23:34:44 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,81 +0,0 @@
-/*
-
-Copyright (C) 1994-2015 John W. Eaton
-Copyright (C) 2008-2009 Jaroslav Hajek
-
-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/>.
-
-*/
-
-// updating/downdating by Jaroslav Hajek 2008
-
-#if ! defined (octave_fCmplxQR_h)
-#define octave_fCmplxQR_h 1
-
-#include "octave-config.h"
-
-#include <iosfwd>
-
-#include "fCMatrix.h"
-#include "fCColVector.h"
-#include "fCRowVector.h"
-#include "base-qr.h"
-
-class
-OCTAVE_API
-FloatComplexQR : public base_qr<FloatComplexMatrix>
-{
-public:
-
-  FloatComplexQR (void) : base_qr<FloatComplexMatrix> () { }
-
-  FloatComplexQR (const FloatComplexMatrix&, qr_type_t = qr_type_std);
-
-  FloatComplexQR (const FloatComplexMatrix& qx, const FloatComplexMatrix& rx)
-    : base_qr<FloatComplexMatrix> (qx, rx) { }
-
-  FloatComplexQR (const FloatComplexQR& a) : base_qr<FloatComplexMatrix> (a) { }
-
-  void init (const FloatComplexMatrix&, qr_type_t = qr_type_std);
-
-  void update (const FloatComplexColumnVector& u,
-               const FloatComplexColumnVector& v);
-
-  void update (const FloatComplexMatrix& u, const FloatComplexMatrix& v);
-
-  void insert_col (const FloatComplexColumnVector& u, octave_idx_type j);
-
-  void insert_col (const FloatComplexMatrix& u,
-                   const Array<octave_idx_type>& j);
-
-  void delete_col (octave_idx_type j);
-
-  void delete_col (const Array<octave_idx_type>& j);
-
-  void insert_row (const FloatComplexRowVector& u, octave_idx_type j);
-
-  void delete_row (octave_idx_type j);
-
-  void shift_cols (octave_idx_type i, octave_idx_type j);
-
-protected:
-
-  void form (octave_idx_type n, FloatComplexMatrix& afact,
-             FloatComplex *tau, qr_type_t qr_type);
-};
-
-#endif
diff -r 48d82f74243e -r eb1524b07fe3 liboctave/numeric/fCmplxQRP.cc
--- a/liboctave/numeric/fCmplxQRP.cc	Tue Feb 16 23:34:44 2016 -0500
+++ b/liboctave/numeric/fCmplxQRP.cc	Wed Feb 17 02:57:21 2016 -0500
@@ -41,19 +41,18 @@
                              const octave_idx_type&, float*, octave_idx_type&);
 }
 
-// It would be best to share some of this code with FloatComplexQR class...
+// It would be best to share some of this code with qr<FloatComplexMatrix> class...
 
-FloatComplexQRP::FloatComplexQRP (const FloatComplexMatrix& a,
-                                  qr_type_t qr_type)
-  : FloatComplexQR (), p ()
+FloatComplexQRP::FloatComplexQRP (const FloatComplexMatrix& a, type qr_type)
+  : qr<FloatComplexMatrix> (), p ()
 {
   init (a, qr_type);
 }
 
 void
-FloatComplexQRP::init (const FloatComplexMatrix& a, qr_type_t qr_type)
+FloatComplexQRP::init (const FloatComplexMatrix& a, type qr_type)
 {
-  assert (qr_type != qr_type_raw);
+  assert (qr_type != qr<FloatComplexMatrix>::raw);
 
   octave_idx_type m = a.rows ();
   octave_idx_type n = a.cols ();
@@ -64,7 +63,7 @@
   octave_idx_type info = 0;
 
   FloatComplexMatrix afact = a;
-  if (m > n && qr_type == qr_type_std)
+  if (m > n && qr_type == qr<FloatComplexMatrix>::std)
     afact.resize (m, m);
 
   MArray<octave_idx_type> jpvt (dim_vector (n, 1), 0);
diff -r 48d82f74243e -r eb1524b07fe3 liboctave/numeric/fCmplxQRP.h
--- a/liboctave/numeric/fCmplxQRP.h	Tue Feb 16 23:34:44 2016 -0500
+++ b/liboctave/numeric/fCmplxQRP.h	Wed Feb 17 02:57:21 2016 -0500
@@ -27,27 +27,31 @@
 
 #include <iosfwd>
 
-#include "fCmplxQR.h"
 #include "PermMatrix.h"
-#include "fColVector.h"
+#include "fCColVector.h"
+#include "fCMatrix.h"
+#include "fCRowVector.h"
+#include "qr.h"
 
 class
 OCTAVE_API
-FloatComplexQRP : public FloatComplexQR
+FloatComplexQRP : public qr<FloatComplexMatrix>
 {
 public:
 
-  FloatComplexQRP (void) : FloatComplexQR (), p () { }
+  typedef qr<FloatComplexMatrix>::type type;
+
+  FloatComplexQRP (void) : qr<FloatComplexMatrix> (), p () { }
 
-  FloatComplexQRP (const FloatComplexMatrix&, qr_type_t = qr_type_std);
+  FloatComplexQRP (const FloatComplexMatrix&, type = qr<FloatComplexMatrix>::std);
 
-  FloatComplexQRP (const FloatComplexQRP& a) : FloatComplexQR (a), p (a.p) { }
+  FloatComplexQRP (const FloatComplexQRP& a) : qr<FloatComplexMatrix> (a), p (a.p) { }
 
   FloatComplexQRP& operator = (const FloatComplexQRP& a)
   {
     if (this != &a)
       {
-        FloatComplexQR::operator = (a);
+        qr<FloatComplexMatrix>::operator = (a);
         p = a.p;
       }
     return *this;
@@ -55,7 +59,7 @@
 
   ~FloatComplexQRP (void) { }
 
-  void init (const FloatComplexMatrix&, qr_type_t = qr_type_std);
+  void init (const FloatComplexMatrix&, type = qr<FloatComplexMatrix>::std);
 
   PermMatrix P (void) const { return p; }
 
diff -r 48d82f74243e -r eb1524b07fe3 liboctave/numeric/floatQR.cc
--- a/liboctave/numeric/floatQR.cc	Tue Feb 16 23:34:44 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,677 +0,0 @@
-/*
-
-Copyright (C) 1994-2015 John W. Eaton
-Copyright (C) 2008-2009 Jaroslav Hajek
-Copyright (C) 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 "floatQR.h"
-#include "f77-fcn.h"
-#include "lo-error.h"
-#include "Range.h"
-#include "idx-vector.h"
-#include "oct-locbuf.h"
-
-#include "base-qr.cc"
-
-template class base_qr<FloatMatrix>;
-
-extern "C"
-{
-  F77_RET_T
-  F77_FUNC (sgeqrf, SGEQRF) (const octave_idx_type&, const octave_idx_type&,
-                             float*, const octave_idx_type&, float*, float*,
-                             const octave_idx_type&, octave_idx_type&);
-
-  F77_RET_T
-  F77_FUNC (sorgqr, SORGQR) (const octave_idx_type&, const octave_idx_type&,
-                             const octave_idx_type&, float*,
-                             const octave_idx_type&, float*, float*,
-                             const octave_idx_type&, octave_idx_type&);
-
-#ifdef HAVE_QRUPDATE
-
-  F77_RET_T
-  F77_FUNC (sqr1up, SQR1UP) (const octave_idx_type&, const octave_idx_type&,
-                             const octave_idx_type&, float*,
-                             const octave_idx_type&, float*,
-                             const octave_idx_type&, float*, float*, float*);
-
-  F77_RET_T
-  F77_FUNC (sqrinc, SQRINC) (const octave_idx_type&, const octave_idx_type&,
-                             const octave_idx_type&, float*,
-                             const octave_idx_type&, float*,
-                             const octave_idx_type&,
-                             const octave_idx_type&, const float*, float*);
-
-  F77_RET_T
-  F77_FUNC (sqrdec, SQRDEC) (const octave_idx_type&, const octave_idx_type&,
-                             const octave_idx_type&, float*,
-                             const octave_idx_type&, float*,
-                             const octave_idx_type&,
-                             const octave_idx_type&, float*);
-
-  F77_RET_T
-  F77_FUNC (sqrinr, SQRINR) (const octave_idx_type&, const octave_idx_type&,
-                             float*, const octave_idx_type&,
-                             float*, const octave_idx_type&,
-                             const octave_idx_type&, const float*, float*);
-
-  F77_RET_T
-  F77_FUNC (sqrder, SQRDER) (const octave_idx_type&, const octave_idx_type&,
-                             float*, const octave_idx_type&,
-                             float*, const octave_idx_type&,
-                             const octave_idx_type&, float*);
-
-  F77_RET_T
-  F77_FUNC (sqrshc, SQRSHC) (const octave_idx_type&, const octave_idx_type&,
-                             const octave_idx_type&, float*,
-                             const octave_idx_type&, float*,
-                             const octave_idx_type&,
-                             const octave_idx_type&, const octave_idx_type&,
-                             float*);
-
-#endif
-}
-
-FloatQR::FloatQR (const FloatMatrix& a, qr_type_t qr_type)
-{
-  init (a, qr_type);
-}
-
-void
-FloatQR::init (const FloatMatrix& a, qr_type_t qr_type)
-{
-  octave_idx_type m = a.rows ();
-  octave_idx_type n = a.cols ();
-
-  octave_idx_type min_mn = m < n ? m : n;
-  OCTAVE_LOCAL_BUFFER (float, tau, min_mn);
-
-  octave_idx_type info = 0;
-
-  FloatMatrix afact = a;
-  if (m > n && qr_type == qr_type_std)
-    afact.resize (m, m);
-
-  if (m > 0)
-    {
-      // workspace query.
-      float rlwork;
-      F77_XFCN (sgeqrf, SGEQRF, (m, n, afact.fortran_vec (), m, tau,
-                                 &rlwork, -1, info));
-
-      // allocate buffer and do the job.
-      octave_idx_type lwork = rlwork;
-      lwork = std::max (lwork, static_cast<octave_idx_type> (1));
-      OCTAVE_LOCAL_BUFFER (float, work, lwork);
-      F77_XFCN (sgeqrf, SGEQRF, (m, n, afact.fortran_vec (), m, tau,
-                                 work, lwork, info));
-    }
-
-  form (n, afact, tau, qr_type);
-}
-
-void FloatQR::form (octave_idx_type n, FloatMatrix& afact,
-                    float *tau, qr_type_t qr_type)
-{
-  octave_idx_type m = afact.rows ();
-  octave_idx_type min_mn = std::min (m, n);
-  octave_idx_type info;
-
-  if (qr_type == qr_type_raw)
-    {
-      for (octave_idx_type j = 0; j < min_mn; j++)
-        {
-          octave_idx_type limit = j < min_mn - 1 ? j : min_mn - 1;
-          for (octave_idx_type i = limit + 1; i < m; i++)
-            afact.elem (i, j) *= tau[j];
-        }
-
-      r = afact;
-    }
-  else
-    {
-      // Attempt to minimize copying.
-      if (m >= n)
-        {
-          // afact will become q.
-          q = afact;
-          octave_idx_type k = qr_type == qr_type_economy ? n : m;
-          r = FloatMatrix (k, n);
-          for (octave_idx_type j = 0; j < n; j++)
-            {
-              octave_idx_type i = 0;
-              for (; i <= j; i++)
-                r.xelem (i, j) = afact.xelem (i, j);
-              for (; i < k; i++)
-                r.xelem (i, j) = 0;
-            }
-          afact = FloatMatrix (); // optimize memory
-        }
-      else
-        {
-          // afact will become r.
-          q = FloatMatrix (m, m);
-          for (octave_idx_type j = 0; j < m; j++)
-            for (octave_idx_type i = j + 1; i < m; i++)
-              {
-                q.xelem (i, j) = afact.xelem (i, j);
-                afact.xelem (i, j) = 0;
-              }
-          r = afact;
-        }
-
-
-      if (m > 0)
-        {
-          octave_idx_type k = q.columns ();
-          // workspace query.
-          float rlwork;
-          F77_XFCN (sorgqr, SORGQR, (m, k, min_mn, q.fortran_vec (), m, tau,
-                                     &rlwork, -1, info));
-
-          // allocate buffer and do the job.
-          octave_idx_type lwork = rlwork;
-          lwork = std::max (lwork, static_cast<octave_idx_type> (1));
-          OCTAVE_LOCAL_BUFFER (float, work, lwork);
-          F77_XFCN (sorgqr, SORGQR, (m, k, min_mn, q.fortran_vec (), m, tau,
-                                     work, lwork, info));
-        }
-    }
-}
-
-#ifdef HAVE_QRUPDATE
-
-void
-FloatQR::update (const FloatColumnVector& u, const FloatColumnVector& v)
-{
-  octave_idx_type m = q.rows ();
-  octave_idx_type n = r.columns ();
-  octave_idx_type k = q.columns ();
-
-  if (u.numel () != m || v.numel () != n)
-    (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
-
-  FloatColumnVector utmp = u;
-  FloatColumnVector vtmp = v;
-  OCTAVE_LOCAL_BUFFER (float, w, 2*k);
-  F77_XFCN (sqr1up, SQR1UP, (m, n, k, q.fortran_vec (),
-                             m, r.fortran_vec (), k,
-                             utmp.fortran_vec (), vtmp.fortran_vec (), w));
-}
-
-void
-FloatQR::update (const FloatMatrix& u, const FloatMatrix& v)
-{
-  octave_idx_type m = q.rows ();
-  octave_idx_type n = r.columns ();
-  octave_idx_type k = q.columns ();
-
-  if (u.rows () != m || v.rows () != n || u.cols () != v.cols ())
-    (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
-
-  OCTAVE_LOCAL_BUFFER (float, w, 2*k);
-  for (volatile octave_idx_type i = 0; i < u.cols (); i++)
-    {
-      FloatColumnVector utmp = u.column (i);
-      FloatColumnVector vtmp = v.column (i);
-      F77_XFCN (sqr1up, SQR1UP, (m, n, k, q.fortran_vec (),
-                                 m, r.fortran_vec (), k,
-                                 utmp.fortran_vec (), vtmp.fortran_vec (),
-                                 w));
-    }
-}
-
-void
-FloatQR::insert_col (const FloatColumnVector& u, octave_idx_type j)
-{
-  octave_idx_type m = q.rows ();
-  octave_idx_type n = r.columns ();
-  octave_idx_type k = q.columns ();
-
-  if (u.numel () != m)
-    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
-  if (j < 0 || j > n)
-    (*current_liboctave_error_handler) ("qrinsert: index out of range");
-
-  if (k < m)
-    {
-      q.resize (m, k+1);
-      r.resize (k+1, n+1);
-    }
-  else
-    {
-      r.resize (k, n+1);
-    }
-
-  FloatColumnVector utmp = u;
-  OCTAVE_LOCAL_BUFFER (float, w, k);
-  F77_XFCN (sqrinc, SQRINC, (m, n, k, q.fortran_vec (), q.rows (),
-                             r.fortran_vec (), r.rows (), j + 1,
-                             utmp.data (), w));
-}
-
-void
-FloatQR::insert_col (const FloatMatrix& u, const Array<octave_idx_type>& j)
-{
-  octave_idx_type m = q.rows ();
-  octave_idx_type n = r.columns ();
-  octave_idx_type k = q.columns ();
-
-  Array<octave_idx_type> jsi;
-  Array<octave_idx_type> js = j.sort (jsi, 0, ASCENDING);
-  octave_idx_type nj = js.numel ();
-  bool dups = false;
-  for (octave_idx_type i = 0; i < nj - 1; i++)
-    dups = dups && js(i) == js(i+1);
-
-  if (dups)
-    (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
-  if (u.numel () != m || u.columns () != nj)
-    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
-  if (nj > 0 && (js(0) < 0 || js(nj-1) > n))
-    (*current_liboctave_error_handler) ("qrinsert: index out of range");
-
-  if (nj > 0)
-    {
-      octave_idx_type kmax = std::min (k + nj, m);
-      if (k < m)
-        {
-          q.resize (m, kmax);
-          r.resize (kmax, n + nj);
-        }
-      else
-        {
-          r.resize (k, n + nj);
-        }
-
-      OCTAVE_LOCAL_BUFFER (float, w, kmax);
-      for (volatile octave_idx_type i = 0; i < js.numel (); i++)
-        {
-          octave_idx_type ii = i;
-          FloatColumnVector utmp = u.column (jsi(i));
-          F77_XFCN (sqrinc, SQRINC, (m, n + ii, std::min (kmax, k + ii),
-                                     q.fortran_vec (), q.rows (),
-                                     r.fortran_vec (), r.rows (), js(ii) + 1,
-                                     utmp.data (), w));
-        }
-    }
-}
-
-void
-FloatQR::delete_col (octave_idx_type j)
-{
-  octave_idx_type m = q.rows ();
-  octave_idx_type k = r.rows ();
-  octave_idx_type n = r.columns ();
-
-  if (j < 0 || j > n-1)
-    (*current_liboctave_error_handler) ("qrdelete: index out of range");
-
-  OCTAVE_LOCAL_BUFFER (float, w, k);
-  F77_XFCN (sqrdec, SQRDEC, (m, n, k, q.fortran_vec (), q.rows (),
-                             r.fortran_vec (), r.rows (), j + 1, w));
-
-  if (k < m)
-    {
-      q.resize (m, k-1);
-      r.resize (k-1, n-1);
-    }
-  else
-    {
-      r.resize (k, n-1);
-    }
-}
-
-void
-FloatQR::delete_col (const Array<octave_idx_type>& j)
-{
-  octave_idx_type m = q.rows ();
-  octave_idx_type n = r.columns ();
-  octave_idx_type k = q.columns ();
-
-  Array<octave_idx_type> jsi;
-  Array<octave_idx_type> js = j.sort (jsi, 0, DESCENDING);
-  octave_idx_type nj = js.numel ();
-  bool dups = false;
-  for (octave_idx_type i = 0; i < nj - 1; i++)
-    dups = dups && js(i) == js(i+1);
-
-  if (dups)
-    (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
-  if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0))
-    (*current_liboctave_error_handler) ("qrinsert: index out of range");
-
-  if (nj > 0)
-    {
-      OCTAVE_LOCAL_BUFFER (float, w, k);
-      for (volatile octave_idx_type i = 0; i < js.numel (); i++)
-        {
-          octave_idx_type ii = i;
-          F77_XFCN (sqrdec, SQRDEC, (m, n - ii, k == m ? k : k - ii,
-                                     q.fortran_vec (), q.rows (),
-                                     r.fortran_vec (), r.rows (),
-                                     js(ii) + 1, w));
-        }
-      if (k < m)
-        {
-          q.resize (m, k - nj);
-          r.resize (k - nj, n - nj);
-        }
-      else
-        {
-          r.resize (k, n - nj);
-        }
-
-    }
-}
-
-void
-FloatQR::insert_row (const FloatRowVector& u, octave_idx_type j)
-{
-  octave_idx_type m = r.rows ();
-  octave_idx_type n = r.columns ();
-  octave_idx_type k = std::min (m, n);
-
-  if (! q.is_square () || u.numel () != n)
-    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
-  if (j < 0 || j > m)
-    (*current_liboctave_error_handler) ("qrinsert: index out of range");
-
-  q.resize (m + 1, m + 1);
-  r.resize (m + 1, n);
-  FloatRowVector utmp = u;
-  OCTAVE_LOCAL_BUFFER (float, w, k);
-  F77_XFCN (sqrinr, SQRINR, (m, n, q.fortran_vec (), q.rows (),
-                             r.fortran_vec (), r.rows (),
-                             j + 1, utmp.fortran_vec (), w));
-
-}
-
-void
-FloatQR::delete_row (octave_idx_type j)
-{
-  octave_idx_type m = r.rows ();
-  octave_idx_type n = r.columns ();
-
-  if (! q.is_square ())
-    (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch");
-  if (j < 0 || j > m-1)
-    (*current_liboctave_error_handler) ("qrdelete: index out of range");
-
-  OCTAVE_LOCAL_BUFFER (float, w, 2*m);
-  F77_XFCN (sqrder, SQRDER, (m, n, q.fortran_vec (), q.rows (),
-                             r.fortran_vec (), r.rows (), j + 1,
-                             w));
-
-  q.resize (m - 1, m - 1);
-  r.resize (m - 1, n);
-}
-
-void
-FloatQR::shift_cols (octave_idx_type i, octave_idx_type j)
-{
-  octave_idx_type m = q.rows ();
-  octave_idx_type k = r.rows ();
-  octave_idx_type n = r.columns ();
-
-  if (i < 0 || i > n-1 || j < 0 || j > n-1)
-    (*current_liboctave_error_handler) ("qrshift: index out of range");
-
-  OCTAVE_LOCAL_BUFFER (float, w, 2*k);
-  F77_XFCN (sqrshc, SQRSHC, (m, n, k,
-                             q.fortran_vec (), q.rows (),
-                             r.fortran_vec (), r.rows (),
-                             i + 1, j + 1, w));
-}
-
-#else
-
-// Replacement update methods.
-
-void
-FloatQR::update (const FloatColumnVector& u, const FloatColumnVector& v)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type m = q.rows ();
-  octave_idx_type n = r.columns ();
-
-  if (u.numel () != m || v.numel () != n)
-    (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
-
-  init (q*r + FloatMatrix (u) * FloatMatrix (v).transpose (), get_type ());
-}
-
-void
-FloatQR::update (const FloatMatrix& u, const FloatMatrix& v)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type m = q.rows ();
-  octave_idx_type n = r.columns ();
-
-  if (u.rows () != m || v.rows () != n || u.cols () != v.cols ())
-    (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
-
-  init (q*r + u * v.transpose (), get_type ());
-}
-
-static
-FloatMatrix insert_col (const FloatMatrix& a, octave_idx_type i,
-                        const FloatColumnVector& x)
-{
-  FloatMatrix retval (a.rows (), a.columns () + 1);
-  retval.assign (idx_vector::colon, idx_vector (0, i),
-                 a.index (idx_vector::colon, idx_vector (0, i)));
-  retval.assign (idx_vector::colon, idx_vector (i), x);
-  retval.assign (idx_vector::colon, idx_vector (i+1, retval.columns ()),
-                 a.index (idx_vector::colon, idx_vector (i, a.columns ())));
-  return retval;
-}
-
-static
-FloatMatrix insert_row (const FloatMatrix& a, octave_idx_type i,
-                        const FloatRowVector& x)
-{
-  FloatMatrix retval (a.rows () + 1, a.columns ());
-  retval.assign (idx_vector (0, i), idx_vector::colon,
-                 a.index (idx_vector (0, i), idx_vector::colon));
-  retval.assign (idx_vector (i), idx_vector::colon, x);
-  retval.assign (idx_vector (i+1, retval.rows ()), idx_vector::colon,
-                 a.index (idx_vector (i, a.rows ()), idx_vector::colon));
-  return retval;
-}
-
-static
-FloatMatrix delete_col (const FloatMatrix& a, octave_idx_type i)
-{
-  FloatMatrix retval = a;
-  retval.delete_elements (1, idx_vector (i));
-  return retval;
-}
-
-static
-FloatMatrix delete_row (const FloatMatrix& a, octave_idx_type i)
-{
-  FloatMatrix retval = a;
-  retval.delete_elements (0, idx_vector (i));
-  return retval;
-}
-
-static
-FloatMatrix shift_cols (const FloatMatrix& a,
-                        octave_idx_type i, octave_idx_type j)
-{
-  octave_idx_type n = a.columns ();
-  Array<octave_idx_type> p (dim_vector (n, 1));
-  for (octave_idx_type k = 0; k < n; k++) p(k) = k;
-  if (i < j)
-    {
-      for (octave_idx_type k = i; k < j; k++) p(k) = k+1;
-      p(j) = i;
-    }
-  else if (j < i)
-    {
-      p(j) = i;
-      for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1;
-    }
-
-  return a.index (idx_vector::colon, idx_vector (p));
-}
-
-void
-FloatQR::insert_col (const FloatColumnVector& u, octave_idx_type j)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type m = q.rows ();
-  octave_idx_type n = r.columns ();
-
-  if (u.numel () != m)
-    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
-  if (j < 0 || j > n)
-    (*current_liboctave_error_handler) ("qrinsert: index out of range");
-
-  init (::insert_col (q*r, j, u), get_type ());
-}
-
-void
-FloatQR::insert_col (const FloatMatrix& u, const Array<octave_idx_type>& j)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type m = q.rows ();
-  octave_idx_type n = r.columns ();
-
-  Array<octave_idx_type> jsi;
-  Array<octave_idx_type> js = j.sort (jsi, 0, ASCENDING);
-  octave_idx_type nj = js.numel ();
-  bool dups = false;
-  for (octave_idx_type i = 0; i < nj - 1; i++)
-    dups = dups && js(i) == js(i+1);
-
-  if (dups)
-    (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
-  if (u.numel () != m || u.columns () != nj)
-    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
-  if (nj > 0 && (js(0) < 0 || js(nj-1) > n))
-    (*current_liboctave_error_handler) ("qrinsert: index out of range");
-
-  if (nj > 0)
-    {
-      FloatMatrix a = q*r;
-      for (octave_idx_type i = 0; i < js.numel (); i++)
-        a = ::insert_col (a, js(i), u.column (i));
-      init (a, get_type ());
-    }
-}
-
-void
-FloatQR::delete_col (octave_idx_type j)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type n = r.columns ();
-
-  if (j < 0 || j > n-1)
-    (*current_liboctave_error_handler) ("qrdelete: index out of range");
-
-  init (::delete_col (q*r, j), get_type ());
-}
-
-void
-FloatQR::delete_col (const Array<octave_idx_type>& j)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type n = r.columns ();
-
-  Array<octave_idx_type> jsi;
-  Array<octave_idx_type> js = j.sort (jsi, 0, DESCENDING);
-  octave_idx_type nj = js.numel ();
-  bool dups = false;
-  for (octave_idx_type i = 0; i < nj - 1; i++)
-    dups = dups && js(i) == js(i+1);
-
-  if (dups)
-    (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
-  if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0))
-    (*current_liboctave_error_handler) ("qrinsert: index out of range");
-
-  if (nj > 0)
-    {
-      FloatMatrix a = q*r;
-      for (octave_idx_type i = 0; i < js.numel (); i++)
-        a = ::delete_col (a, js(i));
-      init (a, get_type ());
-    }
-}
-
-void
-FloatQR::insert_row (const FloatRowVector& u, octave_idx_type j)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type m = r.rows ();
-  octave_idx_type n = r.columns ();
-
-  if (! q.is_square () || u.numel () != n)
-    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
-  if (j < 0 || j > m)
-    (*current_liboctave_error_handler) ("qrinsert: index out of range");
-
-  init (::insert_row (q*r, j, u), get_type ());
-}
-
-void
-FloatQR::delete_row (octave_idx_type j)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type m = r.rows ();
-
-  if (! q.is_square ())
-    (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch");
-  if (j < 0 || j > m-1)
-    (*current_liboctave_error_handler) ("qrdelete: index out of range");
-
-  init (::delete_row (q*r, j), get_type ());
-}
-
-void
-FloatQR::shift_cols (octave_idx_type i, octave_idx_type j)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type n = r.columns ();
-
-  if (i < 0 || i > n-1 || j < 0 || j > n-1)
-    (*current_liboctave_error_handler) ("qrshift: index out of range");
-
-  init (::shift_cols (q*r, i, j), get_type ());
-}
-
-#endif
diff -r 48d82f74243e -r eb1524b07fe3 liboctave/numeric/floatQR.h
--- a/liboctave/numeric/floatQR.h	Tue Feb 16 23:34:44 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,77 +0,0 @@
-/*
-
-Copyright (C) 1994-2015 John W. Eaton
-Copyright (C) 2008-2009 Jaroslav Hajek
-
-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/>.
-
-*/
-
-#if ! defined (octave_floatQR_h)
-#define octave_floatQR_h 1
-
-#include "octave-config.h"
-
-#include <iosfwd>
-
-#include "fMatrix.h"
-#include "fColVector.h"
-#include "fRowVector.h"
-#include "base-qr.h"
-
-class
-OCTAVE_API
-FloatQR : public base_qr<FloatMatrix>
-{
-public:
-
-  FloatQR (void) : base_qr<FloatMatrix> () { }
-
-  FloatQR (const FloatMatrix&, qr_type_t = qr_type_std);
-
-  FloatQR (const FloatMatrix& qx, const FloatMatrix& rx)
-    : base_qr<FloatMatrix> (qx, rx) { }
-
-  FloatQR (const FloatQR& a) : base_qr<FloatMatrix> (a) { }
-
-  void init (const FloatMatrix&, qr_type_t);
-
-  void update (const FloatColumnVector& u, const FloatColumnVector& v);
-
-  void update (const FloatMatrix& u, const FloatMatrix& v);
-
-  void insert_col (const FloatColumnVector& u, octave_idx_type j);
-
-  void insert_col (const FloatMatrix& u, const Array<octave_idx_type>& j);
-
-  void delete_col (octave_idx_type j);
-
-  void delete_col (const Array<octave_idx_type>& j);
-
-  void insert_row (const FloatRowVector& u, octave_idx_type j);
-
-  void delete_row (octave_idx_type j);
-
-  void shift_cols (octave_idx_type i, octave_idx_type j);
-
-protected:
-
-  void form (octave_idx_type n, FloatMatrix& afact,
-             float *tau, qr_type_t qr_type);
-};
-
-#endif
diff -r 48d82f74243e -r eb1524b07fe3 liboctave/numeric/floatQRP.cc
--- a/liboctave/numeric/floatQRP.cc	Tue Feb 16 23:34:44 2016 -0500
+++ b/liboctave/numeric/floatQRP.cc	Wed Feb 17 02:57:21 2016 -0500
@@ -43,16 +43,16 @@
 
 // It would be best to share some of this code with QR class...
 
-FloatQRP::FloatQRP (const FloatMatrix& a, qr_type_t qr_type)
-  : FloatQR (), p ()
+FloatQRP::FloatQRP (const FloatMatrix& a, type qr_type)
+  : qr<FloatMatrix> (), p ()
 {
   init (a, qr_type);
 }
 
 void
-FloatQRP::init (const FloatMatrix& a, qr_type_t qr_type)
+FloatQRP::init (const FloatMatrix& a, type qr_type)
 {
-  assert (qr_type != qr_type_raw);
+  assert (qr_type != qr<FloatMatrix>::raw);
 
   octave_idx_type m = a.rows ();
   octave_idx_type n = a.cols ();
@@ -63,7 +63,7 @@
   octave_idx_type info = 0;
 
   FloatMatrix afact = a;
-  if (m > n && qr_type == qr_type_std)
+  if (m > n && qr_type == qr<FloatMatrix>::std)
     afact.resize (m, m);
 
   MArray<octave_idx_type> jpvt (dim_vector (n, 1), 0);
diff -r 48d82f74243e -r eb1524b07fe3 liboctave/numeric/floatQRP.h
--- a/liboctave/numeric/floatQRP.h	Tue Feb 16 23:34:44 2016 -0500
+++ b/liboctave/numeric/floatQRP.h	Wed Feb 17 02:57:21 2016 -0500
@@ -27,27 +27,31 @@
 
 #include <iosfwd>
 
-#include "floatQR.h"
 #include "PermMatrix.h"
 #include "fColVector.h"
+#include "fMatrix.h"
+#include "fRowVector.h"
+#include "qr.h"
 
 class
 OCTAVE_API
-FloatQRP : public FloatQR
+FloatQRP : public qr<FloatMatrix>
 {
 public:
 
-  FloatQRP (void) : FloatQR (), p () { }
+  typedef qr<FloatMatrix>::type type;
+
+  FloatQRP (void) : qr<FloatMatrix> (), p () { }
 
-  FloatQRP (const FloatMatrix&, qr_type_t = qr_type_std);
+  FloatQRP (const FloatMatrix&, type = qr<FloatMatrix>::std);
 
-  FloatQRP (const FloatQRP& a) : FloatQR (a), p (a.p) { }
+  FloatQRP (const FloatQRP& a) : qr<FloatMatrix> (a), p (a.p) { }
 
   FloatQRP& operator = (const FloatQRP& a)
   {
     if (this != &a)
       {
-        FloatQR::operator = (a);
+        qr<FloatMatrix>::operator = (a);
         p = a.p;
       }
 
@@ -56,7 +60,7 @@
 
   ~FloatQRP (void) { }
 
-  void init (const FloatMatrix&, qr_type_t = qr_type_std);
+  void init (const FloatMatrix&, type = qr<FloatMatrix>::std);
 
   PermMatrix P (void) const { return p; }
 
diff -r 48d82f74243e -r eb1524b07fe3 liboctave/numeric/module.mk
--- a/liboctave/numeric/module.mk	Tue Feb 16 23:34:44 2016 -0500
+++ b/liboctave/numeric/module.mk	Wed Feb 17 02:57:21 2016 -0500
@@ -8,7 +8,6 @@
 LIBOCTAVE_OPT_IN = $(LIBOCTAVE_OPT_INC:.h=.in)
 
 NUMERIC_INC = \
-  liboctave/numeric/CmplxQR.h \
   liboctave/numeric/CmplxQRP.h \
   liboctave/numeric/CollocWt.h \
   liboctave/numeric/DAE.h \
@@ -30,17 +29,13 @@
   liboctave/numeric/base-dae.h \
   liboctave/numeric/base-de.h \
   liboctave/numeric/base-min.h \
-  liboctave/numeric/base-qr.h \
   liboctave/numeric/bsxfun-decl.h \
   liboctave/numeric/bsxfun.h \
   liboctave/numeric/chol.h \
-  liboctave/numeric/dbleQR.h \
   liboctave/numeric/dbleQRP.h \
   liboctave/numeric/eigs-base.h \
-  liboctave/numeric/fCmplxQR.h \
   liboctave/numeric/fCmplxQRP.h \
   liboctave/numeric/fEIG.h \
-  liboctave/numeric/floatQR.h \
   liboctave/numeric/floatQRP.h \
   liboctave/numeric/gepbalance.h \
   liboctave/numeric/hess.h \
@@ -52,6 +47,7 @@
   liboctave/numeric/oct-norm.h \
   liboctave/numeric/oct-rand.h \
   liboctave/numeric/oct-spparms.h \
+  liboctave/numeric/qr.h \
   liboctave/numeric/randgamma.h \
   liboctave/numeric/randmtzig.h \
   liboctave/numeric/randpoisson.h \
@@ -68,7 +64,6 @@
   liboctave/numeric/randpoisson.c
 
 NUMERIC_SRC = \
-  liboctave/numeric/CmplxQR.cc \
   liboctave/numeric/CmplxQRP.cc \
   liboctave/numeric/CollocWt.cc \
   liboctave/numeric/DASPK.cc \
@@ -80,13 +75,10 @@
   liboctave/numeric/Quad.cc \
   liboctave/numeric/aepbalance.cc \
   liboctave/numeric/chol.cc \
-  liboctave/numeric/dbleQR.cc \
   liboctave/numeric/dbleQRP.cc \
   liboctave/numeric/eigs-base.cc \
-  liboctave/numeric/fCmplxQR.cc \
   liboctave/numeric/fCmplxQRP.cc \
   liboctave/numeric/fEIG.cc \
-  liboctave/numeric/floatQR.cc \
   liboctave/numeric/floatQRP.cc \
   liboctave/numeric/gepbalance.cc \
   liboctave/numeric/hess.cc \
@@ -98,6 +90,7 @@
   liboctave/numeric/oct-norm.cc \
   liboctave/numeric/oct-rand.cc \
   liboctave/numeric/oct-spparms.cc \
+  liboctave/numeric/qr.cc \
   liboctave/numeric/schur.cc \
   liboctave/numeric/sparse-chol.cc \
   liboctave/numeric/sparse-dmsolve.cc \
@@ -107,7 +100,6 @@
   $(NUMERIC_C_SRC)
 
 LIBOCTAVE_TEMPLATE_SRC += \
-  liboctave/numeric/base-qr.cc \
   liboctave/numeric/bsxfun-defs.cc
 
 ## Special rules for sources which must be built before rest of compilation.
diff -r 48d82f74243e -r eb1524b07fe3 liboctave/numeric/qr.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/numeric/qr.cc	Wed Feb 17 02:57:21 2016 -0500
@@ -0,0 +1,2049 @@
+/*
+
+Copyright (C) 1994-2015 John W. Eaton
+Copyright (C) 2008-2009 Jaroslav Hajek
+Copyright (C) 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 "CColVector.h"
+#include "CMatrix.h"
+#include "CRowVector.h"
+#include "dColVector.h"
+#include "dMatrix.h"
+#include "dRowVector.h"
+#include "f77-fcn.h"
+#include "fCColVector.h"
+#include "fCMatrix.h"
+#include "fCRowVector.h"
+#include "fColVector.h"
+#include "fMatrix.h"
+#include "fRowVector.h"
+#include "idx-vector.h"
+#include "lo-error.h"
+#include "oct-locbuf.h"
+#include "qr.h"
+
+extern "C"
+{
+  F77_RET_T
+  F77_FUNC (dgeqrf, DGEQRF) (const octave_idx_type&, const octave_idx_type&,
+                             double*, const octave_idx_type&, double*,
+                             double*, const octave_idx_type&,
+                             octave_idx_type&);
+
+  F77_RET_T
+  F77_FUNC (dorgqr, DORGQR) (const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, double*,
+                             const octave_idx_type&, double*, double*,
+                             const octave_idx_type&, octave_idx_type&);
+
+#if defined (HAVE_QRUPDATE)
+
+  F77_RET_T
+  F77_FUNC (dqr1up, DQR1UP) (const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, double*,
+                             const octave_idx_type&, double*,
+                             const octave_idx_type&, double*, double*, double*);
+
+  F77_RET_T
+  F77_FUNC (dqrinc, DQRINC) (const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, double*,
+                             const octave_idx_type&, double*,
+                             const octave_idx_type&, const octave_idx_type&,
+                             const double*, double*);
+
+  F77_RET_T
+  F77_FUNC (dqrdec, DQRDEC) (const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, double*,
+                             const octave_idx_type&, double*,
+                             const octave_idx_type&, const octave_idx_type&,
+                             double*);
+
+  F77_RET_T
+  F77_FUNC (dqrinr, DQRINR) (const octave_idx_type&, const octave_idx_type&,
+                             double*, const octave_idx_type&, double*,
+                             const octave_idx_type&, const octave_idx_type&,
+                             const double*, double*);
+
+  F77_RET_T
+  F77_FUNC (dqrder, DQRDER) (const octave_idx_type&, const octave_idx_type&,
+                             double*, const octave_idx_type&, double*,
+                             const octave_idx_type&, const octave_idx_type&,
+                             double*);
+
+  F77_RET_T
+  F77_FUNC (dqrshc, DQRSHC) (const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, double*,
+                             const octave_idx_type&, double*,
+                             const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, double*);
+
+#endif
+
+  F77_RET_T
+  F77_FUNC (sgeqrf, SGEQRF) (const octave_idx_type&, const octave_idx_type&,
+                             float*, const octave_idx_type&, float*, float*,
+                             const octave_idx_type&, octave_idx_type&);
+
+  F77_RET_T
+  F77_FUNC (sorgqr, SORGQR) (const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, float*,
+                             const octave_idx_type&, float*, float*,
+                             const octave_idx_type&, octave_idx_type&);
+
+#if defined (HAVE_QRUPDATE)
+
+  F77_RET_T
+  F77_FUNC (sqr1up, SQR1UP) (const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, float*,
+                             const octave_idx_type&, float*,
+                             const octave_idx_type&, float*, float*, float*);
+
+  F77_RET_T
+  F77_FUNC (sqrinc, SQRINC) (const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, float*,
+                             const octave_idx_type&, float*,
+                             const octave_idx_type&,
+                             const octave_idx_type&, const float*, float*);
+
+  F77_RET_T
+  F77_FUNC (sqrdec, SQRDEC) (const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, float*,
+                             const octave_idx_type&, float*,
+                             const octave_idx_type&,
+                             const octave_idx_type&, float*);
+
+  F77_RET_T
+  F77_FUNC (sqrinr, SQRINR) (const octave_idx_type&, const octave_idx_type&,
+                             float*, const octave_idx_type&,
+                             float*, const octave_idx_type&,
+                             const octave_idx_type&, const float*, float*);
+
+  F77_RET_T
+  F77_FUNC (sqrder, SQRDER) (const octave_idx_type&, const octave_idx_type&,
+                             float*, const octave_idx_type&,
+                             float*, const octave_idx_type&,
+                             const octave_idx_type&, float*);
+
+  F77_RET_T
+  F77_FUNC (sqrshc, SQRSHC) (const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, float*,
+                             const octave_idx_type&, float*,
+                             const octave_idx_type&,
+                             const octave_idx_type&, const octave_idx_type&,
+                             float*);
+
+#endif
+
+  F77_RET_T
+  F77_FUNC (zgeqrf, ZGEQRF) (const octave_idx_type&, const octave_idx_type&,
+                             Complex*, const octave_idx_type&, Complex*,
+                             Complex*, const octave_idx_type&,
+                             octave_idx_type&);
+
+  F77_RET_T
+  F77_FUNC (zungqr, ZUNGQR) (const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, Complex*,
+                             const octave_idx_type&, Complex*, Complex*,
+                             const octave_idx_type&, octave_idx_type&);
+
+#if defined (HAVE_QRUPDATE)
+
+  F77_RET_T
+  F77_FUNC (zqr1up, ZQR1UP) (const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, Complex*,
+                             const octave_idx_type&, Complex*,
+                             const octave_idx_type&, Complex*,
+                             Complex*, Complex*, double*);
+
+  F77_RET_T
+  F77_FUNC (zqrinc, ZQRINC) (const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, Complex*,
+                             const octave_idx_type&, Complex*,
+                             const octave_idx_type&, const octave_idx_type&,
+                             const Complex*, double*);
+
+  F77_RET_T
+  F77_FUNC (zqrdec, ZQRDEC) (const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, Complex*,
+                             const octave_idx_type&, Complex*,
+                             const octave_idx_type&, const octave_idx_type&,
+                             double*);
+
+  F77_RET_T
+  F77_FUNC (zqrinr, ZQRINR) (const octave_idx_type&, const octave_idx_type&,
+                             Complex*, const octave_idx_type&, Complex*,
+                             const octave_idx_type&, const octave_idx_type&,
+                             const Complex*, double*);
+
+  F77_RET_T
+  F77_FUNC (zqrder, ZQRDER) (const octave_idx_type&, const octave_idx_type&,
+                             Complex*, const octave_idx_type&, Complex*,
+                             const octave_idx_type&, const octave_idx_type&,
+                             Complex*, double*);
+
+  F77_RET_T
+  F77_FUNC (zqrshc, ZQRSHC) (const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, Complex*,
+                             const octave_idx_type&, Complex*,
+                             const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, Complex*, double*);
+
+#endif
+
+  F77_RET_T
+  F77_FUNC (cgeqrf, CGEQRF) (const octave_idx_type&, const octave_idx_type&,
+                             FloatComplex*, const octave_idx_type&,
+                             FloatComplex*, FloatComplex*,
+                             const octave_idx_type&, octave_idx_type&);
+
+  F77_RET_T
+  F77_FUNC (cungqr, CUNGQR) (const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, FloatComplex*,
+                             const octave_idx_type&, FloatComplex*,
+                             FloatComplex*, const octave_idx_type&,
+                             octave_idx_type&);
+
+#ifdef HAVE_QRUPDATE
+
+  F77_RET_T
+  F77_FUNC (cqr1up, CQR1UP) (const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, FloatComplex*,
+                             const octave_idx_type&, FloatComplex*,
+                             const octave_idx_type&, FloatComplex*,
+                             FloatComplex*, FloatComplex*, float*);
+
+  F77_RET_T
+  F77_FUNC (cqrinc, CQRINC) (const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, FloatComplex*,
+                             const octave_idx_type&, FloatComplex*,
+                             const octave_idx_type&,const octave_idx_type&,
+                             const FloatComplex*, float*);
+
+  F77_RET_T
+  F77_FUNC (cqrdec, CQRDEC) (const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, FloatComplex*,
+                             const octave_idx_type&, FloatComplex*,
+                             const octave_idx_type&, const octave_idx_type&,
+                             float*);
+
+  F77_RET_T
+  F77_FUNC (cqrinr, CQRINR) (const octave_idx_type&, const octave_idx_type&,
+                             FloatComplex*, const octave_idx_type&,
+                             FloatComplex*, const octave_idx_type&,
+                             const octave_idx_type&, const FloatComplex*,
+                             float*);
+
+  F77_RET_T
+  F77_FUNC (cqrder, CQRDER) (const octave_idx_type&, const octave_idx_type&,
+                             FloatComplex*, const octave_idx_type&,
+                             FloatComplex*, const octave_idx_type&,
+                             const octave_idx_type&, FloatComplex*, float*);
+
+  F77_RET_T
+  F77_FUNC (cqrshc, CQRSHC) (const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, FloatComplex*,
+                             const octave_idx_type&, FloatComplex*,
+                             const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, FloatComplex*,
+                             float*);
+
+#endif
+}
+
+template <typename T>
+qr<T>::qr (const T& q_arg, const T& r_arg)
+  : q (q_arg), r (r_arg)
+{
+  octave_idx_type q_nr = q.rows ();
+  octave_idx_type q_nc = q.columns ();
+
+  octave_idx_type r_nr = r.rows ();
+  octave_idx_type r_nc = r.columns ();
+
+  if (! (q_nc == r_nr && (q_nr == q_nc || (q_nr > q_nc && r_nr == r_nc))))
+    (*current_liboctave_error_handler) ("QR dimensions mismatch");
+}
+
+template <typename T>
+typename qr<T>::type
+qr<T>::get_type (void) const
+{
+  type retval;
+
+  if (! q.is_empty () && q.is_square ())
+    retval = qr<T>::std;
+  else if (q.rows () > q.columns () && r.is_square ())
+    retval = qr<T>::economy;
+  else
+    retval = qr<T>::raw;
+
+  return retval;
+}
+
+template <typename T>
+bool
+qr<T>::regular (void) const
+{
+  bool retval = true;
+
+  octave_idx_type k = std::min (r.rows (), r.columns ());
+
+  for (octave_idx_type i = 0; i < k; i++)
+    {
+      if (r(i, i) == ELT_T ())
+        {
+          retval = false;
+          break;
+        }
+    }
+
+  return retval;
+}
+
+#if ! defined (HAVE_QRUPDATE)
+
+// Replacement update methods.
+
+static
+void
+warn_qrupdate_once (void)
+{
+  static bool warned = false;
+
+  if (! warned)
+    {
+      (*current_liboctave_warning_with_id_handler)
+        ("Octave:missing-dependency",
+         "In this version of Octave, QR & Cholesky updating routines "
+         "simply update the matrix and recalculate factorizations. "
+         "To use fast algorithms, link Octave with the qrupdate library. "
+         "See <http://sourceforge.net/projects/qrupdate>.");
+
+      warned = true;
+    }
+}
+
+template <typename T>
+void
+qr<T>::update (const CV_T& u, const CV_T& v)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (u.numel () != m || v.numel () != n)
+    (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
+
+  init (q*r + T (u) * T (v).hermitian (), get_type ());
+}
+
+template <typename T>
+void
+qr<T>::update (const T& u, const T& v)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (u.rows () != m || v.rows () != n || u.cols () != v.cols ())
+    (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
+
+  init (q*r + u * v.hermitian (), get_type ());
+}
+
+template <typename T>
+static
+T
+insert_col (const T& a, octave_idx_type i, const CV_T& x)
+{
+  T retval (a.rows (), a.columns () + 1);
+  retval.assign (idx_vector::colon, idx_vector (0, i),
+                 a.index (idx_vector::colon, idx_vector (0, i)));
+  retval.assign (idx_vector::colon, idx_vector (i), x);
+  retval.assign (idx_vector::colon, idx_vector (i+1, retval.columns ()),
+                 a.index (idx_vector::colon, idx_vector (i, a.columns ())));
+  return retval;
+}
+
+template <typename T>
+static
+T
+insert_row (const T& a, octave_idx_type i, const RV_T& x)
+{
+  T retval (a.rows () + 1, a.columns ());
+  retval.assign (idx_vector (0, i), idx_vector::colon,
+                 a.index (idx_vector (0, i), idx_vector::colon));
+  retval.assign (idx_vector (i), idx_vector::colon, x);
+  retval.assign (idx_vector (i+1, retval.rows ()), idx_vector::colon,
+                 a.index (idx_vector (i, a.rows ()), idx_vector::colon));
+  return retval;
+}
+
+template <typename T>
+static
+T
+delete_col (const T& a, octave_idx_type i)
+{
+  T retval = a;
+  retval.delete_elements (1, idx_vector (i));
+  return retval;
+}
+
+template <typename T>
+static
+T
+delete_row (const T& a, octave_idx_type i)
+{
+  T retval = a;
+  retval.delete_elements (0, idx_vector (i));
+  return retval;
+}
+
+template <typename T>
+static
+T
+shift_cols (const T& a, octave_idx_type i, octave_idx_type j)
+{
+  octave_idx_type n = a.columns ();
+  Array<octave_idx_type> p (dim_vector (n, 1));
+  for (octave_idx_type k = 0; k < n; k++) p(k) = k;
+  if (i < j)
+    {
+      for (octave_idx_type k = i; k < j; k++) p(k) = k+1;
+      p(j) = i;
+    }
+  else if (j < i)
+    {
+      p(j) = i;
+      for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1;
+    }
+
+  return a.index (idx_vector::colon, idx_vector (p));
+}
+
+template <typename T>
+void
+qr<T>::insert_col (const CV_T& u, octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (u.numel () != m)
+    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
+  if (j < 0 || j > n)
+    (*current_liboctave_error_handler) ("qrinsert: index out of range");
+
+  init (::insert_col (q*r, j, u), get_type ());
+}
+
+template <typename T>
+void
+qr<T>::insert_col (const T& u, const Array<octave_idx_type>& j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+
+  Array<octave_idx_type> jsi;
+  Array<octave_idx_type> js = j.sort (jsi, 0, ASCENDING);
+  octave_idx_type nj = js.numel ();
+  bool dups = false;
+  for (octave_idx_type i = 0; i < nj - 1; i++)
+    dups = dups && js(i) == js(i+1);
+
+  if (dups)
+    (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
+  if (u.numel () != m || u.columns () != nj)
+    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
+  if (nj > 0 && (js(0) < 0 || js(nj-1) > n))
+    (*current_liboctave_error_handler) ("qrinsert: index out of range");
+
+  if (nj > 0)
+    {
+      T a = q*r;
+      for (octave_idx_type i = 0; i < js.numel (); i++)
+        a = ::insert_col (a, js(i), u.column (i));
+      init (a, get_type ());
+    }
+}
+
+template <typename T>
+void
+qr<T>::delete_col (octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type n = r.columns ();
+
+  if (j < 0 || j > n-1)
+    (*current_liboctave_error_handler) ("qrdelete: index out of range");
+
+  init (::delete_col (q*r, j), get_type ());
+}
+
+template <typename T>
+void
+qr<T>::delete_col (const Array<octave_idx_type>& j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type n = r.columns ();
+
+  Array<octave_idx_type> jsi;
+  Array<octave_idx_type> js = j.sort (jsi, 0, DESCENDING);
+  octave_idx_type nj = js.numel ();
+  bool dups = false;
+  for (octave_idx_type i = 0; i < nj - 1; i++)
+    dups = dups && js(i) == js(i+1);
+
+  if (dups)
+    (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
+  if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0))
+    (*current_liboctave_error_handler) ("qrinsert: index out of range");
+
+  if (nj > 0)
+    {
+      T a = q*r;
+      for (octave_idx_type i = 0; i < js.numel (); i++)
+        a = ::delete_col (a, js(i));
+      init (a, get_type ());
+    }
+}
+
+template <typename T>
+void
+qr<T>::insert_row (const RV_T& u, octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = r.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (! q.is_square () || u.numel () != n)
+    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
+  if (j < 0 || j > m)
+    (*current_liboctave_error_handler) ("qrinsert: index out of range");
+
+  init (::insert_row (q*r, j, u), get_type ());
+}
+
+template <typename T>
+void
+qr<T>::delete_row (octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type m = r.rows ();
+
+  if (! q.is_square ())
+    (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch");
+  if (j < 0 || j > m-1)
+    (*current_liboctave_error_handler) ("qrdelete: index out of range");
+
+  init (::delete_row (q*r, j), get_type ());
+}
+
+template <typename T>
+void
+qr<T>::shift_cols (octave_idx_type i, octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type n = r.columns ();
+
+  if (i < 0 || i > n-1 || j < 0 || j > n-1)
+    (*current_liboctave_error_handler) ("qrshift: index out of range");
+
+  init (::shift_cols (q*r, i, j), get_type ());
+}
+
+#endif
+
+// Specializations.
+
+template <>
+void
+qr<Matrix>::form (octave_idx_type n, Matrix& afact, double *tau, type qr_type)
+{
+  octave_idx_type m = afact.rows ();
+  octave_idx_type min_mn = std::min (m, n);
+  octave_idx_type info;
+
+  if (qr_type == qr<Matrix>::raw)
+    {
+      for (octave_idx_type j = 0; j < min_mn; j++)
+        {
+          octave_idx_type limit = j < min_mn - 1 ? j : min_mn - 1;
+          for (octave_idx_type i = limit + 1; i < m; i++)
+            afact.elem (i, j) *= tau[j];
+        }
+
+      r = afact;
+    }
+  else
+    {
+      // Attempt to minimize copying.
+      if (m >= n)
+        {
+          // afact will become q.
+          q = afact;
+          octave_idx_type k = qr_type == qr<Matrix>::economy ? n : m;
+          r = Matrix (k, n);
+          for (octave_idx_type j = 0; j < n; j++)
+            {
+              octave_idx_type i = 0;
+              for (; i <= j; i++)
+                r.xelem (i, j) = afact.xelem (i, j);
+              for (; i < k; i++)
+                r.xelem (i, j) = 0;
+            }
+          afact = Matrix (); // optimize memory
+        }
+      else
+        {
+          // afact will become r.
+          q = Matrix (m, m);
+          for (octave_idx_type j = 0; j < m; j++)
+            for (octave_idx_type i = j + 1; i < m; i++)
+              {
+                q.xelem (i, j) = afact.xelem (i, j);
+                afact.xelem (i, j) = 0;
+              }
+          r = afact;
+        }
+
+
+      if (m > 0)
+        {
+          octave_idx_type k = q.columns ();
+          // workspace query.
+          double rlwork;
+          F77_XFCN (dorgqr, DORGQR, (m, k, min_mn, q.fortran_vec (), m, tau,
+                                     &rlwork, -1, info));
+
+          // allocate buffer and do the job.
+          octave_idx_type lwork = rlwork;
+          lwork = std::max (lwork, static_cast<octave_idx_type> (1));
+          OCTAVE_LOCAL_BUFFER (double, work, lwork);
+          F77_XFCN (dorgqr, DORGQR, (m, k, min_mn, q.fortran_vec (), m, tau,
+                                     work, lwork, info));
+        }
+    }
+}
+
+template <>
+void
+qr<Matrix>::init (const Matrix& a, type qr_type)
+{
+  octave_idx_type m = a.rows ();
+  octave_idx_type n = a.cols ();
+
+  octave_idx_type min_mn = m < n ? m : n;
+  OCTAVE_LOCAL_BUFFER (double, tau, min_mn);
+
+  octave_idx_type info = 0;
+
+  Matrix afact = a;
+  if (m > n && qr_type == qr<Matrix>::std)
+    afact.resize (m, m);
+
+  if (m > 0)
+    {
+      // workspace query.
+      double rlwork;
+      F77_XFCN (dgeqrf, DGEQRF, (m, n, afact.fortran_vec (), m, tau,
+                                 &rlwork, -1, info));
+
+      // allocate buffer and do the job.
+      octave_idx_type lwork = rlwork;
+      lwork = std::max (lwork, static_cast<octave_idx_type> (1));
+      OCTAVE_LOCAL_BUFFER (double, work, lwork);
+      F77_XFCN (dgeqrf, DGEQRF, (m, n, afact.fortran_vec (), m, tau,
+                                 work, lwork, info));
+    }
+
+  form (n, afact, tau, qr_type);
+}
+
+#if defined (HAVE_QRUPDATE)
+
+template <>
+void
+qr<Matrix>::update (const ColumnVector& u, const ColumnVector& v)
+{
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+  octave_idx_type k = q.columns ();
+
+  if (u.numel () != m || v.numel () != n)
+    (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
+
+  ColumnVector utmp = u;
+  ColumnVector vtmp = v;
+  OCTAVE_LOCAL_BUFFER (double, w, 2*k);
+  F77_XFCN (dqr1up, DQR1UP, (m, n, k, q.fortran_vec (),
+                             m, r.fortran_vec (), k,
+                             utmp.fortran_vec (), vtmp.fortran_vec (), w));
+}
+
+template <>
+void
+qr<Matrix>::update (const Matrix& u, const Matrix& v)
+{
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+  octave_idx_type k = q.columns ();
+
+  if (u.rows () != m || v.rows () != n || u.cols () != v.cols ())
+    (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
+
+  OCTAVE_LOCAL_BUFFER (double, w, 2*k);
+  for (volatile octave_idx_type i = 0; i < u.cols (); i++)
+    {
+      ColumnVector utmp = u.column (i);
+      ColumnVector vtmp = v.column (i);
+      F77_XFCN (dqr1up, DQR1UP, (m, n, k, q.fortran_vec (),
+                                 m, r.fortran_vec (), k,
+                                 utmp.fortran_vec (), vtmp.fortran_vec (),
+                                 w));
+    }
+}
+
+template <>
+void
+qr<Matrix>::insert_col (const ColumnVector& u, octave_idx_type j)
+{
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+  octave_idx_type k = q.columns ();
+
+  if (u.numel () != m)
+    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
+  if (j < 0 || j > n)
+    (*current_liboctave_error_handler) ("qrinsert: index out of range");
+
+  if (k < m)
+    {
+      q.resize (m, k+1);
+      r.resize (k+1, n+1);
+    }
+  else
+    {
+      r.resize (k, n+1);
+    }
+
+  ColumnVector utmp = u;
+  OCTAVE_LOCAL_BUFFER (double, w, k);
+  F77_XFCN (dqrinc, DQRINC, (m, n, k, q.fortran_vec (), q.rows (),
+                             r.fortran_vec (), r.rows (), j + 1,
+                             utmp.data (), w));
+}
+
+template <>
+void
+qr<Matrix>::insert_col (const Matrix& u, const Array<octave_idx_type>& j)
+{
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+  octave_idx_type k = q.columns ();
+
+  Array<octave_idx_type> jsi;
+  Array<octave_idx_type> js = j.sort (jsi, 0, ASCENDING);
+  octave_idx_type nj = js.numel ();
+  bool dups = false;
+  for (octave_idx_type i = 0; i < nj - 1; i++)
+    dups = dups && js(i) == js(i+1);
+
+  if (dups)
+    (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
+  if (u.numel () != m || u.columns () != nj)
+    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
+  if (nj > 0 && (js(0) < 0 || js(nj-1) > n))
+    (*current_liboctave_error_handler) ("qrinsert: index out of range");
+
+  if (nj > 0)
+    {
+      octave_idx_type kmax = std::min (k + nj, m);
+      if (k < m)
+        {
+          q.resize (m, kmax);
+          r.resize (kmax, n + nj);
+        }
+      else
+        {
+          r.resize (k, n + nj);
+        }
+
+      OCTAVE_LOCAL_BUFFER (double, w, kmax);
+      for (volatile octave_idx_type i = 0; i < js.numel (); i++)
+        {
+          octave_idx_type ii = i;
+          ColumnVector utmp = u.column (jsi(i));
+          F77_XFCN (dqrinc, DQRINC, (m, n + ii, std::min (kmax, k + ii),
+                                     q.fortran_vec (), q.rows (),
+                                     r.fortran_vec (), r.rows (), js(ii) + 1,
+                                     utmp.data (), w));
+        }
+    }
+}
+
+template <>
+void
+qr<Matrix>::delete_col (octave_idx_type j)
+{
+  octave_idx_type m = q.rows ();
+  octave_idx_type k = r.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (j < 0 || j > n-1)
+    (*current_liboctave_error_handler) ("qrdelete: index out of range");
+
+  OCTAVE_LOCAL_BUFFER (double, w, k);
+  F77_XFCN (dqrdec, DQRDEC, (m, n, k, q.fortran_vec (), q.rows (),
+                             r.fortran_vec (), r.rows (), j + 1, w));
+
+  if (k < m)
+    {
+      q.resize (m, k-1);
+      r.resize (k-1, n-1);
+    }
+  else
+    {
+      r.resize (k, n-1);
+    }
+}
+
+template <>
+void
+qr<Matrix>::delete_col (const Array<octave_idx_type>& j)
+{
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+  octave_idx_type k = q.columns ();
+
+  Array<octave_idx_type> jsi;
+  Array<octave_idx_type> js = j.sort (jsi, 0, DESCENDING);
+  octave_idx_type nj = js.numel ();
+  bool dups = false;
+  for (octave_idx_type i = 0; i < nj - 1; i++)
+    dups = dups && js(i) == js(i+1);
+
+  if (dups)
+    (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
+  if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0))
+    (*current_liboctave_error_handler) ("qrinsert: index out of range");
+
+  if (nj > 0)
+    {
+      OCTAVE_LOCAL_BUFFER (double, w, k);
+      for (volatile octave_idx_type i = 0; i < js.numel (); i++)
+        {
+          octave_idx_type ii = i;
+          F77_XFCN (dqrdec, DQRDEC, (m, n - ii, k == m ? k : k - ii,
+                                     q.fortran_vec (), q.rows (),
+                                     r.fortran_vec (), r.rows (),
+                                     js(ii) + 1, w));
+        }
+      if (k < m)
+        {
+          q.resize (m, k - nj);
+          r.resize (k - nj, n - nj);
+        }
+      else
+        {
+          r.resize (k, n - nj);
+        }
+
+    }
+}
+
+template <>
+void
+qr<Matrix>::insert_row (const RowVector& u, octave_idx_type j)
+{
+  octave_idx_type m = r.rows ();
+  octave_idx_type n = r.columns ();
+  octave_idx_type k = std::min (m, n);
+
+  if (! q.is_square () || u.numel () != n)
+    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
+  if (j < 0 || j > m)
+    (*current_liboctave_error_handler) ("qrinsert: index out of range");
+
+  q.resize (m + 1, m + 1);
+  r.resize (m + 1, n);
+  RowVector utmp = u;
+  OCTAVE_LOCAL_BUFFER (double, w, k);
+  F77_XFCN (dqrinr, DQRINR, (m, n, q.fortran_vec (), q.rows (),
+                             r.fortran_vec (), r.rows (),
+                             j + 1, utmp.fortran_vec (), w));
+
+}
+
+template <>
+void
+qr<Matrix>::delete_row (octave_idx_type j)
+{
+  octave_idx_type m = r.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (! q.is_square ())
+    (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch");
+  if (j < 0 || j > m-1)
+    (*current_liboctave_error_handler) ("qrdelete: index out of range");
+
+  OCTAVE_LOCAL_BUFFER (double, w, 2*m);
+  F77_XFCN (dqrder, DQRDER, (m, n, q.fortran_vec (), q.rows (),
+                             r.fortran_vec (), r.rows (), j + 1,
+                             w));
+
+  q.resize (m - 1, m - 1);
+  r.resize (m - 1, n);
+}
+
+template <>
+void
+qr<Matrix>::shift_cols (octave_idx_type i, octave_idx_type j)
+{
+  octave_idx_type m = q.rows ();
+  octave_idx_type k = r.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (i < 0 || i > n-1 || j < 0 || j > n-1)
+    (*current_liboctave_error_handler) ("qrshift: index out of range");
+
+  OCTAVE_LOCAL_BUFFER (double, w, 2*k);
+  F77_XFCN (dqrshc, DQRSHC, (m, n, k,
+                             q.fortran_vec (), q.rows (),
+                             r.fortran_vec (), r.rows (),
+                             i + 1, j + 1, w));
+}
+
+#endif
+
+template <>
+void
+qr<FloatMatrix>::form (octave_idx_type n, FloatMatrix& afact, float *tau, type qr_type)
+{
+  octave_idx_type m = afact.rows ();
+  octave_idx_type min_mn = std::min (m, n);
+  octave_idx_type info;
+
+  if (qr_type == qr<FloatMatrix>::raw)
+    {
+      for (octave_idx_type j = 0; j < min_mn; j++)
+        {
+          octave_idx_type limit = j < min_mn - 1 ? j : min_mn - 1;
+          for (octave_idx_type i = limit + 1; i < m; i++)
+            afact.elem (i, j) *= tau[j];
+        }
+
+      r = afact;
+    }
+  else
+    {
+      // Attempt to minimize copying.
+      if (m >= n)
+        {
+          // afact will become q.
+          q = afact;
+          octave_idx_type k = qr_type == qr<FloatMatrix>::economy ? n : m;
+          r = FloatMatrix (k, n);
+          for (octave_idx_type j = 0; j < n; j++)
+            {
+              octave_idx_type i = 0;
+              for (; i <= j; i++)
+                r.xelem (i, j) = afact.xelem (i, j);
+              for (; i < k; i++)
+                r.xelem (i, j) = 0;
+            }
+          afact = FloatMatrix (); // optimize memory
+        }
+      else
+        {
+          // afact will become r.
+          q = FloatMatrix (m, m);
+          for (octave_idx_type j = 0; j < m; j++)
+            for (octave_idx_type i = j + 1; i < m; i++)
+              {
+                q.xelem (i, j) = afact.xelem (i, j);
+                afact.xelem (i, j) = 0;
+              }
+          r = afact;
+        }
+
+
+      if (m > 0)
+        {
+          octave_idx_type k = q.columns ();
+          // workspace query.
+          float rlwork;
+          F77_XFCN (sorgqr, SORGQR, (m, k, min_mn, q.fortran_vec (), m, tau,
+                                     &rlwork, -1, info));
+
+          // allocate buffer and do the job.
+          octave_idx_type lwork = rlwork;
+          lwork = std::max (lwork, static_cast<octave_idx_type> (1));
+          OCTAVE_LOCAL_BUFFER (float, work, lwork);
+          F77_XFCN (sorgqr, SORGQR, (m, k, min_mn, q.fortran_vec (), m, tau,
+                                     work, lwork, info));
+        }
+    }
+}
+
+template <>
+void
+qr<FloatMatrix>::init (const FloatMatrix& a, type qr_type)
+{
+  octave_idx_type m = a.rows ();
+  octave_idx_type n = a.cols ();
+
+  octave_idx_type min_mn = m < n ? m : n;
+  OCTAVE_LOCAL_BUFFER (float, tau, min_mn);
+
+  octave_idx_type info = 0;
+
+  FloatMatrix afact = a;
+  if (m > n && qr_type == qr<FloatMatrix>::std)
+    afact.resize (m, m);
+
+  if (m > 0)
+    {
+      // workspace query.
+      float rlwork;
+      F77_XFCN (sgeqrf, SGEQRF, (m, n, afact.fortran_vec (), m, tau,
+                                 &rlwork, -1, info));
+
+      // allocate buffer and do the job.
+      octave_idx_type lwork = rlwork;
+      lwork = std::max (lwork, static_cast<octave_idx_type> (1));
+      OCTAVE_LOCAL_BUFFER (float, work, lwork);
+      F77_XFCN (sgeqrf, SGEQRF, (m, n, afact.fortran_vec (), m, tau,
+                                 work, lwork, info));
+    }
+
+  form (n, afact, tau, qr_type);
+}
+
+#ifdef HAVE_QRUPDATE
+
+template <>
+void
+qr<FloatMatrix>::update (const FloatColumnVector& u, const FloatColumnVector& v)
+{
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+  octave_idx_type k = q.columns ();
+
+  if (u.numel () != m || v.numel () != n)
+    (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
+
+  FloatColumnVector utmp = u;
+  FloatColumnVector vtmp = v;
+  OCTAVE_LOCAL_BUFFER (float, w, 2*k);
+  F77_XFCN (sqr1up, SQR1UP, (m, n, k, q.fortran_vec (),
+                             m, r.fortran_vec (), k,
+                             utmp.fortran_vec (), vtmp.fortran_vec (), w));
+}
+
+template <>
+void
+qr<FloatMatrix>::update (const FloatMatrix& u, const FloatMatrix& v)
+{
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+  octave_idx_type k = q.columns ();
+
+  if (u.rows () != m || v.rows () != n || u.cols () != v.cols ())
+    (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
+
+  OCTAVE_LOCAL_BUFFER (float, w, 2*k);
+  for (volatile octave_idx_type i = 0; i < u.cols (); i++)
+    {
+      FloatColumnVector utmp = u.column (i);
+      FloatColumnVector vtmp = v.column (i);
+      F77_XFCN (sqr1up, SQR1UP, (m, n, k, q.fortran_vec (),
+                                 m, r.fortran_vec (), k,
+                                 utmp.fortran_vec (), vtmp.fortran_vec (),
+                                 w));
+    }
+}
+
+template <>
+void
+qr<FloatMatrix>::insert_col (const FloatColumnVector& u, octave_idx_type j)
+{
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+  octave_idx_type k = q.columns ();
+
+  if (u.numel () != m)
+    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
+  if (j < 0 || j > n)
+    (*current_liboctave_error_handler) ("qrinsert: index out of range");
+
+  if (k < m)
+    {
+      q.resize (m, k+1);
+      r.resize (k+1, n+1);
+    }
+  else
+    {
+      r.resize (k, n+1);
+    }
+
+  FloatColumnVector utmp = u;
+  OCTAVE_LOCAL_BUFFER (float, w, k);
+  F77_XFCN (sqrinc, SQRINC, (m, n, k, q.fortran_vec (), q.rows (),
+                             r.fortran_vec (), r.rows (), j + 1,
+                             utmp.data (), w));
+}
+
+template <>
+void
+qr<FloatMatrix>::insert_col (const FloatMatrix& u, const Array<octave_idx_type>& j)
+{
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+  octave_idx_type k = q.columns ();
+
+  Array<octave_idx_type> jsi;
+  Array<octave_idx_type> js = j.sort (jsi, 0, ASCENDING);
+  octave_idx_type nj = js.numel ();
+  bool dups = false;
+  for (octave_idx_type i = 0; i < nj - 1; i++)
+    dups = dups && js(i) == js(i+1);
+
+  if (dups)
+    (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
+  if (u.numel () != m || u.columns () != nj)
+    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
+  if (nj > 0 && (js(0) < 0 || js(nj-1) > n))
+    (*current_liboctave_error_handler) ("qrinsert: index out of range");
+
+  if (nj > 0)
+    {
+      octave_idx_type kmax = std::min (k + nj, m);
+      if (k < m)
+        {
+          q.resize (m, kmax);
+          r.resize (kmax, n + nj);
+        }
+      else
+        {
+          r.resize (k, n + nj);
+        }
+
+      OCTAVE_LOCAL_BUFFER (float, w, kmax);
+      for (volatile octave_idx_type i = 0; i < js.numel (); i++)
+        {
+          octave_idx_type ii = i;
+          FloatColumnVector utmp = u.column (jsi(i));
+          F77_XFCN (sqrinc, SQRINC, (m, n + ii, std::min (kmax, k + ii),
+                                     q.fortran_vec (), q.rows (),
+                                     r.fortran_vec (), r.rows (), js(ii) + 1,
+                                     utmp.data (), w));
+        }
+    }
+}
+
+template <>
+void
+qr<FloatMatrix>::delete_col (octave_idx_type j)
+{
+  octave_idx_type m = q.rows ();
+  octave_idx_type k = r.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (j < 0 || j > n-1)
+    (*current_liboctave_error_handler) ("qrdelete: index out of range");
+
+  OCTAVE_LOCAL_BUFFER (float, w, k);
+  F77_XFCN (sqrdec, SQRDEC, (m, n, k, q.fortran_vec (), q.rows (),
+                             r.fortran_vec (), r.rows (), j + 1, w));
+
+  if (k < m)
+    {
+      q.resize (m, k-1);
+      r.resize (k-1, n-1);
+    }
+  else
+    {
+      r.resize (k, n-1);
+    }
+}
+
+template <>
+void
+qr<FloatMatrix>::delete_col (const Array<octave_idx_type>& j)
+{
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+  octave_idx_type k = q.columns ();
+
+  Array<octave_idx_type> jsi;
+  Array<octave_idx_type> js = j.sort (jsi, 0, DESCENDING);
+  octave_idx_type nj = js.numel ();
+  bool dups = false;
+  for (octave_idx_type i = 0; i < nj - 1; i++)
+    dups = dups && js(i) == js(i+1);
+
+  if (dups)
+    (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
+  if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0))
+    (*current_liboctave_error_handler) ("qrinsert: index out of range");
+
+  if (nj > 0)
+    {
+      OCTAVE_LOCAL_BUFFER (float, w, k);
+      for (volatile octave_idx_type i = 0; i < js.numel (); i++)
+        {
+          octave_idx_type ii = i;
+          F77_XFCN (sqrdec, SQRDEC, (m, n - ii, k == m ? k : k - ii,
+                                     q.fortran_vec (), q.rows (),
+                                     r.fortran_vec (), r.rows (),
+                                     js(ii) + 1, w));
+        }
+      if (k < m)
+        {
+          q.resize (m, k - nj);
+          r.resize (k - nj, n - nj);
+        }
+      else
+        {
+          r.resize (k, n - nj);
+        }
+
+    }
+}
+
+template <>
+void
+qr<FloatMatrix>::insert_row (const FloatRowVector& u, octave_idx_type j)
+{
+  octave_idx_type m = r.rows ();
+  octave_idx_type n = r.columns ();
+  octave_idx_type k = std::min (m, n);
+
+  if (! q.is_square () || u.numel () != n)
+    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
+  if (j < 0 || j > m)
+    (*current_liboctave_error_handler) ("qrinsert: index out of range");
+
+  q.resize (m + 1, m + 1);
+  r.resize (m + 1, n);
+  FloatRowVector utmp = u;
+  OCTAVE_LOCAL_BUFFER (float, w, k);
+  F77_XFCN (sqrinr, SQRINR, (m, n, q.fortran_vec (), q.rows (),
+                             r.fortran_vec (), r.rows (),
+                             j + 1, utmp.fortran_vec (), w));
+
+}
+
+template <>
+void
+qr<FloatMatrix>::delete_row (octave_idx_type j)
+{
+  octave_idx_type m = r.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (! q.is_square ())
+    (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch");
+  if (j < 0 || j > m-1)
+    (*current_liboctave_error_handler) ("qrdelete: index out of range");
+
+  OCTAVE_LOCAL_BUFFER (float, w, 2*m);
+  F77_XFCN (sqrder, SQRDER, (m, n, q.fortran_vec (), q.rows (),
+                             r.fortran_vec (), r.rows (), j + 1,
+                             w));
+
+  q.resize (m - 1, m - 1);
+  r.resize (m - 1, n);
+}
+
+template <>
+void
+qr<FloatMatrix>::shift_cols (octave_idx_type i, octave_idx_type j)
+{
+  octave_idx_type m = q.rows ();
+  octave_idx_type k = r.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (i < 0 || i > n-1 || j < 0 || j > n-1)
+    (*current_liboctave_error_handler) ("qrshift: index out of range");
+
+  OCTAVE_LOCAL_BUFFER (float, w, 2*k);
+  F77_XFCN (sqrshc, SQRSHC, (m, n, k,
+                             q.fortran_vec (), q.rows (),
+                             r.fortran_vec (), r.rows (),
+                             i + 1, j + 1, w));
+}
+
+#endif
+
+template <>
+void
+qr<ComplexMatrix>::form (octave_idx_type n, ComplexMatrix& afact,
+                 Complex *tau, type qr_type)
+{
+  octave_idx_type m = afact.rows ();
+  octave_idx_type min_mn = std::min (m, n);
+  octave_idx_type info;
+
+  if (qr_type == qr<ComplexMatrix>::raw)
+    {
+      for (octave_idx_type j = 0; j < min_mn; j++)
+        {
+          octave_idx_type limit = j < min_mn - 1 ? j : min_mn - 1;
+          for (octave_idx_type i = limit + 1; i < m; i++)
+            afact.elem (i, j) *= tau[j];
+        }
+
+      r = afact;
+    }
+  else
+    {
+      // Attempt to minimize copying.
+      if (m >= n)
+        {
+          // afact will become q.
+          q = afact;
+          octave_idx_type k = qr_type == qr<ComplexMatrix>::economy ? n : m;
+          r = ComplexMatrix (k, n);
+          for (octave_idx_type j = 0; j < n; j++)
+            {
+              octave_idx_type i = 0;
+              for (; i <= j; i++)
+                r.xelem (i, j) = afact.xelem (i, j);
+              for (; i < k; i++)
+                r.xelem (i, j) = 0;
+            }
+          afact = ComplexMatrix (); // optimize memory
+        }
+      else
+        {
+          // afact will become r.
+          q = ComplexMatrix (m, m);
+          for (octave_idx_type j = 0; j < m; j++)
+            for (octave_idx_type i = j + 1; i < m; i++)
+              {
+                q.xelem (i, j) = afact.xelem (i, j);
+                afact.xelem (i, j) = 0;
+              }
+          r = afact;
+        }
+
+
+      if (m > 0)
+        {
+          octave_idx_type k = q.columns ();
+          // workspace query.
+          Complex clwork;
+          F77_XFCN (zungqr, ZUNGQR, (m, k, min_mn, q.fortran_vec (), m, tau,
+                                     &clwork, -1, info));
+
+          // allocate buffer and do the job.
+          octave_idx_type lwork = clwork.real ();
+          lwork = std::max (lwork, static_cast<octave_idx_type> (1));
+          OCTAVE_LOCAL_BUFFER (Complex, work, lwork);
+          F77_XFCN (zungqr, ZUNGQR, (m, k, min_mn, q.fortran_vec (), m, tau,
+                                     work, lwork, info));
+        }
+    }
+}
+
+template <>
+void
+qr<ComplexMatrix>::init (const ComplexMatrix& a, type qr_type)
+{
+  octave_idx_type m = a.rows ();
+  octave_idx_type n = a.cols ();
+
+  octave_idx_type min_mn = m < n ? m : n;
+  OCTAVE_LOCAL_BUFFER (Complex, tau, min_mn);
+
+  octave_idx_type info = 0;
+
+  ComplexMatrix afact = a;
+  if (m > n && qr_type == qr<ComplexMatrix>::std)
+    afact.resize (m, m);
+
+  if (m > 0)
+    {
+      // workspace query.
+      Complex clwork;
+      F77_XFCN (zgeqrf, ZGEQRF, (m, n, afact.fortran_vec (), m, tau,
+                                 &clwork, -1, info));
+
+      // allocate buffer and do the job.
+      octave_idx_type lwork = clwork.real ();
+      lwork = std::max (lwork, static_cast<octave_idx_type> (1));
+      OCTAVE_LOCAL_BUFFER (Complex, work, lwork);
+      F77_XFCN (zgeqrf, ZGEQRF, (m, n, afact.fortran_vec (), m, tau,
+                                 work, lwork, info));
+    }
+
+  form (n, afact, tau, qr_type);
+}
+
+#if defined (HAVE_QRUPDATE)
+
+template <>
+void
+qr<ComplexMatrix>::update (const ComplexColumnVector& u, const ComplexColumnVector& v)
+{
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+  octave_idx_type k = q.columns ();
+
+  if (u.numel () != m || v.numel () != n)
+    (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
+
+  ComplexColumnVector utmp = u;
+  ComplexColumnVector vtmp = v;
+  OCTAVE_LOCAL_BUFFER (Complex, w, k);
+  OCTAVE_LOCAL_BUFFER (double, rw, k);
+  F77_XFCN (zqr1up, ZQR1UP, (m, n, k, q.fortran_vec (),
+                             m, r.fortran_vec (), k,
+                             utmp.fortran_vec (), vtmp.fortran_vec (),
+                             w, rw));
+}
+
+template <>
+void
+qr<ComplexMatrix>::update (const ComplexMatrix& u, const ComplexMatrix& v)
+{
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+  octave_idx_type k = q.columns ();
+
+  if (u.rows () != m || v.rows () != n || u.cols () != v.cols ())
+    (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
+
+  OCTAVE_LOCAL_BUFFER (Complex, w, k);
+  OCTAVE_LOCAL_BUFFER (double, rw, k);
+  for (volatile octave_idx_type i = 0; i < u.cols (); i++)
+    {
+      ComplexColumnVector utmp = u.column (i);
+      ComplexColumnVector vtmp = v.column (i);
+      F77_XFCN (zqr1up, ZQR1UP, (m, n, k, q.fortran_vec (),
+                                 m, r.fortran_vec (), k,
+                                 utmp.fortran_vec (), vtmp.fortran_vec (),
+                                 w, rw));
+    }
+}
+
+template <>
+void
+qr<ComplexMatrix>::insert_col (const ComplexColumnVector& u, octave_idx_type j)
+{
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+  octave_idx_type k = q.columns ();
+
+  if (u.numel () != m)
+    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
+  if (j < 0 || j > n)
+    (*current_liboctave_error_handler) ("qrinsert: index out of range");
+
+  if (k < m)
+    {
+      q.resize (m, k+1);
+      r.resize (k+1, n+1);
+    }
+  else
+    {
+      r.resize (k, n+1);
+    }
+
+  ComplexColumnVector utmp = u;
+  OCTAVE_LOCAL_BUFFER (double, rw, k);
+  F77_XFCN (zqrinc, ZQRINC, (m, n, k, q.fortran_vec (), q.rows (),
+                             r.fortran_vec (), r.rows (), j + 1,
+                             utmp.data (), rw));
+}
+
+template <>
+void
+qr<ComplexMatrix>::insert_col (const ComplexMatrix& u, const Array<octave_idx_type>& j)
+{
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+  octave_idx_type k = q.columns ();
+
+  Array<octave_idx_type> jsi;
+  Array<octave_idx_type> js = j.sort (jsi, 0, ASCENDING);
+  octave_idx_type nj = js.numel ();
+  bool dups = false;
+  for (octave_idx_type i = 0; i < nj - 1; i++)
+    dups = dups && js(i) == js(i+1);
+
+  if (dups)
+    (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
+  if (u.numel () != m || u.columns () != nj)
+    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
+  if (nj > 0 && (js(0) < 0 || js(nj-1) > n))
+    (*current_liboctave_error_handler) ("qrinsert: index out of range");
+
+  if (nj > 0)
+    {
+      octave_idx_type kmax = std::min (k + nj, m);
+      if (k < m)
+        {
+          q.resize (m, kmax);
+          r.resize (kmax, n + nj);
+        }
+      else
+        {
+          r.resize (k, n + nj);
+        }
+
+      OCTAVE_LOCAL_BUFFER (double, rw, kmax);
+      for (volatile octave_idx_type i = 0; i < js.numel (); i++)
+        {
+          octave_idx_type ii = i;
+          ComplexColumnVector utmp = u.column (jsi(i));
+          F77_XFCN (zqrinc, ZQRINC, (m, n + ii, std::min (kmax, k + ii),
+                                     q.fortran_vec (), q.rows (),
+                                     r.fortran_vec (), r.rows (), js(ii) + 1,
+                                     utmp.data (), rw));
+        }
+    }
+}
+
+template <>
+void
+qr<ComplexMatrix>::delete_col (octave_idx_type j)
+{
+  octave_idx_type m = q.rows ();
+  octave_idx_type k = r.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (j < 0 || j > n-1)
+    (*current_liboctave_error_handler) ("qrdelete: index out of range");
+
+  OCTAVE_LOCAL_BUFFER (double, rw, k);
+  F77_XFCN (zqrdec, ZQRDEC, (m, n, k, q.fortran_vec (), q.rows (),
+                             r.fortran_vec (), r.rows (), j + 1, rw));
+
+  if (k < m)
+    {
+      q.resize (m, k-1);
+      r.resize (k-1, n-1);
+    }
+  else
+    {
+      r.resize (k, n-1);
+    }
+}
+
+template <>
+void
+qr<ComplexMatrix>::delete_col (const Array<octave_idx_type>& j)
+{
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+  octave_idx_type k = q.columns ();
+
+  Array<octave_idx_type> jsi;
+  Array<octave_idx_type> js = j.sort (jsi, 0, DESCENDING);
+  octave_idx_type nj = js.numel ();
+  bool dups = false;
+  for (octave_idx_type i = 0; i < nj - 1; i++)
+    dups = dups && js(i) == js(i+1);
+
+  if (dups)
+    (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
+  if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0))
+    (*current_liboctave_error_handler) ("qrinsert: index out of range");
+
+  if (nj > 0)
+    {
+      OCTAVE_LOCAL_BUFFER (double, rw, k);
+      for (volatile octave_idx_type i = 0; i < js.numel (); i++)
+        {
+          octave_idx_type ii = i;
+          F77_XFCN (zqrdec, ZQRDEC, (m, n - ii, k == m ? k : k - ii,
+                                     q.fortran_vec (), q.rows (),
+                                     r.fortran_vec (), r.rows (),
+                                     js(ii) + 1, rw));
+        }
+      if (k < m)
+        {
+          q.resize (m, k - nj);
+          r.resize (k - nj, n - nj);
+        }
+      else
+        {
+          r.resize (k, n - nj);
+        }
+
+    }
+}
+
+template <>
+void
+qr<ComplexMatrix>::insert_row (const ComplexRowVector& u, octave_idx_type j)
+{
+  octave_idx_type m = r.rows ();
+  octave_idx_type n = r.columns ();
+  octave_idx_type k = std::min (m, n);
+
+  if (! q.is_square () || u.numel () != n)
+    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
+  if (j < 0 || j > m)
+    (*current_liboctave_error_handler) ("qrinsert: index out of range");
+
+  q.resize (m + 1, m + 1);
+  r.resize (m + 1, n);
+  ComplexRowVector utmp = u;
+  OCTAVE_LOCAL_BUFFER (double, rw, k);
+  F77_XFCN (zqrinr, ZQRINR, (m, n, q.fortran_vec (), q.rows (),
+                             r.fortran_vec (), r.rows (),
+                             j + 1, utmp.fortran_vec (), rw));
+
+}
+
+template <>
+void
+qr<ComplexMatrix>::delete_row (octave_idx_type j)
+{
+  octave_idx_type m = r.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (! q.is_square ())
+    (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch");
+  if (j < 0 || j > m-1)
+    (*current_liboctave_error_handler) ("qrdelete: index out of range");
+
+  OCTAVE_LOCAL_BUFFER (Complex, w, m);
+  OCTAVE_LOCAL_BUFFER (double, rw, m);
+  F77_XFCN (zqrder, ZQRDER, (m, n, q.fortran_vec (), q.rows (),
+                             r.fortran_vec (), r.rows (), j + 1,
+                             w, rw));
+
+  q.resize (m - 1, m - 1);
+  r.resize (m - 1, n);
+}
+
+template <>
+void
+qr<ComplexMatrix>::shift_cols (octave_idx_type i, octave_idx_type j)
+{
+  octave_idx_type m = q.rows ();
+  octave_idx_type k = r.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (i < 0 || i > n-1 || j < 0 || j > n-1)
+    (*current_liboctave_error_handler) ("qrshift: index out of range");
+
+  OCTAVE_LOCAL_BUFFER (Complex, w, k);
+  OCTAVE_LOCAL_BUFFER (double, rw, k);
+  F77_XFCN (zqrshc, ZQRSHC, (m, n, k,
+                             q.fortran_vec (), q.rows (),
+                             r.fortran_vec (), r.rows (),
+                             i + 1, j + 1, w, rw));
+}
+
+#endif
+
+template <>
+void
+qr<FloatComplexMatrix>::form (octave_idx_type n, FloatComplexMatrix& afact, FloatComplex *tau, type qr_type)
+{
+  octave_idx_type m = afact.rows ();
+  octave_idx_type min_mn = std::min (m, n);
+  octave_idx_type info;
+
+  if (qr_type == qr<FloatComplexMatrix>::raw)
+    {
+      for (octave_idx_type j = 0; j < min_mn; j++)
+        {
+          octave_idx_type limit = j < min_mn - 1 ? j : min_mn - 1;
+          for (octave_idx_type i = limit + 1; i < m; i++)
+            afact.elem (i, j) *= tau[j];
+        }
+
+      r = afact;
+    }
+  else
+    {
+      // Attempt to minimize copying.
+      if (m >= n)
+        {
+          // afact will become q.
+          q = afact;
+          octave_idx_type k = qr_type == qr<FloatComplexMatrix>::economy ? n : m;
+          r = FloatComplexMatrix (k, n);
+          for (octave_idx_type j = 0; j < n; j++)
+            {
+              octave_idx_type i = 0;
+              for (; i <= j; i++)
+                r.xelem (i, j) = afact.xelem (i, j);
+              for (; i < k; i++)
+                r.xelem (i, j) = 0;
+            }
+          afact = FloatComplexMatrix (); // optimize memory
+        }
+      else
+        {
+          // afact will become r.
+          q = FloatComplexMatrix (m, m);
+          for (octave_idx_type j = 0; j < m; j++)
+            for (octave_idx_type i = j + 1; i < m; i++)
+              {
+                q.xelem (i, j) = afact.xelem (i, j);
+                afact.xelem (i, j) = 0;
+              }
+          r = afact;
+        }
+
+
+      if (m > 0)
+        {
+          octave_idx_type k = q.columns ();
+          // workspace query.
+          FloatComplex clwork;
+          F77_XFCN (cungqr, CUNGQR, (m, k, min_mn, q.fortran_vec (), m, tau,
+                                     &clwork, -1, info));
+
+          // allocate buffer and do the job.
+          octave_idx_type lwork = clwork.real ();
+          lwork = std::max (lwork, static_cast<octave_idx_type> (1));
+          OCTAVE_LOCAL_BUFFER (FloatComplex, work, lwork);
+          F77_XFCN (cungqr, CUNGQR, (m, k, min_mn, q.fortran_vec (), m, tau,
+                                     work, lwork, info));
+        }
+    }
+}
+
+template <>
+void
+qr<FloatComplexMatrix>::init (const FloatComplexMatrix& a, type qr_type)
+{
+  octave_idx_type m = a.rows ();
+  octave_idx_type n = a.cols ();
+
+  octave_idx_type min_mn = m < n ? m : n;
+  OCTAVE_LOCAL_BUFFER (FloatComplex, tau, min_mn);
+
+  octave_idx_type info = 0;
+
+  FloatComplexMatrix afact = a;
+  if (m > n && qr_type == qr<FloatComplexMatrix>::std)
+    afact.resize (m, m);
+
+  if (m > 0)
+    {
+      // workspace query.
+      FloatComplex clwork;
+      F77_XFCN (cgeqrf, CGEQRF, (m, n, afact.fortran_vec (), m, tau,
+                                 &clwork, -1, info));
+
+      // allocate buffer and do the job.
+      octave_idx_type lwork = clwork.real ();
+      lwork = std::max (lwork, static_cast<octave_idx_type> (1));
+      OCTAVE_LOCAL_BUFFER (FloatComplex, work, lwork);
+      F77_XFCN (cgeqrf, CGEQRF, (m, n, afact.fortran_vec (), m, tau,
+                                 work, lwork, info));
+    }
+
+  form (n, afact, tau, qr_type);
+}
+
+#ifdef HAVE_QRUPDATE
+
+template <>
+void
+qr<FloatComplexMatrix>::update (const FloatComplexColumnVector& u, const FloatComplexColumnVector& v)
+{
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+  octave_idx_type k = q.columns ();
+
+  if (u.numel () != m || v.numel () != n)
+    (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
+
+  FloatComplexColumnVector utmp = u;
+  FloatComplexColumnVector vtmp = v;
+  OCTAVE_LOCAL_BUFFER (FloatComplex, w, k);
+  OCTAVE_LOCAL_BUFFER (float, rw, k);
+  F77_XFCN (cqr1up, CQR1UP, (m, n, k, q.fortran_vec (),
+                             m, r.fortran_vec (), k,
+                             utmp.fortran_vec (), vtmp.fortran_vec (),
+                             w, rw));
+}
+
+template <>
+void
+qr<FloatComplexMatrix>::update (const FloatComplexMatrix& u, const FloatComplexMatrix& v)
+{
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+  octave_idx_type k = q.columns ();
+
+  if (u.rows () != m || v.rows () != n || u.cols () != v.cols ())
+    (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
+
+  OCTAVE_LOCAL_BUFFER (FloatComplex, w, k);
+  OCTAVE_LOCAL_BUFFER (float, rw, k);
+  for (volatile octave_idx_type i = 0; i < u.cols (); i++)
+    {
+      FloatComplexColumnVector utmp = u.column (i);
+      FloatComplexColumnVector vtmp = v.column (i);
+      F77_XFCN (cqr1up, CQR1UP, (m, n, k, q.fortran_vec (),
+                                 m, r.fortran_vec (), k,
+                                 utmp.fortran_vec (), vtmp.fortran_vec (),
+                                 w, rw));
+    }
+}
+
+template <>
+void
+qr<FloatComplexMatrix>::insert_col (const FloatComplexColumnVector& u, octave_idx_type j)
+{
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+  octave_idx_type k = q.columns ();
+
+  if (u.numel () != m)
+    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
+  if (j < 0 || j > n)
+    (*current_liboctave_error_handler) ("qrinsert: index out of range");
+
+  if (k < m)
+    {
+      q.resize (m, k+1);
+      r.resize (k+1, n+1);
+    }
+  else
+    {
+      r.resize (k, n+1);
+    }
+
+  FloatComplexColumnVector utmp = u;
+  OCTAVE_LOCAL_BUFFER (float, rw, k);
+  F77_XFCN (cqrinc, CQRINC, (m, n, k, q.fortran_vec (), q.rows (),
+                             r.fortran_vec (), r.rows (), j + 1,
+                             utmp.data (), rw));
+}
+
+template <>
+void
+qr<FloatComplexMatrix>::insert_col (const FloatComplexMatrix& u, const Array<octave_idx_type>& j)
+{
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+  octave_idx_type k = q.columns ();
+
+  Array<octave_idx_type> jsi;
+  Array<octave_idx_type> js = j.sort (jsi, 0, ASCENDING);
+  octave_idx_type nj = js.numel ();
+  bool dups = false;
+  for (octave_idx_type i = 0; i < nj - 1; i++)
+    dups = dups && js(i) == js(i+1);
+
+  if (dups)
+    (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
+  if (u.numel () != m || u.columns () != nj)
+    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
+  if (nj > 0 && (js(0) < 0 || js(nj-1) > n))
+    (*current_liboctave_error_handler) ("qrinsert: index out of range");
+
+  if (nj > 0)
+    {
+      octave_idx_type kmax = std::min (k + nj, m);
+      if (k < m)
+        {
+          q.resize (m, kmax);
+          r.resize (kmax, n + nj);
+        }
+      else
+        {
+          r.resize (k, n + nj);
+        }
+
+      OCTAVE_LOCAL_BUFFER (float, rw, kmax);
+      for (volatile octave_idx_type i = 0; i < js.numel (); i++)
+        {
+          octave_idx_type ii = i;
+          F77_XFCN (cqrinc, CQRINC, (m, n + ii, std::min (kmax, k + ii),
+                                     q.fortran_vec (), q.rows (),
+                                     r.fortran_vec (), r.rows (), js(ii) + 1,
+                                     u.column (jsi(i)).data (), rw));
+        }
+    }
+}
+
+template <>
+void
+qr<FloatComplexMatrix>::delete_col (octave_idx_type j)
+{
+  octave_idx_type m = q.rows ();
+  octave_idx_type k = r.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (j < 0 || j > n-1)
+    (*current_liboctave_error_handler) ("qrdelete: index out of range");
+
+  OCTAVE_LOCAL_BUFFER (float, rw, k);
+  F77_XFCN (cqrdec, CQRDEC, (m, n, k, q.fortran_vec (), q.rows (),
+                             r.fortran_vec (), r.rows (), j + 1, rw));
+
+  if (k < m)
+    {
+      q.resize (m, k-1);
+      r.resize (k-1, n-1);
+    }
+  else
+    {
+      r.resize (k, n-1);
+    }
+}
+
+template <>
+void
+qr<FloatComplexMatrix>::delete_col (const Array<octave_idx_type>& j)
+{
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+  octave_idx_type k = q.columns ();
+
+  Array<octave_idx_type> jsi;
+  Array<octave_idx_type> js = j.sort (jsi, 0, DESCENDING);
+  octave_idx_type nj = js.numel ();
+  bool dups = false;
+  for (octave_idx_type i = 0; i < nj - 1; i++)
+    dups = dups && js(i) == js(i+1);
+
+  if (dups)
+    (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
+  if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0))
+    (*current_liboctave_error_handler) ("qrinsert: index out of range");
+
+  if (nj > 0)
+    {
+      OCTAVE_LOCAL_BUFFER (float, rw, k);
+      for (volatile octave_idx_type i = 0; i < js.numel (); i++)
+        {
+          octave_idx_type ii = i;
+          F77_XFCN (cqrdec, CQRDEC, (m, n - ii, k == m ? k : k - ii,
+                                     q.fortran_vec (), q.rows (),
+                                     r.fortran_vec (), r.rows (),
+                                     js(ii) + 1, rw));
+        }
+      if (k < m)
+        {
+          q.resize (m, k - nj);
+          r.resize (k - nj, n - nj);
+        }
+      else
+        {
+          r.resize (k, n - nj);
+        }
+
+    }
+}
+
+template <>
+void
+qr<FloatComplexMatrix>::insert_row (const FloatComplexRowVector& u, octave_idx_type j)
+{
+  octave_idx_type m = r.rows ();
+  octave_idx_type n = r.columns ();
+  octave_idx_type k = std::min (m, n);
+
+  if (! q.is_square () || u.numel () != n)
+    (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
+  if (j < 0 || j > m)
+    (*current_liboctave_error_handler) ("qrinsert: index out of range");
+
+  q.resize (m + 1, m + 1);
+  r.resize (m + 1, n);
+  FloatComplexRowVector utmp = u;
+  OCTAVE_LOCAL_BUFFER (float, rw, k);
+  F77_XFCN (cqrinr, CQRINR, (m, n, q.fortran_vec (), q.rows (),
+                             r.fortran_vec (), r.rows (),
+                             j + 1, utmp.fortran_vec (), rw));
+
+}
+
+template <>
+void
+qr<FloatComplexMatrix>::delete_row (octave_idx_type j)
+{
+  octave_idx_type m = r.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (! q.is_square ())
+    (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch");
+  if (j < 0 || j > m-1)
+    (*current_liboctave_error_handler) ("qrdelete: index out of range");
+
+  OCTAVE_LOCAL_BUFFER (FloatComplex, w, m);
+  OCTAVE_LOCAL_BUFFER (float, rw, m);
+  F77_XFCN (cqrder, CQRDER, (m, n, q.fortran_vec (), q.rows (),
+                             r.fortran_vec (), r.rows (), j + 1,
+                             w, rw));
+
+  q.resize (m - 1, m - 1);
+  r.resize (m - 1, n);
+}
+
+template <>
+void
+qr<FloatComplexMatrix>::shift_cols (octave_idx_type i, octave_idx_type j)
+{
+  octave_idx_type m = q.rows ();
+  octave_idx_type k = r.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (i < 0 || i > n-1 || j < 0 || j > n-1)
+    (*current_liboctave_error_handler) ("qrshift: index out of range");
+
+  OCTAVE_LOCAL_BUFFER (FloatComplex, w, k);
+  OCTAVE_LOCAL_BUFFER (float, rw, k);
+  F77_XFCN (cqrshc, CQRSHC, (m, n, k,
+                             q.fortran_vec (), q.rows (),
+                             r.fortran_vec (), r.rows (),
+                             i + 1, j + 1, w, rw));
+}
+
+#endif
+
+// Instantiations we need.
+
+template class qr<Matrix>;
+
+template class qr<FloatMatrix>;
+
+template class qr<ComplexMatrix>;
+
+template class qr<FloatComplexMatrix>;
diff -r 48d82f74243e -r eb1524b07fe3 liboctave/numeric/qr.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/numeric/qr.h	Wed Feb 17 02:57:21 2016 -0500
@@ -0,0 +1,112 @@
+/*
+
+Copyright (C) 1994-2015 John W. Eaton
+Copyright (C) 2008-2009 Jaroslav Hajek
+Copyright (C) 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/>.
+
+*/
+
+#if ! defined (octave_qr_h)
+#define octave_qr_h 1
+
+#include "octave-config.h"
+
+#include "Array.h"
+
+template <typename T>
+class
+qr
+{
+public:
+
+  typedef typename T::element_type ELT_T;
+  typedef typename T::row_vector_type RV_T;
+  typedef typename T::column_vector_type CV_T;
+
+  enum type
+  {
+    std,
+    raw,
+    economy
+  };
+
+  qr (void) : q (), r () { }
+
+  qr (const T& a, type qr_type = qr::std)
+    : q (), r ()
+  {
+    init (a, qr_type);
+  }
+
+  qr (const T& q, const T& r);
+
+  qr (const qr& a) : q (a.q), r (a.r) { }
+
+  qr& operator = (const qr& a)
+  {
+    if (this != &a)
+      {
+        q = a.q;
+        r = a.r;
+      }
+
+    return *this;
+  }
+
+  virtual ~qr (void) { }
+
+  T Q (void) const { return q; }
+
+  T R (void) const { return r; }
+
+  type get_type (void) const;
+
+  bool regular (void) const;
+
+  void init (const T& a, type qr_type);
+
+  void update (const CV_T& u, const CV_T& v);
+
+  void update (const T& u, const T& v);
+
+  void insert_col (const CV_T& u, octave_idx_type j);
+
+  void insert_col (const T& u, const Array<octave_idx_type>& j);
+
+  void delete_col (octave_idx_type j);
+
+  void delete_col (const Array<octave_idx_type>& j);
+
+  void insert_row (const RV_T& u, octave_idx_type j);
+
+  void delete_row (octave_idx_type j);
+
+  void shift_cols (octave_idx_type i, octave_idx_type j);
+
+protected:
+
+  T q;
+  T r;
+
+  void form (octave_idx_type n, T& afact, ELT_T *tau, type qr_type);
+};
+
+extern void warn_qrupdate_once (void);
+
+#endif
diff -r 48d82f74243e -r eb1524b07fe3 liboctave/operators/mx-defs.h
--- a/liboctave/operators/mx-defs.h	Tue Feb 16 23:34:44 2016 -0500
+++ b/liboctave/operators/mx-defs.h	Wed Feb 17 02:57:21 2016 -0500
@@ -74,10 +74,7 @@
 
 template <typename T> class lu;
 
-class QR;
-class ComplexQR;
-class FloatQR;
-class FloatComplexQR;
+template <typename T> class qr;
 
 class QRP;
 class ComplexQRP;
diff -r 48d82f74243e -r eb1524b07fe3 liboctave/operators/mx-ext.h
--- a/liboctave/operators/mx-ext.h	Tue Feb 16 23:34:44 2016 -0500
+++ b/liboctave/operators/mx-ext.h	Wed Feb 17 02:57:21 2016 -0500
@@ -63,10 +63,6 @@
 
 // Result of a QR decomposition.
 
-#include "dbleQR.h"
-#include "CmplxQR.h"
-
-#include "dbleQRP.h"
-#include "CmplxQRP.h"
+#include "qr.h"
 
 #endif