text:
//Filters each vector in X along dim.
//Filter IIR coefficients are given in vector A with length P+1,
//where P is the IIR filter order (P=0 means only a0).
//Filter FIR coefficients are given in vector B with length Q+1,
//where Q is the FIR filter order (Q=0 means only b0).

//For each vector: a0*Y[t] = (b0*X[t] + b1*X[t-1] + ... + bQ*X[t-Q]) - (a1*Y[t-1] + ... + aP*Y[t-P])

//The calling program must ensure that the sizes are correct, the filter is stable, etc.

#include <stdio.h>
#include <math.h>
#include <cblas.h>

#ifdef __cplusplus
namespace codee {
extern "C" {
#endif

int filter_s (float *Y, const float *X, float *A, float *B, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t P, const size_t Q, const size_t dim);
int filter_d (double *Y, const double *X, double *A, double *B, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t P, const size_t Q, const size_t dim);
int filter_c (float *Y, const float *X, float *A, float *B, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t P, const size_t Q, const size_t dim);
int filter_z (double *Y, const double *X, double *A, double *B, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t P, const size_t Q, const size_t dim);

int filter_s (float *Y, const float *X, float *A, float *B, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t P, const size_t Q, const size_t dim)
{
    if (dim>3u) { fprintf(stderr,"error in filter_s: dim must be in [0 3]\n"); return 1; }

    const size_t N = R*C*S*H;
    const size_t L = (dim==0u) ? R : (dim==1u) ? C : (dim==2u) ? S : H;

    //Deal with a0 (usually 1) and b0; also initialize Y
    float a = *A++, b = *B++;
    if (a==1.0f)
    {
        for (size_t p=P; p>0u; --p, ++A) { *A = -*A; }
        for (size_t n=N; n>0u; --n, ++X, ++Y) { *Y = b * *X; }
    }
    else
    {
        const float b_a = b / a;
        for (size_t p=P; p>0u; --p, ++A) { *A = -*A / a; }
        for (size_t q=Q; q>0u; --q, ++B) { *B /= a; }
        for (size_t n=N; n>0u; --n, ++X, ++Y) { *Y = b_a * *X; }
        B -= Q;
    }
    A -= P; X -= N; Y -= N;

    if (N==0u || L==1u) {}
    else if (L==N)
    {
        //FIR
        if (L<3000u)
        {
            for (size_t q=1u; q<=Q; ++q, ++B, X-=L-q+1u, Y-=L-q)
            {
                b = *B;
                for (size_t l=q; l<L; ++l, ++X) { ++Y; *Y = fmaf(b,*X,*Y); }
            }
        }
        else
        {
            for (size_t q=1u; q<=Q; ++q, ++B) { ++Y; cblas_saxpy((int)(L-q),*B,X,1,Y,1); }
        }
        Y -= Q;

        //IIR
        for (size_t l=0u; l<P; A-=l, ++l, ++Y)
        {
            for (size_t p=1u; p<=l; ++p, ++A) { *Y = fmaf(*A,*(Y-p),*Y); }
        }
        for (size_t l=P; l<L; ++l, A-=P, ++Y)
        {
            for (size_t p=1u; p<=P; ++p, ++A) { *Y = fmaf(*A,*(Y-p),*Y); }
        }
    }
    else
    {
        const size_t K = (iscolmajor) ? ((dim==0u) ? 1u : (dim==1u) ? R : (dim==2u) ? R*C : R*C*S) : ((dim==0u) ? C*S*H : (dim==1u) ? S*H : (dim==2u) ? H : 1u);
        const size_t BS = (iscolmajor && dim==0u) ? C*S*H : K;
        const size_t V = N/L, G = V/BS;

        if (K==1u && (G==1u || BS==1u))
        {
            for (size_t v=V; v>0u; --v)
            {
                //FIR
                if (L<30000u)
                {
                    for (size_t q=1u; q<=Q; ++q, ++B)
                    {
                        b = *B;
                        for (size_t l=q; l<L; ++l, ++X) { ++Y; *Y = fmaf(b,*X,*Y); }
                        if (q<Q) { X -= L-q; Y -= L-q-1u; }
                        else { X += Q; Y -= L-1u; }
                    }
                }
                else
                {
                    for (size_t q=1u; q<=Q; ++q, ++B) { ++Y; cblas_saxpy((int)(L-q),*B,X,1,Y,1); }
                    X += L; Y -= Q;
                }
                B -= Q;

                //IIR
                for (size_t l=0u; l<P; A-=l, ++l, ++Y)
                {
                    for (size_t p=1u; p<=l; ++p, ++A) { *Y = fmaf(*A,*(Y-p),*Y); }
                }
                for (size_t l=P; l<L; ++l, A-=P, ++Y)
                {
                    for (size_t p=1u; p<=P; ++p, ++A) { *Y = fmaf(*A,*(Y-p),*Y); }
                }
            }
        }
        else
        {
            for (size_t g=G; g>0u; --g, X+=BS*(L-1u), Y+=BS*(L-1u))
            {
                for (size_t bs=BS; bs>0u; --bs, ++X, Y-=K*L-1u)
                {
                    //FIR
                    for (size_t q=1u; q<=Q; ++q, ++B) { Y+=K; cblas_saxpy((int)(L-q),*B,X,(int)K,Y,(int)K); }
                    Y -= Q*K; B -= Q;

                    //IIR
                    for (size_t l=0u; l<P; A-=l, ++l, Y+=K)
                    {
                        for (size_t p=1u; p<=l; ++p, ++A) { *Y = fmaf(*A,*(Y-K*p),*Y); }
                    }
                    for (size_t l=P; l<L; ++l, A-=P, Y+=K)
                    {
                        for (size_t p=1u; p<=P; ++p, ++A) { *Y = fmaf(*A,*(Y-K*p),*Y); }
                    }
                }
            }
        }
    }

    return 0;
}

int filter_d (double *Y, const double *X, double *A, double *B, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t P, const size_t Q, const size_t dim)
{
    if (dim>3u) { fprintf(stderr,"error in filter_d: dim must be in [0 3]\n"); return 1; }

    const size_t N = R*C*S*H;
    const size_t L = (dim==0u) ? R : (dim==1u) ? C : (dim==2u) ? S : H;

    //Deal with a0 (usually 1) and b0; also initialize Y
    double a = *A++, b = *B++;
    if (a==1.0)
    {
        for (size_t p=P; p>0u; --p, ++A) { *A = -*A; }
        for (size_t n=N; n>0u; --n, ++X, ++Y) { *Y = b * *X; }
    }
    else
    {
        const double b_a = b / a;
        for (size_t p=P; p>0u; --p, ++A) { *A = -*A / a; }
        for (size_t q=Q; q>0u; --q, ++B) { *B /= a; }
        for (size_t n=N; n>0u; --n, ++X, ++Y) { *Y = b_a * *X; }
        B -= Q;
    }
    A -= P; X -= N; Y -= N;

    if (N==0u || L==1u) {}
    else if (L==N)
    {
        //FIR
        if (L<3000u)
        {
            for (size_t q=1u; q<=Q; ++q, ++B, X-=L-q+1u, Y-=L-q)
            {
                b = *B;
                for (size_t l=q; l<L; ++l, ++X) { ++Y; *Y = fma(b,*X,*Y); }
            }
        }
        else
        {
            for (size_t q=1u; q<=Q; ++q, ++B) { ++Y; cblas_daxpy((int)(L-q),*B,X,1,Y,1); }
        }
        Y -= Q;

        //IIR
        for (size_t l=0u; l<P; A-=l, ++l, ++Y)
        {
            for (size_t p=1u; p<=l; ++p, ++A) { *Y = fma(*A,*(Y-p),*Y); }
        }
        for (size_t l=P; l<L; ++l, A-=P, ++Y)
        {
            for (size_t p=1u; p<=P; ++p, ++A) { *Y = fma(*A,*(Y-p),*Y); }
        }
    }
    else
    {
        const size_t K = (iscolmajor) ? ((dim==0u) ? 1u : (dim==1u) ? R : (dim==2u) ? R*C : R*C*S) : ((dim==0u) ? C*S*H : (dim==1u) ? S*H : (dim==2u) ? H : 1u);
        const size_t BS = (iscolmajor && dim==0u) ? C*S*H : K;
        const size_t V = N/L, G = V/BS;

        if (K==1u && (G==1u || BS==1u))
        {
            for (size_t v=V; v>0u; --v)
            {
                //FIR
                if (L<30000u)
                {
                    for (size_t q=1u; q<=Q; ++q, ++B)
                    {
                        b = *B;
                        for (size_t l=q; l<L; ++l, ++X) { ++Y; *Y = fma(b,*X,*Y); }
                        if (q<Q) { X -= L-q; Y -= L-q-1u; }
                        else { X += Q; Y -= L-1u; }
                    }
                }
                else
                {
                    for (size_t q=1u; q<=Q; ++q, ++B) { ++Y; cblas_daxpy((int)(L-q),*B,X,1,Y,1); }
                    X += L; Y -= Q;
                }
                B -= Q;

                //IIR
                for (size_t l=0u; l<P; A-=l, ++l, ++Y)
                {
                    for (size_t p=1u; p<=l; ++p, ++A) { *Y = fma(*A,*(Y-p),*Y); }
                }
                for (size_t l=P; l<L; ++l, A-=P, ++Y)
                {
                    for (size_t p=1u; p<=P; ++p, ++A) { *Y = fma(*A,*(Y-p),*Y); }
                }
            }
        }
        else
        {
            for (size_t g=G; g>0u; --g, X+=BS*(L-1u), Y+=BS*(L-1u))
            {
                for (size_t bs=BS; bs>0u; --bs, ++X, Y-=K*L-1u)
                {
                    //FIR
                    for (size_t q=1u; q<=Q; ++q, ++B) { Y+=K; cblas_daxpy((int)(L-q),*B,X,(int)K,Y,(int)K); }
                    Y -= Q*K; B -= Q;

                    //IIR
                    for (size_t l=0u; l<P; A-=l, ++l, Y+=K)
                    {
                        for (size_t p=1u; p<=l; ++p, ++A) { *Y = fma(*A,*(Y-K*p),*Y); }
                    }
                    for (size_t l=P; l<L; ++l, A-=P, Y+=K)
                    {
                        for (size_t p=1u; p<=P; ++p, ++A) { *Y = fma(*A,*(Y-K*p),*Y); }
                    }
                }
            }
        }
    }

    return 0;
}

int filter_c (float *Y, const float *X, float *A, float *B, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t P, const size_t Q, const size_t dim)
{
    if (dim>3u) { fprintf(stderr,"error in filter_c: dim must be in [0 3]\n"); return 1; }

    const size_t N = R*C*S*H;
    const size_t L = (dim==0u) ? R : (dim==1u) ? C : (dim==2u) ? S : H;
    float ar, ai, br, bi, xr, xi, yr, yi;

    //Deal with a0 (usually 1) and negate a1 to aP
    //Also, initialize Y
    const float b0r = *B++, b0i = *B++;
    const float a0r = *A++, a0i = *A++, a02 = a0r*a0r + a0i*a0i;
    const float b_ar = (b0r*a0r + b0i*a0i) / a02;
    const float b_ai = (b0i*a0r - b0r*a0i) / a02;
    for (size_t p=P; p>0u; --p)
    {
        ar = -*A; ai = -*(A+1u);
        *A++ = (ar*a0r + ai*a0i) / a02;
        *A++ = (ai*a0r - ar*a0i) / a02;
    }
    for (size_t q=Q; q>0u; --q)
    {
        br = *B++; bi = *B--;
        *B++ = (br*a0r + bi*a0i) / a02;
        *B++ = (bi*a0r - br*a0i) / a02;
    }
    for (size_t n=N; n>0u; --n)
    {
        xr = *X++; xi = *X++;
        *Y++ = b_ar*xr - b_ai*xi;
        *Y++ = b_ar*xi + b_ai*xr;
    }
    A -= 2u*P; B -= 2u*Q; X -= 2u*N; Y -= 2u*N;

    if (N==0u || L==1u) {}
    else if (L==N)
    {
        //FIR
        for (size_t q=1u; q<=Q; ++q, X-=2u*(L-q+1u), Y-=2u*(L-q))
        {
            br = *B++; bi = *B++;
            for (size_t l=q; l<L; ++l)
            {
                Y += 2u;
                xr = *X++; xi = *X++;
                *Y++ += br*xr - bi*xi;
                *Y-- += br*xi + bi*xr;
            }
        }
        Y -= 2u*Q;

        //IIR
        for (size_t l=0u; l<P; ++l, A-=2u*l-2u, Y+=2)
        {
            for (size_t p=1u; p<=l; ++p)
            {
                ar = *A++; ai = *A++;
                yr = *(Y-2u*p); yi = *(Y-2u*p+1u);
                *Y += ar*yr - ai*yi;
                *(Y+1u) += ar*yi + ai*yr;
            }
        }
        for (size_t l=P; l<L; ++l, A-=2u*P, Y+=2)
        {
            for (size_t p=1u; p<=P; ++p)
            {
                ar = *A++; ai = *A++;
                yr = *(Y-2u*p); yi = *(Y-2u*p+1u);
                *Y += ar*yr - ai*yi;
                *(Y+1u) += ar*yi + ai*yr;
            }
        }
    }
    else
    {
        const size_t K = (iscolmajor) ? ((dim==0u) ? 1u : (dim==1u) ? R : (dim==2u) ? R*C : R*C*S) : ((dim==0u) ? C*S*H : (dim==1u) ? S*H : (dim==2u) ? H : 1u);
        const size_t BS = (iscolmajor && dim==0u) ? C*S*H : K;
        const size_t V = N/L, G = V/BS;

        if (K==1u && (G==1u || BS==1u))
        {
            for (size_t v=V; v>0u; --v)
            {
                //FIR
                for (size_t q=1u; q<=Q; ++q)
                {
                    br = *B++; bi = *B++;
                    for (size_t l=q; l<L; ++l)
                    {
                        Y += 2u;
                        xr = *X++; xi = *X++;
                        *Y++ += br*xr - bi*xi;
                        *Y-- += br*xi + bi*xr;
                    }
                    if (q<Q) { X -= 2u*(L-q); Y -= 2u*(L-q-1u); }
                    else { X += 2u*Q; Y -= 2u*L - 2u; }
                }
                B -= 2u*Q;

                //IIR
                for (size_t l=0u; l<P; ++l, A-=2u*l-2u, Y+=2)
                {
                    for (size_t p=1u; p<=l; ++p)
                    {
                        ar = *A++; ai = *A++;
                        yr = *(Y-2u*p); yi = *(Y-2u*p+1u);
                        *Y += ar*yr - ai*yi;
                        *(Y+1u) += ar*yi + ai*yr;
                    }
                }
                for (size_t l=P; l<L; ++l, A-=2u*P, Y+=2)
                {
                    for (size_t p=1u; p<=P; ++p)
                    {
                        ar = *A++; ai = *A++;
                        yr = *(Y-2u*p); yi = *(Y-2u*p+1u);
                        *Y += ar*yr - ai*yi;
                        *(Y+1u) += ar*yi + ai*yr;
                    }
                }
            }
        }
        else
        {
            for (size_t g=G; g>0u; --g, X+=2u*BS*(L-1u), Y+=2u*BS*(L-1u))
            {
                for (size_t bs=BS; bs>0u; --bs, X-=2u*K*L-2u, Y-=2u*K*L-2u)
                {
                    //FIR
                    //for (size_t q=1u; q<=Q; X-=2u*K*(L-q), ++q)
                    for (size_t q=1u; q<=Q; ++q)
                    {
                        br = *B++; bi = *B++;
                        for (size_t l=q; l<L; ++l, X+=2u*K)
                        {
                            Y += 2u*K;
                            xr = *X; xi = *(X+1);
                            *Y++ += br*xr - bi*xi;
                            *Y-- += br*xi + bi*xr;
                        }
                        if (q<Q) { X -= 2u*K*(L-q); Y -= 2u*K*(L-q-1u); }
                        else { X += 2u*K*Q; Y -= 2u*K*(L-1u); }
                    }
                    B -= 2u*Q;

                    //IIR
                    for (size_t l=0u; l<P; ++l, A-=2u*l-2u, Y+=2u*K)
                    {
                        for (size_t p=1u; p<=l; ++p)
                        {
                            ar = *A++; ai = *A++;
                            yr = *(Y-2u*K*p); yi = *(Y-2u*K*p+1u);
                            *Y += ar*yr - ai*yi;
                            *(Y+1u) += ar*yi + ai*yr;
                        }
                    }
                    for (size_t l=P; l<L; ++l, A-=2u*P, Y+=2u*K)
                    {
                        for (size_t p=1u; p<=P; ++p)
                        {
                            ar = *A++; ai = *A++;
                            yr = *(Y-2u*K*p); yi = *(Y-2u*K*p+1u);
                            *Y += ar*yr - ai*yi;
                            *(Y+1u) += ar*yi + ai*yr;
                        }
                    }
                }
            }
        }
    }

    return 0;
}

int filter_z (double *Y, const double *X, double *A, double *B, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t P, const size_t Q, const size_t dim)
{
    if (dim>3u) { fprintf(stderr,"error in filter_z: dim must be in [0 3]\n"); return 1; }

    const size_t N = R*C*S*H;
    const size_t L = (dim==0u) ? R : (dim==1u) ? C : (dim==2u) ? S : H;
    double ar, ai, br, bi, xr, xi, yr, yi;

    //Deal with a0 (usually 1) and negate a1 to aP
    //Also, initialize Y
    const double b0r = *B++, b0i = *B++;
    const double a0r = *A++, a0i = *A++, a02 = a0r*a0r + a0i*a0i;
    const double b_ar = (b0r*a0r + b0i*a0i) / a02;
    const double b_ai = (b0i*a0r - b0r*a0i) / a02;
    for (size_t p=P; p>0u; --p)
    {
        ar = -*A; ai = -*(A+1u);
        *A++ = (ar*a0r + ai*a0i) / a02;
        *A++ = (ai*a0r - ar*a0i) / a02;
    }
    for (size_t q=Q; q>0u; --q)
    {
        br = *B++; bi = *B--;
        *B++ = (br*a0r + bi*a0i) / a02;
        *B++ = (bi*a0r - br*a0i) / a02;
    }
    for (size_t n=N; n>0u; --n)
    {
        xr = *X++; xi = *X++;
        *Y++ = b_ar*xr - b_ai*xi;
        *Y++ = b_ar*xi + b_ai*xr;
    }
    A -= 2u*P; B -= 2u*Q; X -= 2u*N; Y -= 2u*N;

    if (N==0u || L==1u) {}
    else if (L==N)
    {
        //FIR
        for (size_t q=1u; q<=Q; ++q, X-=2u*(L-q+1u), Y-=2u*(L-q))
        {
            br = *B++; bi = *B++;
            for (size_t l=q; l<L; ++l)
            {
                Y += 2u;
                xr = *X++; xi = *X++;
                *Y++ += br*xr - bi*xi;
                *Y-- += br*xi + bi*xr;
            }
        }
        Y -= 2u*Q;

        //IIR
        for (size_t l=0u; l<P; ++l, A-=2u*l-2u, Y+=2)
        {
            for (size_t p=1u; p<=l; ++p)
            {
                ar = *A++; ai = *A++;
                yr = *(Y-2u*p); yi = *(Y-2u*p+1u);
                *Y += ar*yr - ai*yi;
                *(Y+1u) += ar*yi + ai*yr;
            }
        }
        for (size_t l=P; l<L; ++l, A-=2u*P, Y+=2)
        {
            for (size_t p=1u; p<=P; ++p)
            {
                ar = *A++; ai = *A++;
                yr = *(Y-2u*p); yi = *(Y-2u*p+1u);
                *Y += ar*yr - ai*yi;
                *(Y+1u) += ar*yi + ai*yr;
            }
        }
    }
    else
    {
        const size_t K = (iscolmajor) ? ((dim==0u) ? 1u : (dim==1u) ? R : (dim==2u) ? R*C : R*C*S) : ((dim==0u) ? C*S*H : (dim==1u) ? S*H : (dim==2u) ? H : 1u);
        const size_t BS = (iscolmajor && dim==0u) ? C*S*H : K;
        const size_t V = N/L, G = V/BS;

        if (K==1u && (G==1u || BS==1u))
        {
            for (size_t v=V; v>0u; --v)
            {
                //FIR
                for (size_t q=1u; q<=Q; ++q)
                {
                    br = *B++; bi = *B++;
                    for (size_t l=q; l<L; ++l)
                    {
                        Y += 2u;
                        xr = *X++; xi = *X++;
                        *Y++ += br*xr - bi*xi;
                        *Y-- += br*xi + bi*xr;
                    }
                    if (q<Q) { X -= 2u*(L-q); Y -= 2u*(L-q-1u); }
                    else { X += 2u*Q; Y -= 2u*L - 2u; }
                }
                B -= 2u*Q;

                //IIR
                for (size_t l=0u; l<P; ++l, A-=2u*l-2u, Y+=2)
                {
                    for (size_t p=1u; p<=l; ++p)
                    {
                        ar = *A++; ai = *A++;
                        yr = *(Y-2u*p); yi = *(Y-2u*p+1u);
                        *Y += ar*yr - ai*yi;
                        *(Y+1u) += ar*yi + ai*yr;
                    }
                }
                for (size_t l=P; l<L; ++l, A-=2u*P, Y+=2)
                {
                    for (size_t p=1u; p<=P; ++p)
                    {
                        ar = *A++; ai = *A++;
                        yr = *(Y-2u*p); yi = *(Y-2u*p+1u);
                        *Y += ar*yr - ai*yi;
                        *(Y+1u) += ar*yi + ai*yr;
                    }
                }
            }
        }
        else
        {
            for (size_t g=G; g>0u; --g, X+=2u*BS*(L-1u), Y+=2u*BS*(L-1u))
            {
                for (size_t bs=BS; bs>0u; --bs, X-=2u*K*L-2u, Y-=2u*K*L-2u)
                {
                    //FIR
                    for (size_t q=1u; q<=Q; ++q)
                    {
                        br = *B++; bi = *B++;
                        for (size_t l=q; l<L; ++l, X+=2u*K)
                        {
                            Y += 2u*K;
                            xr = *X; xi = *(X+1);
                            *Y++ += br*xr - bi*xi;
                            *Y-- += br*xi + bi*xr;
                        }
                        if (q<Q) { X -= 2u*K*(L-q); Y -= 2u*K*(L-q-1u); }
                        else { X += 2u*K*Q; Y -= 2u*K*(L-1u); }
                    }
                    B -= 2u*Q;

                    //IIR
                    for (size_t l=0u; l<P; ++l, A-=2u*l-2u, Y+=2u*K)
                    {
                        for (size_t p=1u; p<=l; ++p)
                        {
                            ar = *A++; ai = *A++;
                            yr = *(Y-2u*K*p); yi = *(Y-2u*K*p+1u);
                            *Y += ar*yr - ai*yi;
                            *(Y+1u) += ar*yi + ai*yr;
                        }
                    }
                    for (size_t l=P; l<L; ++l, A-=2u*P, Y+=2u*K)
                    {
                        for (size_t p=1u; p<=P; ++p)
                        {
                            ar = *A++; ai = *A++;
                            yr = *(Y-2u*K*p); yi = *(Y-2u*K*p+1u);
                            *Y += ar*yr - ai*yi;
                            *(Y+1u) += ar*yi + ai*yr;
                        }
                    }
                }
            }
        }
    }

    return 0;
}

#ifdef __cplusplus
}
}
#endif


text:
#pragma once

#include <iostream>

#include <gsl/gsl.h>

namespace util{
template<typename T>
void print(T const&x){
    std::cerr <<x << " ";
}

template<typename T, int64_t M>
void print(gsl::span<T,M> const span){
    for(auto x : span) print(x);
    std::cerr << " :1D"<<std::endl;
}

template<typename T, int64_t M,int64_t N>
void print(gsl::span<T,M,N> const span){
    ////ranged for iterates flattened span.
    //for(auto x : span) print(x);
    for(int64_t i=0; i<M; ++i) print(span[i]);
    std::cerr << " :2D static"<<std::endl;
}

template<typename T, int64_t M,int64_t N>
void print(gsl::span<const T,gsl::dynamic_range,M,N> span){
    for(int64_t i=0; i<M; ++i) print(span[i]);
    //for(auto it=span.cbegin(); it<span.cend(); ++it) print_span(*it);
    std::cerr << " :2D dynamic_range"<<std::endl;
}
// template<typename T>
// void print(std::vector<T> const & vec){
//     print(gsl::as_span(vec));
// }
// template<typename T, size_t M>
// void print(std::array<T,M> arr){
//     print(gsl::as_span(arr));
// }

}//namespcae util

ext: h

max_stars_repo_name: uphere-co/nlp-prototype

max_issues_repo_name: uphere-co/nlp-prototype

max_forks_repo_name: uphere-co/nlp-prototype


text:
#pragma once

#include <gsl/span>
#include <atomic>

namespace Halley {
	// This is a lock-free, wait-free ring buffer
	// It is only thread safe for one producer and one consumer at any given time
	// Loosely based on implementation from Game Audio Programming 2

    template <typename T>
    class RingBuffer {
    public:
    	explicit RingBuffer(size_t capacity)
	        : numEntries(0)
    	{
            entries.resize(capacity);
    	}

    	size_t availableToRead() const
    	{
            return numEntries.load();
    	}

        size_t availableToWrite() const
        {
            return entries.size() - numEntries.load();
        }

    	bool canWrite(size_t n) const
    	{
            return numEntries.load() + n <= entries.size();
    	}

    	bool canRead(size_t n) const
    	{
            return numEntries.load() >= n;
    	}

    	bool empty() const
    	{
            return numEntries.load() == 0;
    	}

        void writeOne(T e)
        {
            Expects(canWrite(1));
            entries[writePos] = std::move(e);
            writePos = (writePos + 1) % entries.size();
            ++numEntries;
        }

        void write(gsl::span<T> es)
        {
            const size_t numToWrite = size_t(es.size());
            Expects(canWrite(numToWrite));
            const size_t spaceToEnd = entries.size() - writePos;
            const size_t nToWrite1 = std::min(spaceToEnd, numToWrite);

            for (size_t i = 0; i < nToWrite1; ++i) {
                entries[writePos + i] = std::move(es[i]);
            }

            const size_t nToWrite2 = numToWrite - nToWrite1;
            for (size_t i = 0; i < nToWrite2; ++i) {
                entries[i] = std::move(es[i + nToWrite1]);
            }

            writePos = (writePos + numToWrite) % entries.size();
            numEntries.fetch_add(numToWrite);
        }

    	void write(gsl::span<const T> es)
    	{
            const size_t numToWrite = size_t(es.size());
            Expects(canWrite(numToWrite));
            const size_t spaceToEnd = entries.size() - writePos;
            const size_t nToWrite1 = std::min(spaceToEnd, numToWrite);

    		for (size_t i = 0; i < nToWrite1; ++i) {
                entries[writePos + i] = es[i];
    		}

            const size_t nToWrite2 = numToWrite - nToWrite1;
    		for (size_t i = 0; i < nToWrite2; ++i) {
                entries[i] = es[i + nToWrite1];
    		}

            writePos = (writePos + numToWrite) % entries.size();
            numEntries.fetch_add(numToWrite);
    	}

    	T readOne()
    	{
            Expects(canRead(1));
            T v = std::move(entries[readPos]);
    		entries[readPos] = T();
            readPos = (readPos + 1) % entries.size();
            --numEntries;
            return v;
    	}

    	void read(gsl::span<T> es)
    	{
            const size_t numToRead = size_t(es.size());
            Expects(canRead(numToRead));
            const size_t spaceToEnd = entries.size() - readPos;
            const size_t nToRead1 = std::min(spaceToEnd, numToRead);

            for (size_t i = 0; i < nToRead1; ++i) {
                es[i] = std::move(entries[readPos + i]);
            	entries[readPos + i] = T();
            }

            const size_t nToRead2 = numToRead - nToRead1;
            for (size_t i = 0; i < nToRead2; ++i) {
                es[i + nToRead1] = std::move(entries[i]);
            	entries[i] = T();
            }

            readPos = (readPos + numToRead) % entries.size();
            numEntries.fetch_sub(numToRead);
    	}

    private:
        size_t readPos = 0;
        size_t writePos = 0;
        std::atomic<size_t> numEntries;
        Vector<T> entries;
    };
}

max_stars_repo_name: amrezzd/halley

max_issues_repo_name: amrezzd/halley

max_forks_repo_name: amrezzd/halley


text:
/*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*
**                                                                                **
** This file forms part of the Underworld geophysics modelling application.         **
**                                                                                **
** For full license and copyright information, please refer to the LICENSE.md file  **
** located at the project root, or contact the authors.                             **
**                                                                                **
**~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*/
#include <petsc.h>
#include <petscmat.h>
#include <petscvec.h>
#include <petscksp.h>
#include <petscpc.h>
#include <petscversion.h>
#if ( (PETSC_VERSION_MAJOR >= 3) && (PETSC_VERSION_MINOR >=3) )
  #if (PETSC_VERSION_MINOR >=6)
     #include <petsc/private/kspimpl.h>
  #else
     #include <petsc-private/kspimpl.h>
  #endif
#else
  #include <private/kspimpl.h>
#endif
#include <StGermain/StGermain.h>
#include <StgDomain/StgDomain.h>
#include <StgFEM/StgFEM.h>
#include <PICellerator/PICellerator.h>
#include <Underworld/Underworld.h>
#include "Solvers/SLE/SLE.h" /* to give the AugLagStokes_SLE type */
#include "Solvers/KSPSolvers/KSPSolvers.h"
#include "BSSCR.h"
#include "stokes_block_scaling.h"
#include "stokes_mvblock_scaling.h"
#include "writeMatVec.h"
#include "createK2.h"
#include "mg.h"
#include "summary.h"
#include "ksp_pressure_nullspace.h"

#define BSSCR_GetPetscMatrix( matrix ) ( (Mat)(matrix) )
#define BSSCR_GetPetscVector( vector ) ( (Vec)(vector) )

typedef struct {
    PetscTruth initialrtol;    /* default relative residual decrease is computing from initial residual, not rhs */
    PetscTruth mininitialrtol; /* default relative residual decrease is computing from min of initial residual and rhs */
    Vec        work;
    PetscReal  pointwise_max;
} KSPPWConvergedCtx;

PetscErrorCode BSSCR_KSPNormInfMonitor( KSP ksp, PetscInt iteration, PetscReal residualNorm, void *dummy);
PetscErrorCode BSSCR_KSPNormInfToNorm2Monitor( KSP ksp, PetscInt iteration, PetscReal residualNorm, void *dummy);

PetscErrorCode BSSCR_KSPPWConvergedCreate(void **ctx);
PetscErrorCode BSSCR_KSPPWConvergedDestroy(void *ctx);
PetscErrorCode BSSCR_KSPNormInfConverged(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *ctx);
PetscErrorCode BSSCR_KSPNormInfToNorm2Monitor(KSP ksp,PetscInt n,PetscReal rnorm, void *dummy);
PetscErrorCode BSSCR_KSPNormInfMonitor(KSP ksp,PetscInt n,PetscReal rnorm, void *dummy);
PetscErrorCode BSSCR_KSPNorm2RawMonitor(KSP ksp,PetscInt n,PetscReal rnorm, void *dummy);
PetscErrorCode BSSCR_KSPSetNormInfConvergenceTest(KSP ksp);

#undef __FUNCT__
#define __FUNCT__ "BSSCR_DRIVER_auglag"
PetscErrorCode BSSCR_DRIVER_auglag( KSP ksp, Mat stokes_A, Vec stokes_x, Vec stokes_b, Mat approxS,
                                          MatStokesBlockScaling BA, PetscTruth sym, KSP_BSSCR * bsscrp_self )
{
    PetscTruth uzawastyle, KisJustK=PETSC_TRUE, restorek, change_A11rhspresolve;
    PetscTruth usePreviousGuess, useNormInfStoppingConditions, useNormInfMonitor, found, forcecorrection;
    PetscTruth change_backsolve, mg_active, get_flops;
    PetscErrorCode ierr;
    KSPConvergedReason reason;

    PetscInt max_it,min_it;
    KSP ksp_inner, ksp_S, ksp_new_inner;
    PC pc_S, pcInner;
    Mat K,G,D,C, S, K2;
    Vec u,p,f,f2=0,f3=0,h, h_hat,t;
    Vec f_tmp;

    MGContext mgCtx;
    double mgSetupTime, problemBuildTime, scrSolveTime, RHSSolveTime, a11SingleSolveTime, penaltyNumber;// hFactor;
    double backsolveSetupTime, scrSetupTime, RHSSetupTime;
    int been_here = bsscrp_self->been_here;

    char name[PETSC_MAX_PATH_LEN];
    char suffix[PETSC_MAX_PATH_LEN];
    char str[PETSC_MAX_PATH_LEN];
    PetscTruth flg, extractMats;
    PetscLogDouble flopsA,flopsB;

    /***************************************************************************************************************/
    /***************************************************************************************************************/
    //if( bsscrp_self->solver->st_sle->context->loadFromCheckPoint ){
    //    been_here=1;
    //}
    /* get sub matrix / vector objects */
    /* note that here, the matrix D should always exist. It is set up in  _StokesBlockKSPInterface_Solve in StokesBlockKSPInterface.c */
    /* now extract K,G etc from a MatNest object */
    MatNestGetSubMat( stokes_A, 0,0, &K );
    MatNestGetSubMat( stokes_A, 0,1, &G );
    MatNestGetSubMat( stokes_A, 1,0, &D );if(!D){ PetscPrintf( PETSC_COMM_WORLD, "D does not exist but should!!\n"); exit(1); }
    MatNestGetSubMat( stokes_A, 1,1, &C );
    VecNestGetSubVec( stokes_x, 0, &u );
    VecNestGetSubVec( stokes_x, 1, &p );
    VecNestGetSubVec( stokes_b, 0, &f );
    VecNestGetSubVec( stokes_b, 1, &h );

    // Try this ...
    // MatMPIAIJSetPreallocation(K, 375 ,PETSC_NULL, 375, PETSC_NULL);

    PetscPrintf( PETSC_COMM_WORLD, "AUGMENTED LAGRANGIAN K2 METHOD " );
    PetscPrintf( PETSC_COMM_WORLD, "- Penalty = %f\n\n", bsscrp_self->solver->penaltyNumber );
    sprintf(suffix,"%s","x");

    found = PETSC_FALSE;
    get_flops = PETSC_FALSE;
    PetscOptionsGetTruth( PETSC_NULL, "-get_flops", &get_flops, &found);

    PetscOptionsGetString( PETSC_NULL, "-matsuffix", suffix, PETSC_MAX_PATH_LEN-1, &extractMats );

    flg=0;
    PetscOptionsGetString( PETSC_NULL, "-matdumpdir", name, PETSC_MAX_PATH_LEN-1, &flg );
    if(flg){
        sprintf(str,"%s/K%s",name, suffix);
        bsscr_writeMat( K, str, "Writing K matrix in al Solver");
        sprintf(str,"%s/G%s",name, suffix);
        bsscr_writeMat( G, str, "Writing G matrix in al Solver");
        sprintf(str,"%s/D%s",name, suffix);
        bsscr_writeMat( D, str, "Writing D matrix in al Solver");
        sprintf(str,"%s/f%s",name, suffix);
        bsscr_writeVec( f, str, "Writing f vector in al Solver");
        sprintf(str,"%s/h%s",name, suffix);
        bsscr_writeVec( h, str, "Writing h vector in al Solver");
        sprintf(str,"%s/Shat%s",name,suffix);
        bsscr_writeMat( approxS, str, "Writing Shat matrix in al Solver");
        if(C){
          sprintf(str,"%s/C%s",name,suffix);
          bsscr_writeMat( C, str, "Writing C matrix in al Solver");
        }
    }

    mg_active=PETSC_TRUE;
    PetscOptionsGetTruth( PETSC_NULL ,"-A11_mg_active", &mg_active, &found );
    bsscrp_self->solver->mg_active=mg_active;

    penaltyNumber = bsscrp_self->solver->penaltyNumber;

    /***************************************************************************************************************/
    /***************************************************************************************************************/
    /******  GET K2   ****************************************************************************************/

    if(penaltyNumber > 1e-10 && bsscrp_self->k2type){
        flg=0;
        PetscOptionsGetString( PETSC_NULL, "-matdumpdir", name, PETSC_MAX_PATH_LEN-1, &flg );
        if(flg){
            sprintf(str,"%s/K2%s",name, suffix);
            bsscr_writeMat( bsscrp_self->K2, str, "Writing K2 matrix in al Solver");
        }

        K2=bsscrp_self->K2;

        scrSolveTime = MPI_Wtime();
        ierr=MatAXPY(K,penaltyNumber,K2,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);/* Computes K = penaltyNumber*K2 + K */
        scrSolveTime =  MPI_Wtime() - scrSolveTime;

        PetscPrintf( PETSC_COMM_WORLD, "\n\t* K+p*K2 in time: %lf seconds\n\n", scrSolveTime);
        KisJustK=PETSC_FALSE;
        forcecorrection=PETSC_TRUE;
        PetscOptionsGetTruth( PETSC_NULL ,"-force_correction", &forcecorrection, &found );
        if(forcecorrection){
            if(bsscrp_self->f2 && forcecorrection){
                f2=bsscrp_self->f2;
                ierr=VecAXPY(f,penaltyNumber,f2);/*  f <- f +a*f2 */
            }
            else {
                switch (bsscrp_self->k2type) {
                case (K2_GG):
                {
                    VecDuplicate( u, &f3 );
                    MatMult( G, h, f3);   ierr=VecAXPY(f,penaltyNumber,f3);/*  f <- f +a*f2 */
                }
                break;
                case (K2_GMG):
                case (K2_DGMGD):
                {   /* THIS IS THE ONE THAT WE KNOW WORKS PROPERLY */
                    Mat M;
                    Vec Mdiag;
                    VecDuplicate( u, &f3 );
                    M = bsscrp_self->solver->mStiffMat->matrix;
                    MatGetVecs( M, &Mdiag, PETSC_NULL );
                    MatGetDiagonal( M, Mdiag );
                    VecReciprocal(Mdiag);
                    VecPointwiseMult(Mdiag,Mdiag,h);
                    MatMult( G, Mdiag, f3);   ierr=VecAXPY(f,penaltyNumber,f3);/*  f <- f +a*f2 */
                    Stg_VecDestroy(&Mdiag);
                }
                break;
                case (K2_NULL):
                case (K2_SLE):
                {
                    ;
                }
                break;
                }
            }
        }
    }

    /* Create Schur complement matrix */
    MatCreateSchurComplement(K,K,G,D,C, &S);
    MatSchurComplementGetKSP( S, &ksp_inner);
    KSPGetPC( ksp_inner, &pcInner );

    /***************************************************************************************************************/
    /***************************************************************************************************************/
    /*********    SET PREFIX FOR INNER/VELOCITY KSP    *************************************************************/

    KSPSetOptionsPrefix( ksp_inner, "A11_" );
    KSPSetFromOptions( ksp_inner );
    Stg_KSPSetOperators(ksp_inner, K, K, DIFFERENT_NONZERO_PATTERN);

    useNormInfStoppingConditions = PETSC_FALSE;
    PetscOptionsGetTruth( PETSC_NULL ,"-A11_use_norm_inf_stopping_condition", &useNormInfStoppingConditions, &found );
    if(useNormInfStoppingConditions)
        BSSCR_KSPSetNormInfConvergenceTest( ksp_inner );

    useNormInfMonitor = PETSC_FALSE;
    PetscOptionsGetTruth( PETSC_NULL, "-A11_ksp_norm_inf_monitor", &useNormInfMonitor, &found );
    if(useNormInfMonitor)  KSPMonitorSet( ksp_inner, BSSCR_KSPNormInfMonitor, PETSC_NULL, PETSC_NULL );

    useNormInfMonitor = PETSC_FALSE;
    PetscOptionsGetTruth( PETSC_NULL, "-A11_ksp_norm_inf_to_norm_2_monitor", &useNormInfMonitor, &found );
    if(useNormInfMonitor)  KSPMonitorSet( ksp_inner, BSSCR_KSPNormInfToNorm2Monitor, PETSC_NULL, PETSC_NULL );

    usePreviousGuess = PETSC_FALSE;
//    if(been_here)
//        PetscOptionsGetTruth( PETSC_NULL, "-scr_use_previous_guess", &usePreviousGuess, &found );
    /***************************************************************************************************************/

    /***************************************************************************************************************/
    /***************************************************************************************************************/
    /* If multigrid is enabled, set it now. */
    change_A11rhspresolve = PETSC_FALSE;
    PetscOptionsGetTruth( PETSC_NULL, "-change_A11rhspresolve", &change_A11rhspresolve, &found );

    if(bsscrp_self->solver->mg_active && !change_A11rhspresolve) { mgSetupTime=setupMG( bsscrp_self, ksp_inner, pcInner, K, &mgCtx ); }
    /***************************************************************************************************************/
    /***************************************************************************************************************/
    /* create right hand side */

    if(change_A11rhspresolve){
      KSPCreate(PETSC_COMM_WORLD, &ksp_new_inner);
      Stg_KSPSetOperators(ksp_new_inner, K, K, DIFFERENT_NONZERO_PATTERN);
      KSPSetOptionsPrefix(ksp_new_inner, "rhsA11_");
      MatSchurComplementSetKSP( S, ksp_new_inner );/* this call destroys the ksp_inner that is already set on S */
      ksp_inner=ksp_new_inner;
      KSPGetPC( ksp_inner, &pcInner );
      KSPSetFromOptions(ksp_inner); /* make sure we are setting up our solver how we want it */
    }

    MatGetVecs( S, PETSC_NULL, &h_hat );
    /* It may be the case that the current velocity solution might not be bad guess for f_tmp? <-- maybe not */
    MatGetVecs( K, PETSC_NULL, &f_tmp );

    /*************************************/
    /*************************************/

//    if (!been_here) {
//        VecDuplicate( u, &uStar );
//        VecSet(uStar,0.0);
//    }
//
//    if(usePreviousGuess) {
//        VecCopy( uStar, f_tmp);
//        KSPSetInitialGuessNonzero( ksp_inner, PETSC_TRUE );
//    }
//    else {
        KSPSetInitialGuessNonzero( ksp_inner, PETSC_FALSE );
//    }

    RHSSetupTime = MPI_Wtime();
    KSPSetUp(ksp_inner);
    RHSSetupTime = MPI_Wtime() - RHSSetupTime;
    bsscrp_self->solver->stats.velocity_presolve_setup_time = RHSSetupTime;

    RHSSolveTime = MPI_Wtime();
    KSPSolve(ksp_inner, f, f_tmp);
    KSPGetConvergedReason( ksp_inner, &reason ); {if (reason < 0) bsscrp_self->solver->fhat_reason=(int)reason; }
    KSPGetIterationNumber( ksp_inner, &bsscrp_self->solver->stats.velocity_presolve_its );
    RHSSolveTime =  MPI_Wtime() - RHSSolveTime;
    bsscrp_self->solver->stats.velocity_presolve_time = RHSSolveTime;

    // VecCopy( f_tmp, uStar);
    // scrSolveTime =  RHSsolveTime;

    MatMult(D, f_tmp, h_hat);
    VecAYPX(h_hat, -1.0, h); /* Computes y = x + alpha y.  h_hat -> h - Gt*K^(-1)*f*/
    Stg_VecDestroy(&f_tmp);

    /*************************************/
    /*************************************/

    if(bsscrp_self->solver->mg_active && change_A11rhspresolve) {
      //Stg_KSPDestroy(&ksp_inner );
      KSPCreate(PETSC_COMM_WORLD, &ksp_new_inner);
      Stg_KSPSetOperators(ksp_new_inner, K, K, DIFFERENT_NONZERO_PATTERN);
      KSPSetOptionsPrefix( ksp_new_inner, "A11_" );
      MatSchurComplementSetKSP( S, ksp_new_inner );
      ksp_inner=ksp_new_inner;
      //MatSchurSetKSP( S, ksp_inner );
      KSPGetPC( ksp_inner, &pcInner );
      KSPSetFromOptions( ksp_inner );
      mgSetupTime=setupMG( bsscrp_self, ksp_inner, pcInner, K, &mgCtx );
    }

    /* create solver for S p = h_hat */
    KSPCreate( PETSC_COMM_WORLD, &ksp_S );
    KSPSetOptionsPrefix( ksp_S, "scr_");

    /* By default use the UW approxS Schur preconditioner -- same as the one used by the Uzawa solver */
    /* Note that if scaling is activated then the approxS matrix has been scaled already */
    /* so no need to rebuild in the case of scaling as we have been doing */
    if(!approxS){
        PetscPrintf( PETSC_COMM_WORLD,  "WARNING approxS is NULL\n");
    }

    Stg_KSPSetOperators( ksp_S, S, S, SAME_NONZERO_PATTERN );
    KSPSetType( ksp_S, "cg" );
    KSPGetPC( ksp_S, &pc_S );
    BSSCR_BSSCR_StokesCreatePCSchur2( K,G,D,C,approxS, pc_S, sym, bsscrp_self );

    flg=0;
    PetscOptionsGetString( PETSC_NULL, "-NN", name, PETSC_MAX_PATH_LEN-1, &flg );
    if(flg){
      Mat Smat, Pmat;
      //MatStructure mstruct;
      Stg_PCGetOperators( pc_S, &Smat, &Pmat, NULL );
      sprintf(str,"%s/Pmat%s",name, suffix);
      bsscr_writeMat( Pmat, str, "Writing Pmat matrix in al Solver");
    }

    uzawastyle=PETSC_FALSE;
    PetscOptionsGetTruth( PETSC_NULL, "-uzawa_style", &uzawastyle, &found );
    if(uzawastyle){
        /* now want to set up the ksp_S->pc to be of type ksp (gmres) by default to match Uzawa */
        KSP pc_ksp;
        KSPGetPC( ksp_S, &pc_S );
        PCSetType(pc_S,PCKSP);
        PCKSPGetKSP( pc_S, &pc_ksp);
        KSPSetType(pc_ksp, "gmres" );
        KSPSetOptionsPrefix( pc_ksp, "scrPCKSP_");
        KSPSetFromOptions( pc_ksp );
    }

    KSPSetFromOptions( ksp_S );
    /* Set specific monitor test */
    KSPGetTolerances( ksp_S, PETSC_NULL, PETSC_NULL, PETSC_NULL, &max_it );

    if(usePreviousGuess) {   /* Note this should actually look at checkpoint information */
        KSPSetInitialGuessNonzero( ksp_S, PETSC_TRUE ); }
    else {
        KSPSetInitialGuessNonzero( ksp_S, PETSC_FALSE ); }

    /***************************************************************************************************************/
    /* Pressure / Velocity Solve */
    /***************************************************************************************************************/

    /***************************************************************************************************************/
    /*******     SET CONVERGENCE TESTS     *************************************************************************/

    useNormInfStoppingConditions = PETSC_FALSE;
    PetscOptionsGetTruth( PETSC_NULL ,"-scr_use_norm_inf_stopping_condition", &useNormInfStoppingConditions, &found );
    if(useNormInfStoppingConditions)
        BSSCR_KSPSetNormInfConvergenceTest(ksp_S);

    useNormInfMonitor = PETSC_FALSE;
    PetscOptionsGetTruth( PETSC_NULL, "-scr_ksp_norm_inf_monitor", &useNormInfMonitor, &found );
    if(useNormInfMonitor)
        KSPMonitorSet( ksp_S, BSSCR_KSPNormInfToNorm2Monitor, PETSC_NULL, PETSC_NULL );

    /***************************************************************************************************************/
    /***************************************************************************************************************/
    /*******   PRESSURE SOLVE   ************************************************************************************/

    //PetscPrintf( PETSC_COMM_WORLD, "\t* KSPSolve( ksp_S, h_hat, p )\n");
    /* if h_hat needs to be fixed up ..take out any nullspace vectors here */
    /* we want to check that there is no "noise" in the null-space in the h vector */
    /* this causes problems when we are trying to solve a Jacobian system when the Residual is almost converged */

    if(bsscrp_self->check_pressureNS){
        bsscrp_self->buildPNS(ksp);/* build and set nullspace vectors on bsscr - which is on ksp (function pointer is set in KSPSetUp_BSSCR) */
    }

    PetscScalar hnorm, gnorm;
    MatNorm(G,NORM_INFINITY,&gnorm);
    VecNorm(h_hat, NORM_2, &hnorm);
    hnorm=hnorm/gnorm;

    if((hnorm < 1e-6) && (hnorm > 1e-20)){
        VecScale(h_hat,1.0/hnorm);
    }
    /* test to see if v or t are in nullspace of G and orthogonalize wrt h_hat if needed */
    KSPRemovePressureNullspace_BSSCR(ksp, h_hat);
    /* set convergence test to use min_it */
    found = PETSC_FALSE;
    min_it = 0;
    PetscOptionsGetInt( PETSC_NULL,"-scr_ksp_set_min_it_converge", &min_it, &found);
    if(found && min_it > 0){
        BSSCR_KSPSetConvergenceMinIts(ksp_S, min_it, bsscrp_self);
    }

    /** Pressure Setup **/
    scrSetupTime = MPI_Wtime();
    KSPSetUp(ksp_S);
    scrSetupTime = MPI_Wtime() - scrSetupTime;
    bsscrp_self->solver->stats.velocity_pressuresolve_setup_time = scrSetupTime;

    /** Pressure Solve **/
    if(get_flops) PetscGetFlops(&flopsA);
    scrSolveTime = MPI_Wtime();
    KSPSolve( ksp_S, h_hat, p );
    scrSolveTime =  MPI_Wtime() - scrSolveTime;

    KSPGetConvergedReason( ksp_S, &reason ); {if (reason < 0) bsscrp_self->solver->outer_reason=(int)reason; }
    /*************************************/
    /*************************************/
#if ( (PETSC_VERSION_MAJOR >= 3) && (PETSC_VERSION_MINOR >= 6 ) && (PETSC_VERSION_SUBMINOR >= 1 ))
    KSPGetTotalIterations(ksp_inner, &bsscrp_self->solver->stats.velocity_pressuresolve_its);
    if(!change_A11rhspresolve){
      bsscrp_self->solver->stats.velocity_pressuresolve_its -= bsscrp_self->solver->stats.velocity_presolve_its;
    }
#else
    bsscrp_self->solver->stats.velocity_pressuresolve_its=-1;
#endif

    if(get_flops) {
      PetscGetFlops(&flopsB);
      bsscrp_self->solver->stats.pressure_flops=(double)(flopsB-flopsA); }
    //PetscPrintf( PETSC_COMM_WORLD, "\n\t* KSPSolve( ksp_S, h_hat, p )  Solve Finished in time: %lf seconds\n\n", scrSolveTime);
    bsscrp_self->solver->stats.pressure_time=scrSolveTime;

    /***************************************/
    if((hnorm < 1e-6) && (hnorm > 1e-20)){
        VecScale(h_hat,hnorm);
        VecScale(p,hnorm);
    }
    KSPRemovePressureNullspace_BSSCR(ksp, p);

    /***************************************************************************************************************/
    /***************************************************************************************************************/
    /* restore K and f for the Velocity back solve */

    found = PETSC_FALSE;
    restorek = PETSC_FALSE;
    PetscOptionsGetTruth( PETSC_NULL, "-restore_K", &restorek, &found);
    //PetscOptionsGetString( PETSC_NULL, "-restore_K", name, PETSC_MAX_PATH_LEN-1, &flg );
    if(penaltyNumber > 1e-10 && bsscrp_self->k2type){
        if(restorek){
            penaltyNumber = -penaltyNumber;
            if(f2) {
                ierr=VecAXPY(f,penaltyNumber,f2);
            }/*  f <- f +a*f2 */
            if(f3) {
                ierr=VecAXPY(f,penaltyNumber,f3);
            }/*  f <- f +a*f3 */
            ierr=MatAXPY(K,penaltyNumber,K2,DIFFERENT_NONZERO_PATTERN);
            CHKERRQ(ierr);/* Computes K = penaltyNumber*K2 + K */
            Stg_KSPSetOperators(ksp_inner, K, K, DIFFERENT_NONZERO_PATTERN);
            KisJustK=PETSC_TRUE;
        }
    }
    if(f3){ Stg_VecDestroy(&f3 ); } /* always destroy this local vector if was created */

    /* obtain solution for u */

    VecDuplicate( u, &t );
    MatMult( G, p, t);
    VecAYPX( t, -1.0, f ); /*** t <- -t + f   = f - G*p  ***/

    MatSchurComplementGetKSP( S, &ksp_inner );
    a11SingleSolveTime = MPI_Wtime();           /* ----------------------------------  Final V Solve */
    if(usePreviousGuess)
        KSPSetInitialGuessNonzero( ksp_inner, PETSC_TRUE );

    /***************************************************************************************************************/
    /***************************************************************************************************************/
    /*******   VELOCITY SOLVE   ************************************************************************************/

    /** Easier to just create a new KSP here if we want to do backsolve diffferently. (getting petsc errors now when switching from fgmres) */

    change_backsolve=PETSC_FALSE;
    PetscOptionsGetTruth( PETSC_NULL, "-change_backsolve", &change_backsolve, &found );
    if(change_backsolve){
      //Stg_KSPDestroy(&ksp_inner );
      KSPCreate(PETSC_COMM_WORLD, &ksp_new_inner);
      Stg_KSPSetOperators(ksp_new_inner, K, K, DIFFERENT_NONZERO_PATTERN);
      KSPSetOptionsPrefix(ksp_new_inner, "backsolveA11_");
      KSPSetFromOptions(ksp_new_inner); /* make sure we are setting up our solver how we want it */
      MatSchurComplementSetKSP( S, ksp_new_inner );/* need to give the Schur it's inner ksp back for when we destroy it at end */
      ksp_inner=ksp_new_inner;
    }

    if(get_flops) PetscGetFlops(&flopsA);
    /*************************************/
    /*************************************/

    backsolveSetupTime = MPI_Wtime();
    KSPSetUp(ksp_inner);
    backsolveSetupTime = MPI_Wtime() - backsolveSetupTime;
    bsscrp_self->solver->stats.velocity_backsolve_setup_time = backsolveSetupTime;

    KSPSolve(ksp_inner, t, u);
    KSPGetConvergedReason(ksp_inner, &reason ); {if (reason < 0) bsscrp_self->solver->backsolve_reason=(int)reason; }
    /*************************************/
    /*************************************/
    KSPGetIterationNumber( ksp_inner, &bsscrp_self->solver->stats.velocity_backsolve_its );
    if(get_flops) {
      PetscGetFlops(&flopsB);
      bsscrp_self->solver->stats.velocity_backsolve_flops=(double)(flopsB-flopsA); }
    a11SingleSolveTime = MPI_Wtime() - a11SingleSolveTime;              /* ------------------ Final V Solve */

    bsscrp_self->solver->stats.velocity_presolve_time=RHSSolveTime;
    bsscrp_self->solver->stats.velocity_backsolve_time=a11SingleSolveTime;
    bsscrp_self->solver->stats.velocity_total_its = bsscrp_self->solver->stats.velocity_backsolve_its +
                                                    bsscrp_self->solver->stats.velocity_presolve_its +
                                                    bsscrp_self->solver->stats.velocity_pressuresolve_its;

    flg=0;

    PetscOptionsGetString( PETSC_NULL, "-solutiondumpdir", name, PETSC_MAX_PATH_LEN-1, &flg );
    if(flg){
        sprintf(str,"%s/p%s",name,suffix);
        bsscr_writeVec( p, str, "Writing p vector in al Solver");
        sprintf(str,"%s/u%s",name, suffix);
        bsscr_writeVec( u, str, "Writing u vector in al Solver");
        sprintf(str,"%s/h_hat%s",name,suffix);
        bsscr_writeVec( h_hat, str, "Writing h_hat vector in al Solver");

    }

    found = PETSC_FALSE;
    restorek = PETSC_FALSE;
    PetscOptionsGetTruth( PETSC_NULL, "-restore_K_after_solve", &restorek, &found);
    if(penaltyNumber > 1e-10 && bsscrp_self->k2type){
        if(restorek){
            penaltyNumber = -penaltyNumber;
            ierr=MatAXPY(K,penaltyNumber,K2,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);/* Computes K = penaltyNumber*K2 + K */
            KisJustK=PETSC_TRUE;
        }
    }
    /***************************************************************************************************************/
    /***************************************************************************************************************/
    /******     SOLUTION SUMMARY      ******************************************************************************/
    bsscr_summary(bsscrp_self,ksp_S,ksp_inner,K,K2,D,G,C,u,p,f,h,t,penaltyNumber,KisJustK,
                  mgSetupTime, RHSSolveTime, scrSolveTime, a11SingleSolveTime);

    /***************************************************************************************************************/
    /***************************************************************************************************************/
    Stg_VecDestroy(&t );
    Stg_KSPDestroy(&ksp_S );
    Stg_VecDestroy(&h_hat );
    Stg_MatDestroy(&S );//This will destroy ksp_inner: also.. pcInner == pc_MG and is destroyed when ksp_inner is
    been_here = 1;
    PetscFunctionReturn(0);
}

ext: c

max_stars_repo_name: StuartRClark/mantle

max_issues_repo_name: StuartRClark/mantle

max_forks_repo_name: StuartRClark/mantle


text:
/* sample.c
 * written by Sanghyuk Moon
 * Last update: 2020 July 24
 *
 * DESCRIPTION:
 *   Generate 3D particle distribution that is consistent with prescribed density.
 *   Density is assumed to be axisymmetric and normalized such that rho(R,z) ~ [0,1].
 *   R, z coordinates are generated by rejection sampling using the Halton's pseudo-random
 *   sequence, while phi coordinates are uniformly spaced between [0,2pi].
 *
 * RETURN:
 *   Particle positions in Cartesian coordinates (x,y,z).
 *
 * REQUIREMENTS:
 *   GNU Scientific Library ($sudo apt install libgsl-dev)
 *
 * USAGE:
 *   $gcc sample.c -lgsl -lm
 *   $a.out > out.txt
 */

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_qrng.h>
#include <gsl/gsl_rng.h>

const unsigned int dim=2;

double density(double R, double z) {
  double rho, scaleh, Rmin;
  Rmin = 0.1;
  scaleh = 0.1; // Gaussian scale height of the disk
  rho = Rmin/R*exp(-z*z/scaleh/scaleh/2);
  return rho;
}

void sample_particles(long int ngas, double x1min, double x1max, double x3min, double x3max) {
  long id=0;
  double pos[dim], u4;
  double x1,x2,x3; // particle position in cylindrical coordinates
  double den; // desired gas density
  gsl_qrng *qrng;
  gsl_rng *rng;

  // GSL random number generator setup
  qrng = gsl_qrng_alloc(gsl_qrng_halton, dim);
  rng = gsl_rng_alloc(gsl_rng_ranlxd2);

  while (id < ngas) {
    gsl_qrng_get(qrng, pos);
    x1 = x1min + (x1max-x1min)*pos[0];
    x3 = x3min + (x3max-x3min)*pos[1];
    u4 = x1max*gsl_rng_uniform(rng); // uniform random number in [0, x1max]
    den = density(x1, x3);

    // accept if u4 < R*density.
    // additional weighting by R compensates the R*d\phi volume factor.
    if (u4 <= x1*den) {
      x2 = 2.0*M_PI*(double)id/(double)ngas;
      printf("%le %le %le\n", x1*cos(x2), x1*sin(x2), x3);
      id++;
    }
  }
  gsl_qrng_free(qrng);
  gsl_rng_free(rng);
}

int main() {
  long int ngas; // number of particles
  double x1min, x1max, x3min, x3max;

  // set domain info
  ngas = pow(10,6);
  x1min = 0.1;
  x1max = 1.0;
  x3min = -0.5;
  x3max = 0.5;

  // run rejection samplig
  sample_particles(ngas, x1min, x1max, x3min, x3max);
  return 0;
}


text:
#ifndef ALLVARS_H
#define ALLVARS_H

#include <stdio.h>
#include <stdbool.h>
#include <gsl/gsl_rng.h>
#include "core_simulation.h"
#include <stdint.h>

#define ABORT(sigterm)                                                  \
do {                                                                \
  printf("Error in file: %s\tfunc: %s\tline: %i\n", __FILE__, __FUNCTION__, __LINE__); \
  myexit(sigterm);                                                \
} while(0)

#define  STEPS 10       // Number of integration intervals between two snapshots
#define  MAXGALFAC 1
#define  ALLOCPARAMETER 10.0
#define  MAX_NODE_NAME_LEN 50
#define  ABSOLUTEMAXSNAPS 1000
#define STELLAR_TRACKING_TIME 100 // This is the time (Myr) we wish to track the stellar history for.

#define  MAXLEN      1024
#define  MAX_STRING_LEN      1024
#define  GRAVITY     6.672e-8
#define  SOLAR_MASS  1.989e33
#define  SOLAR_LUM   3.826e33
#define  RAD_CONST   7.565e-15
#define  AVOGADRO    6.0222e23
#define  BOLTZMANN   1.3806e-16
#define  GAS_CONST   8.31425e7
#define  C           2.9979e10
#define  PLANCK      6.6262e-27
#define  CM_PER_MPC  3.085678e24
#define  PROTONMASS  1.6726e-24
#define  HUBBLE      3.2407789e-18   /* in h/sec */

#define  SEC_PER_MEGAYEAR   3.155e13
#define  SEC_PER_YEAR       3.155e7

#define	 CUBE(x) (x*x*x)

#ifdef NDEBUG
#define XASSERT(EXP, ...)                                do{} while(0)
#else
#define XASSERT(EXP, ...)                                              \
    do { if (!(EXP)) {                                                  \
            printf("Error in file: %s\tfunc: %s\tline: %d with expression `"#EXP"'\n", __FILE__, __FUNCTION__, __LINE__); \
            printf(__VA_ARGS__);                                        \
            fflush(stdout);                                             \
            ABORT(EXIT_FAILURE);                                         \
        } \
    } while (0)
#endif

#ifdef NDEBUG
#define XPRINT(EXP, ...)                                do{} while(0)
#else
#define XPRINT(EXP, ...)                                              \
    do { if (!(EXP)) {                                                  \
            printf("Warning in file: %s\tfunc: %s\tline: %d with expression `"#EXP"'\n", __FILE__, __FUNCTION__, __LINE__); \
            printf(__VA_ARGS__);                                        \
            fflush(stdout);                                             \
        } \
    } while (0)
#endif

// This structure contains the properties used within the code
struct GALAXY
{

  int32_t FileNr;

  int   SnapNum;
  int   Type;

  int   GalaxyNr;
  int   CentralGal;
  int   HaloNr;
  int   TreeNr;
  long long  MostBoundID;

  int   mergeType;  //0=none; 1=minor merger; 2=major merger; 3=disk instability; 4=disrupt to ICS
  int   mergeIntoID;
  int   mergeIntoSnapNum;
  float   dT;
  float Age;
  // (sub)halo properties
  float Pos[3];
  float Vel[3];
  int   Len;
  float Mvir;
  float deltaMvir;
  float CentralMvir;
  float Rvir;
  float Vvir;
  float Vmax;

  // baryonic reservoirs
  float ColdGas;
  float StellarMass;
  float BulgeMass;
  float HotGas;
  float EjectedMass;
  float EjectedMassSN;
  float EjectedMassQSO;
  float BlackHoleMass;
  float ICS;

  // metals
  float MetalsColdGas;
  float MetalsStellarMass;
  float MetalsBulgeMass;
  float MetalsHotGas;
  float MetalsEjectedMass;
  float MetalsICS;

  // to calculate magnitudes
  float SfrDisk[STEPS];
  float SfrBulge[STEPS];
  float SfrDiskColdGas[STEPS];
  float SfrDiskColdGasMetals[STEPS];
  float SfrBulgeColdGas[STEPS];
  float SfrBulgeColdGasMetals[STEPS];

  // misc
  float DiskScaleRadius;
  float MergTime;
  double Cooling;
  double Heating;
  float r_heat;
  float QuasarModeBHaccretionMass;
  float TimeOfLastMajorMerger;
  float TimeOfLastMinorMerger;
  float OutflowRate;
	float TotalSatelliteBaryons;

  // infall properties
  float infallMvir;
  float infallVvir;
  float infallVmax;

  int32_t *GridType;
  int32_t *GridFoFHaloNr;
  int32_t *GridHistory;
  float   *GridColdGas;
  float   *GridHotGas;
  float   *GridEjectedMass;
  float   *GridDustColdGas;
  float   *GridDustHotGas;
  float   *GridDustEjectedMass;
  float   *GridStellarMass;
  float   *GridBHMass;
  float   *GridSFR;
  float   *GridZ;
  float   *GridFoFMass;
  float   *GridHaloMass;
  float   *EjectedFraction;
  float   *EjectedFractionSN;
  float   *EjectedFractionQSO;
  int32_t *LenHistory;
  float   *GridOutflowRate;
  float   *GridInfallRate;
  int32_t *QuasarActivity;
  int32_t *QuasarSubstep;
  float   *DynamicalTime;
  int32_t *LenMergerGal;
  float   *GridReionMod;
  float   *GridNgamma_HI;
  float   *Gridfesc;
  float   *ColdCrit;

  int IsMerged;
  double GrandSum;
  int IsMalloced;

  // Delayed SN Properties
  float reheated_mass;
  float ejected_mass;
  float mass_stars_recycled;
  float mass_metals_new;
  float NSN; // Number of supernova within a time step.  Can be fractional.
  double *SN_Stars;

  double StellarAge_Numerator;
  double StellarAge_Denominator;
  double Total_SN_SF_Time;
  double Total_SN_Stars;

  // Dust Properties
  float DustColdGas;
  float DustHotGas;
  float DustEjectedMass;

  // Quasar Activity Tracking

  int32_t QuasarActivityToggle;
  float TargetQuasarTime;
  float QuasarBoostActiveTime;
  float QuasarFractionalPhotons;

  // Stellar Age Tracking To Determine Nion.
  float *Stellar_Stars; // Keep this separate from the Stars of the delayed SN.
                        // Its possible tat we want to do tracking of stellar ages to determine Nion but not do delayed SN (or vice versa).
  double Total_Stellar_SF_Time;
  double Total_Stellar_Stars;

  // UV Magnitudes
  float *MUV;

}
*Gal, *HaloGal, *MergedGal;

// auxiliary halo data
struct halo_aux_data
{
  int DoneFlag;
  int HaloFlag;
  int NGalaxies;
  int FirstGalaxy;
}
*HaloAux;

// For the fully self consistent model we read in a list of reionization modifiers for each halo.
// These lists only contain those halos with a reionization modifier less than 0.9999.
// We will read if a list for each of the previous snapshots below the one we are currently on.
// E.g., if we are applying reionization feedback for snapshot 40, we will read in lists for snapshot 0 to 39.
// There is one list for each tree file.
struct REIONMOD_STRUCT
{

  int32_t NumLists; // Will be reionization snapshot number.
  struct REIONMOD_LIST
  {
    int32_t NHalos_Found; // Number of halos that have been matched.
    int32_t NHalos_Ionized; // Number of halos in the list.
    int64_t *HaloID; // This is a unique ID for each tree file.  The most significant 32 bits (left-most) is the tree number, and the least significant (right-most) bits is the halonr within the tree.
    float *ReionMod; // Reionization Modifer for each halo.
  }*ReionMod_List;
}*ReionList;

// Grid for doing photoionization feedback.
// This will contain multiple grids, one for each redshift.
struct GRID_STRUCT
{

  int32_t GridSize;
  uint64_t NumCellsTotal;
  int32_t NumGrids;

  double *ReionRedshift; // This is the redshift the grid cell was ionized at.
  struct PHOTO_GRID
  {
    double *PhotoRate; // Photoionization Rate (s^-1).
    int32_t valid_grid; // This will control whether data has been read in for this snapshot. Useful for the early snapshots where we don't have any ionization.
  }*PhotoGrid;
}*Grid;

struct SELFCON_GRID_STRUCT
{

  int32_t GridSize;
  int32_t NumCellsTotal;
  float Nion_HI_Total;

  float *Nion_HI;
  int32_t *GalCount;

}*SelfConGrid;

extern int    FirstFile;    // first and last file for processing
extern int    LastFile;

extern int count_onehalo;

extern int    Ntrees;      // number of trees in current file
extern int    NumGals;     // Total number of galaxies stored for current tree
extern int    MaxGals;     // Maximum number of galaxies allowed for current tree
extern int    FoF_MaxGals;

extern int    GalaxyCounter;     // unique galaxy ID for main progenitor line in tree

extern int    LastSnapShotNr;

extern char   OutputDir[512];
extern char   GalaxyOutputDir[512];
extern char   GridOutputDir[512];
extern char   RunPrefix[512];
extern char   TreeName[512];
extern char   TreeExtension[512];
extern char   SimulationDir[512];
extern char   IonizationDir[512];
extern char   FileWithSnapList[512];
extern char   PhotoionDir[512];
extern char   PhotoionName[512];
extern char   ReionRedshiftName[512];
extern int    ReionSnap;

extern int    TotHalos;
extern int    TotGalaxies;
extern int    *TreeNgals;

extern int    *FirstHaloInSnap;

extern int    *TreeNHalos;
extern int    *TreeFirstHalo;

#ifdef MPI
extern int ThisTask, NTask, nodeNameLen;
extern char *ThisNode;
#else
extern int ThisTask, NTask;
#endif

extern int IMF;

extern double Omega;
extern double OmegaLambda;
extern double PartMass;
extern double BoxSize;
extern int GridSize;
extern int self_consistent;
extern double Hubble_h;
extern double EnergySNcode, EnergySN;

// recipe flags
extern int    QuasarRecipeOn;
extern int    ReionizationOn;
extern int    SupernovaRecipeOn;
extern int    DiskInstabilityOn;
extern int    AGNrecipeOn;
extern int    SFprescription;

// recipe parameters
extern double RecycleFraction;
extern double Yield;
extern double ReIncorporationFactor;
extern double ThreshMajorMerger;
extern double BaryonFrac;
extern double SfrEfficiency;
extern double FeedbackReheatingEpsilon;
extern double FeedbackEjectionEfficiency;
extern double RadioModeEfficiency;
extern double QuasarModeEfficiency;
extern double BlackHoleGrowthRate;
extern double Reionization_z0;
extern double Reionization_zr;
extern double ThresholdSatDisruption;

// Parameters for the gridding with self_consistent
extern int fescPrescription;

extern double alpha;
extern double beta;
extern double delta;

extern double fescMH_alpha;
extern double fescMH_beta;

extern double quasar_baseline;
extern double quasar_boosted;
extern double N_dyntime;

extern double MH_low;
extern double fesc_low;
extern double MH_high;
extern double fesc_high;

extern double UnitLength_in_cm,
  UnitTime_in_s,
  UnitVelocity_in_cm_per_s,
  UnitMass_in_g,
  RhoCrit,
  UnitPressure_in_cgs,
  UnitDensity_in_cgs,
  UnitCoolingRate_in_cgs,
  UnitEnergy_in_cgs,
  UnitTime_in_Megayears,
  sage_G,
  sage_Hubble,
  a0, ar;

extern int    ListOutputSnaps[ABSOLUTEMAXSNAPS];

extern double ZZ[ABSOLUTEMAXSNAPS];
extern double AA[ABSOLUTEMAXSNAPS];
extern double Age[ABSOLUTEMAXSNAPS];

extern int    MAXSNAPS;
extern int    Snaplistlen;

extern gsl_rng *random_generator;

extern int TreeID;
extern int FileNum;
extern int neg_cell;
extern int pos_cell;

extern double ejectedmass_total;
extern double metalsejectedmass_total;

extern int MergedNr;
extern int TotMerged;
extern int *TreeNMergedgals;
extern int MaxMergedGals;     // Maximum number of galaxies allowed for current tree

extern int zeromass_count;
extern int suppression_count;
extern int previous_tree;
extern int lowmass_halo;

extern double smallest_mass;

extern double IMF_norm;
extern double IMF_slope;
extern double Eta_SNII;
extern double m_SNII;
extern int IRA;
extern int TimeResolutionSN;
extern int SN_Array_Len;
extern int Time_SFH;

extern double alpha_energy;
extern double V_energy;
extern double beta_energy;
extern double alpha_mass;
extern double V_mass;
extern double beta_mass;
extern double epsilon_mass_max;

extern int RescaleSN;
extern int mergedgal_mallocs;
extern int gal_mallocs;
extern int ismerged_mallocs;

extern int mergedgal_frees;
extern int gal_frees;
extern int ismerged_frees;

extern int *gal_to_free;

extern int outside_box;
extern int inside_box;

extern int count_Mvir;
extern int count_Len;

extern long count_gal;

extern int32_t LowSnap;
extern int32_t HighSnap;

extern float Ngamma_HI_Total;

// Delayed Supernovae Lookup Tables

extern float *IMF_massgrid_eta;
extern float *IMF_massgrid_m;
extern float m_IMFbins_low;
extern float m_IMFbins_high;
extern float m_IMFbins_delta;
extern float N_massbins;

extern float *coreburning_times;
extern float coreburning_tbins_low;
extern float coreburning_tbins_high;
extern float coreburning_tbins_delta;
extern float N_tbins;

// Ionizing Photon Lookup Tables

extern int32_t PhotonPrescription;
extern float *stars_tbins;
extern float *stars_Ngamma;
extern double TimeResolutionStellar;
extern int32_t StellarTracking_Len;
extern int32_t HaloPartCut;

// UV Magnitude

extern int32_t calcUVmag;
extern float *stars_LUV;

#endif  // #ifndef ALLVARS_H

max_stars_repo_name: jacobseiler/rsage

max_issues_repo_name: jacobseiler/rsage

max_forks_repo_name: jacobseiler/rsage


text:
static char help[] = "Solve the p-Laplacian equation in 2D using Q^1 FEM.\n"
"Implements an objective function and a residual (gradient) function, but\n"
"no Jacobian.  Defaults to p=4 and quadrature degree n=2.  Run as one of:\n"
"   ./plap -snes_fd_color             [default]\n"
"   ./plap -snes_mf\n"
"   ./plap -snes_fd                   [does not scale]\n"
"   ./plap -snes_fd_function -snes_fd [does not scale]\n"
"Uses a manufactured solution.\n"
"This is NOT a recommended example for further work because of the weird way\n"
"it handles the Dirichlet boundary values.\n\n";

#include <petsc.h>
#include "../../quadrature.h"

#define COMM PETSC_COMM_WORLD

/* cool dendritic failure:
    ./plap -snes_fd_color -snes_converged_reason -ksp_converged_reason -pc_type mg -plap_p 10.0 -da_refine 6 -snes_monitor_solution draw
succeeds in 11 snes iterations (on fine grid) with grid sequencing:
    ./plap -snes_fd_color -snes_converged_reason -ksp_converged_reason -pc_type mg -plap_p 10.0 -snes_grid_sequence 6 -snes_monitor_solution draw
*/

//STARTCTX
typedef struct {
    double    p, eps, alpha;
    int       quaddegree;
    PetscBool no_residual;
} PLapCtx;
//ENDCTX

PetscErrorCode ConfigureCtx(PLapCtx *user) {
    PetscErrorCode ierr;
    user->p = 4.0;
    user->eps = 0.0;
    user->alpha = 1.0;
    user->quaddegree = 2;
    user->no_residual = PETSC_FALSE;
    ierr = PetscOptionsBegin(COMM,"plap_","p-laplacian solver options",""); CHKERRQ(ierr);
    ierr = PetscOptionsReal("-p","exponent p with  1 <= p < infty",
                      "plap.c",user->p,&(user->p),NULL); CHKERRQ(ierr);
    if (user->p < 1.0) { SETERRQ(COMM,1,"p >= 1 required"); }
    ierr = PetscOptionsReal("-eps","regularization parameter eps",
                      "plap.c",user->eps,&(user->eps),NULL); CHKERRQ(ierr);
    ierr = PetscOptionsReal("-alpha","parameter alpha in exact solution",
                      "plap.c",user->alpha,&(user->alpha),NULL); CHKERRQ(ierr);
    ierr = PetscOptionsInt("-quaddegree","quadrature degree n (= 1,2,3 only)",
                     "plap.c",user->quaddegree,&(user->quaddegree),NULL); CHKERRQ(ierr);
    if ((user->quaddegree < 1) || (user->quaddegree > 3)) {
        SETERRQ(COMM,2,"quadrature degree n=1,2,3 only"); }
    ierr = PetscOptionsBool("-no_residual","do not set the residual evaluation function",
                      "plap.c",user->no_residual,&(user->no_residual),NULL);CHKERRQ(ierr);
    ierr = PetscOptionsEnd(); CHKERRQ(ierr);
    return 0;
}

//STARTEXACT
double Uexact(double x, double y, double alpha) {
    return 0.5 * (x+alpha)*(x+alpha) * (y+alpha)*(y+alpha);
}

double Frhs(double x, double y, PLapCtx *user) {
    const double alf = user->alpha,
                 XX = (x+alf)*(x+alf),
                 YY = (y+alf)*(y+alf),
                 D2 = XX + YY,
                 C = PetscPowScalar(XX * YY * D2, (user->p - 2.0) / 2.0),
                 gamma1 = 1.0/(x+alf) + (x+alf)/D2,
                 gamma2 = 1.0/(y+alf) + (y+alf)/D2;
    return - (user->p - 2.0) * C * (gamma1*(x+alf)*YY + gamma2*XX*(y+alf))
           - C * D2;
}
//ENDEXACT

//STARTINITIALITERATE
PetscErrorCode InitialIterateLocal(DMDALocalInfo *info, Vec u, PLapCtx *user) {
    PetscErrorCode ierr;
    const double hx = 1.0 / (info->mx+1), hy = 1.0 / (info->my+1);
    double       x, y, **au;
    int          i, j;
    ierr = DMDAVecGetArray(info->da,u,&au); CHKERRQ(ierr);
    for (j = info->ys; j < info->ys + info->ym; j++) {
        y = hy * (j + 1);
        for (i = info->xs; i < info->xs + info->xm; i++) {
            x = hx * (i + 1);
            au[j][i] = (1.0 - x) * Uexact(0.0,y,user->alpha)
                       + x * Uexact(1.0,y,user->alpha);
        }
    }
    ierr = DMDAVecRestoreArray(info->da,u,&au); CHKERRQ(ierr);
    return 0;
}
//ENDINITIALITERATE

PetscErrorCode GetUexactLocal(DMDALocalInfo *info, Vec uex, PLapCtx *user) {
    PetscErrorCode ierr;
    const double hx = 1.0 / (info->mx+1), hy = 1.0 / (info->my+1);
    double       x, y, **auex;
    int          i, j;
    ierr = DMDAVecGetArray(info->da,uex,&auex); CHKERRQ(ierr);
    for (j = info->ys; j < info->ys + info->ym; j++) {
        y = hy * (j + 1);
        for (i = info->xs; i < info->xs + info->xm; i++) {
            x = hx * (i + 1);
            auex[j][i] = Uexact(x,y,user->alpha);
        }
    }
    ierr = DMDAVecRestoreArray(info->da,uex,&auex); CHKERRQ(ierr);
    return 0;
}

//STARTFEM
static double xiL[4]  = { 1.0, -1.0, -1.0,  1.0},
              etaL[4] = { 1.0,  1.0, -1.0, -1.0};

double chi(int L, double xi, double eta) {
    return 0.25 * (1.0 + xiL[L] * xi) * (1.0 + etaL[L] * eta);
}

// evaluate v(xi,eta) on reference element using local node numbering
double eval(const double v[4], double xi, double eta) {
    double sum = 0.0;
    int    L;
    for (L=0; L<4; L++)
        sum += v[L] * chi(L,xi,eta);
    return sum;
}

typedef struct {
    double  xi, eta;
} gradRef;

gradRef dchi(int L, double xi, double eta) {
    const gradRef result = {0.25 * xiL[L]  * (1.0 + etaL[L] * eta),
                            0.25 * etaL[L] * (1.0 + xiL[L]  * xi)};
    return result;
}

// evaluate partial derivs of v(xi,eta) on reference element
gradRef deval(const double v[4], double xi, double eta) {
    gradRef sum = {0.0,0.0}, tmp;
    int     L;
    for (L=0; L<4; L++) {
        tmp = dchi(L,xi,eta);
        sum.xi += v[L] * tmp.xi;  sum.eta += v[L] * tmp.eta;
    }
    return sum;
}
//ENDFEM

//STARTTOOLS
void GetUorG(DMDALocalInfo *info, int i, int j, double **au, double *u,
             PLapCtx *user) {
    const double hx = 1.0 / (info->mx+1),  hy = 1.0 / (info->my+1),
                 x = hx * (i + 1),  y = hy * (j + 1);
    u[0] = ((i == info->mx) || (j == info->my))
             ? Uexact(x,   y,   user->alpha) : au[j][i];
    u[1] = ((i == 0)  || (j == info->my))
             ? Uexact(x-hx,y,   user->alpha) : au[j][i-1];
    u[2] = ((i == 0)  || (j == 0))
             ? Uexact(x-hx,y-hy,user->alpha) : au[j-1][i-1];
    u[3] = ((i == info->mx) || (j == 0))
             ? Uexact(x,   y-hy,user->alpha) : au[j-1][i];
}

double GradInnerProd(DMDALocalInfo *info, gradRef du, gradRef dv) {
    const double hx = 1.0 / (info->mx+1),  hy = 1.0 / (info->my+1),
                 cx = 4.0 / (hx * hx),  cy = 4.0 / (hy * hy);
    return cx * du.xi  * dv.xi + cy * du.eta * dv.eta;
}

double GradPow(DMDALocalInfo *info, gradRef du, double P, double eps) {
    return PetscPowScalar(GradInnerProd(info,du,du) + eps * eps, P / 2.0);
}
//ENDTOOLS

//STARTOBJECTIVE
double ObjIntegrandRef(DMDALocalInfo *info,
                       const double f[4], const double u[4],
                       double xi, double eta, double p, double eps) {
    const gradRef du = deval(u,xi,eta);
    return GradPow(info,du,p,eps) / p - eval(f,xi,eta) * eval(u,xi,eta);
}

PetscErrorCode FormObjectiveLocal(DMDALocalInfo *info, double **au,
                                  double *obj, PLapCtx *user) {
  PetscErrorCode ierr;
  const double hx = 1.0 / (info->mx+1),  hy = 1.0 / (info->my+1);
  const Quad1D q = gausslegendre[user->quaddegree-1];
  const int    XE = info->xs + info->xm,  YE = info->ys + info->ym;
  double       x, y, lobj = 0.0, u[4];
  int          i,j,r,s;
  MPI_Comm     com;

  // loop over all elements
  for (j = info->ys; j <= YE; j++) {
      y = hy * (j + 1);
      for (i = info->xs; i <= XE; i++) {
          x = hx * (i + 1);
          // owned elements include "right" and "top" edges of grid
          if (i < XE || j < YE || i == info->mx || j == info->my) {
              const double f[4] = {Frhs(x,   y,   user),
                                   Frhs(x-hx,y,   user),
                                   Frhs(x-hx,y-hy,user),
                                   Frhs(x,   y-hy,user)};
              GetUorG(info,i,j,au,u,user);
              // loop over quadrature points
              for (r = 0; r < q.n; r++) {
                  for (s = 0; s < q.n; s++) {
                      lobj += q.w[r] * q.w[s]
                              * ObjIntegrandRef(info,f,u,q.xi[r],q.xi[s],
                                                user->p,user->eps);
                  }
              }
          }
      }
  }
  lobj *= 0.25 * hx * hy;

  ierr = PetscObjectGetComm((PetscObject)(info->da),&com); CHKERRQ(ierr);
  ierr = MPI_Allreduce(&lobj,obj,1,MPI_DOUBLE,MPI_SUM,com); CHKERRQ(ierr);
  return 0;
}
//ENDOBJECTIVE

//STARTFUNCTION
double FunIntegrandRef(DMDALocalInfo *info, int L,
                       const double f[4], const double u[4],
                       double xi, double eta, double p, double eps) {
  const gradRef du    = deval(u,xi,eta),
                dchiL = dchi(L,xi,eta);
  return GradPow(info,du,p - 2.0,eps) * GradInnerProd(info,du,dchiL)
         - eval(f,xi,eta) * chi(L,xi,eta);
}

PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, double **au,
                                 double **FF, PLapCtx *user) {
  const double hx = 1.0 / (info->mx+1),  hy = 1.0 / (info->my+1);
  const Quad1D q = gausslegendre[user->quaddegree-1];
  const int    XE = info->xs + info->xm,  YE = info->ys + info->ym,
               li[4] = {0,-1,-1,0},  lj[4] = {0,0,-1,-1};
  double       x, y, u[4];
  int          i,j,l,r,s,PP,QQ;

  // clear residuals
  for (j = info->ys; j < YE; j++)
      for (i = info->xs; i < XE; i++)
          FF[j][i] = 0.0;

  // loop over all elements
  for (j = info->ys; j <= YE; j++) {
      y = hy * (j + 1);
      for (i = info->xs; i <= XE; i++) {
          x = hx * (i + 1);
          const double f[4] = {Frhs(x,   y,   user),
                               Frhs(x-hx,y,   user),
                               Frhs(x-hx,y-hy,user),
                               Frhs(x,   y-hy,user)};
          GetUorG(info,i,j,au,u,user);
          // loop over corners of element i,j
          for (l = 0; l < 4; l++) {
              PP = i + li[l];
              QQ = j + lj[l];
              // only update residual if we own node
              if (PP >= info->xs && PP < XE
                  && QQ >= info->ys && QQ < YE) {
                  // loop over quadrature points
                  for (r = 0; r < q.n; r++) {
                      for (s = 0; s < q.n; s++) {
                         FF[QQ][PP]
                             += 0.25 * hx * hy * q.w[r] * q.w[s]
                                * FunIntegrandRef(info,l,f,u,
                                                  q.xi[r],q.xi[s],
                                                  user->p,user->eps);
                      }
                  }
              }
          }
      }
  }
  return 0;
}
//ENDFUNCTION

int main(int argc,char **argv) {
  PetscErrorCode ierr;
  DM             da, da_after;
  SNES           snes;
  Vec            u_initial, u, u_exact;
  PLapCtx        user;
  DMDALocalInfo  info;
  double         err, hx, hy;

  PetscInitialize(&argc,&argv,NULL,help);
  ierr = ConfigureCtx(&user); CHKERRQ(ierr);

  ierr = DMDACreate2d(COMM,
               DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX,
               3,3,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da); CHKERRQ(ierr);
  ierr = DMSetFromOptions(da); CHKERRQ(ierr);
  ierr = DMSetUp(da); CHKERRQ(ierr);
  ierr = DMSetApplicationContext(da,&user);CHKERRQ(ierr);

  ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);
  hx = 1.0 / (info.mx+1);
  hy = 1.0 / (info.my+1);
  ierr = DMDASetUniformCoordinates(da,hx,1.0-hx,hy,1.0-hy,0.0,1.0); CHKERRQ(ierr);

  ierr = SNESCreate(COMM,&snes); CHKERRQ(ierr);
  ierr = SNESSetDM(snes,da); CHKERRQ(ierr);
  ierr = DMDASNESSetObjectiveLocal(da,
             (DMDASNESObjective)FormObjectiveLocal,&user); CHKERRQ(ierr);
  if (!user.no_residual) {
      ierr = DMDASNESSetFunctionLocal(da,INSERT_VALUES,
                 (DMDASNESFunction)FormFunctionLocal,&user); CHKERRQ(ierr);
  }
  ierr = SNESSetFromOptions(snes); CHKERRQ(ierr);

  ierr = DMCreateGlobalVector(da,&u_initial);CHKERRQ(ierr);
  ierr = InitialIterateLocal(&info,u_initial,&user); CHKERRQ(ierr);
  ierr = SNESSolve(snes,NULL,u_initial); CHKERRQ(ierr);
  ierr = VecDestroy(&u_initial); CHKERRQ(ierr);
  ierr = DMDestroy(&da); CHKERRQ(ierr);

  ierr = SNESGetSolution(snes,&u); CHKERRQ(ierr);
  ierr = VecDuplicate(u,&u_exact); CHKERRQ(ierr);
  ierr = SNESGetDM(snes,&da_after); CHKERRQ(ierr);
  ierr = DMDAGetLocalInfo(da_after,&info); CHKERRQ(ierr);
  ierr = PetscPrintf(COMM,
            "grid of %d x %d = %d interior nodes\n",
            info.mx,info.my,info.mx*info.my); CHKERRQ(ierr);
  ierr = GetUexactLocal(&info,u_exact,&user); CHKERRQ(ierr);
  ierr = VecAXPY(u,-1.0,u_exact); CHKERRQ(ierr);    // u <- u + (-1.0) uexact
  ierr = VecNorm(u,NORM_INFINITY,&err); CHKERRQ(ierr);
  ierr = PetscPrintf(COMM,"numerical error:  |u-u_exact|_inf = %.3e\n",
           err); CHKERRQ(ierr);

  VecDestroy(&u_exact);  SNESDestroy(&snes);
  return PetscFinalize();
}


text:
#include <stdio.h>
#include <stdlib.h>
#include <cblas.h>
#include <lapacke.h>
#include <string.h>

void eigs(double *F, double *Es, double *l, int n_eigs, int n_patch)
{
    double *eigenvalue = (double *)malloc(n_patch*sizeof(double));
    lapack_int n = n_patch;
    lapack_int info;

    printf("n_patch = %d\n", n);
    info = LAPACKE_dsyevd(LAPACK_COL_MAJOR, 'V', 'L', n, l, n, eigenvalue);
    printf("info = %d\n", info);
    memcpy(F, l, n_eigs*n_patch*sizeof(double));
    memcpy(Es, eigenvalue, n_eigs*sizeof(double));
    free(eigenvalue);
}


text:
/* Overlap-add method
 *
 * Copyright (C) 2017-2018 Jean-Pierre Cimalando.
 */

#pragma once
#include <jsl/dynarray>
#include <gsl/span>

template <class R>
class overlap_add {
public:
    overlap_add() noexcept {}
    overlap_add(size_t step, gsl::span<const R> window);

    // compute one block of overlap-add
    void process(const R *in/*[winsize]*/, R *out/*[step]*/);

    const size_t step() const
        { return step_; }

    const jsl::dynarray<R> &window() const
        { return window_; }

private:
    size_t step_ = 0;
    jsl::dynarray<R> window_;
    jsl::dynarray<R> buffer_;
};

#include "util/dsp/overlap_add.tcc"

max_stars_repo_name: jpcima/pd-externals

max_issues_repo_name: jpcima/pd-externals

max_forks_repo_name: jpcima/pd-externals


text:
/* specfunc/sincos_pi.c
 *
 * Copyright (C) 2017 Gerard Jungman, Konrad Griessinger (konradg@gmx.net)
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

/* routines for computing sin(pi*x) and cos(pi*x), respectively, with argument reduction */

#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf_trig.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_sincos_pi.h>

/* Any double precision number bigger than this is automatically an even integer. */
#define TWOBIG (2.0 / GSL_DBL_EPSILON)

/* routine computing sin(pi*x) valid for |x| <= 0.25 using a Taylor expansion around the origin and otherwise a rational approximation from the reference below. Spot-checked to give around 2e-16 relative accuracy. */
/* I. Koren and O. Zinaty. Evaluating elementary functions in a numerical
coprocessor based on rational approximations. IEEE Transactions on
Computers, Vol.39, No.8, August 1990, pp 1030-1037. */
/*
static int
sin_pi_koren(const double x, gsl_sf_result *result)
{
  result->val = 0.0;
  result->err = 0.0;
  if (16.0*fabs(x) < 1.0) {
    const double y = M_PI * x;
    const double a = y*y;
    result->val = y*(1.0 - a*(1.0 - a*(1.0 - a*(1.0 - a*(1.0 - a/110.0)/72.0)/42.0)/20.0)/6.0);
  }
  else {
    const double a0 = 1805490264.690988571178600370234394843221;
    const double a1 = -164384678.227499837726129612587952660511;
    const double a2 =    3664210.647581261810227924465160827365;
    const double a3 =     -28904.140246461781357223741935980097;
    const double a4 =         76.568981088717405810132543523682;
    const double b0 = 2298821602.638922662086487520330827251172;
    const double b1 =   27037050.118894436776624866648235591988;
    const double b2 =     155791.388546947693206469423979505671;
    const double b3 =        540.567501261284024767779280700089;
    const double t = 16.0*x*x;
    result->val = 4.0*x*(((( a4*t + a3 )*t + a2 )*t + a1 )*t + a0)/(((( t + b3 )*t + b2 )*t + b1 )*t + b0);
  }

  result->err = GSL_DBL_EPSILON*fabs(result->val);

  return GSL_SUCCESS;
}
*/

/* routine computing cos(pi*x) valid for |x| <= 0.25 using a Taylor expansion around the origin and otherwise a rational approximation from the reference below. Spot-checked to give around 2e-16 relative accuracy. */
/* I. Koren and O. Zinaty. Evaluating elementary functions in a numerical
coprocessor based on rational approximations. IEEE Transactions on
Computers, Vol.39, No.8, August 1990, pp 1030-1037. */
/*
static int
cos_pi_koren(const double x, gsl_sf_result *result)
{
  result->val = 0.0;
  result->err = 0.0;
  if (20.0*fabs(x) < 1.0) {
    const double y = M_PI * x;
    const double a = y*y;
    result->val = 1.0 - 0.5*a*(1.0 - a*(1.0 - a*(1.0 - a*(1.0 - a/90.0)/56.0)/30.0)/12.0);
  }
  else {
    const double a0 = 1090157078.174871420428849017262549038606;
    const double a1 = -321324810.993150712401352959397648541681;
    const double a2 =   12787876.849523878944051885325593878177;
    const double a3 =    -150026.206045948110568310887166405972;
    const double a4 =        538.333564203182661664319151379451;
    const double b0 = 1090157078.174871420428867295670039506886;
    const double b1 =   14907035.776643879767410969509628406502;
    const double b2 =     101855.811943661368302608146695082218;
    const double b3 =        429.772865107391823245671264489311;
    const double t = 16.0*x*x;
    result->val = (((( a4*t + a3 )*t + a2 )*t + a1 )*t + a0)/(((( t + b3 )*t + b2 )*t + b1 )*t + b0);
  }

  result->err = GSL_DBL_EPSILON*fabs(result->val);

  return GSL_SUCCESS;
}
*/

/* routine computing sin(pi*x) using a Taylor expansion around the origin and otherwise the library function. */
static int
sin_pi_taylor(const double x, gsl_sf_result *result)
{
  result->val = 0.0;
  result->err = 0.0;
  if (16.0*fabs(x) < 1.0) {
    const double y = M_PI * x;
    const double a = y*y;
    result->val = y*(1.0 - a*(1.0 - a*(1.0 - a*(1.0 - a*(1.0 - a/110.0)/72.0)/42.0)/20.0)/6.0);
  }
  else {
    result->val = sin(M_PI*x);
  }

  result->err = GSL_DBL_EPSILON*fabs(result->val);

  return GSL_SUCCESS;
}

/* routine computing sin(pi*x) using a Taylor expansion around the origin and otherwise the library function. */
static int
cos_pi_taylor(const double x, gsl_sf_result *result)
{
  result->val = 0.0;
  result->err = 0.0;
  if (20.0*fabs(x) < 1.0) {
    const double y = M_PI * x;
    const double a = y*y;
    result->val = 1.0 - 0.5*a*(1.0 - a*(1.0 - a*(1.0 - a*(1.0 - a/90.0)/56.0)/30.0)/12.0);
  }
  else {
    result->val = cos(M_PI*x);
  }

  result->err = GSL_DBL_EPSILON*fabs(result->val);

  return GSL_SUCCESS;
}

int
gsl_sf_sin_pi_e(const double x, gsl_sf_result *result)
{
  double intx = 0.0, fracx = 0.0;
  long q;
  int sign = 1, status;

  result->val = 0.0;
  result->err = 0.0;
  fracx = modf(x,&intx);
  if (fracx == 0.0) return GSL_SUCCESS;
  if(fabs(intx) >= TWOBIG) return GSL_SUCCESS; /* to be sure. Actually should be covered by the line above */

  q = ( ( (intx >= LONG_MIN) && (intx <= LONG_MAX) ) ? intx : fmod(intx, 2.0) );
  sign = ( q % 2 ? -1 : 1 );

  /* int sign = 1 - 2*((int)round(fmod(fabs(intx),2.0))); */
  if (fabs(fracx) == 0.5) { /* probably unnecessary */
    if (fracx < 0.0) sign = -sign;
    result->val = ( sign != 1 ? -1.0 : 1.0 );
    return GSL_SUCCESS;
  }
  if (fabs(fracx) > 0.5) {
    sign = -sign;
    fracx = ( fracx > 0.0 ? fracx-1.0 : fracx+1.0 );
  }

  status = 0;
  if (fracx > 0.25) {
    status = cos_pi_taylor((fracx-0.5), result);
  }
  else if (fracx < -0.25) {
    status = cos_pi_taylor((fracx+0.5), result);
    sign = -sign;
  }
  else {
    status = sin_pi_taylor(fracx, result);
  }
  if (sign != 1) result->val = -result->val;
  return status;
}

int
gsl_sf_cos_pi_e(const double x, gsl_sf_result *result)
{
  double intx = 0.0, fracx = 0.0;
  long q;
  int sign = 1, status;

  result->val = 0.0;
  result->err = 0.0;
  fracx = modf(x,&intx);
  if (fabs(fracx) == 0.5) return GSL_SUCCESS;

  if(fabs(intx) >= TWOBIG) {
    result->val = 1.0;
    return GSL_SUCCESS;
  }

  q = ( ( (intx >= LONG_MIN) && (intx <= LONG_MAX) ) ? intx : fmod(intx, 2.0) );
  sign = ( q % 2 ? -1 : 1 );

  /* int sign = 1 - 2*((int)round(fmod(fabs(intx),2.0))); */
  if (fracx == 0.0) { /* probably unnecessary */
    result->val = ( sign != 1 ? -1.0 : 1.0 );
    return GSL_SUCCESS;
  }
  if (fabs(fracx) > 0.5) {
    sign = -sign;
    fracx = ( fracx > 0.0 ? fracx-1.0 : fracx+1.0 );
  }

  status = 0;
  if (fracx > 0.25) {
    status = sin_pi_taylor((fracx-0.5), result);
    sign = -sign;
  }
  else if (fracx < -0.25) {
    status = sin_pi_taylor((fracx+0.5), result);
  }
  else {
    status = cos_pi_taylor(fracx, result);
  }
  if (sign != 1) result->val = -result->val;
  return status;
}

/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/

#include "eval.h"

double
gsl_sf_sin_pi(const double x)
{
  EVAL_RESULT(gsl_sf_sin_pi_e(x, &result));
}

double
gsl_sf_cos_pi(const double x)
{
  EVAL_RESULT(gsl_sf_cos_pi_e(x, &result));
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
//====---- Sudoku/Board_Iterators.h                                   ----====//
//
// Location aware iterator implementation for Board
//====--------------------------------------------------------------------====//
//
// The iterator can return a Location object.
// included by Board.h
//
//====--------------------------------------------------------------------====//
#pragma once

#include "Sudoku/Location.h"
#include <gsl/gsl>

#include <iterator>
#include <limits>
#include <type_traits>

#include "Board.fwd.h" // Forward declarations

#include <cassert>

namespace Sudoku
{
template<typename T, int N, bool is_const = false, bool is_reverse = false>
class Board_iterator;

//====--- Aliases --------------------------------------------------------====//
template<typename T, int N>
using const_Board_iterator = Board_iterator<T, N, true, false>;
template<typename T, int N>
using reverse_Board_iterator = Board_iterator<T, N, false, true>;
template<typename T, int N>
using const_reverse_Board_iterator = Board_iterator<T, N, true, true>;

//====--------------------------------------------------------------------====//
template<typename T, int N, bool is_const, bool is_reverse>
class Board_iterator
{
	using owner_type =
		std::conditional_t<is_const, Board<T, N> const, Board<T, N>>;
	using Location = ::Sudoku::Location<N>;

public:
	// member types
	using difference_type = std::ptrdiff_t;
	static_assert(
		full_size<N> < std::numeric_limits<difference_type>::max(),
		"Use std::ptrdiff_t for Board_iterator::difference_type.");
	using value_type        = T;
	using pointer           = std::conditional_t<is_const, T const*, T*>;
	using reference         = std::conditional_t<is_const, T const&, T&>;
	using iterator_category = std::random_access_iterator_tag;

	// Constructors
	constexpr Board_iterator() noexcept
	{ // defaults to [r]begin()
		if constexpr (is_reverse)
			elem_ = full_size<N> - 1;
	}
	explicit constexpr Board_iterator(gsl::not_null<owner_type*> owner) noexcept
		: board_(owner)
	{ // defaults to [r]begin()
		if constexpr (is_reverse)
			elem_ = (full_size<N> - 1);
	}
	explicit constexpr Board_iterator(
		gsl::not_null<owner_type*> owner, Location loc) noexcept
		: board_(owner), elem_(loc.element())
	{
	}
	// Assignment
	constexpr Board_iterator& operator=(Location const loc) noexcept
	{
		assert(is_valid(loc));
		elem_ = loc.element();
		return (*this);
	}

	// [[implicit]] Type Conversion to const_(reverse_)iterator
	// NOLINTNEXTLINE(google-explicit-constructor, hicpp-explicit-conversions)
	constexpr operator Board_iterator<T, N, true, is_reverse>() const noexcept
	{
		static_assert(!is_const);
		return Board_iterator<T, N, true, is_reverse>(
			gsl::not_null<owner_type*>{board_}, Location{elem_});
	}

	//====----------------------------------------------------------------====//
	[[nodiscard]] constexpr Location location() const noexcept
	{
		return Location{gsl::narrow_cast<int>(elem_)};
	}
	// Allows for conversion to Location
	explicit constexpr operator Location() const noexcept { return location(); }

	//====----------------------------------------------------------------====//
	[[nodiscard]] constexpr reference operator*() const noexcept
	{ // Only valid in a dereferenceable location:
		assert(board_ != nullptr && is_valid(Location{elem_}));
		return board_->operator[](location());
	}
	[[nodiscard]] constexpr pointer operator->() const noexcept
	{ // member access; equivalent to (*p).member
		return std::pointer_traits<pointer>::pointer_to(**this);
	}

	//====----------------------------------------------------------------====//
	constexpr Board_iterator& operator++() noexcept;
	constexpr Board_iterator operator++(int) noexcept;
	constexpr Board_iterator& operator--() noexcept;
	constexpr Board_iterator operator--(int) noexcept;
	constexpr Board_iterator& operator+=(const difference_type offset) noexcept;
	constexpr Board_iterator& operator-=(const difference_type offset) noexcept;

	[[nodiscard]] friend constexpr difference_type operator-(
		Board_iterator const left, Board_iterator const right) noexcept
	{ // difference
		assert(is_same_Board(left, right));
		if constexpr (is_reverse)
			return right.elem_ - left.elem_;
		else
			return left.elem_ - right.elem_;
	}

	[[nodiscard]] constexpr reference
		operator[](const difference_type offset) const noexcept
	{
		return (*(*this + offset));
	}

	//====----------------------------------------------------------------====//
	[[nodiscard]] friend constexpr bool operator==(
		Board_iterator const left, Board_iterator const right) noexcept
	{
		assert(is_same_Board(left, right));
		return is_same_Board(left, right) && left.elem_ == right.elem_;
	}
	[[nodiscard]] friend constexpr bool operator<(
		Board_iterator const left, Board_iterator const right) noexcept
	{
		assert(is_same_Board(left, right));
		if constexpr (is_reverse)
			return left.elem_ > right.elem_;
		else
			return left.elem_ < right.elem_;
	}

private:
	owner_type* board_{nullptr};
	gsl::index elem_{0};

	friend constexpr bool
		is_same_Board(Board_iterator const A, Board_iterator const B) noexcept
	{ // compare address of:
		return A.board_ == B.board_;
	}
};

//====--------------------------------------------------------------------====//
template<typename T, int N, bool C, bool is_reverse>
constexpr Board_iterator<T, N, C, is_reverse>&
	Board_iterator<T, N, C, is_reverse>::operator++() noexcept
{ // pre-increment
	assert(board_ != nullptr);
	if constexpr (is_reverse)
	{
		assert(elem_ >= 0);
		--elem_;
	}
	else
	{
		assert(elem_ < full_size<N>);
		++elem_;
	}
	return (*this);
}
template<typename T, int N, bool C, bool R>
constexpr Board_iterator<T, N, C, R> // NOLINTNEXTLINE(readability/casting)
	Board_iterator<T, N, C, R>::operator++(int) noexcept
{ // post-increment
	const Board_iterator pre{*this};
	operator++();
	return pre;
}
template<typename T, int N, bool C, bool is_reverse>
constexpr Board_iterator<T, N, C, is_reverse>&
	Board_iterator<T, N, C, is_reverse>::operator--() noexcept
{ // pre-decrement
	assert(board_ != nullptr);
	if constexpr (is_reverse)
	{
		assert(elem_ < full_size<N> - 1);
		++elem_;
	}
	else
	{
		assert(elem_ > 0);
		--elem_;
	}
	return (*this);
}
template<typename T, int N, bool C, bool R>
constexpr Board_iterator<T, N, C, R> // NOLINTNEXTLINE(readability/casting)
	Board_iterator<T, N, C, R>::operator--(int) noexcept
{ // post-decrement
	const Board_iterator pre{*this};
	operator--();
	return pre;
}
template<typename T, int N, bool C, bool is_reverse>
constexpr Board_iterator<T, N, C, is_reverse>&
	Board_iterator<T, N, C, is_reverse>::operator+=(
		const difference_type offset) noexcept
{
	assert(offset == 0 || board_ != nullptr);
	if constexpr (is_reverse)
	{
		elem_ -= offset;
		assert(elem_ >= -1);
		assert(elem_ < full_size<N>);
	}
	else
	{
		elem_ += offset;
		assert(elem_ >= 0);
		assert(elem_ <= full_size<N>);
	}
	return (*this);
}
template<typename T, int N, bool C, bool is_reverse>
constexpr Board_iterator<T, N, C, is_reverse>&
	Board_iterator<T, N, C, is_reverse>::operator-=(
		const difference_type offset) noexcept
{
	return operator+=(-offset);
}

//====--------------------------------------------------------------------====//
// Free functions (non-friend)
//====--------------------------------------------------------------------====//
template<typename T, int N, bool is_const, bool is_reverse>
[[nodiscard]] constexpr auto operator+(
	Board_iterator<T, N, is_const, is_reverse> left,
	typename Board_iterator<T, N, is_const, is_reverse>::difference_type const
		offset) noexcept
{ // itr + offset
	return (left += offset);
}
template<typename T, int N, bool is_const, bool is_reverse>
[[nodiscard]] inline constexpr auto operator+(
	typename Board_iterator<T, N, is_const, is_reverse>::difference_type const
		offset,
	Board_iterator<T, N, is_const, is_reverse> itr) noexcept
{ // offset + itr
	return (itr += offset);
}
template<typename T, int N, bool is_const, bool is_reverse>
[[nodiscard]] constexpr auto operator-(
	Board_iterator<T, N, is_const, is_reverse> left,
	typename Board_iterator<T, N, is_const, is_reverse>::difference_type const
		offset) noexcept
{ // itr - offset
	return (left += -offset);
}

template<typename T, int N, bool is_const, bool is_reverse>
[[nodiscard]] inline constexpr bool operator!=(
	Board_iterator<T, N, is_const, is_reverse> const left,
	Board_iterator<T, N, is_const, is_reverse> const right) noexcept
{
	return !(left == right);
}
template<typename T, int N, bool is_const, bool is_reverse>
[[nodiscard]] inline constexpr bool operator>(
	Board_iterator<T, N, is_const, is_reverse> const left,
	Board_iterator<T, N, is_const, is_reverse> const right) noexcept
{
	return (right < left);
}
template<typename T, int N, bool is_const, bool is_reverse>
[[nodiscard]] inline constexpr bool operator<=(
	Board_iterator<T, N, is_const, is_reverse> const left,
	Board_iterator<T, N, is_const, is_reverse> const right) noexcept
{
	return !(right < left);
}
template<typename T, int N, bool is_const, bool is_reverse>
[[nodiscard]] inline constexpr bool operator>=(
	Board_iterator<T, N, is_const, is_reverse> const left,
	Board_iterator<T, N, is_const, is_reverse> const right) noexcept
{
	return !(left < right);
}

} // namespace Sudoku

max_stars_repo_name: Farwaykorse/fwkSudoku

max_issues_repo_name: Farwaykorse/fwkSudoku

max_forks_repo_name: Farwaykorse/fwkSudoku