# HG changeset patch # User John W. Eaton # Date 1311973392 14400 # Node ID b9c66f8912819d632f96c6f90ab2646ed93fe908 # Parent f754b65f4bc569fe39ea0f7f0570211cdf090f7a maint: style fixes for sqrtm * sqrtm.cc: Style fixes. diff -r f754b65f4bc5 -r b9c66f891281 src/DLD-FUNCTIONS/sqrtm.cc --- a/src/DLD-FUNCTIONS/sqrtm.cc Fri Jul 29 11:32:38 2011 -0500 +++ b/src/DLD-FUNCTIONS/sqrtm.cc Fri Jul 29 17:03:12 2011 -0400 @@ -49,22 +49,20 @@ bool singular = false; - /* - * 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 - * T(i,j) /= (T(i,i) + T(j,j)); - * 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 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 + // T(i,j) /= (T(i,i) + T(j,j)); + // 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. const octave_idx_type n = T.rows (); element_type *Tp = T.fortran_vec (); @@ -117,38 +115,32 @@ { case MatrixType::Upper: case MatrixType::Diagonal: - { - if (! x.diag ().any_element_is_negative ()) - { - // Do it in real arithmetic. - sqrtm_utri_inplace (x); - retval = x; - retval.matrix_type (mt); - } - else - iscomplex = true; - - break; - } - case MatrixType::Lower: + if (! x.diag ().any_element_is_negative ()) { - if (! x.diag ().any_element_is_negative ()) - { - x = x.transpose (); - sqrtm_utri_inplace (x); - retval = x.transpose (); - retval.matrix_type (mt); - } - else - iscomplex = true; + // Do it in real arithmetic. + sqrtm_utri_inplace (x); + retval = x; + retval.matrix_type (mt); + } + else + iscomplex = true; + break; - break; + case MatrixType::Lower: + if (! x.diag ().any_element_is_negative ()) + { + x = x.transpose (); + sqrtm_utri_inplace (x); + retval = x.transpose (); + retval.matrix_type (mt); } + else + iscomplex = true; + break; + default: - { - iscomplex = true; - break; - } + iscomplex = true; + break; } if (iscomplex) @@ -166,46 +158,41 @@ { case MatrixType::Upper: case MatrixType::Diagonal: - { - sqrtm_utri_inplace (x); - retval = x; - retval.matrix_type (mt); + sqrtm_utri_inplace (x); + retval = x; + retval.matrix_type (mt); + break; - break; - } case MatrixType::Lower: - { - x = x.transpose (); - sqrtm_utri_inplace (x); - retval = x.transpose (); - retval.matrix_type (mt); + x = x.transpose (); + sqrtm_utri_inplace (x); + retval = x.transpose (); + retval.matrix_type (mt); + break; - break; - } default: - { - ComplexMatrix u; + { + ComplexMatrix u; - do - { - ComplexSCHUR schur (x, std::string (), true); - x = schur.schur_matrix (); - u = schur.unitary_matrix (); - } - while (0); // schur no longer needed. + do + { + ComplexSCHUR schur (x, std::string (), true); + x = schur.schur_matrix (); + u = schur.unitary_matrix (); + } + while (0); // schur no longer needed. - sqrtm_utri_inplace (x); + sqrtm_utri_inplace (x); - x = u * x; // original x no longer needed. - ComplexMatrix res = xgemm (x, u, blas_no_trans, blas_conj_trans); + 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) - retval = real (res); - else - retval = res; - - break; - } + if (cutoff > 0 && xnorm (imag (res), one) <= cutoff) + retval = real (res); + else + retval = res; + } + break; } } @@ -246,22 +233,17 @@ } if (arg.is_diag_matrix ()) - { - // sqrtm of a diagonal matrix is just sqrt. - retval(0) = arg.sqrt (); - } + // sqrtm of a diagonal matrix is just sqrt. + retval(0) = arg.sqrt (); else if (arg.is_single_type ()) - { - retval(0) = do_sqrtm (arg); - } + retval(0) = do_sqrtm (arg); else if (arg.is_numeric_type ()) - { - retval(0) = do_sqrtm (arg); - } + retval(0) = do_sqrtm (arg); if (nargout > 1 && ! error_state) { // This corresponds to generic code + // // norm (s*s - x, "fro") / norm (x, "fro"); octave_value s = retval(0);