diff src/DLD-FUNCTIONS/chol.cc @ 8547:d66c9b6e506a

imported patch qrupdate.diff
author Jaroslav Hajek <highegg@gmail.com>
date Tue, 20 Jan 2009 21:16:42 +0100
parents 81d6ab3ac93c
children a6edd5c23cb5
line wrap: on
line diff
--- a/src/DLD-FUNCTIONS/chol.cc	Tue Jan 20 21:15:17 2009 +0100
+++ b/src/DLD-FUNCTIONS/chol.cc	Tue Jan 20 21:16:42 2009 +0100
@@ -2,6 +2,8 @@
 
 Copyright (C) 1996, 1997, 1999, 2000, 2002, 2005, 2006, 2007
               John W. Eaton
+Copyright (C) 2008, 2009 Jaroslav Hajek
+Copyright (C) 2008, 2009 VZLU Prague
 
 This file is part of Octave.
 
@@ -21,9 +23,6 @@
 
 */
 
-// The cholupdate, cholinsert, choldelete and cholshift functions were
-//  written by Jaroslav Hajek <highegg@gmail.com>, Copyright (C) 2008  
-//  VZLU Prague, a.s., Czech Republic.
 
 #ifdef HAVE_CONFIG_H
 #include <config.h>
@@ -577,6 +576,8 @@
   return retval;
 }
 
+#ifdef HAVE_QRUPDATE
+
 DEFUN_DLD (cholupdate, args, nargout,
   "-*- texinfo -*-\n\
 @deftypefn {Loadable Function} {[@var{R1}, @var{info}] =} cholupdate (@var{R}, @var{u}, @var{op})\n\
@@ -628,100 +629,84 @@
       if (down || op == "+")
         if (argr.columns () == n && argu.rows () == n && argu.columns () == 1)
           {
+            int err = 0;
 	    if (argr.is_single_type () || argu.is_single_type ())
 	      {
-		if (argr.is_real_matrix () && argu.is_real_matrix ())
+		if (argr.is_real_type () && argu.is_real_type ())
 		  {
 		    // real case
 		    FloatMatrix R = argr.float_matrix_value ();
-		    FloatMatrix u = argu.float_matrix_value ();
+		    FloatColumnVector u = argu.float_column_vector_value ();
 
 		    FloatCHOL fact;
 		    fact.set (R);
-		    int err = 0;
 
 		    if (down)
 		      err = fact.downdate (u);
 		    else
 		      fact.update (u);
 
-		    if (nargout > 1)
-		      retval(1) = err;
-		    else if (err)
-		      error ("cholupdate: downdate violates positiveness");
-
 		    retval(0) = fact.chol_matrix ();
 		  }
 		else
 		  {
 		    // complex case
 		    FloatComplexMatrix R = argr.float_complex_matrix_value ();
-		    FloatComplexMatrix u = argu.float_complex_matrix_value ();
+		    FloatComplexColumnVector u = argu.float_complex_column_vector_value ();
 
 		    FloatComplexCHOL fact;
 		    fact.set (R);
-		    int err = 0;
 
 		    if (down)
 		      err = fact.downdate (u);
 		    else
 		      fact.update (u);
 
-		    if (nargout > 1)
-		      retval(1) = err;
-		    else if (err)
-		      error ("cholupdate: downdate violates positiveness");
-
 		    retval(0) = fact.chol_matrix ();
 		  }
 	      }
 	    else
 	      {
-		if (argr.is_real_matrix () && argu.is_real_matrix ())
+		if (argr.is_real_type () && argu.is_real_type ())
 		  {
 		    // real case
 		    Matrix R = argr.matrix_value ();
-		    Matrix u = argu.matrix_value ();
+		    ColumnVector u = argu.column_vector_value ();
 
 		    CHOL fact;
 		    fact.set (R);
-		    int err = 0;
 
 		    if (down)
 		      err = fact.downdate (u);
 		    else
 		      fact.update (u);
 
-		    if (nargout > 1)
-		      retval(1) = err;
-		    else if (err)
-		      error ("cholupdate: downdate violates positiveness");
-
 		    retval(0) = fact.chol_matrix ();
 		  }
 		else
 		  {
 		    // complex case
 		    ComplexMatrix R = argr.complex_matrix_value ();
-		    ComplexMatrix u = argu.complex_matrix_value ();
+		    ComplexColumnVector u = argu.complex_column_vector_value ();
 
 		    ComplexCHOL fact;
 		    fact.set (R);
-		    int err = 0;
 
 		    if (down)
 		      err = fact.downdate (u);
 		    else
 		      fact.update (u);
 
-		    if (nargout > 1)
-		      retval(1) = err;
-		    else if (err)
-		      error ("cholupdate: downdate violates positiveness");
-
 		    retval(0) = fact.chol_matrix ();
 		  }
 	      }
+
+            if (nargout > 1)
+              retval(1) = err;
+            else if (err == 1)
+              error ("cholupdate: downdate violates positiveness");
+            else if (err == 2)
+              error ("cholupdate: singular matrix");
           }
         else
           error ("cholupdate: dimension mismatch");
@@ -853,22 +838,18 @@
         {
           if (j > 0 && j <= n+1)
             {
+              int err = 0;
 	      if (argr.is_single_type () || argu.is_single_type ())
 		{
-		  if (argr.is_real_matrix () && argu.is_real_matrix ())
+		  if (argr.is_real_type () && argu.is_real_type ())
 		    {
 		      // real case
 		      FloatMatrix R = argr.float_matrix_value ();
-		      FloatMatrix u = argu.float_matrix_value ();
+		      FloatColumnVector u = argu.float_column_vector_value ();
 
 		      FloatCHOL fact;
 		      fact.set (R);
-		      int err = fact.insert_sym (u, j-1);
-
-		      if (nargout > 1)
-			retval(1) = err;
-		      else if (err)
-			error ("cholinsert: insertion violates positiveness");
+		      err = fact.insert_sym (u, j-1);
 
 		      retval(0) = fact.chol_matrix ();
 		    }
@@ -876,36 +857,26 @@
 		    {
 		      // complex case
 		      FloatComplexMatrix R = argr.float_complex_matrix_value ();
-		      FloatComplexMatrix u = argu.float_complex_matrix_value ();
+		      FloatComplexColumnVector u = argu.float_complex_column_vector_value ();
 
 		      FloatComplexCHOL fact;
 		      fact.set (R);
-		      int err = fact.insert_sym (u, j-1);
-
-		      if (nargout > 1)
-			retval(1) = err;
-		      else if (err)
-			error ("cholinsert: insertion violates positiveness");
+		      err = fact.insert_sym (u, j-1);
 
 		      retval(0) = fact.chol_matrix ();
 		    }
 		}
 	      else
 		{
-		  if (argr.is_real_matrix () && argu.is_real_matrix ())
+		  if (argr.is_real_type () && argu.is_real_type ())
 		    {
 		      // real case
 		      Matrix R = argr.matrix_value ();
-		      Matrix u = argu.matrix_value ();
+		      ColumnVector u = argu.column_vector_value ();
 
 		      CHOL fact;
 		      fact.set (R);
-		      int err = fact.insert_sym (u, j-1);
-
-		      if (nargout > 1)
-			retval(1) = err;
-		      else if (err)
-			error ("cholinsert: insertion violates positiveness");
+		      err = fact.insert_sym (u, j-1);
 
 		      retval(0) = fact.chol_matrix ();
 		    }
@@ -913,20 +884,24 @@
 		    {
 		      // complex case
 		      ComplexMatrix R = argr.complex_matrix_value ();
-		      ComplexMatrix u = argu.complex_matrix_value ();
+		      ComplexColumnVector u = argu.complex_column_vector_value ();
 
 		      ComplexCHOL fact;
 		      fact.set (R);
-		      int err = fact.insert_sym (u, j-1);
-
-		      if (nargout > 1)
-			retval(1) = err;
-		      else if (err)
-			error ("cholinsert: insertion violates positiveness");
+		      err = fact.insert_sym (u, j-1);
 
 		      retval(0) = fact.chol_matrix ();
 		    }
 		}
+
+              if (nargout > 1)
+                retval(1) = err;
+              else if (err == 1)
+                error ("cholinsert: insertion violates positiveness");
+              else if (err == 2)
+                error ("cholinsert: singular matrix");
+              else if (err == 3)
+                error ("cholinsert: diagonal element must be real");
             }
           else
             error ("cholinsert: index out of range");
@@ -1037,7 +1012,7 @@
             {
 	      if (argr.is_single_type ())
 		{
-		  if (argr.is_real_matrix ())
+		  if (argr.is_real_type ())
 		    {
 		      // real case
 		      FloatMatrix R = argr.float_matrix_value ();
@@ -1062,7 +1037,7 @@
 		}
 	      else
 		{
-		  if (argr.is_real_matrix ())
+		  if (argr.is_real_type ())
 		    {
 		      // real case
 		      Matrix R = argr.matrix_value ();
@@ -1178,7 +1153,7 @@
 	      if (argr.is_single_type () && argi.is_single_type () && 
 		  argj.is_single_type ())
 		{
-		  if (argr.is_real_matrix ())
+		  if (argr.is_real_type ())
 		    {
 		      // real case
 		      FloatMatrix R = argr.float_matrix_value ();
@@ -1203,7 +1178,7 @@
 		}
 	      else
 		{
-		  if (argr.is_real_matrix ())
+		  if (argr.is_real_type ())
 		    {
 		      // real case
 		      Matrix R = argr.matrix_value ();
@@ -1301,6 +1276,8 @@
 %! assert(norm(R1'*R1 - single(Ac(p,p)),Inf) < 1e1*eps('single'))
 */
 
+#endif
+
 /*
 ;;; Local Variables: ***
 ;;; mode: C++ ***