diff liboctave/fCmplxQRP.cc @ 8597:c86718093c1b

improve & fix QR classes
author Jaroslav Hajek <highegg@gmail.com>
date Tue, 27 Jan 2009 12:40:06 +0100
parents e3c9102431a9
children 20dfb885f877
line wrap: on
line diff
--- a/liboctave/fCmplxQRP.cc	Tue Jan 27 08:15:08 2009 +0100
+++ b/liboctave/fCmplxQRP.cc	Tue Jan 27 12:40:06 2009 +0100
@@ -30,6 +30,7 @@
 #include "fCmplxQRP.h"
 #include "f77-fcn.h"
 #include "lo-error.h"
+#include "oct-locbuf.h"
 
 extern "C"
 {
@@ -37,11 +38,6 @@
   F77_FUNC (cgeqp3, CGEQP3) (const octave_idx_type&, const octave_idx_type&, FloatComplex*,
 			     const octave_idx_type&, octave_idx_type*, FloatComplex*, FloatComplex*,
 			     const octave_idx_type&, float*, octave_idx_type&);
-
-  F77_RET_T
-  F77_FUNC (cungqr, CUNGQR) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&,
-			     FloatComplex*, const octave_idx_type&, FloatComplex*,
-			     FloatComplex*, const octave_idx_type&, octave_idx_type&);
 }
 
 // It would be best to share some of this code with FloatComplexQR class...
@@ -60,44 +56,34 @@
   octave_idx_type m = a.rows ();
   octave_idx_type n = a.cols ();
 
-  if (m == 0 || n == 0)
-    {
-      (*current_liboctave_error_handler)
-	("FloatComplexQR must have non-empty matrix");
-      return;
-    }
-
   octave_idx_type min_mn = m < n ? m : n;
-  Array<FloatComplex> tau (min_mn);
-  FloatComplex *ptau = tau.fortran_vec ();
+  OCTAVE_LOCAL_BUFFER (FloatComplex, tau, min_mn);
 
   octave_idx_type info = 0;
 
-  FloatComplexMatrix A_fact = a;
-  if (m > n && qr_type != QR::economy)
-    A_fact.resize (m, m, 0.0);
-
-  FloatComplex *tmp_data = A_fact.fortran_vec ();
-
-  Array<float> rwork (2*n);
-  float *prwork = rwork.fortran_vec ();
+  FloatComplexMatrix afact = a;
+  if (m > n && qr_type == QR::std)
+    afact.resize (m, m);
 
   MArray<octave_idx_type> jpvt (n, 0);
-  octave_idx_type *pjpvt = jpvt.fortran_vec ();
 
-  FloatComplex rlwork = 0;
-  // Workspace query...
-  F77_XFCN (cgeqp3, CGEQP3, (m, n, tmp_data, m, pjpvt, ptau, &rlwork,
-			     -1, prwork, info));
+  if (m > 0)
+    {
+      OCTAVE_LOCAL_BUFFER (float, rwork, 2*n);
 
-  octave_idx_type lwork = rlwork.real ();
-  Array<FloatComplex> work (lwork);
-  FloatComplex *pwork = work.fortran_vec ();
+      // workspace query.
+      FloatComplex clwork;
+      F77_XFCN (cgeqp3, CGEQP3, (m, n, afact.fortran_vec (), m, jpvt.fortran_vec (),
+                                 tau, &clwork, -1, rwork, info));
 
-  // Code to enforce a certain permutation could go here...
-
-  F77_XFCN (cgeqp3, CGEQP3, (m, n, tmp_data, m, pjpvt, ptau, pwork,
-			     lwork, prwork, info));
+      // allocate buffer and do the job.
+      octave_idx_type lwork = clwork.real (); lwork = std::max (lwork, 1);
+      OCTAVE_LOCAL_BUFFER (FloatComplex, work, lwork);
+      F77_XFCN (cgeqp3, CGEQP3, (m, n, afact.fortran_vec (), m, jpvt.fortran_vec (),
+                                 tau, work, lwork, rwork, info));
+    }
+  else
+    for (octave_idx_type i = 0; i < n; i++) jpvt(i) = i+1;
 
   // Form Permutation matrix (if economy is requested, return the
   // indices only!)
@@ -105,25 +91,8 @@
   jpvt -= 1;
   p = PermMatrix (jpvt, true);
 
-  octave_idx_type n2 = (qr_type == QR::economy) ? min_mn : m;
 
-  if (qr_type == QR::economy && m > n)
-    r.resize (n, n, 0.0);
-  else
-    r.resize (m, n, 0.0);
-
-  for (octave_idx_type j = 0; j < n; j++)
-    {
-      octave_idx_type limit = j < min_mn-1 ? j : min_mn-1;
-      for (octave_idx_type i = 0; i <= limit; i++)
-	r.elem (i, j) = A_fact.elem (i, j);
-    }
-
-  F77_XFCN (cungqr, CUNGQR, (m, n2, min_mn, tmp_data, m, ptau,
-			     pwork, lwork, info));
-
-  q = A_fact;
-  q.resize (m, n2);
+  form (n, afact, tau, qr_type);
 }
 
 FloatColumnVector