Mercurial > hg > octave-max

changeset 1922:8a57554f3142

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1996-02-11 02:53:20 by jwe]
author jwe
date Sun, 11 Feb 1996 02:53:20 +0000
parents ce0db0ca0729
children 1d4f67de06f4
files liboctave/CmplxQRP.cc liboctave/CmplxQRP.h liboctave/dbleQRP.cc liboctave/dbleQRP.h
diffstat 4 files changed, 176 insertions(+), 101 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/CmplxQRP.cc
+++ b/liboctave/CmplxQRP.cc
@@ -53,6 +53,11 @@
 {
   assert (qr_type != QR::raw);
 
+  tau = 0;
+  work = 0;
+  tmp_data = 0;
+  jpvt = 0;
+
   int m = a.rows ();
   int n = a.cols ();
 
@@ -63,11 +68,12 @@
       return;
     }
 
-  Complex *tmp_data;
   int min_mn = m < n ? m : n;
-  Complex *tau = new Complex[min_mn];
+  tau = new Complex[min_mn];
+
   int lwork = n;
-  Complex *work = new Complex[lwork];
+  work = new Complex[lwork];
+
   int info = 0;
 
   if (m > n)
@@ -78,10 +84,9 @@
   else
     tmp_data = dup (a.data (), a.length ());
 
-
   work = new Complex[n];
-  double *rwork = new double[2*n];
-  int *jpvt = new int[n];
+  rwork = new double[2*n];
+  jpvt = new int[n];
 
   // Clear Pivot vector (code to enforce a certain permutation would
   // go here...)
@@ -89,58 +94,85 @@
   for (int i = 0; i < n; i++)
     jpvt[i] = 0;      
 
-  F77_FCN (zgeqpf, ZGEQPF) (m, n, tmp_data, m, jpvt, tau, work, rwork,
-			    info);
+  F77_XFCN (zgeqpf, ZGEQPF, (m, n, tmp_data, m, jpvt, tau, work, rwork,
+			     info));
 
-  // Form Permutation matrix (if economy is requested, return the
-  // indices only!)
+  delete [] work;
+  work = 0;
 
-  if (qr_type == QR::economy && m > n)
-    {
-      p.resize (1, n, 0.0);
-      for (int j = 0; j < n; j++)
-	p.elem (0, j) = jpvt[j];
-    }
+  delete [] rwork;
+  rwork = 0;
+
+  if (f77_exception_encountered)
+    (*current_liboctave_error_handler) ("unrecoverable error in zgeqpf");
   else
     {
-      p.resize (n, n, 0.0);
+      // Form Permutation matrix (if economy is requested, return the
+      // indices only!)
+
+      if (qr_type == QR::economy && m > n)
+	{
+	  p.resize (1, n, 0.0);
+	  for (int j = 0; j < n; j++)
+	    p.elem (0, j) = jpvt[j];
+	}
+      else
+	{
+	  p.resize (n, n, 0.0);
+	  for (int j = 0; j < n; j++)
+	    p.elem (jpvt[j]-1, j) = 1.0;
+	}
+
+      delete [] jpvt;
+      jpvt = 0;
+
+      volatile int n2;
+
+      if (qr_type == QR::economy && m > n)
+	{
+	  n2 = n;
+	  r.resize (n, n, 0.0);
+	}
+      else
+	{
+	  n2 = m;
+	  r.resize (m, n, 0.0);
+	}
+
       for (int j = 0; j < n; j++)
-	p.elem (jpvt[j]-1, j) = 1.0;
+	{
+	  int limit = j < min_mn-1 ? j : min_mn-1;
+	  for (int i = 0; i <= limit; i++)
+	    r.elem (i, j) = tmp_data[m*j+i];
+	}
+
+      lwork = 32*m;
+      work = new Complex[lwork];
+
+      F77_XFCN (zungqr, ZUNGQR, (m, m, min_mn, tmp_data, m, tau, work,
+				 lwork, info));
+
+      if (f77_exception_encountered)
+	(*current_liboctave_error_handler) ("unrecoverable error in zungqr");
+      else
+	{
+	  q = ComplexMatrix (tmp_data, m, m);
+	  tmp_data = 0;
+	  q.resize (m, n2);
+	}
     }
 
+  delete [] tau;
+  tau = 0;
+
+  delete [] work;
+  work = 0;
+
+  delete [] tmp_data;
+  tmp_data = 0;
+
   delete [] jpvt;
-  delete [] work;
-
-  int n2;
-  if (qr_type == QR::economy && m > n)
-    {
-      n2 = n;
-      r.resize (n, n, 0.0);
-    }
-  else
-    {
-      n2 = m;
-      r.resize (m, n, 0.0);
-    }
-
-  for (int j = 0; j < n; j++)
-    {
-      int limit = j < min_mn-1 ? j : min_mn-1;
-      for (int i = 0; i <= limit; i++)
-	r.elem (i, j) = tmp_data[m*j+i];
-    }
-
-  lwork = 32*m;
-  work = new Complex[lwork];
-
-  F77_FCN (zungqr, ZUNGQR) (m, m, min_mn, tmp_data, m, tau, work,
-			    lwork, info);
-
-  q = ComplexMatrix (tmp_data, m, m);
-  q.resize (m, n2);
-
-  delete [] tau;
-  delete [] work;
+  jpvt = 0;
 }
 
 /*
--- a/liboctave/CmplxQRP.h
+++ b/liboctave/CmplxQRP.h
@@ -37,7 +37,7 @@
 {
 public:
 
-  ComplexQRP (void) : ComplexQR (), p () { }
+  ComplexQRP (void) : ComplexQR (), p (), jpvt (0), rwork (0) { }
 
   ComplexQRP (const ComplexMatrix& A, QR::type qr_type = QR::std);
 
@@ -53,6 +53,12 @@
       return *this;
     }
 
+  ~ComplexQRP (void)
+    {
+      delete [] jpvt;
+      delete [] rwork;
+    }
+
   Matrix P (void) const { return p; }
 
   friend ostream&  operator << (ostream& os, const ComplexQRP& a);
@@ -60,6 +66,9 @@
 private:
 
   Matrix p;
+
+  int *jpvt;
+  double *rwork;
 };
 
 #endif
--- a/liboctave/dbleQRP.cc
+++ b/liboctave/dbleQRP.cc
@@ -53,6 +53,11 @@
 {
   assert (qr_type != QR::raw);
 
+  tau = 0;
+  work = 0;
+  tmp_data = 0;
+  jpvt = 0;
+
   int m = a.rows ();
   int n = a.cols ();
 
@@ -62,11 +67,12 @@
       return;
     }
 
-  double *tmp_data;
   int min_mn = m < n ? m : n;
-  double *tau = new double[min_mn];
+  tau = new double[min_mn];
+
   int lwork = 3*n;
-  double *work = new double[lwork];
+  work = new double[lwork];
+
   int info = 0;
 
   if (m > n)
@@ -77,7 +83,7 @@
   else
     tmp_data = dup (a.data (), a.length ());
 
-  int *jpvt = new int[n];
+  jpvt = new int[n];
 
   // Clear Pivot vector (code to enforce a certain permutation would
   // go here...)
@@ -85,57 +91,81 @@
   for (int i = 0; i < n; i++)
     jpvt[i] = 0;
 
-  F77_FCN (dgeqpf, DGEQPF) (m, n, tmp_data, m, jpvt, tau, work, info);
-
-  // Form Permutation matrix (if economy is requested, return the
-  // indices only!)
+  F77_XFCN (dgeqpf, DGEQPF, (m, n, tmp_data, m, jpvt, tau, work, info));
 
-  if (qr_type == QR::economy && m > n)
-    {
-      p.resize (1, n, 0.0);
-      for (int j = 0; j < n; j++)
-	p.elem (0, j) = jpvt[j];
-    }
+  delete [] work;
+  work = 0;
+
+  if (f77_exception_encountered)
+    (*current_liboctave_error_handler) ("unrecoverable error in dgeqpf");
   else
     {
-      p.resize (n, n, 0.0);
+      // Form Permutation matrix (if economy is requested, return the
+      // indices only!)
+
+      if (qr_type == QR::economy && m > n)
+	{
+	  p.resize (1, n, 0.0);
+	  for (int j = 0; j < n; j++)
+	    p.elem (0, j) = jpvt[j];
+	}
+      else
+	{
+	  p.resize (n, n, 0.0);
+	  for (int j = 0; j < n; j++)
+	    p.elem (jpvt[j]-1, j) = 1.0;
+	}
+
+      delete [] jpvt;
+      jpvt = 0;
+
+      volatile int n2;
+
+      if (qr_type == QR::economy && m > n)
+	{
+	  n2 = n;
+	  r.resize (n, n, 0.0);
+	}
+      else
+	{
+	  n2 = m;
+	  r.resize (m, n, 0.0);
+	}
+
       for (int j = 0; j < n; j++)
-	p.elem (jpvt[j]-1, j) = 1.0;
+	{
+	  int limit = j < min_mn-1 ? j : min_mn-1;
+	  for (int i = 0; i <= limit; i++)
+	    r.elem (i, j) = tmp_data[m*j+i];
+	}
+
+      lwork = 32*m;
+      work = new double[lwork];
+
+      F77_XFCN (dorgqr, DORGQR, (m, m, min_mn, tmp_data, m, tau, work,
+				 lwork, info));
+
+      if (f77_exception_encountered)
+	(*current_liboctave_error_handler) ("unrecoverable error in dorgqr");
+      else
+	{
+	  q = Matrix (tmp_data, m, m);
+	  tmp_data = 0;
+	  q.resize (m, n2);
+	}
     }
 
+  delete [] tau;
+  tau = 0;
+
+  delete [] work;
+  work = 0;
+
+  delete [] tmp_data;
+  tmp_data = 0;
+
   delete [] jpvt;
-  delete [] work;
-
-  int n2;
-  if (qr_type == QR::economy && m > n)
-    {
-      n2 = n;
-      r.resize (n, n, 0.0);
-    }
-  else
-    {
-      n2 = m;
-      r.resize (m, n, 0.0);
-    }
-
-  for (int j = 0; j < n; j++)
-    {
-      int limit = j < min_mn-1 ? j : min_mn-1;
-      for (int i = 0; i <= limit; i++)
-	r.elem (i, j) = tmp_data[m*j+i];
-    }
-
-  lwork = 32*m;
-  work = new double[lwork];
-
-  F77_FCN (dorgqr, DORGQR) (m, m, min_mn, tmp_data, m, tau, work,
-			    lwork, info);
-
-  q = Matrix (tmp_data, m, m);
-  q.resize (m, n2);
-
-  delete [] tau;
-  delete [] work;
+  jpvt = 0;
 }
 
 /*
--- a/liboctave/dbleQRP.h
+++ b/liboctave/dbleQRP.h
@@ -37,7 +37,7 @@
 {
 public:
 
-  QRP (void) : QR (), p () { }
+  QRP (void) : QR (), p (), jpvt (0) { }
 
   QRP (const Matrix& A, QR::type qr_type = QR::std);
 
@@ -54,13 +54,17 @@
       return *this;
     }
 
+  ~QRP (void) { delete [] jpvt; }
+
   Matrix P (void) const { return p; }
 
   friend ostream&  operator << (ostream& os, const QRP& a);
 
-private:
+protected:
 
   Matrix p;
+
+  int *jpvt;
 };
 
 #endif