Mercurial > octave-nkf
changeset 20497:5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
* chol.cc (chol): Send 'L' parameter correctly when chol is called with 'lower'.
* floatCHOL.cc (init): Propagate 'lower' to underlying library function.
* floatCHOL.h: Modify the prototype of methods.
* fMatrix.cc (inverse): Invoke chol with additional parameter.
* dbleCHOL.cc (init): Propagate 'lower' to underlying library function.
* dbleCHOL.h: Modify the prototype of methods.
* dMatrix.cc (inverse): Invoke chol with additional parameter.
* CmplxCHOL.cc (init): Propagate 'lower' to underlying library function.
* CmplxCHOL.h: Modify the prototype of methods.
* CMatrix.cc (inverse): Invoke chol with additional parameter.
| author | PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com> |
|---|---|
| date | Sun, 24 Aug 2014 19:35:06 +0530 |
| parents | 0fe7133da8ce |
| children | 7fbba8c8efd5 |
| files | libinterp/dldfcn/chol.cc liboctave/array/CMatrix.cc liboctave/array/dMatrix.cc liboctave/array/fMatrix.cc liboctave/numeric/CmplxCHOL.cc liboctave/numeric/CmplxCHOL.h liboctave/numeric/dbleCHOL.cc liboctave/numeric/dbleCHOL.h liboctave/numeric/floatCHOL.cc liboctave/numeric/floatCHOL.h |
| diffstat | 10 files changed, 248 insertions(+), 103 deletions(-) [+] |
line wrap: on
line diff
--- a/libinterp/dldfcn/chol.cc Thu Aug 20 14:37:57 2015 -0400 +++ b/libinterp/dldfcn/chol.cc Sun Aug 24 19:35:06 2014 +0530 @@ -46,6 +46,13 @@ template <class CHOLT> static octave_value +get_chol (const CHOLT& fact) +{ + return octave_value (fact.chol_matrix()); +} + +template <class CHOLT> +static octave_value get_chol_r (const CHOLT& fact) { return octave_value (fact.chol_matrix (), @@ -162,10 +169,6 @@ if (tmp.compare ("vector") == 0) vecout = true; else if (tmp.compare ("lower") == 0) - // FIXME: currently the option "lower" is handled by transposing - // the matrix, factorizing it with the lapack function - // DPOTRF ('U', ...) and finally transposing the factor. It would - // be more efficient to use DPOTRF ('L', ...) in this case. LLt = true; else if (tmp.compare ("upper") == 0) LLt = false; @@ -266,18 +269,12 @@ octave_idx_type info; FloatCHOL fact; - if (LLt) - fact = FloatCHOL (m.transpose (), info); - else - fact = FloatCHOL (m, info); + fact = FloatCHOL (m, info, LLt != true); if (nargout == 2 || info == 0) { retval(1) = info; - if (LLt) - retval(0) = get_chol_l (fact); - else - retval(0) = get_chol_r (fact); + retval(0) = get_chol (fact); } else error ("chol: input matrix must be positive definite"); @@ -323,18 +320,12 @@ octave_idx_type info; CHOL fact; - if (LLt) - fact = CHOL (m.transpose (), info); - else - fact = CHOL (m, info); + fact = CHOL (m, info, LLt != true); if (nargout == 2 || info == 0) { retval(1) = info; - if (LLt) - retval(0) = get_chol_l (fact); - else - retval(0) = get_chol_r (fact); + retval(0) = get_chol (fact); } else error ("chol: input matrix must be positive definite"); @@ -349,18 +340,12 @@ octave_idx_type info; ComplexCHOL fact; - if (LLt) - fact = ComplexCHOL (m.transpose (), info); - else - fact = ComplexCHOL (m, info); + fact = ComplexCHOL (m, info, LLt != true); if (nargout == 2 || info == 0) { retval(1) = info; - if (LLt) - retval(0) = get_chol_l (fact); - else - retval(0) = get_chol_r (fact); + retval(0) = get_chol (fact); } else error ("chol: input matrix must be positive definite");
--- a/liboctave/array/CMatrix.cc Thu Aug 20 14:37:57 2015 -0400 +++ b/liboctave/array/CMatrix.cc Sun Aug 24 19:35:06 2014 +0530 @@ -1193,7 +1193,7 @@ { if (mattype.is_hermitian ()) { - ComplexCHOL chol (*this, info, calc_cond); + ComplexCHOL chol (*this, info, true, calc_cond); if (info == 0) { if (calc_cond)
--- a/liboctave/array/dMatrix.cc Thu Aug 20 14:37:57 2015 -0400 +++ b/liboctave/array/dMatrix.cc Sun Aug 24 19:35:06 2014 +0530 @@ -843,7 +843,7 @@ { if (mattype.is_hermitian ()) { - CHOL chol (*this, info, calc_cond); + CHOL chol (*this, info, true, calc_cond); if (info == 0) { if (calc_cond)
--- a/liboctave/array/fMatrix.cc Thu Aug 20 14:37:57 2015 -0400 +++ b/liboctave/array/fMatrix.cc Sun Aug 24 19:35:06 2014 +0530 @@ -850,7 +850,7 @@ { if (mattype.is_hermitian ()) { - FloatCHOL chol (*this, info, calc_cond); + FloatCHOL chol (*this, info, true, calc_cond); if (info == 0) { if (calc_cond)
--- a/liboctave/numeric/CmplxCHOL.cc Thu Aug 20 14:37:57 2015 -0400 +++ b/liboctave/numeric/CmplxCHOL.cc Sun Aug 24 19:35:06 2014 +0530 @@ -86,7 +86,7 @@ } octave_idx_type -ComplexCHOL::init (const ComplexMatrix& a, bool calc_cond) +ComplexCHOL::init (const ComplexMatrix& a, bool upper, bool calc_cond) { octave_idx_type a_nr = a.rows (); octave_idx_type a_nc = a.cols (); @@ -101,13 +101,28 @@ octave_idx_type n = a_nc; octave_idx_type info; + is_upper = upper; + chol_mat.clear (n, n); - for (octave_idx_type j = 0; j < n; j++) + if (is_upper) { - 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; + 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 (); @@ -116,8 +131,18 @@ if (calc_cond) anorm = xnorm (a, 1); - F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1), n, h, n, info - F77_CHAR_ARG_LEN (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) @@ -143,7 +168,7 @@ } static ComplexMatrix -chol2inv_internal (const ComplexMatrix& r) +chol2inv_internal (const ComplexMatrix& r, bool is_upper = true) { ComplexMatrix retval; @@ -157,17 +182,37 @@ ComplexMatrix tmp = r; - F77_XFCN (zpotri, ZPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n, - tmp.fortran_vec (), n, info - F77_CHAR_ARG_LEN (1))); + 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) - 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) = std::conj (tmp.xelem (j, i)); + { + 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; } @@ -181,7 +226,7 @@ ComplexMatrix ComplexCHOL::inverse (void) const { - return chol2inv_internal (chol_mat); + return chol2inv_internal (chol_mat, is_upper); } void
--- a/liboctave/numeric/CmplxCHOL.h Thu Aug 20 14:37:57 2015 -0400 +++ b/liboctave/numeric/CmplxCHOL.h Sun Aug 24 19:35:06 2014 +0530 @@ -37,17 +37,17 @@ ComplexCHOL (void) : chol_mat (), xrcond (0) { } - ComplexCHOL (const ComplexMatrix& a, bool calc_cond = false) + ComplexCHOL (const ComplexMatrix& a, bool upper = true, bool calc_cond = false) : chol_mat (), xrcond (0) { - init (a, calc_cond); + init (a, upper, calc_cond); } - ComplexCHOL (const ComplexMatrix& a, octave_idx_type& info, + ComplexCHOL (const ComplexMatrix& a, octave_idx_type& info, bool upper = true, bool calc_cond = false) : chol_mat (), xrcond (0) { - info = init (a, calc_cond); + info = init (a, upper, calc_cond); } ComplexCHOL (const ComplexCHOL& a) @@ -91,7 +91,9 @@ double xrcond; - octave_idx_type init (const ComplexMatrix& a, bool calc_cond); + bool is_upper; + + octave_idx_type init (const ComplexMatrix& a, bool upper, bool calc_cond); }; ComplexMatrix OCTAVE_API chol2inv (const ComplexMatrix& r);
--- a/liboctave/numeric/dbleCHOL.cc Thu Aug 20 14:37:57 2015 -0400 +++ b/liboctave/numeric/dbleCHOL.cc Sun Aug 24 19:35:06 2014 +0530 @@ -87,7 +87,7 @@ } octave_idx_type -CHOL::init (const Matrix& a, bool calc_cond) +CHOL::init (const Matrix& a, bool upper, bool calc_cond) { octave_idx_type a_nr = a.rows (); octave_idx_type a_nc = a.cols (); @@ -101,13 +101,28 @@ octave_idx_type n = a_nc; octave_idx_type info; + is_upper = upper; + chol_mat.clear (n, n); - for (octave_idx_type j = 0; j < n; j++) + if (is_upper) { - 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; + 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 (); @@ -116,9 +131,18 @@ if (calc_cond) anorm = xnorm (a, 1); - F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1), - n, h, n, info - F77_CHAR_ARG_LEN (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) @@ -132,9 +156,19 @@ double *pz = z.fortran_vec (); Array<octave_idx_type> iz (dim_vector (n, 1)); octave_idx_type *piz = iz.fortran_vec (); - F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h, - n, anorm, xrcond, pz, piz, dpocon_info - F77_CHAR_ARG_LEN (1))); + 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; @@ -144,7 +178,7 @@ } static Matrix -chol2inv_internal (const Matrix& r) +chol2inv_internal (const Matrix& r, bool is_upper = true) { Matrix retval; @@ -161,17 +195,37 @@ if (info == 0) { - F77_XFCN (dpotri, DPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n, - v, n, info - F77_CHAR_ARG_LEN (1))); + 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) - 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); + { + 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; } @@ -186,7 +240,7 @@ Matrix CHOL::inverse (void) const { - return chol2inv_internal (chol_mat); + return chol2inv_internal (chol_mat, is_upper); } void
--- a/liboctave/numeric/dbleCHOL.h Thu Aug 20 14:37:57 2015 -0400 +++ b/liboctave/numeric/dbleCHOL.h Sun Aug 24 19:35:06 2014 +0530 @@ -37,16 +37,16 @@ CHOL (void) : chol_mat (), xrcond (0) { } - CHOL (const Matrix& a, bool calc_cond = false) + CHOL (const Matrix& a, bool upper = true, bool calc_cond = false) : chol_mat (), xrcond (0) { - init (a, calc_cond); + init (a, upper, calc_cond); } - CHOL (const Matrix& a, octave_idx_type& info, bool calc_cond = false) - : chol_mat (), xrcond (0) + CHOL (const Matrix& a, octave_idx_type& info, bool upper = true, + bool calc_cond = false) : chol_mat (), xrcond (0) { - info = init (a, calc_cond); + info = init (a, upper, calc_cond); } CHOL (const CHOL& a) : chol_mat (a.chol_mat), xrcond (a.xrcond) { } @@ -88,7 +88,9 @@ double xrcond; - octave_idx_type init (const Matrix& a, bool calc_cond); + bool is_upper; + + octave_idx_type init (const Matrix& a, bool upper, bool calc_cond); }; Matrix OCTAVE_API chol2inv (const Matrix& r);
--- a/liboctave/numeric/floatCHOL.cc Thu Aug 20 14:37:57 2015 -0400 +++ b/liboctave/numeric/floatCHOL.cc Sun Aug 24 19:35:06 2014 +0530 @@ -87,7 +87,7 @@ } octave_idx_type -FloatCHOL::init (const FloatMatrix& a, bool calc_cond) +FloatCHOL::init (const FloatMatrix& a, bool upper, bool calc_cond) { octave_idx_type a_nr = a.rows (); octave_idx_type a_nc = a.cols (); @@ -101,14 +101,30 @@ octave_idx_type n = a_nc; octave_idx_type info; + is_upper = upper; + chol_mat.clear (n, n); - for (octave_idx_type j = 0; j < n; j++) + if (is_upper) { - 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; + 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+1; 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. @@ -116,9 +132,18 @@ if (calc_cond) anorm = xnorm (a, 1); - F77_XFCN (spotrf, SPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1), - n, h, n, info - F77_CHAR_ARG_LEN (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) @@ -132,9 +157,19 @@ float *pz = z.fortran_vec (); Array<octave_idx_type> iz (dim_vector (n, 1)); octave_idx_type *piz = iz.fortran_vec (); - F77_XFCN (spocon, SPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h, - n, anorm, xrcond, pz, piz, spocon_info - F77_CHAR_ARG_LEN (1))); + 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; @@ -144,7 +179,7 @@ } static FloatMatrix -chol2inv_internal (const FloatMatrix& r) +chol2inv_internal (const FloatMatrix& r, bool is_upper = true) { FloatMatrix retval; @@ -161,17 +196,37 @@ if (info == 0) { - F77_XFCN (spotri, SPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n, - v, n, info - F77_CHAR_ARG_LEN (1))); + 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) - 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); + { + 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; } @@ -186,7 +241,7 @@ FloatMatrix FloatCHOL::inverse (void) const { - return chol2inv_internal (chol_mat); + return chol2inv_internal (chol_mat, is_upper); } void
--- a/liboctave/numeric/floatCHOL.h Thu Aug 20 14:37:57 2015 -0400 +++ b/liboctave/numeric/floatCHOL.h Sun Aug 24 19:35:06 2014 +0530 @@ -37,17 +37,17 @@ FloatCHOL (void) : chol_mat (), xrcond (0) { } - FloatCHOL (const FloatMatrix& a, bool calc_cond = false) + FloatCHOL (const FloatMatrix& a, bool upper = true, bool calc_cond = false) : chol_mat (), xrcond (0) { - init (a, calc_cond); + init (a, upper, calc_cond); } - FloatCHOL (const FloatMatrix& a, octave_idx_type& info, + FloatCHOL (const FloatMatrix& a, octave_idx_type& info, bool upper = true, bool calc_cond = false) : chol_mat (), xrcond (0) { - info = init (a, calc_cond); + info = init (a, upper, calc_cond); } FloatCHOL (const FloatCHOL& a) : chol_mat (a.chol_mat), xrcond (a.xrcond) { } @@ -90,7 +90,9 @@ float xrcond; - octave_idx_type init (const FloatMatrix& a, bool calc_cond); + bool is_upper; + + octave_idx_type init (const FloatMatrix& a, bool upper, bool calc_cond); }; FloatMatrix OCTAVE_API chol2inv (const FloatMatrix& r);