--- a/liboctave/dbleQRP.cc	Tue Jan 27 08:15:08 2009 +0100
+++ b/liboctave/dbleQRP.cc	Tue Jan 27 12:40:06 2009 +0100
@@ -30,6 +30,7 @@
 #include "dbleQRP.h"
 #include "f77-fcn.h"
 #include "lo-error.h"
+#include "oct-locbuf.h"
 
 extern "C"
 {
@@ -37,11 +38,6 @@
   F77_FUNC (dgeqp3, DGEQP3) (const octave_idx_type&, const octave_idx_type&, double*,
 			     const octave_idx_type&, octave_idx_type*, double*, double*,
                              const octave_idx_type&, octave_idx_type&);
-
-  F77_RET_T
-  F77_FUNC (dorgqr, DORGQR) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&,
-			     double*, const octave_idx_type&, double*, double*,
-			     const octave_idx_type&, octave_idx_type&);
 }
 
 // It would be best to share some of this code with QR class...
@@ -60,38 +56,32 @@
   octave_idx_type m = a.rows ();
   octave_idx_type n = a.cols ();
 
-  if (m == 0 || n == 0)
-    {
-      (*current_liboctave_error_handler) ("QR must have non-empty matrix");
-      return;
-    }
-
   octave_idx_type min_mn = m < n ? m : n;
-  Array<double> tau (min_mn);
-  double *ptau = tau.fortran_vec ();
+  OCTAVE_LOCAL_BUFFER (double, tau, min_mn);
 
   octave_idx_type info = 0;
 
-  Matrix A_fact = a;
-  if (m > n && qr_type != QR::economy)
-    A_fact.resize (m, m, 0.0);
-
-  double *tmp_data = A_fact.fortran_vec ();
+  Matrix 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 ();
 
-  double rlwork = 0;
-  // Workspace query...
-  F77_XFCN (dgeqp3, DGEQP3, (m, n, tmp_data, m, pjpvt, ptau, &rlwork, -1, info));
+  if (m > 0)
+    {
+      // workspace query.
+      double rlwork;
+      F77_XFCN (dgeqp3, DGEQP3, (m, n, afact.fortran_vec (), m, jpvt.fortran_vec (),
+                                 tau, &rlwork, -1, info));
 
-  octave_idx_type lwork = rlwork;
-  Array<double> work (lwork);
-  double *pwork = work.fortran_vec ();
-
-  // Code to enforce a certain permutation could go here...
-
-  F77_XFCN (dgeqp3, DGEQP3, (m, n, tmp_data, m, pjpvt, ptau, pwork, lwork, info));
+      // allocate buffer and do the job.
+      octave_idx_type lwork = rlwork; lwork = std::max (lwork, 1);
+      OCTAVE_LOCAL_BUFFER (double, work, lwork);
+      F77_XFCN (dgeqp3, DGEQP3, (m, n, afact.fortran_vec (), m, jpvt.fortran_vec (),
+                                 tau, work, lwork, 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!)
@@ -99,25 +89,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 (dorgqr, DORGQR, (m, n2, min_mn, tmp_data, m, ptau,
-			     pwork, lwork, info));
-
-  q = A_fact;
-  q.resize (m, n2);
+  form (n, afact, tau, qr_type);
 }
 
 ColumnVector