// Test libstd linear algebra matinv(), solve(), svd() externs

#include "libstd.oah"

proc main()
{
    var a, udv, inv, b, x, zero;
    using namespace linalg;

    a = [[7,9,8],[3,4,5],[6,2,1]];
    b = [1,2,3];
    zero = 0 .reshape(3,3);

    // Test int=>float
    io::FltFormatChar = 'F';
    io::FieldWidth = 8;
    io::NumDigits = 5;
    inv = matinv(a);
    "Int inverse:\n", inv, '\n';
    udv = svd(a);
    "Int U:\n", udv[0], '\n';
    "Int D: ", udv[1], '\n';
    "Int V:\n", udv[2], '\n';
    x = solve(a, b);
    "Int X: ", x, '\n';
    inv = matinv(udv);
    "Int UDV inverse:\n", inv, '\n';
    x = solve(udv, b);
    "Int UDV X: ", x, '\n';

    // Next, test Double=>Double
    io::FieldWidth = 15;
    io::NumDigits = 12;
    inv = matinv(Double(a));
    "Double inverse:\n", inv, '\n';
    udv = svd(Double(a));
    "Double U:\n", udv[0], '\n';
    "Double D: ", udv[1], '\n';
    "Double V:\n", udv[2], '\n';
    x = solve(Double(a),Double(b));
    "Double X: ", x, '\n';
    inv = matinv(udv);
    "Double UDV inverse:\n", inv, '\n';
    x = solve(udv, b);
    "Double UDV X: ", x, '\n';

    // Test errors
    try { inv = matinv(zero); } catch (e) { "Zero inv: ", e, '\n'; }
    try { udv = svd(zero); } catch (e) { "Zero svd: ", e, '\n'; }
    try { x = solve(zero,b); } catch (e) { "Zero solve: ", e, '\n'; }
}