Mercurial > octave
changeset 32770:ed5e61fc7a4f bytecode-interpreter
maint: Merge default to bytecode-interpreter
| author | Arun Giridhar <arungiridhar@gmail.com> |
|---|---|
| date | Wed, 17 Jan 2024 14:28:58 -0500 |
| parents | adf0b68212e1 (current diff) d7ff09ba5303 (diff) |
| children | 80f7aec2c1c6 |
| files | |
| diffstat | 2 files changed, 116 insertions(+), 106 deletions(-) [+] |
line wrap: on
line diff
--- a/libinterp/corefcn/sqrtm.cc Tue Jan 16 20:38:58 2024 -0500 +++ b/libinterp/corefcn/sqrtm.cc Wed Jan 17 14:28:58 2024 -0500 @@ -40,51 +40,94 @@ OCTAVE_BEGIN_NAMESPACE(octave) -template <typename Matrix> +template <typename T> static void -sqrtm_utri_inplace (Matrix& T) +sqrtm_utri_inplace (T& m) { - typedef typename Matrix::element_type element_type; + typedef typename T::element_type element_type; + typedef typename T::real_matrix_type real_matrix_type; + typedef typename T::real_elt_type real_elt_type; const element_type zero = element_type (); bool singular = false; + bool diagonal = true; - // The following code is equivalent to this triple loop: - // - // n = rows (T); - // for j = 1:n - // T(j,j) = sqrt (T(j,j)); - // for i = j-1:-1:1 - // if T(i,j) != 0 - // T(i,j) /= (T(i,i) + T(j,j)); - // endif - // k = 1:i-1; - // T(k,j) -= T(k,i) * T(i,j); - // endfor - // endfor - // - // this is an in-place, cache-aligned variant of the code - // given in Higham's paper. + // The Schur matrix of Hermitian matrices is diagonal. + // check for off-diagonal elements above tolerance + const octave_idx_type n = m.rows (); + real_matrix_type abs_m = m.abs (); + + real_elt_type max_abs_diag = 0; + for (octave_idx_type i = 0; i < n; i++) + max_abs_diag = std::max (max_abs_diag, abs_m(i,i)); + + const real_elt_type tol = n * max_abs_diag + * std::numeric_limits<real_elt_type>::epsilon (); + + for (octave_idx_type j = 0; j < n; j++) + { + for (octave_idx_type i = j-1; i >= 0; i--) + { + if (abs_m(i,j) > tol) + { + diagonal = false; + break; + } + } + if (! diagonal) + break; + } - const octave_idx_type n = T.rows (); - element_type *Tp = T.rwdata (); - for (octave_idx_type j = 0; j < n; j++) + element_type *mp = m.fortran_vec (); + if (diagonal) + { + // shortcut for diagonal Schur matrices + for (octave_idx_type i = 0; i < n; i++) + { + octave_idx_type idx_diag = i*(n+1); + if (mp[idx_diag] != zero) + mp[idx_diag] = sqrt (mp[idx_diag]); + else + singular = true; + } + } + else { - element_type *colj = Tp + n*j; - if (colj[j] != zero) - colj[j] = sqrt (colj[j]); - else - singular = true; + // The following code is equivalent to this triple loop: + // + // n = rows (m); + // for j = 1:n + // m(j,j) = sqrt (m(j,j)); + // for i = j-1:-1:1 + // if m(i,j) != 0 + // m(i,j) /= (m(i,i) + m(j,j)); + // endif + // k = 1:i-1; + // m(k,j) -= m(k,i) * m(i,j); + // endfor + // endfor + // + // this is an in-place, cache-aligned variant of the code + // given in Higham's paper. - for (octave_idx_type i = j-1; i >= 0; i--) + for (octave_idx_type j = 0; j < n; j++) { - const element_type *coli = Tp + n*i; - if (colj[i] != zero) - colj[i] /= (coli[i] + colj[j]); - const element_type colji = colj[i]; - for (octave_idx_type k = 0; k < i; k++) - colj[k] -= coli[k] * colji; + element_type *colj = mp + n*j; + if (colj[j] != zero) + colj[j] = sqrt (colj[j]); + else + singular = true; + + for (octave_idx_type i = j-1; i >= 0; i--) + { + const element_type *coli = mp + n*i; + if (colj[i] != zero) + colj[i] /= (coli[i] + colj[j]); + const element_type colji = colj[i]; + for (octave_idx_type k = 0; k < i; k++) + colj[k] -= coli[k] * colji; + } } } @@ -186,11 +229,11 @@ x = schur_fact.schur_matrix (); u = schur_fact.unitary_schur_matrix (); } - while (0); // schur no longer needed. + while (0); // schur_fact no longer needed. sqrtm_utri_inplace (x); - x = u * x; // original x no longer needed. + x = u * x; // original x no longer needed. ComplexMatrix res = xgemm (x, u, blas_no_trans, blas_conj_trans); if (cutoff > 0 && xnorm (imag (res), one) <= cutoff)
--- a/libinterp/parse-tree/oct-parse.yy Tue Jan 16 20:38:58 2024 -0500 +++ b/libinterp/parse-tree/oct-parse.yy Wed Jan 17 14:28:58 2024 -0500 @@ -5498,34 +5498,6 @@ m_lexer.m_allow_command_syntax = false; } - // FIXME: this function partially duplicates do_dbtype in debug.cc. - static std::string - get_file_line (const std::string& name, int line) - { - // NAME should be an absolute file name and the file should exist. - - std::ifstream fs = sys::ifstream (name.c_str (), std::ios::in); - - std::string text; - - if (fs) - { - int i = 1; - - do - { - if (! std::getline (fs, text)) - { - text = ""; - break; - } - } - while (i++ < line); - } - - return text; - } - void base_parser::bison_error (const std::string& str) { @@ -5543,56 +5515,51 @@ { std::ostringstream output_buf; - if (m_lexer.m_reading_fcn_file || m_lexer.m_reading_script_file - || m_lexer.m_reading_classdef_file) - output_buf << "parse error near line " << err_line - << " of file " << m_lexer.m_fcn_file_full_name; - else - output_buf << "parse error:"; - - if (str != "parse error") - output_buf << "\n\n " << str; - - output_buf << "\n\n"; - - std::string curr_line; - - if (m_lexer.m_reading_fcn_file || m_lexer.m_reading_script_file - || m_lexer.m_reading_classdef_file) - curr_line = get_file_line (m_lexer.m_fcn_file_full_name, err_line); - else - curr_line = m_lexer.m_current_input_line; + bool in_file = (m_lexer.m_reading_fcn_file || m_lexer.m_reading_script_file + || m_lexer.m_reading_classdef_file); // Adjust the error column for display because it is 1-based in the // lexer for easier reporting. err_col--; - if (! curr_line.empty ()) + if (in_file) + { + output_buf << str + << " near line " << err_line << ", column " << err_col << "\n" + << "error: called from\n" + << " " << m_lexer.m_fcn_file_name + << " at line " << err_line << " column " << err_col << "\n"; + } + else { - // FIXME: we could do better if we just cached lines from the - // input file in a list. See also functions for managing input - // buffers in lex.ll. - - std::size_t len = curr_line.length (); - - if (curr_line[len-1] == '\n') - curr_line.resize (len-1); - - // Print the line, maybe with a pointer near the error token. - - output_buf << ">>> " << curr_line << "\n"; - - if (err_col == 0) - err_col = len; - - for (int i = 0; i < err_col + 3; i++) - output_buf << " "; - - output_buf << "^"; + // On command line, point directly to error + output_buf << str << "\n\n"; + std::string curr_line = m_lexer.m_current_input_line; + + if (! curr_line.empty ()) + { + // FIXME: we could do better if we just cached lines from the + // input file in a list. See also functions for managing input + // buffers in lex.ll. + std::size_t len = curr_line.length (); + + if (curr_line[len-1] == '\n') + curr_line.resize (len-1); + + // Print the line, maybe with a pointer near the error token. + output_buf << ">>> " << curr_line << "\n"; + + if (err_col == 0) + err_col = len; + + for (int i = 0; i < err_col + 3; i++) + output_buf << " "; + + output_buf << "^" << "\n"; + } + } - output_buf << "\n"; - m_parse_error_msg = output_buf.str (); }