Seems like another addition for Matlab compatibility will be to map new options for the type of eigenvalues to return to their old two-letter codes based on symmetry and real/complex nature of input.

sigma                      Description                                    sigma (R2017a and earlier)
----------------------------------------------------------------------------------------------------

scalar (real or complex)   The eigenvalues closest to the number sigma.   No change

'largestabs' (default)     Largest magnitude.                             'lm'

'smallestabs'              Smallest magnitude. Same as sigma = 0.         'sm'

'largestreal'              Largest real.                                  'lr', 'la'

'smallestreal'             Smallest real.                                 'sr', 'sa'

'bothendsreal'             Both ends, with k/2 values with largest and    'be'
                           smallest real part respectively
                           (one more from high end if k is odd).

Adding David Bateman to CC list since he is deeply familiar with the eigs function and the interface to ARPACK.

The issue of the ordering of singular values depends on the underlying function eigs().  That function does not guarantee a particular ordering.  However, Matlab has changed behavior relatively recently (R2017b).  They now state

Behavior changed in R2017b

    Changes to sorting order of output

    eigs now sorts the output according to the value of sigma. For example, the command eigs(A,k,'largestabs') produces k eigenvalues sorted in descending order by magnitude.

    Previously, the sorting order of the output produced by eigs was not guaranteed.

    Reproducibility

    Calling eigs multiple times in succession now produces the same result. Set 'StartVector' to a random vector to change this behavior.

    Display

    A display value of 2 no longer returns timing information. Instead, eigs treats a value of 2 the same as a value of 1. Also, the messages shown by the 'Display' option have changed. The new messages show the residual in each iteration, instead of the Ritz values.

Re-titling bug report to reflect underlying issue.

I filed a bug report about eigs() not returning strictly real values for symmetric inputs (bug #57196).  That would fix this, but in the meantime, it is appropriate to apply the patch from comment #9.

For the second question, I think we need to check the code.  In svds.m there is

    if (nargout > 1)
      [V, s, flag] = eigs ([sparse(m,m), b; b', sparse(n,n)],
                           b_k, b_sigma, b_opts);
      s = diag (s);
    else
      s = eigs ([sparse(m,m), b; b', sparse(n,n)], b_k, b_sigma, b_opts);
    endif

so we are calling eigs with a different function prototype in each case.  I seem to recall that if only the singular values are requested we may perform a sort on them to ensure uniform ordering.  But, we can't do that arbitrarily when someone has requested the full decomposition because it would also require changing the U and V matrices to ensure that the decomposition

A = U * S * V'

still holds. 

Now this makes sense, same for default hg_id 959a63d1ede6:

x = randn(10, 5) ;
y = randn(10, 5) ;
z = complex(x, y) ;

[~, Dr] = svds(x);
svds_x = [diag(Dr), svds(x), svd(x)]
[~, Dc] = svds(z);
svd_z = [diag(Dc), svds(z), svd(z)]

isreal(Dr)
arrayfun (@isreal, svd_z)

The output

svds_x =

   5.1420   5.1420   5.1420
   3.0993   3.0993   3.0993
   2.7733   2.7733   2.7733
   2.1100   2.1100   2.1100
   1.4318   1.4318   1.4318

svd_z =

   2.6499 + 0.0000i   7.3784 - 0.0000i   7.3784 + 0.0000i
   3.8673 + 0.0000i   5.7718 - 0.0000i   5.7718 + 0.0000i
   4.3697 - 0.0000i   4.3697 - 0.0000i   4.3697 + 0.0000i
   5.7718 - 0.0000i   3.8673 + 0.0000i   3.8673 + 0.0000i
   7.3784 + 0.0000i   2.6499 - 0.0000i   2.6499 + 0.0000i

ans = 1
ans =

  0  0  1
  0  0  1
  0  0  1
  0  0  1
  0  0  1

Thus two problems:

1. The order of the singular values depends on the number of output arguments

2. The singular values are complex values.

Agree to comment #1 and comment #2.  Nothing suspicious using

>> __blas_version__ ()
ans = OpenBLAS (config: OpenBLAS 0.3.6 NO_LAPACKE DYNAMIC_ARCH NO_AFFINITY Haswell MAX_THREADS=48)
>> __lapack_version__ ()
ans = Linear Algebra PACKage Version 3.8.0

on openSUSE 15.1 using the flatpak installation or a custom build.

Same suspicion as comment #3: What is your output for the functions above and did you customize your Octave somehow (please explain how you installed Octave)?

Also works for me.  In case it was something related to random numbers, I tried several thousand runs with

for i = 1:1000
  x = randn(10, 5) ;
  y = randn(10, 5) ;
  z = complex(x, y) ;

  [~, Dr] = svd(x) ; diag(Dr);
  svd(x);
  [~, Dc] = svd(z) ; diag(Dc);
  svd(z);

  assert (isreal(Dr))
  assert (isreal(Dc))
endfor

It always passes.

x = randn(10, 5) ;
y = randn(10, 5) ;
z = complex(x, y) ;

[~, Dr] = svd(x) ; diag(Dr)
svd(x)
[~, Dc] = svd(z) ; diag(Dc)
svd(z)

isreal(Dr)
isreal(Dc)

This shows that in the complex case, the ordering depends on having explicit output (Dz) or not. This is inconsistent and also different from Matlab.

svd should also return real and not complex eigenvalues.