# HG changeset patch
# User jwe
# Date 857182473 0
# Node ID d9d00d7e271e1346e7614764a6b2ea3894296501
# Parent  ac1427f5c9e635d1998f09feb3c315f72490b22b
[project @ 1997-03-01 02:14:33 by jwe]

diff --git a/liboctave/CmplxQR.cc b/liboctave/CmplxQR.cc
--- a/liboctave/CmplxQR.cc
+++ b/liboctave/CmplxQR.cc
@@ -45,6 +45,13 @@
 }
 
 ComplexQR::ComplexQR (const ComplexMatrix& a, QR::type qr_type)
+  : q (), r ()
+{
+  init (a, qr_type);
+}
+
+void
+ComplexQR::init (const ComplexMatrix& a, QR::type qr_type)
 {
   int m = a.rows ();
   int n = a.cols ();
diff --git a/liboctave/CmplxQR.h b/liboctave/CmplxQR.h
--- a/liboctave/CmplxQR.h
+++ b/liboctave/CmplxQR.h
@@ -39,7 +39,7 @@
 
   ComplexQR (void) : q (), r () { }
 
-  ComplexQR (const ComplexMatrix& A, QR::type qr_type = QR::std);
+  ComplexQR (const ComplexMatrix&, QR::type = QR::std);
 
   ComplexQR (const ComplexQR& a) : q (a.q), r (a.r) { }
 
@@ -55,10 +55,13 @@
 
   ~ComplexQR (void) { }
 
+  void init (const ComplexMatrix&, QR::type = QR::std);
+
   ComplexMatrix Q (void) const { return q; }
+
   ComplexMatrix R (void) const { return r; }
 
-  friend ostream&  operator << (ostream& os, const ComplexQR& a);
+  friend ostream&  operator << (ostream&, const ComplexQR&);
 
 protected:
 
diff --git a/liboctave/CmplxQRP.cc b/liboctave/CmplxQRP.cc
--- a/liboctave/CmplxQRP.cc
+++ b/liboctave/CmplxQRP.cc
@@ -49,6 +49,13 @@
 // It would be best to share some of this code with ComplexQR class...
 
 ComplexQRP::ComplexQRP (const ComplexMatrix& a, QR::type qr_type)
+  : ComplexQR (), p ()
+{
+  init (a, qr_type);
+}
+
+void
+ComplexQRP::init (const ComplexMatrix& a, QR::type qr_type)
 {
   assert (qr_type != QR::raw);
 
@@ -62,24 +69,18 @@
       return;
     }
 
-  int min_mn = m < n ? m : n;
-  Array<Complex> tau (min_mn);
+  Array<Complex> tau (m < n ? m : n);
   Complex *ptau = tau.fortran_vec ();
 
-  int lwork = n;
+  int lwork = 3*n > 32*m ? 3*n : 32*m;
   Array<Complex> work (lwork);
   Complex *pwork = work.fortran_vec ();
 
   int info = 0;
 
-  ComplexMatrix A_fact;
+  ComplexMatrix A_fact = a;
   if (m > n)
-    {
-      A_fact.resize (m, m);
-      A_fact.insert (a, 0, 0);
-    }
-  else
-    A_fact = a;
+    A_fact.resize (m, m, 0.0);
 
   Complex *tmp_data = A_fact.fortran_vec ();
 
@@ -101,7 +102,7 @@
       // Form Permutation matrix (if economy is requested, return the
       // indices only!)
 
-      if (qr_type == QR::economy && m > n)
+      if (qr_type == QR::economy)
 	{
 	  p.resize (1, n, 0.0);
 	  for (int j = 0; j < n; j++)
@@ -114,18 +115,12 @@
 	    p.elem (jpvt.elem (j) - 1, j) = 1.0;
 	}
 
-      volatile int n2;
+      if (qr_type == QR::economy && m > n)
+	r.resize (n, n, 0.0);
+      else
+	r.resize (m, n, 0.0);
 
-      if (qr_type == QR::economy && m > n)
-	{
-	  n2 = n;
-	  r.resize (n, n, 0.0);
-	}
-      else
-	{
-	  n2 = m;
-	  r.resize (m, n, 0.0);
-	}
+      int min_mn = m < n ? m : n;
 
       for (int j = 0; j < n; j++)
 	{
@@ -134,11 +129,11 @@
 	    r.elem (i, j) = A_fact.elem (i, j);
 	}
 
-      lwork = 32*m;
-      work.resize (lwork);
-      Complex *pwork = work.fortran_vec ();
+      int n2 = m;
+      if (qr_type == QR::economy)
+	n2 = min_mn;
 
-      F77_XFCN (zungqr, ZUNGQR, (m, m, min_mn, tmp_data, m, ptau,
+      F77_XFCN (zungqr, ZUNGQR, (m, n2, min_mn, tmp_data, m, ptau,
 				 pwork, lwork, info));
 
       if (f77_exception_encountered)
diff --git a/liboctave/CmplxQRP.h b/liboctave/CmplxQRP.h
--- a/liboctave/CmplxQRP.h
+++ b/liboctave/CmplxQRP.h
@@ -38,7 +38,7 @@
 
   ComplexQRP (void) : ComplexQR (), p () { }
 
-  ComplexQRP (const ComplexMatrix& A, QR::type qr_type = QR::std);
+  ComplexQRP (const ComplexMatrix&, QR::type = QR::std);
 
   ComplexQRP (const ComplexQRP& a) : ComplexQR (a), p (a.p) { }
 
@@ -54,9 +54,11 @@
 
   ~ComplexQRP (void) { }
 
+  void init (const ComplexMatrix&, QR::type = QR::std);
+
   Matrix P (void) const { return p; }
 
-  friend ostream&  operator << (ostream& os, const ComplexQRP& a);
+  friend ostream&  operator << (ostream&, const ComplexQRP&);
 
 private:
 
diff --git a/liboctave/dbleQR.cc b/liboctave/dbleQR.cc
--- a/liboctave/dbleQR.cc
+++ b/liboctave/dbleQR.cc
@@ -45,6 +45,13 @@
 }
 
 QR::QR (const Matrix& a, QR::type qr_type)
+  : q (), r ()
+{
+  init (a, qr_type);
+}
+
+void
+QR::init (const Matrix& a, QR::type qr_type)
 {
   int m = a.rows ();
   int n = a.cols ();
diff --git a/liboctave/dbleQR.h b/liboctave/dbleQR.h
--- a/liboctave/dbleQR.h
+++ b/liboctave/dbleQR.h
@@ -45,7 +45,7 @@
 
   QR (void) : q (), r () { }
 
-  QR (const Matrix& A, type qr_type = QR::std);
+  QR (const Matrix&, QR::type = QR::std);
 
   QR (const QR& a) : q (a.q), r (a.r) { }
 
@@ -61,11 +61,13 @@
 
   ~QR (void) { }
 
+  void init (const Matrix&, QR::type);
+
   Matrix Q (void) const { return q; }
 
   Matrix R (void) const { return r; }
 
-  friend ostream&  operator << (ostream& os, const QR& a);
+  friend ostream&  operator << (ostream&, const QR&);
 
 protected:
 
diff --git a/liboctave/dbleQRP.cc b/liboctave/dbleQRP.cc
--- a/liboctave/dbleQRP.cc
+++ b/liboctave/dbleQRP.cc
@@ -49,6 +49,13 @@
 // It would be best to share some of this code with QR class...
 
 QRP::QRP (const Matrix& a, QR::type qr_type)
+  : QR (), p ()
+{
+  init (a, qr_type);
+}
+
+void
+QRP::init (const Matrix& a, QR::type qr_type)
 {
   assert (qr_type != QR::raw);
 
@@ -61,24 +68,18 @@
       return;
     }
 
-  int min_mn = m < n ? m : n;
-  Array<double> tau (min_mn);
+  Array<double> tau (m < n ? m : n);
   double *ptau = tau.fortran_vec ();
 
-  int lwork = 3*n;
+  int lwork = 3*n > 32*m ? 3*n : 32*m;
   Array<double> work (lwork);
   double *pwork = work.fortran_vec ();
 
   int info = 0;
 
-  Matrix A_fact;
+  Matrix A_fact = a;
   if (m > n)
-    {
-      A_fact.resize (m, m);
-      A_fact.insert (a, 0, 0);
-    }
-  else
-    A_fact = a;
+    A_fact.resize (m, m, 0.0);
 
   double *tmp_data = A_fact.fortran_vec ();
 
@@ -96,7 +97,7 @@
       // Form Permutation matrix (if economy is requested, return the
       // indices only!)
 
-      if (qr_type == QR::economy && m > n)
+      if (qr_type == QR::economy)
 	{
 	  p.resize (1, n, 0.0);
 	  for (int j = 0; j < n; j++)
@@ -109,18 +110,12 @@
 	    p.elem (jpvt.elem (j) - 1, j) = 1.0;
 	}
 
-      volatile int n2;
+      if (qr_type == QR::economy && m > n)
+	r.resize (n, n, 0.0);
+      else
+	r.resize (m, n, 0.0);
 
-      if (qr_type == QR::economy && m > n)
-	{
-	  n2 = n;
-	  r.resize (n, n, 0.0);
-	}
-      else
-	{
-	  n2 = m;
-	  r.resize (m, n, 0.0);
-	}
+      int min_mn = m < n ? m : n;
 
       for (int j = 0; j < n; j++)
 	{
@@ -129,11 +124,11 @@
 	    r.elem (i, j) = A_fact.elem (i, j);
 	}
 
-      lwork = 32*m;
-      work.resize (lwork);
-      double *pwork = work.fortran_vec ();
+      int n2 = m;
+      if (qr_type == QR::economy)
+	n2 = min_mn;
 
-      F77_XFCN (dorgqr, DORGQR, (m, m, min_mn, tmp_data, m, ptau,
+      F77_XFCN (dorgqr, DORGQR, (m, n2, min_mn, tmp_data, m, ptau,
 				 pwork, lwork, info));
 
       if (f77_exception_encountered)
diff --git a/liboctave/dbleQRP.h b/liboctave/dbleQRP.h
--- a/liboctave/dbleQRP.h
+++ b/liboctave/dbleQRP.h
@@ -38,7 +38,7 @@
 
   QRP (void) : QR (), p () { }
 
-  QRP (const Matrix& A, QR::type qr_type = QR::std);
+  QRP (const Matrix&, QR::type = QR::std);
 
   QRP (const QRP& a) : QR (a), p (a.p) { }
 
@@ -55,9 +55,11 @@
 
   ~QRP (void) { }
 
+  void init (const Matrix&, QR::type = QR::std);
+
   Matrix P (void) const { return p; }
 
-  friend ostream&  operator << (ostream& os, const QRP& a);
+  friend ostream&  operator << (ostream&, const QRP&);
 
 protected: