--- a/doc/interpreter/basics.txi
+++ b/doc/interpreter/basics.txi
@@ -247,11 +247,13 @@
 beep_on_error                   = true
 confirm_recursive_rmdir         = false
 crash_dumps_octave_core         = false
-save_default_options            = "-mat-binary"
+disable_diagonal_matrix         = true
+disable_permutation_matrix      = true
 fixed_point_format              = true
 history_timestamp_format_string = "%%-- %D %I:%M %p --%%"
 page_screen_output              = false
 print_empty_dimensions          = false
+save_default_options            = "-mat-binary"
 struct_levels_to_print          = 0
 @end group
 @end example

--- a/doc/interpreter/diagperm.txi
+++ b/doc/interpreter/diagperm.txi
@@ -59,7 +59,13 @@
 Octave provides special treatment of real and complex rectangular diagonal
 matrices, as well as permutation matrices.  They are stored as special objects,
 using efficient storage and algorithms, facilitating writing both readable and
-efficient matrix algebra expressions in the Octave language.
+efficient matrix algebra expressions in the Octave language.  The
+special treatment may be disabled by using the functions
+@dfn{disable_permutation_matrix} and @dfn{disable_diagonal_matrix}
+
+@DOCSTRING(disable_permutation_matrix)
+
+@DOCSTRING(disable_diagonal_matrix)
 
 @menu
 * Creating Diagonal Matrices::

--- a/libinterp/octave-value/ov.cc
+++ b/libinterp/octave-value/ov.cc
@@ -95,6 +95,14 @@
 // make the grow_size large.
 DEFINE_OCTAVE_ALLOCATOR2(octave_value, 1024);
 
+// If TRUE, don't create special diagonal matrix objects.
+
+static bool Vdisable_diagonal_matrix = false;
+
+// If TRUE, don't create special permutation matrix objects.
+
+static bool Vdisable_permutation_matrix = false;
+
 // FIXME
 
 // Octave's value type.
@@ -674,37 +682,49 @@
 }
 
 octave_value::octave_value (const DiagArray2<double>& d)
-  : rep (new octave_diag_matrix (d))
+  : rep (Vdisable_diagonal_matrix
+         ? dynamic_cast<octave_base_value *> (new octave_matrix (Matrix (d)))
+         : dynamic_cast<octave_base_value *> (new octave_diag_matrix (d)))
 {
   maybe_mutate ();
 }
 
 octave_value::octave_value (const DiagArray2<float>& d)
-  : rep (new octave_float_diag_matrix (d))
+  : rep (Vdisable_diagonal_matrix
+         ? dynamic_cast<octave_base_value *> (new octave_float_matrix (FloatMatrix (d)))
+         : dynamic_cast<octave_base_value *> (new octave_float_diag_matrix (d)))
 {
   maybe_mutate ();
 }
 
 octave_value::octave_value (const DiagArray2<Complex>& d)
-  : rep (new octave_complex_diag_matrix (d))
+  : rep (Vdisable_diagonal_matrix
+         ? dynamic_cast<octave_base_value *> (new octave_complex_matrix (ComplexMatrix (d)))
+         : dynamic_cast<octave_base_value *> (new octave_complex_diag_matrix (d)))
 {
   maybe_mutate ();
 }
 
 octave_value::octave_value (const DiagArray2<FloatComplex>& d)
-  : rep (new octave_float_complex_diag_matrix (d))
+  : rep (Vdisable_diagonal_matrix
+         ? dynamic_cast<octave_base_value *> (new octave_float_complex_matrix (FloatComplexMatrix (d)))
+         : dynamic_cast<octave_base_value *> (new octave_float_complex_diag_matrix (d)))
 {
   maybe_mutate ();
 }
 
 octave_value::octave_value (const DiagMatrix& d)
-  : rep (new octave_diag_matrix (d))
+  : rep (Vdisable_diagonal_matrix
+         ? dynamic_cast<octave_base_value *> (new octave_matrix (Matrix (d)))
+         : dynamic_cast<octave_base_value *> (new octave_diag_matrix (d)))
 {
   maybe_mutate ();
 }
 
 octave_value::octave_value (const FloatDiagMatrix& d)
-  : rep (new octave_float_diag_matrix (d))
+  : rep (Vdisable_diagonal_matrix
+         ? dynamic_cast<octave_base_value *> (new octave_float_matrix (FloatMatrix (d)))
+         : dynamic_cast<octave_base_value *> (new octave_float_diag_matrix (d)))
 {
   maybe_mutate ();
 }
@@ -782,13 +802,17 @@
 }
 
 octave_value::octave_value (const ComplexDiagMatrix& d)
-  : rep (new octave_complex_diag_matrix (d))
+  : rep (Vdisable_diagonal_matrix
+         ? dynamic_cast<octave_base_value *> (new octave_complex_matrix (ComplexMatrix (d)))
+         : dynamic_cast<octave_base_value *> (new octave_complex_diag_matrix (d)))
 {
   maybe_mutate ();
 }
 
 octave_value::octave_value (const FloatComplexDiagMatrix& d)
-  : rep (new octave_float_complex_diag_matrix (d))
+  : rep (Vdisable_diagonal_matrix
+         ? dynamic_cast<octave_base_value *> (new octave_float_complex_matrix (FloatComplexMatrix (d)))
+         : dynamic_cast<octave_base_value *> (new octave_float_complex_diag_matrix (d)))
 {
   maybe_mutate ();
 }
@@ -818,7 +842,9 @@
 }
 
 octave_value::octave_value (const PermMatrix& p)
-  : rep (new octave_perm_matrix (p))
+  : rep (Vdisable_permutation_matrix
+         ? dynamic_cast<octave_base_value *> (new octave_matrix (Matrix (p)))
+         : dynamic_cast<octave_base_value *> (new octave_perm_matrix (p)))
 {
   maybe_mutate ();
 }
@@ -3148,3 +3174,73 @@
 %!error is_dq_string ()
 %!error is_dq_string ("foo", 2)
 */
+
+DEFUN (disable_permutation_matrix, args, nargout,
+       "-*- texinfo -*-\n\
+@deftypefn  {Built-in Function} {@var{val} =} disable_permutation_matrix ()\n\
+@deftypefnx {Built-in Function} {@var{old_val} =} disable_permutation_matrix (@var{new_val})\n\
+@deftypefnx {Built-in Function} {} disable_permutation_matrix (@var{new_val}, \"local\")\n\
+Query or set the internal variable that controls whether permutation\n\
+matrices are stored in a special space-efficient format.  The default\n\
+value is true.  If this option is disabled Octave will store permutation\n\
+matrices as full matrices.\n\
+\n\
+When called from inside a function with the @qcode{\"local\"} option, the\n\
+variable is changed locally for the function and any subroutines it calls.\n\
+The original variable value is restored when exiting the function.\n\
+@end deftypefn")
+{
+  return SET_INTERNAL_VARIABLE (disable_permutation_matrix);
+}
+
+/*
+%!function p = __test_dpm__ (dpm)
+%!  disable_permutation_matrix (dpm, "local");
+%!  [~, ~, p] = lu ([1,2;3,4]);
+%!endfunction
+
+%!assert (typeinfo (__test_dpm__ (false)), "permutation matrix");
+%!assert (typeinfo (__test_dpm__ (true)), "matrix");
+*/
+
+DEFUN (disable_diagonal_matrix, args, nargout,
+       "-*- texinfo -*-\n\
+@deftypefn  {Built-in Function} {@var{val} =} disable_diagonal_matrix ()\n\
+@deftypefnx {Built-in Function} {@var{old_val} =} disable_diagonal_matrix (@var{new_val})\n\
+@deftypefnx {Built-in Function} {} disable_diagonal_matrix (@var{new_val}, \"local\")\n\
+Query or set the internal variable that controls whether diagonal\n\
+matrices are stored in a special space-efficient format.  The default\n\
+value is true.  If this option is disabled Octave will store diagonal\n\
+matrices as full matrices.\n\
+\n\
+When called from inside a function with the @qcode{\"local\"} option, the\n\
+variable is changed locally for the function and any subroutines it calls.\n\
+The original variable value is restored when exiting the function.\n\
+@end deftypefn")
+{
+  return SET_INTERNAL_VARIABLE (disable_diagonal_matrix);
+}
+
+/*
+%!function [x, xi, fx, fxi] = __test_ddm__ (ddm)
+%!  disable_diagonal_matrix (ddm, "local");
+%!  x = eye (2);
+%!  xi = x*i;
+%!  fx = single (x);
+%!  fxi = single (xi);
+%!endfunction
+
+%!shared x, xi, fx, fxi
+%!  [x, xi, fx, fxi] = __test_ddm__ (false);
+%!assert (typeinfo (x), "diagonal matrix");
+%!assert (typeinfo (xi), "complex diagonal matrix");
+%!assert (typeinfo (fx), "float diagonal matrix");
+%!assert (typeinfo (fxi), "float complex diagonal matrix");
+
+%!shared x, xi, fx, fxi
+%!  [x, xi, fx, fxi] = __test_ddm__ (true);
+%!assert (typeinfo (x), "matrix");
+%!assert (typeinfo (xi), "complex matrix");
+%!assert (typeinfo (fx), "float matrix");
+%!assert (typeinfo (fxi), "float complex matrix");
+*/

--- a/libinterp/octave.cc
+++ b/libinterp/octave.cc
@@ -486,11 +486,13 @@
   Fbeep_on_error (octave_value (true));
   Fconfirm_recursive_rmdir (octave_value (false));
   Fcrash_dumps_octave_core (octave_value (false));
-  Fsave_default_options (octave_value ("-mat-binary"));
+  Fdisable_diagonal_matrix (octave_value (true));
+  Fdisable_permutation_matrix (octave_value (true));
   Ffixed_point_format (octave_value (true));
   Fhistory_timestamp_format_string (octave_value ("%%-- %D %I:%M %p --%%"));
   Fpage_screen_output (octave_value (false));
   Fprint_empty_dimensions (octave_value (false));
+  Fsave_default_options (octave_value ("-mat-binary"));
   Fstruct_levels_to_print (octave_value (0));
 
   disable_warning ("Octave:abbreviated-property-match");

--- a/liboctave/array/CMatrix.cc
+++ b/liboctave/array/CMatrix.cc
@@ -289,6 +289,20 @@
     elem (i, i) = a.elem (i, i);
 }
 
+ComplexMatrix::ComplexMatrix (const MDiagArray2<double>& a)
+  : ComplexNDArray (a.dims (), 0.0)
+{
+  for (octave_idx_type i = 0; i < a.length (); i++)
+    elem (i, i) = a.elem (i, i);
+}
+
+ComplexMatrix::ComplexMatrix (const DiagArray2<double>& a)
+  : ComplexNDArray (a.dims (), 0.0)
+{
+  for (octave_idx_type i = 0; i < a.length (); i++)
+    elem (i, i) = a.elem (i, i);
+}
+
 ComplexMatrix::ComplexMatrix (const ComplexRowVector& rv)
   : ComplexNDArray (rv)
 {
@@ -306,6 +320,20 @@
     elem (i, i) = a.elem (i, i);
 }
 
+ComplexMatrix::ComplexMatrix (const MDiagArray2<Complex>& a)
+  : ComplexNDArray (a.dims (), 0.0)
+{
+  for (octave_idx_type i = 0; i < a.length (); i++)
+    elem (i, i) = a.elem (i, i);
+}
+
+ComplexMatrix::ComplexMatrix (const DiagArray2<Complex>& a)
+  : ComplexNDArray (a.dims (), 0.0)
+{
+  for (octave_idx_type i = 0; i < a.length (); i++)
+    elem (i, i) = a.elem (i, i);
+}
+
 // FIXME: could we use a templated mixed-type copy function here?
 
 ComplexMatrix::ComplexMatrix (const boolMatrix& a)

--- a/liboctave/array/CMatrix.h
+++ b/liboctave/array/CMatrix.h
@@ -75,12 +75,20 @@
 
   explicit ComplexMatrix (const DiagMatrix& a);
 
+  explicit ComplexMatrix (const MDiagArray2<double>& a);
+
+  explicit ComplexMatrix (const DiagArray2<double>& a);
+
   explicit ComplexMatrix (const ComplexRowVector& rv);
 
   explicit ComplexMatrix (const ComplexColumnVector& cv);
 
   explicit ComplexMatrix (const ComplexDiagMatrix& a);
 
+  explicit ComplexMatrix (const MDiagArray2<Complex>& a);
+
+  explicit ComplexMatrix (const DiagArray2<Complex>& a);
+
   explicit ComplexMatrix (const boolMatrix& a);
 
   explicit ComplexMatrix (const charMatrix& a);

--- a/liboctave/array/dMatrix.cc
+++ b/liboctave/array/dMatrix.cc
@@ -262,6 +262,20 @@
     elem (i, i) = a.elem (i, i);
 }
 
+Matrix::Matrix (const MDiagArray2<double>& a)
+  : NDArray (a.dims (), 0.0)
+{
+  for (octave_idx_type i = 0; i < a.length (); i++)
+    elem (i, i) = a.elem (i, i);
+}
+
+Matrix::Matrix (const DiagArray2<double>& a)
+  : NDArray (a.dims (), 0.0)
+{
+  for (octave_idx_type i = 0; i < a.length (); i++)
+    elem (i, i) = a.elem (i, i);
+}
+
 Matrix::Matrix (const PermMatrix& a)
   : NDArray (a.dims (), 0.0)
 {

--- a/liboctave/array/dMatrix.h
+++ b/liboctave/array/dMatrix.h
@@ -70,6 +70,10 @@
 
   explicit Matrix (const DiagMatrix& a);
 
+  explicit Matrix (const MDiagArray2<double>& a);
+
+  explicit Matrix (const DiagArray2<double>& a);
+
   explicit Matrix (const PermMatrix& a);
 
   explicit Matrix (const boolMatrix& a);

--- a/liboctave/array/fCMatrix.cc
+++ b/liboctave/array/fCMatrix.cc
@@ -290,6 +290,20 @@
     elem (i, i) = a.elem (i, i);
 }
 
+FloatComplexMatrix::FloatComplexMatrix (const MDiagArray2<float>& a)
+  : FloatComplexNDArray (a.dims (), 0.0)
+{
+  for (octave_idx_type i = 0; i < a.length (); i++)
+    elem (i, i) = a.elem (i, i);
+}
+
+FloatComplexMatrix::FloatComplexMatrix (const DiagArray2<float>& a)
+  : FloatComplexNDArray (a.dims (), 0.0)
+{
+  for (octave_idx_type i = 0; i < a.length (); i++)
+    elem (i, i) = a.elem (i, i);
+}
+
 FloatComplexMatrix::FloatComplexMatrix (const FloatComplexRowVector& rv)
   : FloatComplexNDArray (rv)
 {
@@ -307,6 +321,20 @@
     elem (i, i) = a.elem (i, i);
 }
 
+FloatComplexMatrix::FloatComplexMatrix (const MDiagArray2<FloatComplex>& a)
+  : FloatComplexNDArray (a.dims (), 0.0)
+{
+  for (octave_idx_type i = 0; i < a.length (); i++)
+    elem (i, i) = a.elem (i, i);
+}
+
+FloatComplexMatrix::FloatComplexMatrix (const DiagArray2<FloatComplex>& a)
+  : FloatComplexNDArray (a.dims (), 0.0)
+{
+  for (octave_idx_type i = 0; i < a.length (); i++)
+    elem (i, i) = a.elem (i, i);
+}
+
 // FIXME: could we use a templated mixed-type copy function
 // here?
 

--- a/liboctave/array/fCMatrix.h
+++ b/liboctave/array/fCMatrix.h
@@ -78,12 +78,20 @@
 
   explicit FloatComplexMatrix (const FloatDiagMatrix& a);
 
+  explicit FloatComplexMatrix (const MDiagArray2<float>& a);
+
+  explicit FloatComplexMatrix (const DiagArray2<float>& a);
+
   explicit FloatComplexMatrix (const FloatComplexRowVector& rv);
 
   explicit FloatComplexMatrix (const FloatComplexColumnVector& cv);
 
   explicit FloatComplexMatrix (const FloatComplexDiagMatrix& a);
 
+  explicit FloatComplexMatrix (const MDiagArray2<FloatComplex>& a);
+
+  explicit FloatComplexMatrix (const DiagArray2<FloatComplex>& a);
+
   explicit FloatComplexMatrix (const boolMatrix& a);
 
   explicit FloatComplexMatrix (const charMatrix& a);

--- a/liboctave/array/fMatrix.cc
+++ b/liboctave/array/fMatrix.cc
@@ -265,6 +265,20 @@
     elem (i, i) = a.elem (i, i);
 }
 
+FloatMatrix::FloatMatrix (const MDiagArray2<float>& a)
+  : FloatNDArray (a.dims (), 0.0)
+{
+  for (octave_idx_type i = 0; i < a.length (); i++)
+    elem (i, i) = a.elem (i, i);
+}
+
+FloatMatrix::FloatMatrix (const DiagArray2<float>& a)
+  : FloatNDArray (a.dims (), 0.0)
+{
+  for (octave_idx_type i = 0; i < a.length (); i++)
+    elem (i, i) = a.elem (i, i);
+}
+
 FloatMatrix::FloatMatrix (const PermMatrix& a)
   : FloatNDArray (a.dims (), 0.0)
 {

--- a/liboctave/array/fMatrix.h
+++ b/liboctave/array/fMatrix.h
@@ -70,6 +70,10 @@
 
   explicit FloatMatrix (const FloatDiagMatrix& a);
 
+  explicit FloatMatrix (const MDiagArray2<float>& a);
+
+  explicit FloatMatrix (const DiagArray2<float>& a);
+
   explicit FloatMatrix (const PermMatrix& a);
 
   explicit FloatMatrix (const boolMatrix& a);