test.d

/+dub.sdl:
dependency "mir-lapack" version="~>1.2.1"
dependency "mir-algorithm" version="~>3.5.0-alpha12"
+/
import mir.ndslice.slice;
import mir.rc.array;
import mir.ndslice.allocation;
import std.traits: Unqual;
import mir.internal.utility: isComplex, realType;
import mir.ndslice.dynamic: transposed;
import mir.ndslice.topology: as, canonical;
import mir.lapack;
import mir.utility: min;
import mir.algorithm.iteration: equal;
import std.math:sqrt,pow;

size_t mldivide (T, SliceKind kindA, SliceKind kindB)(
    Slice!(const(T)*, 2, kindA) a,
    Slice!(const(T)*, 2, kindB) b,
    ref Slice!(T*, 2) x
 )
    if (isComplex!T)
in
{
    assert(a.length!0 == b.length!0);
}
out
{
    assert(x.length!0 == a.length!1);
    assert(x.length!1 == b.length!1);
}
do
{
    import std.stdio:writeln;
    writeln("Complex");
    auto rcat = a.transposed.as!T.slice;
    // writeln(rcat);
    auto at = rcat.canonical;
    // writeln(at);
    auto rcbt = b.transposed.as!T.slice;
    auto bt = rcbt.canonical;
    // writeln(bt);
    size_t info;
    writeln("correct");
    size_t liwork = void;
    size_t lrwork = void;
    writeln("at ", at.length!0, ", ", at.length!1);
    writeln("bt ", bt.length!0, ", ", bt.length!1);
    assert(at.length!1 == bt.length!1);
    auto lwork = gelsd_wq(at, bt, liwork, lrwork);
    writeln(liwork, " ", lrwork, " ", lwork);
    auto s = min(at.length!0, at.length!1).mininitRcslice!(realType!T);
    auto work = lwork.slice!T;
    auto work_lightScope = work.lightScope;
    auto iwork = liwork.slice!lapackint;
    auto iwork_lightScope = iwork.lightScope;
    auto rwork = lrwork.slice!(realType!T);
    auto rwork_lightScope = rwork.lightScope;
    size_t rank = void;
    realType!T rcond = -1;

    info = gelsd!T(at.lightScope, bt.lightScope, s.lightScope, rcond, rank, work_lightScope, rwork_lightScope, iwork_lightScope);
    writeln("a ", at);
    writeln("b ", bt);
    writeln("s ", s);
    writeln("rank ", rank);
    writeln("info ", info);
    //info > 0 means that many components failed to converge
        
    writeln("exit static if");
    writeln(at.length!0, ", ", at.length!1);
    bt = bt[0 .. $, 0 .. at.length!0];
    writeln("bt crop ", bt);
    writeln("exit if");
    // x = mininitRcslice!T(bt.length!1, bt.length!0);
    writeln("x alloc");
    x = bt.transposed.slice;
    // writeln(x);
    return info;
}

void main() {
    import std.stdio:writeln;
    auto ca = [
        -0.57 + 0.0i,  -1.28 + 0.0i,  -0.39 + 0.0i,   0.25 + 0.0i,
        -1.93 + 0.0i,   1.08 + 0.0i,  -0.31 + 0.0i,  -2.14 + 0.0i,
        2.30 + 0.0i,   0.24 + 0.0i,   0.40 + 0.0i,  -0.35 + 0.0i,
        -1.93 + 0.0i,   0.64 + 0.0i,  -0.66 + 0.0i,   0.08 + 0.0i,
        0.15 + 0.0i,   0.30 + 0.0i,   0.15 + 0.0i,  -2.13 + 0.0i,
        -0.02 + 0.0i,   1.03 + 0.0i,  -1.43 + 0.0i,   0.50 + 0.0i,
    ].sliced(6, 4); //.rcslice;

    auto cx = [
        1.5339 + 0.0i,
        1.8707 + 0.0i,
        -1.5241 + 0.0i,
        0.0392 + 0.0i
    ].sliced(4, 1); //.rcslice;

    auto cb = [
        -2.67 + 0.0i,
        -0.55 + 0.0i,
        3.34 + 0.0i,
        -0.77 + 0.0i,
        0.48 + 0.0i,
        4.10 + 0.0i,
    ].sliced(6, 1); //.rcslice;

    Slice!(cdouble*, 2) cres;

    auto cinfo = mldivide!(cdouble)(ca.lightConst, cb.lightConst, cres);
    writeln("cinfo", cinfo);
    assert(cinfo == 0);
    assert(equal!((a,b)=>norm!cdouble(a-b)<5e-5)(cres, cx));
}


realType!T norm(T)(in T z)
{
    return sqrt(pow(z.re,2.0)+pow(z.im,2.0));
}