# HG changeset patch
# User John W. Eaton <jwe@octave.org>
# Date 1455608849 18000
# Node ID 3c8a3d35661a6422812b35416970181128a5d12c
# Parent  f08ae27289e40c4732fff014d825d3ff95a86011
better use of templates for Cholesky factorization

* liboctave/numeric/chol.h, liboctave/numeric/chol.cc: New files
generated from CmplxCHOL.cc, fCmplxCHOL.cc, floatCHOL.cc, CmplxCHOL.h,
dbleCHOL.cc, dbleCHOL.h, fCmplxCHOL.h, and floatCHOL.h and converted
to templates.

* liboctave/numeric/module.mk: Update.

* __qp__.cc, chol.cc, CMatrix.cc, CMatrix.h, dMatrix.cc, dMatrix.h,
fCMatrix.cc, fCMatrix.h, fMatrix.cc, fMatrix.h, eigs-base.cc,
mx-defs.h, mx-ext.h: Use new classes.

diff -r f08ae27289e4 -r 3c8a3d35661a libinterp/corefcn/__qp__.cc
--- a/libinterp/corefcn/__qp__.cc	Tue Feb 16 00:32:29 2016 -0500
+++ b/libinterp/corefcn/__qp__.cc	Tue Feb 16 02:47:29 2016 -0500
@@ -26,7 +26,7 @@
 
 #include <cfloat>
 
-#include "dbleCHOL.h"
+#include "chol.h"
 #include "dbleSVD.h"
 #include "mx-m-dm.h"
 #include "EIG.h"
@@ -191,7 +191,7 @@
               // factorization since the Hessian is positive
               // definite.
 
-              CHOL cholH (H);
+              chol<Matrix> cholH (H);
 
               R = cholH.chol_matrix ();
 
@@ -250,7 +250,7 @@
               // Computing the Cholesky factorization (pR = 0 means
               // that the reduced Hessian was positive definite).
 
-              CHOL cholrH (rH, pR);
+              chol<Matrix> cholrH (rH, pR);
               Matrix tR = cholrH.chol_matrix ();
               if (pR == 0)
                 R = tR;
diff -r f08ae27289e4 -r 3c8a3d35661a libinterp/dldfcn/chol.cc
--- a/libinterp/dldfcn/chol.cc	Tue Feb 16 00:32:29 2016 -0500
+++ b/libinterp/dldfcn/chol.cc	Tue Feb 16 02:47:29 2016 -0500
@@ -27,10 +27,7 @@
 #  include <config.h>
 #endif
 
-#include "CmplxCHOL.h"
-#include "dbleCHOL.h"
-#include "fCmplxCHOL.h"
-#include "floatCHOL.h"
+#include "chol.h"
 #include "sparse-chol.h"
 #include "oct-spparms.h"
 #include "sparse-util.h"
@@ -250,8 +247,7 @@
 
           octave_idx_type info;
 
-          FloatCHOL fact;
-          fact = FloatCHOL (m, info, LLt != true);
+          chol<FloatMatrix> fact (m, info, LLt != true);
 
           if (nargout == 2 || info == 0)
             retval = ovl (get_chol (fact), info);
@@ -264,8 +260,7 @@
 
           octave_idx_type info;
 
-          FloatComplexCHOL fact;
-          fact = FloatComplexCHOL (m, info, LLt != true);
+          chol<FloatComplexMatrix> fact (m, info, LLt != true);
 
           if (nargout == 2 || info == 0)
             retval = ovl (get_chol (fact), info);
@@ -283,8 +278,7 @@
 
           octave_idx_type info;
 
-          CHOL fact;
-          fact = CHOL (m, info, LLt != true);
+          chol<Matrix> fact (m, info, LLt != true);
 
           if (nargout == 2 || info == 0)
             retval = ovl (get_chol (fact), info);
@@ -297,8 +291,7 @@
 
           octave_idx_type info;
 
-          ComplexCHOL fact;
-          fact = ComplexCHOL (m, info, LLt != true);
+          chol<ComplexMatrix> fact (m, info, LLt != true);
 
           if (nargout == 2 || info == 0)
             retval = ovl (get_chol (fact), info);
@@ -385,7 +378,7 @@
               FloatMatrix m = arg.float_matrix_value ();
 
               octave_idx_type info;
-              FloatCHOL chol (m, info);
+              chol<FloatMatrix> chol (m, info);
               if (info == 0)
                 retval = chol.inverse ();
               else
@@ -396,7 +389,7 @@
               FloatComplexMatrix m = arg.float_complex_matrix_value ();
 
               octave_idx_type info;
-              FloatComplexCHOL chol (m, info);
+              chol<FloatComplexMatrix> chol (m, info);
               if (info == 0)
                 retval = chol.inverse ();
               else
@@ -412,7 +405,7 @@
               Matrix m = arg.matrix_value ();
 
               octave_idx_type info;
-              CHOL chol (m, info);
+              chol<Matrix> chol (m, info);
               if (info == 0)
                 retval = chol.inverse ();
               else
@@ -423,7 +416,7 @@
               ComplexMatrix m = arg.complex_matrix_value ();
 
               octave_idx_type info;
-              ComplexCHOL chol (m, info);
+              chol<ComplexMatrix> chol (m, info);
               if (info == 0)
                 retval = chol.inverse ();
               else
@@ -602,7 +595,7 @@
           FloatMatrix R = argr.float_matrix_value ();
           FloatColumnVector u = argu.float_column_vector_value ();
 
-          FloatCHOL fact;
+          chol<FloatMatrix> fact;
           fact.set (R);
 
           if (down)
@@ -619,7 +612,7 @@
           FloatComplexColumnVector u =
             argu.float_complex_column_vector_value ();
 
-          FloatComplexCHOL fact;
+          chol<FloatComplexMatrix> fact;
           fact.set (R);
 
           if (down)
@@ -638,7 +631,7 @@
           Matrix R = argr.matrix_value ();
           ColumnVector u = argu.column_vector_value ();
 
-          CHOL fact;
+          chol<Matrix> fact;
           fact.set (R);
 
           if (down)
@@ -654,7 +647,7 @@
           ComplexMatrix R = argr.complex_matrix_value ();
           ComplexColumnVector u = argu.complex_column_vector_value ();
 
-          ComplexCHOL fact;
+          chol<ComplexMatrix> fact;
           fact.set (R);
 
           if (down)
@@ -793,7 +786,7 @@
           FloatMatrix R = argr.float_matrix_value ();
           FloatColumnVector u = argu.float_column_vector_value ();
 
-          FloatCHOL fact;
+          chol<FloatMatrix> fact;
           fact.set (R);
           err = fact.insert_sym (u, j-1);
 
@@ -806,7 +799,7 @@
           FloatComplexColumnVector u =
             argu.float_complex_column_vector_value ();
 
-          FloatComplexCHOL fact;
+          chol<FloatComplexMatrix> fact;
           fact.set (R);
           err = fact.insert_sym (u, j-1);
 
@@ -821,7 +814,7 @@
           Matrix R = argr.matrix_value ();
           ColumnVector u = argu.column_vector_value ();
 
-          CHOL fact;
+          chol<Matrix> fact;
           fact.set (R);
           err = fact.insert_sym (u, j-1);
 
@@ -834,7 +827,7 @@
           ComplexColumnVector u =
             argu.complex_column_vector_value ();
 
-          ComplexCHOL fact;
+          chol<ComplexMatrix> fact;
           fact.set (R);
           err = fact.insert_sym (u, j-1);
 
@@ -1028,7 +1021,7 @@
           // real case
           FloatMatrix R = argr.float_matrix_value ();
 
-          FloatCHOL fact;
+          chol<FloatMatrix> fact;
           fact.set (R);
           fact.delete_sym (j-1);
 
@@ -1039,7 +1032,7 @@
           // complex case
           FloatComplexMatrix R = argr.float_complex_matrix_value ();
 
-          FloatComplexCHOL fact;
+          chol<FloatComplexMatrix> fact;
           fact.set (R);
           fact.delete_sym (j-1);
 
@@ -1053,7 +1046,7 @@
           // real case
           Matrix R = argr.matrix_value ();
 
-          CHOL fact;
+          chol<Matrix> fact;
           fact.set (R);
           fact.delete_sym (j-1);
 
@@ -1064,7 +1057,7 @@
           // complex case
           ComplexMatrix R = argr.complex_matrix_value ();
 
-          ComplexCHOL fact;
+          chol<ComplexMatrix> fact;
           fact.set (R);
           fact.delete_sym (j-1);
 
@@ -1158,7 +1151,7 @@
           // real case
           FloatMatrix R = argr.float_matrix_value ();
 
-          FloatCHOL fact;
+          chol<FloatMatrix> fact;
           fact.set (R);
           fact.shift_sym (i-1, j-1);
 
@@ -1169,7 +1162,7 @@
           // complex case
           FloatComplexMatrix R = argr.float_complex_matrix_value ();
 
-          FloatComplexCHOL fact;
+          chol<FloatComplexMatrix> fact;
           fact.set (R);
           fact.shift_sym (i-1, j-1);
 
@@ -1183,7 +1176,7 @@
           // real case
           Matrix R = argr.matrix_value ();
 
-          CHOL fact;
+          chol<Matrix> fact;
           fact.set (R);
           fact.shift_sym (i-1, j-1);
 
@@ -1194,7 +1187,7 @@
           // complex case
           ComplexMatrix R = argr.complex_matrix_value ();
 
-          ComplexCHOL fact;
+          chol<ComplexMatrix> fact;
           fact.set (R);
           fact.shift_sym (i-1, j-1);
 
diff -r f08ae27289e4 -r 3c8a3d35661a liboctave/array/CMatrix.cc
--- a/liboctave/array/CMatrix.cc	Tue Feb 16 00:32:29 2016 -0500
+++ b/liboctave/array/CMatrix.cc	Tue Feb 16 02:47:29 2016 -0500
@@ -38,6 +38,7 @@
 #include "Array-util.h"
 #include "boolMatrix.h"
 #include "chMatrix.h"
+#include "chol.h"
 #include "dMatrix.h"
 #include "CMatrix.h"
 #include "CNDArray.h"
@@ -45,7 +46,6 @@
 #include "dRowVector.h"
 #include "CDiagMatrix.h"
 #include "dDiagMatrix.h"
-#include "CmplxCHOL.h"
 #include "schur.h"
 #include "CmplxSVD.h"
 #include "DET.h"
@@ -1109,7 +1109,7 @@
     {
       if (mattype.is_hermitian ())
         {
-          ComplexCHOL chol (*this, info, true, calc_cond);
+          chol<ComplexMatrix> chol (*this, info, true, calc_cond);
           if (info == 0)
             {
               if (calc_cond)
diff -r f08ae27289e4 -r 3c8a3d35661a liboctave/array/CMatrix.h
--- a/liboctave/array/CMatrix.h	Tue Feb 16 00:32:29 2016 -0500
+++ b/liboctave/array/CMatrix.h	Tue Feb 16 02:47:29 2016 -0500
@@ -50,6 +50,9 @@
   typedef Matrix real_matrix_type;
   typedef ComplexMatrix complex_matrix_type;
 
+  typedef double real_elt_type;
+  typedef Complex complex_elt_type;
+
   typedef void (*solve_singularity_handler) (double rcon);
 
   ComplexMatrix (void) : ComplexNDArray () { }
diff -r f08ae27289e4 -r 3c8a3d35661a liboctave/array/dMatrix.cc
--- a/liboctave/array/dMatrix.cc	Tue Feb 16 00:32:29 2016 -0500
+++ b/liboctave/array/dMatrix.cc	Tue Feb 16 02:47:29 2016 -0500
@@ -36,6 +36,7 @@
 #include "byte-swap.h"
 #include "boolMatrix.h"
 #include "chMatrix.h"
+#include "chol.h"
 #include "dMatrix.h"
 #include "dDiagMatrix.h"
 #include "CMatrix.h"
@@ -46,7 +47,6 @@
 #include "DET.h"
 #include "schur.h"
 #include "dbleSVD.h"
-#include "dbleCHOL.h"
 #include "f77-fcn.h"
 #include "functor.h"
 #include "lo-error.h"
@@ -799,7 +799,7 @@
     {
       if (mattype.is_hermitian ())
         {
-          CHOL chol (*this, info, true, calc_cond);
+          chol<Matrix> chol (*this, info, true, calc_cond);
           if (info == 0)
             {
               if (calc_cond)
diff -r f08ae27289e4 -r 3c8a3d35661a liboctave/array/dMatrix.h
--- a/liboctave/array/dMatrix.h	Tue Feb 16 00:32:29 2016 -0500
+++ b/liboctave/array/dMatrix.h	Tue Feb 16 02:47:29 2016 -0500
@@ -49,6 +49,9 @@
   typedef Matrix real_matrix_type;
   typedef ComplexMatrix complex_matrix_type;
 
+  typedef double real_elt_type;
+  typedef Complex complex_elt_type;
+
   typedef void (*solve_singularity_handler) (double rcon);
 
   Matrix (void) : NDArray () { }
diff -r f08ae27289e4 -r 3c8a3d35661a liboctave/array/fCMatrix.cc
--- a/liboctave/array/fCMatrix.cc	Tue Feb 16 00:32:29 2016 -0500
+++ b/liboctave/array/fCMatrix.cc	Tue Feb 16 02:47:29 2016 -0500
@@ -40,12 +40,12 @@
 #include "f77-fcn.h"
 #include "boolMatrix.h"
 #include "chMatrix.h"
+#include "chol.h"
 #include "fCMatrix.h"
 #include "fCNDArray.h"
 #include "fCDiagMatrix.h"
 #include "fCColVector.h"
 #include "fCRowVector.h"
-#include "fCmplxCHOL.h"
 #include "schur.h"
 #include "fCmplxSVD.h"
 #include "functor.h"
@@ -1112,7 +1112,7 @@
     {
       if (mattype.is_hermitian ())
         {
-          FloatComplexCHOL chol (*this, info, true, calc_cond);
+          chol<FloatComplexMatrix> chol (*this, info, true, calc_cond);
           if (info == 0)
             {
               if (calc_cond)
diff -r f08ae27289e4 -r 3c8a3d35661a liboctave/array/fCMatrix.h
--- a/liboctave/array/fCMatrix.h	Tue Feb 16 00:32:29 2016 -0500
+++ b/liboctave/array/fCMatrix.h	Tue Feb 16 02:47:29 2016 -0500
@@ -50,6 +50,9 @@
   typedef FloatMatrix real_matrix_type;
   typedef FloatComplexMatrix complex_matrix_type;
 
+  typedef float real_elt_type;
+  typedef FloatComplex complex_elt_type;
+
   typedef void (*solve_singularity_handler) (float rcon);
 
   FloatComplexMatrix (void) : FloatComplexNDArray () { }
diff -r f08ae27289e4 -r 3c8a3d35661a liboctave/array/fMatrix.cc
--- a/liboctave/array/fMatrix.cc	Tue Feb 16 00:32:29 2016 -0500
+++ b/liboctave/array/fMatrix.cc	Tue Feb 16 02:47:29 2016 -0500
@@ -36,6 +36,7 @@
 #include "Array-util.h"
 #include "boolMatrix.h"
 #include "chMatrix.h"
+#include "chol.h"
 #include "fMatrix.h"
 #include "fDiagMatrix.h"
 #include "fCMatrix.h"
@@ -47,7 +48,6 @@
 #include "byte-swap.h"
 #include "f77-fcn.h"
 #include "fMatrix.h"
-#include "floatCHOL.h"
 #include "schur.h"
 #include "floatSVD.h"
 #include "functor.h"
@@ -806,7 +806,7 @@
     {
       if (mattype.is_hermitian ())
         {
-          FloatCHOL chol (*this, info, true, calc_cond);
+          chol<FloatMatrix> chol (*this, info, true, calc_cond);
           if (info == 0)
             {
               if (calc_cond)
diff -r f08ae27289e4 -r 3c8a3d35661a liboctave/array/fMatrix.h
--- a/liboctave/array/fMatrix.h	Tue Feb 16 00:32:29 2016 -0500
+++ b/liboctave/array/fMatrix.h	Tue Feb 16 02:47:29 2016 -0500
@@ -49,6 +49,9 @@
   typedef FloatMatrix real_matrix_type;
   typedef FloatComplexMatrix complex_matrix_type;
 
+  typedef float real_elt_type;
+  typedef FloatComplex complex_elt_type;
+
   typedef void (*solve_singularity_handler) (float rcon);
 
   FloatMatrix (void) : FloatNDArray () { }
diff -r f08ae27289e4 -r 3c8a3d35661a liboctave/numeric/CmplxCHOL.cc
--- a/liboctave/numeric/CmplxCHOL.cc	Tue Feb 16 00:32:29 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,447 +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/>.
-
-*/
-
-#ifdef HAVE_CONFIG_H
-#  include <config.h>
-#endif
-
-#include <vector>
-
-#include "dMatrix.h"
-#include "dRowVector.h"
-#include "CmplxCHOL.h"
-#include "f77-fcn.h"
-#include "lo-error.h"
-#include "oct-locbuf.h"
-#include "oct-norm.h"
-#ifndef HAVE_QRUPDATE
-#  include "dbleQR.h"
-#endif
-
-extern "C"
-{
-  F77_RET_T
-  F77_FUNC (zpotrf, ZPOTRF) (F77_CONST_CHAR_ARG_DECL,
-                             const octave_idx_type&, Complex*,
-                             const octave_idx_type&, octave_idx_type&
-                             F77_CHAR_ARG_LEN_DECL);
-  F77_RET_T
-  F77_FUNC (zpotri, ZPOTRI) (F77_CONST_CHAR_ARG_DECL,
-                             const octave_idx_type&, Complex*,
-                             const octave_idx_type&, octave_idx_type&
-                             F77_CHAR_ARG_LEN_DECL);
-
-  F77_RET_T
-  F77_FUNC (zpocon, ZPOCON) (F77_CONST_CHAR_ARG_DECL,
-                             const octave_idx_type&, Complex*,
-                             const octave_idx_type&, const double&,
-                             double&, Complex*, double*, octave_idx_type&
-                             F77_CHAR_ARG_LEN_DECL);
-#ifdef HAVE_QRUPDATE
-
-  F77_RET_T
-  F77_FUNC (zch1up, ZCH1UP) (const octave_idx_type&, Complex*,
-                             const octave_idx_type&, Complex*, double*);
-
-  F77_RET_T
-  F77_FUNC (zch1dn, ZCH1DN) (const octave_idx_type&, Complex*,
-                             const octave_idx_type&, Complex*, double*,
-                             octave_idx_type&);
-
-  F77_RET_T
-  F77_FUNC (zchinx, ZCHINX) (const octave_idx_type&, Complex*,
-                             const octave_idx_type&, const octave_idx_type&,
-                             Complex*, double*, octave_idx_type&);
-
-  F77_RET_T
-  F77_FUNC (zchdex, ZCHDEX) (const octave_idx_type&, Complex*,
-                             const octave_idx_type&, const octave_idx_type&,
-                             double*);
-
-  F77_RET_T
-  F77_FUNC (zchshx, ZCHSHX) (const octave_idx_type&, Complex*,
-                             const octave_idx_type&, const octave_idx_type&,
-                             const octave_idx_type&, Complex*, double*);
-#endif
-}
-
-octave_idx_type
-ComplexCHOL::init (const ComplexMatrix& a, bool upper, bool calc_cond)
-{
-  octave_idx_type a_nr = a.rows ();
-  octave_idx_type a_nc = a.cols ();
-
-  if (a_nr != a_nc)
-    (*current_liboctave_error_handler) ("ComplexCHOL requires square matrix");
-
-  octave_idx_type n = a_nc;
-  octave_idx_type info;
-
-  is_upper = upper;
-
-  chol_mat.clear (n, n);
-  if (is_upper)
-    for (octave_idx_type j = 0; j < n; j++)
-      {
-        for (octave_idx_type i = 0; i <= j; i++)
-          chol_mat.xelem (i, j) = a(i, j);
-        for (octave_idx_type i = j+1; i < n; i++)
-          chol_mat.xelem (i, j) = 0.0;
-      }
-  else
-    for (octave_idx_type j = 0; j < n; j++)
-      {
-        for (octave_idx_type i = 0; i < j; i++)
-          chol_mat.xelem (i, j) = 0.0;
-        for (octave_idx_type i = j; i < n; i++)
-          chol_mat.xelem (i, j) = a(i, j);
-      }
-  Complex *h = chol_mat.fortran_vec ();
-
-  // Calculate the norm of the matrix, for later use.
-  double anorm = 0;
-  if (calc_cond)
-    anorm = xnorm (a, 1);
-
-  if (is_upper)
-    F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1), n, h, n, info
-                               F77_CHAR_ARG_LEN (1)));
-  else
-    F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), n, h, n, info
-                               F77_CHAR_ARG_LEN (1)));
-
-  xrcond = 0.0;
-  if (info > 0)
-    chol_mat.resize (info - 1, info - 1);
-  else if (calc_cond)
-    {
-      octave_idx_type zpocon_info = 0;
-
-      // Now calculate the condition number for non-singular matrix.
-      Array<Complex> z (dim_vector (2*n, 1));
-      Complex *pz = z.fortran_vec ();
-      Array<double> rz (dim_vector (n, 1));
-      double *prz = rz.fortran_vec ();
-      F77_XFCN (zpocon, ZPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h,
-                                 n, anorm, xrcond, pz, prz, zpocon_info
-                                 F77_CHAR_ARG_LEN (1)));
-
-      if (zpocon_info != 0)
-        info = -1;
-    }
-
-  return info;
-}
-
-static ComplexMatrix
-chol2inv_internal (const ComplexMatrix& r, bool is_upper = true)
-{
-  ComplexMatrix retval;
-
-  octave_idx_type r_nr = r.rows ();
-  octave_idx_type r_nc = r.cols ();
-
-  if (r_nr != r_nc)
-    (*current_liboctave_error_handler) ("chol2inv requires square matrix");
-
-  octave_idx_type n = r_nc;
-  octave_idx_type info;
-
-  ComplexMatrix tmp = r;
-
-  if (is_upper)
-    F77_XFCN (zpotri, ZPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n,
-                               tmp.fortran_vec (), n, info
-                               F77_CHAR_ARG_LEN (1)));
-  else
-    F77_XFCN (zpotri, ZPOTRI, (F77_CONST_CHAR_ARG2 ("L", 1), n,
-                               tmp.fortran_vec (), n, info
-                               F77_CHAR_ARG_LEN (1)));
-
-  // If someone thinks of a more graceful way of doing this (or
-  // faster for that matter :-)), please let me know!
-
-  if (n > 1)
-    {
-      if (is_upper)
-        for (octave_idx_type j = 0; j < r_nc; j++)
-          for (octave_idx_type i = j+1; i < r_nr; i++)
-            tmp.xelem (i, j) = tmp.xelem (j, i);
-      else
-        for (octave_idx_type j = 0; j < r_nc; j++)
-          for (octave_idx_type i = j+1; i < r_nr; i++)
-            tmp.xelem (j, i) = tmp.xelem (i, j);
-    }
-
-  retval = tmp;
-
-  return retval;
-}
-
-// Compute the inverse of a matrix using the Cholesky factorization.
-ComplexMatrix
-ComplexCHOL::inverse (void) const
-{
-  return chol2inv_internal (chol_mat, is_upper);
-}
-
-void
-ComplexCHOL::set (const ComplexMatrix& R)
-{
-  if (R.is_square ())
-    chol_mat = R;
-  else
-    (*current_liboctave_error_handler) ("CHOL requires square matrix");
-}
-
-#ifdef HAVE_QRUPDATE
-
-void
-ComplexCHOL::update (const ComplexColumnVector& u)
-{
-  octave_idx_type n = chol_mat.rows ();
-
-  if (u.numel () != n)
-    (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
-
-  ComplexColumnVector utmp = u;
-
-  OCTAVE_LOCAL_BUFFER (double, rw, n);
-
-  F77_XFCN (zch1up, ZCH1UP, (n, chol_mat.fortran_vec (), chol_mat.rows (),
-                             utmp.fortran_vec (), rw));
-}
-
-octave_idx_type
-ComplexCHOL::downdate (const ComplexColumnVector& u)
-{
-  octave_idx_type info = -1;
-
-  octave_idx_type n = chol_mat.rows ();
-
-  if (u.numel () != n)
-    (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
-
-  ComplexColumnVector utmp = u;
-
-  OCTAVE_LOCAL_BUFFER (double, rw, n);
-
-  F77_XFCN (zch1dn, ZCH1DN, (n, chol_mat.fortran_vec (), chol_mat.rows (),
-                             utmp.fortran_vec (), rw, info));
-
-  return info;
-}
-
-octave_idx_type
-ComplexCHOL::insert_sym (const ComplexColumnVector& u, octave_idx_type j)
-{
-  octave_idx_type info = -1;
-
-  octave_idx_type n = chol_mat.rows ();
-
-  if (u.numel () != n + 1)
-    (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
-  if (j < 0 || j > n)
-    (*current_liboctave_error_handler) ("cholinsert: index out of range");
-
-  ComplexColumnVector utmp = u;
-
-  OCTAVE_LOCAL_BUFFER (double, rw, n);
-
-  chol_mat.resize (n+1, n+1);
-
-  F77_XFCN (zchinx, ZCHINX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
-                             j + 1, utmp.fortran_vec (), rw, info));
-
-  return info;
-}
-
-void
-ComplexCHOL::delete_sym (octave_idx_type j)
-{
-  octave_idx_type n = chol_mat.rows ();
-
-  if (j < 0 || j > n-1)
-    (*current_liboctave_error_handler) ("choldelete: index out of range");
-
-  OCTAVE_LOCAL_BUFFER (double, rw, n);
-
-  F77_XFCN (zchdex, ZCHDEX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
-                             j + 1, rw));
-
-  chol_mat.resize (n-1, n-1);
-}
-
-void
-ComplexCHOL::shift_sym (octave_idx_type i, octave_idx_type j)
-{
-  octave_idx_type n = chol_mat.rows ();
-
-  if (i < 0 || i > n-1 || j < 0 || j > n-1)
-    (*current_liboctave_error_handler) ("cholshift: index out of range");
-
-  OCTAVE_LOCAL_BUFFER (Complex, w, n);
-  OCTAVE_LOCAL_BUFFER (double, rw, n);
-
-  F77_XFCN (zchshx, ZCHSHX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
-                             i + 1, j + 1, w, rw));
-}
-
-#else
-
-void
-ComplexCHOL::update (const ComplexColumnVector& u)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type n = chol_mat.rows ();
-
-  if (u.numel () != n)
-    (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
-
-  init (chol_mat.hermitian () * chol_mat
-        + ComplexMatrix (u) * ComplexMatrix (u).hermitian (), true, false);
-}
-
-static bool
-singular (const ComplexMatrix& a)
-{
-  for (octave_idx_type i = 0; i < a.rows (); i++)
-    if (a(i,i) == 0.0) return true;
-  return false;
-}
-
-octave_idx_type
-ComplexCHOL::downdate (const ComplexColumnVector& u)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type info = -1;
-
-  octave_idx_type n = chol_mat.rows ();
-
-  if (u.numel () != n)
-    (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
-
-  if (singular (chol_mat))
-    info = 2;
-  else
-    {
-      info = init (chol_mat.hermitian () * chol_mat
-                   - ComplexMatrix (u) * ComplexMatrix (u).hermitian (),
-                   true, false);
-      if (info) info = 1;
-    }
-
-  return info;
-}
-
-octave_idx_type
-ComplexCHOL::insert_sym (const ComplexColumnVector& u, octave_idx_type j)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type info = -1;
-
-  octave_idx_type n = chol_mat.rows ();
-
-  if (u.numel () != n + 1)
-    (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
-  if (j < 0 || j > n)
-    (*current_liboctave_error_handler) ("cholinsert: index out of range");
-
-  if (singular (chol_mat))
-    info = 2;
-  else if (u(j).imag () != 0.0)
-    info = 3;
-  else
-    {
-      ComplexMatrix a = chol_mat.hermitian () * chol_mat;
-      ComplexMatrix a1 (n+1, n+1);
-      for (octave_idx_type k = 0; k < n+1; k++)
-        for (octave_idx_type l = 0; l < n+1; l++)
-          {
-            if (l == j)
-              a1(k, l) = u(k);
-            else if (k == j)
-              a1(k, l) = std::conj (u(l));
-            else
-              a1(k, l) = a(k < j ? k : k-1, l < j ? l : l-1);
-          }
-      info = init (a1, true, false);
-      if (info) info = 1;
-    }
-
-  return info;
-}
-
-void
-ComplexCHOL::delete_sym (octave_idx_type j)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type n = chol_mat.rows ();
-
-  if (j < 0 || j > n-1)
-    (*current_liboctave_error_handler) ("choldelete: index out of range");
-
-  ComplexMatrix a = chol_mat.hermitian () * chol_mat;
-  a.delete_elements (1, idx_vector (j));
-  a.delete_elements (0, idx_vector (j));
-  init (a, true, false);
-}
-
-void
-ComplexCHOL::shift_sym (octave_idx_type i, octave_idx_type j)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type n = chol_mat.rows ();
-
-  if (i < 0 || i > n-1 || j < 0 || j > n-1)
-    (*current_liboctave_error_handler) ("cholshift: index out of range");
-
-  ComplexMatrix a = chol_mat.hermitian () * chol_mat;
-  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;
-    }
-
-  init (a.index (idx_vector (p), idx_vector (p)), true, false);
-}
-
-#endif
-
-ComplexMatrix
-chol2inv (const ComplexMatrix& r)
-{
-  return chol2inv_internal (r);
-}
diff -r f08ae27289e4 -r 3c8a3d35661a liboctave/numeric/CmplxCHOL.h
--- a/liboctave/numeric/CmplxCHOL.h	Tue Feb 16 00:32:29 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,103 +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_CmplxCHOL_h)
-#define octave_CmplxCHOL_h 1
-
-#include "octave-config.h"
-
-#include <iosfwd>
-
-#include "CMatrix.h"
-#include "CColVector.h"
-
-class
-OCTAVE_API
-ComplexCHOL
-{
-public:
-
-  ComplexCHOL (void) : chol_mat (), xrcond (0) { }
-
-  ComplexCHOL (const ComplexMatrix& a, bool upper = true, bool calc_cond = false)
-    : chol_mat (), xrcond (0)
-  {
-    init (a, upper, calc_cond);
-  }
-
-  ComplexCHOL (const ComplexMatrix& a, octave_idx_type& info, bool upper = true,
-               bool calc_cond = false)
-    : chol_mat (), xrcond (0)
-  {
-    info = init (a, upper, calc_cond);
-  }
-
-  ComplexCHOL (const ComplexCHOL& a)
-    : chol_mat (a.chol_mat), xrcond (a.xrcond) { }
-
-  ComplexCHOL& operator = (const ComplexCHOL& a)
-  {
-    if (this != &a)
-      {
-        chol_mat = a.chol_mat;
-        xrcond = a.xrcond;
-      }
-
-    return *this;
-  }
-
-  ComplexMatrix chol_matrix (void) const { return chol_mat; }
-
-  double rcond (void) const { return xrcond; }
-
-  ComplexMatrix inverse (void) const;
-
-  void set (const ComplexMatrix& R);
-
-  void update (const ComplexColumnVector& u);
-
-  octave_idx_type downdate (const ComplexColumnVector& u);
-
-  octave_idx_type insert_sym (const ComplexColumnVector& u, octave_idx_type j);
-
-  void delete_sym (octave_idx_type j);
-
-  void shift_sym (octave_idx_type i, octave_idx_type j);
-
-  friend OCTAVE_API std::ostream& operator << (std::ostream& os,
-                                               const ComplexCHOL& a);
-
-private:
-
-  ComplexMatrix chol_mat;
-
-  double xrcond;
-
-  bool is_upper;
-
-  octave_idx_type init (const ComplexMatrix& a, bool upper, bool calc_cond);
-};
-
-ComplexMatrix OCTAVE_API chol2inv (const ComplexMatrix& r);
-
-#endif
diff -r f08ae27289e4 -r 3c8a3d35661a liboctave/numeric/chol.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/numeric/chol.cc	Tue Feb 16 02:47:29 2016 -0500
@@ -0,0 +1,1298 @@
+/*
+
+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/>.
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#  include <config.h>
+#endif
+
+#include <vector>
+
+
+#include "CColVector.h"
+#include "CMatrix.h"
+#include "chol.h"
+#include "dColVector.h"
+#include "dMatrix.h"
+#include "f77-fcn.h"
+#include "fCColVector.h"
+#include "fCMatrix.h"
+#include "fColVector.h"
+#include "fMatrix.h"
+#include "lo-error.h"
+#include "oct-locbuf.h"
+#include "oct-norm.h"
+
+#if ! defined (HAVE_QRUPDATE)
+#  include "CmplxQR.h"
+#  include "dbleQR.h"
+#  include "fCmplxQR.h"
+#  include "floatQR.h"
+#endif
+
+extern "C"
+{
+  F77_RET_T
+  F77_FUNC (dpotrf, DPOTRF) (F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type&, double*,
+                             const octave_idx_type&, octave_idx_type&
+                             F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (dpotri, DPOTRI) (F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type&, double*,
+                             const octave_idx_type&, octave_idx_type&
+                             F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (dpocon, DPOCON) (F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type&, double*,
+                             const octave_idx_type&, const double&,
+                             double&, double*, octave_idx_type*,
+                             octave_idx_type&
+                             F77_CHAR_ARG_LEN_DECL);
+#ifdef HAVE_QRUPDATE
+
+  F77_RET_T
+  F77_FUNC (dch1up, DCH1UP) (const octave_idx_type&, double*,
+                             const octave_idx_type&, double*, double*);
+
+  F77_RET_T
+  F77_FUNC (dch1dn, DCH1DN) (const octave_idx_type&, double*,
+                             const octave_idx_type&, double*, double*,
+                             octave_idx_type&);
+
+  F77_RET_T
+  F77_FUNC (dchinx, DCHINX) (const octave_idx_type&, double*,
+                             const octave_idx_type&, const octave_idx_type&,
+                             double*, double*, octave_idx_type&);
+
+  F77_RET_T
+  F77_FUNC (dchdex, DCHDEX) (const octave_idx_type&, double*,
+                             const octave_idx_type&, const octave_idx_type&,
+                             double*);
+
+  F77_RET_T
+  F77_FUNC (dchshx, DCHSHX) (const octave_idx_type&, double*,
+                             const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, double*);
+#endif
+
+  F77_RET_T
+  F77_FUNC (spotrf, SPOTRF) (F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type&, float*,
+                             const octave_idx_type&, octave_idx_type&
+                             F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (spotri, SPOTRI) (F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type&, float*,
+                             const octave_idx_type&, octave_idx_type&
+                             F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (spocon, SPOCON) (F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type&, float*,
+                             const octave_idx_type&, const float&,
+                             float&, float*, octave_idx_type*,
+                             octave_idx_type&
+                             F77_CHAR_ARG_LEN_DECL);
+#ifdef HAVE_QRUPDATE
+
+  F77_RET_T
+  F77_FUNC (sch1up, SCH1UP) (const octave_idx_type&, float*,
+                             const octave_idx_type&, float*, float*);
+
+  F77_RET_T
+  F77_FUNC (sch1dn, SCH1DN) (const octave_idx_type&, float*,
+                             const octave_idx_type&, float*, float*,
+                             octave_idx_type&);
+
+  F77_RET_T
+  F77_FUNC (schinx, SCHINX) (const octave_idx_type&, float*,
+                             const octave_idx_type&, const octave_idx_type&,
+                             float*, float*, octave_idx_type&);
+
+  F77_RET_T
+  F77_FUNC (schdex, SCHDEX) (const octave_idx_type&, float*,
+                             const octave_idx_type&, const octave_idx_type&,
+                             float*);
+
+  F77_RET_T
+  F77_FUNC (schshx, SCHSHX) (const octave_idx_type&, float*,
+                             const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, float*);
+#endif
+
+  F77_RET_T
+  F77_FUNC (zpotrf, ZPOTRF) (F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type&, Complex*,
+                             const octave_idx_type&, octave_idx_type&
+                             F77_CHAR_ARG_LEN_DECL);
+  F77_RET_T
+  F77_FUNC (zpotri, ZPOTRI) (F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type&, Complex*,
+                             const octave_idx_type&, octave_idx_type&
+                             F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (zpocon, ZPOCON) (F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type&, Complex*,
+                             const octave_idx_type&, const double&,
+                             double&, Complex*, double*, octave_idx_type&
+                             F77_CHAR_ARG_LEN_DECL);
+#ifdef HAVE_QRUPDATE
+
+  F77_RET_T
+  F77_FUNC (zch1up, ZCH1UP) (const octave_idx_type&, Complex*,
+                             const octave_idx_type&, Complex*, double*);
+
+  F77_RET_T
+  F77_FUNC (zch1dn, ZCH1DN) (const octave_idx_type&, Complex*,
+                             const octave_idx_type&, Complex*, double*,
+                             octave_idx_type&);
+
+  F77_RET_T
+  F77_FUNC (zchinx, ZCHINX) (const octave_idx_type&, Complex*,
+                             const octave_idx_type&, const octave_idx_type&,
+                             Complex*, double*, octave_idx_type&);
+
+  F77_RET_T
+  F77_FUNC (zchdex, ZCHDEX) (const octave_idx_type&, Complex*,
+                             const octave_idx_type&, const octave_idx_type&,
+                             double*);
+
+  F77_RET_T
+  F77_FUNC (zchshx, ZCHSHX) (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 (cpotrf, CPOTRF) (F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type&, FloatComplex*,
+                             const octave_idx_type&, octave_idx_type&
+                             F77_CHAR_ARG_LEN_DECL);
+  F77_RET_T
+  F77_FUNC (cpotri, CPOTRI) (F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type&, FloatComplex*,
+                             const octave_idx_type&, octave_idx_type&
+                             F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (cpocon, CPOCON) (F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type&, FloatComplex*,
+                             const octave_idx_type&, const float&,
+                             float&, FloatComplex*, float*, octave_idx_type&
+                             F77_CHAR_ARG_LEN_DECL);
+#ifdef HAVE_QRUPDATE
+
+  F77_RET_T
+  F77_FUNC (cch1up, CCH1UP) (const octave_idx_type&, FloatComplex*,
+                             const octave_idx_type&, FloatComplex*, float*);
+
+  F77_RET_T
+  F77_FUNC (cch1dn, CCH1DN) (const octave_idx_type&, FloatComplex*,
+                             const octave_idx_type&, FloatComplex*,
+                             float*, octave_idx_type&);
+
+  F77_RET_T
+  F77_FUNC (cchinx, CCHINX) (const octave_idx_type&, FloatComplex*,
+                             const octave_idx_type&, const octave_idx_type&,
+                             FloatComplex*, float*, octave_idx_type&);
+
+  F77_RET_T
+  F77_FUNC (cchdex, CCHDEX) (const octave_idx_type&, FloatComplex*,
+                             const octave_idx_type&, const octave_idx_type&,
+                             float*);
+
+  F77_RET_T
+  F77_FUNC (cchshx, CCHSHX) (const octave_idx_type&, FloatComplex*,
+                             const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, FloatComplex*, float*);
+#endif
+}
+
+static Matrix
+chol2inv_internal (const Matrix& r, bool is_upper = true)
+{
+  Matrix retval;
+
+  octave_idx_type r_nr = r.rows ();
+  octave_idx_type r_nc = r.cols ();
+
+  if (r_nr != r_nc)
+    (*current_liboctave_error_handler) ("chol2inv requires square matrix");
+
+  octave_idx_type n = r_nc;
+  octave_idx_type info = 0;
+
+  Matrix tmp = r;
+  double *v = tmp.fortran_vec ();
+
+  if (info == 0)
+    {
+      if (is_upper)
+        F77_XFCN (dpotri, DPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n,
+                                   v, n, info
+                                   F77_CHAR_ARG_LEN (1)));
+      else
+        F77_XFCN (dpotri, DPOTRI, (F77_CONST_CHAR_ARG2 ("L", 1), n,
+                                   v, n, info
+                                   F77_CHAR_ARG_LEN (1)));
+
+      // If someone thinks of a more graceful way of doing this (or
+      // faster for that matter :-)), please let me know!
+
+      if (n > 1)
+        {
+          if (is_upper)
+            for (octave_idx_type j = 0; j < r_nc; j++)
+              for (octave_idx_type i = j+1; i < r_nr; i++)
+                tmp.xelem (i, j) = tmp.xelem (j, i);
+          else
+            for (octave_idx_type j = 0; j < r_nc; j++)
+              for (octave_idx_type i = j+1; i < r_nr; i++)
+                tmp.xelem (j, i) = tmp.xelem (i, j);
+        }
+
+      retval = tmp;
+    }
+
+  return retval;
+}
+
+static FloatMatrix
+chol2inv_internal (const FloatMatrix& r, bool is_upper = true)
+{
+  FloatMatrix retval;
+
+  octave_idx_type r_nr = r.rows ();
+  octave_idx_type r_nc = r.cols ();
+
+  if (r_nr != r_nc)
+    (*current_liboctave_error_handler) ("chol2inv requires square matrix");
+
+  octave_idx_type n = r_nc;
+  octave_idx_type info = 0;
+
+  FloatMatrix tmp = r;
+  float *v = tmp.fortran_vec ();
+
+  if (info == 0)
+    {
+      if (is_upper)
+        F77_XFCN (spotri, SPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n,
+                                   v, n, info
+                                   F77_CHAR_ARG_LEN (1)));
+      else
+        F77_XFCN (spotri, SPOTRI, (F77_CONST_CHAR_ARG2 ("L", 1), n,
+                                   v, n, info
+                                   F77_CHAR_ARG_LEN (1)));
+
+      // If someone thinks of a more graceful way of doing this (or
+      // faster for that matter :-)), please let me know!
+
+      if (n > 1)
+        {
+          if (is_upper)
+            for (octave_idx_type j = 0; j < r_nc; j++)
+              for (octave_idx_type i = j+1; i < r_nr; i++)
+                tmp.xelem (i, j) = tmp.xelem (j, i);
+          else
+            for (octave_idx_type j = 0; j < r_nc; j++)
+              for (octave_idx_type i = j+1; i < r_nr; i++)
+                tmp.xelem (j, i) = tmp.xelem (i, j);
+        }
+
+      retval = tmp;
+    }
+
+  return retval;
+}
+
+static ComplexMatrix
+chol2inv_internal (const ComplexMatrix& r, bool is_upper = true)
+{
+  ComplexMatrix retval;
+
+  octave_idx_type r_nr = r.rows ();
+  octave_idx_type r_nc = r.cols ();
+
+  if (r_nr != r_nc)
+    (*current_liboctave_error_handler) ("chol2inv requires square matrix");
+
+  octave_idx_type n = r_nc;
+  octave_idx_type info;
+
+  ComplexMatrix tmp = r;
+
+  if (is_upper)
+    F77_XFCN (zpotri, ZPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n,
+                               tmp.fortran_vec (), n, info
+                               F77_CHAR_ARG_LEN (1)));
+  else
+    F77_XFCN (zpotri, ZPOTRI, (F77_CONST_CHAR_ARG2 ("L", 1), n,
+                               tmp.fortran_vec (), n, info
+                               F77_CHAR_ARG_LEN (1)));
+
+  // If someone thinks of a more graceful way of doing this (or
+  // faster for that matter :-)), please let me know!
+
+  if (n > 1)
+    {
+      if (is_upper)
+        for (octave_idx_type j = 0; j < r_nc; j++)
+          for (octave_idx_type i = j+1; i < r_nr; i++)
+            tmp.xelem (i, j) = tmp.xelem (j, i);
+      else
+        for (octave_idx_type j = 0; j < r_nc; j++)
+          for (octave_idx_type i = j+1; i < r_nr; i++)
+            tmp.xelem (j, i) = tmp.xelem (i, j);
+    }
+
+  retval = tmp;
+
+  return retval;
+}
+
+static FloatComplexMatrix
+chol2inv_internal (const FloatComplexMatrix& r, bool is_upper = true)
+{
+  FloatComplexMatrix retval;
+
+  octave_idx_type r_nr = r.rows ();
+  octave_idx_type r_nc = r.cols ();
+
+  if (r_nr != r_nc)
+    (*current_liboctave_error_handler) ("chol2inv requires square matrix");
+
+  octave_idx_type n = r_nc;
+  octave_idx_type info;
+
+  FloatComplexMatrix tmp = r;
+
+  if (is_upper)
+    F77_XFCN (cpotri, CPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n,
+                               tmp.fortran_vec (), n, info
+                               F77_CHAR_ARG_LEN (1)));
+  else
+    F77_XFCN (cpotri, CPOTRI, (F77_CONST_CHAR_ARG2 ("L", 1), n,
+                               tmp.fortran_vec (), n, info
+                               F77_CHAR_ARG_LEN (1)));
+
+  // If someone thinks of a more graceful way of doing this (or
+  // faster for that matter :-)), please let me know!
+
+  if (n > 1)
+    {
+      if (is_upper)
+        for (octave_idx_type j = 0; j < r_nc; j++)
+          for (octave_idx_type i = j+1; i < r_nr; i++)
+            tmp.xelem (i, j) = tmp.xelem (j, i);
+      else
+        for (octave_idx_type j = 0; j < r_nc; j++)
+          for (octave_idx_type i = j+1; i < r_nr; i++)
+            tmp.xelem (j, i) = tmp.xelem (i, j);
+    }
+
+  retval = tmp;
+
+  return retval;
+}
+
+template <typename T>
+T
+chol2inv (const T& r)
+{
+  return chol2inv_internal (r);
+}
+
+// Compute the inverse of a matrix using the Cholesky factorization.
+template <typename T>
+T
+chol<T>::inverse (void) const
+{
+  return chol2inv_internal (chol_mat, is_upper);
+}
+
+template <typename T>
+void
+chol<T>::set (const T& R)
+{
+  if (! R.is_square ())
+    (*current_liboctave_error_handler) ("chol: requires square matrix");
+
+  chol_mat = R;
+}
+
+#if ! defined (HAVE_QRUPDATE)
+
+template <typename T>
+void
+chol<T>::update (const T::VT& u)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type n = chol_mat.rows ();
+
+  if (u.numel () != n)
+    (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
+
+  init (chol_mat.hermitian () * chol_mat + T (u) * T (u).hermitian (),
+        true, false);
+}
+
+template <typename T>
+static bool
+singular (const T& a)
+{
+  static typename T::element_type zero ();
+  for (octave_idx_type i = 0; i < a.rows (); i++)
+    if (a(i,i) == zero) return true;
+  return false;
+}
+
+template <typename T>
+octave_idx_type
+chol<T>::downdate (const T::VT& u)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type info = -1;
+
+  octave_idx_type n = chol_mat.rows ();
+
+  if (u.numel () != n)
+    (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
+
+  if (singular (chol_mat))
+    info = 2;
+  else
+    {
+      info = init (chol_mat.hermitian () * chol_mat
+                   - T (u) * T (u).hermitian (), true, false);
+      if (info) info = 1;
+    }
+
+  return info;
+}
+
+template <typename T>
+octave_idx_type
+chol<T>::insert_sym (const T::VT& u, octave_idx_type j)
+{
+  static typename T::element_type zero ();
+
+  warn_qrupdate_once ();
+
+  octave_idx_type info = -1;
+
+  octave_idx_type n = chol_mat.rows ();
+
+  if (u.numel () != n + 1)
+    (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
+  if (j < 0 || j > n)
+    (*current_liboctave_error_handler) ("cholinsert: index out of range");
+
+  if (singular (chol_mat))
+    info = 2;
+  else if (ximag (u(j)) != zero)
+    info = 3;
+  else
+    {
+      T a = chol_mat.hermitian () * chol_mat;
+      T a1 (n+1, n+1);
+      for (octave_idx_type k = 0; k < n+1; k++)
+        for (octave_idx_type l = 0; l < n+1; l++)
+          {
+            if (l == j)
+              a1(k, l) = u(k);
+            else if (k == j)
+              a1(k, l) = xconj (u(l));
+            else
+              a1(k, l) = a(k < j ? k : k-1, l < j ? l : l-1);
+          }
+      info = init (a1, true, false);
+      if (info) info = 1;
+    }
+
+  return info;
+}
+
+template <typename T>
+void
+chol<T>::delete_sym (octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type n = chol_mat.rows ();
+
+  if (j < 0 || j > n-1)
+    (*current_liboctave_error_handler) ("choldelete: index out of range");
+
+  T a = chol_mat.hermitian () * chol_mat;
+  a.delete_elements (1, idx_vector (j));
+  a.delete_elements (0, idx_vector (j));
+  init (a, true, false);
+}
+
+template <typename T>
+void
+chol<T>::shift_sym (octave_idx_type i, octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type n = chol_mat.rows ();
+
+  if (i < 0 || i > n-1 || j < 0 || j > n-1)
+    (*current_liboctave_error_handler) ("cholshift: index out of range");
+
+  T a = chol_mat.hermitian () * chol_mat;
+  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;
+    }
+
+  init (a.index (idx_vector (p), idx_vector (p)), true, false);
+}
+
+#endif
+
+// Specializations.
+
+template <>
+octave_idx_type
+chol<Matrix>::init (const Matrix& a, bool upper, bool calc_cond)
+{
+  octave_idx_type a_nr = a.rows ();
+  octave_idx_type a_nc = a.cols ();
+
+  if (a_nr != a_nc)
+    (*current_liboctave_error_handler) ("chol: requires square matrix");
+
+  octave_idx_type n = a_nc;
+  octave_idx_type info;
+
+  is_upper = upper;
+
+  chol_mat.clear (n, n);
+  if (is_upper)
+    for (octave_idx_type j = 0; j < n; j++)
+      {
+        for (octave_idx_type i = 0; i <= j; i++)
+          chol_mat.xelem (i, j) = a(i, j);
+        for (octave_idx_type i = j+1; i < n; i++)
+          chol_mat.xelem (i, j) = 0.0;
+      }
+  else
+    for (octave_idx_type j = 0; j < n; j++)
+      {
+        for (octave_idx_type i = 0; i < j; i++)
+          chol_mat.xelem (i, j) = 0.0;
+        for (octave_idx_type i = j; i < n; i++)
+          chol_mat.xelem (i, j) = a(i, j);
+      }
+  double *h = chol_mat.fortran_vec ();
+
+  // Calculate the norm of the matrix, for later use.
+  double anorm = 0;
+  if (calc_cond)
+    anorm = xnorm (a, 1);
+
+  if (is_upper)
+    F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1), n, h, n, info
+                               F77_CHAR_ARG_LEN (1)));
+  else
+    F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), n, h, n, info
+                               F77_CHAR_ARG_LEN (1)));
+
+  xrcond = 0.0;
+  if (info > 0)
+    chol_mat.resize (info - 1, info - 1);
+  else if (calc_cond)
+    {
+      octave_idx_type dpocon_info = 0;
+
+      // Now calculate the condition number for non-singular matrix.
+      Array<double> z (dim_vector (3*n, 1));
+      double *pz = z.fortran_vec ();
+      Array<octave_idx_type> iz (dim_vector (n, 1));
+      octave_idx_type *piz = iz.fortran_vec ();
+      if (is_upper)
+        F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h,
+                                   n, anorm, xrcond, pz, piz, dpocon_info
+                                   F77_CHAR_ARG_LEN (1)));
+      else
+        F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 ("L", 1), n, h,
+                                   n, anorm, xrcond, pz, piz, dpocon_info
+                                   F77_CHAR_ARG_LEN (1)));
+
+      if (dpocon_info != 0)
+        info = -1;
+    }
+
+  return info;
+}
+
+#if defined (HAVE_QRUPDATE)
+
+template <>
+void
+chol<Matrix>::update (const ColumnVector& u)
+{
+  octave_idx_type n = chol_mat.rows ();
+
+  if (u.numel () != n)
+    (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
+
+  ColumnVector utmp = u;
+
+  OCTAVE_LOCAL_BUFFER (double, w, n);
+
+  F77_XFCN (dch1up, DCH1UP, (n, chol_mat.fortran_vec (), chol_mat.rows (),
+                             utmp.fortran_vec (), w));
+}
+
+template <>
+octave_idx_type
+chol<Matrix>::downdate (const ColumnVector& u)
+{
+  octave_idx_type info = -1;
+
+  octave_idx_type n = chol_mat.rows ();
+
+  if (u.numel () != n)
+    (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
+
+  ColumnVector utmp = u;
+
+  OCTAVE_LOCAL_BUFFER (double, w, n);
+
+  F77_XFCN (dch1dn, DCH1DN, (n, chol_mat.fortran_vec (), chol_mat.rows (),
+                             utmp.fortran_vec (), w, info));
+
+  return info;
+}
+
+template <>
+octave_idx_type
+chol<Matrix>::insert_sym (const ColumnVector& u, octave_idx_type j)
+{
+  octave_idx_type info = -1;
+
+  octave_idx_type n = chol_mat.rows ();
+
+  if (u.numel () != n + 1)
+    (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
+  if (j < 0 || j > n)
+    (*current_liboctave_error_handler) ("cholinsert: index out of range");
+
+  ColumnVector utmp = u;
+
+  OCTAVE_LOCAL_BUFFER (double, w, n);
+
+  chol_mat.resize (n+1, n+1);
+
+  F77_XFCN (dchinx, DCHINX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
+                             j + 1, utmp.fortran_vec (), w, info));
+
+  return info;
+}
+
+template <>
+void
+chol<Matrix>::delete_sym (octave_idx_type j)
+{
+  octave_idx_type n = chol_mat.rows ();
+
+  if (j < 0 || j > n-1)
+    (*current_liboctave_error_handler) ("choldelete: index out of range");
+
+  OCTAVE_LOCAL_BUFFER (double, w, n);
+
+  F77_XFCN (dchdex, DCHDEX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
+                             j + 1, w));
+
+  chol_mat.resize (n-1, n-1);
+}
+
+template <>
+void
+chol<Matrix>::shift_sym (octave_idx_type i, octave_idx_type j)
+{
+  octave_idx_type n = chol_mat.rows ();
+
+  if (i < 0 || i > n-1 || j < 0 || j > n-1)
+    (*current_liboctave_error_handler) ("cholshift: index out of range");
+
+  OCTAVE_LOCAL_BUFFER (double, w, 2*n);
+
+  F77_XFCN (dchshx, DCHSHX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
+                             i + 1, j + 1, w));
+}
+
+#endif
+
+template <>
+octave_idx_type
+chol<FloatMatrix>::init (const FloatMatrix& a, bool upper, bool calc_cond)
+{
+  octave_idx_type a_nr = a.rows ();
+  octave_idx_type a_nc = a.cols ();
+
+  if (a_nr != a_nc)
+    (*current_liboctave_error_handler) ("chol: requires square matrix");
+
+  octave_idx_type n = a_nc;
+  octave_idx_type info;
+
+  is_upper = upper;
+
+  chol_mat.clear (n, n);
+  if (is_upper)
+    for (octave_idx_type j = 0; j < n; j++)
+      {
+        for (octave_idx_type i = 0; i <= j; i++)
+          chol_mat.xelem (i, j) = a(i, j);
+        for (octave_idx_type i = j+1; i < n; i++)
+          chol_mat.xelem (i, j) = 0.0f;
+      }
+  else
+    for (octave_idx_type j = 0; j < n; j++)
+      {
+        for (octave_idx_type i = 0; i < j; i++)
+          chol_mat.xelem (i, j) = 0.0f;
+        for (octave_idx_type i = j; i < n; i++)
+          chol_mat.xelem (i, j) = a(i, j);
+      }
+  float *h = chol_mat.fortran_vec ();
+
+  // Calculate the norm of the matrix, for later use.
+  float anorm = 0;
+  if (calc_cond)
+    anorm = xnorm (a, 1);
+
+  if (is_upper)
+    F77_XFCN (spotrf, SPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1), n, h, n, info
+                               F77_CHAR_ARG_LEN (1)));
+  else
+    F77_XFCN (spotrf, SPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), n, h, n, info
+                               F77_CHAR_ARG_LEN (1)));
+
+  xrcond = 0.0;
+  if (info > 0)
+    chol_mat.resize (info - 1, info - 1);
+  else if (calc_cond)
+    {
+      octave_idx_type spocon_info = 0;
+
+      // Now calculate the condition number for non-singular matrix.
+      Array<float> z (dim_vector (3*n, 1));
+      float *pz = z.fortran_vec ();
+      Array<octave_idx_type> iz (dim_vector (n, 1));
+      octave_idx_type *piz = iz.fortran_vec ();
+      if (is_upper)
+        F77_XFCN (spocon, SPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h,
+                                   n, anorm, xrcond, pz, piz, spocon_info
+                                   F77_CHAR_ARG_LEN (1)));
+      else
+        F77_XFCN (spocon, SPOCON, (F77_CONST_CHAR_ARG2 ("L", 1), n, h,
+                                   n, anorm, xrcond, pz, piz, spocon_info
+                                   F77_CHAR_ARG_LEN (1)));
+
+      if (spocon_info != 0)
+        info = -1;
+    }
+
+  return info;
+}
+
+#ifdef HAVE_QRUPDATE
+
+template <>
+void
+chol<FloatMatrix>::update (const FloatColumnVector& u)
+{
+  octave_idx_type n = chol_mat.rows ();
+
+  if (u.numel () != n)
+    (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
+
+  FloatColumnVector utmp = u;
+
+  OCTAVE_LOCAL_BUFFER (float, w, n);
+
+  F77_XFCN (sch1up, SCH1UP, (n, chol_mat.fortran_vec (), chol_mat.rows (),
+                             utmp.fortran_vec (), w));
+}
+
+template <>
+octave_idx_type
+chol<FloatMatrix>::downdate (const FloatColumnVector& u)
+{
+  octave_idx_type info = -1;
+
+  octave_idx_type n = chol_mat.rows ();
+
+  if (u.numel () != n)
+    (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
+
+  FloatColumnVector utmp = u;
+
+  OCTAVE_LOCAL_BUFFER (float, w, n);
+
+  F77_XFCN (sch1dn, SCH1DN, (n, chol_mat.fortran_vec (), chol_mat.rows (),
+                             utmp.fortran_vec (), w, info));
+
+  return info;
+}
+
+template <>
+octave_idx_type
+chol<FloatMatrix>::insert_sym (const FloatColumnVector& u, octave_idx_type j)
+{
+  octave_idx_type info = -1;
+
+  octave_idx_type n = chol_mat.rows ();
+
+  if (u.numel () != n + 1)
+    (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
+  if (j < 0 || j > n)
+    (*current_liboctave_error_handler) ("cholinsert: index out of range");
+
+  FloatColumnVector utmp = u;
+
+  OCTAVE_LOCAL_BUFFER (float, w, n);
+
+  chol_mat.resize (n+1, n+1);
+
+  F77_XFCN (schinx, SCHINX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
+                             j + 1, utmp.fortran_vec (), w, info));
+
+  return info;
+}
+
+template <>
+void
+chol<FloatMatrix>::delete_sym (octave_idx_type j)
+{
+  octave_idx_type n = chol_mat.rows ();
+
+  if (j < 0 || j > n-1)
+    (*current_liboctave_error_handler) ("choldelete: index out of range");
+
+  OCTAVE_LOCAL_BUFFER (float, w, n);
+
+  F77_XFCN (schdex, SCHDEX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
+                             j + 1, w));
+
+  chol_mat.resize (n-1, n-1);
+}
+
+template <>
+void
+chol<FloatMatrix>::shift_sym (octave_idx_type i, octave_idx_type j)
+{
+  octave_idx_type n = chol_mat.rows ();
+
+  if (i < 0 || i > n-1 || j < 0 || j > n-1)
+    (*current_liboctave_error_handler) ("cholshift: index out of range");
+
+  OCTAVE_LOCAL_BUFFER (float, w, 2*n);
+
+  F77_XFCN (schshx, SCHSHX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
+                             i + 1, j + 1, w));
+}
+
+#endif
+
+template <>
+octave_idx_type
+chol<ComplexMatrix>::init (const ComplexMatrix& a, bool upper, bool calc_cond)
+{
+  octave_idx_type a_nr = a.rows ();
+  octave_idx_type a_nc = a.cols ();
+
+  if (a_nr != a_nc)
+    (*current_liboctave_error_handler) ("chol: requires square matrix");
+
+  octave_idx_type n = a_nc;
+  octave_idx_type info;
+
+  is_upper = upper;
+
+  chol_mat.clear (n, n);
+  if (is_upper)
+    for (octave_idx_type j = 0; j < n; j++)
+      {
+        for (octave_idx_type i = 0; i <= j; i++)
+          chol_mat.xelem (i, j) = a(i, j);
+        for (octave_idx_type i = j+1; i < n; i++)
+          chol_mat.xelem (i, j) = 0.0;
+      }
+  else
+    for (octave_idx_type j = 0; j < n; j++)
+      {
+        for (octave_idx_type i = 0; i < j; i++)
+          chol_mat.xelem (i, j) = 0.0;
+        for (octave_idx_type i = j; i < n; i++)
+          chol_mat.xelem (i, j) = a(i, j);
+      }
+  Complex *h = chol_mat.fortran_vec ();
+
+  // Calculate the norm of the matrix, for later use.
+  double anorm = 0;
+  if (calc_cond)
+    anorm = xnorm (a, 1);
+
+  if (is_upper)
+    F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1), n, h, n, info
+                               F77_CHAR_ARG_LEN (1)));
+  else
+    F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), n, h, n, info
+                               F77_CHAR_ARG_LEN (1)));
+
+  xrcond = 0.0;
+  if (info > 0)
+    chol_mat.resize (info - 1, info - 1);
+  else if (calc_cond)
+    {
+      octave_idx_type zpocon_info = 0;
+
+      // Now calculate the condition number for non-singular matrix.
+      Array<Complex> z (dim_vector (2*n, 1));
+      Complex *pz = z.fortran_vec ();
+      Array<double> rz (dim_vector (n, 1));
+      double *prz = rz.fortran_vec ();
+      F77_XFCN (zpocon, ZPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h,
+                                 n, anorm, xrcond, pz, prz, zpocon_info
+                                 F77_CHAR_ARG_LEN (1)));
+
+      if (zpocon_info != 0)
+        info = -1;
+    }
+
+  return info;
+}
+
+#ifdef HAVE_QRUPDATE
+
+template <>
+void
+chol<ComplexMatrix>::update (const ComplexColumnVector& u)
+{
+  octave_idx_type n = chol_mat.rows ();
+
+  if (u.numel () != n)
+    (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
+
+  ComplexColumnVector utmp = u;
+
+  OCTAVE_LOCAL_BUFFER (double, rw, n);
+
+  F77_XFCN (zch1up, ZCH1UP, (n, chol_mat.fortran_vec (), chol_mat.rows (),
+                             utmp.fortran_vec (), rw));
+}
+
+template <>
+octave_idx_type
+chol<ComplexMatrix>::downdate (const ComplexColumnVector& u)
+{
+  octave_idx_type info = -1;
+
+  octave_idx_type n = chol_mat.rows ();
+
+  if (u.numel () != n)
+    (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
+
+  ComplexColumnVector utmp = u;
+
+  OCTAVE_LOCAL_BUFFER (double, rw, n);
+
+  F77_XFCN (zch1dn, ZCH1DN, (n, chol_mat.fortran_vec (), chol_mat.rows (),
+                             utmp.fortran_vec (), rw, info));
+
+  return info;
+}
+
+template <>
+octave_idx_type
+chol<ComplexMatrix>::insert_sym (const ComplexColumnVector& u,
+                                 octave_idx_type j)
+{
+  octave_idx_type info = -1;
+
+  octave_idx_type n = chol_mat.rows ();
+
+  if (u.numel () != n + 1)
+    (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
+  if (j < 0 || j > n)
+    (*current_liboctave_error_handler) ("cholinsert: index out of range");
+
+  ComplexColumnVector utmp = u;
+
+  OCTAVE_LOCAL_BUFFER (double, rw, n);
+
+  chol_mat.resize (n+1, n+1);
+
+  F77_XFCN (zchinx, ZCHINX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
+                             j + 1, utmp.fortran_vec (), rw, info));
+
+  return info;
+}
+
+template <>
+void
+chol<ComplexMatrix>::delete_sym (octave_idx_type j)
+{
+  octave_idx_type n = chol_mat.rows ();
+
+  if (j < 0 || j > n-1)
+    (*current_liboctave_error_handler) ("choldelete: index out of range");
+
+  OCTAVE_LOCAL_BUFFER (double, rw, n);
+
+  F77_XFCN (zchdex, ZCHDEX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
+                             j + 1, rw));
+
+  chol_mat.resize (n-1, n-1);
+}
+
+template <>
+void
+chol<ComplexMatrix>::shift_sym (octave_idx_type i, octave_idx_type j)
+{
+  octave_idx_type n = chol_mat.rows ();
+
+  if (i < 0 || i > n-1 || j < 0 || j > n-1)
+    (*current_liboctave_error_handler) ("cholshift: index out of range");
+
+  OCTAVE_LOCAL_BUFFER (Complex, w, n);
+  OCTAVE_LOCAL_BUFFER (double, rw, n);
+
+  F77_XFCN (zchshx, ZCHSHX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
+                             i + 1, j + 1, w, rw));
+}
+
+#endif
+
+template <>
+octave_idx_type
+chol<FloatComplexMatrix>::init (const FloatComplexMatrix& a, bool upper,
+                                bool calc_cond)
+{
+  octave_idx_type a_nr = a.rows ();
+  octave_idx_type a_nc = a.cols ();
+
+  if (a_nr != a_nc)
+    (*current_liboctave_error_handler) ("chol: requires square matrix");
+
+  octave_idx_type n = a_nc;
+  octave_idx_type info;
+
+  is_upper = upper;
+
+  chol_mat.clear (n, n);
+  if (is_upper)
+    for (octave_idx_type j = 0; j < n; j++)
+      {
+        for (octave_idx_type i = 0; i <= j; i++)
+          chol_mat.xelem (i, j) = a(i, j);
+        for (octave_idx_type i = j+1; i < n; i++)
+          chol_mat.xelem (i, j) = 0.0f;
+      }
+  else
+    for (octave_idx_type j = 0; j < n; j++)
+      {
+        for (octave_idx_type i = 0; i < j; i++)
+          chol_mat.xelem (i, j) = 0.0f;
+        for (octave_idx_type i = j; i < n; i++)
+          chol_mat.xelem (i, j) = a(i, j);
+      }
+  FloatComplex *h = chol_mat.fortran_vec ();
+
+  // Calculate the norm of the matrix, for later use.
+  float anorm = 0;
+  if (calc_cond)
+    anorm = xnorm (a, 1);
+
+  if (is_upper)
+    F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1), n, h, n, info
+                               F77_CHAR_ARG_LEN (1)));
+  else
+    F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), n, h, n, info
+                               F77_CHAR_ARG_LEN (1)));
+
+  xrcond = 0.0;
+  if (info > 0)
+    chol_mat.resize (info - 1, info - 1);
+  else if (calc_cond)
+    {
+      octave_idx_type cpocon_info = 0;
+
+      // Now calculate the condition number for non-singular matrix.
+      Array<FloatComplex> z (dim_vector (2*n, 1));
+      FloatComplex *pz = z.fortran_vec ();
+      Array<float> rz (dim_vector (n, 1));
+      float *prz = rz.fortran_vec ();
+      F77_XFCN (cpocon, CPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h,
+                                 n, anorm, xrcond, pz, prz, cpocon_info
+                                 F77_CHAR_ARG_LEN (1)));
+
+      if (cpocon_info != 0)
+        info = -1;
+    }
+
+  return info;
+}
+
+#ifdef HAVE_QRUPDATE
+
+template <>
+void
+chol<FloatComplexMatrix>::update (const FloatComplexColumnVector& u)
+{
+  octave_idx_type n = chol_mat.rows ();
+
+  if (u.numel () != n)
+    (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
+
+  FloatComplexColumnVector utmp = u;
+
+  OCTAVE_LOCAL_BUFFER (float, rw, n);
+
+  F77_XFCN (cch1up, CCH1UP, (n, chol_mat.fortran_vec (), chol_mat.rows (),
+                             utmp.fortran_vec (), rw));
+}
+
+template <>
+octave_idx_type
+chol<FloatComplexMatrix>::downdate (const FloatComplexColumnVector& u)
+{
+  octave_idx_type info = -1;
+
+  octave_idx_type n = chol_mat.rows ();
+
+  if (u.numel () != n)
+    (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
+
+  FloatComplexColumnVector utmp = u;
+
+  OCTAVE_LOCAL_BUFFER (float, rw, n);
+
+  F77_XFCN (cch1dn, CCH1DN, (n, chol_mat.fortran_vec (), chol_mat.rows (),
+                             utmp.fortran_vec (), rw, info));
+
+  return info;
+}
+
+template <>
+octave_idx_type
+chol<FloatComplexMatrix>::insert_sym (const FloatComplexColumnVector& u,
+                                      octave_idx_type j)
+{
+  octave_idx_type info = -1;
+
+  octave_idx_type n = chol_mat.rows ();
+
+  if (u.numel () != n + 1)
+    (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
+  if (j < 0 || j > n)
+    (*current_liboctave_error_handler) ("cholinsert: index out of range");
+
+  FloatComplexColumnVector utmp = u;
+
+  OCTAVE_LOCAL_BUFFER (float, rw, n);
+
+  chol_mat.resize (n+1, n+1);
+
+  F77_XFCN (cchinx, CCHINX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
+                             j + 1, utmp.fortran_vec (), rw, info));
+
+  return info;
+}
+
+template <>
+void
+chol<FloatComplexMatrix>::delete_sym (octave_idx_type j)
+{
+  octave_idx_type n = chol_mat.rows ();
+
+  if (j < 0 || j > n-1)
+    (*current_liboctave_error_handler) ("choldelete: index out of range");
+
+  OCTAVE_LOCAL_BUFFER (float, rw, n);
+
+  F77_XFCN (cchdex, CCHDEX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
+                             j + 1, rw));
+
+  chol_mat.resize (n-1, n-1);
+}
+
+template <>
+void
+chol<FloatComplexMatrix>::shift_sym (octave_idx_type i, octave_idx_type j)
+{
+  octave_idx_type n = chol_mat.rows ();
+
+  if (i < 0 || i > n-1 || j < 0 || j > n-1)
+    (*current_liboctave_error_handler) ("cholshift: index out of range");
+
+  OCTAVE_LOCAL_BUFFER (FloatComplex, w, n);
+  OCTAVE_LOCAL_BUFFER (float, rw, n);
+
+  F77_XFCN (cchshx, CCHSHX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
+                             i + 1, j + 1, w, rw));
+}
+
+#endif
+
+// Instantiations we need.
+
+template class chol<Matrix>;
+
+template class chol<FloatMatrix>;
+
+template class chol<ComplexMatrix>;
+
+template class chol<FloatComplexMatrix>;
+
+template Matrix
+chol2inv<Matrix> (const Matrix& r);
+
+template ComplexMatrix
+chol2inv<ComplexMatrix> (const ComplexMatrix& r);
+
+template FloatMatrix
+chol2inv<FloatMatrix> (const FloatMatrix& r);
+
+template FloatComplexMatrix
+chol2inv<FloatComplexMatrix> (const FloatComplexMatrix& r);
diff -r f08ae27289e4 -r 3c8a3d35661a liboctave/numeric/chol.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/numeric/chol.h	Tue Feb 16 02:47:29 2016 -0500
@@ -0,0 +1,101 @@
+/*
+
+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_chol_h)
+#define octave_chol_h 1
+
+#include "octave-config.h"
+
+template <typename T>
+class
+chol
+{
+public:
+
+  typedef typename T::column_vector_type VT;
+  typedef typename T::real_elt_type COND_T;
+
+  chol (void) : chol_mat (), xrcond (0) { }
+
+  chol (const T& a, bool upper = true, bool calc_cond = false)
+    : chol_mat (), xrcond (0)
+  {
+    init (a, upper, calc_cond);
+  }
+
+  chol (const T& a, octave_idx_type& info, bool upper = true,
+        bool calc_cond = false)
+    : chol_mat (), xrcond (0)
+  {
+    info = init (a, upper, calc_cond);
+  }
+
+  chol (const chol& a)
+    : chol_mat (a.chol_mat), xrcond (a.xrcond) { }
+
+  chol& operator = (const chol& a)
+  {
+    if (this != &a)
+      {
+        chol_mat = a.chol_mat;
+        xrcond = a.xrcond;
+      }
+
+    return *this;
+  }
+
+  T chol_matrix (void) const { return chol_mat; }
+
+  COND_T rcond (void) const { return xrcond; }
+
+  // Compute the inverse of a matrix using the Cholesky factorization.
+  T inverse (void) const;
+
+  void set (const T& R);
+
+  void update (const VT& u);
+
+  octave_idx_type downdate (const VT& u);
+
+  octave_idx_type insert_sym (const VT& u, octave_idx_type j);
+
+  void delete_sym (octave_idx_type j);
+
+  void shift_sym (octave_idx_type i, octave_idx_type j);
+
+private:
+
+  T chol_mat;
+
+  COND_T xrcond;
+
+  bool is_upper;
+
+  octave_idx_type init (const T& a, bool upper, bool calc_cond);
+};
+
+template <typename T>
+T
+chol2inv (const T& r);
+
+#endif
diff -r f08ae27289e4 -r 3c8a3d35661a liboctave/numeric/dbleCHOL.cc
--- a/liboctave/numeric/dbleCHOL.cc	Tue Feb 16 00:32:29 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,453 +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/>.
-
-*/
-
-#ifdef HAVE_CONFIG_H
-#  include <config.h>
-#endif
-
-#include <vector>
-
-#include "dRowVector.h"
-#include "dbleCHOL.h"
-#include "f77-fcn.h"
-#include "lo-error.h"
-#include "oct-locbuf.h"
-#include "oct-norm.h"
-#ifndef HAVE_QRUPDATE
-#  include "dbleQR.h"
-#endif
-
-extern "C"
-{
-  F77_RET_T
-  F77_FUNC (dpotrf, DPOTRF) (F77_CONST_CHAR_ARG_DECL,
-                             const octave_idx_type&, double*,
-                             const octave_idx_type&, octave_idx_type&
-                             F77_CHAR_ARG_LEN_DECL);
-
-  F77_RET_T
-  F77_FUNC (dpotri, DPOTRI) (F77_CONST_CHAR_ARG_DECL,
-                             const octave_idx_type&, double*,
-                             const octave_idx_type&, octave_idx_type&
-                             F77_CHAR_ARG_LEN_DECL);
-
-  F77_RET_T
-  F77_FUNC (dpocon, DPOCON) (F77_CONST_CHAR_ARG_DECL,
-                             const octave_idx_type&, double*,
-                             const octave_idx_type&, const double&,
-                             double&, double*, octave_idx_type*,
-                             octave_idx_type&
-                             F77_CHAR_ARG_LEN_DECL);
-#ifdef HAVE_QRUPDATE
-
-  F77_RET_T
-  F77_FUNC (dch1up, DCH1UP) (const octave_idx_type&, double*,
-                             const octave_idx_type&, double*, double*);
-
-  F77_RET_T
-  F77_FUNC (dch1dn, DCH1DN) (const octave_idx_type&, double*,
-                             const octave_idx_type&, double*, double*,
-                             octave_idx_type&);
-
-  F77_RET_T
-  F77_FUNC (dchinx, DCHINX) (const octave_idx_type&, double*,
-                             const octave_idx_type&, const octave_idx_type&,
-                             double*, double*, octave_idx_type&);
-
-  F77_RET_T
-  F77_FUNC (dchdex, DCHDEX) (const octave_idx_type&, double*,
-                             const octave_idx_type&, const octave_idx_type&,
-                             double*);
-
-  F77_RET_T
-  F77_FUNC (dchshx, DCHSHX) (const octave_idx_type&, double*,
-                             const octave_idx_type&, const octave_idx_type&,
-                             const octave_idx_type&, double*);
-#endif
-}
-
-octave_idx_type
-CHOL::init (const Matrix& a, bool upper, bool calc_cond)
-{
-  octave_idx_type a_nr = a.rows ();
-  octave_idx_type a_nc = a.cols ();
-
-  if (a_nr != a_nc)
-    (*current_liboctave_error_handler) ("CHOL requires square matrix");
-
-  octave_idx_type n = a_nc;
-  octave_idx_type info;
-
-  is_upper = upper;
-
-  chol_mat.clear (n, n);
-  if (is_upper)
-    for (octave_idx_type j = 0; j < n; j++)
-      {
-        for (octave_idx_type i = 0; i <= j; i++)
-          chol_mat.xelem (i, j) = a(i, j);
-        for (octave_idx_type i = j+1; i < n; i++)
-          chol_mat.xelem (i, j) = 0.0;
-      }
-  else
-    for (octave_idx_type j = 0; j < n; j++)
-      {
-        for (octave_idx_type i = 0; i < j; i++)
-          chol_mat.xelem (i, j) = 0.0;
-        for (octave_idx_type i = j; i < n; i++)
-          chol_mat.xelem (i, j) = a(i, j);
-      }
-  double *h = chol_mat.fortran_vec ();
-
-  // Calculate the norm of the matrix, for later use.
-  double anorm = 0;
-  if (calc_cond)
-    anorm = xnorm (a, 1);
-
-  if (is_upper)
-    F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1), n, h, n, info
-                               F77_CHAR_ARG_LEN (1)));
-  else
-    F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), n, h, n, info
-                               F77_CHAR_ARG_LEN (1)));
-
-  xrcond = 0.0;
-  if (info > 0)
-    chol_mat.resize (info - 1, info - 1);
-  else if (calc_cond)
-    {
-      octave_idx_type dpocon_info = 0;
-
-      // Now calculate the condition number for non-singular matrix.
-      Array<double> z (dim_vector (3*n, 1));
-      double *pz = z.fortran_vec ();
-      Array<octave_idx_type> iz (dim_vector (n, 1));
-      octave_idx_type *piz = iz.fortran_vec ();
-      if (is_upper)
-        F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h,
-                                   n, anorm, xrcond, pz, piz, dpocon_info
-                                   F77_CHAR_ARG_LEN (1)));
-      else
-        F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 ("L", 1), n, h,
-                                   n, anorm, xrcond, pz, piz, dpocon_info
-                                   F77_CHAR_ARG_LEN (1)));
-
-      if (dpocon_info != 0)
-        info = -1;
-    }
-
-  return info;
-}
-
-static Matrix
-chol2inv_internal (const Matrix& r, bool is_upper = true)
-{
-  Matrix retval;
-
-  octave_idx_type r_nr = r.rows ();
-  octave_idx_type r_nc = r.cols ();
-
-  if (r_nr != r_nc)
-    (*current_liboctave_error_handler) ("chol2inv requires square matrix");
-
-  octave_idx_type n = r_nc;
-  octave_idx_type info = 0;
-
-  Matrix tmp = r;
-  double *v = tmp.fortran_vec ();
-
-  if (info == 0)
-    {
-      if (is_upper)
-        F77_XFCN (dpotri, DPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n,
-                                   v, n, info
-                                   F77_CHAR_ARG_LEN (1)));
-      else
-        F77_XFCN (dpotri, DPOTRI, (F77_CONST_CHAR_ARG2 ("L", 1), n,
-                                   v, n, info
-                                   F77_CHAR_ARG_LEN (1)));
-
-      // If someone thinks of a more graceful way of doing this (or
-      // faster for that matter :-)), please let me know!
-
-      if (n > 1)
-        {
-          if (is_upper)
-            for (octave_idx_type j = 0; j < r_nc; j++)
-              for (octave_idx_type i = j+1; i < r_nr; i++)
-                tmp.xelem (i, j) = tmp.xelem (j, i);
-          else
-            for (octave_idx_type j = 0; j < r_nc; j++)
-              for (octave_idx_type i = j+1; i < r_nr; i++)
-                tmp.xelem (j, i) = tmp.xelem (i, j);
-        }
-
-      retval = tmp;
-    }
-
-  return retval;
-}
-
-// Compute the inverse of a matrix using the Cholesky factorization.
-Matrix
-CHOL::inverse (void) const
-{
-  return chol2inv_internal (chol_mat, is_upper);
-}
-
-void
-CHOL::set (const Matrix& R)
-{
-  if (! R.is_square ())
-    (*current_liboctave_error_handler) ("CHOL requires square matrix");
-
-  chol_mat = R;
-}
-
-#ifdef HAVE_QRUPDATE
-
-void
-CHOL::update (const ColumnVector& u)
-{
-  octave_idx_type n = chol_mat.rows ();
-
-  if (u.numel () != n)
-    (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
-
-  ColumnVector utmp = u;
-
-  OCTAVE_LOCAL_BUFFER (double, w, n);
-
-  F77_XFCN (dch1up, DCH1UP, (n, chol_mat.fortran_vec (), chol_mat.rows (),
-                             utmp.fortran_vec (), w));
-}
-
-octave_idx_type
-CHOL::downdate (const ColumnVector& u)
-{
-  octave_idx_type info = -1;
-
-  octave_idx_type n = chol_mat.rows ();
-
-  if (u.numel () != n)
-    (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
-
-  ColumnVector utmp = u;
-
-  OCTAVE_LOCAL_BUFFER (double, w, n);
-
-  F77_XFCN (dch1dn, DCH1DN, (n, chol_mat.fortran_vec (), chol_mat.rows (),
-                             utmp.fortran_vec (), w, info));
-
-  return info;
-}
-
-octave_idx_type
-CHOL::insert_sym (const ColumnVector& u, octave_idx_type j)
-{
-  octave_idx_type info = -1;
-
-  octave_idx_type n = chol_mat.rows ();
-
-  if (u.numel () != n + 1)
-    (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
-  if (j < 0 || j > n)
-    (*current_liboctave_error_handler) ("cholinsert: index out of range");
-
-  ColumnVector utmp = u;
-
-  OCTAVE_LOCAL_BUFFER (double, w, n);
-
-  chol_mat.resize (n+1, n+1);
-
-  F77_XFCN (dchinx, DCHINX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
-                             j + 1, utmp.fortran_vec (), w, info));
-
-  return info;
-}
-
-void
-CHOL::delete_sym (octave_idx_type j)
-{
-  octave_idx_type n = chol_mat.rows ();
-
-  if (j < 0 || j > n-1)
-    (*current_liboctave_error_handler) ("choldelete: index out of range");
-
-  OCTAVE_LOCAL_BUFFER (double, w, n);
-
-  F77_XFCN (dchdex, DCHDEX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
-                             j + 1, w));
-
-  chol_mat.resize (n-1, n-1);
-}
-
-void
-CHOL::shift_sym (octave_idx_type i, octave_idx_type j)
-{
-  octave_idx_type n = chol_mat.rows ();
-
-  if (i < 0 || i > n-1 || j < 0 || j > n-1)
-    (*current_liboctave_error_handler) ("cholshift: index out of range");
-
-  OCTAVE_LOCAL_BUFFER (double, w, 2*n);
-
-  F77_XFCN (dchshx, DCHSHX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
-                             i + 1, j + 1, w));
-}
-
-#else
-
-void
-CHOL::update (const ColumnVector& u)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type n = chol_mat.rows ();
-
-  if (u.numel () != n)
-    (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
-
-  init (chol_mat.transpose () * chol_mat
-        + Matrix (u) * Matrix (u).transpose (), true, false);
-}
-
-static bool
-singular (const Matrix& a)
-{
-  for (octave_idx_type i = 0; i < a.rows (); i++)
-    if (a(i,i) == 0.0) return true;
-  return false;
-}
-
-octave_idx_type
-CHOL::downdate (const ColumnVector& u)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type info = -1;
-
-  octave_idx_type n = chol_mat.rows ();
-
-  if (u.numel () != n)
-    (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
-
-  if (singular (chol_mat))
-    info = 2;
-  else
-    {
-      info = init (chol_mat.transpose () * chol_mat
-                   - Matrix (u) * Matrix (u).transpose (), true, false);
-      if (info) info = 1;
-    }
-
-  return info;
-}
-
-octave_idx_type
-CHOL::insert_sym (const ColumnVector& u, octave_idx_type j)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type info = -1;
-
-  octave_idx_type n = chol_mat.rows ();
-
-  if (u.numel () != n + 1)
-    (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
-  if (j < 0 || j > n)
-    (*current_liboctave_error_handler) ("cholinsert: index out of range");
-
-  if (singular (chol_mat))
-    info = 2;
-  else
-    {
-      Matrix a = chol_mat.transpose () * chol_mat;
-      Matrix a1 (n+1, n+1);
-      for (octave_idx_type k = 0; k < n+1; k++)
-        for (octave_idx_type l = 0; l < n+1; l++)
-          {
-            if (l == j)
-              a1(k, l) = u(k);
-            else if (k == j)
-              a1(k, l) = u(l);
-            else
-              a1(k, l) = a(k < j ? k : k-1, l < j ? l : l-1);
-          }
-      info = init (a1, true, false);
-      if (info) info = 1;
-    }
-
-  return info;
-}
-
-void
-CHOL::delete_sym (octave_idx_type j)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type n = chol_mat.rows ();
-
-  if (j < 0 || j > n-1)
-    (*current_liboctave_error_handler) ("choldelete: index out of range");
-
-  Matrix a = chol_mat.transpose () * chol_mat;
-  a.delete_elements (1, idx_vector (j));
-  a.delete_elements (0, idx_vector (j));
-  init (a, true, false);
-}
-
-void
-CHOL::shift_sym (octave_idx_type i, octave_idx_type j)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type n = chol_mat.rows ();
-
-  if (i < 0 || i > n-1 || j < 0 || j > n-1)
-    (*current_liboctave_error_handler) ("cholshift: index out of range");
-
-  Matrix a = chol_mat.transpose () * chol_mat;
-  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;
-    }
-
-  init (a.index (idx_vector (p), idx_vector (p)), true, false);
-}
-
-#endif
-
-Matrix
-chol2inv (const Matrix& r)
-{
-  return chol2inv_internal (r);
-}
diff -r f08ae27289e4 -r 3c8a3d35661a liboctave/numeric/dbleCHOL.h
--- a/liboctave/numeric/dbleCHOL.h	Tue Feb 16 00:32:29 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,100 +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_dbleCHOL_h)
-#define octave_dbleCHOL_h 1
-
-#include "octave-config.h"
-
-#include <iosfwd>
-
-#include "dMatrix.h"
-#include "dColVector.h"
-
-class
-OCTAVE_API
-CHOL
-{
-public:
-
-  CHOL (void) : chol_mat (), xrcond (0) { }
-
-  CHOL (const Matrix& a, bool upper = true, bool calc_cond = false)
-    : chol_mat (), xrcond (0)
-  {
-    init (a, upper, calc_cond);
-  }
-
-  CHOL (const Matrix& a, octave_idx_type& info, bool upper = true,
-        bool calc_cond = false) : chol_mat (), xrcond (0)
-  {
-    info = init (a, upper, calc_cond);
-  }
-
-  CHOL (const CHOL& a) : chol_mat (a.chol_mat), xrcond (a.xrcond) { }
-
-  CHOL& operator = (const CHOL& a)
-  {
-    if (this != &a)
-      {
-        chol_mat = a.chol_mat;
-        xrcond = a.xrcond;
-      }
-    return *this;
-  }
-
-  Matrix chol_matrix (void) const { return chol_mat; }
-
-  double rcond (void) const { return xrcond; }
-
-  // Compute the inverse of a matrix using the Cholesky factorization.
-  Matrix inverse (void) const;
-
-  void set (const Matrix& R);
-
-  void update (const ColumnVector& u);
-
-  octave_idx_type downdate (const ColumnVector& u);
-
-  octave_idx_type insert_sym (const ColumnVector& u, octave_idx_type j);
-
-  void delete_sym (octave_idx_type j);
-
-  void shift_sym (octave_idx_type i, octave_idx_type j);
-
-  friend OCTAVE_API std::ostream& operator << (std::ostream& os, const CHOL& a);
-
-private:
-
-  Matrix chol_mat;
-
-  double xrcond;
-
-  bool is_upper;
-
-  octave_idx_type init (const Matrix& a, bool upper, bool calc_cond);
-};
-
-Matrix OCTAVE_API chol2inv (const Matrix& r);
-
-#endif
diff -r f08ae27289e4 -r 3c8a3d35661a liboctave/numeric/eigs-base.cc
--- a/liboctave/numeric/eigs-base.cc	Tue Feb 16 00:32:29 2016 -0500
+++ b/liboctave/numeric/eigs-base.cc	Tue Feb 16 02:47:29 2016 -0500
@@ -30,11 +30,10 @@
 #include <iostream>
 
 #include "CSparse.h"
-#include "CmplxCHOL.h"
 #include "CmplxLU.h"
 #include "MatrixType.h"
+#include "chol.h"
 #include "dSparse.h"
-#include "dbleCHOL.h"
 #include "dbleLU.h"
 #include "eigs-base.h"
 #include "f77-fcn.h"
@@ -305,7 +304,7 @@
 make_cholb (Matrix& b, Matrix& bt, ColumnVector& permB)
 {
   octave_idx_type info;
-  CHOL fact (b, info);
+  chol<Matrix> fact (b, info);
   octave_idx_type n = b.cols ();
 
   if (info != 0)
@@ -342,7 +341,7 @@
 make_cholb (ComplexMatrix& b, ComplexMatrix& bt, ColumnVector& permB)
 {
   octave_idx_type info;
-  ComplexCHOL fact (b, info);
+  chol<ComplexMatrix> fact (b, info);
   octave_idx_type n = b.cols ();
 
   if (info != 0)
diff -r f08ae27289e4 -r 3c8a3d35661a liboctave/numeric/fCmplxCHOL.cc
--- a/liboctave/numeric/fCmplxCHOL.cc	Tue Feb 16 00:32:29 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,452 +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/>.
-
-*/
-
-#ifdef HAVE_CONFIG_H
-#  include <config.h>
-#endif
-
-#include <vector>
-
-#include "fMatrix.h"
-#include "fRowVector.h"
-#include "fCmplxCHOL.h"
-#include "f77-fcn.h"
-#include "lo-error.h"
-#include "oct-locbuf.h"
-#include "oct-norm.h"
-#ifndef HAVE_QRUPDATE
-#  include "dbleQR.h"
-#endif
-
-extern "C"
-{
-  F77_RET_T
-  F77_FUNC (cpotrf, CPOTRF) (F77_CONST_CHAR_ARG_DECL,
-                             const octave_idx_type&, FloatComplex*,
-                             const octave_idx_type&, octave_idx_type&
-                             F77_CHAR_ARG_LEN_DECL);
-  F77_RET_T
-  F77_FUNC (cpotri, CPOTRI) (F77_CONST_CHAR_ARG_DECL,
-                             const octave_idx_type&, FloatComplex*,
-                             const octave_idx_type&, octave_idx_type&
-                             F77_CHAR_ARG_LEN_DECL);
-
-  F77_RET_T
-  F77_FUNC (cpocon, CPOCON) (F77_CONST_CHAR_ARG_DECL,
-                             const octave_idx_type&, FloatComplex*,
-                             const octave_idx_type&, const float&,
-                             float&, FloatComplex*, float*, octave_idx_type&
-                             F77_CHAR_ARG_LEN_DECL);
-#ifdef HAVE_QRUPDATE
-
-  F77_RET_T
-  F77_FUNC (cch1up, CCH1UP) (const octave_idx_type&, FloatComplex*,
-                             const octave_idx_type&, FloatComplex*, float*);
-
-  F77_RET_T
-  F77_FUNC (cch1dn, CCH1DN) (const octave_idx_type&, FloatComplex*,
-                             const octave_idx_type&, FloatComplex*,
-                             float*, octave_idx_type&);
-
-  F77_RET_T
-  F77_FUNC (cchinx, CCHINX) (const octave_idx_type&, FloatComplex*,
-                             const octave_idx_type&, const octave_idx_type&,
-                             FloatComplex*, float*, octave_idx_type&);
-
-  F77_RET_T
-  F77_FUNC (cchdex, CCHDEX) (const octave_idx_type&, FloatComplex*,
-                             const octave_idx_type&, const octave_idx_type&,
-                             float*);
-
-  F77_RET_T
-  F77_FUNC (cchshx, CCHSHX) (const octave_idx_type&, FloatComplex*,
-                             const octave_idx_type&, const octave_idx_type&,
-                             const octave_idx_type&, FloatComplex*, float*);
-#endif
-}
-
-octave_idx_type
-FloatComplexCHOL::init (const FloatComplexMatrix& a, bool upper, bool calc_cond)
-{
-  octave_idx_type a_nr = a.rows ();
-  octave_idx_type a_nc = a.cols ();
-
-  if (a_nr != a_nc)
-    (*current_liboctave_error_handler)
-      ("FloatComplexCHOL requires square matrix");
-
-  octave_idx_type n = a_nc;
-  octave_idx_type info;
-
-  is_upper = upper;
-
-  chol_mat.clear (n, n);
-  if (is_upper)
-    for (octave_idx_type j = 0; j < n; j++)
-      {
-        for (octave_idx_type i = 0; i <= j; i++)
-          chol_mat.xelem (i, j) = a(i, j);
-        for (octave_idx_type i = j+1; i < n; i++)
-          chol_mat.xelem (i, j) = 0.0f;
-      }
-  else
-    for (octave_idx_type j = 0; j < n; j++)
-      {
-        for (octave_idx_type i = 0; i < j; i++)
-          chol_mat.xelem (i, j) = 0.0f;
-        for (octave_idx_type i = j; i < n; i++)
-          chol_mat.xelem (i, j) = a(i, j);
-      }
-  FloatComplex *h = chol_mat.fortran_vec ();
-
-  // Calculate the norm of the matrix, for later use.
-  float anorm = 0;
-  if (calc_cond)
-    anorm = xnorm (a, 1);
-
-  if (is_upper)
-    F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1), n, h, n, info
-                               F77_CHAR_ARG_LEN (1)));
-  else
-    F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), n, h, n, info
-                               F77_CHAR_ARG_LEN (1)));
-
-  xrcond = 0.0;
-  if (info > 0)
-    chol_mat.resize (info - 1, info - 1);
-  else if (calc_cond)
-    {
-      octave_idx_type cpocon_info = 0;
-
-      // Now calculate the condition number for non-singular matrix.
-      Array<FloatComplex> z (dim_vector (2*n, 1));
-      FloatComplex *pz = z.fortran_vec ();
-      Array<float> rz (dim_vector (n, 1));
-      float *prz = rz.fortran_vec ();
-      F77_XFCN (cpocon, CPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h,
-                                 n, anorm, xrcond, pz, prz, cpocon_info
-                                 F77_CHAR_ARG_LEN (1)));
-
-      if (cpocon_info != 0)
-        info = -1;
-    }
-
-  return info;
-}
-
-static FloatComplexMatrix
-chol2inv_internal (const FloatComplexMatrix& r, bool is_upper = true)
-{
-  FloatComplexMatrix retval;
-
-  octave_idx_type r_nr = r.rows ();
-  octave_idx_type r_nc = r.cols ();
-
-  if (r_nr != r_nc)
-    (*current_liboctave_error_handler) ("chol2inv requires square matrix");
-
-  octave_idx_type n = r_nc;
-  octave_idx_type info;
-
-  FloatComplexMatrix tmp = r;
-
-  if (is_upper)
-    F77_XFCN (cpotri, CPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n,
-                               tmp.fortran_vec (), n, info
-                               F77_CHAR_ARG_LEN (1)));
-  else
-    F77_XFCN (cpotri, CPOTRI, (F77_CONST_CHAR_ARG2 ("L", 1), n,
-                               tmp.fortran_vec (), n, info
-                               F77_CHAR_ARG_LEN (1)));
-
-  // If someone thinks of a more graceful way of doing this (or
-  // faster for that matter :-)), please let me know!
-
-  if (n > 1)
-    {
-      if (is_upper)
-        for (octave_idx_type j = 0; j < r_nc; j++)
-          for (octave_idx_type i = j+1; i < r_nr; i++)
-            tmp.xelem (i, j) = tmp.xelem (j, i);
-      else
-        for (octave_idx_type j = 0; j < r_nc; j++)
-          for (octave_idx_type i = j+1; i < r_nr; i++)
-            tmp.xelem (j, i) = tmp.xelem (i, j);
-    }
-
-  retval = tmp;
-
-  return retval;
-}
-
-// Compute the inverse of a matrix using the Cholesky factorization.
-FloatComplexMatrix
-FloatComplexCHOL::inverse (void) const
-{
-  return chol2inv_internal (chol_mat, is_upper);
-}
-
-void
-FloatComplexCHOL::set (const FloatComplexMatrix& R)
-{
-  if (! R.is_square ())
-    (*current_liboctave_error_handler) ("CHOL requires square matrix");
-
-  chol_mat = R;
-}
-
-#ifdef HAVE_QRUPDATE
-
-void
-FloatComplexCHOL::update (const FloatComplexColumnVector& u)
-{
-  octave_idx_type n = chol_mat.rows ();
-
-  if (u.numel () != n)
-    (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
-
-  FloatComplexColumnVector utmp = u;
-
-  OCTAVE_LOCAL_BUFFER (float, rw, n);
-
-  F77_XFCN (cch1up, CCH1UP, (n, chol_mat.fortran_vec (), chol_mat.rows (),
-                             utmp.fortran_vec (), rw));
-}
-
-octave_idx_type
-FloatComplexCHOL::downdate (const FloatComplexColumnVector& u)
-{
-  octave_idx_type info = -1;
-
-  octave_idx_type n = chol_mat.rows ();
-
-  if (u.numel () != n)
-    (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
-
-  FloatComplexColumnVector utmp = u;
-
-  OCTAVE_LOCAL_BUFFER (float, rw, n);
-
-  F77_XFCN (cch1dn, CCH1DN, (n, chol_mat.fortran_vec (), chol_mat.rows (),
-                             utmp.fortran_vec (), rw, info));
-
-  return info;
-}
-
-octave_idx_type
-FloatComplexCHOL::insert_sym (const FloatComplexColumnVector& u,
-                              octave_idx_type j)
-{
-  octave_idx_type info = -1;
-
-  octave_idx_type n = chol_mat.rows ();
-
-  if (u.numel () != n + 1)
-    (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
-  if (j < 0 || j > n)
-    (*current_liboctave_error_handler) ("cholinsert: index out of range");
-
-  FloatComplexColumnVector utmp = u;
-
-  OCTAVE_LOCAL_BUFFER (float, rw, n);
-
-  chol_mat.resize (n+1, n+1);
-
-  F77_XFCN (cchinx, CCHINX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
-                             j + 1, utmp.fortran_vec (), rw, info));
-
-  return info;
-}
-
-void
-FloatComplexCHOL::delete_sym (octave_idx_type j)
-{
-  octave_idx_type n = chol_mat.rows ();
-
-  if (j < 0 || j > n-1)
-    (*current_liboctave_error_handler) ("choldelete: index out of range");
-
-  OCTAVE_LOCAL_BUFFER (float, rw, n);
-
-  F77_XFCN (cchdex, CCHDEX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
-                             j + 1, rw));
-
-  chol_mat.resize (n-1, n-1);
-}
-
-void
-FloatComplexCHOL::shift_sym (octave_idx_type i, octave_idx_type j)
-{
-  octave_idx_type n = chol_mat.rows ();
-
-  if (i < 0 || i > n-1 || j < 0 || j > n-1)
-    (*current_liboctave_error_handler) ("cholshift: index out of range");
-
-  OCTAVE_LOCAL_BUFFER (FloatComplex, w, n);
-  OCTAVE_LOCAL_BUFFER (float, rw, n);
-
-  F77_XFCN (cchshx, CCHSHX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
-                             i + 1, j + 1, w, rw));
-}
-
-#else
-
-void
-FloatComplexCHOL::update (const FloatComplexColumnVector& u)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type n = chol_mat.rows ();
-
-  if (u.numel () != n)
-    (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
-
-  init (chol_mat.hermitian () * chol_mat
-        + FloatComplexMatrix (u) * FloatComplexMatrix (u).hermitian (),
-        true, false);
-}
-
-static bool
-singular (const FloatComplexMatrix& a)
-{
-  for (octave_idx_type i = 0; i < a.rows (); i++)
-    if (a(i,i) == 0.0f) return true;
-  return false;
-}
-
-octave_idx_type
-FloatComplexCHOL::downdate (const FloatComplexColumnVector& u)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type info = -1;
-
-  octave_idx_type n = chol_mat.rows ();
-
-  if (u.numel () != n)
-    (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
-
-  if (singular (chol_mat))
-    info = 2;
-  else
-    {
-      info = init (chol_mat.hermitian () * chol_mat
-                   - FloatComplexMatrix (u)
-                   * FloatComplexMatrix (u).hermitian (),
-                   true, false);
-      if (info) info = 1;
-    }
-
-  return info;
-}
-
-octave_idx_type
-FloatComplexCHOL::insert_sym (const FloatComplexColumnVector& u,
-                              octave_idx_type j)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type info = -1;
-
-  octave_idx_type n = chol_mat.rows ();
-
-  if (u.numel () != n + 1)
-    (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
-  if (j < 0 || j > n)
-    (*current_liboctave_error_handler) ("cholinsert: index out of range");
-
-  if (singular (chol_mat))
-    info = 2;
-  else if (u(j).imag () != 0.0)
-    info = 3;
-  else
-    {
-      FloatComplexMatrix a = chol_mat.hermitian () * chol_mat;
-      FloatComplexMatrix a1 (n+1, n+1);
-      for (octave_idx_type k = 0; k < n+1; k++)
-        for (octave_idx_type l = 0; l < n+1; l++)
-          {
-            if (l == j)
-              a1(k, l) = u(k);
-            else if (k == j)
-              a1(k, l) = std::conj (u(l));
-            else
-              a1(k, l) = a(k < j ? k : k-1, l < j ? l : l-1);
-          }
-      info = init (a1, true, false);
-      if (info) info = 1;
-    }
-
-  return info;
-}
-
-void
-FloatComplexCHOL::delete_sym (octave_idx_type j)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type n = chol_mat.rows ();
-
-  if (j < 0 || j > n-1)
-    (*current_liboctave_error_handler) ("choldelete: index out of range");
-
-  FloatComplexMatrix a = chol_mat.hermitian () * chol_mat;
-  a.delete_elements (1, idx_vector (j));
-  a.delete_elements (0, idx_vector (j));
-  init (a, true, false);
-}
-
-void
-FloatComplexCHOL::shift_sym (octave_idx_type i, octave_idx_type j)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type n = chol_mat.rows ();
-
-  if (i < 0 || i > n-1 || j < 0 || j > n-1)
-    (*current_liboctave_error_handler) ("cholshift: index out of range");
-
-  FloatComplexMatrix a = chol_mat.hermitian () * chol_mat;
-  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;
-    }
-
-  init (a.index (idx_vector (p), idx_vector (p)), true, false);
-}
-
-#endif
-
-FloatComplexMatrix
-chol2inv (const FloatComplexMatrix& r)
-{
-  return chol2inv_internal (r);
-}
diff -r f08ae27289e4 -r 3c8a3d35661a liboctave/numeric/fCmplxCHOL.h
--- a/liboctave/numeric/fCmplxCHOL.h	Tue Feb 16 00:32:29 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,105 +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_fCmplxCHOL_h)
-#define octave_fCmplxCHOL_h 1
-
-#include "octave-config.h"
-
-#include <iosfwd>
-
-#include "fCMatrix.h"
-#include "fCColVector.h"
-
-class
-OCTAVE_API
-FloatComplexCHOL
-{
-public:
-
-  FloatComplexCHOL (void) : chol_mat (), xrcond (0) { }
-
-  FloatComplexCHOL (const FloatComplexMatrix& a, bool upper = true,
-                    bool calc_cond = false)
-    : chol_mat (), xrcond (0)
-  {
-    init (a, upper, calc_cond);
-  }
-
-  FloatComplexCHOL (const FloatComplexMatrix& a, octave_idx_type& info,
-                    bool upper = true, bool calc_cond = false)
-    : chol_mat (), xrcond (0)
-  {
-    info = init (a, upper, calc_cond);
-  }
-
-  FloatComplexCHOL (const FloatComplexCHOL& a)
-    : chol_mat (a.chol_mat), xrcond (a.xrcond) { }
-
-  FloatComplexCHOL& operator = (const FloatComplexCHOL& a)
-  {
-    if (this != &a)
-      {
-        chol_mat = a.chol_mat;
-        xrcond = a.xrcond;
-      }
-
-    return *this;
-  }
-
-  FloatComplexMatrix chol_matrix (void) const { return chol_mat; }
-
-  float rcond (void) const { return xrcond; }
-
-  FloatComplexMatrix inverse (void) const;
-
-  void set (const FloatComplexMatrix& R);
-
-  void update (const FloatComplexColumnVector& u);
-
-  octave_idx_type downdate (const FloatComplexColumnVector& u);
-
-  octave_idx_type insert_sym (const FloatComplexColumnVector& u,
-                              octave_idx_type j);
-
-  void delete_sym (octave_idx_type j);
-
-  void shift_sym (octave_idx_type i, octave_idx_type j);
-
-  friend OCTAVE_API std::ostream& operator << (std::ostream& os,
-                                               const FloatComplexCHOL& a);
-
-private:
-
-  FloatComplexMatrix chol_mat;
-
-  float xrcond;
-
-  bool is_upper;
-
-  octave_idx_type init (const FloatComplexMatrix& a, bool upper, bool calc_cond);
-};
-
-FloatComplexMatrix OCTAVE_API chol2inv (const FloatComplexMatrix& r);
-
-#endif
diff -r f08ae27289e4 -r 3c8a3d35661a liboctave/numeric/floatCHOL.cc
--- a/liboctave/numeric/floatCHOL.cc	Tue Feb 16 00:32:29 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,454 +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/>.
-
-*/
-
-#ifdef HAVE_CONFIG_H
-#  include <config.h>
-#endif
-
-#include <vector>
-
-#include "fRowVector.h"
-#include "floatCHOL.h"
-#include "f77-fcn.h"
-#include "lo-error.h"
-#include "oct-locbuf.h"
-#include "oct-norm.h"
-#ifndef HAVE_QRUPDATE
-#  include "dbleQR.h"
-#endif
-
-extern "C"
-{
-  F77_RET_T
-  F77_FUNC (spotrf, SPOTRF) (F77_CONST_CHAR_ARG_DECL,
-                             const octave_idx_type&, float*,
-                             const octave_idx_type&, octave_idx_type&
-                             F77_CHAR_ARG_LEN_DECL);
-
-  F77_RET_T
-  F77_FUNC (spotri, SPOTRI) (F77_CONST_CHAR_ARG_DECL,
-                             const octave_idx_type&, float*,
-                             const octave_idx_type&, octave_idx_type&
-                             F77_CHAR_ARG_LEN_DECL);
-
-  F77_RET_T
-  F77_FUNC (spocon, SPOCON) (F77_CONST_CHAR_ARG_DECL,
-                             const octave_idx_type&, float*,
-                             const octave_idx_type&, const float&,
-                             float&, float*, octave_idx_type*,
-                             octave_idx_type&
-                             F77_CHAR_ARG_LEN_DECL);
-#ifdef HAVE_QRUPDATE
-
-  F77_RET_T
-  F77_FUNC (sch1up, SCH1UP) (const octave_idx_type&, float*,
-                             const octave_idx_type&, float*, float*);
-
-  F77_RET_T
-  F77_FUNC (sch1dn, SCH1DN) (const octave_idx_type&, float*,
-                             const octave_idx_type&, float*, float*,
-                             octave_idx_type&);
-
-  F77_RET_T
-  F77_FUNC (schinx, SCHINX) (const octave_idx_type&, float*,
-                             const octave_idx_type&, const octave_idx_type&,
-                             float*, float*, octave_idx_type&);
-
-  F77_RET_T
-  F77_FUNC (schdex, SCHDEX) (const octave_idx_type&, float*,
-                             const octave_idx_type&, const octave_idx_type&,
-                             float*);
-
-  F77_RET_T
-  F77_FUNC (schshx, SCHSHX) (const octave_idx_type&, float*,
-                             const octave_idx_type&, const octave_idx_type&,
-                             const octave_idx_type&, float*);
-#endif
-}
-
-octave_idx_type
-FloatCHOL::init (const FloatMatrix& a, bool upper, bool calc_cond)
-{
-  octave_idx_type a_nr = a.rows ();
-  octave_idx_type a_nc = a.cols ();
-
-  if (a_nr != a_nc)
-    (*current_liboctave_error_handler) ("FloatCHOL requires square matrix");
-
-  octave_idx_type n = a_nc;
-  octave_idx_type info;
-
-  is_upper = upper;
-
-  chol_mat.clear (n, n);
-  if (is_upper)
-    for (octave_idx_type j = 0; j < n; j++)
-      {
-        for (octave_idx_type i = 0; i <= j; i++)
-          chol_mat.xelem (i, j) = a(i, j);
-        for (octave_idx_type i = j+1; i < n; i++)
-          chol_mat.xelem (i, j) = 0.0f;
-      }
-  else
-    for (octave_idx_type j = 0; j < n; j++)
-      {
-        for (octave_idx_type i = 0; i < j; i++)
-          chol_mat.xelem (i, j) = 0.0f;
-        for (octave_idx_type i = j; i < n; i++)
-          chol_mat.xelem (i, j) = a(i, j);
-      }
-  float *h = chol_mat.fortran_vec ();
-
-  // Calculate the norm of the matrix, for later use.
-  float anorm = 0;
-  if (calc_cond)
-    anorm = xnorm (a, 1);
-
-  if (is_upper)
-    F77_XFCN (spotrf, SPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1), n, h, n, info
-                               F77_CHAR_ARG_LEN (1)));
-  else
-    F77_XFCN (spotrf, SPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), n, h, n, info
-                               F77_CHAR_ARG_LEN (1)));
-
-  xrcond = 0.0;
-  if (info > 0)
-    chol_mat.resize (info - 1, info - 1);
-  else if (calc_cond)
-    {
-      octave_idx_type spocon_info = 0;
-
-      // Now calculate the condition number for non-singular matrix.
-      Array<float> z (dim_vector (3*n, 1));
-      float *pz = z.fortran_vec ();
-      Array<octave_idx_type> iz (dim_vector (n, 1));
-      octave_idx_type *piz = iz.fortran_vec ();
-      if (is_upper)
-        F77_XFCN (spocon, SPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h,
-                                   n, anorm, xrcond, pz, piz, spocon_info
-                                   F77_CHAR_ARG_LEN (1)));
-      else
-        F77_XFCN (spocon, SPOCON, (F77_CONST_CHAR_ARG2 ("L", 1), n, h,
-                                   n, anorm, xrcond, pz, piz, spocon_info
-                                   F77_CHAR_ARG_LEN (1)));
-
-      if (spocon_info != 0)
-        info = -1;
-    }
-
-  return info;
-}
-
-static FloatMatrix
-chol2inv_internal (const FloatMatrix& r, bool is_upper = true)
-{
-  FloatMatrix retval;
-
-  octave_idx_type r_nr = r.rows ();
-  octave_idx_type r_nc = r.cols ();
-
-  if (r_nr != r_nc)
-    (*current_liboctave_error_handler) ("chol2inv requires square matrix");
-
-  octave_idx_type n = r_nc;
-  octave_idx_type info = 0;
-
-  FloatMatrix tmp = r;
-  float *v = tmp.fortran_vec ();
-
-  if (info == 0)
-    {
-      if (is_upper)
-        F77_XFCN (spotri, SPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n,
-                                   v, n, info
-                                   F77_CHAR_ARG_LEN (1)));
-      else
-        F77_XFCN (spotri, SPOTRI, (F77_CONST_CHAR_ARG2 ("L", 1), n,
-                                   v, n, info
-                                   F77_CHAR_ARG_LEN (1)));
-
-      // If someone thinks of a more graceful way of doing this (or
-      // faster for that matter :-)), please let me know!
-
-      if (n > 1)
-        {
-          if (is_upper)
-            for (octave_idx_type j = 0; j < r_nc; j++)
-              for (octave_idx_type i = j+1; i < r_nr; i++)
-                tmp.xelem (i, j) = tmp.xelem (j, i);
-          else
-            for (octave_idx_type j = 0; j < r_nc; j++)
-              for (octave_idx_type i = j+1; i < r_nr; i++)
-                tmp.xelem (j, i) = tmp.xelem (i, j);
-        }
-
-      retval = tmp;
-    }
-
-  return retval;
-}
-
-// Compute the inverse of a matrix using the Cholesky factorization.
-FloatMatrix
-FloatCHOL::inverse (void) const
-{
-  return chol2inv_internal (chol_mat, is_upper);
-}
-
-void
-FloatCHOL::set (const FloatMatrix& R)
-{
-  if (! R.is_square ())
-    (*current_liboctave_error_handler) ("FloatCHOL requires square matrix");
-
-  chol_mat = R;
-}
-
-#ifdef HAVE_QRUPDATE
-
-void
-FloatCHOL::update (const FloatColumnVector& u)
-{
-  octave_idx_type n = chol_mat.rows ();
-
-  if (u.numel () != n)
-    (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
-
-  FloatColumnVector utmp = u;
-
-  OCTAVE_LOCAL_BUFFER (float, w, n);
-
-  F77_XFCN (sch1up, SCH1UP, (n, chol_mat.fortran_vec (), chol_mat.rows (),
-                             utmp.fortran_vec (), w));
-}
-
-octave_idx_type
-FloatCHOL::downdate (const FloatColumnVector& u)
-{
-  octave_idx_type info = -1;
-
-  octave_idx_type n = chol_mat.rows ();
-
-  if (u.numel () != n)
-    (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
-
-  FloatColumnVector utmp = u;
-
-  OCTAVE_LOCAL_BUFFER (float, w, n);
-
-  F77_XFCN (sch1dn, SCH1DN, (n, chol_mat.fortran_vec (), chol_mat.rows (),
-                             utmp.fortran_vec (), w, info));
-
-  return info;
-}
-
-octave_idx_type
-FloatCHOL::insert_sym (const FloatColumnVector& u, octave_idx_type j)
-{
-  octave_idx_type info = -1;
-
-  octave_idx_type n = chol_mat.rows ();
-
-  if (u.numel () != n + 1)
-    (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
-  if (j < 0 || j > n)
-    (*current_liboctave_error_handler) ("cholinsert: index out of range");
-
-  FloatColumnVector utmp = u;
-
-  OCTAVE_LOCAL_BUFFER (float, w, n);
-
-  chol_mat.resize (n+1, n+1);
-
-  F77_XFCN (schinx, SCHINX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
-                             j + 1, utmp.fortran_vec (), w, info));
-
-  return info;
-}
-
-void
-FloatCHOL::delete_sym (octave_idx_type j)
-{
-  octave_idx_type n = chol_mat.rows ();
-
-  if (j < 0 || j > n-1)
-    (*current_liboctave_error_handler) ("choldelete: index out of range");
-
-  OCTAVE_LOCAL_BUFFER (float, w, n);
-
-  F77_XFCN (schdex, SCHDEX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
-                             j + 1, w));
-
-  chol_mat.resize (n-1, n-1);
-}
-
-void
-FloatCHOL::shift_sym (octave_idx_type i, octave_idx_type j)
-{
-  octave_idx_type n = chol_mat.rows ();
-
-  if (i < 0 || i > n-1 || j < 0 || j > n-1)
-    (*current_liboctave_error_handler) ("cholshift: index out of range");
-
-  OCTAVE_LOCAL_BUFFER (float, w, 2*n);
-
-  F77_XFCN (schshx, SCHSHX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
-                             i + 1, j + 1, w));
-}
-
-#else
-
-void
-FloatCHOL::update (const FloatColumnVector& u)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type n = chol_mat.rows ();
-
-  if (u.numel () != n)
-    (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
-
-  init (chol_mat.transpose () * chol_mat
-        + FloatMatrix (u) * FloatMatrix (u).transpose (), true, false);
-}
-
-static bool
-singular (const FloatMatrix& a)
-{
-  for (octave_idx_type i = 0; i < a.rows (); i++)
-    if (a(i,i) == 0.0f) return true;
-  return false;
-}
-
-octave_idx_type
-FloatCHOL::downdate (const FloatColumnVector& u)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type info = -1;
-
-  octave_idx_type n = chol_mat.rows ();
-
-  if (u.numel () != n)
-    (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
-
-  if (singular (chol_mat))
-    info = 2;
-  else
-    {
-      info = init (chol_mat.transpose () * chol_mat
-                   - FloatMatrix (u) * FloatMatrix (u).transpose (), true,
-                   false);
-      if (info) info = 1;
-    }
-
-  return info;
-}
-
-octave_idx_type
-FloatCHOL::insert_sym (const FloatColumnVector& u, octave_idx_type j)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type info = -1;
-
-  octave_idx_type n = chol_mat.rows ();
-
-  if (u.numel () != n + 1)
-    (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
-  if (j < 0 || j > n)
-    (*current_liboctave_error_handler) ("cholinsert: index out of range");
-
-  if (singular (chol_mat))
-    info = 2;
-  else
-    {
-      FloatMatrix a = chol_mat.transpose () * chol_mat;
-      FloatMatrix a1 (n+1, n+1);
-      for (octave_idx_type k = 0; k < n+1; k++)
-        for (octave_idx_type l = 0; l < n+1; l++)
-          {
-            if (l == j)
-              a1(k, l) = u(k);
-            else if (k == j)
-              a1(k, l) = u(l);
-            else
-              a1(k, l) = a(k < j ? k : k-1, l < j ? l : l-1);
-          }
-      info = init (a1, true, false);
-      if (info) info = 1;
-    }
-
-  return info;
-}
-
-void
-FloatCHOL::delete_sym (octave_idx_type j)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type n = chol_mat.rows ();
-
-  if (j < 0 || j > n-1)
-    (*current_liboctave_error_handler) ("choldelete: index out of range");
-
-  FloatMatrix a = chol_mat.transpose () * chol_mat;
-  a.delete_elements (1, idx_vector (j));
-  a.delete_elements (0, idx_vector (j));
-  init (a, true, false);
-}
-
-void
-FloatCHOL::shift_sym (octave_idx_type i, octave_idx_type j)
-{
-  warn_qrupdate_once ();
-
-  octave_idx_type n = chol_mat.rows ();
-
-  if (i < 0 || i > n-1 || j < 0 || j > n-1)
-    (*current_liboctave_error_handler) ("cholshift: index out of range");
-
-  FloatMatrix a = chol_mat.transpose () * chol_mat;
-  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;
-    }
-
-  init (a.index (idx_vector (p), idx_vector (p)), true, false);
-}
-
-#endif
-
-FloatMatrix
-chol2inv (const FloatMatrix& r)
-{
-  return chol2inv_internal (r);
-}
diff -r f08ae27289e4 -r 3c8a3d35661a liboctave/numeric/floatCHOL.h
--- a/liboctave/numeric/floatCHOL.h	Tue Feb 16 00:32:29 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,102 +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_floatCHOL_h)
-#define octave_floatCHOL_h 1
-
-#include "octave-config.h"
-
-#include <iosfwd>
-
-#include "fMatrix.h"
-#include "fColVector.h"
-
-class
-OCTAVE_API
-FloatCHOL
-{
-public:
-
-  FloatCHOL (void) : chol_mat (), xrcond (0) { }
-
-  FloatCHOL (const FloatMatrix& a, bool upper = true, bool calc_cond = false)
-    : chol_mat (), xrcond (0)
-  {
-    init (a, upper, calc_cond);
-  }
-
-  FloatCHOL (const FloatMatrix& a, octave_idx_type& info, bool upper = true,
-             bool calc_cond = false)
-    : chol_mat (), xrcond (0)
-  {
-    info = init (a, upper, calc_cond);
-  }
-
-  FloatCHOL (const FloatCHOL& a) : chol_mat (a.chol_mat), xrcond (a.xrcond) { }
-
-  FloatCHOL& operator = (const FloatCHOL& a)
-  {
-    if (this != &a)
-      {
-        chol_mat = a.chol_mat;
-        xrcond = a.xrcond;
-      }
-    return *this;
-  }
-
-  FloatMatrix chol_matrix (void) const { return chol_mat; }
-
-  float rcond (void) const { return xrcond; }
-
-  // Compute the inverse of a matrix using the Cholesky factorization.
-  FloatMatrix inverse (void) const;
-
-  void set (const FloatMatrix& R);
-
-  void update (const FloatColumnVector& u);
-
-  octave_idx_type downdate (const FloatColumnVector& u);
-
-  octave_idx_type insert_sym (const FloatColumnVector& u, octave_idx_type j);
-
-  void delete_sym (octave_idx_type j);
-
-  void shift_sym (octave_idx_type i, octave_idx_type j);
-
-  friend OCTAVE_API std::ostream& operator << (std::ostream& os,
-                                               const FloatCHOL& a);
-
-private:
-
-  FloatMatrix chol_mat;
-
-  float xrcond;
-
-  bool is_upper;
-
-  octave_idx_type init (const FloatMatrix& a, bool upper, bool calc_cond);
-};
-
-FloatMatrix OCTAVE_API chol2inv (const FloatMatrix& r);
-
-#endif
diff -r f08ae27289e4 -r 3c8a3d35661a liboctave/numeric/module.mk
--- a/liboctave/numeric/module.mk	Tue Feb 16 00:32:29 2016 -0500
+++ b/liboctave/numeric/module.mk	Tue Feb 16 02:47:29 2016 -0500
@@ -8,7 +8,6 @@
 LIBOCTAVE_OPT_IN = $(LIBOCTAVE_OPT_INC:.h=.in)
 
 NUMERIC_INC = \
-  liboctave/numeric/CmplxCHOL.h \
   liboctave/numeric/CmplxLU.h \
   liboctave/numeric/CmplxQR.h \
   liboctave/numeric/CmplxQRP.h \
@@ -37,19 +36,17 @@
   liboctave/numeric/base-qr.h \
   liboctave/numeric/bsxfun-decl.h \
   liboctave/numeric/bsxfun.h \
-  liboctave/numeric/dbleCHOL.h \
+  liboctave/numeric/chol.h \
   liboctave/numeric/dbleLU.h \
   liboctave/numeric/dbleQR.h \
   liboctave/numeric/dbleQRP.h \
   liboctave/numeric/dbleSVD.h \
   liboctave/numeric/eigs-base.h \
-  liboctave/numeric/fCmplxCHOL.h \
   liboctave/numeric/fCmplxLU.h \
   liboctave/numeric/fCmplxQR.h \
   liboctave/numeric/fCmplxQRP.h \
   liboctave/numeric/fCmplxSVD.h \
   liboctave/numeric/fEIG.h \
-  liboctave/numeric/floatCHOL.h \
   liboctave/numeric/floatLU.h \
   liboctave/numeric/floatQR.h \
   liboctave/numeric/floatQRP.h \
@@ -78,7 +75,6 @@
   liboctave/numeric/randpoisson.c
 
 NUMERIC_SRC = \
-  liboctave/numeric/CmplxCHOL.cc \
   liboctave/numeric/CmplxLU.cc \
   liboctave/numeric/CmplxQR.cc \
   liboctave/numeric/CmplxQRP.cc \
@@ -92,19 +88,17 @@
   liboctave/numeric/ODES.cc \
   liboctave/numeric/Quad.cc \
   liboctave/numeric/aepbalance.cc \
-  liboctave/numeric/dbleCHOL.cc \
+  liboctave/numeric/chol.cc \
   liboctave/numeric/dbleLU.cc \
   liboctave/numeric/dbleQR.cc \
   liboctave/numeric/dbleQRP.cc \
   liboctave/numeric/dbleSVD.cc \
   liboctave/numeric/eigs-base.cc \
-  liboctave/numeric/fCmplxCHOL.cc \
   liboctave/numeric/fCmplxLU.cc \
   liboctave/numeric/fCmplxQR.cc \
   liboctave/numeric/fCmplxQRP.cc \
   liboctave/numeric/fCmplxSVD.cc \
   liboctave/numeric/fEIG.cc \
-  liboctave/numeric/floatCHOL.cc \
   liboctave/numeric/floatLU.cc \
   liboctave/numeric/floatQR.cc \
   liboctave/numeric/floatQRP.cc \
diff -r f08ae27289e4 -r 3c8a3d35661a liboctave/operators/mx-defs.h
--- a/liboctave/operators/mx-defs.h	Tue Feb 16 00:32:29 2016 -0500
+++ b/liboctave/operators/mx-defs.h	Tue Feb 16 02:47:29 2016 -0500
@@ -62,10 +62,7 @@
 
 template <typename T> class gepbalance;
 
-class CHOL;
-class ComplexCHOL;
-class FloatCHOL;
-class FloatComplexCHOL;
+template <typename T> class chol;
 
 class EIG;
 
diff -r f08ae27289e4 -r 3c8a3d35661a liboctave/operators/mx-ext.h
--- a/liboctave/operators/mx-ext.h	Tue Feb 16 00:32:29 2016 -0500
+++ b/liboctave/operators/mx-ext.h	Tue Feb 16 02:47:29 2016 -0500
@@ -39,10 +39,7 @@
 
 // Result of a Cholesky Factorization
 
-#include "dbleCHOL.h"
-#include "CmplxCHOL.h"
-#include "floatCHOL.h"
-#include "fCmplxCHOL.h"
+#include "chol.h"
 
 // Result of a Hessenberg Decomposition