fixes a mistake and improves on the previous benchmarking script

benchmarking_nmatrix_solve_2.rb

require 'benchmark'
require 'nmatrix'

puts "------------ Without nmatrix-lapacke -----------"

n = 1000 
m = 100

puts "(A is #{n} x #{n}, B is #{n} x #{m})"

rand_mat = NMatrix.random([n, n], dtype: :float64)
# make the matrix well-conditioned
rand_mat = rand_mat + NMatrix.identity([n,n], dtype: :float64)
b = NMatrix.random([n, m], dtype: :float64)

# variables to hold the maximal errors of the solvers
err_general = -999.999
err_special = -999.999

puts "Solving an upper triangular linear system AX = B:"
a = rand_mat.triu
Benchmark.bm(10) do |bm|
  bm.report('general') do
    x = a.solve(b)
    err_general = (a.dot(x) - b).abs.to_flat_a.max
  end
  bm.report('upper_tri') do
    x = a.solve(b, form: :upper_tri)
    err_special = (a.dot(x) - b).abs.to_flat_a.max
  end
  puts "(Max. error: #{err_general} and #{err_special})"
end

puts "Solving a lower triangular linear system AX = B:"
a = rand_mat.tril
Benchmark.bm(10) do |bm|
  bm.report('general') do
    x = a.solve(b)
    err_general = (a.dot(x) - b).abs.to_flat_a.max
  end
  bm.report('lower_tri') do
    x = a.solve(b, form: :lower_tri)
    err_special = (a.dot(x) - b).abs.to_flat_a.max
  end
  puts "(Max. error: #{err_general} and #{err_special})"
end
 
puts "------------ With nmatrix-lapacke -----------"

n = 10000
m = 1

puts "(A is #{n} x #{n}, B is #{n} x #{m})"

require 'nmatrix/lapacke'

rand_mat = NMatrix.random([n, n], dtype: :float64)
# make the matrix well-conditioned
rand_mat = rand_mat + NMatrix.identity([n,n], dtype: :float64) * 4
b = NMatrix.random([n, m], dtype: :float64)

puts "Solving an upper triangular linear system AX = B:"
a = rand_mat.triu
Benchmark.bm(10) do |bm|
  bm.report('general') do
    x = a.solve(b)
    err_general = (a.dot(x) - b).abs.to_flat_a.max
  end
  bm.report('upper_tri') do
    x = a.solve(b, form: :upper_triangular)
    err_special = (a.dot(x) - b).abs.to_flat_a.max
  end
  puts "(Max. error: #{err_general} and #{err_special})"
end

puts "Solving a lower triangular linear system AX = B:"
a = rand_mat.tril
Benchmark.bm(10) do |bm|
  bm.report('general') do
    x = a.solve(b)
    err_general = (a.dot(x) - b).abs.to_flat_a.max
  end
  bm.report('lower_tri') do
    x = a.solve(b, form: :lower_triangular)
    err_special = (a.dot(x) - b).abs.to_flat_a.max
  end
  puts "(Max. error: #{err_general} and #{err_special})"
end

puts "Solving a linear system AX = B with positive definite A:"
a = rand_mat.transpose.dot rand_mat
Benchmark.bm(10) do |bm|
  bm.report('general') do
    x = a.solve(b)
    err_general = (a.dot(x) - b).abs.to_flat_a.max
  end
  bm.report('pos_def') do
    x = a.solve(b, form: :positive_definite)
    err_special = (a.dot(x) - b).abs.to_flat_a.max
  end
  puts "(Max. error: #{err_general} and #{err_special})"
end