--- 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;
 }
 
 /*