text:
/*
 * Dsp.c
 * Hamilton Kibbe
 * Copyright 2013 Hamilton Kibbe
 */

#include "Dsp.h"
#include "Utilities.h"
#include <string.h>
#include <float.h>
#include <math.h>

#ifdef __APPLE__
#include <Accelerate/Accelerate.h>
#elif defined(USE_BLAS)
#include <cblas.h>
#endif

/*******************************************************************************
 FloatBufferToInt16 */
Error_t
FloatBufferToInt16(signed short* dest, const float* src, unsigned length)
{
#ifdef __APPLE__
    // Use the Accelerate framework if we have it
    float scale = (float)INT16_MAX;
    float temp[length];
    vDSP_vsmul(src, 1, &scale, temp, 1, length);
    vDSP_vfix16(temp,1,dest,1,length);

#else
    // Otherwise do it manually
    unsigned i;
    const unsigned end = 4 * (length / 4);
    for (i = 0; i < end; i+=4)
    {
        dest[i] = floatToInt16(*src++);
        dest[i + 1] = floatToInt16(*src++);
        dest[i + 2] = floatToInt16(*src++);
        dest[i + 3] = floatToInt16(*src++);
    }
    for (i = end; i < length; ++i)
    {
        dest[i] = floatToInt16(*src++);
    }
#endif
    return NOERR;
}

/*******************************************************************************
 DoubleBufferToInt16 */
Error_t
DoubleBufferToInt16(signed short* dest, const double* src, unsigned length)
{
#ifdef __APPLE__
    // Use the Accelerate framework if we have it
    double scale = (float)INT16_MAX;
    double temp[length];
    vDSP_vsmulD(src, 1, &scale, temp, 1, length);
    vDSP_vfix16D(temp,1,dest,1,length);

#else
    // Otherwise do it manually
    unsigned i;
    const unsigned end = 4 * (length / 4);
    for (i = 0; i < end; i+=4)
    {
        dest[i] = floatToInt16(*src++);
        dest[i + 1] = floatToInt16(*src++);
        dest[i + 2] = floatToInt16(*src++);
        dest[i + 3] = floatToInt16(*src++);
    }
    for (i = end; i < length; ++i)
    {
        dest[i] = floatToInt16(*src++);
    }
#endif
    return NOERR;
}

/*******************************************************************************
 Int16BufferToFloat */
Error_t
Int16BufferToFloat(float* dest, const signed short* src, unsigned length)
{
#ifdef __APPLE__
    // Use the Accelerate framework if we have it
    float temp[length];
    float scale = 1.0 / (float)INT16_MAX;
    vDSP_vflt16(src,1,temp,1,length);
    vDSP_vsmul(temp, 1, &scale, dest, 1, length);

#else
    // Otherwise do it manually
    unsigned i;
    const unsigned end = 4 * (length / 4);
    for (i = 0; i < end; i+=4)
    {
        dest[i] = int16ToFloat(*src++);
        dest[i + 1] = int16ToFloat(*src++);
        dest[i + 2] = int16ToFloat(*src++);
        dest[i + 3] = int16ToFloat(*src++);
    }
    for (i = end; i < length; ++i)
    {
        dest[i] = int16ToFloat(*src++);
    }
#endif
    return NOERR;
}

/*******************************************************************************
 Int16BufferToDouble */
Error_t
Int16BufferToDouble(double* dest, const signed short* src, unsigned length)
{
#ifdef __APPLE__
    // Use the Accelerate framework if we have it
    double temp[length];
    double scale = 1.0 / (double)INT16_MAX;
    vDSP_vflt16D(src,1,temp,1,length);
    vDSP_vsmulD(temp, 1, &scale, dest, 1, length);

#else
    // Otherwise do it manually
    unsigned i;
    const unsigned end = 4 * (length / 4);
    for (i = 0; i < end; i+=4)
    {
        dest[i] = int16ToFloat(*src++);
        dest[i + 1] = int16ToFloat(*src++);
        dest[i + 2] = int16ToFloat(*src++);
        dest[i + 3] = int16ToFloat(*src++);
    }
    for (i = end; i < length; ++i)
    {
        dest[i] = int16ToFloat(*src++);
    }
#endif
    return NOERR;
}

/*******************************************************************************
 DoubleToFloat */
Error_t
DoubleToFloat(float* dest, const double* src, unsigned length)
{
#ifdef __APPLE__
    vDSP_vdpsp(src, 1, dest, 1, length);
#else
    // Otherwise do it manually
    unsigned i;
    const unsigned end = 4 * (length / 4);
    for (i = 0; i < end; i+=4)
    {
        dest[i] = (float)src[i];
        dest[i + 1] = (float)src[i + 1];
        dest[i + 2] = (float)src[i + 2];
        dest[i + 3] = (float)src[i + 3];
    }
    for (i = end; i < length; ++i)
    {
        dest[i] = (float)src[i];
    }
#endif
    return NOERR;
}

/*******************************************************************************
 Float To Double */
Error_t
FloatToDouble(double* dest, const float* src, unsigned length)
{
#ifdef __APPLE__
    vDSP_vspdp(src, 1, dest, 1, length);
#else
    // Otherwise do it manually
    unsigned i;
    const unsigned end = 4 * (length / 4);
    for (i = 0; i < end; i+=4)
    {
        dest[i] = (double)(src[i]);
        dest[i + 1] = (double)(src[i + 1]);
        dest[i + 2] = (double)(src[i + 2]);
        dest[i + 3] = (double)(src[i + 3]);
    }
    for (i = end; i < length; ++i)
    {
        dest[i] = (double)(src[i]);
    }
#endif
    return NOERR;
}

/*******************************************************************************
 FillBuffer */
Error_t
FillBuffer(float *dest, unsigned length, float value)
{
#ifdef __APPLE__
    // Use the Accelerate framework if we have it
    vDSP_vfill(&value, dest, 1, length);

#else
    // Otherwise do it manually
    unsigned i = 0;
    const unsigned end = 4 * (length / 4);
    for (i = 0; i < end; i+=4)
    {
        dest[i] = value;
        dest[i + 1] = value;
        dest[i + 2] = value;
        dest[i + 3] = value;
    }
    for (i = end; i < length; ++i)
    {
        dest[i] = value;
    }
#endif
    return NOERR;
}

/*******************************************************************************
 FillBufferD */
Error_t
FillBufferD(double *dest, unsigned length, double value)
{
#ifdef __APPLE__
    // Use the Accelerate framework if we have it
    vDSP_vfillD(&value, dest, 1, length);

#else
    // Otherwise do it manually
    unsigned i = 0;
    const unsigned end = 4 * (length / 4);
    for (i = 0; i < end; i+=4)
    {
        dest[i] = value;
        dest[i + 1] = value;
        dest[i + 2] = value;
        dest[i + 3] = value;
    }
    for (i = end; i < length; ++i)
    {
        dest[i] = value;
    }
#endif
    return NOERR;
}

/*******************************************************************************
 ClearBuffer */
Error_t
ClearBuffer(float *dest, unsigned length)
{
#ifdef __APPLE__
    // Use the Accelerate framework if we have it
    vDSP_vclr(dest, 1, length);

#else
    // Otherwise do it manually. Yes this works for floats
    memset(dest, 0, length * sizeof(float));
#endif
    return NOERR;
}

/*******************************************************************************
 ClearBufferD */
Error_t
ClearBufferD(double *dest, unsigned length)
{
#ifdef __APPLE__
    // Use the Accelerate framework if we have it
    vDSP_vclrD(dest, 1, length);

#else
    // Otherwise do it manually. Yes this works for doubles
    memset(dest, 0, length * sizeof(double));
#endif
    return NOERR;
}

/*******************************************************************************
 CopyBuffer */
Error_t
CopyBuffer(float* dest, const float* src, unsigned length)
{
    if (src != dest)
    {
#if defined(__APPLE__) || defined(USE_BLAS)
        // Use the Accelerate framework if we have it
        cblas_scopy(length, src, 1, dest, 1);
#else
        // Do it the boring way. Prefer memmove to memcpy in case src and dest
        // overlap
        memmove(dest, src, length * sizeof(float));
#endif
    }
    return NOERR;
}

/*******************************************************************************
 CopyBufferD */
Error_t
CopyBufferD(double* dest, const double* src, unsigned length)
{
    if (src != dest)
    {
#if defined(__APPLE__) || defined(USE_BLAS)
        // Use the Accelerate framework if we have it
        cblas_dcopy(length, src, 1, dest, 1);
#else
        // Do it the boring way. Prefer memmove to memcpy in case src and dest
        // overlap
        memmove(dest, src, length * sizeof(double));
#endif
    }
    return NOERR;
}

/*******************************************************************************
 CopyBufferStride */
Error_t
CopyBufferStride(float*         dest,
                 unsigned       dest_stride,
                 const float*   src,
                 unsigned       src_stride,
                 unsigned       length)
{
#if defined(__APPLE__) || defined(USE_BLAS)
    // Use the Accelerate framework if we have it
    cblas_scopy(length, src, src_stride, dest, dest_stride);
#else
    for (unsigned i = 0; i < length; ++i)
    {
        dest[i * dest_stride] = src[i * src_stride];
    }
#endif
    return NOERR;
}

/*******************************************************************************
 CopyBufferStrideD */
Error_t
CopyBufferStrideD(double*       dest,
                  unsigned      dest_stride,
                  const double* src,
                  unsigned      src_stride,
                  unsigned      length)
{
#if defined(__APPLE__) || defined(USE_BLAS)
    // Use the Accelerate framework if we have it
    cblas_dcopy(length, src, src_stride, dest, dest_stride);
#else
    for (unsigned i = 0; i < length; ++i)
    {
        dest[i * dest_stride] = src[i * src_stride];
    }
#endif
    return NOERR;
}

/*******************************************************************************
 SplitToInterleaved */
Error_t
SplitToInterleaved(float* dest, const float* real, const float* imag, unsigned length)
{
#if defined(__APPLE__)
    DSPSplitComplex in = {.realp = (float*)real, .imagp = (float*)imag};
    vDSP_ztoc(&in, 1, (DSPComplex*)dest, 2, length);

#elif defined(USE_BLAS)
    cblas_scopy(length, real, 1, dest, 2);
    cblas_scopy(length, imag, 1, dest + 1, 2);
#else
    unsigned i;
    unsigned i2;
    const unsigned end = 4 * (length / 4);
    for (i = 0; i < end; i+=4)
    {
        i2 = i * 2;
        dest[i2] = real[i];
        dest[i2 + 2] = real[i + 1];
        dest[i2 + 4] = real[i + 2];
        dest[i2 + 6] = real[i + 3];

        dest[i2 + 1] = imag[i];
        dest[i2 + 3] = imag[i + 1];
        dest[i2 + 5] = imag[i + 2];
        dest[i2 + 7] = imag[i + 3];
    }
    for (i = end; i < length; ++i)
    {
        i2 = i * 2;
        dest[i2] = real[i];
        dest[i2 + 1] = imag[i];
    }
#endif
    return NOERR;
}

Error_t
SplitToInterleavedD(double* dest, const double* real, const double* imag, unsigned length)
{
#if defined(__APPLE__)
    DSPDoubleSplitComplex in = {.realp = (double*)real, .imagp = (double*)imag};
    vDSP_ztocD(&in, 1, (DSPDoubleComplex*)dest, 2, length);

#elif defined(USE_BLAS)
    cblas_dcopy(length, real, 1, dest, 2);
    cblas_dcopy(length, imag, 1, dest + 1, 2);
#else
    unsigned i;
    unsigned i2;
    const unsigned end = 4 * (length / 4);
    for (i = 0; i < end; i+=4)
    {
        i2 = i * 2;
        dest[i2] = real[i];
        dest[i2 + 2] = real[i + 1];
        dest[i2 + 4] = real[i + 2];
        dest[i2 + 6] = real[i + 3];

        dest[i2 + 1] = imag[i];
        dest[i2 + 3] = imag[i + 1];
        dest[i2 + 5] = imag[i + 2];
        dest[i2 + 7] = imag[i + 3];
    }
    for (i = end; i < length; ++i)
    {
        i2 = i * 2;
        dest[i2] = real[i];
        dest[i2 + 1] = imag[i];
    }
#endif
    return NOERR;
}

/*******************************************************************************
 InterleavedToSplit */
Error_t
InterleavedToSplit(float*       real,
                   float*       imag,
                   const float* input,
                   unsigned     length)
{
#if defined(__APPLE__)
    DSPSplitComplex out = {.realp = real, .imagp = imag};
    vDSP_ctoz((const DSPComplex*)input, 2, &out, 1, length);

#elif defined(USE_BLAS)
    cblas_scopy(length, input, 2, real, 1);
    cblas_scopy(length, input + 1, 2, imag, 1);
#else
    unsigned i;
    unsigned i2;
    const unsigned end = 4 * (length / 4);
    for (i = 0; i < end; i+=4)
    {
        i2 = i * 2;
        real[i] = input[i2];
        real[i + 1] = input[i2 + 2];
        real[i + 2] = input[i2 + 4];
        real[i + 3] = input[i2 + 6];

        imag[i] = input[i2 + 1];
        imag[i + 1] = input[i2 + 3];
        imag[i + 2] = input[i2 + 5];
        imag[i + 3] = input[i2 + 7];
    }
    for (i = end; i < length; ++i)
    {
        i2 = i * 2;
        real[i] = input[i2];
        imag[i] = input[i2 + 1];
    }

#endif
    return NOERR;
}

Error_t
InterleavedToSplitD(double*         real,
                    double*         imag,
                    const double*   input,
                    unsigned        length)
{
#if defined(__APPLE__)
    DSPDoubleSplitComplex out = {.realp = real, .imagp = imag};
    vDSP_ctozD((const DSPDoubleComplex*)input, 2, &out, 1, length);

#elif defined(USE_BLAS)
    cblas_dcopy(length, input, 2, real, 1);
    cblas_dcopy(length, input + 1, 2, imag, 1);
#else

    unsigned i;
    unsigned i2;
    const unsigned end = 4 * (length / 4);
    for (i = 0; i < end; i+=4)
    {
        i2 = i * 2;
        real[i] = input[i2];
        real[i + 1] = input[i2 + 2];
        real[i + 2] = input[i2 + 4];
        real[i + 3] = input[i2 + 6];

        imag[i] = input[i2 + 1];
        imag[i + 1] = input[i2 + 3];
        imag[i + 2] = input[i2 + 5];
        imag[i + 3] = input[i2 + 7];
    }
    for (i = end; i < length; ++i)
    {
        i2 = i * 2;
        real[i] = input[i2];
        imag[i] = input[i2 + 1];
    }

#endif
    return NOERR;
}

/*******************************************************************************
 VectorAbs */
Error_t
VectorAbs(float *dest, const float *in, unsigned length)
{
#ifdef __APPLE__
    vDSP_vabs(in, 1, dest, 1, length);
#else
    // Otherwise do it manually
    unsigned i = 0;
    const unsigned end = 4 * (length / 4);
    for (i = 0; i < end; i+=4)
    {
        dest[i] = fabsf(in[i]);
        dest[i + 1] = fabsf(in[i + 1]);
        dest[i + 2] = fabsf(in[i + 2]);
        dest[i + 3] = fabsf(in[i + 3]);
    }
    for (unsigned i = end; i < length; ++i)

    {
        dest[i] = fabsf(in[i]);
    }
#endif
    return NOERR;
}

/*******************************************************************************
 VectorAbsD */
Error_t
VectorAbsD(double *dest, const double *in, unsigned length)
{
#ifdef __APPLE__
    vDSP_vabsD(in, 1, dest, 1, length);
#else
    // Otherwise do it manually
    unsigned i = 0;
    const unsigned end = 4 * (length / 4);
    for (i = 0; i < end; i+=4)
    {
        dest[i] = fabs(in[i]);
        dest[i + 1] = fabs(in[i + 1]);
        dest[i + 2] = fabs(in[i + 2]);
        dest[i + 3] = fabs(in[i + 3]);
    }
    for (unsigned i = end; i < length; ++i)

    {
        dest[i] = fabs(in[i]);
    }
#endif
    return NOERR;
}

/*******************************************************************************
 VectorNegate */
Error_t
VectorNegate(float          *dest,
             const float    *in,
             unsigned       length)
{
#if defined(__APPLE__)
    // Use the Accelerate framework if we have it
    vDSP_vneg(in, 1, dest, 1, length);
#else
    // Otherwise do it manually
    unsigned i;
    const unsigned end = 4 * (length / 4);
    for (i = 0; i < end; i+=4)
    {
        dest[i] = -in[i];
        dest[i + 1] = -in[i + 1];
        dest[i + 2] = -in[i + 2];
        dest[i + 3] = -in[i + 3];
    }
    for (i = end; i < length; ++i)
    {
        dest[i] = -in[i];
    }

#endif
    return NOERR;
}

/*******************************************************************************
 VectorNegateD */
Error_t
VectorNegateD(double          *dest,
              const double    *in,
              unsigned       length)
{
#ifdef __APPLE__
    // Use the Accelerate framework if we have it
    vDSP_vnegD(in, 1, dest, 1, length);

#else
    // Otherwise do it manually
    unsigned i;
    const unsigned end = 4 * (length / 4);
    for (i = 0; i < end; i+=4)
    {
        dest[i] = -in[i];
        dest[i + 1] = -in[i + 1];
        dest[i + 2] = -in[i + 2];
        dest[i + 3] = -in[i + 3];
    }
    for (i = end; i < length; ++i)
    {
        dest[i] = -in[i];
    }
#endif
    return NOERR;
}

/*******************************************************************************
 VectorSum */
float
VectorSum(const float* src, unsigned length)
{
    float res = 0.0;
#ifdef __APPLE__
    // Use the Accelerate framework if we have it
    vDSP_sve(src, 1, &res, length);
#else
    for (unsigned i = 0; i < length; ++i)
    {
        res += src[i];
    }
#endif
    return res;
}

/*******************************************************************************
 VectorSumD */
double
VectorSumD(const double* src, unsigned length)
{
    double res = 0.0;
#ifdef __APPLE__
    // Use the Accelerate framework if we have it
    vDSP_sveD(src, 1, &res, length);
#else
    for (unsigned i = 0; i < length; ++i)
    {
        res += src[i];
    }
#endif
    return res;
}

/*******************************************************************************
 VectorMax */
float
VectorMax(const float* src, unsigned length)
{
    float res = FLT_MIN;
#ifdef __APPLE__
    // Use the Accelerate framework if we have it
    vDSP_maxv(src, 1, &res, length);
#else
    for (unsigned i = 0; i < length; ++i)
    {
        if (src[i] > res)
        {
            res = src[i];
        }
    }
#endif
    return res;
}

/*******************************************************************************
 VectorMax */
double
VectorMaxD(const double* src, unsigned length)
{
    double res = DBL_MIN;
#ifdef __APPLE__
    // Use the Accelerate framework if we have it
    vDSP_maxvD(src, 1, &res, length);
#else
    for (unsigned i = 0; i < length; ++i)
    {
        if (src[i] > res)
        {
            res = src[i];
        }
    }
#endif
    return res;
}

/*******************************************************************************
 VectorMaxVI */
Error_t
VectorMaxVI(float* value, unsigned* index, const float* src, unsigned length)
{
#ifdef __APPLE__
    // Use the Accelerate framework if we have it
    vDSP_maxvi(src, 1, value, (unsigned long*)index, length);
#else
    float res = FLT_MIN;
    for (unsigned i = 0; i < length; ++i)
    {
        if (src[i] > res)
        {
            *value = res = src[i];
            *index = i;
        }
    }
#endif
    return NOERR;
}

/*******************************************************************************
 VectorMaxVID*/
Error_t
VectorMaxVID(double* value, unsigned* index, const double* src, unsigned length)
{
#ifdef __APPLE__
    // Use the Accelerate framework if we have it
    vDSP_maxviD(src, 1, value, (unsigned long*)index, length);
#else
    double res = src[0];
    for (unsigned i = 0; i < length; ++i)
    {
        if (src[i] > res)
        {
            *value = res = src[i];
            *index = i;
        }
    }
#endif
    return NOERR;
}

/*******************************************************************************
 VectorMin */
float
VectorMin(const float* src, unsigned length)
{
    float res = FLT_MAX;
#ifdef __APPLE__
    // Use the Accelerate framework if we have it
    vDSP_minv(src, 1, &res, length);
#else
    for (unsigned i = 0; i < length; ++i)
    {
        if (src[i] < res)
        {
            res = src[i];
        }
    }
#endif
    return res;
}

/*******************************************************************************
 VectorMinD */
double
VectorMinD(const double* src, unsigned length)
{
    double res = DBL_MAX;
#ifdef __APPLE__
    // Use the Accelerate framework if we have it
    vDSP_minvD(src, 1, &res, length);
#else
    for (unsigned i = 0; i < length; ++i)
    {
        if (src[i] < res)
        {
            res = src[i];
        }
    }
#endif
    return res;
}

/*******************************************************************************
 VectorMinVI */
Error_t
VectorMinVI(float* value, unsigned* index, const float* src, unsigned length)
{

#ifdef __APPLE__
    // Use the Accelerate framework if we have it
    vDSP_minvi(src, 1, value, (unsigned long*)index, length);
#else
    float res = src[0];
    for (unsigned i = 0; i < length; ++i)
    {
        if (src[i] < res)
        {
            *value = res = src[i];
            *index = i;
        }
    }
#endif
    return NOERR;
}

/*******************************************************************************
 VectorMinVID */
Error_t
VectorMinVID(double* value, unsigned* index, const double* src, unsigned length)
{

#ifdef __APPLE__
    // Use the Accelerate framework if we have it
    vDSP_minviD(src, 1, value, (unsigned long*)index, length);
#else
    double res = src[0];
    for (unsigned i = 0; i < length; ++i)
    {
        if (src[i] < res)
        {
            *value = res = src[i];
            *index = i;
        }
    }
#endif
    return NOERR;
}

/*******************************************************************************
 VectorVectorAdd */
Error_t
VectorVectorAdd(float         *dest,
                const float   *in1,
                const float   *in2,
                unsigned      length)
{
#ifdef __APPLE__
    // Use the Accelerate framework if we have it
    vDSP_vadd(in1, 1, in2, 1, dest, 1, length);

#else
    // Otherwise do it manually
    unsigned i;
    const unsigned end = 4 * (length / 4);
    for (i = 0; i < end; i+=4)
    {
        dest[i] = in1[i] + in2[i];
        dest[i + 1] = in1[i + 1] + in2[i + 1];
        dest[i + 2] = in1[i + 2] + in2[i + 2];
        dest[i + 3] = in1[i + 3] + in2[i + 3];
    }
    for (i = end; i < length; ++i)
    {
        dest[i] = in1[i] + in2[i];
    }

#endif
    return NOERR;
}

/*******************************************************************************
 VectorVectorAddD */
Error_t
VectorVectorAddD(double           *dest,
                 const double     *in1,
                 const double     *in2,
                 unsigned         length)
{
#ifdef __APPLE__
    // Use the Accelerate framework if we have it
    vDSP_vaddD(in1, 1, in2, 1, dest, 1, length);

#else
    // Otherwise do it manually
    unsigned i;
    const unsigned end = 4 * (length / 4);
    for (i = 0; i < end; i+=4)
    {
        dest[i] = in1[i] + in2[i];
        dest[i + 1] = in1[i + 1] + in2[i + 1];
        dest[i + 2] = in1[i + 2] + in2[i + 2];
        dest[i + 3] = in1[i + 3] + in2[i + 3];
    }
    for (i = end; i < length; ++i)
    {
        dest[i] = in1[i] + in2[i];
    }

#endif
    return NOERR;
}

/*******************************************************************************
 VectorVectorSub */
Error_t
VectorVectorSub(float         *dest,
                const float   *in1,
                const float   *in2,
                unsigned      length)
{
#ifdef __APPLE__
    // Use the Accelerate framework if we have it
    vDSP_vsub(in1, 1, in2, 1, dest, 1, length);

#else
    // Otherwise do it manually
    unsigned i;
    const unsigned end = 4 * (length / 4);
    for (i = 0; i < end; i+=4)
    {
        dest[i] = in2[i] - in1[i];
        dest[i + 1] = in2[i + 1] - in1[i + 1];
        dest[i + 2] = in2[i + 2] - in1[i + 2];
        dest[i + 3] = in2[i + 3] - in1[i + 3];
    }
    for (i = end; i < length; ++i)
    {
        dest[i] = in2[i] - in1[i];
    }

#endif
    return NOERR;
}

/*******************************************************************************
 VectorVectorSubD */
Error_t
VectorVectorSubD(double           *dest,
                 const double     *in1,
                 const double     *in2,
                 unsigned         length)
{
#ifdef __APPLE__
    // Use the Accelerate framework if we have it
    vDSP_vsubD(in1, 1, in2, 1, dest, 1, length);

#else
    // Otherwise do it manually
    unsigned i;
    const unsigned end = 4 * (length / 4);
    for (i = 0; i < end; i+=4)
    {
        dest[i] = in2[i] - in1[i];
        dest[i + 1] = in2[i + 1] - in1[i + 1];
        dest[i + 2] = in2[i + 2] - in1[i + 2];
        dest[i + 3] = in2[i + 3] - in1[i + 3];
    }
    for (i = end; i < length; ++i)
    {
        dest[i] = in2[i] - in1[i];
    }

#endif
    return NOERR;
}

/*******************************************************************************
 VectorScalarAdd */
Error_t
VectorScalarAdd(float           *dest,
                const float     *in1,
                const float     scalar,
                unsigned        length)
{
#ifdef __APPLE__
    // Use the Accelerate framework if we have it
    vDSP_vsadd(in1, 1, &scalar, dest, 1, length);
#else
    // Otherwise do it manually
    unsigned i;
    const unsigned end = 4 * (length / 4);
    for (i = 0; i < end; i+=4)
    {
        dest[i] = in1[i] + scalar;
        dest[i + 1] = in1[i + 1] + scalar;
        dest[i + 2] = in1[i + 2] + scalar;
        dest[i + 3] = in1[i + 3] + scalar;
    }
    for (i = end; i < length; ++i)
    {
        dest[i] = in1[i] + scalar;
    }

#endif
    return NOERR;
}

/*******************************************************************************
 VectorScalarAddD */
Error_t
VectorScalarAddD(double         *dest,
                 const double   *in1,
                 const double   scalar,
                 unsigned       length)
{
#ifdef __APPLE__
    // Use the Accelerate framework if we have it
    vDSP_vsaddD(in1, 1, &scalar, dest, 1, length);
#else
    // Otherwise do it manually
    unsigned i;
    const unsigned end = 4 * (length / 4);
    for (i = 0; i < end; i+=4)
    {
        dest[i] = in1[i] + scalar;
        dest[i + 1] = in1[i + 1] + scalar;
        dest[i + 2] = in1[i + 2] + scalar;
        dest[i + 3] = in1[i + 3] + scalar;
    }
    for (i = end; i < length; ++i)
    {
        dest[i] = in1[i] + scalar;
    }

#endif
    return NOERR;
}

/*******************************************************************************
 VectorVectorMultiply */
Error_t
VectorVectorMultiply(float          *dest,
                     const float    *in1,
                     const float    *in2,
                     unsigned       length)
{
#ifdef __APPLE__
    // Use the Accelerate framework if we have it
    vDSP_vmul(in1, 1, in2, 1, dest, 1, length);

#else
    // Otherwise do it manually
    unsigned i;
    const unsigned end = 4 * (length / 4);
    for (i = 0; i < end; i+=4)
    {
        dest[i] = in1[i] * in2[i];
        dest[i + 1] = in1[i + 1] * in2[i + 1];
        dest[i + 2] = in1[i + 2] * in2[i + 2];
        dest[i + 3] = in1[i + 3] * in2[i + 3];
    }
    for (i = end; i < length; ++i)
    {
        dest[i] = in1[i] * in2[i];
    }

#endif
    return NOERR;
}

/*******************************************************************************
 VectorVectorMultiplyD */
Error_t
VectorVectorMultiplyD(double        *dest,
                      const double  *in1,
                      const double  *in2,
                      unsigned      length)
{
#ifdef __APPLE__
    // Use the Accelerate framework if we have it
    vDSP_vmulD(in1, 1, in2, 1, dest, 1, length);

#else
    // Otherwise do it manually
    unsigned i;
    const unsigned end = 4 * (length / 4);
    for (i = 0; i < end; i+=4)
    {
        dest[i] = in1[i] * in2[i];
        dest[i + 1] = in1[i + 1] * in2[i + 1];
        dest[i + 2] = in1[i + 2] * in2[i + 2];
        dest[i + 3] = in1[i + 3] * in2[i + 3];
    }
    for (i = end; i < length; ++i)
    {
        dest[i] = in1[i] * in2[i];
    }

#endif
    return NOERR;
}

/*******************************************************************************
 VectorScalarMultiply */
Error_t
VectorScalarMultiply(float          *dest,
                     const float    *in1,
                     const float    scalar,
                     unsigned       length)
{
#ifdef __APPLE__
    // Use the Accelerate framework if we have it
    vDSP_vsmul(in1, 1, &scalar,dest, 1, length);

#else
    // Otherwise do it manually
    unsigned i;
    const unsigned end = 4 * (length / 4);
    for (i = 0; i < end; i+=4)
    {
        dest[i] = in1[i] * scalar;
        dest[i + 1] = in1[i + 1] * scalar;
        dest[i + 2] = in1[i + 2] * scalar;
        dest[i + 3] = in1[i + 3] * scalar;
    }
    for (i = end; i < length; ++i)
    {
        dest[i] = in1[i] * scalar;
    }
#endif
    return NOERR;
}

/*******************************************************************************
 VectorScalarMultiplyD */
Error_t
VectorScalarMultiplyD(double        *dest,
                      const double  *in1,
                      const double  scalar,
                      unsigned      length)
{
#ifdef __APPLE__
    // Use the Accelerate framework if we have it
    vDSP_vsmulD(in1, 1, &scalar,dest, 1, length);

#else
    // Otherwise do it manually
    unsigned i;
    const unsigned end = 4 * (length / 4);
    for (i = 0; i < end; i+=4)
    {
        dest[i] = in1[i] * scalar;
        dest[i + 1] = in1[i + 1] * scalar;
        dest[i + 2] = in1[i + 2] * scalar;
        dest[i + 3] = in1[i + 3] * scalar;
    }
    for (i = end; i < length; ++i)
    {
        dest[i] = in1[i] * scalar;
    }

#endif
    return NOERR;
}

/*******************************************************************************
 VectorVectorMix */

Error_t
VectorVectorMix(float        *dest,
                const float  *in1,
                const float  *scalar1,
                const float  *in2,
                const float  *scalar2,
                unsigned     length)
{
#ifdef __APPLE__
  // Use the Accelerate framework if we have it
  vDSP_vsmsma(in1, 1, scalar1, in2, 1, scalar2, dest, 1, length);
#else
  unsigned i;
  for (i = 0; i < length; ++i)
  {
    dest[i] = (in1[i] * (*scalar1)) + (in2[i] * (*scalar2));
  }
#endif
  return NOERR;
}

/*******************************************************************************
 VectorVectorMixD */

Error_t
VectorVectorMixD(double        *dest,
                 const double  *in1,
                 const double  *scalar1,
                 const double  *in2,
                 const double  *scalar2,
                 unsigned     length)
{
#ifdef __APPLE__
  // Use the Accelerate framework if we have it
  vDSP_vsmsmaD(in1, 1, scalar1, in2, 1, scalar2, dest, 1, length);
#else
  unsigned i;
  for (i = 0; i < length; ++i)
  {
    dest[i] = (in1[i] * (*scalar1)) + (in2[i] * (*scalar2));
  }
#endif
  return NOERR;
}

/*******************************************************************************
 VectorVectorSumScale */

Error_t
VectorVectorSumScale(float        *dest,
                     const float  *in1,
                     const float  *in2,
                     const float  *scalar,
                     unsigned     length)
{
#ifdef __APPLE__
// Use the Accelerate framework if we have it
vDSP_vasm(in1, 1, in2, 1, scalar, dest, 1, length);

#else
unsigned i;
for (i = 0; i < length; ++i)
{
  dest[i] = (in1[i] + in2[i]) * (*scalar);
}
#endif
return NOERR;
}

/*******************************************************************************
 VectorVectorSumScaleD */
Error_t
VectorVectorSumScaleD(double        *dest,
                      const double  *in1,
                      const double  *in2,
                      const double  *scalar,
                      unsigned     length)
{
#ifdef __APPLE__
  // Use the Accelerate framework if we have it
  vDSP_vasmD(in1, 1, in2, 1, scalar, dest, 1, length);

#else
  unsigned i;
  for (i = 0; i < length; ++i)
  {
    dest[i] = (in1[i] + in2[i]) * (*scalar);
  }
#endif
  return NOERR;
}

/*******************************************************************************
 VectorPower */
Error_t
VectorPower(float* dest, const float* in, float power, unsigned length)
{
#ifdef __APPLE__
    if (power == 2.0)
    {
        vDSP_vsq(in, 1, dest, 1, length);
    }
    else
    {
        unsigned i;
        const unsigned end = 4 * (length / 4);
        for (i = 0; i < end; i+=4)
        {
            dest[i] = powf(in[i], power);
            dest[i + 1] = powf(in[i + 1], power);
            dest[i + 2] = powf(in[i + 2], power);
            dest[i + 3] = powf(in[i + 3], power);
        }
        for (i = end; i < length; ++i)
        {
            dest[i] = powf(in[i], power);
        }
    }
#else
    unsigned i;
    const unsigned end = 4 * (length / 4);
    for (i = 0; i < end; i+=4)
    {
        dest[i] = powf(in[i], power);
        dest[i + 1] = powf(in[i + 1], power);
        dest[i + 2] = powf(in[i + 2], power);
        dest[i + 3] = powf(in[i + 3], power);
    }
    for (i = end; i < length; ++i)
    {
        dest[i] = powf(in[i], power);
    }
#endif
    return NOERR;

}

/*******************************************************************************
 VectorPowerD */
Error_t
VectorPowerD(double* dest, const double* in, double power, unsigned length)
{
#ifdef __APPLE__
    if (power == 2.0)
    {
        vDSP_vsqD(in, 1, dest, 1, length);
    }
    else
    {
        unsigned i;
        const unsigned end = 4 * (length / 4);
        for (i = 0; i < end; i+=4)
        {
            dest[i] = pow(in[i], power);
            dest[i + 1] = pow(in[i + 1], power);
            dest[i + 2] = pow(in[i + 2], power);
            dest[i + 3] = pow(in[i + 3], power);
        }
        for (i = end; i < length; ++i)
        {
            dest[i] = pow(in[i], power);
        }
    }
#else
    unsigned i;
    const unsigned end = 4 * (length / 4);
    for (i = 0; i < end; i+=4)
    {
        dest[i] = pow(in[i], power);
        dest[i + 1] = pow(in[i + 1], power);
        dest[i + 2] = pow(in[i + 2], power);
        dest[i + 3] = pow(in[i + 3], power);
    }
    for (i = end; i < length; ++i)
    {
        dest[i] = pow(in[i], power);
    }
#endif
    return NOERR;

}

/*******************************************************************************
 Convolve */
Error_t
Convolve(float       *in1,
         unsigned    in1_length,
         float       *in2,
         unsigned    in2_length,
         float       *dest)
{

    unsigned resultLength = in1_length + (in2_length - 1);
#ifdef __APPLE__
    //Use Native vectorized convolution function if available
    float    *in2_end = in2 + (in2_length - 1);
    unsigned signalLength = (in2_length + resultLength);

    float padded[signalLength];

    //float zero = 0.0;
    ClearBuffer(padded, signalLength);

    // Pad the input signal with (filter_length - 1) zeros.
    cblas_scopy(in1_length, in1, 1, (padded + (in2_length - 1)), 1);
    vDSP_conv(padded, 1, in2_end, -1, dest, 1, resultLength, in2_length);

#else
    // Use (boring, slow) canonical implementation
    unsigned i;
    for (i = 0; i <resultLength; ++i)
    {
        unsigned kmin, kmax, k;
        dest[i] = 0;

        kmin = (i >= (in2_length - 1)) ? i - (in2_length - 1) : 0;
        kmax = (i < in1_length - 1) ? i : in1_length - 1;
        for (k = kmin; k <= kmax; k++)
        {
            dest[i] += in1[k] * in2[i - k];
        }
    }

#endif
    return NOERR;
}

/*******************************************************************************
 ConvolveD */
Error_t
ConvolveD(double    *in1,
          unsigned  in1_length,
          double    *in2,
          unsigned  in2_length,
          double    *dest)
{

    unsigned resultLength = in1_length + (in2_length - 1);

#ifdef __APPLE__
    //Use Native vectorized convolution function if available
    double    *in2_end = in2 + (in2_length - 1);
    double signalLength = (in2_length + resultLength);

    // So there's some hella weird requirement that the signal input to
    //vDSP_conv has to be larger than (result_length + filter_length - 1),
    // (the output vector length) and it has to be zero-padded. What. The. Fuck!
    double padded[(unsigned)ceil(signalLength)];

    //float zero = 0.0;
    FillBufferD(padded, signalLength, 0.0);

    // Pad the input signal with (filter_length - 1) zeros.
    cblas_dcopy(in1_length, in1, 1, (padded + (in2_length - 1)), 1);
    vDSP_convD(padded, 1, in2_end, -1, dest, 1, resultLength, in2_length);

#else
    // Use (boring, slow) canonical implementation
    unsigned i;
    for (i = 0; i <resultLength; ++i)
    {
        unsigned kmin, kmax, k;
        dest[i] = 0;

        kmin = (i >= (in2_length - 1)) ? i - (in2_length - 1) : 0;
        kmax = (i < in1_length - 1) ? i : in1_length - 1;
        for (k = kmin; k <= kmax; k++)
        {
            dest[i] += in1[k] * in2[i - k];
        }
    }

#endif
    return NOERR;
}

/*******************************************************************************
 VectorDbConvert */
Error_t
VectorDbConvert(float* dest,
                const float* in,
                unsigned length)
{
#ifdef __APPLE__
    float one = 1.0;
    vDSP_vdbcon(in, 1, &one, dest, 1, length, 1);
#else
    for (unsigned i = 0; i < length; ++i)
    {
        dest[i] = AMP_TO_DB(in[i]);
    }
#endif
    return NOERR;
}

/*******************************************************************************
 VectorDbConvertD */
Error_t
VectorDbConvertD(double*        dest,
                 const double*  in,
                 unsigned       length)
{
#ifdef __APPLE__
    double one = 1.0;
    vDSP_vdbconD(in, 1, &one, dest, 1, length, 1);
#else
    for (unsigned i = 0; i < length; ++i)

    {
        dest[i] = AMP_TO_DBD(in[i]);
    }
#endif
    return NOERR;
}

/*******************************************************************************
 ComplexMultiply */
Error_t
ComplexMultiply(float*          re,
                float*          im,
                const float*    re1,
                const float*    im1,
                const float*    re2,
                const float*    im2,
                unsigned        length)
{
#if defined(__APPLE__)
    DSPSplitComplex in1 = {.realp = (float*)re1, .imagp = (float*)im1};
    DSPSplitComplex in2 = {.realp = (float*)re2, .imagp = (float*)im2};
    DSPSplitComplex out = {.realp = re, .imagp = im};
    vDSP_zvmul(&in1, 1, &in2, 1, &out, 1, length, 1);
#else

    for (unsigned i = 0; i < length; ++i)
    {
        float ire1 = re1[i];
        float iim1 = im1[i];
        float ire2 = re2[i];
        float iim2 = im2[i];
        re[i] = (ire1 * ire2 - iim1 * iim2);
        im[i] = (iim1 * ire2 + iim2 * ire1);
    }
#endif
    return NOERR;
}

/*******************************************************************************
 ComplexMultiplyD */
Error_t
ComplexMultiplyD(double*        re,
                 double*        im,
                 const double*  re1,
                 const double*  im1,
                 const double*  re2,
                 const double*  im2,
                 unsigned       length)
{
#if defined(__APPLE__)
    DSPDoubleSplitComplex in1 = {.realp = (double*)re1, .imagp = (double*)im1};
    DSPDoubleSplitComplex in2 = {.realp = (double*)re2, .imagp = (double*)im2};
    DSPDoubleSplitComplex out = {.realp = re, .imagp = im};
    vDSP_zvmulD(&in1, 1, &in2, 1, &out, 1, length, 1);
#else
    for (unsigned i = 0; i < length; ++i)
    {
        double ire1 = re1[i];
        double iim1 = im1[i];
        double ire2 = re2[i];
        double iim2 = im2[i];
        re[i] = (ire1 * ire2 - iim1 * iim2);
        im[i] = (iim1 * ire2 + iim2 * ire1);
    }
#endif
    return NOERR;
}

/*******************************************************************************
 VectorRectToPolar */
Error_t
VectorRectToPolar(float*        magnitude,
                  float*        phase,
                  const float*  real,
                  const float*  imaginary,
                  unsigned      length)
{
#ifdef __APPLE__
    float dest[2*length];
    SplitToInterleaved(dest, real, imaginary, length);
    vDSP_polar(dest, 2, dest, 2, length);
    InterleavedToSplit(magnitude, phase, dest, length);
#else
    unsigned i;
    const unsigned end = 4 * (length / 4);
    for (i = 0; i < end; i+=4)
    {
        RectToPolar(real[i], imaginary[i], &magnitude[i], &phase[i]);
        RectToPolar(real[i+1], imaginary[i+1], &magnitude[i+1], &phase[i+1]);
        RectToPolar(real[i+2], imaginary[i+2], &magnitude[i+2], &phase[i+2]);
        RectToPolar(real[i+3], imaginary[i+3], &magnitude[i+3], &phase[i+3]);
    }
    for (i = end; i < length; ++i)
    {
        RectToPolar(real[i], imaginary[i], &magnitude[i], &phase[i]);
    }
#endif
    return NOERR;
}

/*******************************************************************************
 VectorRectToPolarD */
Error_t
VectorRectToPolarD(double*          magnitude,
                   double*          phase,
                   const double*    real,
                   const double*    imaginary,
                   unsigned         length)
{
#ifdef __APPLE__
    double dest[2*length];
    SplitToInterleavedD(dest, real, imaginary, length);
    vDSP_polarD(dest, 2, dest, 2, length);
    InterleavedToSplitD(magnitude, phase, dest, length);
#else
    unsigned i;
    const unsigned end = 4 * (length / 4);
    for (i = 0; i < end; i+=4)
    {
        RectToPolarD(real[i], imaginary[i], &magnitude[i], &phase[i]);
        RectToPolarD(real[i+1], imaginary[i+1], &magnitude[i+1], &phase[i+1]);
        RectToPolarD(real[i+2], imaginary[i+2], &magnitude[i+2], &phase[i+2]);
        RectToPolarD(real[i+3], imaginary[i+3], &magnitude[i+3], &phase[i+3]);
    }
    for (i = end; i < length; ++i)
    {
        RectToPolarD(real[i], imaginary[i], &magnitude[i], &phase[i]);
    }
#endif
    return NOERR;
}

/*******************************************************************************
 MeanSquare */
float
MeanSquare(const float* data, unsigned length)
{
    float result = 0.;
#ifdef __APPLE__
    vDSP_measqv(data, 1, &result, length);
#else
    float scratch[length];
    VectorPower(scratch, data, 2, length);
    result = VectorSum(scratch, length) / length;

#endif
    return result;
}

/*******************************************************************************
 MeanSquareD */
double
MeanSquareD(const double* data, unsigned length)
{
    double result = 0.;
#ifdef __APPLE__
    vDSP_measqvD(data, 1, &result, length);
#else
    double scratch[length];
    VectorPowerD(scratch, data, 2, length);
    result = VectorSumD(scratch, length) / length;
#endif
    return result;
}

max_stars_repo_name: akshaynannaware/FxDSP

max_issues_repo_name: hamiltonkibbe/FutureTechAudio

max_forks_repo_name: hamiltonkibbe/FutureTechAudio


text:
/// @file
/// @brief shared pointers to gsl types
/// @details We use gsl for linear algebra routines. However, the interface is C-like and is based on
/// raw pointers. The file contains code to leverage on boost shared pointer and use it with gsl types
/// such as vector and matrix.
///
/// @copyright (c) 2008 CSIRO
/// Australia Telescope National Facility (ATNF)
/// Commonwealth Scientific and Industrial Research Organisation (CSIRO)
/// PO Box 76, Epping NSW 1710, Australia
/// atnf-enquiries@csiro.au
///
/// This file is part of the ASKAP software distribution.
///
/// The ASKAP software distribution 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 2 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
///
/// @author Max Voronkov <maxim.voronkov@csiro.au>

#ifndef ASKAP_SCIMATH_SHARED_GSL_TYPES_H
#define ASKAP_SCIMATH_SHARED_GSL_TYPES_H

// boost includes
#include <boost/shared_ptr.hpp>
#include <boost/static_assert.hpp>

// gsl includes
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_eigen.h>

// own includes
#include "askap/AskapError.h"

namespace askap {

namespace utility {

// typedefs for shared pointers, we can add more here when needed

/// @brief shared pointer to gsl vector
/// @ingroup utils
typedef boost::shared_ptr<gsl_vector> SharedGSLVector;

/// @brief shared pointer to gsl matrix
/// @ingroup utils
typedef boost::shared_ptr<gsl_matrix> SharedGSLMatrix;

/// @brief custom deleter for gsl types
/// @details This is a custom deleter class to allow deallocation of
/// a gsl object. Exact operation should be present in specialised types
/// @ingroup utils
template<typename GSLType> struct CustomGSLDeleter {
    /// @brief this method frees the object
    /// @param[in] obj pointer to gsl object
    void operator()(GSLType *obj) const { BOOST_STATIC_ASSERT_MSG(sizeof(GSLType) == 0,
                            "Generic version is not supposed to be used - specialise it for your type"); }
};

/// @brief specialised deleter for a gsl vector
/// @ingroup utils
template<> struct CustomGSLDeleter<gsl_vector> {
    /// @brief this method frees the object
    /// @param[in] obj pointer to gsl vector
    void operator()(gsl_vector *obj) const;
};

/// @brief specialised deleter for a gsl matrix
/// @ingroup utils
template<> struct CustomGSLDeleter<gsl_matrix> {
    /// @brief this method frees the object
    /// @param[in] obj pointer to gsl matrix
    void operator()(gsl_matrix *obj) const;
};

/// @brief specialised deleter for a gsl eigen_symmv_workspace
/// @ingroup utils
template<> struct CustomGSLDeleter<gsl_eigen_symmv_workspace> {
    /// @brief this method frees the object
    /// @param[in] obj pointer to gsl eigen_symmv_workspace
    void operator()(gsl_eigen_symmv_workspace *obj) const;
};

/// @brief templated helper method to wrap newly allocated gsl object
/// @details This method processes a pointer to a gsl object and transfers the ownership (i.e.
/// responsibility for deallocation to the boost shared pointer. Templates are used to automatically
/// deduce the object type and add an appropriate deleter.
/// @note the pointer passed as a parameter is tested to be non-zero. Otherwise an exception is thrown.
/// @param[in] obj a pointer to gsl object
/// @return boost shared pointer to the object
template<typename GSLType>
inline boost::shared_ptr<GSLType> createGSLObject(GSLType *obj)
   { ASKAPASSERT(obj); return boost::shared_ptr<GSLType>(obj, CustomGSLDeleter<GSLType>()); }

/// @brief helper method to allocate gsl vector
/// @details This method allocates a gsl vector of requested length and returns a
/// shared pointer.
/// @param[in] size size of the vector
/// @return shared pointer to gsl vector
inline SharedGSLVector createGSLVector(size_t size) { return createGSLObject(gsl_vector_alloc(size));}

/// @brief helper method to allocate gsl matrix
/// @details This method allocates a gsl matrix of requested shape and returns a
/// shared pointer.
/// @param[in] nrow number of rows
/// @param[in] ncol number of columns
/// @return shared pointer to gsl matrix
inline SharedGSLMatrix createGSLMatrix(size_t nrow, size_t ncol) { return createGSLObject(gsl_matrix_alloc(nrow,ncol));}

} // namespace utility

} // namespace askap

#endif // #ifndef ASKAP_SCIMATH_SHARED_GSL_TYPES_H

ext: h

max_stars_repo_name: rtobar/askapsoft

max_issues_repo_name: ATNF/askapsoft

max_forks_repo_name: ATNF/askapsoft


text:
/*
 * VGMTrans (c) 2002-2019
 * Licensed under the zlib license,
 * refer to the included LICENSE.txt file
 */

#pragma once

#include <string>
#include <sstream>
#include <cassert>
#include <cmath>
#include <cstring>
#include <iomanip>
#include <climits>
#include <algorithm>
#include <vector>
#include <gsl-lite.hpp>

#include "helper.h"

#define VERSION "1.0.3"

/* Type aliases to save some typing */
using size_t = std::size_t;

using s8 = std::int8_t;
using s16 = std::int16_t;
using s32 = std::int32_t;
using s64 = std::int64_t;
using sptr = std::intptr_t;

using u8 = std::uint8_t;
using u16 = std::uint16_t;
using u32 = std::uint32_t;
using u64 = std::uint64_t;
using uptr = std::uintptr_t;

std::string removeExtFromPath(const std::string &s);

std::string StringToUpper(std::string myString);
std::string StringToLower(std::string myString);

uint32_t StringToHex(const std::string &str);

std::string ConvertToSafeFileName(const std::string &str);

inline int CountBytesOfVal(uint8_t *buf, uint32_t numBytes, uint8_t val) {
    int count = 0;
    for (uint32_t i = 0; i < numBytes; i++)
        if (buf[i] == val)
            count++;
    return count;
}

struct SizeOffsetPair {
    uint32_t size;
    uint32_t offset;

    SizeOffsetPair() : size(0), offset(0) {}

    SizeOffsetPair(uint32_t offset_, uint32_t size_) : size(size_), offset(offset_) {}
};

char *GetFileWithBase(const char *f, const char *newfile);

std::vector<char> zdecompress(gsl::span<const char> data);
std::vector<char> zdecompress(std::vector<char> &data);

max_stars_repo_name: sykhro/vgmtrans-qt

max_issues_repo_name: sykhro/vgmtrans-qt

max_forks_repo_name: sykhro/vgmtrans-qt


text:
///
/// @file
///
/// @author Angelika Schwarz (angies@cs.umu.se), Umeå University
///
/// @internal LICENSE
///
/// Copyright (c) 2019-2020, Umeå Universitet
///
/// Redistribution and use in source and binary forms, with or without
/// modification, are permitted provided that the following conditions are met:
///
/// 1. Redistributions of source code must retain the above copyright notice,
///    this list of conditions and the following disclaimer.
///
/// 2. Redistributions in binary form must reproduce the above copyright notice,
///    this list of conditions and the following disclaimer in the documentation
///    and/or other materials provided with the distribution.
///
/// 3. Neither the name of the copyright holder nor the names of its
///    contributors may be used to endorse or promote products derived from this
///    software without specific prior written permission.
///
/// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
/// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
/// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
/// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
/// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
/// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
/// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
/// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
/// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
/// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
/// POSSIBILITY OF SUCH DAMAGE.
///

#include <starneig_config.h>
#include <starneig/configuration.h>

#include "typedefs.h"
#include "cpu.h"
// #include "tasks.h"  // <--- moved to misc/
#include "robust.h"

#include "../../common/common.h"
#include "../../common/tiles.h"
#include <starpu.h>
#include <cblas.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <float.h>

static double vec_real_infnorm(int n, const double *x)
{
    double norm = 0.0;
    for (int i = 0; i < n; ++i) {
        double abs = fabs(x[i]);
        if (abs > norm)
            norm = abs;
    }

    return norm;
}

static double vec_cmplx_infnorm(int n, const double *x_re, const double *x_im)
{
    double norm = 0.0;
    for (int i = 0; i < n; ++i) {
        // Compute len = sqrt(x_re[i]*x_re[i]+x_im[i]*x_im[i]) robustly.
        double maxabs = MAX(fabs(x_re[i]), fabs(x_im[i]));
        double len = maxabs*sqrt(  (x_re[i]/maxabs)*(x_re[i]/maxabs)
                                 + (x_im[i]/maxabs)*(x_im[i]/maxabs));

        if (len > norm)
            norm = len;
    }

    return norm;
}

static double mat_infnorm(int m, int n, const double *A, int ldA)
{
#define A(i,j) A[(i) + (j) * (size_t)ldA]

    double *rowsums = calloc(m, sizeof(double));

    for (int j = 0; j < n; ++j)
        for (int i = 0; i < m; i++)
            rowsums[i] += fabs(A(i,j));

    double norm = rowsums[0];
    for (int i = 1; i < m; i++)
        if (rowsums[i] > norm)
            norm = rowsums[i];

    free(rowsums);

    return norm;

#undef A
}

static void find_max(
    int num_rows, int num_selected, int n,
    const double *restrict const X, int ldX,
    const int *restrict const lambda_type, const int *restrict const selected,
    double *restrict emax)
{
    int si = num_selected-1;
    for (int ki = n-1; ki >= 0; ki--) {
        if (!selected[ki]) {
            // Proceed with the next eigenvalue.
            if (lambda_type[ki] == 1) { // CMPLX
                // A complex conjugate pair of eigenvalues is not selected, so
                // skip the next diagonal entry.
                ki--;
            }
        }
        else { // ki-th eigenvalue is selected
            if (lambda_type[ki] == 0) { // REAL
                // Locate eigenvector to ki-th eigenvalue.
               const double *X_re = X+si*(size_t)ldX;

               // Find max entry.
               emax[si] = vec_real_infnorm(num_rows, X_re);
               si--;
            }
            else { // CMPLX
                // Locate real and imaginary part of complex eigenvector for
                // complex conjugate pair of eigenvalues (ki, ki-1).
                const double *X_re = X+(si-1)*(size_t)ldX;
                const double *X_im = X+si*(size_t)ldX;

                // Find max entry.
                emax[si] = 0.0;
                for (int i = 0; i < num_rows; i++)
                    if (fabs(X_re[i])+fabs(X_im[i]) > emax[si])
                        emax[si] = fabs(X_re[i])+fabs(X_im[i]);

                // Duplicate max entry for real and imaginary column.
                emax[si-1] = emax[si];

                ki--;
                si -= 2;
            }
        }
    }
}

void starneig_eigvec_std_unify_scaling(int num_tiles, int *first_row, int *first_col,
    scaling_t *restrict scales,
    double *restrict X, int ldX,
    const int *restrict lambda_type, const int *restrict selected)
{
    // TODO: Is it possible to add a #pragma parallel for here?

    const int num_selected = first_col[num_tiles];

    //
    // Compute the most constraining scaling factor.
    //
    scaling_t *smin = (scaling_t *) malloc(num_selected*sizeof(scaling_t));
    starneig_eigvec_std_init_scaling_factor(num_selected, smin);

    starneig_eigvec_std_find_smallest_scaling(num_tiles, num_selected, scales, smin);

#ifndef STARNEIG_ENABLE_INTEGER_SCALING

    //
    // Check if the range of double-precision scaling factors was sufficient
    // and replace flushed entries with 1/Omega to avoid NaNs.
    //

    for (int i = 0; i < num_selected; i++) {
        if (smin[i] == 0.0) {
            if (lambda_type[i] == 1) { // CMPLX
                starneig_warning("Scaling factor of complex eigenvector "
                                 "X(:,%d:%d) was flushed. Rerun with integer "
                                 "scaling factors.", i, i+1);
                smin[i] = DBL_MIN;
                smin[i+1] = DBL_MIN;
                i++;
            }
            else { // REAL
                starneig_warning("Scaling factor of real eigenvector X(:,%d) "
                                 "was flushed. Rerun with integer "
                                 "scaling factors.", i);
                smin[i] = DBL_MIN;

            }
        }
    }

    for (int i = 0; i < (size_t)num_tiles*num_selected;i++)
        if (scales[i] == 0.0)
            scales[i] = DBL_MIN;

#endif

    double *tmp = (double *) malloc(num_selected * sizeof(double));
    double *emax = (double *) malloc(num_selected * sizeof(double));
    memset(emax, 0.0, num_selected * sizeof(double));

#define scales(col, tilerow) scales[(col) + (tilerow) * (size_t)num_selected]

    for (int blkj = 0; blkj < num_tiles; blkj++) {
        for (int blki = 0; blki <= blkj; blki++) {
            const int num_rows = first_row[blki+1]-first_row[blki];
            const int num_sel = first_col[blkj+1]-first_col[blkj];
            const int num_cols = first_row[blkj+1]-first_row[blkj];

            double *tile = X+(size_t)first_col[blkj]*ldX+first_row[blki];
            const int *lambda_type_tile = lambda_type+first_row[blkj];
            const int *selected_tile = selected+first_row[blkj];
            double *tmp_tile = tmp+first_col[blkj];
            memset(tmp_tile, 0.0, num_sel*sizeof(double));

            find_max(num_rows, num_sel, num_cols,
                tile, ldX, lambda_type_tile, selected_tile, tmp_tile);

            // Reduce to maximum normalization factor.
            for (int j = first_col[blkj]; j < first_col[blkj+1]; j++) {
                // Compute normalization factor simulating consistent scaling.
               double s = starneig_eigvec_std_compute_upscaling(smin[j], scales(j, blki));
               emax[j] = MAX(s * tmp[j], emax[j]);
            }
        }

        // Apply scaling.
        for (int blki = 0; blki <= blkj; blki++) {
            const int num_rows = first_row[blki+1]-first_row[blki];
            const int num_cols = first_col[blkj+1]-first_col[blkj];

            // Scale column.
            for (int j = 0; j < num_cols; j++) {
                // The current column index.
                size_t col = (size_t)first_col[blkj]+j;

                // The current column.
                double *x = X+col*ldX+first_row[blki];
                double s =
                    starneig_eigvec_std_compute_upscaling(smin[col], scales(col, blki));

                // Avoid oo.
                if (isinf(s))
                    s = DBL_MIN;

                for (int i = 0; i < num_rows; i++)
                    x[i] = (s*x[i])/emax[col];
            }
        }
    }

    free(tmp);
    free(emax);
#undef scales
}

void starneig_eigvec_std_cpu_bound_DM(void *buffers[], void *cl_args)
{
    struct packing_info packing_info;
    starpu_codelet_unpack_args(cl_args, &packing_info);

    // extract tile dimensions from the packing information struct
    int m = packing_info.rend - packing_info.rbegin;
    int n = packing_info.cend - packing_info.cbegin;

    int k = 0;

    double *norm = (double *) STARPU_VARIABLE_GET_PTR(buffers[k]);
    k++;

    double *T = (double *) STARPU_MATRIX_GET_PTR(buffers[k]);
    int ldT = STARPU_MATRIX_GET_LD(buffers[k]);
    k++;

    struct starpu_matrix_interface **A_i =
        (struct starpu_matrix_interface **)buffers + k;
    k += packing_info.handles;

    starneig_join_diag_window(&packing_info, ldT, A_i, T, 0);

    *norm = mat_infnorm(m, n, T, ldT);
}

void starneig_eigvec_std_cpu_bound(void *buffers[], void *cl_args)
{
    double *T = (double *) STARPU_MATRIX_GET_PTR(buffers[0]);
    int ldT = STARPU_MATRIX_GET_LD(buffers[0]);
    int m = STARPU_MATRIX_GET_NX(buffers[0]);
    int n = STARPU_MATRIX_GET_NY(buffers[0]);

    double *tnorm = (double *) STARPU_VARIABLE_GET_PTR(buffers[1]);
    double ub = mat_infnorm(m, n, T, ldT);
    *tnorm = ub;
}

static void backsolve(
    int n, const double *restrict T, int ldT, double tnorm,
    double *restrict X, int ldX, scaling_t *restrict const scales,
    double *restrict Xnorms, const int *restrict lambda_type,
    const int *restrict selected, int num_selected,
    double smlnum, int *restrict infos)
{
#define T(i,j) T[(i) + (j) * (size_t)ldT]

    // Solve T * X = X * \Lambda.

    // The i-th selected eigenvalue.
    int si = num_selected-1;

    // Loop over eigenvalues from bottom to top.
    for (int ki = n-1; ki >= 0; ki--) {
        if (!selected[ki]) {
            // Proceed with the next eigenvalue.
            if (lambda_type[ki] == 1) { // CMPLX
                // A complex conjugate pair of eigenvalues is not selected,
                // so skip the next diagonal entry.
                ki--;
            }
        }
        else { // ki-th eigenvalue is selected

            // Error flag.
            int info = 0;

            if (lambda_type[ki] == 0) {
                // Compute a real right eigenvector.
                double lambda = T(ki,ki);

                // Locate eigenvector to ki-th eigenvalue.
                double *X_re = X+si*(size_t)ldX;

                // Locate corresponding scaling factor.
                scaling_t *beta = scales+si;

                // Critical threshold to detect unsafe divisions.
                const double smin = MAX(DBL_EPSILON/2*fabs(lambda), smlnum);

                // Form right-hand side.
                X_re[ki] = 1.0;
                for (int i = 0; i < ki; i++)
                    X_re[i] = -T(i, ki);
                for (int i = ki+1; i < n; i++)
                    X_re[i] = 0.0;

                // Compute norm of entire vector.
                double norm = vec_real_infnorm(ki+1, X_re);

                // Solve the upper quasi-triangular system.
                // (T(0:ki-1,0:ki-1) - lambda I) \ X.
                for (int j = ki-1; j >= 0; j--) {
                    // if next block is 1-by-1 diagonal block:
                    if (lambda_type[j] == 0) { // REAL
                        scaling_t phi;
                        starneig_eigvec_std_init_scaling_factor(1, &phi);
                        info |= starneig_eigvec_std_solve_1x1_real_system(
                            smin, T(j,j), lambda, X_re+j, &phi);
                        starneig_eigvec_std_update_global_scaling(beta, phi);

                        // Scale remaining parts of the vector.
                        starneig_eigvec_std_scale(j, X_re, &phi);
                        starneig_eigvec_std_scale(ki-(j+1), X_re+(j+1), &phi);
                        starneig_eigvec_std_update_norm(&norm, phi);

                        // Protect against overflow in the linear update.
                        phi =
                            starneig_eigvec_std_protect_update(tnorm, fabs(X_re[j]), norm);
                        starneig_eigvec_std_update_global_scaling(beta, phi);

                        // Apply the scaling to the whole eigenvector.
                        starneig_eigvec_std_scale(ki+1, X_re, &phi);

                        // Now it is safe to execute the linear update.
                        for (int i = 0; i < j; i++)
                            X_re[i] = X_re[i]-T(i,j)*X_re[j];

                        // Recompute norm excluding the entries j:n.
                        norm = vec_real_infnorm(j, X_re);
                    }
                    // if next block is 2-by-2 diagonal block:
                    else {
                        scaling_t phi;
                        starneig_eigvec_std_init_scaling_factor(1, &phi);
                        info |= starneig_eigvec_std_solve_2x2_real_system(smin,
                            &T(j-1,j-1), ldT, lambda, &X_re[j-1], &phi);
                        starneig_eigvec_std_update_global_scaling(beta, phi);

                        // Scale remaining parts of vector.
                        starneig_eigvec_std_scale(j-1, X_re, &phi);
                        starneig_eigvec_std_scale(ki-(j+1), X_re+(j+1), &phi);
                        starneig_eigvec_std_update_norm(&norm, phi);

                        // Protect against overflow in the first linear update.
                        phi = starneig_eigvec_std_protect_update(
                            tnorm, fabs(X_re[j-1]), norm);
                        starneig_eigvec_std_update_global_scaling(beta, phi);

                        // Apply the scaling to the whole eigenvector.
                        starneig_eigvec_std_scale(ki+1, X_re, &phi);

                        // Now it is safe to execute the first linear update.
                        for (int i = 0; i < j-1; i++)
                            X_re[i] = X_re[i]-T(i,j-1)*X_re[j-1];

                        // Recompute norm excluding the entries j:n.
                        norm = vec_real_infnorm(j, X_re);

                        // Protect against overflow in the second linear update.
                        phi =
                            starneig_eigvec_std_protect_update(tnorm, fabs(X_re[j]), norm);
                        starneig_eigvec_std_update_global_scaling(beta, phi);

                        // Apply the scaling to the whole eigenvector.
                        starneig_eigvec_std_scale(ki+1, X_re, &phi);

                        // Now it is safe to execute the second linear update.
                        for (int i = 0; i < j-1; i++)
                            X_re[i] = X_re[i]-T(i,j)*X_re[j];

                        // Recompute norm excluding the entries j-1:n.
                        norm = vec_real_infnorm(j-1, X_re);

                        // We processed a 2-by-2 block, so skip the next entry.
                        j--;
                    }
                }

                // The ki-th real eigenvector has been computed. Recompute norm.
                Xnorms[si] = vec_real_infnorm(ki+1, X_re);

                // Record error flag.
                infos[si] = info;

                // This eigenvector spans 1 column. Update selected counter.
                si--;
            }
            else { // lambda_type[ki] == CMPLX
                // Locate real and imaginary part of complex eigenvector for
                // complex conjugate pair of eigenvalues (ki, ki-1).
                double *X_re = X+(si-1)*(size_t)ldX;
                double *X_im = X+si*(size_t)ldX;

                // Locate corresponding scaling factor, one per complex pair.
                scaling_t *beta = scales+si;

                // Compute eigenvalue as lambda = lambda_re + i * lambda_im.
                // A 2-by-2 block in canonical Schur form has the shape
                // [ T(ki-1,ki-1) T(ki-1,ki) ] = [ a b ]  or [ a -b ]
                // [ T(ki, ki-1)  T(ki, ki)  ]   [-c a ]     [ c  a ].
                double lambda_re = T(ki,ki);
                double lambda_im =
                    sqrt(fabs(T(ki,ki-1)))*sqrt(fabs(T(ki-1,ki)));

                // Critical threshold to detect unsafe divisions.
                const double smin = MAX(
                    DBL_EPSILON/2*(fabs(lambda_re)+fabs(lambda_im)), smlnum);

                // Form right-hand side.
                if (fabs(T(ki-1,ki)) >= fabs(T(ki,ki-1))) {
                    X_re[ki-1] = 1.0;
                    X_im[ki] = lambda_im/T(ki-1,ki);
                }
                else {
                    X_re[ki-1] = -lambda_im/T(ki,ki-1);
                    X_im[ki] = 1.0;
                }
                X_re[ki] = 0.0;
                X_im[ki-1] = 0.0;
                for (int i = 0; i < ki-1; i++) {
                    X_re[i] = -X_re[ki-1]*T(i,ki-1);
                    X_im[i] = -X_im[ki]*T(i,ki);
                }
                for (int i = ki+1; i < n; i++) {
                    X_re[i] = 0.0;
                    X_im[i] = 0.0;
                }

                // Compute norm of the entire vector.
                double norm = vec_cmplx_infnorm(ki+1, X_re, X_im);

                // Solve the upper quasi-triangular system
                // (T(0:ki-2,0:ki-2) - (lambda * I) \ (X_re + i * X_im).

                // Loop over triangular matrix above the eigenvalue. Note ki-2!
                for (int j = ki-2; j >= 0; j--) {
                    // If next block is 1-by-1 diagonal bock:
                    if (lambda_type[j] == 0) { // REAL
                        scaling_t phi;
                        starneig_eigvec_std_init_scaling_factor(1, &phi);
                        info |= starneig_eigvec_std_solve_1x1_cmplx_system(
                            smin, T(j,j), lambda_re, lambda_im,
                            X_re+j, X_im+j, &phi);
                        starneig_eigvec_std_update_global_scaling(beta, phi);

                        // Scale the remaining parts of the vector.
                        starneig_eigvec_std_scale(j, X_re, &phi);
                        starneig_eigvec_std_scale(ki-(j+1), X_re+(j+1), &phi);
                        starneig_eigvec_std_scale(j, X_im, &phi);
                        starneig_eigvec_std_scale(ki-(j+1), X_im+(j+1), &phi);
                        starneig_eigvec_std_update_norm(&norm, phi);

                        // Protect against overflow in the linear update.
                        double absmax = MAX(fabs(X_re[j]), fabs(X_im[j]));
                        phi = starneig_eigvec_std_protect_update(tnorm, absmax, norm);
                        starneig_eigvec_std_update_global_scaling(beta, phi);

                        // Apply scaling to the whole eigenvector.
                        starneig_eigvec_std_scale(ki+1, X_re, &phi);
                        starneig_eigvec_std_scale(ki+1, X_im, &phi);

                        // Now it is safe to execute the linear update.
                        for (int i = 0; i < j; i++) {
                            X_re[i] = X_re[i]-T(i,j)*X_re[j];
                            X_im[i] = X_im[i]-T(i,j)*X_im[j];
                        }

                        // Recompute norm excluding the entries j:n.
                        norm = vec_cmplx_infnorm(j, X_re, X_im);
                    }
                    // If next block is 2-by-2 diagonal block:
                    else {
                        scaling_t phi;
                        starneig_eigvec_std_init_scaling_factor(1, &phi);
                        info |= starneig_eigvec_std_solve_2x2_cmplx_system(
                            smin, &T(j-1,j-1), ldT, lambda_re, lambda_im,
                            X_re+j-1, X_im+j-1, &phi);
                        starneig_eigvec_std_update_global_scaling(beta, phi);

                        // Scale the remaining parts of the vector.
                        starneig_eigvec_std_scale(j-1, X_re, &phi);
                        starneig_eigvec_std_scale(ki-j, X_re+(j+1), &phi);
                        starneig_eigvec_std_scale(j-1, X_im, &phi);
                        starneig_eigvec_std_scale(ki-j, X_im+(j+1), &phi);
                        starneig_eigvec_std_update_norm(&norm, phi);

                        // Protect against overflow in the first linear update.
                        double absmax = MAX(fabs(X_re[j-1]), fabs(X_im[j-1]));
                        phi = starneig_eigvec_std_protect_update(tnorm, absmax, norm);
                        starneig_eigvec_std_update_global_scaling(beta, phi);

                        // Apply scaling to the whole eigenvector.
                        starneig_eigvec_std_scale(ki+1, X_re, &phi);
                        starneig_eigvec_std_scale(ki+1, X_im, &phi);

                        // Now it is safe to execute the first linear update.
                        for (int i = 0; i < j-1; i++) {
                            X_re[i] = X_re[i]-T(i,j-1)*X_re[j-1];
                            X_im[i] = X_im[i]-T(i,j-1)*X_im[j-1];
                        }

                        // Recompute norm excluding the entries j+1:n.
                        norm = vec_cmplx_infnorm(j+1, X_re, X_im);

                        // Protect against overflow in the second linear update.
                        absmax = MAX(fabs(X_re[j]), fabs(X_im[j]));
                        phi = starneig_eigvec_std_protect_update(tnorm, absmax, norm);
                        starneig_eigvec_std_update_global_scaling(beta, phi);

                        // Apply scaling to the whole eigenvector.
                        starneig_eigvec_std_scale(ki+1, X_re, &phi);
                        starneig_eigvec_std_scale(ki+1, X_im, &phi);

                        // Now it is safe to execute the second linear update.
                        for (int i = 0; i < j-1; i++) {
                            X_re[i] = X_re[i]-T(i,j)*X_re[j];
                            X_im[i] = X_im[i]-T(i,j)*X_im[j];
                        }

                        // Recompute norm excluding the entries j:n.
                        norm = vec_cmplx_infnorm(j-1, X_re, X_im);

                        // We processed a 2-by-2 block, so skip the next entry.
                        j--;
                    }
                }

                // Note that the second column of a complex conjugate pair is
                // never allocated or computed. Obtaining it is left to the
                // user. If the positions si - 1, si mark a 2-by-2 block, then
                // the eigenvector corresponding to lambda = alpha + i beta
                // is X(:, si - 1) + i * X(:, si).
                // The complex conjugate eigenvector corresponding to
                // lambda = alpha - i beta can be derived as
                // conj(X) := X(:, si -1) - i * X(:, si).

                // The real part and the imaginary part are scaled alike.
                scales[si-1] = scales[si];

                // The ki-th complex eigenvector has been computed. Recompute
                // norm of the complex eigenvector.
                Xnorms[si] = vec_cmplx_infnorm(ki+1, X_re, X_im);
                Xnorms[si-1] = Xnorms[si];

                // Record error flag.
                infos[si] = info;
                infos[si-1] = info;

                // We processed a complex conjugate pair of eigenvalues,
                // so skip the next entry.
                ki--;

                // This eigenvector spans 2 columns. Update selected counter.
                si -= 2;
            }
        }
    }

    // Note that the eigenvectors are not normalized.

#undef T

}

void starneig_eigvec_std_cpu_backsolve(void *buffers[], void *cl_args)
{
    double *T = (double *) STARPU_MATRIX_GET_PTR(buffers[0]);
    int n = STARPU_MATRIX_GET_NX(buffers[0]);
    int ldT = STARPU_MATRIX_GET_LD(buffers[0]);

    double *ubT = (double *) STARPU_VARIABLE_GET_PTR(buffers[1]);
    double tnorm = *ubT;

    double *X = (double *) STARPU_MATRIX_GET_PTR(buffers[2]);
    int num_selected = STARPU_MATRIX_GET_NY(buffers[2]);
    int ldX = STARPU_MATRIX_GET_LD(buffers[2]);

    scaling_t *scales = (scaling_t *) STARPU_VECTOR_GET_PTR(buffers[3]);

    double *Xnorms = (double *) STARPU_VECTOR_GET_PTR(buffers[4]);

    int *lambda_type = (int *) STARPU_VECTOR_GET_PTR(buffers[5]);

    int *selected = (int *) STARPU_VECTOR_GET_PTR(buffers[6]);

    int *infos = (int *) STARPU_VECTOR_GET_PTR(buffers[7]);

    double smlnum;
    starpu_codelet_unpack_args(cl_args, &smlnum);

    backsolve(n, T, ldT, tnorm, X, ldX, scales, Xnorms, lambda_type,
        selected, num_selected, smlnum, infos);
}

void starneig_eigvec_std_cpu_solve(void *buffers[], void *cl_args)
{
    double *T = (double *) STARPU_MATRIX_GET_PTR(buffers[0]);
    int n = STARPU_MATRIX_GET_NX(buffers[0]);
    int ldT = STARPU_MATRIX_GET_LD(buffers[0]);

    double *ubT = (double *) STARPU_VARIABLE_GET_PTR(buffers[1]);
    double tnorm = *ubT;

    double *X = (double *) STARPU_MATRIX_GET_PTR(buffers[2]);
    int num_selected = STARPU_MATRIX_GET_NY(buffers[2]);
    int ldX = STARPU_MATRIX_GET_LD(buffers[2]);

    scaling_t *scales = (scaling_t *) STARPU_VECTOR_GET_PTR(buffers[3]);

    double *Xnorms = (double *) STARPU_VECTOR_GET_PTR(buffers[4]);

    double *lambda = (double *) STARPU_VECTOR_GET_PTR(buffers[5]);
    int num_rhs = STARPU_VECTOR_GET_NX(buffers[5]);

    int *lambda_type = (int *) STARPU_VECTOR_GET_PTR(buffers[6]);

    int *selected = (int *) STARPU_VECTOR_GET_PTR(buffers[7]);

    int *diag_type = (int *) STARPU_VECTOR_GET_PTR(buffers[8]);

    int *infos = (int *) STARPU_VECTOR_GET_PTR(buffers[9]);

    double smlnum;
    starpu_codelet_unpack_args(cl_args, &smlnum);

#define T(i,j) T[(i) + (j) * (size_t)ldT]
    // The i-th selected eigenvalue.
    int si = num_selected-1;

    // Loop over eigenvalues.
    for (int k = num_rhs-1; k >= 0; k--) {
        if (!selected[k]) {
            // Proceed with the next eigenvalue.
            if (lambda_type[k] == 1) { // CMPLX
                // A complex conjugate pair of eigenvalues is not selected,
                // so skip the next diagonal entry.
                k--;
            }
        }
        else { // k-th eigenvalue is selected

            // Error flag.
            int info = 0;

            if (lambda_type[k] == 0) { // REAL
                // Locate eigenvector to k-th eigenvalue.
                double *X_re = X+si*(size_t)ldX;

                // Locate corresponding scaling factor.
                scaling_t *beta = scales+si;

                // Compute norm of entire vector.
                double norm = vec_real_infnorm(n, X_re);

                // Critical threshold to detect unsafe divisions.
                const double smin = MAX(DBL_EPSILON/2*fabs(lambda[k]), smlnum);

                for (int j = n-1; j >= 0; j--) {
                    // if next block is 1-by-1 diagonal block:
                    if (diag_type[j] == 0) { // REAL
                        scaling_t phi;
                        starneig_eigvec_std_init_scaling_factor(1, &phi);
                        info |= starneig_eigvec_std_solve_1x1_real_system(
                            smin, T(j,j), lambda[k], X_re+j, &phi);
                        starneig_eigvec_std_update_global_scaling(beta, phi);

                        // Scale remaining parts of vector.
                        starneig_eigvec_std_scale(j, X_re, &phi);
                        starneig_eigvec_std_scale(n-(j+1), X_re+(j+1), &phi);
                        starneig_eigvec_std_update_norm(&norm, phi);

                        // Protect against overflow in the linear update.
                        phi =
                            starneig_eigvec_std_protect_update(tnorm, fabs(X_re[j]), norm);
                        starneig_eigvec_std_update_global_scaling(beta, phi);

                        // Apply the scaling to the whole eigenvector.
                        starneig_eigvec_std_scale(n, X_re, &phi);

                        // Now it is safe to execute the linear update.
                        for (int i = 0; i < j; i++)
                            X_re[i] = X_re[i]-T(i,j)*X_re[j];

                        // Recompute norm excluding the entries j:n.
                        norm = vec_real_infnorm(j, X_re);
                    }
                    // if next block is 2-by-2 block:
                    else {
                        scaling_t phi;
                        starneig_eigvec_std_init_scaling_factor(1, &phi);
                        info |= starneig_eigvec_std_solve_2x2_real_system(smin,
                            &T(j-1,j-1), ldT, lambda[k], &X_re[j-1], &phi);
                        starneig_eigvec_std_update_global_scaling(beta, phi);

                        // Scale remaining parts of vector.
                        starneig_eigvec_std_scale(j-1, X_re, &phi);
                        starneig_eigvec_std_scale(n-(j+1), X_re+(j+1), &phi);
                        starneig_eigvec_std_update_norm(&norm, phi);

                        // Protect first linear update against overflow.
                        phi = starneig_eigvec_std_protect_update(
                            tnorm, fabs(X_re[j-1]), norm);
                        starneig_eigvec_std_update_global_scaling(beta, phi);

                        // Apply the scaling to the whole eigenvector.
                        starneig_eigvec_std_scale(n, X_re, &phi);

                        // Now it is safe to execute the linear udpate.
                        for (int i = 0; i < j-1; i++)
                            X_re[i] = X_re[i]-T(i,j-1)*X_re[j-1];

                        // Recompute norm excluding the entries j:n.
                        norm = vec_real_infnorm(j, X_re);

                        // Protect second linear update against overflow.
                        phi =
                            starneig_eigvec_std_protect_update(tnorm, fabs(X_re[j]), norm);
                        starneig_eigvec_std_update_global_scaling(beta, phi);

                        // Apply the scaling to the whole eigenvector.
                        starneig_eigvec_std_scale(n, X_re, &phi);

                        // Now it is safe to execute the linear update.
                        for (int i = 0; i < j - 1; i++)
                            X_re[i] = X_re[i]-T(i,j)*X_re[j];

                        // Recompute norm excluding the entries j-1:n.
                        norm = vec_real_infnorm(j-1, X_re);

                        // We processed a 2-by-2 block, so skip the next
                        // diagonal entry.
                        j--;
                    }
                }

                // The k-th real eigenvector has been computed.
                // Recompute norm.
                Xnorms[si] = vec_real_infnorm(n, X_re);

                // Record error flag.
                infos[si] = info;

                // This eigenvalue spans 1 column. Update selected counter.
                si--;
            }
            else { // lambda_type[k] == CMPLX
                // Locate real and imaginary part of complex eigenvector
                // for complex conjugate pair of eigenvalues (k, k-1).
                double *X_re = X+(si-1)*(size_t)ldX;
                double *X_im = X+si*(size_t)ldX;

                // The eigenvalue is lambda = lambda_re+i*lambda_im.
                double lambda_re = lambda[k-1];
                double lambda_im = fabs(lambda[k]);

                // Locate corresponding scaling factor.
                scaling_t *beta = scales+si;

                // Compute norm of entire vector.
                double norm = vec_cmplx_infnorm(n, X_re, X_im);

                // Critical threshold to detect unsafe divisions.
                const double smin = MAX(
                    DBL_EPSILON/2*(fabs(lambda_re)+fabs(lambda_im)), smlnum);

                for (int j = n-1; j >= 0; j--) {
                    // if the next block is 1-by-1 diagonal block:
                    if (diag_type[j] == 0) { // REAL
                        scaling_t phi;
                        starneig_eigvec_std_init_scaling_factor(1, &phi);
                        info |= starneig_eigvec_std_solve_1x1_cmplx_system(smin, T(j,j),
                            lambda_re, lambda_im, X_re+j, X_im+j, &phi);
                        starneig_eigvec_std_update_global_scaling(beta, phi);

                        // Scale the remaining parts of the 2 columns.
                        starneig_eigvec_std_scale(j, X_re, &phi);
                        starneig_eigvec_std_scale(n-(j+1), X_re+(j+1), &phi);
                        starneig_eigvec_std_scale(j, X_im, &phi);
                        starneig_eigvec_std_scale(n-(j+1), X_im+(j+1), &phi);
                        starneig_eigvec_std_update_norm(&norm, phi);

                        // Protect against overflow in the linear update.
                        double absmax = MAX(fabs(X_re[j]), fabs(X_im[j]));
                        phi = starneig_eigvec_std_protect_update(tnorm, absmax, norm);
                        starneig_eigvec_std_update_global_scaling(beta, phi);

                        // Apply scaling to the whole eigenvector.
                        starneig_eigvec_std_scale(n, X_re, &phi);
                        starneig_eigvec_std_scale(n, X_im, &phi);

                        // Now it is safe to execute the linear update.
                        for (int i = 0; i < j; i++) {
                            X_re[i] = X_re[i]-T(i,j)*X_re[j];
                            X_im[i] = X_im[i]-T(i,j)*X_im[j];
                        }

                        // Recompute norm excluding the entries j:n.
                        norm = vec_cmplx_infnorm(j, X_re, X_im);
                    }
                    // if next block is 2-by-2 diagonal block:
                    else {
                        scaling_t phi;
                        starneig_eigvec_std_init_scaling_factor(1, &phi);
                        info |= starneig_eigvec_std_solve_2x2_cmplx_system(smin,
                            &T(j-1,j-1), ldT, lambda_re, lambda_im,
                            X_re+j-1, X_im+j-1, &phi);
                        starneig_eigvec_std_update_global_scaling(beta, phi);

                        // Scale remaining parts of vector.
                        starneig_eigvec_std_scale(j-1, X_re, &phi);
                        starneig_eigvec_std_scale(n-(j+1), X_re+(j+1), &phi);
                        starneig_eigvec_std_scale(j-1, X_im, &phi);
                        starneig_eigvec_std_scale(n-(j+1), X_im+(j+1), &phi);
                        starneig_eigvec_std_update_norm(&norm, phi);

                        // Protect against overflow in the first linear update.
                        double absmax = MAX(fabs(X_re[j-1]), fabs(X_im[j-1]));
                        phi = starneig_eigvec_std_protect_update(tnorm, absmax, norm);
                        starneig_eigvec_std_update_global_scaling(beta, phi);

                        // Apply scaling to the whole eigenvector.
                        starneig_eigvec_std_scale(n, X_re, &phi);
                        starneig_eigvec_std_scale(n, X_im, &phi);

                        // Now it is safe to execute the first linear update.
                        for (int i = 0; i < j-1; i++) {
                            X_re[i] = X_re[i]-T(i,j-1)*X_re[j-1];
                            X_im[i] = X_im[i]-T(i,j-1)*X_im[j-1];
                        }

                        // Recompute norm excluding the entries j+1:n.
                        norm = vec_cmplx_infnorm(j+1, X_re, X_im);

                        // Protect against overflow in the second linear update.
                        absmax = MAX(fabs(X_re[j]), fabs(X_im[j]));
                        phi = starneig_eigvec_std_protect_update(tnorm, absmax, norm);
                        starneig_eigvec_std_update_global_scaling(beta, phi);

                        // Apply scaling to the whole eigenvector.
                        starneig_eigvec_std_scale(n, X_re, &phi);
                        starneig_eigvec_std_scale(n, X_im, &phi);

                        // Now it is safe to execute the second linear update.
                        for (int i = 0; i < j-1; i++) {
                            X_re[i] = X_re[i]-T(i,j)*X_re[j];
                            X_im[i] = X_im[i]-T(i,j)*X_im[j];
                        }

                        // Recompute norm excluding the entries j-1:n.
                        norm = vec_cmplx_infnorm(j-1, X_re, X_im);

                        // We processed a 2-by-2 block, so skip the next
                        // diagonal entry.
                        j--;
                    }
                }

                // Note that the second column of a complex conjugate pair is
                // never allocated or computed. Obtaining it is left to the
                // user. If the positions si - 1, si mark a 2-by-2 block, then
                // the eigenvector corresponding to lambda = alpha + i beta
                // is X(:, si - 1) + i * X(:, si).
                // The complex conjugate eigenvector corresponding to
                // lambda = alpha - i beta can be derived as
                // conj(X) := X(:, si -1) - i * X(:, si).

                // The real part and the imaginary part are scaled alike.
                scales[si-1] = scales[si];

                // The ki-th complex eigenvector has been computed.
                // Recompute norm of the complex eigenvector.
                Xnorms[si] = vec_cmplx_infnorm(n, X_re, X_im);
                Xnorms[si-1] = Xnorms[si];

                // Record error flag.
                infos[si] = info;
                infos[si-1] = info;

                // We processed a complex conjugate pair of eigenvalues,
                // so skip the next entry.
                k--;

                // This eigenvector spans 2 cols. Update selected counter.
                si -= 2;
            }
        }
    }

#undef T
}

void starneig_eigvec_std_cpu_update(void *buffers[], void *cl_args)
{
    // T is n-by-m.
    double *T = (double *) STARPU_MATRIX_GET_PTR(buffers[0]);
    int ldT = STARPU_MATRIX_GET_LD(buffers[0]);
    int n = STARPU_MATRIX_GET_NX(buffers[0]);
    int m = STARPU_MATRIX_GET_NY(buffers[0]);

    double *ubT = (double *) STARPU_VARIABLE_GET_PTR(buffers[1]);
    double tnorm = *ubT;

    double *Xin = (double *) STARPU_MATRIX_GET_PTR(buffers[2]);
    int num_rhs = STARPU_MATRIX_GET_NY(buffers[2]);
    int ldX = STARPU_MATRIX_GET_LD(buffers[2]);

    scaling_t *Xscales = (scaling_t *) STARPU_VECTOR_GET_PTR(buffers[3]);

    double *Xinnorms = (double *) STARPU_VECTOR_GET_PTR(buffers[4]);

    double *Y = (double *) STARPU_MATRIX_GET_PTR(buffers[5]);
    int ldY = STARPU_MATRIX_GET_LD(buffers[5]);

    scaling_t *Yscales = (scaling_t *) STARPU_VECTOR_GET_PTR(buffers[6]);

    double *Ynorms = (double *) STARPU_VECTOR_GET_PTR(buffers[7]);

    int *lambda_type = (int *) STARPU_VECTOR_GET_PTR(buffers[8]);

    // Pointer to X - either a copy or the original memory.
    double *X;

    // Pointer to norms of X - either a copy or the original memory.
    double *Xnorms;

    // Status flag if X has to be rescaled.
    int rescale_X = 0;

    // Status flag if Y has to be rescaled.
    int rescale_Y = 0;

    // Indicator if xinnorms has to be rescaled. Note that this only affects
    // consistency scaling and not overflow protection.
    int rescale_xnorms = 0;

    // Workspace to store locally computed scaling factors.
    scaling_t tmp_scales[num_rhs];

    //
    // Compute scaling factor.
    //

    for (int k = 0; k < num_rhs; k++) {
        if (Yscales[k] < Xscales[k]) {
            rescale_xnorms = 1;
            break;
        }
    }

    if (rescale_xnorms) {
        // As X is read-only, copy xnorms.
        Xnorms = (double *) malloc(num_rhs * sizeof(double));
        memcpy(Xnorms, Xinnorms, num_rhs * sizeof(double));

        // Simulate the consistency scaling.
        for (int k = 0; k < num_rhs; k++) {
            if (Yscales[k] < Xscales[k]) {
                // The common scaling factor is Yscales[k].
                const double s =
                    starneig_eigvec_std_compute_upscaling(Yscales[k], Xscales[k]);

                // Mark X for scaling. Physical rescaling is deferred.
                rescale_X = 1;

                // Update norm.
                Xnorms[k] = s*Xinnorms[k];
            }
            else if (Xscales[k] < Yscales[k]) {
                // The common scaling factor is Xscales[k].
                const double s =
                    starneig_eigvec_std_compute_upscaling(Xscales[k], Yscales[k]);

                // Mark Y for scaling. Physical rescaling is deferred.
                rescale_Y = 1;

                // Update norm: norm(s * Y) = s * norm(Y).
                Ynorms[k] = s*Ynorms[k];
            }
        }
    }
    else { // !rescale_xnorms
        // No changes to Xinnorms necessary. Operate on original memory.
        Xnorms = Xinnorms;

        // Xnorms does not need scaling, but Ynorms may do.
        for (int k = 0; k < num_rhs; k++) {
            if (Xscales[k] < Yscales[k]) {
                // The common scaling factor is Xscales[k].
                const double s =
                    starneig_eigvec_std_compute_upscaling(Xscales[k], Yscales[k]);

                // Mark Y for scaling. Phyiscal rescaling is deferred.
                rescale_Y = 1;

                // Update norm: norm(s * Y) = s * norm(Y).
                Ynorms[k] = s*Ynorms[k];
            }
        }
    }

    //
    // Apply scaling.
    //

    // Status flag if update is safe to execute or if rescaling is required.
    int rescale;

    // Compute scaling factors needed to survive the linear update.
    rescale = starneig_eigvec_std_protect_multi_rhs_update(
        Xnorms, num_rhs, tnorm, Ynorms, lambda_type, tmp_scales);

    if (rescale) {
        rescale_X = 1;
        rescale_Y = 1;
    }

    // If X has to be rescaled, take a copy of X and do scaling on the copy.
    if (rescale_X) {
        X = (double *) malloc((size_t)ldX * num_rhs * sizeof(double));

        for (int k = 0; k < num_rhs; k++) {
            if (Yscales[k] < Xscales[k]) {
                // Copy X and simultaneously rescale.
                const double s = starneig_eigvec_std_compute_combined_upscaling(
                    Yscales[k], Xscales[k], tmp_scales[k]);
                for (int i = 0; i < m; i++)
                    X[(size_t)ldX*k+i] = s*Xin[(size_t)ldX*k+i];
            }
            else if (Xscales[k] < Yscales[k]) {
                // Copy X and simultaneously rescale with robust update factor.
                const double s = starneig_eigvec_std_convert_scaling(tmp_scales[k]);
                for (int i = 0; i < m; i++)
                    X[(size_t)ldX*k+i] = s*Xin[(size_t)ldX*k+i];
            }
            else {
                // Xscales[k] == Yscales[k].

                // Copy X and simultaneously rescale with robust update factor.
                const double s = starneig_eigvec_std_convert_scaling(tmp_scales[k]);
                for (int i = 0; i < m; i++)
                    X[(size_t)ldX*k+i] = s*Xin[(size_t)ldX*k+i];
            }
        }
    }
    else { // !rescale_X
        // No changes to X necessary. Operate on original memory.
        X = Xin;
    }

    // If Y has to be rescaled, directly modify Y.
    if (rescale_Y) {
        for (int k = 0; k < num_rhs; k++) {
            if (Yscales[k] < Xscales[k]) {
                // The common scaling factor is Yscales[k]. Rescale Y with
                // robust update factor, if necessary.
                starneig_eigvec_std_scale(n, Y+(size_t)ldY*k, tmp_scales+k);
            }
            else if (Xscales[k] < Yscales[k]) {
                // The common scaling factor is Xscales[k]. Combine with
                // robust update scaling factor.
                const double s = starneig_eigvec_std_compute_combined_upscaling(
                    Xscales[k], Yscales[k], tmp_scales[k]);
                for (int i = 0; i < n; i++)
                    Y[(size_t)ldY*k+i] = s*Y[(size_t)ldY*k+i];
            }
            else {
                // Xscales[k] == Yscales[k].

                // Rescale Y with robust update factor, if necessary.
                starneig_eigvec_std_scale(n, Y+(size_t)ldY*k, tmp_scales+k);
            }
        }
    }

    //
    // Update global scaling of Y.
    //

#ifdef STARNEIG_ENABLE_INTEGER_SCALING
    for (int k = 0; k < num_rhs; k++)
        Yscales[k] = MIN(Yscales[k], Xscales[k])+tmp_scales[k];
#else
    for (int k = 0; k < num_rhs; k++)
        Yscales[k] = MIN(Yscales[k], Xscales[k])*tmp_scales[k];
#endif

    //
    // Compute update.
    //

    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
        n, num_rhs, m, -1.0, T, ldT, X, ldX, 1.0, Y, ldY);

    //
    // Recompute norms.
    //

    for (int k = 0; k < num_rhs; k++) {
        if (lambda_type[k] == 1) { // CMPLX
            // We store only one scaling factor per complex eigenvector pair.
            // So interpret columns as real and imaginary part.
            const double *Y_re = Y+(size_t)k*ldY;
            const double *Y_im = Y+(size_t)(k+1)*ldY;

            Ynorms[k] = vec_cmplx_infnorm(n, Y_re, Y_im);

            // Duplicate norm for real and imaginary column.
            Ynorms[k+1] = Ynorms[k];

            k++;
        }
        else { // lambda_type[k] == REAL
            Ynorms[k] = vec_real_infnorm(n, Y+(size_t)k*ldY);
        }
    }

    //
    // Clean up.
    //

    if (rescale_X)
        free(X);

    if (rescale_xnorms)
        free(Xnorms);
}

void starneig_eigvec_std_cpu_backtransform(void *buffers[], void *cl_args)
{
    double *Q = (double *) STARPU_MATRIX_GET_PTR(buffers[0]);
    int ldQ = STARPU_MATRIX_GET_LD(buffers[0]);

    double *X = (double *) STARPU_MATRIX_GET_PTR(buffers[1]);
    int ldX = STARPU_MATRIX_GET_LD(buffers[1]);

    double *Y = (double *) STARPU_MATRIX_GET_PTR(buffers[2]);
    int ldY = STARPU_MATRIX_GET_LD(buffers[2]);
    int m = STARPU_MATRIX_GET_NX(buffers[2]);
    int n = STARPU_MATRIX_GET_NY(buffers[2]);

    int k;
    starpu_codelet_unpack_args(cl_args, &k);

    //   Yij  :=   Qi:  *   X:j
    // (m x n)   (m x k)  (k x n)

    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
        m, n, k, 1.0, Q, ldQ, X, ldX, 0.0, Y, ldY);

}

ext: c

max_stars_repo_name: NLAFET/StarNEig

max_issues_repo_name: NLAFET/StarNEig

max_forks_repo_name: NLAFET/StarNEig


text:
#if !defined(ISPLAINASCII_H_INCLUDED)
#define ISPLAINASCII_H_INCLUDED

#include "FileEnumerator.h"
#include "Utils.h"

#include <filesystem>
#include <gsl/gsl>
#include <iosfwd>
#include <string>

class IsPlainAscii
{
public:
	static int usage(::std::ostream& strm, const ::std::string& progName,
		const char* pMsg);

	IsPlainAscii(::gsl::span<const char*const> args);
	int run() const;

	IsPlainAscii(const IsPlainAscii&) = delete;
	IsPlainAscii& operator=(const IsPlainAscii&) = delete;
	IsPlainAscii(IsPlainAscii&&) = delete;
	IsPlainAscii& operator=(IsPlainAscii&&) = delete;

PRIVATE_EXCEPT_IN_TEST:
	using Path = ::std::filesystem::path;

	static void scanFile(const Path& filePath);
	static void scanFile2(const Path& filePath, ::std::istream& in,
		::std::ostream& out);
	static void reportNonAsciiRun(const Path& filePath, unsigned long lineNum,
		unsigned long approxColNum, ::std::string& nonAsciiRun, ::std::ostream& out);

	FileEnumerator m_fileEnumerator;
};

#endif // ISPLAINASCII_H_INCLUDED

max_stars_repo_name: IanEmmons/CmdLineUtil

max_issues_repo_name: IanEmmons/CmdLineUtil

max_forks_repo_name: IanEmmons/CmdLineUtil