diff liboctave/CmplxCHOL.cc @ 8562:a6edd5c23cb5

use replacement methods if qrupdate is not available
author Jaroslav Hajek <highegg@gmail.com>
date Thu, 22 Jan 2009 11:10:47 +0100
parents d66c9b6e506a
children c0aeedd8fb86
line wrap: on
line diff
--- a/liboctave/CmplxCHOL.cc	Wed Jan 21 21:48:24 2009 -0500
+++ b/liboctave/CmplxCHOL.cc	Thu Jan 22 11:10:47 2009 +0100
@@ -34,6 +34,9 @@
 #include "f77-fcn.h"
 #include "lo-error.h"
 #include "oct-locbuf.h"
+#ifndef HAVE_QRUPDATE
+#include "dbleQR.h"
+#endif
 
 extern "C"
 {
@@ -291,6 +294,146 @@
     }
 }
 
+#else
+
+void
+ComplexCHOL::update (const ComplexColumnVector& u)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type n = chol_mat.rows ();
+
+  if (u.length () == n)
+    {
+      init (chol_mat.hermitian () * chol_mat 
+            + ComplexMatrix (u) * ComplexMatrix (u).hermitian (), false);
+    }
+  else
+    (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
+}
+
+static bool
+singular (const ComplexMatrix& a)
+{
+  for (octave_idx_type i = 0; i < a.rows (); i++)
+    if (a(i,i) == 0.0) return true;
+  return false;
+}
+
+octave_idx_type
+ComplexCHOL::downdate (const ComplexColumnVector& u)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type info = -1;
+
+  octave_idx_type n = chol_mat.rows ();
+
+  if (u.length () == n)
+    {
+      if (singular (chol_mat))
+        info = 2;
+      else
+        {
+          info = init (chol_mat.hermitian () * chol_mat 
+                       - ComplexMatrix (u) * ComplexMatrix (u).hermitian (), false);
+          if (info) info = 1;
+        }
+    }
+  else
+    (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
+
+  return info;
+}
+
+octave_idx_type
+ComplexCHOL::insert_sym (const ComplexColumnVector& u, octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type info = -1;
+
+  octave_idx_type n = chol_mat.rows ();
+  
+  if (u.length () != n + 1)
+    (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
+  else if (j < 0 || j > n)
+    (*current_liboctave_error_handler) ("cholinsert: index out of range");
+  else
+    {
+      if (singular (chol_mat))
+        info = 2;
+      else if (u(j).imag () != 0.0)
+        info = 3;
+      else
+        {
+          ComplexMatrix a = chol_mat.hermitian () * chol_mat;
+          ComplexMatrix a1 (n+1, n+1);
+          for (octave_idx_type k = 0; k < n+1; k++)
+            for (octave_idx_type l = 0; l < n+1; l++)
+              {
+                if (l == j)
+                  a1(k, l) = u(k);
+                else if (k == j)
+                  a1(k, l) = std::conj (u(l));
+                else
+                  a1(k, l) = a(k < j ? k : k-1, l < j ? l : l-1);
+              }
+          info = init (a1, false);
+          if (info) info = 1;
+        }
+    }
+
+  return info;
+}
+
+void
+ComplexCHOL::delete_sym (octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type n = chol_mat.rows ();
+  
+  if (j < 0 || j > n-1)
+    (*current_liboctave_error_handler) ("choldelete: index out of range");
+  else
+    {
+      ComplexMatrix a = chol_mat.hermitian () * chol_mat;
+      a.delete_elements (1, idx_vector (j));
+      a.delete_elements (0, idx_vector (j));
+      init (a, false);
+    }
+}
+
+void
+ComplexCHOL::shift_sym (octave_idx_type i, octave_idx_type j)
+{
+  warn_qrupdate_once ();
+
+  octave_idx_type n = chol_mat.rows ();
+  
+  if (i < 0 || i > n-1 || j < 0 || j > n-1) 
+    (*current_liboctave_error_handler) ("cholshift: index out of range");
+  else
+    {
+      ComplexMatrix a = chol_mat.hermitian () * chol_mat;
+      Array<octave_idx_type> p (n);
+      for (octave_idx_type k = 0; k < n; k++) p(k) = k;
+      if (i < j)
+        {
+          for (octave_idx_type k = i; k < j; k++) p(k) = k+1;
+          p(j) = i;
+        }
+      else if (j < i)
+        {
+          p(j) = i;
+          for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1;
+        }
+
+      init (a.index (idx_vector (p), idx_vector (p)), false);
+    }
+}
+
 #endif
 
 ComplexMatrix