Mercurial > octave
changeset 21269:3c8a3d35661a
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.
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Tue, 16 Feb 2016 02:47:29 -0500 |
| parents | f08ae27289e4 |
| children | 230e186e292d |
| files | libinterp/corefcn/__qp__.cc libinterp/dldfcn/chol.cc liboctave/array/CMatrix.cc liboctave/array/CMatrix.h liboctave/array/dMatrix.cc liboctave/array/dMatrix.h liboctave/array/fCMatrix.cc liboctave/array/fCMatrix.h liboctave/array/fMatrix.cc liboctave/array/fMatrix.h liboctave/numeric/CmplxCHOL.cc liboctave/numeric/CmplxCHOL.h liboctave/numeric/chol.cc liboctave/numeric/chol.h liboctave/numeric/dbleCHOL.cc liboctave/numeric/dbleCHOL.h liboctave/numeric/eigs-base.cc liboctave/numeric/fCmplxCHOL.cc liboctave/numeric/fCmplxCHOL.h liboctave/numeric/floatCHOL.cc liboctave/numeric/floatCHOL.h liboctave/numeric/module.mk liboctave/operators/mx-defs.h liboctave/operators/mx-ext.h |
| diffstat | 24 files changed, 1454 insertions(+), 2279 deletions(-) [+] |
line wrap: on
line diff
--- 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;
--- 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);
--- 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)
--- 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 () { }
--- 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)
--- 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 () { }
--- 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)
--- 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 () { }
--- 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)
--- 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 () { }
--- 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); -}
--- 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
--- /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);
--- /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
--- 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); -}
--- 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
--- 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)
--- 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); -}
--- 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
--- 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); -}
--- 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
--- 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 \
--- 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;
--- 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