# HG changeset patch
# User Júlio Hoffimann <julio.hoffimann@gmail.com>
# Date 1362951535 10800
# Node ID 9a78540d8b38128ea9fd3646a0efdad95fa96340
# Parent  c31fd42f9000493ecc36050485ae7b98615b9992
Added trivial OpenMP constructs to liboctave/array classes as an experiment.

Affected files: Array.cc, CColVector.cc, CDiagMatrix.cc, CMatrix.cc,
CNDArray.cc, CRowVector.cc, CSparse.cc, MSparse.cc, Range.cc, Sparse.cc,
boolSparse.cc, dColVector.cc, dDiagMatrix.cc, dMatrix.cc, dNDArray.cc,
dRowVector.cc, dSparse.cc, fCColVector.cc, fCDiagMatrix.c, fCMatrix.cc,
fCNDArray.cc, fCRowVector.cc, fColVector.cc, fDiagMatrix.cc, fMatrix.cc,
fRowVector.cc, idx-vector.cc

diff -r c31fd42f9000 -r 9a78540d8b38 liboctave/array/Array.cc
--- a/liboctave/array/Array.cc	Sat Mar 09 18:42:43 2013 +0000
+++ b/liboctave/array/Array.cc	Sun Mar 10 18:38:55 2013 -0300
@@ -1427,8 +1427,15 @@
           octave_idx_type nd = n + l - u, dl = 1, du = 1;
           dim_vector rdv = dimensions;
           rdv(dim) = nd;
-          for (int k = 0; k < dim; k++) dl *= dimensions(k);
-          for (int k = dim + 1; k < ndims (); k++) du *= dimensions(k);
+          #pragma omp parallel
+          {
+            #pragma omp for reduction(*:dl)
+            for (int k = 0; k < dim; k++)
+              dl *= dimensions(k);
+            #pragma omp for reduction(*:du)
+            for (int k = dim + 1; k < ndims (); k++)
+              du *= dimensions(k);
+          }
 
           // Special case deleting a contiguous range.
           Array<T> tmp = Array<T> (rdv);
@@ -1715,6 +1722,7 @@
   octave_idx_type iter = dv.numel () / ns;
   octave_idx_type stride = 1;
 
+  #pragma omp parallel for reduction(*:stride)
   for (int i = 0; i < dim; i++)
     stride *= dv(i);
 
@@ -1834,6 +1842,7 @@
   octave_idx_type iter = dv.numel () / ns;
   octave_idx_type stride = 1;
 
+  #pragma omp parallel for reduction(*:stride)
   for (int i = 0; i < dim; i++)
     stride *= dv(i);
 
@@ -2173,6 +2182,7 @@
       // We want all elements, which means we'll almost surely need
       // to resize. So count first, then allocate array of exact size.
       octave_idx_type cnt = 0;
+      #pragma omp parallel for reduction(+:cnt)
       for (octave_idx_type i = 0; i < nel; i++)
         cnt += src[i] != zero;
 
@@ -2313,6 +2323,7 @@
   octave_idx_type iter = numel () / ns;
   octave_idx_type stride = 1;
 
+  #pragma omp parallel for reduction(*:stride)
   for (int i = 0; i < dim; i++)
     stride *= dv(i);
 
diff -r c31fd42f9000 -r 9a78540d8b38 liboctave/array/CColVector.cc
--- a/liboctave/array/CColVector.cc	Sat Mar 09 18:42:43 2013 +0000
+++ b/liboctave/array/CColVector.cc	Sun Mar 10 18:38:55 2013 -0300
@@ -394,6 +394,7 @@
 
   ComplexColumnVector result (nr);
 
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < a_len; i++)
     result.elem (i) = a.elem (i) * m.elem (i, i);
 
@@ -422,6 +423,7 @@
 
   ComplexColumnVector result (nr);
 
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < a_len; i++)
     result.elem (i) = a.elem (i) * m.elem (i, i);
 
@@ -450,6 +452,7 @@
 
   ComplexColumnVector result (nr);
 
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < a_len; i++)
     result.elem (i) = a.elem (i) * m.elem (i, i);
 
@@ -506,9 +509,8 @@
 std::ostream&
 operator << (std::ostream& os, const ComplexColumnVector& a)
 {
-//  int field_width = os.precision () + 7;
   for (octave_idx_type i = 0; i < a.length (); i++)
-    os << /* setw (field_width) << */ a.elem (i) << "\n";
+    os << a.elem (i) << "\n";
   return os;
 }
 
diff -r c31fd42f9000 -r 9a78540d8b38 liboctave/array/CDiagMatrix.cc
--- a/liboctave/array/CDiagMatrix.cc	Sat Mar 09 18:42:43 2013 +0000
+++ b/liboctave/array/CDiagMatrix.cc	Sun Mar 10 18:38:55 2013 -0300
@@ -389,6 +389,7 @@
 
   ComplexDiagMatrix retval (c, r);
 
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < len; i++)
     {
       if (elem (i, i) != 0.0)
@@ -448,8 +449,10 @@
 
   octave_idx_type len = c.length (), lenm = len < a_nc ? len : a_nc;
 
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < lenm; i++)
     c.dgxelem (i) = a.dgelem (i) * b.dgelem (i);
+
   for (octave_idx_type i = lenm; i < len; i++)
     c.dgxelem (i) = 0.0;
 
@@ -478,6 +481,7 @@
 
   octave_idx_type len = a_nr < b_nc ? a_nr : b_nc;
 
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < len; i++)
     {
       double a_element = a.elem (i, i);
@@ -511,6 +515,7 @@
 
   octave_idx_type len = a_nr < b_nc ? a_nr : b_nc;
 
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < len; i++)
     {
       Complex a_element = a.elem (i, i);
@@ -557,17 +562,18 @@
 operator << (std::ostream& os, const ComplexDiagMatrix& a)
 {
   Complex ZERO (0.0);
-//  int field_width = os.precision () + 7;
+
   for (octave_idx_type i = 0; i < a.rows (); i++)
     {
       for (octave_idx_type j = 0; j < a.cols (); j++)
         {
           if (i == j)
-            os << " " /* setw (field_width) */ << a.elem (i, i);
+            os << " " << a.elem (i, i);
           else
-            os << " " /* setw (field_width) */ << ZERO;
+            os << " " << ZERO;
         }
       os << "\n";
     }
+
   return os;
 }
diff -r c31fd42f9000 -r 9a78540d8b38 liboctave/array/CMatrix.cc
--- a/liboctave/array/CMatrix.cc	Sat Mar 09 18:42:43 2013 +0000
+++ b/liboctave/array/CMatrix.cc	Sun Mar 10 18:38:55 2013 -0300
@@ -1399,6 +1399,7 @@
       F77_FUNC (zfftb, ZFFTB) (npts, &tmp_data[npts*j], pwsave);
     }
 
+  #pragma omp parallel for
   for (octave_idx_type j = 0; j < npts*nsamples; j++)
     tmp_data[j] = tmp_data[j] / static_cast<double> (npts);
 
@@ -1509,6 +1510,7 @@
       F77_FUNC (zfftb, ZFFTB) (npts, &tmp_data[npts*j], pwsave);
     }
 
+  #pragma omp parallel for
   for (octave_idx_type j = 0; j < npts*nsamples; j++)
     tmp_data[j] = tmp_data[j] / static_cast<double> (npts);
 
@@ -2958,6 +2960,7 @@
       return *this;
     }
 
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < a.length (); i++)
     elem (i, i) += a.elem (i, i);
 
@@ -2979,6 +2982,7 @@
       return *this;
     }
 
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < a.length (); i++)
     elem (i, i) -= a.elem (i, i);
 
@@ -3000,6 +3004,7 @@
       return *this;
     }
 
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < a.length (); i++)
     elem (i, i) += a.elem (i, i);
 
@@ -3021,6 +3026,7 @@
       return *this;
     }
 
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < a.length (); i++)
     elem (i, i) -= a.elem (i, i);
 
@@ -3781,6 +3787,7 @@
                                        a.data (), lda, 0.0, c, a_nr
                                        F77_CHAR_ARG_LEN (1)
                                        F77_CHAR_ARG_LEN (1)));
+              #pragma omp parallel for
               for (octave_idx_type j = 0; j < a_nr; j++)
                 for (octave_idx_type i = 0; i < j; i++)
                   retval.xelem (j,i) = std::conj (retval.xelem (i,j));
@@ -3793,6 +3800,7 @@
                                        a.data (), lda, 0.0, c, a_nr
                                        F77_CHAR_ARG_LEN (1)
                                        F77_CHAR_ARG_LEN (1)));
+              #pragma omp parallel for
               for (octave_idx_type j = 0; j < a_nr; j++)
                 for (octave_idx_type i = 0; i < j; i++)
                   retval.xelem (j,i) = retval.xelem (i,j);
@@ -4065,15 +4073,21 @@
 
   // The last column is not needed while using delta.
   Complex *delta = &retval(0, n-1);
-  for (octave_idx_type i = 0; i < m; i++)
-    delta[i] = (x2(i) - x1(i)) / (n - 1.0);
-
-  for (octave_idx_type j = 1; j < n-1; j++)
+  #pragma omp parallel
+  {
+    #pragma omp for
     for (octave_idx_type i = 0; i < m; i++)
-      retval(i, j) = x1(i) + static_cast<double> (j)*delta[i];
-
-  for (octave_idx_type i = 0; i < m; i++)
-    retval(i, n-1) = x2(i);
+      delta[i] = (x2(i) - x1(i)) / (n - 1.0);
+
+    #pragma omp for collapse(2)
+    for (octave_idx_type j = 1; j < n-1; j++)
+      for (octave_idx_type i = 0; i < m; i++)
+        retval(i, j) = x1(i) + static_cast<double> (j)*delta[i];
+
+    #pragma omp for
+    for (octave_idx_type i = 0; i < m; i++)
+      retval(i, n-1) = x2(i);
+  }
 
   return retval;
 }
diff -r c31fd42f9000 -r 9a78540d8b38 liboctave/array/CNDArray.cc
--- a/liboctave/array/CNDArray.cc	Sat Mar 09 18:42:43 2013 +0000
+++ b/liboctave/array/CNDArray.cc	Sun Mar 10 18:38:55 2013 -0300
@@ -65,6 +65,7 @@
   octave_idx_type stride = 1;
   octave_idx_type n = dv(dim);
 
+  #pragma omp parallel for reduction(*:stride)
   for (int i = 0; i < dim; i++)
     stride *= dv(i);
 
@@ -96,6 +97,7 @@
   octave_idx_type stride = 1;
   octave_idx_type n = dv(dim);
 
+  #pragma omp parallel for reduction(*:stride)
   for (int i = 0; i < dim; i++)
     stride *= dv(i);
 
@@ -223,6 +225,7 @@
 
   octave_idx_type stride = 1;
 
+  #pragma omp parallel for reduction(*:stride)
   for (int i = 0; i < dim; i++)
     stride *= dv(i);
 
@@ -270,6 +273,7 @@
 
   octave_idx_type stride = 1;
 
+  #pragma omp parallel for reduction(*:stride)
   for (int i = 0; i < dim; i++)
     stride *= dv(i);
 
@@ -776,7 +780,7 @@
 
       octave_idx_type n_elt = a.numel ();
 
-      // IS make_unique () NECCESSARY HERE??
+      // FIXME: is make_unique () neccessary here?
 
       for (octave_idx_type i = 0; i < n_elt; i++)
         {
diff -r c31fd42f9000 -r 9a78540d8b38 liboctave/array/CRowVector.cc
--- a/liboctave/array/CRowVector.cc	Sat Mar 09 18:42:43 2013 +0000
+++ b/liboctave/array/CRowVector.cc	Sun Mar 10 18:38:55 2013 -0300
@@ -402,9 +402,8 @@
 std::ostream&
 operator << (std::ostream& os, const ComplexRowVector& a)
 {
-//  int field_width = os.precision () + 7;
   for (octave_idx_type i = 0; i < a.length (); i++)
-    os << " " /* setw (field_width) */ << a.elem (i);
+    os << " " << a.elem (i);
   return os;
 }
 
@@ -467,8 +466,11 @@
 
   Complex delta = (x2 - x1) / (n - 1.0);
   retval(0) = x1;
+
+  #pragma omp parallel for
   for (octave_idx_type i = 1; i < n-1; i++)
     retval(i) = x1 + static_cast<double> (i)*delta;
+
   retval(n-1) = x2;
 
   return retval;
diff -r c31fd42f9000 -r 9a78540d8b38 liboctave/array/CSparse.cc
--- a/liboctave/array/CSparse.cc	Sat Mar 09 18:42:43 2013 +0000
+++ b/liboctave/array/CSparse.cc	Sun Mar 10 18:38:55 2013 -0300
@@ -150,6 +150,7 @@
   for (octave_idx_type i = 0; i < nc + 1; i++)
     cidx (i) = a.cidx (i);
 
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < nz; i++)
     {
       data (i) = Complex (a.data (i));
@@ -687,6 +688,7 @@
   for (octave_idx_type i = 0; i < nc + 1; i++)
     retval.cidx (i) = a.cidx (i);
 
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < nz; i++)
     {
       retval.data (i) = conj (a.data (i));
@@ -764,6 +766,7 @@
               rcond = dmin / dmax;
             }
 
+          #pragma omp parallel for
           for (octave_idx_type i = 0; i < nr; i++)
             v[i] = 1.0 / v[i];
         }
@@ -805,6 +808,7 @@
               for (octave_idx_type j = 0; j < nr; j++)
                 {
                   double atmp = 0.;
+                  #pragma omp parallel for reduction(+:atmp)
                   for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
                     atmp += std::abs (data (i));
                   if (atmp > anorm)
@@ -913,6 +917,7 @@
                     }
 
                   if (pivot != 1.0)
+                    #pragma omp parallel for
                     for (octave_idx_type j = cx_colstart; j < cx; j++)
                       retval.xdata (j) /= pivot;
                 }
@@ -1032,8 +1037,8 @@
               for (octave_idx_type j = 0; j < nr; j++)
                 {
                   double atmp = 0.;
-                  for (octave_idx_type i = retval.cidx (j);
-                       i < retval.cidx (j+1); i++)
+                  #pragma omp parallel for reduction(+:atmp)
+                  for (octave_idx_type i = retval.cidx (j); i < retval.cidx (j+1); i++)
                     atmp += std::abs (retval.data (i));
                   if (atmp > ainvnorm)
                     ainvnorm = atmp;
@@ -1281,10 +1286,12 @@
         {
           retval.resize (nc, b.cols (), Complex (0.,0.));
           if (typ == MatrixType::Diagonal)
+            #pragma omp parallel for collapse(2)
             for (octave_idx_type j = 0; j < b.cols (); j++)
                 for (octave_idx_type i = 0; i < nm; i++)
                   retval(i,j) = b(i,j) / data (i);
           else
+            #pragma omp parallel for collapse(2)
             for (octave_idx_type j = 0; j < b.cols (); j++)
               for (octave_idx_type k = 0; k < nc; k++)
                 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
@@ -1433,10 +1440,12 @@
         {
           retval.resize (nc, b.cols (), Complex (0.,0.));
           if (typ == MatrixType::Diagonal)
+            #pragma omp parallel for collapse(2)
             for (octave_idx_type j = 0; j < b.cols (); j++)
               for (octave_idx_type i = 0; i < nm; i++)
                 retval(i,j) = b(i,j) / data (i);
           else
+            #pragma omp parallel for collapse(2)
             for (octave_idx_type j = 0; j < b.cols (); j++)
               for (octave_idx_type k = 0; k < nc; k++)
                 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
@@ -1689,6 +1698,7 @@
                 {
                   for (octave_idx_type i = 0; i < nr; i++)
                     work[i] = b(i,j);
+
                   for (octave_idx_type i = nr; i < nc; i++)
                     work[i] = 0.;
 
@@ -1705,11 +1715,9 @@
 
                           Complex tmp = work[k] / data (cidx (k+1)-1);
                           work[k] = tmp;
+                          #pragma omp parallel for
                           for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
-                            {
-                              octave_idx_type iidx = ridx (i);
-                              work[iidx] = work[iidx] - tmp * data (i);
-                            }
+                            work[ridx (i)] -= tmp * data (i);
                         }
                     }
 
@@ -1742,6 +1750,7 @@
                             }
                         }
                       double atmp = 0;
+                      #pragma omp parallel for reduction(+:atmp)
                       for (octave_idx_type i = 0; i < j+1; i++)
                         {
                           atmp += std::abs (work[i]);
@@ -1829,6 +1838,7 @@
               for (octave_idx_type j = 0; j < nc; j++)
                 {
                   double atmp = 0.;
+                  #pragma omp parallel for reduction(+:atmp)
                   for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
                     atmp += std::abs (data (i));
                   if (atmp > anorm)
@@ -1874,12 +1884,9 @@
 
                           Complex tmp = work[k] / data (cidx (kidx+1)-1);
                           work[k] = tmp;
-                          for (octave_idx_type i = cidx (kidx);
-                               i < cidx (kidx+1)-1; i++)
-                            {
-                              octave_idx_type iidx = ridx (i);
-                              work[iidx] = work[iidx] - tmp * data (i);
-                            }
+                          #pragma omp parallel for
+                          for (octave_idx_type i = cidx (kidx); i < cidx (kidx+1)-1; i++)
+                            work[ridx (i)] -= tmp * data (i);
                         }
                     }
 
@@ -1936,6 +1943,7 @@
                             }
                         }
                       double atmp = 0;
+                      #pragma omp parallel for reduction(+:atmp)
                       for (octave_idx_type i = 0; i < j+1; i++)
                         {
                           atmp += std::abs (work[i]);
@@ -1971,11 +1979,9 @@
 
                           Complex tmp = work[k] / data (cidx (k+1)-1);
                           work[k] = tmp;
+                          #pragma omp parallel for
                           for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
-                            {
-                              octave_idx_type iidx = ridx (i);
-                              work[iidx] = work[iidx] - tmp * data (i);
-                            }
+                            work[ridx (i)] -= tmp * data (i);
                         }
                     }
 
@@ -2021,15 +2027,13 @@
                             {
                               Complex tmp = work[k] / data (cidx (k+1)-1);
                               work[k] = tmp;
-                              for (octave_idx_type i = cidx (k);
-                                   i < cidx (k+1)-1; i++)
-                                {
-                                  octave_idx_type iidx = ridx (i);
-                                  work[iidx] = work[iidx] - tmp * data (i);
-                                }
+                              #pragma omp parallel for
+                              for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
+                                work[ridx (i)] -= tmp * data (i);
                             }
                         }
                       double atmp = 0;
+                      #pragma omp parallel for reduction(+:atmp)
                       for (octave_idx_type i = 0; i < j+1; i++)
                         {
                           atmp += std::abs (work[i]);
@@ -2117,6 +2121,7 @@
               for (octave_idx_type j = 0; j < nc; j++)
                 {
                   double atmp = 0.;
+                  #pragma omp parallel for reduction(+:atmp)
                   for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
                     atmp += std::abs (data (i));
                   if (atmp > anorm)
@@ -2152,12 +2157,9 @@
 
                           Complex tmp = work[k] / data (cidx (kidx+1)-1);
                           work[k] = tmp;
-                          for (octave_idx_type i = cidx (kidx);
-                               i < cidx (kidx+1)-1; i++)
-                            {
-                              octave_idx_type iidx = ridx (i);
-                              work[iidx] = work[iidx] - tmp * data (i);
-                            }
+                          #pragma omp parallel for
+                          for (octave_idx_type i = cidx (kidx); i < cidx (kidx+1)-1; i++)
+                            work[ridx (i)] -= tmp * data (i);
                         }
                     }
 
@@ -2183,15 +2185,13 @@
                             {
                               Complex tmp = work[k] / data (cidx (iidx+1)-1);
                               work[k] = tmp;
-                              for (octave_idx_type i = cidx (iidx);
-                                   i < cidx (iidx+1)-1; i++)
-                                {
-                                  octave_idx_type idx2 = ridx (i);
-                                  work[idx2] = work[idx2] - tmp * data (i);
-                                }
+                              #pragma omp parallel for
+                              for (octave_idx_type i = cidx (iidx); i < cidx (iidx+1)-1; i++)
+                                work[ridx (i)] -= tmp * data (i);
                             }
                         }
                       double atmp = 0;
+                      #pragma omp parallel for reduction(+:atmp)
                       for (octave_idx_type i = 0; i < j+1; i++)
                         {
                           atmp += std::abs (work[i]);
@@ -2228,11 +2228,9 @@
 
                           Complex tmp = work[k] / data (cidx (k+1)-1);
                           work[k] = tmp;
+                          #pragma omp parallel for
                           for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
-                            {
-                              octave_idx_type iidx = ridx (i);
-                              work[iidx] = work[iidx] - tmp * data (i);
-                            }
+                            work[ridx (i)] -= tmp * data (i);
                         }
                     }
 
@@ -2256,15 +2254,13 @@
                             {
                               Complex tmp = work[k] / data (cidx (k+1)-1);
                               work[k] = tmp;
-                              for (octave_idx_type i = cidx (k);
-                                   i < cidx (k+1)-1; i++)
-                                {
-                                  octave_idx_type iidx = ridx (i);
-                                  work[iidx] = work[iidx] - tmp * data (i);
-                                }
+                              #pragma omp parallel for
+                              for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
+                                work[ridx (i)] -= tmp * data (i);
                             }
                         }
                       double atmp = 0;
+                      #pragma omp parallel for reduction(+:atmp)
                       for (octave_idx_type i = 0; i < j+1; i++)
                         {
                           atmp += std::abs (work[i]);
@@ -2352,6 +2348,7 @@
               for (octave_idx_type j = 0; j < nc; j++)
                 {
                   double atmp = 0.;
+                  #pragma omp parallel for reduction(+:atmp)
                   for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
                     atmp += std::abs (data (i));
                   if (atmp > anorm)
@@ -2397,12 +2394,9 @@
 
                           Complex tmp = work[k] / data (cidx (kidx+1)-1);
                           work[k] = tmp;
-                          for (octave_idx_type i = cidx (kidx);
-                               i < cidx (kidx+1)-1; i++)
-                            {
-                              octave_idx_type iidx = ridx (i);
-                              work[iidx] = work[iidx] - tmp * data (i);
-                            }
+                          #pragma omp parallel for
+                          for (octave_idx_type i = cidx (kidx); i < cidx (kidx+1)-1; i++)
+                            work[ridx (i)] -= tmp * data (i);
                         }
                     }
 
@@ -2450,15 +2444,13 @@
                             {
                               Complex tmp = work[k] / data (cidx (iidx+1)-1);
                               work[k] = tmp;
-                              for (octave_idx_type i = cidx (iidx);
-                                   i < cidx (iidx+1)-1; i++)
-                                {
-                                  octave_idx_type idx2 = ridx (i);
-                                  work[idx2] = work[idx2] - tmp * data (i);
-                                }
+                              #pragma omp parallel for
+                              for (octave_idx_type i = cidx (iidx); i < cidx (iidx+1)-1; i++)
+                                work[ridx (i)] -= tmp * data (i);
                             }
                         }
                       double atmp = 0;
+                      #pragma omp parallel for reduction(+:atmp)
                       for (octave_idx_type i = 0; i < j+1; i++)
                         {
                           atmp += std::abs (work[i]);
@@ -2494,11 +2486,9 @@
 
                           Complex tmp = work[k] / data (cidx (k+1)-1);
                           work[k] = tmp;
+                          #pragma omp parallel for
                           for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
-                            {
-                              octave_idx_type iidx = ridx (i);
-                              work[iidx] = work[iidx] - tmp * data (i);
-                            }
+                            work[ridx (i)] -= tmp * data (i);
                         }
                     }
 
@@ -2544,15 +2534,13 @@
                             {
                               Complex tmp = work[k] / data (cidx (k+1)-1);
                               work[k] = tmp;
-                              for (octave_idx_type i = cidx (k);
-                                   i < cidx (k+1)-1; i++)
-                                {
-                                  octave_idx_type iidx = ridx (i);
-                                  work[iidx] = work[iidx] - tmp * data (i);
-                                }
+                              #pragma omp parallel for
+                              for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
+                                work[ridx (i)] -= tmp * data (i);
                             }
                         }
                       double atmp = 0;
+                      #pragma omp parallel for reduction(+:atmp)
                       for (octave_idx_type i = 0; i < j+1; i++)
                         {
                           atmp += std::abs (work[i]);
@@ -2641,6 +2629,7 @@
               for (octave_idx_type j = 0; j < nc; j++)
                 {
                   double atmp = 0.;
+                  #pragma omp parallel for reduction(+:atmp)
                   for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
                     atmp += std::abs (data (i));
                   if (atmp > anorm)
@@ -2738,6 +2727,7 @@
                         }
 
                       double atmp = 0;
+                      #pragma omp parallel for reduction(+:atmp)
                       for (octave_idx_type i = j; i < nc; i++)
                         {
                           atmp += std::abs (work[i]);
@@ -2773,11 +2763,9 @@
 
                           Complex tmp = work[k] / data (cidx (k));
                           work[k] = tmp;
+                          #pragma omp parallel for
                           for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
-                            {
-                              octave_idx_type iidx = ridx (i);
-                              work[iidx] = work[iidx] - tmp * data (i);
-                            }
+                            work[ridx (i)] -= tmp * data (i);
                         }
                     }
                   for (octave_idx_type i = 0; i < nc; i++)
@@ -2801,15 +2789,13 @@
                             {
                               Complex tmp = work[k] / data (cidx (k));
                               work[k] = tmp;
-                              for (octave_idx_type i = cidx (k)+1;
-                                   i < cidx (k+1); i++)
-                                {
-                                  octave_idx_type iidx = ridx (i);
-                                  work[iidx] = work[iidx] - tmp * data (i);
-                                }
+                              #pragma omp parallel for
+                              for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
+                                work[ridx (i)] -= tmp * data (i);
                             }
                         }
                       double atmp = 0;
+                      #pragma omp parallel for reduction(+:atmp)
                       for (octave_idx_type i = j; i < nc; i++)
                         {
                           atmp += std::abs (work[i]);
@@ -2897,6 +2883,7 @@
               for (octave_idx_type j = 0; j < nc; j++)
                 {
                   double atmp = 0.;
+                  #pragma omp parallel for reduction(+:atmp)
                   for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
                     atmp += std::abs (data (i));
                   if (atmp > anorm)
@@ -3022,6 +3009,7 @@
                         }
 
                       double atmp = 0;
+                      #pragma omp parallel for reduction(+:atmp)
                       for (octave_idx_type i = j; i < nc; i++)
                         {
                           atmp += std::abs (work[i]);
@@ -3057,11 +3045,9 @@
 
                           Complex tmp = work[k] / data (cidx (k));
                           work[k] = tmp;
+                          #pragma omp parallel for
                           for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
-                            {
-                              octave_idx_type iidx = ridx (i);
-                              work[iidx] = work[iidx] - tmp * data (i);
-                            }
+                            work[ridx (i)] -= tmp * data (i);
                         }
                     }
 
@@ -3108,15 +3094,13 @@
                             {
                               Complex tmp = work[k] / data (cidx (k));
                               work[k] = tmp;
-                              for (octave_idx_type i = cidx (k)+1;
-                                   i < cidx (k+1); i++)
-                                {
-                                  octave_idx_type iidx = ridx (i);
-                                  work[iidx] = work[iidx] - tmp * data (i);
-                                }
+                              #pragma omp parallel for
+                              for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
+                                work[ridx (i)] -= tmp * data (i);
                             }
                         }
                       double atmp = 0;
+                      #pragma omp parallel for reduction(+:atmp)
                       for (octave_idx_type i = j; i < nc; i++)
                         {
                           atmp += std::abs (work[i]);
@@ -3205,6 +3189,7 @@
               for (octave_idx_type j = 0; j < nc; j++)
                 {
                   double atmp = 0.;
+                  #pragma omp parallel for reduction(+:atmp)
                   for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
                     atmp += std::abs (data (i));
                   if (atmp > anorm)
@@ -3302,6 +3287,7 @@
                         }
 
                       double atmp = 0;
+                      #pragma omp parallel for reduction(+:atmp)
                       for (octave_idx_type i = j; i < nc; i++)
                         {
                           atmp += std::abs (work[i]);
@@ -3339,11 +3325,9 @@
 
                           Complex tmp = work[k] / data (cidx (k));
                           work[k] = tmp;
+                          #pragma omp parallel for
                           for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
-                            {
-                              octave_idx_type iidx = ridx (i);
-                              work[iidx] = work[iidx] - tmp * data (i);
-                            }
+                            work[ridx (i)] -= tmp * data (i);
                         }
                     }
 
@@ -3368,15 +3352,13 @@
                             {
                               Complex tmp = work[k] / data (cidx (k));
                               work[k] = tmp;
-                              for (octave_idx_type i = cidx (k)+1;
-                                   i < cidx (k+1); i++)
-                                {
-                                  octave_idx_type iidx = ridx (i);
-                                  work[iidx] = work[iidx] - tmp * data (i);
-                                }
+                              #pragma omp parallel for
+                              for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
+                                work[ridx (i)] -= tmp * data (i);
                             }
                         }
                       double atmp = 0;
+                      #pragma omp parallel for reduction(+:atmp)
                       for (octave_idx_type i = j; i < nc; i++)
                         {
                           atmp += std::abs (work[i]);
@@ -3464,6 +3446,7 @@
               for (octave_idx_type j = 0; j < nc; j++)
                 {
                   double atmp = 0.;
+                  #pragma omp parallel for reduction(+:atmp)
                   for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
                     atmp += std::abs (data (i));
                   if (atmp > anorm)
@@ -3589,6 +3572,7 @@
                         }
 
                       double atmp = 0;
+                      #pragma omp parallel for reduction(+:atmp)
                       for (octave_idx_type i = j; i < nc; i++)
                         {
                           atmp += std::abs (work[i]);
@@ -3624,11 +3608,9 @@
 
                           Complex tmp = work[k] / data (cidx (k));
                           work[k] = tmp;
+                          #pragma omp parallel for
                           for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
-                            {
-                              octave_idx_type iidx = ridx (i);
-                              work[iidx] = work[iidx] - tmp * data (i);
-                            }
+                            work[ridx (i)] -= tmp * data (i);
                         }
                     }
 
@@ -3675,15 +3657,13 @@
                             {
                               Complex tmp = work[k] / data (cidx (k));
                               work[k] = tmp;
-                              for (octave_idx_type i = cidx (k)+1;
-                                   i < cidx (k+1); i++)
-                                {
-                                  octave_idx_type iidx = ridx (i);
-                                  work[iidx] = work[iidx] - tmp * data (i);
-                                }
+                              #pragma omp parallel for
+                              for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
+                                work[ridx (i)] -= tmp * data (i);
                             }
                         }
                       double atmp = 0;
+                      #pragma omp parallel for reduction(+:atmp)
                       for (octave_idx_type i = j; i < nc; i++)
                         {
                           atmp += std::abs (work[i]);
@@ -4500,6 +4480,7 @@
               for (octave_idx_type j = 0; j < nr; j++)
                 {
                   double atmp = 0.;
+                  #pragma omp parallel for reduction(+:atmp)
                   for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
                     atmp += std::abs (data (i));
                   if (atmp > anorm)
@@ -4782,6 +4763,7 @@
               for (octave_idx_type j = 0; j < nr; j++)
                 {
                   double atmp = 0.;
+                  #pragma omp parallel for reduction(+:atmp)
                   for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
                     atmp += std::abs (data (i));
                   if (atmp > anorm)
@@ -5067,6 +5049,7 @@
               for (octave_idx_type j = 0; j < nr; j++)
                 {
                   double atmp = 0.;
+                  #pragma omp parallel for reduction(+:atmp)
                   for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
                     atmp += std::abs (data (i));
                   if (atmp > anorm)
@@ -5354,6 +5337,7 @@
               for (octave_idx_type j = 0; j < nr; j++)
                 {
                   double atmp = 0.;
+                  #pragma omp parallel for reduction(+:atmp)
                   for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
                     atmp += std::abs (data (i));
                   if (atmp > anorm)
@@ -7375,6 +7359,7 @@
   for (octave_idx_type i = 0; i < nc + 1; i++)
     retval.cidx (i) = cidx (i);
 
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < nz; i++)
     {
       retval.data (i) = std::abs (data (i));
diff -r c31fd42f9000 -r 9a78540d8b38 liboctave/array/MSparse.cc
--- a/liboctave/array/MSparse.cc	Sat Mar 09 18:42:43 2013 +0000
+++ b/liboctave/array/MSparse.cc	Sun Mar 10 18:38:55 2013 -0300
@@ -169,6 +169,7 @@
 
   MSparse<T> r (nr, nc, nz);
 
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < nz; i++)
     {
       r.data (i) = op (a.data (i), s);
@@ -517,6 +518,7 @@
         {
           octave_idx_type b_nnz = b.nnz ();
           r = MSparse<T> (b);
+          #pragma omp parallel for
           for (octave_idx_type i = 0 ; i < b_nnz ; i++)
             r.data (i) = val / r.data (i);
           r.maybe_compress ();
diff -r c31fd42f9000 -r 9a78540d8b38 liboctave/array/Range.cc
--- a/liboctave/array/Range.cc	Sat Mar 09 18:42:43 2013 +0000
+++ b/liboctave/array/Range.cc	Sun Mar 10 18:38:55 2013 -0300
@@ -60,7 +60,8 @@
         {
           // The first element must always be *exactly* the base.
           // E.g, -0 would otherwise become +0 in the loop (-0 + 0*increment).
-          cache(0) = b; 
+          cache(0) = b;
+          #pragma omp parallel for
           for (octave_idx_type i = 1; i < rng_nelem; i++)
             cache(i) = b + i * increment;
         }
@@ -88,7 +89,7 @@
 
   if (i == 0)
     return rng_base;
-  else if (i < rng_nelem - 1) 
+  else if (i < rng_nelem - 1)
     return rng_base + i * rng_inc;
   else
     {
@@ -109,7 +110,7 @@
 #else
   if (i == 0)
     return rng_base;
-  else if (i < rng_nelem - 1) 
+  else if (i < rng_nelem - 1)
     return rng_base + i * rng_inc;
   else
     {
@@ -134,7 +135,7 @@
     {
       if (i == 0)
         *array++ = base;
-      else if (i < nmax) 
+      else if (i < nmax)
         *array++ = base + i * inc;
       else
         {
diff -r c31fd42f9000 -r 9a78540d8b38 liboctave/array/Sparse.cc
--- a/liboctave/array/Sparse.cc	Sat Mar 09 18:42:43 2013 +0000
+++ b/liboctave/array/Sparse.cc	Sun Mar 10 18:38:55 2013 -0300
@@ -111,6 +111,7 @@
             }
         }
 
+      #pragma omp parallel for
       for (octave_idx_type j = _c + 1; j < ncols + 1; j++)
         c[j] = c[j] + 1;
 
@@ -327,6 +328,7 @@
           xridx(0) = r(0);
           xdata(0) = a(0);
 
+          #pragma omp parallel for
           for (octave_idx_type j = 0; j < nc; j++)
             xcidx (j+1) = j >= c(0);
         }
@@ -2202,6 +2204,7 @@
         }
       for (octave_idx_type k = 0; k < i; k++)
         mridx[k] = k;
+      #pragma omp parallel for
       for (octave_idx_type k = i; k < ns; k++)
         mridx[k] = k - ns + nr;
 
diff -r c31fd42f9000 -r 9a78540d8b38 liboctave/array/boolSparse.cc
--- a/liboctave/array/boolSparse.cc	Sat Mar 09 18:42:43 2013 +0000
+++ b/liboctave/array/boolSparse.cc	Sun Mar 10 18:38:55 2013 -0300
@@ -213,6 +213,7 @@
         {
           // We can use O(nr) memory.
           Array<double> tmp (dim_vector (nr, 1), 0);
+          #pragma omp parallel for
           for (octave_idx_type i = 0; i < nz; i++)
             tmp.xelem (ridx (i)) += 1.0;
           retval = tmp;
diff -r c31fd42f9000 -r 9a78540d8b38 liboctave/array/dColVector.cc
--- a/liboctave/array/dColVector.cc	Sat Mar 09 18:42:43 2013 +0000
+++ b/liboctave/array/dColVector.cc	Sun Mar 10 18:38:55 2013 -0300
@@ -251,6 +251,7 @@
         {
           retval.resize (nr);
 
+          #pragma omp parallel for
           for (octave_idx_type i = 0; i < a_len; i++)
             retval.elem (i) = a.elem (i) * m.elem (i, i);
 
@@ -299,9 +300,8 @@
 std::ostream&
 operator << (std::ostream& os, const ColumnVector& a)
 {
-//  int field_width = os.precision () + 7;
   for (octave_idx_type i = 0; i < a.length (); i++)
-    os << /* setw (field_width) << */ a.elem (i) << "\n";
+    os << a.elem (i) << "\n";
   return os;
 }
 
diff -r c31fd42f9000 -r 9a78540d8b38 liboctave/array/dDiagMatrix.cc
--- a/liboctave/array/dDiagMatrix.cc	Sat Mar 09 18:42:43 2013 +0000
+++ b/liboctave/array/dDiagMatrix.cc	Sun Mar 10 18:38:55 2013 -0300
@@ -299,6 +299,7 @@
 
   DiagMatrix retval (c, r);
 
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < len; i++)
     {
       if (elem (i, i) != 0.0)
@@ -330,8 +331,10 @@
 
   octave_idx_type len = c.length (), lenm = len < a_nc ? len : a_nc;
 
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < lenm; i++)
     c.dgxelem (i) = a.dgelem (i) * b.dgelem (i);
+
   for (octave_idx_type i = lenm; i < len; i++)
     c.dgxelem (i) = 0.0;
 
@@ -370,16 +373,14 @@
 std::ostream&
 operator << (std::ostream& os, const DiagMatrix& a)
 {
-//  int field_width = os.precision () + 7;
-
   for (octave_idx_type i = 0; i < a.rows (); i++)
     {
       for (octave_idx_type j = 0; j < a.cols (); j++)
         {
           if (i == j)
-            os << " " /* setw (field_width) */ << a.elem (i, i);
+            os << " " << a.elem (i, i);
           else
-            os << " " /* setw (field_width) */ << 0.0;
+            os << " " << 0.0;
         }
       os << "\n";
     }
diff -r c31fd42f9000 -r 9a78540d8b38 liboctave/array/dMatrix.cc
--- a/liboctave/array/dMatrix.cc	Sat Mar 09 18:42:43 2013 +0000
+++ b/liboctave/array/dMatrix.cc	Sun Mar 10 18:38:55 2013 -0300
@@ -1068,6 +1068,7 @@
       F77_FUNC (zfftb, ZFFTB) (npts, &tmp_data[npts*j], pwsave);
     }
 
+  #pragma omp parallel for
   for (octave_idx_type j = 0; j < npts*nsamples; j++)
     tmp_data[j] = tmp_data[j] / static_cast<double> (npts);
 
@@ -1178,6 +1179,7 @@
       F77_FUNC (zfftb, ZFFTB) (npts, &tmp_data[npts*j], pwsave);
     }
 
+  #pragma omp parallel for
   for (octave_idx_type j = 0; j < npts*nsamples; j++)
     tmp_data[j] = tmp_data[j] / static_cast<double> (npts);
 
@@ -1202,6 +1204,7 @@
 
       F77_FUNC (zfftb, ZFFTB) (npts, prow, pwsave);
 
+      #pragma omp parallel for
       for (octave_idx_type i = 0; i < npts; i++)
         tmp_data[i*nr + j] = prow[i] / static_cast<double> (npts);
     }
@@ -2012,6 +2015,7 @@
   Matrix retval (m, 2*n);
   const Complex *cmd = cm.data ();
   double *rd = retval.fortran_vec ();
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < nel; i++)
     {
       rd[i] = std::real (cmd[i]);
@@ -2027,6 +2031,7 @@
   ComplexMatrix retval (m, n);
   const double *smd = sm.data ();
   Complex *rd = retval.fortran_vec ();
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < nel; i++)
     rd[i] = Complex (smd[i], smd[nel+i]);
   return retval;
@@ -2590,6 +2595,7 @@
       return *this;
     }
 
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < a.length (); i++)
     elem (i, i) += a.elem (i, i);
 
@@ -2611,6 +2617,7 @@
       return *this;
     }
 
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < a.length (); i++)
     elem (i, i) -= a.elem (i, i);
 
@@ -3422,15 +3429,21 @@
 
   // The last column is not needed while using delta.
   double *delta = &retval(0, n-1);
-  for (octave_idx_type i = 0; i < m; i++)
-    delta[i] = (x2(i) - x1(i)) / (n - 1);
-
-  for (octave_idx_type j = 1; j < n-1; j++)
+  #pragma omp parallel
+  {
+    #pragma omp for
     for (octave_idx_type i = 0; i < m; i++)
-      retval(i, j) = x1(i) + j*delta[i];
-
-  for (octave_idx_type i = 0; i < m; i++)
-    retval(i, n-1) = x2(i);
+      delta[i] = (x2(i) - x1(i)) / (n - 1);
+
+    #pragma omp for
+    for (octave_idx_type j = 1; j < n-1; j++)
+      for (octave_idx_type i = 0; i < m; i++)
+        retval(i, j) = x1(i) + j*delta[i];
+
+    #pragma omp for
+    for (octave_idx_type i = 0; i < m; i++)
+      retval(i, n-1) = x2(i);
+  }
 
   return retval;
 }
diff -r c31fd42f9000 -r 9a78540d8b38 liboctave/array/dNDArray.cc
--- a/liboctave/array/dNDArray.cc	Sat Mar 09 18:42:43 2013 +0000
+++ b/liboctave/array/dNDArray.cc	Sun Mar 10 18:38:55 2013 -0300
@@ -77,10 +77,12 @@
   else
     {
       if (zero_based)
+        #pragma omp parallel for
         for (octave_idx_type i = 0; i < a.numel (); i++)
           ptmp[i] = static_cast<double>
             (pa[i] + static_cast<octave_idx_type> (1));
       else
+        #pragma omp parallel for
         for (octave_idx_type i = 0; i < a.numel (); i++)
           ptmp[i] = static_cast<double> (pa[i]);
     }
@@ -90,6 +92,7 @@
   : MArray<double> (a.dims ())
 {
   octave_idx_type n = a.numel ();
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < n; i++)
     xelem (i) = static_cast<unsigned char> (a(i));
 }
@@ -107,6 +110,7 @@
   octave_idx_type stride = 1;
   octave_idx_type n = dv(dim);
 
+  #pragma omp parallel for reduction(*:stride)
   for (int i = 0; i < dim; i++)
     stride *= dv(i);
 
@@ -138,6 +142,7 @@
   octave_idx_type stride = 1;
   octave_idx_type n = dv(dim);
 
+  #pragma omp parallel for reduction(*:stride)
   for (int i = 0; i < dim; i++)
     stride *= dv(i);
 
@@ -264,6 +269,7 @@
 
   octave_idx_type stride = 1;
 
+  #pragma omp parallel for reduction(*:stride)
   for (int i = 0; i < dim; i++)
     stride *= dv(i);
 
@@ -311,6 +317,7 @@
 
   octave_idx_type stride = 1;
 
+  #pragma omp parallel for reduction(*:stride)
   for (int i = 0; i < dim; i++)
     stride *= dv(i);
 
@@ -332,6 +339,7 @@
 
           F77_FUNC (zfftb, ZFFTB) (npts, tmp, pwsave);
 
+          #pragma omp parallel for
           for (octave_idx_type i = 0; i < npts; i++)
             retval((i + k*npts)*stride + j*dist) = tmp[i] /
               static_cast<double> (npts);
@@ -426,6 +434,7 @@
 
               F77_FUNC (zfftb, ZFFTB) (npts, prow, pwsave);
 
+              #pragma omp parallel for
               for (octave_idx_type l = 0; l < npts; l++)
                 retval((l + k*npts)*stride + j*dist) = prow[l] /
                   static_cast<double> (npts);
@@ -521,6 +530,7 @@
 
               F77_FUNC (zfftb, ZFFTB) (npts, prow, pwsave);
 
+              #pragma omp parallel for
               for (octave_idx_type l = 0; l < npts; l++)
                 retval((l + k*npts)*stride + j*dist) = prow[l] /
                   static_cast<double> (npts);
diff -r c31fd42f9000 -r 9a78540d8b38 liboctave/array/dRowVector.cc
--- a/liboctave/array/dRowVector.cc	Sat Mar 09 18:42:43 2013 +0000
+++ b/liboctave/array/dRowVector.cc	Sun Mar 10 18:38:55 2013 -0300
@@ -261,10 +261,8 @@
 std::ostream&
 operator << (std::ostream& os, const RowVector& a)
 {
-//  int field_width = os.precision () + 7;
-
   for (octave_idx_type i = 0; i < a.length (); i++)
-    os << " " /* setw (field_width) */ << a.elem (i);
+    os << " " << a.elem (i);
   return os;
 }
 
@@ -299,8 +297,10 @@
 
   double delta = (x2 - x1) / (n - 1);
   retval(0) = x1;
+  #pragma omp parallel for
   for (octave_idx_type i = 1; i < n-1; i++)
     retval(i) = x1 + i*delta;
+
   retval(n-1) = x2;
 
   return retval;
diff -r c31fd42f9000 -r 9a78540d8b38 liboctave/array/dSparse.cc
--- a/liboctave/array/dSparse.cc	Sat Mar 09 18:42:43 2013 +0000
+++ b/liboctave/array/dSparse.cc	Sun Mar 10 18:38:55 2013 -0300
@@ -622,6 +622,7 @@
   for (octave_idx_type i = 0; i < nc +1; i++)
     r.cidx (i) = a.cidx (i);
 
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < nz; i++)
     {
       r.data (i) = std::real (a.data (i));
@@ -642,6 +643,7 @@
   for (octave_idx_type i = 0; i < nc +1; i++)
     r.cidx (i) = a.cidx (i);
 
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < nz; i++)
     {
       r.data (i) = std::imag (a.data (i));
@@ -858,6 +860,7 @@
               rcond = dmin / dmax;
             }
 
+          #pragma omp parallel for
           for (octave_idx_type i = 0; i < nr; i++)
             v[i] = 1.0 / v[i];
         }
@@ -1006,6 +1009,7 @@
                     }
 
                   if (pivot != 1.0)
+                    #pragma omp parallel for
                     for (octave_idx_type j = cx_colstart; j < cx; j++)
                       retval.xdata (j) /= pivot;
                 }
@@ -1125,8 +1129,8 @@
               for (octave_idx_type j = 0; j < nr; j++)
                 {
                   double atmp = 0.;
-                  for (octave_idx_type i = retval.cidx (j);
-                       i < retval.cidx (j+1); i++)
+                  #pragma omp parallel for reduction(+:atmp)
+                  for (octave_idx_type i = retval.cidx (j); i < retval.cidx (j+1); i++)
                     atmp += fabs (retval.data (i));
                   if (atmp > ainvnorm)
                     ainvnorm = atmp;
@@ -1368,10 +1372,12 @@
         {
           retval.resize (nc, b.cols (), 0.);
           if (typ == MatrixType::Diagonal)
+            #pragma omp parallel for collapse(2)
             for (octave_idx_type j = 0; j < b.cols (); j++)
               for (octave_idx_type i = 0; i < nm; i++)
                 retval(i,j) = b(i,j) / data (i);
           else
+            #pragma omp parallel for collapse(2)
             for (octave_idx_type j = 0; j < b.cols (); j++)
               for (octave_idx_type k = 0; k < nc; k++)
                 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
@@ -1518,10 +1524,12 @@
         {
           retval.resize (nc, b.cols (), 0);
           if (typ == MatrixType::Diagonal)
+            #pragma omp parallel for collapse(2)
             for (octave_idx_type j = 0; j < b.cols (); j++)
                 for (octave_idx_type i = 0; i < nm; i++)
                   retval(i,j) = b(i,j) / data (i);
           else
+            #pragma omp parallel for collapse(2)
             for (octave_idx_type j = 0; j < b.cols (); j++)
               for (octave_idx_type k = 0; k < nc; k++)
                 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
@@ -1678,6 +1686,7 @@
               for (octave_idx_type j = 0; j < nc; j++)
                 {
                   double atmp = 0.;
+                  #pragma omp parallel for reduction(+:atmp)
                   for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
                     atmp += fabs (data (i));
                   if (atmp > anorm)
@@ -1713,12 +1722,9 @@
 
                           double tmp = work[k] / data (cidx (kidx+1)-1);
                           work[k] = tmp;
-                          for (octave_idx_type i = cidx (kidx);
-                               i < cidx (kidx+1)-1; i++)
-                            {
-                              octave_idx_type iidx = ridx (i);
-                              work[iidx] = work[iidx] - tmp * data (i);
-                            }
+                          #pragma omp parallel for
+                          for (octave_idx_type i = cidx (kidx); i < cidx (kidx+1)-1; i++)
+                            work[ridx (i)] -= tmp * data (i);
                         }
                     }
 
@@ -1744,15 +1750,13 @@
                             {
                               double tmp = work[k] / data (cidx (iidx+1)-1);
                               work[k] = tmp;
-                              for (octave_idx_type i = cidx (iidx);
-                                   i < cidx (iidx+1)-1; i++)
-                                {
-                                  octave_idx_type idx2 = ridx (i);
-                                  work[idx2] = work[idx2] - tmp * data (i);
-                                }
+                              #pragma omp parallel for
+                              for (octave_idx_type i = cidx (iidx); i < cidx (iidx+1)-1; i++)
+                                work[ridx (i)] -= tmp * data (i);
                             }
                         }
                       double atmp = 0;
+                      #pragma omp parallel for reduction(+:atmp)
                       for (octave_idx_type i = 0; i < j+1; i++)
                         {
                           atmp += fabs (work[i]);
@@ -1789,11 +1793,9 @@
 
                           double tmp = work[k] / data (cidx (k+1)-1);
                           work[k] = tmp;
+                          #pragma omp parallel for
                           for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
-                            {
-                              octave_idx_type iidx = ridx (i);
-                              work[iidx] = work[iidx] - tmp * data (i);
-                            }
+                            work[ridx (i)] -= tmp * data (i);
                         }
                     }
 
@@ -1817,14 +1819,13 @@
                             {
                               double tmp = work[k] / data (cidx (k+1)-1);
                               work[k] = tmp;
+                              #pragma omp parallel for
                               for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
-                                {
-                                  octave_idx_type iidx = ridx (i);
-                                  work[iidx] = work[iidx] - tmp * data (i);
-                                }
+                                work[ridx (i)] -= tmp * data (i);
                             }
                         }
                       double atmp = 0;
+                      #pragma omp parallel for reduction(+:atmp)
                       for (octave_idx_type i = 0; i < j+1; i++)
                         {
                           atmp += fabs (work[i]);
@@ -1912,6 +1913,7 @@
               for (octave_idx_type j = 0; j < nc; j++)
                 {
                   double atmp = 0.;
+                  #pragma omp parallel for reduction(+:atmp)
                   for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
                     atmp += fabs (data (i));
                   if (atmp > anorm)
@@ -1957,12 +1959,9 @@
 
                           double tmp = work[k] / data (cidx (kidx+1)-1);
                           work[k] = tmp;
-                          for (octave_idx_type i = cidx (kidx);
-                               i < cidx (kidx+1)-1; i++)
-                            {
-                              octave_idx_type iidx = ridx (i);
-                              work[iidx] = work[iidx] - tmp * data (i);
-                            }
+                          #pragma omp parallel for
+                          for (octave_idx_type i = cidx (kidx); i < cidx (kidx+1)-1; i++)
+                            work[ridx (i)] -= tmp * data (i);
                         }
                     }
 
@@ -2010,15 +2009,13 @@
                             {
                               double tmp = work[k] / data (cidx (iidx+1)-1);
                               work[k] = tmp;
-                              for (octave_idx_type i = cidx (iidx);
-                                   i < cidx (iidx+1)-1; i++)
-                                {
-                                  octave_idx_type idx2 = ridx (i);
-                                  work[idx2] = work[idx2] - tmp * data (i);
-                                }
+                              #pragma omp parallel for
+                              for (octave_idx_type i = cidx (iidx); i < cidx (iidx+1)-1; i++)
+                                work[ridx (i)] -= tmp * data (i);
                             }
                         }
                       double atmp = 0;
+                      #pragma omp parallel for reduction(+:atmp)
                       for (octave_idx_type i = 0; i < j+1; i++)
                         {
                           atmp += fabs (work[i]);
@@ -2054,11 +2051,9 @@
 
                           double tmp = work[k] / data (cidx (k+1)-1);
                           work[k] = tmp;
+                          #pragma omp parallel for
                           for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
-                            {
-                              octave_idx_type iidx = ridx (i);
-                              work[iidx] = work[iidx] - tmp * data (i);
-                            }
+                            work[ridx (i)] -= tmp * data (i);
                         }
                     }
 
@@ -2104,15 +2099,13 @@
                             {
                               double tmp = work[k] / data (cidx (k+1)-1);
                               work[k] = tmp;
-                              for (octave_idx_type i = cidx (k);
-                                   i < cidx (k+1)-1; i++)
-                                {
-                                  octave_idx_type iidx = ridx (i);
-                                  work[iidx] = work[iidx] - tmp * data (i);
-                                }
+                              #pragma omp parallel for
+                              for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
+                                work[ridx (i)] -= tmp * data (i);
                             }
                         }
                       double atmp = 0;
+                      #pragma omp parallel for reduction(+:atmp)
                       for (octave_idx_type i = 0; i < j+1; i++)
                         {
                           atmp += fabs (work[i]);
@@ -2200,6 +2193,7 @@
               for (octave_idx_type j = 0; j < nc; j++)
                 {
                   double atmp = 0.;
+                  #pragma omp parallel for reduction(+:atmp)
                   for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
                     atmp += fabs (data (i));
                   if (atmp > anorm)
@@ -2235,12 +2229,9 @@
 
                           Complex tmp = cwork[k] / data (cidx (kidx+1)-1);
                           cwork[k] = tmp;
-                          for (octave_idx_type i = cidx (kidx);
-                               i < cidx (kidx+1)-1; i++)
-                            {
-                              octave_idx_type iidx = ridx (i);
-                              cwork[iidx] = cwork[iidx] - tmp * data (i);
-                            }
+                          #pragma omp parallel for
+                          for (octave_idx_type i = cidx (kidx); i < cidx (kidx+1)-1; i++)
+                            cwork[ridx (i)] -= tmp * data (i);
                         }
                     }
 
@@ -2267,15 +2258,13 @@
                             {
                               double tmp = work[k] / data (cidx (iidx+1)-1);
                               work[k] = tmp;
-                              for (octave_idx_type i = cidx (iidx);
-                                   i < cidx (iidx+1)-1; i++)
-                                {
-                                  octave_idx_type idx2 = ridx (i);
-                                  work[idx2] = work[idx2] - tmp * data (i);
-                                }
+                              #pragma omp parallel for
+                              for (octave_idx_type i = cidx (iidx); i < cidx (iidx+1)-1; i++)
+                                work[ridx (i)] -= tmp * data (i);
                             }
                         }
                       double atmp = 0;
+                      #pragma omp parallel for reduction(+:atmp)
                       for (octave_idx_type i = 0; i < j+1; i++)
                         {
                           atmp += fabs (work[i]);
@@ -2312,11 +2301,9 @@
 
                           Complex tmp = cwork[k] / data (cidx (k+1)-1);
                           cwork[k] = tmp;
+                          #pragma omp parallel for
                           for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
-                            {
-                              octave_idx_type iidx = ridx (i);
-                              cwork[iidx] = cwork[iidx] - tmp  * data (i);
-                            }
+                            cwork[ridx (i)] -= tmp  * data (i);
                         }
                     }
 
@@ -2341,15 +2328,13 @@
                             {
                               double tmp = work[k] / data (cidx (k+1)-1);
                               work[k] = tmp;
-                              for (octave_idx_type i = cidx (k);
-                                   i < cidx (k+1)-1; i++)
-                                {
-                                  octave_idx_type iidx = ridx (i);
-                                  work[iidx] = work[iidx] - tmp * data (i);
-                                }
+                              #pragma omp parallel for
+                              for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
+                                work[ridx (i)] -= tmp * data (i);
                             }
                         }
                       double atmp = 0;
+                      #pragma omp parallel for reduction(+:atmp)
                       for (octave_idx_type i = 0; i < j+1; i++)
                         {
                           atmp += fabs (work[i]);
@@ -2482,12 +2467,9 @@
 
                           Complex tmp = cwork[k] / data (cidx (kidx+1)-1);
                           cwork[k] = tmp;
-                          for (octave_idx_type i = cidx (kidx);
-                               i < cidx (kidx+1)-1; i++)
-                            {
-                              octave_idx_type iidx = ridx (i);
-                              cwork[iidx] = cwork[iidx] - tmp * data (i);
-                            }
+                          #pragma omp parallel for
+                          for (octave_idx_type i = cidx (kidx); i < cidx (kidx+1)-1; i++)
+                            cwork[ridx (i)] -= tmp * data (i);
                         }
                     }
 
@@ -2536,15 +2518,13 @@
                             {
                               double tmp = work[k] / data (cidx (iidx+1)-1);
                               work[k] = tmp;
-                              for (octave_idx_type i = cidx (iidx);
-                                   i < cidx (iidx+1)-1; i++)
-                                {
-                                  octave_idx_type idx2 = ridx (i);
-                                  work[idx2] = work[idx2] - tmp * data (i);
-                                }
+                              #pragma omp parallel for
+                              for (octave_idx_type i = cidx (iidx); i < cidx (iidx+1)-1; i++)
+                                work[ridx (i)] -= tmp * data (i);
                             }
                         }
                       double atmp = 0;
+                      #pragma omp parallel for reduction(+:atmp)
                       for (octave_idx_type i = 0; i < j+1; i++)
                         {
                           atmp += fabs (work[i]);
@@ -2580,11 +2560,9 @@
 
                           Complex tmp = cwork[k] / data (cidx (k+1)-1);
                           cwork[k] = tmp;
+                          #pragma omp parallel for
                           for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
-                            {
-                              octave_idx_type iidx = ridx (i);
-                              cwork[iidx] = cwork[iidx] - tmp * data (i);
-                            }
+                            cwork[ridx (i)] -= tmp * data (i);
                         }
                     }
 
@@ -2631,15 +2609,13 @@
                             {
                               double tmp = work[k] / data (cidx (k+1)-1);
                               work[k] = tmp;
-                              for (octave_idx_type i = cidx (k);
-                                   i < cidx (k+1)-1; i++)
-                                {
-                                  octave_idx_type iidx = ridx (i);
-                                  work[iidx] = work[iidx] - tmp * data (i);
-                                }
+                              #pragma omp parallel for
+                              for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
+                                work[ridx (i)] -= tmp * data (i);
                             }
                         }
                       double atmp = 0;
+                      #pragma omp parallel for reduction(+:atmp)
                       for (octave_idx_type i = 0; i < j+1; i++)
                         {
                           atmp += fabs (work[i]);
@@ -2728,6 +2704,7 @@
               for (octave_idx_type j = 0; j < nc; j++)
                 {
                   double atmp = 0.;
+                  #pragma omp parallel for reduction(+:atmp)
                   for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
                     atmp += fabs (data (i));
                   if (atmp > anorm)
@@ -2771,13 +2748,14 @@
 
                           double tmp = work[k] / data (mini);
                           work[k] = tmp;
+                          #pragma omp parallel for
                           for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
                             {
                               if (i == mini)
                                 continue;
 
                               octave_idx_type iidx = perm[ridx (i)];
-                              work[iidx] = work[iidx] - tmp * data (i);
+                              work[iidx] -= tmp * data (i);
                             }
                         }
                     }
@@ -2813,19 +2791,20 @@
 
                               double tmp = work[k] / data (mini);
                               work[k] = tmp;
-                              for (octave_idx_type i = cidx (k);
-                                   i < cidx (k+1); i++)
+                              #pragma omp parallel for
+                              for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
                                 {
                                   if (i == mini)
                                     continue;
 
                                   octave_idx_type iidx = perm[ridx (i)];
-                                  work[iidx] = work[iidx] - tmp * data (i);
+                                  work[iidx] -= tmp * data (i);
                                 }
                             }
                         }
 
                       double atmp = 0;
+                      #pragma omp parallel for reduction(+:atmp)
                       for (octave_idx_type i = j; i < nc; i++)
                         {
                           atmp += fabs (work[i]);
@@ -2861,12 +2840,9 @@
 
                           double tmp = work[k] / data (cidx (k));
                           work[k] = tmp;
-                          for (octave_idx_type i = cidx (k)+1;
-                               i < cidx (k+1); i++)
-                            {
-                              octave_idx_type iidx = ridx (i);
-                              work[iidx] = work[iidx] - tmp * data (i);
-                            }
+                          #pragma omp parallel for
+                          for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
+                            work[ridx (i)] -= tmp * data (i);
                         }
                     }
 
@@ -2891,15 +2867,13 @@
                             {
                               double tmp = work[k] / data (cidx (k));
                               work[k] = tmp;
-                              for (octave_idx_type i = cidx (k)+1;
-                                   i < cidx (k+1); i++)
-                                {
-                                  octave_idx_type iidx = ridx (i);
-                                  work[iidx] = work[iidx] - tmp * data (i);
-                                }
+                              #pragma omp parallel for
+                              for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
+                                work[ridx (i)] -= tmp * data (i);
                             }
                         }
                       double atmp = 0;
+                      #pragma omp parallel for reduction(+:atmp)
                       for (octave_idx_type i = j; i < nc; i++)
                         {
                           atmp += fabs (work[i]);
@@ -2987,6 +2961,7 @@
               for (octave_idx_type j = 0; j < nc; j++)
                 {
                   double atmp = 0.;
+                  #pragma omp parallel for reduction(+:atmp)
                   for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
                     atmp += fabs (data (i));
                   if (atmp > anorm)
@@ -3035,13 +3010,14 @@
 
                           double tmp = work[k] / data (mini);
                           work[k] = tmp;
+                          #pragma omp parallel for
                           for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
                             {
                               if (i == mini)
                                 continue;
 
                               octave_idx_type iidx = perm[ridx (i)];
-                              work[iidx] = work[iidx] - tmp * data (i);
+                              work[iidx] -= tmp * data (i);
                             }
                         }
                     }
@@ -3099,19 +3075,20 @@
 
                               double tmp = work[k] / data (mini);
                               work[k] = tmp;
-                              for (octave_idx_type i = cidx (k);
-                                   i < cidx (k+1); i++)
+                              #pragma omp parallel for
+                              for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
                                 {
                                   if (i == mini)
                                     continue;
 
                                   octave_idx_type iidx = perm[ridx (i)];
-                                  work[iidx] = work[iidx] - tmp * data (i);
+                                  work[iidx] -= tmp * data (i);
                                 }
                             }
                         }
 
                       double atmp = 0;
+                      #pragma omp parallel for reduction(+:atmp)
                       for (octave_idx_type i = j; i < nr; i++)
                         {
                           atmp += fabs (work[i]);
@@ -3147,11 +3124,9 @@
 
                           double tmp = work[k] / data (cidx (k));
                           work[k] = tmp;
+                          #pragma omp parallel for
                           for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
-                            {
-                              octave_idx_type iidx = ridx (i);
-                              work[iidx] = work[iidx] - tmp * data (i);
-                            }
+                            work[ridx (i)] -= tmp * data (i);
                         }
                     }
 
@@ -3198,15 +3173,13 @@
                             {
                               double tmp = work[k] / data (cidx (k));
                               work[k] = tmp;
-                              for (octave_idx_type i = cidx (k)+1;
-                                   i < cidx (k+1); i++)
-                                {
-                                  octave_idx_type iidx = ridx (i);
-                                  work[iidx] = work[iidx] - tmp * data (i);
-                                }
+                              #pragma omp parallel for
+                              for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
+                                work[ridx (i)] -= tmp * data (i);
                             }
                         }
                       double atmp = 0;
+                      #pragma omp parallel for reduction(+:atmp)
                       for (octave_idx_type i = j; i < nc; i++)
                         {
                           atmp += fabs (work[i]);
@@ -3295,6 +3268,7 @@
               for (octave_idx_type j = 0; j < nc; j++)
                 {
                   double atmp = 0.;
+                  #pragma omp parallel for reduction(+:atmp)
                   for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
                     atmp += fabs (data (i));
                   if (atmp > anorm)
@@ -3337,13 +3311,14 @@
 
                           Complex tmp = cwork[k] / data (mini);
                           cwork[k] = tmp;
+                          #pragma omp parallel for
                           for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
                             {
                               if (i == mini)
                                 continue;
 
                               octave_idx_type iidx = perm[ridx (i)];
-                              cwork[iidx] = cwork[iidx] - tmp * data (i);
+                              cwork[iidx] -= tmp * data (i);
                             }
                         }
                     }
@@ -3380,19 +3355,20 @@
 
                               double tmp = work[k] / data (mini);
                               work[k] = tmp;
-                              for (octave_idx_type i = cidx (k);
-                                   i < cidx (k+1); i++)
+                              #pragma omp parallel for
+                              for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
                                 {
                                   if (i == mini)
                                     continue;
 
                                   octave_idx_type iidx = perm[ridx (i)];
-                                  work[iidx] = work[iidx] - tmp * data (i);
+                                  work[iidx] -= tmp * data (i);
                                 }
                             }
                         }
 
                       double atmp = 0;
+                      #pragma omp parallel for reduction(+:atmp)
                       for (octave_idx_type i = j; i < nc; i++)
                         {
                           atmp += fabs (work[i]);
@@ -3429,11 +3405,9 @@
 
                           Complex tmp = cwork[k] / data (cidx (k));
                           cwork[k] = tmp;
+                          #pragma omp parallel for
                           for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
-                            {
-                              octave_idx_type iidx = ridx (i);
-                              cwork[iidx] = cwork[iidx] - tmp * data (i);
-                            }
+                            cwork[ridx (i)] -= tmp * data (i);
                         }
                     }
 
@@ -3459,15 +3433,13 @@
                             {
                               double tmp = work[k] / data (cidx (k));
                               work[k] = tmp;
-                              for (octave_idx_type i = cidx (k)+1;
-                                   i < cidx (k+1); i++)
-                                {
-                                  octave_idx_type iidx = ridx (i);
-                                  work[iidx] = work[iidx] - tmp * data (i);
-                                }
+                              #pragma omp parallel for
+                              for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
+                                work[ridx (i)] -= tmp * data (i);
                             }
                         }
                       double atmp = 0;
+                      #pragma omp parallel for reduction(+:atmp)
                       for (octave_idx_type i = j; i < nc; i++)
                         {
                           atmp += fabs (work[i]);
@@ -3555,6 +3527,7 @@
               for (octave_idx_type j = 0; j < nc; j++)
                 {
                   double atmp = 0.;
+                  #pragma omp parallel for reduction(+:atmp)
                   for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
                     atmp += fabs (data (i));
                   if (atmp > anorm)
@@ -3603,13 +3576,14 @@
 
                           Complex tmp = cwork[k] / data (mini);
                           cwork[k] = tmp;
+                          #pragma omp parallel for
                           for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
                             {
                               if (i == mini)
                                 continue;
 
                               octave_idx_type iidx = perm[ridx (i)];
-                              cwork[iidx] = cwork[iidx] - tmp * data (i);
+                              cwork[iidx] -= tmp * data (i);
                             }
                         }
                     }
@@ -3668,19 +3642,20 @@
 
                               double tmp = work[k] / data (mini);
                               work[k] = tmp;
-                              for (octave_idx_type i = cidx (k);
-                                   i < cidx (k+1); i++)
+                              #pragma omp parallel for
+                              for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
                                 {
                                   if (i == mini)
                                     continue;
 
                                   octave_idx_type iidx = perm[ridx (i)];
-                                  work[iidx] = work[iidx] - tmp * data (i);
+                                  work[iidx] -= tmp * data (i);
                                 }
                             }
                         }
 
                       double atmp = 0;
+                      #pragma omp parallel for reduction(+:atmp)
                       for (octave_idx_type i = j; i < nc; i++)
                         {
                           atmp += fabs (work[i]);
@@ -3716,11 +3691,9 @@
 
                           Complex tmp = cwork[k] / data (cidx (k));
                           cwork[k] = tmp;
+                          #pragma omp parallel for
                           for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
-                            {
-                              octave_idx_type iidx = ridx (i);
-                              cwork[iidx] = cwork[iidx] - tmp * data (i);
-                            }
+                            cwork[ridx (i)] -= tmp * data (i);
                         }
                     }
 
@@ -3768,15 +3741,13 @@
                             {
                               double tmp = work[k] / data (cidx (k));
                               work[k] = tmp;
-                              for (octave_idx_type i = cidx (k)+1;
-                                   i < cidx (k+1); i++)
-                                {
-                                  octave_idx_type iidx = ridx (i);
-                                  work[iidx] = work[iidx] - tmp * data (i);
-                                }
+                              #pragma omp parallel for
+                              for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
+                                work[ridx (i)] -= tmp * data (i);
                             }
                         }
                       double atmp = 0;
+                      #pragma omp parallel for reduction(+:atmp)
                       for (octave_idx_type i = j; i < nc; i++)
                         {
                           atmp += fabs (work[i]);
@@ -4383,6 +4354,7 @@
               for (volatile octave_idx_type j = 0; j < b_nc; j++)
                 {
 
+                  #pragma omp parallel for
                   for (octave_idx_type i = 0; i < b_nr; i++)
                     {
                       Complex c = b (i,j);
@@ -4613,6 +4585,7 @@
               for (octave_idx_type j = 0; j < nr; j++)
                 {
                   double atmp = 0.;
+                  #pragma omp parallel for reduction(+:atmp)
                   for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
                     atmp += fabs (data (i));
                   if (atmp > anorm)
@@ -4896,6 +4869,7 @@
               for (octave_idx_type j = 0; j < nr; j++)
                 {
                   double atmp = 0.;
+                  #pragma omp parallel for reduction(+:atmp)
                   for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
                     atmp += fabs (data (i));
                   if (atmp > anorm)
@@ -5140,6 +5114,7 @@
 
                   for (volatile octave_idx_type j = 0; j < b_nc; j++)
                     {
+                      #pragma omp parallel for
                       for (octave_idx_type i = 0; i < b_nr; i++)
                         {
                           Complex c = b (i,j);
@@ -5175,6 +5150,7 @@
                           break;
                         }
 
+                      #pragma omp parallel for
                       for (octave_idx_type i = 0; i < b_nr; i++)
                         retval(i, j) = Complex (Bx[i], Bz[i]);
                     }
@@ -5212,6 +5188,7 @@
               for (octave_idx_type j = 0; j < nr; j++)
                 {
                   double atmp = 0.;
+                  #pragma omp parallel for reduction(+:atmp)
                   for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
                     atmp += fabs (data (i));
                   if (atmp > anorm)
@@ -5290,6 +5267,7 @@
 
                   for (volatile octave_idx_type j = 0; j < b_nc; j++)
                     {
+                      #pragma omp parallel for
                       for (octave_idx_type i = 0; i < nr; i++)
                         {
                           Complex c = b (i, j);
@@ -5309,6 +5287,7 @@
                                  ldm, pipvt, Bz, b.rows (), err
                                  F77_CHAR_ARG_LEN (1)));
 
+                      #pragma omp parallel for
                       for (octave_idx_type i = 0; i < nr; i++)
                         retval(i, j) = Complex (Bx[i], Bz[i]);
                     }
@@ -5444,7 +5423,7 @@
                   retval.xcidx (0) = 0;
                   for (volatile octave_idx_type j = 0; j < b_nc; j++)
                     {
-
+                      #pragma omp parallel for
                       for (octave_idx_type i = 0; i < b_nr; i++)
                         {
                           Complex c = b (i,j);
@@ -5542,6 +5521,7 @@
               for (octave_idx_type j = 0; j < nr; j++)
                 {
                   double atmp = 0.;
+                  #pragma omp parallel for reduction(+:atmp)
                   for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
                     atmp += fabs (data (i));
                   if (atmp > anorm)
@@ -5628,8 +5608,8 @@
                           Bx[i] = 0.;
                           Bz[i] = 0.;
                         }
-                      for (octave_idx_type i = b.cidx (j);
-                           i < b.cidx (j+1); i++)
+                      #pragma omp parallel for
+                      for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
                         {
                           Complex c = b.data (i);
                           Bx[b.ridx (i)] = std::real (c);
@@ -6460,6 +6440,7 @@
 
               for (octave_idx_type j = 0; j < b_nc; j++)
                 {
+                  #pragma omp parallel for
                   for (octave_idx_type i = 0; i < b_nr; i++)
                     {
                       Complex c = b (i,j);
@@ -6486,6 +6467,7 @@
                       break;
                     }
 
+                  #pragma omp parallel for
                   for (octave_idx_type i = 0; i < b_nr; i++)
                     retval(i, j) = Complex (Xx[i], Xz[i]);
                 }
@@ -6715,6 +6697,7 @@
               retval.xcidx (0) = 0;
               for (octave_idx_type j = 0; j < b_nc; j++)
                 {
+                  #pragma omp parallel for
                   for (octave_idx_type i = 0; i < b_nr; i++)
                     {
                       Complex c = b (i,j);
@@ -7525,6 +7508,7 @@
 
   SparseMatrix retval (*this);
 
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < nz; i++)
     retval.data (i) = fabs (retval.data (i));
 
diff -r c31fd42f9000 -r 9a78540d8b38 liboctave/array/fCColVector.cc
--- a/liboctave/array/fCColVector.cc	Sat Mar 09 18:42:43 2013 +0000
+++ b/liboctave/array/fCColVector.cc	Sun Mar 10 18:38:55 2013 -0300
@@ -393,6 +393,7 @@
 
   FloatComplexColumnVector result (nr);
 
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < a_len; i++)
     result.elem (i) = a.elem (i) * m.elem (i, i);
 
@@ -421,6 +422,7 @@
 
   FloatComplexColumnVector result (nr);
 
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < a_len; i++)
     result.elem (i) = a.elem (i) * m.elem (i, i);
 
@@ -449,6 +451,7 @@
 
   FloatComplexColumnVector result (nr);
 
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < a_len; i++)
     result.elem (i) = a.elem (i) * m.elem (i, i);
 
@@ -505,9 +508,8 @@
 std::ostream&
 operator << (std::ostream& os, const FloatComplexColumnVector& a)
 {
-//  int field_width = os.precision () + 7;
   for (octave_idx_type i = 0; i < a.length (); i++)
-    os << /* setw (field_width) << */ a.elem (i) << "\n";
+    os << a.elem (i) << "\n";
   return os;
 }
 
diff -r c31fd42f9000 -r 9a78540d8b38 liboctave/array/fCDiagMatrix.cc
--- a/liboctave/array/fCDiagMatrix.cc	Sat Mar 09 18:42:43 2013 +0000
+++ b/liboctave/array/fCDiagMatrix.cc	Sun Mar 10 18:38:55 2013 -0300
@@ -389,6 +389,7 @@
 
   FloatComplexDiagMatrix retval (c, r);
 
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < len; i++)
     {
       if (elem (i, i) != 0.0f)
@@ -448,8 +449,10 @@
 
   octave_idx_type len = c.length (), lenm = len < a_nc ? len : a_nc;
 
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < lenm; i++)
     c.dgxelem (i) = a.dgelem (i) * b.dgelem (i);
+
   for (octave_idx_type i = lenm; i < len; i++)
     c.dgxelem (i) = 0.0f;
 
@@ -478,6 +481,7 @@
 
   octave_idx_type len = a_nr < b_nc ? a_nr : b_nc;
 
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < len; i++)
     {
       float a_element = a.elem (i, i);
@@ -511,6 +515,7 @@
 
   octave_idx_type len = a_nr < b_nc ? a_nr : b_nc;
 
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < len; i++)
     {
       FloatComplex a_element = a.elem (i, i);
@@ -557,15 +562,15 @@
 operator << (std::ostream& os, const FloatComplexDiagMatrix& a)
 {
   FloatComplex ZERO (0.0);
-//  int field_width = os.precision () + 7;
+
   for (octave_idx_type i = 0; i < a.rows (); i++)
     {
       for (octave_idx_type j = 0; j < a.cols (); j++)
         {
           if (i == j)
-            os << " " /* setw (field_width) */ << a.elem (i, i);
+            os << " " << a.elem (i, i);
           else
-            os << " " /* setw (field_width) */ << ZERO;
+            os << " " << ZERO;
         }
       os << "\n";
     }
diff -r c31fd42f9000 -r 9a78540d8b38 liboctave/array/fCMatrix.cc
--- a/liboctave/array/fCMatrix.cc	Sat Mar 09 18:42:43 2013 +0000
+++ b/liboctave/array/fCMatrix.cc	Sun Mar 10 18:38:55 2013 -0300
@@ -325,6 +325,8 @@
     (*current_liboctave_error_handler) ("complex: internal error");
 
   octave_idx_type nel = numel ();
+
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < nel; i++)
     xelem (i) = FloatComplex (re(i), im(i));
 }
@@ -1396,6 +1398,7 @@
       F77_FUNC (cfftb, CFFTB) (npts, &tmp_data[npts*j], pwsave);
     }
 
+  #pragma omp parallel for
   for (octave_idx_type j = 0; j < npts*nsamples; j++)
     tmp_data[j] = tmp_data[j] / static_cast<float> (npts);
 
@@ -1506,6 +1509,7 @@
       F77_FUNC (cfftb, CFFTB) (npts, &tmp_data[npts*j], pwsave);
     }
 
+  #pragma omp parallel for
   for (octave_idx_type j = 0; j < npts*nsamples; j++)
     tmp_data[j] = tmp_data[j] / static_cast<float> (npts);
 
@@ -1530,6 +1534,7 @@
 
       F77_FUNC (cfftb, CFFTB) (npts, prow, pwsave);
 
+      #pragma omp parallel for
       for (octave_idx_type i = 0; i < npts; i++)
         tmp_data[i*nr + j] = prow[i] / static_cast<float> (npts);
     }
@@ -2954,6 +2959,7 @@
       return *this;
     }
 
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < a.length (); i++)
     elem (i, i) += a.elem (i, i);
 
@@ -2975,6 +2981,7 @@
       return *this;
     }
 
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < a.length (); i++)
     elem (i, i) -= a.elem (i, i);
 
@@ -2996,6 +3003,7 @@
       return *this;
     }
 
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < a.length (); i++)
     elem (i, i) += a.elem (i, i);
 
@@ -3017,6 +3025,7 @@
       return *this;
     }
 
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < a.length (); i++)
     elem (i, i) -= a.elem (i, i);
 
@@ -3934,6 +3943,7 @@
 
       if (columns_are_real_only)
         {
+          #pragma omp parallel for
           for (octave_idx_type i = 0; i < nr; i++)
             result (i, j) = xmin (std::real (a (i, j)), std::real (b (i, j)));
         }
@@ -4061,9 +4071,11 @@
 
   // The last column is not needed while using delta.
   FloatComplex *delta = &retval(0, n-1);
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < m; i++)
     delta[i] = (x2(i) - x1(i)) / (n - 1.0f);
 
+  #pragma omp parallel for collapse(2)
   for (octave_idx_type j = 1; j < n-1; j++)
     for (octave_idx_type i = 0; i < m; i++)
       retval(i, j) = x1(i) + static_cast<float> (j)*delta[i];
diff -r c31fd42f9000 -r 9a78540d8b38 liboctave/array/fCNDArray.cc
--- a/liboctave/array/fCNDArray.cc	Sat Mar 09 18:42:43 2013 +0000
+++ b/liboctave/array/fCNDArray.cc	Sun Mar 10 18:38:55 2013 -0300
@@ -65,6 +65,7 @@
   octave_idx_type stride = 1;
   octave_idx_type n = dv(dim);
 
+  #pragma omp parallel for reduction(*:stride)
   for (int i = 0; i < dim; i++)
     stride *= dv(i);
 
@@ -96,6 +97,7 @@
   octave_idx_type stride = 1;
   octave_idx_type n = dv(dim);
 
+  #pragma omp parallel for reduction(*:stride)
   for (int i = 0; i < dim; i++)
     stride *= dv(i);
 
@@ -220,6 +222,7 @@
 
   octave_idx_type stride = 1;
 
+  #pragma omp parallel for reduction(*:stride)
   for (int i = 0; i < dim; i++)
     stride *= dv(i);
 
@@ -267,6 +270,7 @@
 
   octave_idx_type stride = 1;
 
+  #pragma omp parallel for reduction(*:stride)
   for (int i = 0; i < dim; i++)
     stride *= dv(i);
 
diff -r c31fd42f9000 -r 9a78540d8b38 liboctave/array/fCRowVector.cc
--- a/liboctave/array/fCRowVector.cc	Sat Mar 09 18:42:43 2013 +0000
+++ b/liboctave/array/fCRowVector.cc	Sun Mar 10 18:38:55 2013 -0300
@@ -402,9 +402,8 @@
 std::ostream&
 operator << (std::ostream& os, const FloatComplexRowVector& a)
 {
-//  int field_width = os.precision () + 7;
   for (octave_idx_type i = 0; i < a.length (); i++)
-    os << " " /* setw (field_width) */ << a.elem (i);
+    os << " " << a.elem (i);
   return os;
 }
 
diff -r c31fd42f9000 -r 9a78540d8b38 liboctave/array/fColVector.cc
--- a/liboctave/array/fColVector.cc	Sat Mar 09 18:42:43 2013 +0000
+++ b/liboctave/array/fColVector.cc	Sun Mar 10 18:38:55 2013 -0300
@@ -250,6 +250,7 @@
         {
           retval.resize (nr);
 
+          #pragma omp parallel for
           for (octave_idx_type i = 0; i < a_len; i++)
             retval.elem (i) = a.elem (i) * m.elem (i, i);
 
@@ -298,9 +299,8 @@
 std::ostream&
 operator << (std::ostream& os, const FloatColumnVector& a)
 {
-//  int field_width = os.precision () + 7;
   for (octave_idx_type i = 0; i < a.length (); i++)
-    os << /* setw (field_width) << */ a.elem (i) << "\n";
+    os << a.elem (i) << "\n";
   return os;
 }
 
diff -r c31fd42f9000 -r 9a78540d8b38 liboctave/array/fDiagMatrix.cc
--- a/liboctave/array/fDiagMatrix.cc	Sat Mar 09 18:42:43 2013 +0000
+++ b/liboctave/array/fDiagMatrix.cc	Sun Mar 10 18:38:55 2013 -0300
@@ -299,6 +299,7 @@
 
   FloatDiagMatrix retval (c, r);
 
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < len; i++)
     {
       if (elem (i, i) != 0.0f)
@@ -330,8 +331,10 @@
 
   octave_idx_type len = c.length (), lenm = len < a_nc ? len : a_nc;
 
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < lenm; i++)
     c.dgxelem (i) = a.dgelem (i) * b.dgelem (i);
+
   for (octave_idx_type i = lenm; i < len; i++)
     c.dgxelem (i) = 0.0f;
 
@@ -370,16 +373,14 @@
 std::ostream&
 operator << (std::ostream& os, const FloatDiagMatrix& a)
 {
-//  int field_width = os.precision () + 7;
-
   for (octave_idx_type i = 0; i < a.rows (); i++)
     {
       for (octave_idx_type j = 0; j < a.cols (); j++)
         {
           if (i == j)
-            os << " " /* setw (field_width) */ << a.elem (i, i);
+            os << " " << a.elem (i, i);
           else
-            os << " " /* setw (field_width) */ << 0.0;
+            os << " " << 0.0;
         }
       os << "\n";
     }
diff -r c31fd42f9000 -r 9a78540d8b38 liboctave/array/fMatrix.cc
--- a/liboctave/array/fMatrix.cc	Sat Mar 09 18:42:43 2013 +0000
+++ b/liboctave/array/fMatrix.cc	Sun Mar 10 18:38:55 2013 -0300
@@ -1068,6 +1068,7 @@
       F77_FUNC (cfftb, CFFTB) (npts, &tmp_data[npts*j], pwsave);
     }
 
+  #pragma omp parallel for
   for (octave_idx_type j = 0; j < npts*nsamples; j++)
     tmp_data[j] = tmp_data[j] / static_cast<float> (npts);
 
@@ -1178,6 +1179,7 @@
       F77_FUNC (cfftb, CFFTB) (npts, &tmp_data[npts*j], pwsave);
     }
 
+  #pragma omp parallel for
   for (octave_idx_type j = 0; j < npts*nsamples; j++)
     tmp_data[j] = tmp_data[j] / static_cast<float> (npts);
 
@@ -2012,6 +2014,7 @@
   FloatMatrix retval (m, 2*n);
   const FloatComplex *cmd = cm.data ();
   float *rd = retval.fortran_vec ();
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < nel; i++)
     {
       rd[i] = std::real (cmd[i]);
@@ -2027,6 +2030,7 @@
   FloatComplexMatrix retval (m, n);
   const float *smd = sm.data ();
   FloatComplex *rd = retval.fortran_vec ();
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < nel; i++)
     rd[i] = FloatComplex (smd[i], smd[nel+i]);
   return retval;
@@ -2590,6 +2594,7 @@
       return *this;
     }
 
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < a.length (); i++)
     elem (i, i) += a.elem (i, i);
 
@@ -2611,6 +2616,7 @@
       return *this;
     }
 
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < a.length (); i++)
     elem (i, i) -= a.elem (i, i);
 
@@ -3417,15 +3423,21 @@
 
   // The last column is not needed while using delta.
   float *delta = &retval(0, n-1);
-  for (octave_idx_type i = 0; i < m; i++)
-    delta[i] = (x2(i) - x1(i)) / (n - 1);
-
-  for (octave_idx_type j = 1; j < n-1; j++)
+  #pragma omp parallel
+  {
+    #pragma omp for
     for (octave_idx_type i = 0; i < m; i++)
-      retval(i, j) = x1(i) + j*delta[i];
-
-  for (octave_idx_type i = 0; i < m; i++)
-    retval(i, n-1) = x2(i);
+      delta[i] = (x2(i) - x1(i)) / (n - 1);
+
+    #pragma omp for
+    for (octave_idx_type j = 1; j < n-1; j++)
+      for (octave_idx_type i = 0; i < m; i++)
+        retval(i, j) = x1(i) + j*delta[i];
+
+    #pragma omp for
+    for (octave_idx_type i = 0; i < m; i++)
+      retval(i, n-1) = x2(i);
+  }
 
   return retval;
 }
diff -r c31fd42f9000 -r 9a78540d8b38 liboctave/array/fRowVector.cc
--- a/liboctave/array/fRowVector.cc	Sat Mar 09 18:42:43 2013 +0000
+++ b/liboctave/array/fRowVector.cc	Sun Mar 10 18:38:55 2013 -0300
@@ -261,10 +261,8 @@
 std::ostream&
 operator << (std::ostream& os, const FloatRowVector& a)
 {
-//  int field_width = os.precision () + 7;
-
   for (octave_idx_type i = 0; i < a.length (); i++)
-    os << " " /* setw (field_width) */ << a.elem (i);
+    os << " " << a.elem (i);
   return os;
 }
 
@@ -299,8 +297,10 @@
 
   float delta = (x2 - x1) / (n - 1);
   retval(0) = x1;
+  #pragma omp parallel for
   for (octave_idx_type i = 1; i < n-1; i++)
     retval(i) = x1 + i*delta;
+
   retval(n-1) = x2;
 
   return retval;
diff -r c31fd42f9000 -r 9a78540d8b38 liboctave/array/idx-vector.cc
--- a/liboctave/array/idx-vector.cc	Sat Mar 09 18:42:43 2013 +0000
+++ b/liboctave/array/idx-vector.cc	Sun Mar 10 18:38:55 2013 -0300
@@ -180,8 +180,11 @@
   if (step < 0 && len > 0)
     {
       idx.clear (1, len);
+
+      #pragma omp parallel for
       for (octave_idx_type i = 0; i < len; i++)
         idx.xelem (i) = len - 1 - i;
+
       return new idx_range_rep (start + (len - 1)*step, len, -step, DIRECT);
     }
   else
@@ -212,6 +215,8 @@
 idx_vector::idx_range_rep::as_array (void)
 {
   Array<octave_idx_type> retval (dim_vector (1, len));
+
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < len; i++)
     retval.xelem (i) = start + i*step;
 
@@ -515,6 +520,7 @@
         has[data[i]] = true;
 
       octave_idx_type new_len = 0;
+      #pragma omp parallel for reduction(+:new_len)
       for (octave_idx_type i = 0; i < ext; i++)
         new_len += has[i];
 
@@ -624,6 +630,7 @@
 idx_vector::idx_vector_rep::unconvert (void) const
 {
   Array<double> retval (orig_dims);
+  #pragma omp parallel for
   for (octave_idx_type i = 0; i < len; i++)
     retval.xelem (i) = data[i] + 1;
   return retval;
