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;