octave-maintainers
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Help with eigs


From: Rik
Subject: Help with eigs
Date: Mon, 7 Jan 2019 13:48:45 -0800

Marco,

I've been working through the static analyzer warnings about the Octave
code base and I ran in to a bit of trouble with eigs.

A sample warning is

__eigs__.cc (369) V550     An odd precise comparison: tmp.double_value() !=
0.0. It's probably better to use a comparison with defined precision:
fabs(A - B) > Epsilon.

for which the code is

tmp = map.getfield ("cholB");
if (tmp.is_defined ())
  cholB = tmp.double_value () != 0.0;

Checking the documentation for eigs

'cholB'
     Flag if 'chol (B)' is passed rather than B.  The default is
     false.

So, I thought the input was just a true/false and that I could change the
code to

  cholB = tmp.bool_value ();

But now when I run the BIST tests I get

octave:1> test eigs
***** testif HAVE_ARPACK, HAVE_UMFPACK
 A = toeplitz (sparse (1:10));
 B = toeplitz (sparse ([1, 1], [1, 2], [2, 1], 1, 10));
 R = chol (B);
 opts.cholB = R;
 [v, d] = eigs (A, R, 4, "lm", opts);
 for i = 1:4
   assert (A * v(:,i), d(i, i) * B * v(:,i), 1e-12)
 endfor
!!!!! test failed
octave_base_value::bool_value(): wrong type argument 'sparse matrix'

If I execute the commands individually I find that opts.cholB is not
true/false, but a sparse matrix.

Compressed Column Sparse (rows = 10, cols = 10, nnz = 19 [19%])

  (1, 1) ->  1.4142
  (1, 2) ->  0.70711
  (2, 2) ->  1.2247
  (2, 3) ->  0.81650
  (3, 3) ->  1.1547
  (3, 4) ->  0.86603
  (4, 4) ->  1.1180
  (4, 5) ->  0.89443
  (5, 5) ->  1.0954
  (5, 6) ->  0.91287
  (6, 6) ->  1.0801
  (6, 7) ->  0.92582
  (7, 7) ->  1.0690
  (7, 8) ->  0.93541
  (8, 8) ->  1.0607
  (8, 9) ->  0.94281
  (9, 9) ->  1.0541
  (9, 10) ->  0.94868
  (10, 10) ->  1.0488

Can you describe what is supposed to take place here and how to modify the
BIST tests?  Is it as simple as "opts.cholB = true;"?

--Rik





reply via email to

[Prev in Thread] Current Thread [Next in Thread]