diff liboctave/CmplxQRP.cc @ 2763:d9d00d7e271e

[project @ 1997-03-01 02:14:33 by jwe]
author jwe
date Sat, 01 Mar 1997 02:14:33 +0000
parents 3db30620918e
children eedc2f3f61f7
line wrap: on
line diff
--- a/liboctave/CmplxQRP.cc	Fri Feb 28 08:26:43 1997 +0000
+++ b/liboctave/CmplxQRP.cc	Sat Mar 01 02:14:33 1997 +0000
@@ -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)