diff liboctave/fCmplxLU.cc @ 9708:6f3ffe11d926

implement luupdate
author Jaroslav Hajek <highegg@gmail.com>
date Thu, 08 Oct 2009 16:05:53 +0200
parents 50db3c5175b5
children 4c0cdbe0acca
line wrap: on
line diff
--- a/liboctave/fCmplxLU.cc
+++ b/liboctave/fCmplxLU.cc
@@ -2,6 +2,7 @@
 
 Copyright (C) 1994, 1995, 1996, 1997, 1999, 2002, 2003, 2004, 2005,
               2007, 2008 John W. Eaton
+Copyright (C) 2009 VZLU Prague, a.s.
 
 This file is part of Octave.
 
@@ -28,6 +29,8 @@
 #include "fCmplxLU.h"
 #include "f77-fcn.h"
 #include "lo-error.h"
+#include "oct-locbuf.h"
+#include "fCColVector.h"
 
 // Instantiate the base LU class for the types we need.
 
@@ -43,6 +46,21 @@
   F77_RET_T
   F77_FUNC (cgetrf, CGETRF) (const octave_idx_type&, const octave_idx_type&, FloatComplex*,
 			     const octave_idx_type&, octave_idx_type*, octave_idx_type&);
+
+#ifdef HAVE_QRUPDATE_LUU
+  F77_RET_T
+  F77_FUNC (clu1up, CLU1UP) (const octave_idx_type&, const octave_idx_type&,
+                             FloatComplex *, const octave_idx_type&,
+                             FloatComplex *, const octave_idx_type&,
+                             FloatComplex *, FloatComplex *);
+
+  F77_RET_T
+  F77_FUNC (clup1up, CLUP1UP) (const octave_idx_type&, const octave_idx_type&,
+                               FloatComplex *, const octave_idx_type&,
+                               FloatComplex *, const octave_idx_type&,
+                               octave_idx_type *, const FloatComplex *,
+                               const FloatComplex *, FloatComplex *);
+#endif
 }
 
 FloatComplexLU::FloatComplexLU (const FloatComplexMatrix& a)
@@ -65,6 +83,132 @@
     pipvt[i] -= 1;
 }
 
+#ifdef HAVE_QRUPDATE_LUU
+
+void FloatComplexLU::update (const FloatComplexColumnVector& u, const FloatComplexColumnVector& v)
+{
+  if (packed ())
+    unpack ();
+
+  FloatComplexMatrix& l = l_fact;
+  FloatComplexMatrix& r = a_fact;
+
+  octave_idx_type m = l.rows ();
+  octave_idx_type n = r.columns ();
+  octave_idx_type k = l.columns ();
+
+  if (u.length () == m && v.length () == n)
+    {
+      FloatComplexColumnVector utmp = u, vtmp = v;
+      F77_XFCN (clu1up, CLU1UP, (m, n, l.fortran_vec (), m, r.fortran_vec (), k,
+                                 utmp.fortran_vec (), vtmp.fortran_vec ()));
+    }
+  else
+    (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
+}
+
+void FloatComplexLU::update (const FloatComplexMatrix& u, const FloatComplexMatrix& v)
+{
+  if (packed ())
+    unpack ();
+
+  FloatComplexMatrix& l = l_fact;
+  FloatComplexMatrix& r = a_fact;
+
+  octave_idx_type m = l.rows ();
+  octave_idx_type n = r.columns ();
+  octave_idx_type k = l.columns ();
+
+  if (u.rows () == m && v.rows () == n && u.cols () == v.cols ())
+    {
+      for (volatile octave_idx_type i = 0; i < u.cols (); i++)
+        {
+          FloatComplexColumnVector utmp = u.column (i), vtmp = v.column (i);
+          F77_XFCN (clu1up, CLU1UP, (m, n, l.fortran_vec (), m, r.fortran_vec (), k,
+                                     utmp.fortran_vec (), vtmp.fortran_vec ()));
+        }
+    }
+  else
+    (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
+}
+
+void FloatComplexLU::update_piv (const FloatComplexColumnVector& u, const FloatComplexColumnVector& v)
+{
+  if (packed ())
+    unpack ();
+
+  FloatComplexMatrix& l = l_fact;
+  FloatComplexMatrix& r = a_fact;
+
+  octave_idx_type m = l.rows ();
+  octave_idx_type n = r.columns ();
+  octave_idx_type k = l.columns ();
+
+  if (u.length () == m && v.length () == n)
+    {
+      FloatComplexColumnVector utmp = u, vtmp = v;
+      OCTAVE_LOCAL_BUFFER (FloatComplex, w, m);
+      for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment
+      F77_XFCN (clup1up, CLUP1UP, (m, n, l.fortran_vec (), m, r.fortran_vec (), k,
+                                   ipvt.fortran_vec (), utmp.data (), vtmp.data (), w));
+      for (octave_idx_type i = 0; i < m; i++) ipvt(i) -= 1; // decrement
+    }
+  else
+    (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
+}
+
+void FloatComplexLU::update_piv (const FloatComplexMatrix& u, const FloatComplexMatrix& v)
+{
+  if (packed ())
+    unpack ();
+
+  FloatComplexMatrix& l = l_fact;
+  FloatComplexMatrix& r = a_fact;
+
+  octave_idx_type m = l.rows ();
+  octave_idx_type n = r.columns ();
+  octave_idx_type k = l.columns ();
+
+  if (u.rows () == m && v.rows () == n && u.cols () == v.cols ())
+    {
+      OCTAVE_LOCAL_BUFFER (FloatComplex, w, m);
+      for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment
+      for (volatile octave_idx_type i = 0; i < u.cols (); i++)
+        {
+          FloatComplexColumnVector utmp = u.column (i), vtmp = v.column (i);
+          F77_XFCN (clup1up, CLUP1UP, (m, n, l.fortran_vec (), m, r.fortran_vec (), k,
+                                       ipvt.fortran_vec (), utmp.data (), vtmp.data (), w));
+        }
+      for (octave_idx_type i = 0; i < m; i++) ipvt(i) -= 1; // decrement
+    }
+  else
+    (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
+}
+
+#else
+
+void FloatComplexLU::update (const FloatComplexColumnVector&, const FloatComplexColumnVector&)
+{
+  (*current_liboctave_error_handler) ("luupdate: not available in this version");
+}
+
+void FloatComplexLU::update (const FloatComplexMatrix&, const FloatComplexMatrix&)
+{
+  (*current_liboctave_error_handler) ("luupdate: not available in this version");
+}
+
+void FloatComplexLU::update_piv (const FloatComplexColumnVector&, const FloatComplexColumnVector&)
+{
+  (*current_liboctave_error_handler) ("luupdate: not available in this version");
+}
+
+void FloatComplexLU::update_piv (const FloatComplexMatrix&, const FloatComplexMatrix&)
+{
+  (*current_liboctave_error_handler) ("luupdate: not available in this version");
+}
+
+#endif
+
 /*
 ;;; Local Variables: ***
 ;;; mode: C++ ***