Mercurial > hg > octave-nkf

changeset 1919:b582f7f765e0

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1996-02-11 01:58:47 by jwe]
author jwe
date Sun, 11 Feb 1996 02:02:55 +0000
parents adde795927b5
children 6fdb56b2cb33
files liboctave/CmplxLU.cc liboctave/CmplxLU.h liboctave/dbleLU.cc liboctave/dbleLU.h
diffstat 4 files changed, 112 insertions(+), 68 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/CmplxLU.cc
+++ b/liboctave/CmplxLU.cc
@@ -43,8 +43,12 @@
 
 ComplexLU::ComplexLU (const ComplexMatrix& a)
 {
+  ipvt = 0;
+  pvt = 0;
+
   int a_nr = a.rows ();
   int a_nc = a.cols ();
+
   if (a_nr == 0 || a_nc == 0 || a_nr != a_nc)
     {
       (*current_liboctave_error_handler) ("ComplexLU requires square matrix");
@@ -53,51 +57,60 @@
 
   int n = a_nr;
 
-  int *ipvt = new int [n];
-  int *pvt = new int [n];
-  Complex *tmp_data = dup (a.data (), a.length ());
+  ipvt = new int [n];
+  pvt = new int [n];
+
+  ComplexMatrix A_fact = a;
+  Complex *tmp_data = A_fact.fortran_vec ();
+
   int info = 0;
   Complex *dummy = 0;
 
-  F77_FCN (zgesv, ZGESV) (n, 0, tmp_data, n, ipvt, dummy, n, info);
+  F77_XFCN (zgesv, ZGESV) (n, 0, tmp_data, n, ipvt, dummy, n, info);
 
-  ComplexMatrix A_fact (tmp_data, n, n);
-
-  for (int i = 0; i < n; i++)
+  if (f77_exception_encountered)
+    (*current_liboctave_error_handler) ("unrecoverable error in zgesv");
+  else
     {
-      ipvt[i] -= 1;
-      pvt[i] = i;
-    }
+      for (int i = 0; i < n; i++)
+	{
+	  ipvt[i] -= 1;
+	  pvt[i] = i;
+	}
 
-  for (int i = 0; i < n - 1; i++)
-    {
-      int k = ipvt[i];
-      if (k != i)
+      for (int i = 0; i < n - 1; i++)
 	{
-	  int tmp = pvt[k];
-	  pvt[k] = pvt[i];
-	  pvt[i] = tmp;
+	  int k = ipvt[i];
+	  if (k != i)
+	    {
+	      int tmp = pvt[k];
+	      pvt[k] = pvt[i];
+	      pvt[i] = tmp;
+	    }
+	}
+
+      l.resize (n, n, 0.0);
+      u.resize (n, n, 0.0);
+      p.resize (n, n, 0.0);
+
+      for (int i = 0; i < n; i++)
+	{
+	  p.elem (i, pvt[i]) = 1.0;
+
+	  l.elem (i, i) = 1.0;
+	  for (int j = 0; j < i; j++)
+	    l.elem (i, j) = A_fact.elem (i, j);
+
+	  for (int j = i; j < n; j++)
+	    u.elem (i, j) = A_fact.elem (i, j);
 	}
     }
 
-  l.resize (n, n, 0.0);
-  u.resize (n, n, 0.0);
-  p.resize (n, n, 0.0);
-
-  for (int i = 0; i < n; i++)
-    {
-      p.elem (i, pvt[i]) = 1.0;
+  delete [] ipvt;
+  ipvt = 0;
 
-      l.elem (i, i) = 1.0;
-      for (int j = 0; j < i; j++)
-	l.elem (i, j) = A_fact.elem (i, j);
-
-      for (int j = i; j < n; j++)
-	u.elem (i, j) = A_fact.elem (i, j);
-    }
-
-  delete [] ipvt;
   delete [] pvt;
+  pvt = 0;
 }
 
 /*
--- a/liboctave/CmplxLU.h
+++ b/liboctave/CmplxLU.h
@@ -38,7 +38,7 @@
 {
 public:
 
-  ComplexLU (void) : l (), u (), p () { }
+  ComplexLU (void) : l (), u (), p (), ipvt (0), pvt (0) { }
 
   ComplexLU (const ComplexMatrix& a);
 
@@ -55,6 +55,12 @@
       return *this;
     }
 
+  ~ComplexLU (void)
+    {
+      delete [] ipvt;
+      delete [] pvt;
+    }
+
   ComplexMatrix L (void) const { return l; }
   ComplexMatrix U (void) const { return u; }
 
@@ -67,6 +73,9 @@
   ComplexMatrix l;
   ComplexMatrix u;
   Matrix p;
+
+  int *ipvt;
+  int *pvt;
 };
 
 #endif
--- a/liboctave/dbleLU.cc
+++ b/liboctave/dbleLU.cc
@@ -43,8 +43,12 @@
 
 LU::LU (const Matrix& a)
 {
+  ipvt = 0;
+  pvt = 0;
+
   int a_nr = a.rows ();
   int a_nc = a.cols ();
+
   if (a_nr == 0 || a_nc == 0 || a_nr != a_nc)
     {
       (*current_liboctave_error_handler) ("LU requires square matrix");
@@ -53,51 +57,60 @@
 
   int n = a_nr;
 
-  int *ipvt = new int [n];
-  int *pvt = new int [n];
-  double *tmp_data = dup (a.data (), a.length ());
+  ipvt = new int [n];
+  pvt = new int [n];
+
+  Matrix A_fact = a;
+  double *tmp_data = A_fact.fortran_vec ();
+
   int info = 0;
   double dummy = 0;
 
-  F77_FCN (dgesv, DGESV) (n, 0, tmp_data, n, ipvt, dummy, n, info);
+  F77_XFCN (dgesv, DGESV) (n, 0, tmp_data, n, ipvt, dummy, n, info);
 
-  Matrix A_fact (tmp_data, n, n);
-
-  for (int i = 0; i < n; i++)
+  if (f77_exception_encountered)
+    (*current_liboctave_error_handler) ("unrecoverable error in dgesv");
+  else
     {
-      ipvt[i] -= 1;
-      pvt[i] = i;
-    }
+      for (int i = 0; i < n; i++)
+	{
+	  ipvt[i] -= 1;
+	  pvt[i] = i;
+	}
 
-  for (int i = 0; i < n - 1; i++)
-    {
-      int k = ipvt[i];
-      if (k != i)
+      for (int i = 0; i < n - 1; i++)
 	{
-	  int tmp = pvt[k];
-	  pvt[k] = pvt[i];
-	  pvt[i] = tmp;
+	  int k = ipvt[i];
+	  if (k != i)
+	    {
+	      int tmp = pvt[k];
+	      pvt[k] = pvt[i];
+	      pvt[i] = tmp;
+	    }
+	}
+
+      l.resize (n, n, 0.0);
+      u.resize (n, n, 0.0);
+      p.resize (n, n, 0.0);
+
+      for (int i = 0; i < n; i++)
+	{
+	  p.elem (i, pvt[i]) = 1.0;
+
+	  l.elem (i, i) = 1.0;
+	  for (int j = 0; j < i; j++)
+	    l.elem (i, j) = A_fact.elem (i, j);
+
+	  for (int j = i; j < n; j++)
+	    u.elem (i, j) = A_fact.elem (i, j);
 	}
     }
 
-  l.resize (n, n, 0.0);
-  u.resize (n, n, 0.0);
-  p.resize (n, n, 0.0);
-
-  for (int i = 0; i < n; i++)
-    {
-      p.elem (i, pvt[i]) = 1.0;
+  delete [] ipvt;
+  ipvt = 0;
 
-      l.elem (i, i) = 1.0;
-      for (int j = 0; j < i; j++)
-	l.elem (i, j) = A_fact.elem (i, j);
-
-      for (int j = i; j < n; j++)
-	u.elem (i, j) = A_fact.elem (i, j);
-    }
-
-  delete [] ipvt;
   delete [] pvt;
+  pvt = 0;
 }
 
 /*
--- a/liboctave/dbleLU.h
+++ b/liboctave/dbleLU.h
@@ -36,7 +36,7 @@
 {
 public:
 
-  LU (void) : l (), u (), p () { }
+  LU (void) : l (), u (), p (), ipvt (0), pvt (0) { }
 
   LU (const Matrix& a);
 
@@ -53,6 +53,12 @@
       return *this;
     }
 
+  ~LU (void)
+    {
+      delete [] ipvt;
+      delete [] pvt;
+    }
+
   Matrix L (void) const { return l; }
 
   Matrix U (void) const { return u; }
@@ -66,6 +72,9 @@
   Matrix l;
   Matrix u;
   Matrix p;
+
+  int *ipvt;
+  int *pvt;
 };
 
 #endif