diff liboctave/dbleCHOL.cc @ 8547:d66c9b6e506a

imported patch qrupdate.diff
author Jaroslav Hajek <highegg@gmail.com>
date Tue, 20 Jan 2009 21:16:42 +0100
parents 25bc2d31e1bf
children a6edd5c23cb5
line wrap: on
line diff
--- a/liboctave/dbleCHOL.cc
+++ b/liboctave/dbleCHOL.cc
@@ -2,6 +2,7 @@
 
 Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 2007
               John W. Eaton
+Copyright (C) 2008, 2009 Jaroslav Hajek
 
 This file is part of Octave.
 
@@ -21,8 +22,6 @@
 
 */
 
-// updating/downdating by Jaroslav Hajek 2008
-
 #ifdef HAVE_CONFIG_H
 #include <config.h>
 #endif
@@ -52,23 +51,30 @@
 			     double*, const octave_idx_type&, const double&,
 			     double&, double*, octave_idx_type*, 
 			     octave_idx_type& F77_CHAR_ARG_LEN_DECL);
+#ifdef HAVE_QRUPDATE
+
   F77_RET_T
-  F77_FUNC (dch1up, DCH1UP) (const octave_idx_type&, double*, double*, double*);
+  F77_FUNC (dch1up, DCH1UP) (const octave_idx_type&, double*, const octave_idx_type&,
+                             double*, double*);
 
   F77_RET_T
-  F77_FUNC (dch1dn, DCH1DN) (const octave_idx_type&, double*, double*, double*, 
+  F77_FUNC (dch1dn, DCH1DN) (const octave_idx_type&, double*, const octave_idx_type&,
+                             double*, double*, octave_idx_type&);
+
+  F77_RET_T
+  F77_FUNC (dchinx, DCHINX) (const octave_idx_type&, double*, const octave_idx_type&,
+                             const octave_idx_type&, double*, double*, 
                              octave_idx_type&);
 
   F77_RET_T
-  F77_FUNC (dqrshc, DQRSHC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&,
-                             double*, double*, const octave_idx_type&, const octave_idx_type&);
+  F77_FUNC (dchdex, DCHDEX) (const octave_idx_type&, double*, const octave_idx_type&,
+                             const octave_idx_type&, double*);
 
   F77_RET_T
-  F77_FUNC (dchinx, DCHINX) (const octave_idx_type&, const double*, double*, const octave_idx_type&,
-                             const double*, octave_idx_type&);
-
-  F77_RET_T
-  F77_FUNC (dchdex, DCHDEX) (const octave_idx_type&, const double*, double*, const octave_idx_type&);
+  F77_FUNC (dchshx, DCHSHX) (const octave_idx_type&, double*, const octave_idx_type&,
+                             const octave_idx_type&, const octave_idx_type&, 
+                             double*);
+#endif
 }
 
 octave_idx_type
@@ -186,26 +192,28 @@
     (*current_liboctave_error_handler) ("CHOL requires square matrix");
 }
 
+#ifdef HAVE_QRUPDATE
+
 void
-CHOL::update (const Matrix& u)
+CHOL::update (const ColumnVector& u)
 {
   octave_idx_type n = chol_mat.rows ();
 
   if (u.length () == n)
     {
-      Matrix tmp = u;
+      ColumnVector utmp = u;
 
       OCTAVE_LOCAL_BUFFER (double, w, n);
 
-      F77_XFCN (dch1up, DCH1UP, (n, chol_mat.fortran_vec (),
-				 tmp.fortran_vec (), w));
+      F77_XFCN (dch1up, DCH1UP, (n, chol_mat.fortran_vec (), chol_mat.rows (),
+                                 utmp.fortran_vec (), w));
     }
   else
-    (*current_liboctave_error_handler) ("CHOL update dimension mismatch");
+    (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
 }
 
 octave_idx_type
-CHOL::downdate (const Matrix& u)
+CHOL::downdate (const ColumnVector& u)
 {
   octave_idx_type info = -1;
 
@@ -213,38 +221,40 @@
 
   if (u.length () == n)
     {
-      Matrix tmp = u;
+      ColumnVector utmp = u;
 
       OCTAVE_LOCAL_BUFFER (double, w, n);
 
-      F77_XFCN (dch1dn, DCH1DN, (n, chol_mat.fortran_vec (),
-				 tmp.fortran_vec (), w, info));
+      F77_XFCN (dch1dn, DCH1DN, (n, chol_mat.fortran_vec (), chol_mat.rows (),
+                                 utmp.fortran_vec (), w, info));
     }
   else
-    (*current_liboctave_error_handler) ("CHOL downdate dimension mismatch");
+    (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
 
   return info;
 }
 
 octave_idx_type
-CHOL::insert_sym (const Matrix& u, octave_idx_type j)
+CHOL::insert_sym (const ColumnVector& u, octave_idx_type j)
 {
   octave_idx_type info = -1;
 
   octave_idx_type n = chol_mat.rows ();
   
-  if (u.length () != n+1)
-    (*current_liboctave_error_handler) ("CHOL insert dimension mismatch");
+  if (u.length () != n + 1)
+    (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
   else if (j < 0 || j > n)
-    (*current_liboctave_error_handler) ("CHOL insert index out of range");
+    (*current_liboctave_error_handler) ("cholinsert: index out of range");
   else
     {
-      Matrix chol_mat1 (n+1, n+1);
+      ColumnVector utmp = u;
+
+      OCTAVE_LOCAL_BUFFER (double, w, n);
 
-      F77_XFCN (dchinx, DCHINX, (n, chol_mat.data (), chol_mat1.fortran_vec (), 
-                                 j+1, u.data (), info));
+      chol_mat.resize (n+1, n+1);
 
-      chol_mat = chol_mat1;
+      F77_XFCN (dchinx, DCHINX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
+                                 j + 1, utmp.fortran_vec (), w, info));
     }
 
   return info;
@@ -256,14 +266,15 @@
   octave_idx_type n = chol_mat.rows ();
   
   if (j < 0 || j > n-1)
-    (*current_liboctave_error_handler) ("CHOL delete index out of range");
+    (*current_liboctave_error_handler) ("choldelete: index out of range");
   else
     {
-      Matrix chol_mat1 (n-1, n-1);
+      OCTAVE_LOCAL_BUFFER (double, w, n);
 
-      F77_XFCN (dchdex, DCHDEX, (n, chol_mat.data (), chol_mat1.fortran_vec (), j+1));
+      F77_XFCN (dchdex, DCHDEX, (n, chol_mat.fortran_vec (), chol_mat.rows (), 
+                                 j + 1, w));
 
-      chol_mat = chol_mat1;
+      chol_mat.resize (n-1, n-1);
     }
 }
 
@@ -271,14 +282,20 @@
 CHOL::shift_sym (octave_idx_type i, octave_idx_type j)
 {
   octave_idx_type n = chol_mat.rows ();
-  double dummy;
   
   if (i < 0 || i > n-1 || j < 0 || j > n-1) 
-    (*current_liboctave_error_handler) ("CHOL shift index out of range");
+    (*current_liboctave_error_handler) ("cholshift: index out of range");
   else
-    F77_XFCN (dqrshc, DQRSHC, (0, n, n, &dummy, chol_mat.fortran_vec (), i+1, j+1));
+    {
+      OCTAVE_LOCAL_BUFFER (double, w, 2*n);
+
+      F77_XFCN (dchshx, DCHSHX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
+                                 i + 1, j + 1, w));
+    }
 }
 
+#endif
+
 Matrix
 chol2inv (const Matrix& r)
 {