// A, B, C = matrix types (B,C may be NO_MATRIX)
// a, b, c are names of corresponding matrices
// s is not indexed (usually a scalar)
// m,n are fixed size (-1 if dynamic)
// #define FUNCTION(aij,bij,cij,s)
// Optionally define TRANS_A, TRANS_B, TRANS_C to access those matrices transposed.
#ifdef TRANS_A
#undef TRANS_A
#define TRANS_A(x,y) y,x
#else
#define TRANS_A(x,y) x,y
#endif

#ifdef TRANS_B
#undef TRANS_B
#define TRANS_B(x,y) y,x
#else
#define TRANS_B(x,y) x,y
#endif

#ifdef TRANS_C
#undef TRANS_C
#define TRANS_C(x,y) y,x
#else
#define TRANS_C(x,y) x,y
#endif

{
    constexpr int ac=access_flags(A::type)&access_flags(B::type)&access_flags(C::type);
    if constexpr(A::type==mat_zero){}
    else if constexpr(A::type==mat_id){}
    else if constexpr(A::type==mat_scale)
    {
        // Careful not to access out of bounds.
        if constexpr(m>0){FUNCTION(a.x,diag(b,0),diag(c,0),s);}
    }
    else if constexpr(A::type==mat_diag)
    {
        // If we are small enough to hard-code.
        if constexpr(m<=size_limit)
        {
            []<int... i>(auto&& a,auto&& b,auto&& c,auto&& s,std::integer_sequence<int,i...>)
            {
                ((FUNCTION(
                        a(i,i),
                        (get<TRANS_B(i,i)>(b)),
                        (get<TRANS_C(i,i)>(c)),s)),...);
            }(a,b,c,s,std::make_integer_sequence<int,m>{});
        }
        else // Too big, generate dynamically.
        {
            for(int i=0;i<m;i++)
                FUNCTION(a(i,i),diag(b,i),diag(c,i),s);
        }
    }
    else if constexpr(A::type==mat_sym)
    {
        // If we are small enough to hard-code.
        if constexpr(m<=size_limit)
        {
            []<int... i>(auto&& a,auto&& b,auto&& c,auto&& s,std::integer_sequence<int,i...>)
            {
                ((FUNCTION(
                        a.x[i],
                        (get<TRANS_B(sym_i(m,i),sym_j(m,i))>(b)),
                        (get<TRANS_C(sym_i(m,i),sym_j(m,i))>(c)),s)),...);
            }(a,b,c,s,std::make_integer_sequence<int,m*(m+1)/2>{});
        }
        else if constexpr(ac&1) // we can use lower
        {
            for(int j=0;j<m;j++)
                for(int i=j;i<m;i++)
                    FUNCTION(a.Element_Lower(TRANS_A(i,j)),lower(b,TRANS_B(i,j)),lower(c,TRANS_C(i,j)),s);
        }
        else // must use diag/lower_
        {
            for(int j=0;j<m;j++)
                for(int i=j+1;i<m;i++)
                    FUNCTION(a.Element_Lower(TRANS_A(i,j)),lower_(b,TRANS_B(i,j)),lower_(c,TRANS_C(i,j)),s);
            for(int i=0;i<m;i++)
                FUNCTION(a.Element_Lower(i,i),diag(b,i),diag(c,i),s);
        }
    }
    else // mat_gen
    {
        // If we are small enough to hard-code.
        if constexpr(m>=0 && m<=size_limit && n<=size_limit)
        {
            []<int... i>(auto&& a,auto&& b,auto&& c,auto&& s,std::integer_sequence<int,i...>)
            {
                ((FUNCTION(
                        (a(TRANS_A(i%m,i/m))),
                        (get<TRANS_B(i%m,i/m)>(b)),
                        (get<TRANS_C(i%m,i/m)>(c)),s)),...);
            }(a,b,c,s,std::make_integer_sequence<int,m*n>{});
        }
        else if constexpr(ac&2) // we can use get
        {
            int dyn_m=a.Rows(),dyn_n=a.Columns();
            for(int j=0;j<dyn_n;j++)
                for(int i=0;i<dyn_m;i++)
                    FUNCTION(a(TRANS_A(i,j)),get(b,TRANS_B(i,j)),get(c,TRANS_C(i,j)),s);
        }
        else if constexpr(ac&1) // we can use lower/upper
        {
            int dyn_m=a.Rows(),dyn_n=a.Columns();
            assert(dyn_m==dyn_n);
            for(int j=0;j<dyn_n;j++)
            {
                for(int i=0;i<j;i++)
                    FUNCTION(a(TRANS_A(i,j)),upper(b,TRANS_B(i,j)),upper(c,TRANS_C(i,j)),s);
                for(int i=j;i<dyn_m;i++)
                    FUNCTION(a(TRANS_A(i,j)),lower(b,TRANS_B(i,j)),lower(c,TRANS_C(i,j)),s);
            }
        }
        else // must use diag/lower_
        {
            int dyn_m=a.Rows(),dyn_n=a.Columns();
            assert(dyn_m==dyn_n);
            for(int j=0;j<dyn_n;j++)
            {
                for(int i=0;i<j;i++)
                    FUNCTION(a(TRANS_A(i,j)),upper_(b,TRANS_B(i,j)),upper_(c,TRANS_C(i,j)),s);
                FUNCTION(a(j,j),diag(b,j),diag(c,j),s);
                for(int i=j+1;i<dyn_m;i++)
                    FUNCTION(a(TRANS_A(i,j)),lower_(b,TRANS_B(i,j)),lower_(c,TRANS_C(i,j)),s);
            }
        }
    }
}
#undef FUNCTION
#undef TRANS_A
#undef TRANS_B
#undef TRANS_C
