# HG changeset patch
# User PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
# Date 1408889106 -19800
# Node ID 5ce959c55cc04638c2c8ee43b30c53199df27c1f
# Parent  0fe7133da8ceb02744b3288841927f586b9a514c
Propagate 'lower' in chol(a, 'lower') to underlying library function.

* chol.cc (chol): Send 'L' parameter correctly when chol is called with 'lower'.
* floatCHOL.cc (init): Propagate 'lower' to underlying library function.
* floatCHOL.h: Modify the prototype of methods.
* fMatrix.cc (inverse): Invoke chol with additional parameter.
* dbleCHOL.cc (init): Propagate 'lower' to underlying library function.
* dbleCHOL.h: Modify the prototype of methods.
* dMatrix.cc (inverse): Invoke chol with additional parameter.
* CmplxCHOL.cc (init): Propagate 'lower' to underlying library function.
* CmplxCHOL.h: Modify the prototype of methods.
* CMatrix.cc (inverse): Invoke chol with additional parameter.

diff -r 0fe7133da8ce -r 5ce959c55cc0 libinterp/dldfcn/chol.cc
--- a/libinterp/dldfcn/chol.cc	Thu Aug 20 14:37:57 2015 -0400
+++ b/libinterp/dldfcn/chol.cc	Sun Aug 24 19:35:06 2014 +0530
@@ -46,6 +46,13 @@
 
 template <class CHOLT>
 static octave_value
+get_chol (const CHOLT& fact)
+{
+  return octave_value (fact.chol_matrix());
+}
+
+template <class CHOLT>
+static octave_value
 get_chol_r (const CHOLT& fact)
 {
   return octave_value (fact.chol_matrix (),
@@ -162,10 +169,6 @@
           if (tmp.compare ("vector") == 0)
             vecout = true;
           else if (tmp.compare ("lower") == 0)
-            // FIXME: currently the option "lower" is handled by transposing
-            //  the matrix, factorizing it with the lapack function
-            //  DPOTRF ('U', ...) and finally transposing the factor.  It would
-            //  be more efficient to use DPOTRF ('L', ...) in this case.
             LLt = true;
           else if (tmp.compare ("upper") == 0)
             LLt = false;
@@ -266,18 +269,12 @@
                   octave_idx_type info;
 
                   FloatCHOL fact;
-                  if (LLt)
-                    fact = FloatCHOL (m.transpose (), info);
-                  else
-                    fact = FloatCHOL (m, info);
+                  fact = FloatCHOL (m, info, LLt != true);
 
                   if (nargout == 2 || info == 0)
                     {
                       retval(1) = info;
-                      if (LLt)
-                        retval(0) = get_chol_l (fact);
-                      else
-                        retval(0) = get_chol_r (fact);
+                      retval(0) = get_chol (fact);
                     }
                   else
                     error ("chol: input matrix must be positive definite");
@@ -323,18 +320,12 @@
                   octave_idx_type info;
 
                   CHOL fact;
-                  if (LLt)
-                    fact = CHOL (m.transpose (), info);
-                  else
-                    fact = CHOL (m, info);
+                  fact = CHOL (m, info, LLt != true);
 
                   if (nargout == 2 || info == 0)
                     {
                       retval(1) = info;
-                      if (LLt)
-                        retval(0) = get_chol_l (fact);
-                      else
-                        retval(0) = get_chol_r (fact);
+                      retval(0) = get_chol (fact);
                     }
                   else
                     error ("chol: input matrix must be positive definite");
@@ -349,18 +340,12 @@
                   octave_idx_type info;
 
                   ComplexCHOL fact;
-                  if (LLt)
-                    fact = ComplexCHOL (m.transpose (), info);
-                  else
-                    fact = ComplexCHOL (m, info);
+                  fact = ComplexCHOL (m, info, LLt != true);
 
                   if (nargout == 2 || info == 0)
                     {
                       retval(1) = info;
-                      if (LLt)
-                        retval(0) = get_chol_l (fact);
-                      else
-                        retval(0) = get_chol_r (fact);
+                      retval(0) = get_chol (fact);
                     }
                   else
                     error ("chol: input matrix must be positive definite");
diff -r 0fe7133da8ce -r 5ce959c55cc0 liboctave/array/CMatrix.cc
--- a/liboctave/array/CMatrix.cc	Thu Aug 20 14:37:57 2015 -0400
+++ b/liboctave/array/CMatrix.cc	Sun Aug 24 19:35:06 2014 +0530
@@ -1193,7 +1193,7 @@
     {
       if (mattype.is_hermitian ())
         {
-          ComplexCHOL chol (*this, info, calc_cond);
+          ComplexCHOL chol (*this, info, true, calc_cond);
           if (info == 0)
             {
               if (calc_cond)
diff -r 0fe7133da8ce -r 5ce959c55cc0 liboctave/array/dMatrix.cc
--- a/liboctave/array/dMatrix.cc	Thu Aug 20 14:37:57 2015 -0400
+++ b/liboctave/array/dMatrix.cc	Sun Aug 24 19:35:06 2014 +0530
@@ -843,7 +843,7 @@
     {
       if (mattype.is_hermitian ())
         {
-          CHOL chol (*this, info, calc_cond);
+          CHOL chol (*this, info, true, calc_cond);
           if (info == 0)
             {
               if (calc_cond)
diff -r 0fe7133da8ce -r 5ce959c55cc0 liboctave/array/fMatrix.cc
--- a/liboctave/array/fMatrix.cc	Thu Aug 20 14:37:57 2015 -0400
+++ b/liboctave/array/fMatrix.cc	Sun Aug 24 19:35:06 2014 +0530
@@ -850,7 +850,7 @@
     {
       if (mattype.is_hermitian ())
         {
-          FloatCHOL chol (*this, info, calc_cond);
+          FloatCHOL chol (*this, info, true, calc_cond);
           if (info == 0)
             {
               if (calc_cond)
diff -r 0fe7133da8ce -r 5ce959c55cc0 liboctave/numeric/CmplxCHOL.cc
--- a/liboctave/numeric/CmplxCHOL.cc	Thu Aug 20 14:37:57 2015 -0400
+++ b/liboctave/numeric/CmplxCHOL.cc	Sun Aug 24 19:35:06 2014 +0530
@@ -86,7 +86,7 @@
 }
 
 octave_idx_type
-ComplexCHOL::init (const ComplexMatrix& a, bool calc_cond)
+ComplexCHOL::init (const ComplexMatrix& a, bool upper, bool calc_cond)
 {
   octave_idx_type a_nr = a.rows ();
   octave_idx_type a_nc = a.cols ();
@@ -101,13 +101,28 @@
   octave_idx_type n = a_nc;
   octave_idx_type info;
 
+  is_upper = upper;
+
   chol_mat.clear (n, n);
-  for (octave_idx_type j = 0; j < n; j++)
+  if (is_upper)
     {
-      for (octave_idx_type i = 0; i <= j; i++)
-        chol_mat.xelem (i, j) = a(i, j);
-      for (octave_idx_type i = j+1; i < n; i++)
-        chol_mat.xelem (i, j) = 0.0;
+      for (octave_idx_type j = 0; j < n; j++)
+        {
+          for (octave_idx_type i = 0; i <= j; i++)
+            chol_mat.xelem (i, j) = a (i, j);
+          for (octave_idx_type i = j + 1; i < n; i++)
+            chol_mat.xelem (i, j) = 0.0;
+        }
+     }
+  else
+    {
+      for (octave_idx_type j = 0; j < n; j++)
+        {
+          for (octave_idx_type i = 0; i < j; i++)
+            chol_mat.xelem (i, j) = 0.0;
+       	  for (octave_idx_type i = j; i < n; i++)
+            chol_mat.xelem (i, j) = a (i, j);
+        }
     }
   Complex *h = chol_mat.fortran_vec ();
 
@@ -116,8 +131,18 @@
   if (calc_cond)
     anorm = xnorm (a, 1);
 
-  F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1), n, h, n, info
-                             F77_CHAR_ARG_LEN (1)));
+  if (is_upper)
+    {
+      F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1),
+                                 n, h, n, info
+                                 F77_CHAR_ARG_LEN (1)));
+    }
+  else
+    {
+      F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1),
+                                 n, h, n, info
+                                 F77_CHAR_ARG_LEN (1)));
+    }
 
   xrcond = 0.0;
   if (info > 0)
@@ -143,7 +168,7 @@
 }
 
 static ComplexMatrix
-chol2inv_internal (const ComplexMatrix& r)
+chol2inv_internal (const ComplexMatrix& r, bool is_upper = true)
 {
   ComplexMatrix retval;
 
@@ -157,17 +182,37 @@
 
       ComplexMatrix tmp = r;
 
-      F77_XFCN (zpotri, ZPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n,
-                                 tmp.fortran_vec (), n, info
-                                 F77_CHAR_ARG_LEN (1)));
+      if (is_upper)
+        {
+          F77_XFCN (zpotri, ZPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n,
+                                     tmp.fortran_vec (), n, info
+                                     F77_CHAR_ARG_LEN (1)));
+        }
+      else
+        {
+          F77_XFCN (zpotri, ZPOTRI, (F77_CONST_CHAR_ARG2 ("L", 1), n,
+                                     tmp.fortran_vec (), n, info
+                                     F77_CHAR_ARG_LEN (1)));
+        }
 
       // If someone thinks of a more graceful way of doing this (or
       // faster for that matter :-)), please let me know!
 
       if (n > 1)
-        for (octave_idx_type j = 0; j < r_nc; j++)
-          for (octave_idx_type i = j+1; i < r_nr; i++)
-            tmp.xelem (i, j) = std::conj (tmp.xelem (j, i));
+        {
+          if (is_upper)
+            {
+              for (octave_idx_type j = 0; j < r_nc; j++)
+                for (octave_idx_type i = j+1; i < r_nr; i++)
+                  tmp.xelem (i, j) = tmp.xelem (j, i);
+            }
+          else
+            {
+              for (octave_idx_type j = 0; j < r_nc; j++)
+                for (octave_idx_type i = j+1; i < r_nr; i++)
+                  tmp.xelem (j, i) = tmp.xelem (i, j);
+            }
+        }
 
       retval = tmp;
     }
@@ -181,7 +226,7 @@
 ComplexMatrix
 ComplexCHOL::inverse (void) const
 {
-  return chol2inv_internal (chol_mat);
+  return chol2inv_internal (chol_mat, is_upper);
 }
 
 void
diff -r 0fe7133da8ce -r 5ce959c55cc0 liboctave/numeric/CmplxCHOL.h
--- a/liboctave/numeric/CmplxCHOL.h	Thu Aug 20 14:37:57 2015 -0400
+++ b/liboctave/numeric/CmplxCHOL.h	Sun Aug 24 19:35:06 2014 +0530
@@ -37,17 +37,17 @@
 
   ComplexCHOL (void) : chol_mat (), xrcond (0) { }
 
-  ComplexCHOL (const ComplexMatrix& a, bool calc_cond = false)
+  ComplexCHOL (const ComplexMatrix& a, bool upper = true, bool calc_cond = false)
     : chol_mat (), xrcond (0)
   {
-    init (a, calc_cond);
+    init (a, upper, calc_cond);
   }
 
-  ComplexCHOL (const ComplexMatrix& a, octave_idx_type& info,
+  ComplexCHOL (const ComplexMatrix& a, octave_idx_type& info, bool upper = true,
                bool calc_cond = false)
     : chol_mat (), xrcond (0)
   {
-    info = init (a, calc_cond);
+    info = init (a, upper, calc_cond);
   }
 
   ComplexCHOL (const ComplexCHOL& a)
@@ -91,7 +91,9 @@
 
   double xrcond;
 
-  octave_idx_type init (const ComplexMatrix& a, bool calc_cond);
+  bool is_upper;
+
+  octave_idx_type init (const ComplexMatrix& a, bool upper, bool calc_cond);
 };
 
 ComplexMatrix OCTAVE_API chol2inv (const ComplexMatrix& r);
diff -r 0fe7133da8ce -r 5ce959c55cc0 liboctave/numeric/dbleCHOL.cc
--- a/liboctave/numeric/dbleCHOL.cc	Thu Aug 20 14:37:57 2015 -0400
+++ b/liboctave/numeric/dbleCHOL.cc	Sun Aug 24 19:35:06 2014 +0530
@@ -87,7 +87,7 @@
 }
 
 octave_idx_type
-CHOL::init (const Matrix& a, bool calc_cond)
+CHOL::init (const Matrix& a, bool upper, bool calc_cond)
 {
   octave_idx_type a_nr = a.rows ();
   octave_idx_type a_nc = a.cols ();
@@ -101,13 +101,28 @@
   octave_idx_type n = a_nc;
   octave_idx_type info;
 
+  is_upper = upper;
+
   chol_mat.clear (n, n);
-  for (octave_idx_type j = 0; j < n; j++)
+  if (is_upper)
     {
-      for (octave_idx_type i = 0; i <= j; i++)
-        chol_mat.xelem (i, j) = a(i, j);
-      for (octave_idx_type i = j+1; i < n; i++)
-        chol_mat.xelem (i, j) = 0.0;
+      for (octave_idx_type j = 0; j < n; j++)
+        {
+          for (octave_idx_type i = 0; i <= j; i++)
+            chol_mat.xelem (i, j) = a (i, j);
+          for (octave_idx_type i = j + 1; i < n; i++)
+            chol_mat.xelem (i, j) = 0.0;
+        }
+     }
+  else
+    {
+      for (octave_idx_type j = 0; j < n; j++)
+        {
+          for (octave_idx_type i = 0; i < j; i++)
+            chol_mat.xelem (i, j) = 0.0;
+       	  for (octave_idx_type i = j; i < n; i++)
+            chol_mat.xelem (i, j) = a (i, j);
+        }
     }
   double *h = chol_mat.fortran_vec ();
 
@@ -116,9 +131,18 @@
   if (calc_cond)
     anorm = xnorm (a, 1);
 
-  F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1),
-                             n, h, n, info
-                             F77_CHAR_ARG_LEN (1)));
+  if (is_upper)
+    {
+      F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1),
+                                 n, h, n, info
+                                 F77_CHAR_ARG_LEN (1)));
+    }
+  else
+    {
+      F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1),
+                                 n, h, n, info
+                                 F77_CHAR_ARG_LEN (1)));
+    }
 
   xrcond = 0.0;
   if (info > 0)
@@ -132,9 +156,19 @@
       double *pz = z.fortran_vec ();
       Array<octave_idx_type> iz (dim_vector (n, 1));
       octave_idx_type *piz = iz.fortran_vec ();
-      F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h,
-                                 n, anorm, xrcond, pz, piz, dpocon_info
-                                 F77_CHAR_ARG_LEN (1)));
+      if (is_upper)
+        {
+          F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h,
+                                     n, anorm, xrcond, pz, piz, dpocon_info
+                                     F77_CHAR_ARG_LEN (1)));
+	}
+      else
+        {
+          F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 ("L", 1), n, h,
+                                     n, anorm, xrcond, pz, piz, dpocon_info
+                                     F77_CHAR_ARG_LEN (1)));
+	}
+      
 
       if (dpocon_info != 0)
         info = -1;
@@ -144,7 +178,7 @@
 }
 
 static Matrix
-chol2inv_internal (const Matrix& r)
+chol2inv_internal (const Matrix& r, bool is_upper = true)
 {
   Matrix retval;
 
@@ -161,17 +195,37 @@
 
       if (info == 0)
         {
-          F77_XFCN (dpotri, DPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n,
-                                     v, n, info
-                                     F77_CHAR_ARG_LEN (1)));
+          if (is_upper)
+            {
+              F77_XFCN (dpotri, DPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n,
+                                         v, n, info
+                                         F77_CHAR_ARG_LEN (1)));
+            }
+          else
+            {
+              F77_XFCN (dpotri, DPOTRI, (F77_CONST_CHAR_ARG2 ("L", 1), n,
+                                         v, n, info
+                                         F77_CHAR_ARG_LEN (1)));
+            }
 
           // If someone thinks of a more graceful way of doing this (or
           // faster for that matter :-)), please let me know!
 
           if (n > 1)
-            for (octave_idx_type j = 0; j < r_nc; j++)
-              for (octave_idx_type i = j+1; i < r_nr; i++)
-                tmp.xelem (i, j) = tmp.xelem (j, i);
+            {
+              if (is_upper)
+                {
+                  for (octave_idx_type j = 0; j < r_nc; j++)
+                    for (octave_idx_type i = j+1; i < r_nr; i++)
+                      tmp.xelem (i, j) = tmp.xelem (j, i);
+                }
+              else
+                {
+                  for (octave_idx_type j = 0; j < r_nc; j++)
+                    for (octave_idx_type i = j+1; i < r_nr; i++)
+                      tmp.xelem (j, i) = tmp.xelem (i, j);
+                }
+            }
 
           retval = tmp;
         }
@@ -186,7 +240,7 @@
 Matrix
 CHOL::inverse (void) const
 {
-  return chol2inv_internal (chol_mat);
+  return chol2inv_internal (chol_mat, is_upper);
 }
 
 void
diff -r 0fe7133da8ce -r 5ce959c55cc0 liboctave/numeric/dbleCHOL.h
--- a/liboctave/numeric/dbleCHOL.h	Thu Aug 20 14:37:57 2015 -0400
+++ b/liboctave/numeric/dbleCHOL.h	Sun Aug 24 19:35:06 2014 +0530
@@ -37,16 +37,16 @@
 
   CHOL (void) : chol_mat (), xrcond (0) { }
 
-  CHOL (const Matrix& a, bool calc_cond = false)
+  CHOL (const Matrix& a, bool upper = true, bool calc_cond = false)
     : chol_mat (), xrcond (0)
   {
-    init (a, calc_cond);
+    init (a, upper, calc_cond);
   }
 
-  CHOL (const Matrix& a, octave_idx_type& info, bool calc_cond = false)
-    : chol_mat (), xrcond (0)
+  CHOL (const Matrix& a, octave_idx_type& info, bool upper = true,
+        bool calc_cond = false) : chol_mat (), xrcond (0)
   {
-    info = init (a, calc_cond);
+    info = init (a, upper, calc_cond);
   }
 
   CHOL (const CHOL& a) : chol_mat (a.chol_mat), xrcond (a.xrcond) { }
@@ -88,7 +88,9 @@
 
   double xrcond;
 
-  octave_idx_type init (const Matrix& a, bool calc_cond);
+  bool is_upper;
+
+  octave_idx_type init (const Matrix& a, bool upper, bool calc_cond);
 };
 
 Matrix OCTAVE_API chol2inv (const Matrix& r);
diff -r 0fe7133da8ce -r 5ce959c55cc0 liboctave/numeric/floatCHOL.cc
--- a/liboctave/numeric/floatCHOL.cc	Thu Aug 20 14:37:57 2015 -0400
+++ b/liboctave/numeric/floatCHOL.cc	Sun Aug 24 19:35:06 2014 +0530
@@ -87,7 +87,7 @@
 }
 
 octave_idx_type
-FloatCHOL::init (const FloatMatrix& a, bool calc_cond)
+FloatCHOL::init (const FloatMatrix& a, bool upper, bool calc_cond)
 {
   octave_idx_type a_nr = a.rows ();
   octave_idx_type a_nc = a.cols ();
@@ -101,14 +101,30 @@
   octave_idx_type n = a_nc;
   octave_idx_type info;
 
+  is_upper = upper;
+
   chol_mat.clear (n, n);
-  for (octave_idx_type j = 0; j < n; j++)
+  if (is_upper)
     {
-      for (octave_idx_type i = 0; i <= j; i++)
-        chol_mat.xelem (i, j) = a(i, j);
-      for (octave_idx_type i = j+1; i < n; i++)
-        chol_mat.xelem (i, j) = 0.0f;
+      for (octave_idx_type j = 0; j < n; j++)
+        {
+          for (octave_idx_type i = 0; i <= j; i++)
+            chol_mat.xelem (i, j) = a(i, j);
+          for (octave_idx_type i = j+1; i < n; i++)
+            chol_mat.xelem (i, j) = 0.0f;
+        }
     }
+  else
+    {
+      for (octave_idx_type j = 0; j < n; j++)
+        {
+          for (octave_idx_type i = 0; i <= j; i++)
+            chol_mat.xelem (i, j) = 0.0f;
+          for (octave_idx_type i = j+1; i < n; i++)
+            chol_mat.xelem (i, j) = a(i, j);
+        }
+    }
+
   float *h = chol_mat.fortran_vec ();
 
   // Calculate the norm of the matrix, for later use.
@@ -116,9 +132,18 @@
   if (calc_cond)
     anorm = xnorm (a, 1);
 
-  F77_XFCN (spotrf, SPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1),
-                             n, h, n, info
-                             F77_CHAR_ARG_LEN (1)));
+  if (is_upper)
+    {
+      F77_XFCN (spotrf, SPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1),
+                                 n, h, n, info
+                                 F77_CHAR_ARG_LEN (1)));   
+    }
+  else
+    {
+      F77_XFCN (spotrf, SPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1),
+                                 n, h, n, info
+                                 F77_CHAR_ARG_LEN (1)));   
+    }
 
   xrcond = 0.0;
   if (info > 0)
@@ -132,9 +157,19 @@
       float *pz = z.fortran_vec ();
       Array<octave_idx_type> iz (dim_vector (n, 1));
       octave_idx_type *piz = iz.fortran_vec ();
-      F77_XFCN (spocon, SPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h,
-                                 n, anorm, xrcond, pz, piz, spocon_info
-                                 F77_CHAR_ARG_LEN (1)));
+      if (is_upper)
+        {
+          F77_XFCN (spocon, SPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h,
+                                     n, anorm, xrcond, pz, piz, spocon_info
+                                     F77_CHAR_ARG_LEN (1)));       
+        }
+      else
+        {
+          F77_XFCN (spocon, SPOCON, (F77_CONST_CHAR_ARG2 ("L", 1), n, h,
+                                     n, anorm, xrcond, pz, piz, spocon_info
+                                     F77_CHAR_ARG_LEN (1)));       
+        }
+
 
       if (spocon_info != 0)
         info = -1;
@@ -144,7 +179,7 @@
 }
 
 static FloatMatrix
-chol2inv_internal (const FloatMatrix& r)
+chol2inv_internal (const FloatMatrix& r, bool is_upper = true)
 {
   FloatMatrix retval;
 
@@ -161,17 +196,37 @@
 
       if (info == 0)
         {
-          F77_XFCN (spotri, SPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n,
-                                     v, n, info
-                                     F77_CHAR_ARG_LEN (1)));
+          if (is_upper)
+            {
+              F77_XFCN (spotri, SPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n,
+                                         v, n, info
+                                         F77_CHAR_ARG_LEN (1)));
+            }
+          else
+            {
+              F77_XFCN (spotri, SPOTRI, (F77_CONST_CHAR_ARG2 ("L", 1), n,
+                                         v, n, info
+                                         F77_CHAR_ARG_LEN (1)));
+            }
 
           // If someone thinks of a more graceful way of doing this (or
           // faster for that matter :-)), please let me know!
 
           if (n > 1)
-            for (octave_idx_type j = 0; j < r_nc; j++)
-              for (octave_idx_type i = j+1; i < r_nr; i++)
-                tmp.xelem (i, j) = tmp.xelem (j, i);
+            {
+              if (is_upper)
+                {
+                  for (octave_idx_type j = 0; j < r_nc; j++)
+                    for (octave_idx_type i = j+1; i < r_nr; i++)
+                      tmp.xelem (i, j) = tmp.xelem (j, i); 
+                }
+              else
+                {
+                  for (octave_idx_type j = 0; j < r_nc; j++)
+                    for (octave_idx_type i = j+1; i < r_nr; i++)
+                      tmp.xelem (j, i) = tmp.xelem (i, j);
+                }
+            }
 
           retval = tmp;
         }
@@ -186,7 +241,7 @@
 FloatMatrix
 FloatCHOL::inverse (void) const
 {
-  return chol2inv_internal (chol_mat);
+  return chol2inv_internal (chol_mat, is_upper);
 }
 
 void
diff -r 0fe7133da8ce -r 5ce959c55cc0 liboctave/numeric/floatCHOL.h
--- a/liboctave/numeric/floatCHOL.h	Thu Aug 20 14:37:57 2015 -0400
+++ b/liboctave/numeric/floatCHOL.h	Sun Aug 24 19:35:06 2014 +0530
@@ -37,17 +37,17 @@
 
   FloatCHOL (void) : chol_mat (), xrcond (0) { }
 
-  FloatCHOL (const FloatMatrix& a, bool calc_cond = false)
+  FloatCHOL (const FloatMatrix& a, bool upper = true, bool calc_cond = false)
     : chol_mat (), xrcond (0)
   {
-    init (a, calc_cond);
+    init (a, upper, calc_cond);
   }
 
-  FloatCHOL (const FloatMatrix& a, octave_idx_type& info,
+  FloatCHOL (const FloatMatrix& a, octave_idx_type& info, bool upper = true,
              bool calc_cond = false)
     : chol_mat (), xrcond (0)
   {
-    info = init (a, calc_cond);
+    info = init (a, upper, calc_cond);
   }
 
   FloatCHOL (const FloatCHOL& a) : chol_mat (a.chol_mat), xrcond (a.xrcond) { }
@@ -90,7 +90,9 @@
 
   float xrcond;
 
-  octave_idx_type init (const FloatMatrix& a, bool calc_cond);
+  bool is_upper;
+
+  octave_idx_type init (const FloatMatrix& a, bool upper, bool calc_cond);
 };
 
 FloatMatrix OCTAVE_API chol2inv (const FloatMatrix& r);