Mercurial > octave
diff src/DLD-FUNCTIONS/sqrtm.cc @ 10154:40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Wed, 20 Jan 2010 17:33:41 -0500 |
parents | 7c02ec148a3c |
children | d0ce5e973937 |
line wrap: on
line diff
--- a/src/DLD-FUNCTIONS/sqrtm.cc Wed Jan 20 17:24:23 2010 -0500 +++ b/src/DLD-FUNCTIONS/sqrtm.cc Wed Jan 20 17:33:41 2010 -0500 @@ -111,23 +111,23 @@ for (octave_idx_type p = 0; p < n-1; p++) { for (octave_idx_type i = 0; i < n-(p+1); i++) - { - const octave_idx_type j = i + p + 1; + { + const octave_idx_type j = i + p + 1; - Complex s = T(i,j); + Complex s = T(i,j); - for (octave_idx_type k = i+1; k < j; k++) - s -= R(i,k) * R(k,j); + for (octave_idx_type k = i+1; k < j; k++) + s -= R(i,k) * R(k,j); - // dividing - // R(i,j) = s/(R(i,i)+R(j,j)); - // screwing around to not / 0 + // dividing + // R(i,j) = s/(R(i,i)+R(j,j)); + // screwing around to not / 0 - const Complex d = R(i,i) + R(j,j) + fudge; - const Complex conjd = conj (d); + const Complex d = R(i,i) + R(j,j) + fudge; + const Complex conjd = conj (d); - R(i,j) = (s*conjd)/(d*conjd); - } + R(i,j) = (s*conjd)/(d*conjd); + } } return U * R * U.hermitian (); @@ -148,23 +148,23 @@ for (octave_idx_type p = 0; p < n-1; p++) { for (octave_idx_type i = 0; i < n-(p+1); i++) - { - const octave_idx_type j = i + p + 1; + { + const octave_idx_type j = i + p + 1; - FloatComplex s = T(i,j); + FloatComplex s = T(i,j); - for (octave_idx_type k = i+1; k < j; k++) - s -= R(i,k) * R(k,j); + for (octave_idx_type k = i+1; k < j; k++) + s -= R(i,k) * R(k,j); - // dividing - // R(i,j) = s/(R(i,i)+R(j,j)); - // screwing around to not / 0 + // dividing + // R(i,j) = s/(R(i,i)+R(j,j)); + // screwing around to not / 0 - const FloatComplex d = R(i,i) + R(j,j) + fudge; - const FloatComplex conjd = conj (d); + const FloatComplex d = R(i,i) + R(j,j) + fudge; + const FloatComplex conjd = conj (d); - R(i,j) = (s*conjd)/(d*conjd); - } + R(i,j) = (s*conjd)/(d*conjd); + } } return U * R * U.hermitian (); @@ -216,244 +216,244 @@ if (arg.is_single_type ()) { if (arg.is_real_scalar ()) - { - float d = arg.float_value (); - if (d > 0.0) - { - retval(0) = sqrt (d); - retval(1) = 0.0; - } - else - { - retval(0) = FloatComplex (0.0, sqrt (d)); - retval(1) = 0.0; - } - } + { + float d = arg.float_value (); + if (d > 0.0) + { + retval(0) = sqrt (d); + retval(1) = 0.0; + } + else + { + retval(0) = FloatComplex (0.0, sqrt (d)); + retval(1) = 0.0; + } + } else if (arg.is_complex_scalar ()) - { - FloatComplex c = arg.float_complex_value (); - retval(0) = sqrt (c); - retval(1) = 0.0; - } + { + FloatComplex c = arg.float_complex_value (); + retval(0) = sqrt (c); + retval(1) = 0.0; + } else if (arg.is_matrix_type ()) - { - float err, minT; + { + float err, minT; - if (arg.is_real_matrix ()) - { - FloatMatrix A = arg.float_matrix_value(); + if (arg.is_real_matrix ()) + { + FloatMatrix A = arg.float_matrix_value(); - if (error_state) - return retval; + if (error_state) + return retval; - // FIXME -- eventually, FloatComplexSCHUR will accept a - // real matrix arg. + // FIXME -- eventually, FloatComplexSCHUR will accept a + // real matrix arg. - FloatComplexMatrix Ac (A); + FloatComplexMatrix Ac (A); - const FloatComplexSCHUR schur (Ac, std::string ()); + const FloatComplexSCHUR schur (Ac, std::string ()); - if (error_state) - return retval; + if (error_state) + return retval; - const FloatComplexMatrix U (schur.unitary_matrix ()); - const FloatComplexMatrix T (schur.schur_matrix ()); - const FloatComplexMatrix X (sqrtm_from_schur (U, T)); + const FloatComplexMatrix U (schur.unitary_matrix ()); + const FloatComplexMatrix T (schur.schur_matrix ()); + const FloatComplexMatrix X (sqrtm_from_schur (U, T)); - // Check for minimal imaginary part - float normX = 0.0; - float imagX = 0.0; - for (octave_idx_type i = 0; i < n; i++) - for (octave_idx_type j = 0; j < n; j++) - { - imagX = getmax (imagX, imag (X(i,j))); - normX = getmax (normX, abs (X(i,j))); - } + // Check for minimal imaginary part + float normX = 0.0; + float imagX = 0.0; + for (octave_idx_type i = 0; i < n; i++) + for (octave_idx_type j = 0; j < n; j++) + { + imagX = getmax (imagX, imag (X(i,j))); + normX = getmax (normX, abs (X(i,j))); + } - if (imagX < normX * 100 * FLT_EPSILON) - retval(0) = real (X); - else - retval(0) = X; + if (imagX < normX * 100 * FLT_EPSILON) + retval(0) = real (X); + else + retval(0) = X; - // Compute error - // FIXME can we estimate the error without doing the - // matrix multiply? + // Compute error + // FIXME can we estimate the error without doing the + // matrix multiply? - err = frobnorm (X*X - FloatComplexMatrix (A)) / frobnorm (A); + err = frobnorm (X*X - FloatComplexMatrix (A)) / frobnorm (A); - if (xisnan (err)) - err = lo_ieee_float_inf_value (); + if (xisnan (err)) + err = lo_ieee_float_inf_value (); - // Find min diagonal - minT = lo_ieee_float_inf_value (); - for (octave_idx_type i=0; i < n; i++) - minT = getmin(minT, abs(T(i,i))); - } - else - { - FloatComplexMatrix A = arg.float_complex_matrix_value (); + // Find min diagonal + minT = lo_ieee_float_inf_value (); + for (octave_idx_type i=0; i < n; i++) + minT = getmin(minT, abs(T(i,i))); + } + else + { + FloatComplexMatrix A = arg.float_complex_matrix_value (); - if (error_state) - return retval; + if (error_state) + return retval; - const FloatComplexSCHUR schur (A, std::string ()); + const FloatComplexSCHUR schur (A, std::string ()); - if (error_state) - return retval; + if (error_state) + return retval; - const FloatComplexMatrix U (schur.unitary_matrix ()); - const FloatComplexMatrix T (schur.schur_matrix ()); - const FloatComplexMatrix X (sqrtm_from_schur (U, T)); + const FloatComplexMatrix U (schur.unitary_matrix ()); + const FloatComplexMatrix T (schur.schur_matrix ()); + const FloatComplexMatrix X (sqrtm_from_schur (U, T)); - retval(0) = X; + retval(0) = X; - err = frobnorm (X*X - A) / frobnorm (A); + err = frobnorm (X*X - A) / frobnorm (A); - if (xisnan (err)) - err = lo_ieee_float_inf_value (); + if (xisnan (err)) + err = lo_ieee_float_inf_value (); - minT = lo_ieee_float_inf_value (); - for (octave_idx_type i = 0; i < n; i++) - minT = getmin (minT, abs (T(i,i))); - } + minT = lo_ieee_float_inf_value (); + for (octave_idx_type i = 0; i < n; i++) + minT = getmin (minT, abs (T(i,i))); + } - retval(1) = err; + retval(1) = err; - if (nargout < 2) - { - if (err > 100*(minT+FLT_EPSILON)*n) - { - if (minT == 0.0) - error ("sqrtm: A is singular, sqrt may not exist"); - else if (minT <= sqrt (FLT_MIN)) - error ("sqrtm: A is nearly singular, failed to find sqrt"); - else - error ("sqrtm: failed to find sqrt"); - } - } - } + if (nargout < 2) + { + if (err > 100*(minT+FLT_EPSILON)*n) + { + if (minT == 0.0) + error ("sqrtm: A is singular, sqrt may not exist"); + else if (minT <= sqrt (FLT_MIN)) + error ("sqrtm: A is nearly singular, failed to find sqrt"); + else + error ("sqrtm: failed to find sqrt"); + } + } + } } else { if (arg.is_real_scalar ()) - { - double d = arg.double_value (); - if (d > 0.0) - { - retval(0) = sqrt (d); - retval(1) = 0.0; - } - else - { - retval(0) = Complex (0.0, sqrt (d)); - retval(1) = 0.0; - } - } + { + double d = arg.double_value (); + if (d > 0.0) + { + retval(0) = sqrt (d); + retval(1) = 0.0; + } + else + { + retval(0) = Complex (0.0, sqrt (d)); + retval(1) = 0.0; + } + } else if (arg.is_complex_scalar ()) - { - Complex c = arg.complex_value (); - retval(0) = sqrt (c); - retval(1) = 0.0; - } + { + Complex c = arg.complex_value (); + retval(0) = sqrt (c); + retval(1) = 0.0; + } else if (arg.is_matrix_type ()) - { - double err, minT; + { + double err, minT; - if (arg.is_real_matrix ()) - { - Matrix A = arg.matrix_value(); + if (arg.is_real_matrix ()) + { + Matrix A = arg.matrix_value(); - if (error_state) - return retval; + if (error_state) + return retval; - // FIXME -- eventually, ComplexSCHUR will accept a - // real matrix arg. + // FIXME -- eventually, ComplexSCHUR will accept a + // real matrix arg. - ComplexMatrix Ac (A); + ComplexMatrix Ac (A); - const ComplexSCHUR schur (Ac, std::string ()); + const ComplexSCHUR schur (Ac, std::string ()); - if (error_state) - return retval; + if (error_state) + return retval; - const ComplexMatrix U (schur.unitary_matrix ()); - const ComplexMatrix T (schur.schur_matrix ()); - const ComplexMatrix X (sqrtm_from_schur (U, T)); + const ComplexMatrix U (schur.unitary_matrix ()); + const ComplexMatrix T (schur.schur_matrix ()); + const ComplexMatrix X (sqrtm_from_schur (U, T)); - // Check for minimal imaginary part - double normX = 0.0; - double imagX = 0.0; - for (octave_idx_type i = 0; i < n; i++) - for (octave_idx_type j = 0; j < n; j++) - { - imagX = getmax (imagX, imag (X(i,j))); - normX = getmax (normX, abs (X(i,j))); - } + // Check for minimal imaginary part + double normX = 0.0; + double imagX = 0.0; + for (octave_idx_type i = 0; i < n; i++) + for (octave_idx_type j = 0; j < n; j++) + { + imagX = getmax (imagX, imag (X(i,j))); + normX = getmax (normX, abs (X(i,j))); + } - if (imagX < normX * 100 * DBL_EPSILON) - retval(0) = real (X); - else - retval(0) = X; + if (imagX < normX * 100 * DBL_EPSILON) + retval(0) = real (X); + else + retval(0) = X; - // Compute error - // FIXME can we estimate the error without doing the - // matrix multiply? + // Compute error + // FIXME can we estimate the error without doing the + // matrix multiply? - err = frobnorm (X*X - ComplexMatrix (A)) / frobnorm (A); + err = frobnorm (X*X - ComplexMatrix (A)) / frobnorm (A); - if (xisnan (err)) - err = lo_ieee_inf_value (); + if (xisnan (err)) + err = lo_ieee_inf_value (); - // Find min diagonal - minT = lo_ieee_inf_value (); - for (octave_idx_type i=0; i < n; i++) - minT = getmin(minT, abs(T(i,i))); - } - else - { - ComplexMatrix A = arg.complex_matrix_value (); + // Find min diagonal + minT = lo_ieee_inf_value (); + for (octave_idx_type i=0; i < n; i++) + minT = getmin(minT, abs(T(i,i))); + } + else + { + ComplexMatrix A = arg.complex_matrix_value (); - if (error_state) - return retval; + if (error_state) + return retval; - const ComplexSCHUR schur (A, std::string ()); + const ComplexSCHUR schur (A, std::string ()); - if (error_state) - return retval; + if (error_state) + return retval; - const ComplexMatrix U (schur.unitary_matrix ()); - const ComplexMatrix T (schur.schur_matrix ()); - const ComplexMatrix X (sqrtm_from_schur (U, T)); + const ComplexMatrix U (schur.unitary_matrix ()); + const ComplexMatrix T (schur.schur_matrix ()); + const ComplexMatrix X (sqrtm_from_schur (U, T)); - retval(0) = X; + retval(0) = X; - err = frobnorm (X*X - A) / frobnorm (A); + err = frobnorm (X*X - A) / frobnorm (A); - if (xisnan (err)) - err = lo_ieee_inf_value (); + if (xisnan (err)) + err = lo_ieee_inf_value (); - minT = lo_ieee_inf_value (); - for (octave_idx_type i = 0; i < n; i++) - minT = getmin (minT, abs (T(i,i))); - } + minT = lo_ieee_inf_value (); + for (octave_idx_type i = 0; i < n; i++) + minT = getmin (minT, abs (T(i,i))); + } - retval(1) = err; + retval(1) = err; - if (nargout < 2) - { - if (err > 100*(minT+DBL_EPSILON)*n) - { - if (minT == 0.0) - error ("sqrtm: A is singular, sqrt may not exist"); - else if (minT <= sqrt (DBL_MIN)) - error ("sqrtm: A is nearly singular, failed to find sqrt"); - else - error ("sqrtm: failed to find sqrt"); - } - } - } + if (nargout < 2) + { + if (err > 100*(minT+DBL_EPSILON)*n) + { + if (minT == 0.0) + error ("sqrtm: A is singular, sqrt may not exist"); + else if (minT <= sqrt (DBL_MIN)) + error ("sqrtm: A is nearly singular, failed to find sqrt"); + else + error ("sqrtm: failed to find sqrt"); + } + } + } else - gripe_wrong_type_arg ("sqrtm", arg); + gripe_wrong_type_arg ("sqrtm", arg); } return retval;