About diagonal matrices

[John W. Eaton]

On 20-Feb-2009, José Luis García Pallero wrote:

| Playing with the new diag() function in octave 3.1.52 I can see some
| extrange behaviour (or I don't understand some particular cases). For
| example, for the matrix
| a = diag([2 3 4])
| 
| This is a diagonal matrix: typeinfo(a) -> diagonal matrix
| But if I add a scalar: typeinfo(4+a) -> matrix
| If I multiply by a scalar: typeinfo(4*a) -> diagonal matrix
| But if I multiply by a scalar in element-by-element form: typeinfo(4.*a) ->
| matrix
| If I divide by a scalar in element-by-element form: typeinfo(4./a) -> matrix
| For a power operation: typeinfo(a^2) -> diagonal matrix
| But element-by-element: typeinfo(a.^2) -> matrix
| For functions: typeinfo(sqrt(a)) -> matrix (the same for sin(), cos(), etc.)
| 
| Has any reason for this behaviour?

A definition is needed in each case for oprators that receive special
treatment.  So I guess some have just not been implemented yet.

But I think there are also some bugs.  Shouldn't the following return
full matrices with the zero elements replaced by NaN (or -0, in the
case of dividing by -Inf)?

  diag ([1,2,3]) / 0
  diag ([1,2,3]) / NaN
  diag ([1,2,3]) / -Inf
  diag ([1,2,3]) * NaN
  diag ([1,2,3]) * Inf

etc.  Maybe the simplest way to implement this would be to do
something like

  tmp = 0 / scalar_value;
  if (tmp != 0 || lo_ieee_signbit (signbit))
    convert to full matrix, replacing zeros with tmp;

in these methods.

If you want to discuss this further, then maybe we should move the
discussion to the maintainers list.

jwe