Mercurial > octave
changeset 20626:9fb8133288e8
chol: Support chol(a, "lower") for complex float types (bug #34850)
* fCmplxCHOL.cc, fCmplxCHOL.h (FloatComplexCHOL::init, chol2inv_internal):
Add is_upper parameter, propagate to underlying library functions.
* fCMatrix.cc (FloatComplexMatrix::inverse): Use new argument.
* chol.cc (Fchol): Likewise.
This adds float complex support that was omitted in cset 5ce959c55cc0.
| author | Mike Miller <mtmiller@octave.org> |
|---|---|
| date | Wed, 14 Oct 2015 21:31:22 -0400 |
| parents | 974b218e7292 |
| children | ec003232f7aa |
| files | libinterp/dldfcn/chol.cc liboctave/array/fCMatrix.cc liboctave/numeric/fCmplxCHOL.cc liboctave/numeric/fCmplxCHOL.h |
| diffstat | 4 files changed, 54 insertions(+), 32 deletions(-) [+] |
line wrap: on
line diff
--- a/libinterp/dldfcn/chol.cc Sat Sep 26 19:25:03 2015 +0200 +++ b/libinterp/dldfcn/chol.cc Wed Oct 14 21:31:22 2015 -0400 @@ -271,18 +271,12 @@ octave_idx_type info; FloatComplexCHOL fact; - if (LLt) - fact = FloatComplexCHOL (m.transpose (), info); - else - fact = FloatComplexCHOL (m, info); + fact = FloatComplexCHOL (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/fCMatrix.cc Sat Sep 26 19:25:03 2015 +0200 +++ b/liboctave/array/fCMatrix.cc Wed Oct 14 21:31:22 2015 -0400 @@ -1199,7 +1199,7 @@ { if (mattype.is_hermitian ()) { - FloatComplexCHOL chol (*this, info, calc_cond); + FloatComplexCHOL chol (*this, info, true, calc_cond); if (info == 0) { if (calc_cond)
--- a/liboctave/numeric/fCmplxCHOL.cc Sat Sep 26 19:25:03 2015 +0200 +++ b/liboctave/numeric/fCmplxCHOL.cc Wed Oct 14 21:31:22 2015 -0400 @@ -86,7 +86,7 @@ } octave_idx_type -FloatComplexCHOL::init (const FloatComplexMatrix& a, bool calc_cond) +FloatComplexCHOL::init (const FloatComplexMatrix& a, bool upper, bool calc_cond) { octave_idx_type a_nr = a.rows (); octave_idx_type a_nc = a.cols (); @@ -101,14 +101,25 @@ 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++) - { - 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; - } + 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. @@ -116,8 +127,12 @@ if (calc_cond) anorm = xnorm (a, 1); - F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1), n, h, n, info - F77_CHAR_ARG_LEN (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) @@ -143,7 +158,7 @@ } static FloatComplexMatrix -chol2inv_internal (const FloatComplexMatrix& r) +chol2inv_internal (const FloatComplexMatrix& r, bool is_upper = true) { FloatComplexMatrix retval; @@ -157,17 +172,27 @@ FloatComplexMatrix tmp = r; - F77_XFCN (cpotri, CPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n, - tmp.fortran_vec (), n, info - F77_CHAR_ARG_LEN (1))); + 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) - 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 +206,7 @@ FloatComplexMatrix FloatComplexCHOL::inverse (void) const { - return chol2inv_internal (chol_mat); + return chol2inv_internal (chol_mat, is_upper); } void
--- a/liboctave/numeric/fCmplxCHOL.h Sat Sep 26 19:25:03 2015 +0200 +++ b/liboctave/numeric/fCmplxCHOL.h Wed Oct 14 21:31:22 2015 -0400 @@ -37,17 +37,18 @@ FloatComplexCHOL (void) : chol_mat (), xrcond (0) { } - FloatComplexCHOL (const FloatComplexMatrix& a, bool calc_cond = false) + FloatComplexCHOL (const FloatComplexMatrix& a, bool upper = true, + bool calc_cond = false) : chol_mat (), xrcond (0) { - init (a, calc_cond); + init (a, upper, calc_cond); } FloatComplexCHOL (const FloatComplexMatrix& a, octave_idx_type& info, - bool calc_cond = false) + bool upper = true, bool calc_cond = false) : chol_mat (), xrcond (0) { - info = init (a, calc_cond); + info = init (a, upper, calc_cond); } FloatComplexCHOL (const FloatComplexCHOL& a) @@ -92,7 +93,9 @@ float xrcond; - octave_idx_type init (const FloatComplexMatrix& a, bool calc_cond); + bool is_upper; + + octave_idx_type init (const FloatComplexMatrix& a, bool upper, bool calc_cond); }; FloatComplexMatrix OCTAVE_API chol2inv (const FloatComplexMatrix& r);