clear all;
rand("seed", 0);

tol = sqrt(eps);

for j=1:2
  for i=1:100
    A = rand(10, 10);
    if j == 1
      B = eye(columns(A));
    else
      B = rand(rows(A), columns(A));
    endif
    A *= A.';
    B *= B.';
    opts.isreal = 1;
    opts.issym = 1;
    [v, lambda] = eigs(A, B, 3, "sm");
    [v2, lambda2] = eigs(@(x) A \ x, columns(A), B, 3, "sm", opts);
    assert(lambda, lambda2, tol*max(max(abs(lambda))));
    for i=1:columns(v)
      assert(A * v(:, i), lambda(i,i) * B * v(:,i), tol * norm(A * v(:,i)));
    endfor
    for i=1:columns(v2)
      assert(A * v2(:, i), lambda2(i,i) * B * v2(:, i), tol * norm(A * v2(:,i)));
    endfor
  endfor
endfor
