text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <getopt.h>
#include <mpi.h>
#include <iniparser.h>
#include "prepmt/prepmt_hpulse96.h"
#include "prepmt/prepmt_event.h"
#ifdef PARMT_USE_INTEL
#include <mkl_cblas.h>
#else
#include <cblas.h>
#endif
#include "sacio.h"
#include "cps.h"
#include "cps_mpi.h"
#include "cps_utils.h"
#include "cps_defaults.h"
#include "iscl/array/array.h"
#include "iscl/memory/memory.h"
#include "iscl/os/os.h"
#include "iscl/string/string.h"

#define PROGRAM_NAME "hspec96"

static void grns2sac(const struct hwave_greens_struct grns,
                     struct sacData_struct sac[10]);
struct precompute_struct
{
    char sourceModel[PATH_MAX]; /*!< Source model */
    char crustDir[PATH_MAX];  /*!< Crust1.0 directory */
    char hspecDistanceFile[PATH_MAX];  /*!< Table of distances for hspec */
    char hspecArchiveFile[PATH_MAX];  /*!< Name of archive file with ff green's
                                           functions. */
    double *depths; /*!< Depths (km) at which to compute greens fns */
    double *dists;  /*!< Distances (degrees) at which to compute
                         surface wave greens fns */
    double slat;    /*!< Source latitude (degrees) */
    double slon;    /*!< Source longitude (degrees) */
    double dt;   /*!< Sampling period */
    int npts;    /*!< Number of samples */
    int ndepths; /*!< Number of depths at which to compute greens fns */
    int ndists;  /*!< Number of distances to compute surface waves */
    bool luseSourceMod; /*!< If true then use the local source model */
    bool luseCrust; /*!< If true then use crust1.0 at teleseismic distances */
    bool luseDistanceTable; /*!< If true then use a distance table.
                                 Otherwise linearly interpolate distances */
};
int prepmt_hspec96_readParameters(const char *iniFile,
                                  const char *section,
                                  struct precompute_struct *precompute);
int prepmt_hspec96_initializeArchive(
    const char *archiveName, const char *model,
    const int ndepths, const double *__restrict__ depths,
    const int ngcarcs,  const double *__restrict__ gcarcs);
int prepmt_hspec96_readHspec96Parameters(
    const char *iniFile,
    struct hprep96_parms_struct *hprepParms,
    struct hspec96_parms_struct *hspecParms,
    struct hpulse96_parms_struct *hpulseParms);
static void printUsage(void);
static int parseArguments(int argc, char *argv[], char iniFile[PATH_MAX]);

/*!
 * @brief Makes the regional fullwaveform fundamental fault Green's
 *        functions archive.
 *
 */
int main(int argc, char **argv)
{
    char fname[PATH_MAX], iniFile[PATH_MAX], groupName[512], *modelName;
    struct prepmtEventParms_struct event;
    struct sacData_struct *sac;
    struct precompute_struct precompute;
    struct hwave_greens_struct *grns;
    struct hpulse96_data_struct *zresp;
    struct hprep96_parms_struct hprepParms;
    struct hspec96_parms_struct hspecParms;
    struct hpulse96_parms_struct hpulseParms;
    struct vmodel_struct recmod, vmodel;
    double *depthr, *depths, *r, *tshift, *vred, dt;
    hid_t groupID, h5fl;
    int idep, idist, ierr, k, l, myid, ndeps, ndists, nprocs, npts, nwaves;
    // Here's the idea: I want to input magnitudes with units of N-m so:
    // N-m -> Dyne-dm (1.e+7)
    // CPS internally scales from Dyne-cm to cm (1.e+20)
    // Finally, I want outputs proprotional to m (1.e-2)
    const double xmom = 1.0;     // no confusing `relative' magnitudes
    const double xcps = 1.e-20;  // convert dyne-cm mt to output cm
    const double cm2m = 1.e-2;   // cm to meters
    const double dcm2nm = 1.e+7; // magnitudes intended to be specified in
                                 // Dyne-cm but I work in N-m
    const char *cfaults[10] = {"ZDS\0", "ZSS\0", "ZDD", "ZEX\0",
                               "RDS\0", "RSS\0", "RDD", "REX\0",
                               "TDS\0", "TSS\0"};
    // Given a M0 in Newton-meters get a seismogram in meters
    const double xscal = xmom*xcps*cm2m*dcm2nm;
    const int master = 0;
    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &myid);
    MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
    sac = NULL;
    grns = NULL;
    zresp = NULL;
    modelName = NULL;
    r = NULL;
    depthr = NULL;
    depths = NULL;
    vred = NULL;
    tshift = NULL;
    memset(&hprepParms, 0, sizeof(struct hprep96_parms_struct));
    memset(&hspecParms, 0, sizeof(struct hspec96_parms_struct));
    memset(&hpulseParms, 0, sizeof(struct hpulse96_parms_struct));
    memset(&vmodel, 0, sizeof(struct vmodel_struct));
    if (myid == master)
    {
        // Parse the input arguments
        ierr = parseArguments(argc, argv, iniFile);
        if (ierr != 0)
        {
            if (ierr !=-1)
            {
                printf("%s: Failed to parse input commands\n", PROGRAM_NAME);
            }
            goto FINISH_EARLY;
        }
        // Read the CPS variables
        ierr = prepmt_hspec96_readHspec96Parameters(iniFile,
                                                    &hprepParms, &hspecParms,
                                                    &hpulseParms);
        if (ierr != 0)
        {
            printf("Error reading hprep parameters\n");
            MPI_Abort(MPI_COMM_WORLD, 30);
        }
        // Read the additional contrived forward modeling variables
        ierr = prepmt_hspec96_readParameters(iniFile,
                                             "precompute", &precompute);
        if (ierr != 0)
        {
            printf("Error reading precompute parameters\n");
            MPI_Abort(MPI_COMM_WORLD, 30);
        }
        ndists = precompute.ndists;
        ndeps = precompute.ndepths;
        // expand to make a grid
        r = (double *) calloc((size_t) (ndists*ndeps), sizeof(double));
        depths = (double *) calloc((size_t) (ndists*ndeps), sizeof(double));
        for (idep=0; idep<ndeps; idep++)
        {
            for (idist=0; idist<ndists; idist++)
            {
                k = idist*ndeps + idep;
                r[k] = precompute.dists[idist]*111.195;
                depths[k] = precompute.depths[idep];
                printf("Wavform %d has distance %f (deg) and depth %f (km)\n",
                       k+1, r[k]/111.195, depths[k]);
            }
        }
        npts = precompute.npts;
        dt = precompute.dt;
//ndists = 31;
//ndeps = 26;
//npts = 1024; //2048;
//dt = 1.0; //0.5;
        printf("%d %f %s\n", npts, dt, hprepParms.mfile);
        // load the model into memory
        if (precompute.luseSourceMod)
        {
            printf("%s: Loading source model %s...\n",
                   PROGRAM_NAME, precompute.sourceModel);
            strcpy(hprepParms.mfile, precompute.sourceModel);
            ierr = cps_getmod(hprepParms.mfile, &vmodel);
            if (ierr != 0)
            {
                printf("Failed to load source model\n");
                goto FINISH_EARLY;
            }
        }
        else if (precompute.luseCrust)
        {
            ierr = prepmt_event_initializeFromIniFile(iniFile, &event);
            if (ierr != 0)
            {
                printf("%s: Error getting event info\n", PROGRAM_NAME);
                goto FINISH_EARLY;
            }
            printf("%s: Loading CPS model (lat,lon)=(%f,%f)\n", PROGRAM_NAME,
                   event.latitude, event.longitude);
            ierr = cps_crust1_getCrust1ForHerrmann(precompute.crustDir, false,
                                              event.latitude, event.longitude,
                                              &event.latitude, &event.longitude,
                                              1, &vmodel, &recmod);
            cps_utils_freeVmodelStruct(&recmod);
        }
        else
        {
            printf("%s: Model is not set\n", PROGRAM_NAME);
            ierr = 1;
            goto FINISH_EARLY;
        }
        // Initialize the archive - remove NULL terminator
        modelName = string_strip(NULL, vmodel.title);
        printf("%s: Initializing archive file %s with model %s\n",
               PROGRAM_NAME, precompute.hspecArchiveFile, modelName);
        ierr = prepmt_hspec96_initializeArchive(precompute.hspecArchiveFile,
                                                modelName,
                                                ndeps, depths,
                                                ndists, precompute.dists);
        if (ierr != 0)
        {
            printf("%s: Failed to initialize archive\n", PROGRAM_NAME);
            goto FINISH_EARLY;
        }
        free(precompute.depths);
        free(precompute.dists);
    }
FINISH_EARLY:;
    MPI_Barrier(MPI_COMM_WORLD);
    MPI_Bcast(&ierr, 1, MPI_INTEGER, master, MPI_COMM_WORLD);
    if (ierr != 0){goto FINISH;}
    cps_broadcast_vmodelStruct(MPI_COMM_WORLD, master, &vmodel);
    cps_broadcast_hprep96ParmsStruct(MPI_COMM_WORLD, master, &hprepParms);
    cps_broadcast_hspec96ParmsStruct(MPI_COMM_WORLD, master, &hspecParms);
    cps_broadcast_hpulse96ParmsStruct(MPI_COMM_WORLD, master, &hpulseParms);
    MPI_Bcast(&ndists, 1, MPI_INTEGER, master, MPI_COMM_WORLD);
    MPI_Bcast(&ndeps,  1, MPI_INTEGER, master, MPI_COMM_WORLD);
    MPI_Bcast(&npts, 1, MPI_INTEGER, master, MPI_COMM_WORLD);
    MPI_Bcast(&dt, 1, MPI_DOUBLE, master, MPI_COMM_WORLD);
    // initialize receiver information
    nwaves = ndists*ndeps;
    depthr = (double *) calloc((size_t) nwaves, sizeof(double));
    if (myid != master)
    {
        depths = (double *) calloc((size_t) nwaves, sizeof(double));
        r = (double *) calloc((size_t) nwaves, sizeof(double));
    }
    tshift = (double *) calloc((size_t) nwaves, sizeof(double));
    vred = (double *) calloc((size_t) nwaves, sizeof(double));
    MPI_Bcast(depths, nwaves, MPI_DOUBLE, master, MPI_COMM_WORLD);
    MPI_Bcast(r, nwaves, MPI_DOUBLE, master, MPI_COMM_WORLD);
    // set the receiver distances
//r[0] = 55.597; //111.195;
//r[1] = 111.195;
//r[2] = 166.792;
//depths[0] = 7.0;
//for (idep=0; idep<ndeps; idep++)
//{
//    for (idist=0;  idist<ndists; idist++)
//    {
//        k = idist*ndeps + idep;
//        r[k] = 55.597*(double) (idist + 1);
//        depths[k] = (double) idep;
//    }
//}
/*
printf("premature end\n");
MPI_Finalize();
return 0;
*/
    // call hspec
    zresp = hspec96_mpi_interface(MPI_COMM_WORLD, nwaves, npts, dt,
                                  r, tshift, vred, depthr, depths,
                                  hprepParms, hspecParms,
                                  vmodel, &ierr);
    MPI_Barrier(MPI_COMM_WORLD);
    if (ierr != 0)
    {
        printf("Error calling hspec\n");
        goto FINISH;
    }
    // convolve STF and convert back to time domain
    grns = hpulse96_mpi_interface(MPI_COMM_WORLD, nwaves, hpulseParms,
                                  zresp, &ierr);
    if (ierr != 0)
    {
        printf("Error computing greens fns\n");
        goto FINISH;
    }
    // dump the results
    if (myid == master)
    {
        printf("Writing results...\n");
        h5fl = H5Fopen(precompute.hspecArchiveFile, H5F_ACC_RDWR, H5P_DEFAULT);
        sac = (struct sacData_struct *)
              calloc(10, sizeof(struct sacData_struct));
        k = 0;
        //for (idist=0; idist<ndists; idist++)
        for (idep=0; idep<ndeps; idep++)
        {
            memset(fname, 0, PATH_MAX*sizeof(char));
            sprintf(fname, "dump/depth_%d", idep);
            os_makedirs(fname);

            memset(groupName, 0, 512*sizeof(char));
            sprintf(groupName, "/%s/Depth_%d", modelName, idep);
            groupID = H5Gopen(h5fl, groupName, H5P_DEFAULT);
            //for (idep=0; idep<ndeps; idep++)
            for (idist=0; idist<ndists; idist++)
            {
                k = idist*ndeps + idep;
                grns2sac(grns[k], sac);
                memset(fname, 0, PATH_MAX*sizeof(char));
                sprintf(fname, "dump/depth_%d/%04d01ZDD.SAC", idep, idist+1);
                sacio_writeTimeSeriesFile(fname, sac[0]);
                sprintf(fname, "dump/depth_%d/%04d02RDD.SAC", idep, idist+1);
                sacio_writeTimeSeriesFile(fname, sac[1]);
                sprintf(fname, "dump/depth_%d/%04d03ZDS.SAC", idep, idist+1);
                sacio_writeTimeSeriesFile(fname, sac[2]);
                sprintf(fname, "dump/depth_%d/%04d04RDS.SAC", idep, idist+1);
                sacio_writeTimeSeriesFile(fname, sac[3]);
                sprintf(fname, "dump/depth_%d/%04d05TDS.SAC", idep, idist+1);
                sacio_writeTimeSeriesFile(fname, sac[4]);
                sprintf(fname, "dump/depth_%d/%04d06ZSS.SAC", idep, idist+1);
                sacio_writeTimeSeriesFile(fname, sac[5]);
                sprintf(fname, "dump/depth_%d/%04d07RSS.SAC", idep, idist+1);
                sacio_writeTimeSeriesFile(fname, sac[6]);
                sprintf(fname, "dump/depth_%d/%04d08TSS.SAC", idep, idist+1);
                sacio_writeTimeSeriesFile(fname, sac[7]);
                sprintf(fname, "dump/depth_%d/%04d09ZEX.SAC", idep, idist+1);
                sacio_writeTimeSeriesFile(fname, sac[8]);
                sprintf(fname, "dump/depth_%d/%04d10REX.SAC", idep, idist+1);
                sacio_writeTimeSeriesFile(fname, sac[9]);
                // Dump to the H5 archive and release memory
                for (l=0; l<10; l++)
                {
                    // Rescale to unit magnitude
                    cblas_dscal(sac[l].npts, xscal, sac[l].data, 1);
                    sprintf(fname, "Greens_%s_%d", cfaults[l], idist);
                    sacioh5_writeTimeSeries2(fname, groupID, sac[l]);
                    sacio_free(&sac[l]);
                }
//break;
            }
            H5Gclose(groupID);
//break;
        }
        if (sac != NULL){free(sac);}
        H5Fclose(h5fl);
        printf("%s: Finishing program\n", PROGRAM_NAME);
    }
    MPI_Barrier(MPI_COMM_WORLD);
FINISH:;
    if (grns != NULL)
    {
        for (k=0; k<nwaves; k++)
        {
            cps_utils_freeHwaveGreensStruct(&grns[k]);
        }
        free(grns);
        grns = NULL;
    }
    if (zresp != NULL)
    {
        for (k=0; k<nwaves; k++)
        {
            cps_utils_freeHpulse96DataStruct(&zresp[k]);
        }
        free(zresp);
        zresp = NULL;
    }
    // release memory
    memory_free8c(&modelName);
    free(depthr);
    free(depths);
    free(r);
    free(tshift);
    free(vred);
    cps_utils_freeHpulse96ParmsStruct(&hpulseParms);
    cps_utils_freeVmodelStruct(&vmodel);
    MPI_Finalize();
    return 0;
}
//============================================================================//
int prepmt_hspec96_initializeArchive(
    const char *archiveName, const char *model,
    const int ndepths, const double *__restrict__ depths,
    const int ngcarcs,  const double *__restrict__ gcarcs)
{
    char modelName[512], depthName[512]; //, varName[512];
    hid_t h5fl, dataSet, dataSpace, depthGroup, modelGroup;
    hsize_t dims[1];
    int idep;
    //------------------------------------------------------------------------//
    h5fl = H5Fcreate(archiveName, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
    memset(modelName, 0, 512*sizeof(char));
    sprintf(modelName, "/%s", model);
    printf("%s: Initializing model: %s\n", __func__, modelName);
    modelGroup = H5Gcreate2(h5fl, modelName, H5P_DEFAULT,
                            H5P_DEFAULT, H5P_DEFAULT);
    // write the depths and distances
    dims[0] = (hsize_t) ndepths;
    dataSpace = H5Screate_simple(1, dims, NULL);
    dataSet = H5Dcreate2(modelGroup, "Depths\0", H5T_NATIVE_DOUBLE,
                         dataSpace,
                         H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
    H5Dwrite(dataSet, H5T_NATIVE_DOUBLE, H5S_ALL,
             H5S_ALL, H5P_DEFAULT, depths);
    H5Dclose(dataSet);
    H5Sclose(dataSpace);
    // Write the distances
    dims[0] = (hsize_t) ngcarcs;
    dataSpace = H5Screate_simple(1, dims, NULL);
    dataSet = H5Dcreate2(modelGroup, "Distances\0", H5T_NATIVE_DOUBLE,
                         dataSpace,
                         H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
    H5Dwrite(dataSet, H5T_NATIVE_DOUBLE, H5S_ALL,
             H5S_ALL, H5P_DEFAULT, gcarcs);
    H5Dclose(dataSet);
    H5Sclose(dataSpace);
    // Create the depths groups
    for (idep=0; idep<ndepths; idep++)
    {
        memset(depthName, 0, 512*sizeof(char));
        sprintf(depthName, "%s/Depth_%d", modelName, idep);
        depthGroup = H5Gcreate2(h5fl, depthName, H5P_DEFAULT,
                                H5P_DEFAULT, H5P_DEFAULT);
        H5Gclose(depthGroup);
    }
    H5Gclose(modelGroup);
    H5Fclose(h5fl);
    return 0;
}
//============================================================================//
int prepmt_hspec96_readParameters(const char *iniFile,
                                  const char *section,
                                  struct precompute_struct *precompute)
{
    const char *fcnm = "prepmt_hspec96_readParameters\0";
    FILE *dfl;
    const char *s;
    char *dirName, cline[256], vname[128];
    double depth0, depth1, dist0, dist1;
    int i, ierr;
    dictionary *ini;
    memset(precompute, 0, sizeof(struct precompute_struct));
    if (!os_path_isfile(iniFile))
    {
        printf("%s: ini file: %s does not exist\n", fcnm, iniFile);
        return -1;
    }
    ini = iniparser_load(iniFile);
    if (ini == NULL)
    {
        printf("%s: Cannot parse ini file\n", fcnm);
        return -1;
    }
    memset(vname, 0, 128*sizeof(char));
    sprintf(vname, "%s:ndepths", section);
    precompute->ndepths = iniparser_getint(ini, vname, 1);

    memset(vname, 0, 128*sizeof(char));
    sprintf(vname, "%s:depthMin", section);
    depth0 = iniparser_getdouble(ini, vname, 1.0);

    memset(vname, 0, 128*sizeof(char));
    sprintf(vname, "%s:depthMax", section);
    depth1 = iniparser_getdouble(ini, vname, 1.0);

    memset(vname, 0, 128*sizeof(char));
    sprintf(vname, "%s:ndists", section);
    precompute->ndists = iniparser_getint(ini, vname, 1);

    memset(vname, 0, 128*sizeof(char));
    sprintf(vname, "%s:minDistance", section);
    dist0 = iniparser_getdouble(ini, vname, 1.0);

    memset(vname, 0, 128*sizeof(char));
    sprintf(vname, "%s:maxDistance", section);
    dist1 = iniparser_getdouble(ini, vname, 1.0);

    memset(vname, 0, 128*sizeof(char));
    sprintf(vname, "%s:dtHspec", section);
    precompute->dt = iniparser_getdouble(ini, vname, 1.0);

    memset(vname, 0, 128*sizeof(char));
    sprintf(vname, "%s:nptsHspec", section);
    precompute->npts = iniparser_getint(ini, vname, 1024);

    memset(vname, 0, 128*sizeof(char));
    sprintf(vname, "%s:luseCrust1", section);
    precompute->luseCrust = iniparser_getboolean(ini, vname, 0);

    memset(vname, 0, 128*sizeof(char));
    sprintf(vname, "%s:luseSourceModel", section);
    precompute->luseSourceMod = iniparser_getboolean(ini, vname, 1);

    memset(vname, 0, 128*sizeof(char));
    sprintf(vname, "%s:sourceModel", section);
    s = iniparser_getstring(ini, vname, NULL);
    if (s != NULL)
    {
        strcpy(precompute->sourceModel, s);
        if (!os_path_isfile(precompute->sourceModel) &&
            precompute->luseSourceMod)
        {
            printf("%s: Source model doesn't %s exist\n", fcnm, s);
            precompute->luseSourceMod = false;
        }
    }
    else
    {
        precompute->luseSourceMod = false;
    }

    memset(vname, 0, 128*sizeof(char));
    sprintf(vname, "%s:hspecArchiveFile", section);
    s = iniparser_getstring(ini, vname, "hspecFFgreens.h5");
    dirName = os_dirname(s, &ierr);
    if (!os_path_isdir(dirName))
    {
        ierr = os_makedirs(dirName);
        if (ierr != 0)
        {
            printf("%s: Failed to make output directory: %s\n", fcnm, dirName);
            return -1;
        }
    }
    memory_free8c(&dirName);
    strcpy(precompute->hspecArchiveFile, s);

    //memset(vname, 0, 128*sizeof(char));
    //sprintf(vname, "%s:modelIdentifier", section);
    //s = iniparser_getstring(ini, vname, "UNKNOWN");
    //strcpy(precompute->modelIdentifer, s);

    memset(vname, 0, 128*sizeof(char));
    sprintf(vname, "%s:crustDir", section);
    s = iniparser_getstring(ini, vname, CPS_DEFAULT_CRUST1_DIRECTORY);
    if (s != NULL)
    {
        strcpy(precompute->crustDir, s);
        if (!os_path_isdir(precompute->crustDir) && precompute->luseCrust)
        {
            printf("%s: crust1.0 directory %s doesn't exist\n", fcnm, s);
            precompute->luseCrust = false;
        }
    }
    else
    {
        precompute->luseCrust = false;
    }

    memset(vname, 0, 128*sizeof(char));
    sprintf(vname, "%s:luseDistanceTable", section);
    precompute->luseDistanceTable = iniparser_getboolean(ini, vname, false);
    if (precompute->luseDistanceTable)
    {
        memset(vname, 0, 128*sizeof(char));
        sprintf(vname, "%s:hspecDistanceFile", section);
        s = iniparser_getstring(ini, vname, NULL);
        if (!os_path_isfile(s))
        {
            printf("%s: files %s doesn't exist\n", fcnm, s);
            return -1;
        }
        strcpy(precompute->hspecDistanceFile, s);
        dfl = fopen(s, "r");
        memset(cline, 0, 256*sizeof(char));
        fgets(cline, 256, dfl);
        memset(cline, 0, 256*sizeof(char));
        fgets(cline, 256, dfl);
        precompute->ndists = atoi(cline);
        if (precompute->ndists < 1)
        {
            printf("%s: Invalid number of distances\n", fcnm);
            return -1;
        }
        precompute->dists
            = (double *) calloc((size_t) precompute->ndists, sizeof(double));
        for (i=0; i<precompute->ndists; i++)
        {
            memset(cline, 0, 256*sizeof(char));
            if (fgets(cline, 256, dfl) == NULL)
            {
                printf("%s: premature end of distance file\n", fcnm);
                return -1;
            }
            precompute->dists[i] = atof(cline);
        }
        fclose(dfl);
    }
    else
    {
        precompute->dists
            = (double *) calloc((size_t) precompute->ndists, sizeof(double));
        ierr = array_linspace64f_work(dist0, dist1,
                                      precompute->ndists, precompute->dists);
    }
    // depths
    precompute->depths
        = (double *) calloc((size_t) precompute->ndepths, sizeof(double));
    ierr = array_linspace64f_work(depth0, depth1,
                                  precompute->ndepths, precompute->depths);
    iniparser_freedict(ini);
    return ierr;
}
//============================================================================//
int prepmt_hspec96_readHspec96Parameters(
    const char *iniFile,
    struct hprep96_parms_struct *hprepParms,
    struct hspec96_parms_struct *hspecParms,
    struct hpulse96_parms_struct *hpulseParms)
{
    const char *fcnm = "prepmt_hspec96_readHspec96Parameters\0";
    int ierr;
    cps_setHprep96Defaults(hprepParms);
    cps_setHspec96Defaults(hspecParms);
    ierr = prepmt_hpulse96_readHpulse96Parameters(iniFile,
                                                  "hpulse96\0", hpulseParms);
    if (ierr != 0)
    {
        printf("%s: Error reading hpulse parameters\n", fcnm);
    }
    return ierr;
}
//============================================================================//
/*!
 * @brief Parses the command line arguments for the ini file.
 *
 * @param[in] argc      Number of arguments input to the program.
 * @param[in] argv      Command line arguments.
 *
 * @param[out] iniFile  If the result is 0 then this is the ini file.
 *
 * @result 0 indicates success. \n
 *        -1 indicates the user inquired about the usage. \n
 *        -2 indicates the command line argument is invalid.
 *
 * @author Ben Baker, ISTI
 *
 */
static int parseArguments(int argc, char *argv[], char iniFile[PATH_MAX])
{
    bool linFile;
    linFile = false;
    memset(iniFile, 0, PATH_MAX*sizeof(char));
    while (true)
    {
        static struct option longOptions[] =
        {
            {"help", no_argument, 0, '?'},
            {"help", no_argument, 0, 'h'},
            {"ini_file", required_argument, 0, 'i'},
            {"section", required_argument, 0, 's'},
            {0, 0, 0, 0}
        };
        int c, optionIndex;
        c = getopt_long(argc, argv, "?hi:",
                        longOptions, &optionIndex);
        if (c ==-1){break;}
        if (c == 'i')
        {
            strcpy(iniFile, (const char *) optarg);
            linFile = true;
        }
        else if (c == 'h' || c == '?')
        {
            printUsage();
            return -2;
        }
        else
        {
            printf("%s: Unknown options: %s\n",
                   PROGRAM_NAME, argv[optionIndex]);
        }
    }
    if (!linFile)
    {
        printf("%s: Error must specify ini file\n\n", PROGRAM_NAME);
        printUsage();
        return -1;
    }
    else
    {
        if (!os_path_isfile(iniFile))
        {
            printf("%s: Error - ini file: %s does not exist\n",
                   PROGRAM_NAME, iniFile);
            return EXIT_FAILURE;
        }
    }
    return 0;
}
//============================================================================//
static void printUsage(void)
{
    printf("Usage:\n   %s -i iniFile\n", PROGRAM_NAME);
    printf("or\n");
    printf("    mpirun -np nProcessors %s -i iniFile\n", PROGRAM_NAME);
    printf("Required arguments:\n");
    printf("    -i ini_file specifies the initialization file\n");
    printf("Optional arguments:\n");
    printf("    -h displays this message\n");
    printf("    -np Number of computer processors to use\n");
    return;
}
//============================================================================//
static void grns2sac(const struct hwave_greens_struct grns,
                     struct sacData_struct sac[10])
{
    const char *kcmpnm[10] = {"ZDD\0", "RDD\0", "ZDS\0", "RDS\0", "TDS\0",
                              "ZSS\0", "RSS\0", "TSS\0", "ZEX\0", "REX\0"};
    double cmpinc, cmpaz;
    int i, k;
    bool lsh;
    for (k=0; k<10; k++)
    {
        memset(&sac[k], 0, sizeof(struct sacData_struct));
        sacio_setDefaultHeader(&sac[k].header);
        lsh = false;
        sac[k].header.npts = grns.npts;
        //sac[k].header.b = grns.b;
        sac[k].npts = grns.npts;
        sac[k].header.b = grns.t0;
        sac[k].header.e = grns.t0 + (double) (grns.npts - 1)*grns.dt;
        sac[k].header.delta = grns.dt;
        sac[k].header.dist = grns.dist;
        sac[k].header.evdp = grns.evdep;
        sac[k].header.stel = grns.stelel;
        sac[k].header.gcarc = grns.dist/111.195;
        sac[k].header.o =-grns.t0;
        strcpy(sac[k].header.kevnm, "SYNTHETIC");
        strcpy(sac[k].header.ko, "O\0");
        sac[k].header.a = grns.timep;
        strcpy(sac[k].header.ka, "P\0");
        sac[k].header.az = 0.0;
        sac[k].header.baz = 180.0;
        sac[k].header.lcalda = 0;
        sac[k].header.norid = 0;
        sac[k].header.nevid = 0;
        sac[k].header.nzyear = 1970;
        sac[k].header.nzjday = 1;
        sac[k].header.nzhour = 0;
        sac[k].header.nzmin = 0;
        sac[k].header.nzsec = 0;
        sac[k].header.nzmsec = 0;
        sac[k].header.lhaveHeader = true;
        strcpy(sac[k].header.kcmpnm, kcmpnm[k]);
        sac[k].data = (double *) calloc((size_t) grns.npts, sizeof(double));
        cmpaz = 0.0;
        cmpinc = 0.0;
        if (k == 0)
        {
            for (i=0; i<grns.npts; i++){sac[k].data[i] = grns.zdd[i];}
        }
        else if (k == 1)
        {
            for (i=0; i<grns.npts; i++){sac[k].data[i] = grns.rdd[i];}
            cmpinc = 90.0;
        }
        else if (k == 2)
        {
            for (i=0; i<grns.npts; i++){sac[k].data[i] = grns.zds[i];}
        }
        else if (k == 3)
        {
            for (i=0; i<grns.npts; i++){sac[k].data[i] = grns.rds[i];}
            cmpinc = 90.0;
        }
        else if (k == 4)
        {
            for (i=0; i<grns.npts; i++){sac[k].data[i] = grns.tds[i];}
            lsh = true;
            cmpaz = 90.0;
            cmpinc = 90.0;
        }
        else if (k == 5)
        {
            for (i=0; i<grns.npts; i++){sac[k].data[i] = grns.zss[i];}
        }
        else if (k == 6)
        {
            for (i=0; i<grns.npts; i++){sac[k].data[i] = grns.rss[i];}
            cmpinc = 90.0;
        }
        else if (k == 7)
        {
            for (i=0; i<grns.npts; i++){sac[k].data[i] = grns.tss[i];}
            lsh = true;
            cmpaz = 90.0;
            cmpinc = 90.0;
        }
        else if (k == 8)
        {
            for (i=0; i<grns.npts; i++){sac[k].data[i] = grns.zex[i];}
        }
        else if (k == 9)
        {
            for (i=0; i<grns.npts; i++){sac[k].data[i] = grns.rex[i];}
            cmpinc = 90.0;
        }
        if (lsh)
        {
            sac[k].header.t1 = grns.timesh;
            strcpy(sac[k].header.kt0, "SH\0");
        }
        else
        {
            sac[k].header.t1 = grns.timesv;
            strcpy(sac[k].header.kt0, "SV\0");
        }
        sac[k].header.cmpaz = cmpaz;
        sac[k].header.cmpinc = cmpinc;
    }
    return;
}


text:
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_statistics_double.h>
#include <math.h>
#include "tz_error.h"
#include "tz_constant.h"
#include "tz_voxel_linked_list.h"
#include "tz_stack_sampling.h"
#include "tz_voxel_graphics.h"
#include "tz_neurotrace.h"
#include "tz_stack_math.h"
#include "tz_stack_utils.h"
#include "tz_objdetect.h"
#include "tz_voxeltrans.h"
#include "tz_stack_stat.h"
#include "tz_stack_draw.h"
#include "tz_geo3d_vector.h"
#include "tz_stack_bwmorph.h"
#include "tz_stack_bwdist.h"

INIT_EXCEPTION_MAIN(e)

int main(int argc, char* argv[])
{
  Stack *seedimg = Read_Stack("../data/fly_neuron_objseed.tif");
  Object_3d_List *objs = Stack_Find_Object(seedimg, 1, 0);

  Stack *stack2 = NULL;
  Object_3d *obj = NULL;
  stack2 = NULL;

  Object_3d_List *objs2 = NULL;

  while (objs != NULL) {
    stack2 = Read_Stack("../data/fly_neuron_bundle.tif");
    obj = Stack_Grow_Object(stack2, 1, objs->data->voxels[0]);
    if (obj->size > 1000) {
      if (objs2 == NULL) {
	objs2 = Object_3d_List_New();
	objs2->data = obj;
      } else {
	Object_3d_List_Add(&objs2, obj);
	Print_Object_3d_Info(obj);
      }
    } else {
      Kill_Object_3d(obj);
    }
    objs = objs->next;
  }

  Print_Object_3d_List_Compact(objs2);

  Stack *sig = Read_Stack("../data/fly_neuron.tif");
  //Stack *sig = Read_Stack("../data/fly_neuron.tif");
  Stack *stack = Translate_Stack(sig, COLOR, 0);
  //Stack *stack = Read_Stack("../data/blobimg2.tif");

  Stack *traced = Make_Stack(GREY, sig->width, sig->height, sig->depth);
  One_Stack(traced);

  Stack *soma = Read_Stack("../data/fly_neuron_soma.tif");
  Stack_Threshold(soma, 1);
  Stack_Binarize(soma);
  Stack_Not(soma, soma);
  Stack_And(traced, soma, traced);
  Kill_Stack(soma);

  Neurochain *chain = New_Neurochain();
  Neurochain *chain2 = NULL;
  while (objs2 != NULL) {
    if (objs2->data->size > 100) {
      Set_Neuroseg(&(chain->seg), 5, 5, 12, -TZ_PI_2, 0.0);
      Set_Position(chain->position, objs2->data->voxels[0][0],
		   objs2->data->voxels[0][1], objs2->data->voxels[0][2]);

#if 1
      Try {
	Fit_Neuroseg(sig, &(chain->seg), chain->position, TRUE);
      } Catch(e) {
	TZ_ERROR(e);
      }
#endif

#if 1
      Try {
	chain2 = Trace_Neuron(sig, chain, BOTH, traced);
      } Catch(e) {
	TZ_ERROR(e);
      }
#endif
      Neurochain_Erase(traced,  chain2);
      chain = Extend_Neurochain(chain2);
    }
    objs2 = objs2->next;
  }

  Neurochain_Label(stack, Neurochain_Head(chain));
  Write_Stack("../data/fly_neuron_trace.tif", stack);

  return 0;
}

ext: c


text:
#ifndef __UTILS_UTILS_H_INCLUDED__
#include <iostream>
#include <vector>
#include <Eigen/Dense>
#include <gsl/gsl>
#include <numeric>
#include <random>
#define __UTILS_UTILS_H_INCLUDED__

template<typename Matrix>
void removeRow(Matrix& matrix, gsl::index rowToRemove)
{
    gsl::index numRows = matrix.rows()-1;
    gsl::index numCols = matrix.cols();

    if( rowToRemove < numRows )
        matrix.block(rowToRemove,0,numRows-rowToRemove,numCols) = matrix.block(rowToRemove+1,0,numRows-rowToRemove,numCols).eval();

    matrix.conservativeResize(numRows,numCols);
}
template<typename Matrix>
void removeColumn(Matrix& matrix, gsl::index colToRemove)
{
    gsl::index numRows = matrix.rows();
    gsl::index numCols = matrix.cols()-1;

    if( colToRemove < numCols )
        matrix.block(0,colToRemove,numRows,numCols-colToRemove) = matrix.block(0,colToRemove+1,numRows,numCols-colToRemove).eval();

    matrix.conservativeResize(numRows,numCols);
}

template <typename T>
void select(std::vector<T>& result, const std::vector<T>& in, const std::vector<typename std::vector<T>::size_type>& s) {
    result.reserve(s.size());
    std::transform(s.begin(), s.end(), std::back_inserter(result),
                   [&in](typename std::vector<T>::size_type idx) {
                       return in.at(idx);
                   });
}

std::tuple<std::vector<gsl::index>,std::vector<gsl::index> > TrainTestSplit(gsl::index nrows,double testRatio=0.3){
    assert(testRatio>=0 && testRatio<=1);
    gsl::index ntest = (gsl::index) nrows*(testRatio);
    gsl::index ntrain = nrows - ntest;
    std::vector<int> idx(nrows);
    std::iota(idx.begin(), idx.end(), 0);

    std::vector<gsl::index> trainIdx, testIdx;
    trainIdx.resize(ntrain);
    testIdx.resize(ntest);
    std::sample(idx.begin(), idx.end(), std::back_inserter(testIdx),
                ntest, std::mt19937{std::random_device{}()});
    std::sort(testIdx.begin(),testIdx.end());
    std::set_difference(idx.begin(), idx.end(), trainIdx.begin(), trainIdx.end(),
                        std::inserter(testIdx, testIdx.begin()));
    return {trainIdx,testIdx};
}

std::tuple<std::vector<gsl::index>,std::vector<gsl::index> > TrainTestSplit(std::vector<gsl::index> inputIdx,double testRatio=0.3){
    assert(testRatio>=0 && testRatio<=1);
    gsl::index nrows = inputIdx.size();
    gsl::index ntest = (gsl::index) nrows*(testRatio);
    gsl::index ntrain = nrows - ntest;
    std::vector<int> idx(nrows);
    std::iota(idx.begin(), idx.end(), 0);

    std::vector<gsl::index> trainIdx, testIdx;
    trainIdx.resize(ntrain);
    testIdx.resize(ntest);
    std::sample(inputIdx.begin(), inputIdx.end(), std::back_inserter(testIdx),
                ntest, std::mt19937{std::random_device{}()});
    std::sort(testIdx.begin(),testIdx.end());
    std::sort(inputIdx.begin(),inputIdx.end());
    std::set_difference(inputIdx.begin(), inputIdx.end(), trainIdx.begin(), trainIdx.end(),
                        std::inserter(testIdx, testIdx.begin()));
    return {trainIdx,testIdx};
}

template<typename T>
std::vector<T> linspace(T start_in, T end_in, int num_in)
{

    std::vector<T> linspaced;

    double start = static_cast<double>(start_in);
    double end = static_cast<double>(end_in);
    double num = static_cast<double>(num_in);

    if (num == 0) { return linspaced; }
    if (num == 1)
    {
        linspaced.push_back(start);
        return linspaced;
    }

    double delta = (end - start) / (num - 1);
    for(int i=0; i < num-1; ++i)
    {
        linspaced.push_back(start + delta * i);
    }
    linspaced.push_back(end);
    return linspaced;
}

template<typename T>
void print_vector(std::vector<T> vec)
{
    std::cout << "size: " << vec.size() << std::endl;
    for (T & d : vec)
        std::cout << d << " ";
    std::cout << std::endl;
}
#endif

ext: h

max_stars_repo_name: zhxinyu/cosan

max_issues_repo_name: zhxinyu/cosan

max_forks_repo_name: zhxinyu/cosan


text:
#ifndef BCT_TEST_H
#define BCT_TEST_H

#include <bct/bct.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <octave/oct.h>
#include <vector>

namespace bct_test {
	Matrix from_gsl(const gsl_vector*, int = 0);
	Matrix from_gsl(const gsl_matrix*, int = 0, int = 0);
	NDArray from_gsl(const std::vector<gsl_matrix*>, int = 0);
	gsl_vector* to_gslv(const Matrix, int = 0);
	gsl_matrix* to_gslm(const Matrix, int = 0, int = 0);
	std::vector<gsl_matrix*> to_gsl(const NDArray, int = 0);
}

#include "bct_test.cpp"

#define MATRIX_TO_SCALAR_FUNCTION(function_name) \
	DEFUN_DLD(function_name##_cpp, args, , "Wrapper for C++ function.") { \
		if (args.length() == 0) { \
			return octave_value_list(); \
		} \
		Matrix m = args(0).matrix_value(); \
		if (!error_state) { \
			gsl_matrix* m_gsl = bct_test::to_gslm(m); \
			octave_value ret = octave_value(bct::function_name(m_gsl)); \
			gsl_matrix_free(m_gsl); \
			return ret; \
		} else { \
			return octave_value_list(); \
		} \
	}

#define MATRIX_TO_VECTOR_FUNCTION(function_name) \
	DEFUN_DLD(function_name##_cpp, args, , "Wrapper for C++ function.") { \
		if (args.length() == 0) { \
			return octave_value_list(); \
		} \
		Matrix m = args(0).matrix_value(); \
		if (!error_state) { \
			gsl_matrix* m_gsl = bct_test::to_gslm(m); \
			gsl_vector* ret_gsl = bct::function_name(m_gsl); \
			octave_value ret = bct_test::from_gsl(ret_gsl); \
			gsl_matrix_free(m_gsl); \
			gsl_vector_free(ret_gsl); \
			return ret; \
		} else { \
			return octave_value_list(); \
		} \
	}

#define MATRIX_TO_MATRIX_FUNCTION(function_name) \
	DEFUN_DLD(function_name##_cpp, args, , "Wrapper for C++ function.") { \
		if (args.length() == 0) { \
			return octave_value_list(); \
		} \
		Matrix m = args(0).matrix_value(); \
		if (!error_state) { \
			gsl_matrix* m_gsl = bct_test::to_gslm(m); \
			gsl_matrix* ret_gsl = bct::function_name(m_gsl); \
			octave_value ret = bct_test::from_gsl(ret_gsl); \
			gsl_matrix_free(m_gsl); \
			gsl_matrix_free(ret_gsl); \
			return ret; \
		} else { \
			return octave_value_list(); \
		} \
	}

#endif

max_stars_repo_name: devuci/bct-cpp

max_issues_repo_name: devuci/bct-cpp

max_forks_repo_name: devuci/bct-cpp


text:
/*
 * Copyright (c) 2010-2018 The University of Tennessee and The University
 *                         of Tennessee Research Foundation.  All rights
 *                         reserved.
 * Copyright (c) 2013      Inria. All rights reserved.
 * $COPYRIGHT
 *
 * @precisions normal z -> s d c
 *
 */
#include "dplasma.h"
#include <math.h>
#include <lapacke.h>
#include "parsec/data_dist/matrix/two_dim_rectangle_cyclic.h"

/**
 *******************************************************************************
 *
 * @ingroup dplasma_complex64_check
 *
 * check_zpotrf - Check the correctness of the Cholesky factorization computed
 * Cholesky functions with the following criteria:
 *
 *    \f[ ||L'L-A||_oo/(||A||_oo.N.eps) < 60. \f]
 *
 *  or
 *
 *    \f[ ||UU'-A||_oo/(||A||_oo.N.eps) < 60. \f]
 *
 *  where A is the original matrix, and L, or U, the result of the Cholesky
 *  factorization.
 *
 *******************************************************************************
 *
 * @param[in,out] parsec
 *          The parsec context of the application that will run the operation.
 *
 * @param[in] loud
 *          The level of verbosity required.
 *
 * @param[in] uplo
 *          = PlasmaUpper: Upper triangle of A and A0 are referenced;
 *          = PlasmaLower: Lower triangle of A and A0 are referenced.
 *
 * @param[in] A
 *          Descriptor of the distributed matrix A result of the Cholesky
 *          factorization. Holds L or U. If uplo == PlasmaUpper, the only the
 *          upper part is referenced, otherwise if uplo == PlasmaLower, the
 *          lower part is referenced.
 *
 * @param[in] A0
 *          Descriptor of the original distributed matrix A before
 *          factorization. If uplo == PlasmaUpper, the only the upper part is
 *          referenced, otherwise if uplo == PlasmaLower, the lower part is
 *          referenced.
 *
 *******************************************************************************
 *
 * @return
 *          \retval 1, if the result is incorrect
 *          \retval 0, if the result is correct
 *
 ******************************************************************************/
int check_zpotrf( parsec_context_t *parsec, int loud,
                  PLASMA_enum uplo,
                  parsec_tiled_matrix_dc_t *A,
                  parsec_tiled_matrix_dc_t *A0 )
{
    two_dim_block_cyclic_t *twodA = (two_dim_block_cyclic_t *)A0;
    two_dim_block_cyclic_t LLt;
    int info_factorization;
    double Rnorm = 0.0;
    double Anorm = 0.0;
    double result = 0.0;
    int M = A->m;
    int N = A->n;
    double eps = LAPACKE_dlamch_work('e');
    PLASMA_enum side;

    two_dim_block_cyclic_init(&LLt, matrix_ComplexDouble, matrix_Tile,
                              A->super.nodes, twodA->grid.rank,
                              A->mb, A->nb, M, N, 0, 0,
                              M, N, twodA->grid.strows, twodA->grid.stcols, twodA->grid.rows);

    LLt.mat = parsec_data_allocate((size_t)LLt.super.nb_local_tiles *
                                  (size_t)LLt.super.bsiz *
                                  (size_t)parsec_datadist_getsizeoftype(LLt.super.mtype));

    dplasma_zlaset( parsec, PlasmaUpperLower, 0., 0.,(parsec_tiled_matrix_dc_t *)&LLt );
    dplasma_zlacpy( parsec, uplo, A, (parsec_tiled_matrix_dc_t *)&LLt );

    /* Compute LL' or U'U  */
    side = (uplo == PlasmaUpper ) ? PlasmaLeft : PlasmaRight;
    dplasma_ztrmm( parsec, side, uplo, PlasmaConjTrans, PlasmaNonUnit, 1.0,
                   A, (parsec_tiled_matrix_dc_t*)&LLt);

    /* compute LL' - A or U'U - A */
    dplasma_ztradd( parsec, uplo, PlasmaNoTrans,
                    -1.0, A0, 1., (parsec_tiled_matrix_dc_t*)&LLt);

    Anorm = dplasma_zlanhe(parsec, PlasmaInfNorm, uplo, A0);
    Rnorm = dplasma_zlanhe(parsec, PlasmaInfNorm, uplo,
                           (parsec_tiled_matrix_dc_t*)&LLt);

    result = Rnorm / ( Anorm * N * eps ) ;

    if ( loud > 2 ) {
        printf("============\n");
        printf("Checking the Cholesky factorization \n");

        if ( loud > 3 )
            printf( "-- ||A||_oo = %e, ||L'L-A||_oo = %e\n", Anorm, Rnorm );

        printf("-- ||L'L-A||_oo/(||A||_oo.N.eps) = %e \n", result);
    }

    if ( isnan(Rnorm)  || isinf(Rnorm)  ||
         isnan(result) || isinf(result) ||
         (result > 60.0) )
    {
        if( loud ) printf("-- Factorization is suspicious ! \n");
        info_factorization = 1;
    }
    else
    {
        if( loud ) printf("-- Factorization is CORRECT ! \n");
        info_factorization = 0;
    }

    parsec_data_free(LLt.mat); LLt.mat = NULL;
    parsec_tiled_matrix_dc_destroy( (parsec_tiled_matrix_dc_t*)&LLt);

    return info_factorization;
}

/**
 *******************************************************************************
 *
 * @ingroup dplasma_complex64_check
 *
 * check_zaxmb - Returns the result of the following test
 *
 *    \f[ (|| A x - b ||_oo / ((||A||_oo * ||x||_oo + ||b||_oo) * N * eps) ) < 60. \f]
 *
 *  where A is the original matrix, b the original right hand side, and x the
 *  solution computed through any factorization.
 *
 *******************************************************************************
 *
 * @param[in,out] parsec
 *          The parsec context of the application that will run the operation.
 *
 * @param[in] loud
 *          The level of verbosity required.
 *
 * @param[in] uplo
 *          = PlasmaUpper: Upper triangle of A is referenced;
 *          = PlasmaLower: Lower triangle of A is referenced.
 *
 * @param[in] A
 *          Descriptor of the distributed matrix A result of the Cholesky
 *          factorization. Holds L or U. If uplo == PlasmaUpper, the only the
 *          upper part is referenced, otherwise if uplo == PlasmaLower, the
 *          lower part is referenced.
 *
 * @param[in,out] b
 *          Descriptor of the original distributed right hand side b.
 *          On exit, b is overwritten by (b - A * x).
 *
 * @param[in] x
 *          Descriptor of the solution to the problem, x.
 *
 *******************************************************************************
 *
 * @return
 *          \retval 1, if the result is incorrect
 *          \retval 0, if the result is correct
 *
 ******************************************************************************/
int check_zaxmb( parsec_context_t *parsec, int loud,
                 PLASMA_enum uplo,
                 parsec_tiled_matrix_dc_t *A,
                 parsec_tiled_matrix_dc_t *b,
                 parsec_tiled_matrix_dc_t *x )
{
    int info_solution;
    double Rnorm = 0.0;
    double Anorm = 0.0;
    double Bnorm = 0.0;
    double Xnorm, result;
    int N = b->m;
    double eps = LAPACKE_dlamch_work('e');

    Anorm = dplasma_zlanhe(parsec, PlasmaInfNorm, uplo, A);
    Bnorm = dplasma_zlange(parsec, PlasmaInfNorm, b);
    Xnorm = dplasma_zlange(parsec, PlasmaInfNorm, x);

    /* Compute b - A*x */
    dplasma_zhemm( parsec, PlasmaLeft, uplo, -1.0, A, x, 1.0, b);

    Rnorm = dplasma_zlange(parsec, PlasmaInfNorm, b);

    result = Rnorm / ( ( Anorm * Xnorm + Bnorm ) * N * eps ) ;

    if ( loud > 2 ) {
        printf("============\n");
        printf("Checking the Residual of the solution \n");
        if ( loud > 3 )
            printf( "-- ||A||_oo = %e, ||X||_oo = %e, ||B||_oo= %e, ||A X - B||_oo = %e\n",
                    Anorm, Xnorm, Bnorm, Rnorm );

        printf("-- ||Ax-B||_oo/((||A||_oo||x||_oo+||B||_oo).N.eps) = %e \n", result);
    }

    if (  isnan(Xnorm) || isinf(Xnorm) || isnan(result) || isinf(result) || (result > 60.0) ) {
        if( loud ) printf("-- Solution is suspicious ! \n");
        info_solution = 1;
    }
    else{
        if( loud ) printf("-- Solution is CORRECT ! \n");
        info_solution = 0;
    }

    return info_solution;
}

/**
 *******************************************************************************
 *
 * @ingroup dplasma_complex64_check
 *
 * check_zpoinv - Returns the result of the following test
 *
 *    \f[ (|| I - A * A^(-1) ||_one / (||A||_one * ||A^(-1)||_one * N * eps) ) < 10. \f]
 *
 *  where A is the original matrix, and Ainv the result of a cholesky inversion.
 *
 *******************************************************************************
 *
 * @param[in,out] parsec
 *          The parsec context of the application that will run the operation.
 *
 * @param[in] loud
 *          The level of verbosity required.
 *
 * @param[in] uplo
 *          = PlasmaUpper: Upper triangle of A is referenced;
 *          = PlasmaLower: Lower triangle of A is referenced.
 *
 * @param[in] A
 *          Descriptor of the distributed original matrix A.
 *          A must be two_dim_block_cyclic and fully generated.
 *
 * @param[in] Ainv
 *          Descriptor of the computed distributed A inverse.
 *
 *******************************************************************************
 *
 * @return
 *          \retval 1, if the result is incorrect
 *          \retval 0, if the result is correct
 *
 ******************************************************************************/
int check_zpoinv( parsec_context_t *parsec, int loud,
                  PLASMA_enum uplo,
                  parsec_tiled_matrix_dc_t *A,
                  parsec_tiled_matrix_dc_t *Ainv )
{
    two_dim_block_cyclic_t *twodA = (two_dim_block_cyclic_t *)A;
    two_dim_block_cyclic_t Id;
    int info_solution;
    double Anorm, Ainvnorm, Rnorm;
    double eps, result;

    eps = LAPACKE_dlamch_work('e');

    two_dim_block_cyclic_init(&Id, matrix_ComplexDouble, matrix_Tile,
                               A->super.nodes, twodA->grid.rank,
                               A->mb, A->nb, A->n, A->n, 0, 0,
                               A->n, A->n, twodA->grid.strows, twodA->grid.stcols, twodA->grid.rows);

    Id.mat = parsec_data_allocate((size_t)Id.super.nb_local_tiles *
                                  (size_t)Id.super.bsiz *
                                  (size_t)parsec_datadist_getsizeoftype(Id.super.mtype));

    dplasma_zlaset( parsec, PlasmaUpperLower, 0., 1., (parsec_tiled_matrix_dc_t *)&Id);

    /* Id - A^-1 * A */
    dplasma_zhemm(parsec, PlasmaLeft, uplo,
                  -1., Ainv, A,
                  1., (parsec_tiled_matrix_dc_t *)&Id );

    Anorm    = dplasma_zlanhe( parsec, PlasmaOneNorm, uplo, A );
    Ainvnorm = dplasma_zlanhe( parsec, PlasmaOneNorm, uplo, Ainv );
    Rnorm    = dplasma_zlange( parsec, PlasmaOneNorm, (parsec_tiled_matrix_dc_t*)&Id );

    result = Rnorm / ( (Anorm*Ainvnorm)*A->n*eps );
    if ( loud > 2 ) {
        printf("  ||A||_one = %e, ||A^(-1)||_one = %e, ||I - A * A^(-1)||_one = %e, result = %e\n",
               Anorm, Ainvnorm, Rnorm, result);
    }

    if ( isinf(Ainvnorm) || isnan(result) || isinf(result) || (result > 10.0) ) {
        info_solution = 1;
    }
    else {
        info_solution = 0;
    }

    parsec_data_free(Id.mat);
    parsec_tiled_matrix_dc_destroy((parsec_tiled_matrix_dc_t*)&Id);

    return info_solution;
}

max_stars_repo_name: NLAFET/ABFT

max_issues_repo_name: NLAFET/ABFT

max_forks_repo_name: NLAFET/ABFT


text:
#include <math.h>
#include <cblas.h>

void matmult_mnk(int M, int N, int K, double **A, double **B, double **C) {
	for (int l = 0; l < M*N; l++) {
		C[0][l] = 0;
	}
    	for (int m = 0; m < M; m++) {
    		for (int n = 0; n < N; n++) {
			for (int k = 0; k < K; k++) {
				C[m][n] += A[m][k] * B[k][n];
		}
    	}
    }
}

void matmult_mkn(int M, int N, int K, double **A, double **B, double **C) {
	for (int l = 0; l < M*N; l++) {
		C[0][l] = 0;
	}
    for (int m = 0; m < M; m++) {
    	for (int k = 0; k < K; k++) {
			for (int n = 0; n < N; n++) {
				C[m][n] += A[m][k] * B[k][n];
			}
    	}
    }
}

void matmult_nmk(int M, int N, int K, double **A, double **B, double **C) {
	for (int l = 0; l < M*N; l++) {
		C[0][l] = 0;
	}
    for (int n = 0; n < N; n++) {
    	for (int m = 0; m < M; m++) {
			for (int k = 0; k < K; k++) {
				C[m][n] += A[m][k] * B[k][n];
			}
    	}
    }
}

void matmult_nkm(int M, int N, int K, double **A, double **B, double **C) {
	for (int l = 0; l < M*N; l++) {
		C[0][l] = 0;
	}
    for (int n = 0; n < N; n++) {
    	for (int k = 0; k < K; k++) {
			for (int m = 0; m < M; m++) {
				C[m][n] += A[m][k] * B[k][n];
			}
    	}
    }
}

void matmult_kmn(int M, int N, int K, double **A, double **B, double **C) {
	for (int l = 0; l < M*N; l++) {
		C[0][l] = 0;
	}
    for (int k = 0; k < K; k++) {
    	for (int m = 0; m < M; m++) {
			for (int n = 0; n < N; n++) {
				C[m][n] += A[m][k] * B[k][n];
			}
    	}
    }
}

void matmult_knm(int M, int N, int K, double **A, double **B, double **C) {
	for (int l = 0; l < M*N; l++) {
		C[0][l] = 0;
	}
    for (int k = 0; k < K; k++) {
    	for (int n = 0; n < N; n++) {
			for (int m = 0; m < M; m++) {
				C[m][n] += A[m][k] * B[k][n];
			}
    	}
    }
}

void matmult_nat(int M, int N, int K, double **A, double **B, double **C) {
	for (int m = 0; m < M; m++) {
		for (int n = 0; n < N; n++) {
			C[m][n] = 0;
			for (int k = 0; k < K; k++) {
				C[m][n] += A[m][k] * B[k][n];
			}
    	}
    }
}

void matmult_lib(int M, int N, int K, double **A, double **B, double **C) {
	int layout, TRANSA, TRANSB, LDA, LDB, LDC;
	double alpha, beta;

    layout = 101; // rowmajor
	TRANSA = 111; // A is not transposed
	TRANSB = 111; // B is not transposed

	LDA = fmax(1,K); // leading dimension of A
	LDB = fmax(1,N); // leading dimension of B
	LDC = fmax(1,N); // leading dimension of C

	alpha = 1.0; // no scaling
	beta = 0.0; //

	cblas_dgemm(layout, TRANSA, TRANSB, M, N, K, alpha, *A, LDA, *B, LDB, beta, *C, LDC);
}

void matmult_blk(int M, int N, int K, double **A, double **B, double **C, int bs)
{
	// MKN

	bs = fmax(1, fmin(bs, K));

	for (int i = 0; i < M*N; i++)
		C[0][i] = 0.0;

	for (int m0 = 0; m0 < M; m0 += bs)
	{
		for (int k0 = 0; k0 < K; k0 += bs)
		{
			for (int n0 = 0; n0 < N; n0 += bs)
			{
				for (int m = m0; m < fmin(m0 + bs, M); m++)
				{
					for (int k = k0; k < fmin(k0 + bs, K); k++)
					{
						for (int n = n0; n < fmin(n0 + bs, N); n++)
							C[m][n] += A[m][k] * B[k][n];
					}
				}
			}
		}
	}
}

max_stars_repo_name: anderslaunerbaek/HPC

max_issues_repo_name: anderslaunerbaek/HPC

max_forks_repo_name: anderslaunerbaek/HPC


text:
#ifndef EIGEN_BENCH_UTIL_H
#define EIGEN_BENCH_UTIL_H

#include <Eigen/Core>
#include "BenchTimer.h"

using namespace std;
using namespace Eigen;

#include <boost/preprocessor/repetition/enum_params.hpp>
#include <boost/preprocessor/repetition.hpp>
#include <boost/preprocessor/seq.hpp>
#include <boost/preprocessor/array.hpp>
#include <boost/preprocessor/arithmetic.hpp>
#include <boost/preprocessor/comparison.hpp>
#include <boost/preprocessor/punctuation.hpp>
#include <boost/preprocessor/punctuation/comma.hpp>
#include <boost/preprocessor/stringize.hpp>

template<typename MatrixType> void initMatrix_random(MatrixType& mat) __attribute__((noinline));
template<typename MatrixType> void initMatrix_random(MatrixType& mat)
{
  mat.setRandom();// = MatrixType::random(mat.rows(), mat.cols());
}

template<typename MatrixType> void initMatrix_identity(MatrixType& mat) __attribute__((noinline));
template<typename MatrixType> void initMatrix_identity(MatrixType& mat)
{
  mat.setIdentity();
}

#ifndef __INTEL_COMPILER
#define DISABLE_SSE_EXCEPTIONS()  { \
  int aux; \
  asm( \
  "stmxcsr   %[aux]           \n\t" \
  "orl       $32832, %[aux]   \n\t" \
  "ldmxcsr   %[aux]           \n\t" \
  : : [aux] "m" (aux)); \
}
#else
#define DISABLE_SSE_EXCEPTIONS()
#endif

#ifdef BENCH_GMM
#include <gmm/gmm.h>
template <typename EigenMatrixType, typename GmmMatrixType>
void eiToGmm(const EigenMatrixType& src, GmmMatrixType& dst)
{
  dst.resize(src.rows(),src.cols());
  for (int j=0; j<src.cols(); ++j)
    for (int i=0; i<src.rows(); ++i)
      dst(i,j) = src.coeff(i,j);
}
#endif

#ifdef BENCH_GSL
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_eigen.h>
template <typename EigenMatrixType>
void eiToGsl(const EigenMatrixType& src, gsl_matrix** dst)
{
  for (int j=0; j<src.cols(); ++j)
    for (int i=0; i<src.rows(); ++i)
      gsl_matrix_set(*dst, i, j, src.coeff(i,j));
}
#endif

#ifdef BENCH_UBLAS
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/vector.hpp>
template <typename EigenMatrixType, typename UblasMatrixType>
void eiToUblas(const EigenMatrixType& src, UblasMatrixType& dst)
{
  dst.resize(src.rows(),src.cols());
  for (int j=0; j<src.cols(); ++j)
    for (int i=0; i<src.rows(); ++i)
      dst(i,j) = src.coeff(i,j);
}
template <typename EigenType, typename UblasType>
void eiToUblasVec(const EigenType& src, UblasType& dst)
{
  dst.resize(src.size());
  for (int j=0; j<src.size(); ++j)
      dst[j] = src.coeff(j);
}
#endif

#endif // EIGEN_BENCH_UTIL_H


text:
#pragma once

#include <gsl\gsl>
#include <winrt\Windows.Foundation.h>
#include <d3d11.h>
#include "DrawableGameComponent.h"
#include "MatrixHelper.h"
#include "DirectionalLight.h"

namespace Library
{
	class ProxyModel;
}

namespace Rendering
{
	class TransparencyMaterial;

	class TransparencyDemo final : public Library::DrawableGameComponent
	{
	public:
		TransparencyDemo(Library::Game& game, const std::shared_ptr<Library::Camera>& camera);
		TransparencyDemo(const TransparencyDemo&) = delete;
		TransparencyDemo(TransparencyDemo&&) = default;
		TransparencyDemo& operator=(const TransparencyDemo&) = default;
		TransparencyDemo& operator=(TransparencyDemo&&) = default;
		~TransparencyDemo();

		float AmbientLightIntensity() const;
		void SetAmbientLightIntensity(float intensity);

		float DirectionalLightIntensity() const;
		void SetDirectionalLightIntensity(float intensity);

		const DirectX::XMFLOAT3& LightDirection() const;
		void RotateDirectionalLight(DirectX::XMFLOAT2 amount);

		float SpecularIntensity() const;
		void SetSpecularIntensity(float intensity);

		float SpecularPower() const;
		void SetSpecularPower(float power);

		float FogStart() const;
		void SetFogStart(float fogStart);

		float FogRange() const;
		void SetFogRange(float fogRange);

		virtual void Initialize() override;
		virtual void Update(const Library::GameTime& gameTime) override;
		virtual void Draw(const Library::GameTime& gameTime) override;

	private:
		inline static const float RotationRate{ DirectX::XM_PI };

		std::shared_ptr<TransparencyMaterial> mMaterial;
		DirectX::XMFLOAT4X4 mWorldMatrix{ Library::MatrixHelper::Identity };
		winrt::com_ptr<ID3D11Buffer> mVertexBuffer;
		std::uint32_t mVertexCount{ 0 };
		Library::DirectionalLight mDirectionalLight;
		std::unique_ptr<Library::ProxyModel> mProxyModel;
		bool mUpdateMaterial{ true };
	};
}


text:
/* siman/siman_test.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Mark Galassi
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

#include <config.h>
#include <math.h>
#include <stdio.h>
#include <string.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_siman.h>
#include <gsl/gsl_ieee_utils.h>

/* A cute example of this program to then look at stuff would be:
   ./siman_test D3 | grep -v "^#" | awk '{print $3}' | sed 's/::/ /gp' > output

  then go into gnuplot and try typing:
  gnuplot> set parametric
  gnuplot> splot 'output' with lines

  Or you can plot from a pipe:
  gnuplot> plot '<../siman/siman_test | grep -v "^#"' using 1:2 with lines
  gnuplot> plot '<../siman/siman_test | grep -v "^#"' using 1:3 with lines
  and so forth.  */

/* set up parameters for this simulated annealing run */

#define N_TRIES 200		/* how many points do we try before stepping */
#define ITERS_FIXED_T 1000	/* how many iterations for each T? */
#define STEP_SIZE 1.0		/* max step size in random walk */
#define K 1.0			/* Boltzmann constant */
#define T_INITIAL 0.008		/* initial temperature */
#define MU_T 1.003		/* damping factor for temperature */
#define T_MIN 2.0e-6

gsl_siman_params_t params = {N_TRIES, ITERS_FIXED_T, STEP_SIZE,
			     K, T_INITIAL, MU_T, T_MIN};

double square(double x);

double test_E_1D(Element x);
void test_step_1D(const gsl_rng * r, Element *x_p, double step_size);
void print_pos_1D(Element x);
double distance_1D(Element x, Element y);
double test_E_2D(Element x);
void test_step_2D(const gsl_rng * r, Element *x_p, double step_size);
void print_pos_2D(Element x);
double distance_2D(Element x, Element y);
double test_E_3D(Element x);
void test_step_3D(const gsl_rng * r, Element *x_p, double step_size);
void print_pos_3D(Element x);
double distance_3D(Element x, Element y);
double E1(void *xp);
double M1(void *xp, void *yp);
void S1(const gsl_rng * r, void *xp, double step_size);
void P1(void *xp);

int main(int argc, char *argv[])
{
  Element x0;			/* initial guess for search */

/*   double x_initial = 2.5; */
  double x_initial = -10.0;

  gsl_rng * r = gsl_rng_alloc (gsl_rng_env_setup()) ;

  gsl_ieee_env_setup ();

  gsl_siman_solve(r, &x_initial, E1, S1, M1, P1, NULL, NULL, NULL,
		  sizeof(double), params);

  return 0;

  if (argc != 2) {
/*     fprintf(stderr, "usage: %s [D1, D2, D3, CA]\n", argv[0]); */
    fprintf(stderr, "usage: %s [D1 | D2 | D3]\n", argv[0]);
    return 1;
  }

  printf("#testing the simulated annealing routines\n");
  if (strcmp(argv[1], "D1") == 0) {
    x0.D1 = 12.0;
    printf("#one dimensional problem, x0 = %f\n", x0.D1);
/*     gsl_siman_Usolve(r, &x0, test_E_1D, test_step_1D, distance_1D, */
/* 		    print_pos_1D, params); */
    return 0;
  }

  if (strcmp(argv[1], "D2") == 0) {
    x0.D2[0] = 12.0;
    x0.D2[1] = 5.5;
    printf("#two dimensional problem, (x0,y0) = (%f,%f)\n",
	   x0.D2[0], x0.D2[1]);
/*     gsl_siman_Usolve(r, &x0, test_E_2D, test_step_2D, distance_2D, */
/* 		    print_pos_2D, params); */
    return 0;
  }

  if (strcmp(argv[1], "D3") == 0) {
    x0.D3[0] = 12.2;
    x0.D3[1] = 5.5;
    x0.D3[2] = -15.5;
    printf("#three dimensional problem, (x0,y0,z0) = (%f,%f,%f)\n",
	   x0.D3[0], x0.D3[1], x0.D3[2]);
/*     gsl_siman_Usolve(r, &x0, test_E_3D, test_step_3D, distance_3D, */
/* 		    print_pos_3D, params); */
  }
/*
  x0.D2[0] = 12.2;
  x0.D2[1] = 5.5;

  gsl_siman_solve(r, &x0, test_E_2D, test_step_2D, distance_2D, print_pos_2D, params);
*/
/*
  x0.D3[0] = 12.2;
  x0.D3[1] = 5.5;
  x0.D3[2] = -15.5;

  gsl_siman_solve(r, &x0, test_E_3D, test_step_3D, distance_3D, print_pos_3D, params);
  */

  return 0;
}

double square(double x)
{
  return x*x;
}

double test_E_1D(Element x)
{
  double val = x.D1;
/*return sin(sin(val*val) - cos(val)) + cos(sin(val) + sin(val)*sin(val));*/
  return exp(-square(val-1))*sin(8*val);
/*  return 1.0/(square(x-1.0) + 1.0)*sin(8.0*x); */
/*  return 1.0/(square(x-1.0) + 1.0)*sin(8.0*x) + x*x/100.0; */
}

/* takes a step for the test function; max distance: step_size.
 * the new point is put in x_p and returned.
 */
void test_step_1D(const gsl_rng * r, Element *x_p, double step_size)
{
  double old_x = x_p->D1;
  double new_x;

  new_x = gsl_rng_uniform(r);
  new_x = new_x*2*step_size;
  new_x = new_x - step_size + old_x;

  x_p->D1 = new_x;
}

/* simple routine to print out a position value */
void print_pos_1D(Element x)
{
  printf("%12g", x.D1);
}

/* a metric for the 2D space */
double distance_1D(Element x, Element y)
{
  return fabs(y.D1 - x.D1);
}

/* a 2-D function to be minimized */
double test_E_2D(Element x)
{
  double old_x = x.D2[0], old_y = x.D2[1];

  return exp(-square(old_x-1) - square(old_y - 0.8))*sin(8*old_x + 8 * old_y);
}

/* takes a step for the test function; max distance: step_size.  the
   new point is put in x_p and returned. */
void test_step_2D(const gsl_rng * r, Element *x_p, double step_size)
{
  double old_x =  x_p->D2[0], old_y = x_p->D2[1], new_x, new_y;

  new_x = gsl_rng_uniform(r);
  new_x = new_x*2*step_size;
  new_x = new_x - step_size + old_x;

  new_y = gsl_rng_uniform(r);
  new_y = new_y*2*step_size;
  new_y = new_y - step_size + old_y;

  x_p->D2[0] = new_x;
  x_p->D2[1] = new_y;
}

/* simple routine to print out a position value */
void print_pos_2D(Element x)
{
  printf("%g::%g", x.D2[0], x.D2[1]);
}

/* a metric for the 2D space */
double distance_2D(Element x, Element y)
{
  return sqrt(square(y.D2[0]-x.D2[0]) + square(y.D2[1]-x.D2[1]));
}

/**********************************************/
/************ 3-dimensional search ************/
/**********************************************/

/* a 3-D function to be minimized */
double test_E_3D(Element x)
{
  return exp(-square(x.D3[0]-1) - square(x.D3[1] - 0.8)
    - square(x.D3[2] - 0.8)) * sin(8*x.D3[0] + 8*x.D3[1] + 8*x.D3[2])
      + (square(x.D3[0]) + square(x.D3[1]) + square(x.D3[2]))/10000.0;
}

/* takes a step for the test function; max distance: step_size.
 * the new point is put in x_p and returned.
 */
void test_step_3D(const gsl_rng * r, Element *x_p, double step_size)
{
  double old_x = x_p->D3[0], old_y = x_p->D3[1], old_z = x_p->D3[2];
  double new_x, new_y, new_z;

  new_x = gsl_rng_uniform(r);
  new_x = new_x*2*step_size;
  new_x = new_x - step_size + old_x;

  new_y = gsl_rng_uniform(r);
  new_y = new_y*2*step_size;
  new_y = new_y - step_size + old_y;

  new_z = gsl_rng_uniform(r);
  new_z = new_z*2*step_size;
  new_z = new_z - step_size + old_z;

  x_p->D3[0] = new_x;
  x_p->D3[1] = new_y;
  x_p->D3[2] = new_z;
}

/* simple routine to print out a position value */
void print_pos_3D(Element x)
{
  printf("%g::%g::%g", x.D3[0], x.D3[1], x.D3[2]);
}

/* a metric for the 2D space */
double distance_3D(Element x, Element y)
{
  return sqrt(square(y.D3[0]-x.D3[0]) + square(y.D3[1]-x.D3[1])
	      + square(y.D3[2]-x.D3[2]));
}

/* now some functions to test in one dimension */
double E1(void *xp)
{
  double x = * ((double *) xp);

  return exp(-square(x-1))*sin(8*x) - exp(-square(x-1000))*0.89;
/*   return exp(-square(x-1))*sin(8*x); */
}

double M1(void *xp, void *yp)
{
  double x = *((double *) xp);
  double y = *((double *) yp);

  return fabs(x - y);
}

void S1(const gsl_rng * r, void *xp, double step_size)
{
  double old_x = *((double *) xp);
  double new_x;

  new_x = gsl_rng_uniform(r)*2*step_size - step_size + old_x;
/*   new_x = new_x*2*step_size; */
/*   new_x = new_x - step_size + old_x; */

  memcpy(xp, &new_x, sizeof(new_x));
}

void P1(void *xp)
{
  printf(" %12g ", *((double *) xp));
}


text:
/* -*- Mode: C++; tab-width: 4; c-basic-offset: 4; indent-tabs-mode: nil -*- */
/*
 *     Copyright 2018 Couchbase, Inc.
 *
 *   Licensed under the Apache License, Version 2.0 (the "License");
 *   you may not use this file except in compliance with the License.
 *   You may obtain a copy of the License at
 *
 *       http://www.apache.org/licenses/LICENSE-2.0
 *
 *   Unless required by applicable law or agreed to in writing, software
 *   distributed under the License is distributed on an "AS IS" BASIS,
 *   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *   See the License for the specific language governing permissions and
 *   limitations under the License.
 */

#pragma once

#include "engine_error.h"
#include "types.h"

#include <gsl/gsl>

namespace cb {
namespace audit {
namespace document {
enum class Operation;
}
} // namespace audit
} // namespace cb

struct ServerDocumentIface {
    virtual ~ServerDocumentIface() = default;

    /**
     * This callback is called from the underlying engine right before
     * it is linked into the list of available documents (it is currently
     * not visible to anyone). The engine should have validated all
     * properties set in the document by the client and the core, and
     * assigned a new CAS number for the document (and sequence number if
     * the underlying engine use those).
     *
     * The callback may at this time do post processing of the document
     * content (it is allowed to modify the content data, but not
     * reallocate or change the size of the data in any way).
     *
     * Given that the engine MAY HOLD LOCKS when calling this function
     * the core is *NOT* allowed to acquire *ANY* locks (except for doing
     * some sort of memory allocation for a temporary buffer).
     *
     * @param cookie The cookie provided to the engine for the storage
     *               command which may (which may hold more context)
     * @param info the items underlying data
     * @return ENGINE_SUCCESS means that the underlying engine should
     *                        proceed to link the item. All other
     *                        error codes means that the engine should
     *                        *NOT* link the item
     */
    virtual ENGINE_ERROR_CODE pre_link(gsl::not_null<const void*> cookie,
                                       item_info& info) = 0;

    /**
     * This callback is called from the underlying engine right before
     * a particular document expires. The callback is responsible examining
     * the value and possibly returning a new and modified value.
     *
     * @param itm_info info pertaining to the item that is to be expired.
     * @return std::string empty if the value required no modification, not
     *         empty then the string contains the modified value. When not empty
     *         the datatype of the new value is datatype xattr only.
     *
     * @throws std::bad_alloc in case of memory allocation failure
     */
    virtual std::string pre_expiry(const item_info& itm_info) = 0;

    /**
     * Add an entry to the audit trail for access to the document specified
     * in the key for this cookie.
     *
     * @param cookie The cookie representing the operation
     * @param operation The type of access for the operation
     */
    virtual void audit_document_access(
            gsl::not_null<const void*> cookie,
            cb::audit::document::Operation operation) = 0;
};


text:
#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include <assert.h>
#include <gsl/gsl_vector.h>

#define MAX(A,B) (((A)>(B))?(A):(B))

void reoshift(gsl_vector *v, int h)
{
  if ( h > 0 ) {
    gsl_vector *temp = gsl_vector_alloc(v->size);
    gsl_vector_view p = gsl_vector_subvector(v, 0, v->size - h);
    gsl_vector_view p1 = gsl_vector_subvector(temp, h, v->size - h);
    gsl_vector_memcpy(&p1.vector, &p.vector);
    p = gsl_vector_subvector(temp, 0, h);
    gsl_vector_set_zero(&p.vector);
    gsl_vector_memcpy(v, temp);
    gsl_vector_free(temp);
  }
}

gsl_vector *poly_long_div(gsl_vector *n, gsl_vector *d, gsl_vector **r)
{
  gsl_vector *nt = NULL, *dt = NULL, *rt = NULL, *d2 = NULL, *q = NULL;
  int gn, gt, gd;

  if ( (n->size >= d->size) && (d->size > 0) && (n->size > 0) ) {
    nt = gsl_vector_alloc(n->size); assert(nt != NULL);
    dt = gsl_vector_alloc(n->size); assert(dt != NULL);
    rt = gsl_vector_alloc(n->size); assert(rt != NULL);
    d2 = gsl_vector_alloc(n->size); assert(d2 != NULL);
    gsl_vector_memcpy(nt, n);
    gsl_vector_set_zero(dt); gsl_vector_set_zero(rt);
    gsl_vector_view p = gsl_vector_subvector(dt, 0, d->size);
    gsl_vector_memcpy(&p.vector, d);
    gsl_vector_memcpy(d2, dt);
    gn = n->size - 1;
    gd = d->size - 1;
    gt = 0;

    while( gsl_vector_get(d, gd) == 0 ) gd--;

    while ( gn >= gd ) {
      reoshift(dt, gn-gd);
      double v = gsl_vector_get(nt, gn)/gsl_vector_get(dt, gn);
      gsl_vector_set(rt, gn-gd, v);
      gsl_vector_scale(dt, v);
      gsl_vector_sub(nt, dt);
      gt = MAX(gt, gn-gd);
      while( (gn>=0) && (gsl_vector_get(nt, gn) == 0.0) ) gn--;
      gsl_vector_memcpy(dt, d2);
    }

    q = gsl_vector_alloc(gt+1); assert(q != NULL);
    p = gsl_vector_subvector(rt, 0, gt+1);
    gsl_vector_memcpy(q, &p.vector);
    if ( r != NULL ) {
      if ( (gn+1) > 0 ) {
	*r = gsl_vector_alloc(gn+1); assert( *r != NULL );
	p = gsl_vector_subvector(nt, 0, gn+1);
	gsl_vector_memcpy(*r, &p.vector);
      } else {
	*r = gsl_vector_alloc(1); assert( *r != NULL );
	gsl_vector_set_zero(*r);
      }
    }
    gsl_vector_free(nt); gsl_vector_free(dt);
    gsl_vector_free(rt); gsl_vector_free(d2);
    return q;
  } else {
    q = gsl_vector_alloc(1); assert( q != NULL );
    gsl_vector_set_zero(q);
    if ( r != NULL ) {
      *r = gsl_vector_alloc(n->size); assert( *r != NULL );
      gsl_vector_memcpy(*r, n);
    }
    return q;
  }
}

void poly_print(gsl_vector *p)
{
  int i;
  for(i=p->size-1; i >= 0; i--) {
    if ( i > 0 )
      printf("%lfx^%d + ",
	     gsl_vector_get(p, i), i);
    else
      printf("%lf\n", gsl_vector_get(p, i));
  }
}

gsl_vector *create_poly(int d, ...)
{
  va_list al;
  int i;
  gsl_vector *r = NULL;

  va_start(al, d);
  r = gsl_vector_alloc(d); assert( r != NULL );

  for(i=0; i < d; i++)
    gsl_vector_set(r, i, va_arg(al, double));

  return r;
}

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
//This computes "the" DCT (discrete cosine transformation) along dim of matrix X.
//This uses a CBLAS matrix multiplication by the DCT-II matrix.

//For complex input X, the output Y is complex, and the DCT is just the DCT of the
//real and imaginary parts separately (following Octave convention).

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

#ifndef M_PI
    #define M_PI 3.14159265358979323846
#endif

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

int dct_cblas_s (float *Y, const float *X, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t dim, const size_t ndct, const int sc);
int dct_cblas_d (double *Y, const double *X, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t dim, const size_t ndct, const int sc);
int dct_cblas_c (float *Y, const float *X, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t dim, const size_t ndct, const int sc);
int dct_cblas_z (double *Y, const double *X, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t dim, const size_t ndct, const int sc);

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

    const size_t N = R*C*S*H;
    const size_t Lx = (dim==0u) ? R : (dim==1u) ? C : (dim==2u) ? S : H;
    if (ndct<Lx) { fprintf(stderr,"error in dct_cblas_s: ndct must be >= Lx (length of vecs in X)\n"); return 1; }

    if (ndct==0u || N==0u) {}
    else if (ndct==1u)
    {
        for (size_t n=N; n>0u; --n, ++X, ++Y) { *Y = *X; }
    }
    else
    {
        //Scaling
        const float s = sc ? 2.0f/sqrtf((float)(2u*ndct)) : 2.0f;
        const float dcsc = sc ? 1.0f/sqrtf((float)ndct) : 2.0f;

        //Initialize DCT-II matrix
        const size_t LN = Lx * ndct;
        const float P_N = (float)(M_PI/(double)ndct);
        float *DCT;
        DCT = (float *)aligned_alloc(sizeof(float),LN*sizeof(float));
        if (!DCT) { fprintf(stderr,"error in dct_cblas_s: problem with aligned_alloc. "); perror("aligned_alloc"); return 1; }
        for (size_t l=0u; l<Lx; ++l, ++DCT) { *DCT = dcsc; }
        for (size_t n=1u; n<ndct; ++n)
        {
            for (size_t l=0u; l<Lx; ++l, ++DCT)
            {
                *DCT = s * cosf(P_N*(0.5f+(float)l)*(float)n);
            }
        }
        DCT -= LN;

        if (Lx==N)
        {
            cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)ndct,(int)Lx,1.0f,DCT,(int)Lx,X,1,0.0f,Y,1);
        }
        else
        {
            const size_t K = (iscolmajor) ? ((dim==0u) ? 1u : (dim==1u) ? R : (dim==2u) ? R*C : R*C*S) : ((dim==0u) ? C*S*H : (dim==1u) ? S*H : (dim==2u) ? H : 1u);
            const size_t B = (iscolmajor && dim==0u) ? C*S*H : K;
            const size_t V = N/Lx, G = V/B;

            if (K==1u && (G==1u || B==1u))
            {
                cblas_sgemm(CblasColMajor,CblasTrans,CblasNoTrans,(int)ndct,(int)V,(int)Lx,1.0f,DCT,(int)Lx,X,(int)Lx,0.0f,Y,(int)ndct);
            }
            else
            {
                for (size_t g=G; g>0u; --g, X+=B*(Lx-1u), Y+=B*(ndct-1u))
                {
                    for (size_t b=B; b>0u; --b, ++X, ++Y)
                    {
                        cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)ndct,(int)Lx,1.0f,DCT,(int)Lx,X,(int)K,0.0f,Y,(int)K);
                    }
                }
            }
        }
        free(DCT);
    }

    return 0;
}

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

    const size_t N = R*C*S*H;
    const size_t Lx = (dim==0u) ? R : (dim==1u) ? C : (dim==2u) ? S : H;
    if (ndct<Lx) { fprintf(stderr,"error in dct_cblas_d: ndct must be >= Lx (length of vecs in X)\n"); return 1; }

    if (ndct==0u || N==0u) {}
    else if (ndct==1u)
    {
        for (size_t n=N; n>0u; --n, ++X, ++Y) { *Y = *X; }
    }
    else
    {
        //Scaling
        const double s = sc ? 2.0/sqrt((double)(2u*ndct)) : 2.0;
        const double dcsc = sc ? 1.0/sqrt((double)ndct) : 2.0;

        //Initialize DCT-II matrix
        const size_t LN = Lx * ndct;
        const double P_N = M_PI/(double)ndct;
        double *DCT;
        DCT = (double *)aligned_alloc(sizeof(double),LN*sizeof(double));
        if (!DCT) { fprintf(stderr,"error in dct_cblas_d: problem with aligned_alloc. "); perror("aligned_alloc"); return 1; }
        for (size_t l=0u; l<Lx; ++l, ++DCT) { *DCT = dcsc; }
        for (size_t n=1u; n<ndct; ++n)
        {
            for (size_t l=0u; l<Lx; ++l, ++DCT)
            {
                *DCT = s * cos(P_N*(0.5+(double)l)*(double)n);
            }
        }
        DCT -= LN;

        if (Lx==N)
        {
            cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)ndct,(int)Lx,1.0,DCT,(int)Lx,X,1,0.0,Y,1);
        }
        else
        {
            const size_t K = (iscolmajor) ? ((dim==0u) ? 1u : (dim==1u) ? R : (dim==2u) ? R*C : R*C*S) : ((dim==0u) ? C*S*H : (dim==1u) ? S*H : (dim==2u) ? H : 1u);
            const size_t B = (iscolmajor && dim==0u) ? C*S*H : K;
            const size_t V = N/Lx, G = V/B;

            if (K==1u && (G==1u || B==1u))
            {
                cblas_dgemm(CblasColMajor,CblasTrans,CblasNoTrans,(int)ndct,(int)V,(int)Lx,1.0,DCT,(int)Lx,X,(int)Lx,0.0,Y,(int)ndct);
            }
            else
            {
                for (size_t g=G; g>0u; --g, X+=B*(Lx-1u), Y+=B*(ndct-1u))
                {
                    for (size_t b=B; b>0u; --b, ++X, ++Y)
                    {
                        cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)ndct,(int)Lx,1.0,DCT,(int)Lx,X,(int)K,0.0,Y,(int)K);
                    }
                }
            }
        }
        free(DCT);
    }

    return 0;
}

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

    const size_t N = R*C*S*H;
    const size_t Lx = (dim==0u) ? R : (dim==1u) ? C : (dim==2u) ? S : H;
    if (ndct<Lx) { fprintf(stderr,"error in dct_cblas_c: ndct must be >= Lx (length of vecs in X)\n"); return 1; }

    if (ndct==0u || N==0u) {}
    else if (ndct==1u)
    {
        for (size_t n=2u*N; n>0u; --n, ++X, ++Y) { *Y = *X; }
    }
    else
    {
        //Scaling
        const float s = sc ? 2.0f/sqrtf((float)(2u*ndct)) : 2.0f;
        const float dcsc = sc ? 1.0f/sqrtf((float)ndct) : 2.0f;

        //Initialize DCT-II matrix
        const size_t LN = Lx * ndct;
        const float P_N = (float)(M_PI/(double)ndct);
        float *DCT;
        DCT = (float *)aligned_alloc(sizeof(float),LN*sizeof(float));
        if (!DCT) { fprintf(stderr,"error in dct_cblas_c: problem with aligned_alloc. "); perror("aligned_alloc"); return 1; }
        for (size_t l=0u; l<Lx; ++l, ++DCT) { *DCT = dcsc; }
        for (size_t n=1u; n<ndct; ++n)
        {
            for (size_t l=0u; l<Lx; ++l, ++DCT)
            {
                *DCT = s * cosf(P_N*(0.5f+(float)l)*(float)n);
            }
        }
        DCT -= LN;

        if (Lx==N)
        {
            cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)ndct,(int)Lx,1.0f,DCT,(int)Lx,X,2,0.0f,Y,2);
            cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)ndct,(int)Lx,1.0f,DCT,(int)Lx,X+1,2,0.0f,Y+1,2);
        }
        else
        {
            const size_t K = (iscolmajor) ? ((dim==0u) ? 1u : (dim==1u) ? R : (dim==2u) ? R*C : R*C*S) : ((dim==0u) ? C*S*H : (dim==1u) ? S*H : (dim==2u) ? H : 1u);
            const size_t B = (iscolmajor && dim==0u) ? C*S*H : K;
            const size_t V = N/Lx, G = V/B;

            if (K==1u && (G==1u || B==1u))
            {
                for (size_t v=0u; v<V; ++v, X+=2u*Lx, Y+=2u*ndct)
                {
                    cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)ndct,(int)Lx,1.0f,DCT,(int)Lx,X,2,0.0f,Y,2);
                    cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)ndct,(int)Lx,1.0f,DCT,(int)Lx,X+1,2,0.0f,Y+1,2);
                }
            }
            else
            {
                for (size_t g=G; g>0u; --g, X+=2u*B*(Lx-1u), Y+=2u*B*(ndct-1u))
                {
                    for (size_t b=B; b>0u; --b, ++X, ++Y)
                    {
                        cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)ndct,(int)Lx,1.0f,DCT,(int)Lx,X,2*(int)K,0.0f,Y,2*(int)K);
                        ++X; ++Y;
                        cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)ndct,(int)Lx,1.0f,DCT,(int)Lx,X,2*(int)K,0.0f,Y,2*(int)K);
                    }
                }
            }
        }
        free(DCT);
    }

    return 0;
}

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

    const size_t N = R*C*S*H;
    const size_t Lx = (dim==0u) ? R : (dim==1u) ? C : (dim==2u) ? S : H;
    if (ndct<Lx) { fprintf(stderr,"error in dct_cblas_z: ndct must be >= Lx (length of vecs in X)\n"); return 1; }

    if (ndct==0u || N==0u) {}
    else if (ndct==1u)
    {
        for (size_t n=2u*N; n>0u; --n, ++X, ++Y) { *Y = *X; }
    }
    else
    {
        //Scaling
        const double s = sc ? 2.0/sqrt((double)(2u*ndct)) : 2.0;
        const double dcsc = sc ? 1.0/sqrt((double)ndct) : 2.0;

        //Initialize DCT-II matrix
        const size_t LN = Lx * ndct;
        const double P_N = M_PI/(double)ndct;
        double *DCT;
        DCT = (double *)aligned_alloc(sizeof(double),LN*sizeof(double));
        if (!DCT) { fprintf(stderr,"error in dct_cblas_z: problem with aligned_alloc. "); perror("aligned_alloc"); return 1; }
        for (size_t l=0u; l<Lx; ++l, ++DCT) { *DCT = dcsc; }
        for (size_t n=1u; n<ndct; ++n)
        {
            for (size_t l=0u; l<Lx; ++l, ++DCT)
            {
                *DCT = s * cos(P_N*(0.5+(double)l)*(double)n);
            }
        }
        DCT -= LN;

        if (Lx==N)
        {
            cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)ndct,(int)Lx,1.0,DCT,(int)Lx,X,2,0.0,Y,2);
            cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)ndct,(int)Lx,1.0,DCT,(int)Lx,X+1,2,0.0,Y+1,2);
        }
        else
        {
            const size_t K = (iscolmajor) ? ((dim==0u) ? 1u : (dim==1u) ? R : (dim==2u) ? R*C : R*C*S) : ((dim==0u) ? C*S*H : (dim==1u) ? S*H : (dim==2u) ? H : 1u);
            const size_t B = (iscolmajor && dim==0u) ? C*S*H : K;
            const size_t V = N/Lx, G = V/B;

            if (K==1u && (G==1u || B==1u))
            {
                for (size_t v=0u; v<V; ++v, X+=2u*Lx, Y+=2u*ndct)
                {
                    cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)ndct,(int)Lx,1.0,DCT,(int)Lx,X,2,0.0,Y,2);
                    cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)ndct,(int)Lx,1.0,DCT,(int)Lx,X+1,2,0.0,Y+1,2);
                }
            }
            else
            {
                for (size_t g=G; g>0u; --g, X+=2u*B*(Lx-1u), Y+=2u*B*(ndct-1u))
                {
                    for (size_t b=B; b>0u; --b, ++X, ++Y)
                    {
                        cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)ndct,(int)Lx,1.0,DCT,(int)Lx,X,2*(int)K,0.0,Y,2*(int)K);
                        ++X; ++Y;
                        cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)ndct,(int)Lx,1.0,DCT,(int)Lx,X,2*(int)K,0.0,Y,2*(int)K);
                    }
                }
            }
        }
        free(DCT);
    }

    return 0;
}

#ifdef __cplusplus
}
}
#endif


text:
/* matrix/gsl_matrix_ushort.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman, Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

#ifndef __GSL_MATRIX_USHORT_H__
#define __GSL_MATRIX_USHORT_H__

#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_vector_ushort.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

typedef struct
{
  size_t size1;
  size_t size2;
  size_t tda;
  unsigned short * data;
  gsl_block_ushort * block;
  int owner;
} gsl_matrix_ushort;

typedef struct
{
  gsl_matrix_ushort matrix;
} _gsl_matrix_ushort_view;

typedef _gsl_matrix_ushort_view gsl_matrix_ushort_view;

typedef struct
{
  gsl_matrix_ushort matrix;
} _gsl_matrix_ushort_const_view;

typedef const _gsl_matrix_ushort_const_view gsl_matrix_ushort_const_view;

/* Allocation */

gsl_matrix_ushort *
gsl_matrix_ushort_alloc (const size_t n1, const size_t n2);

gsl_matrix_ushort *
gsl_matrix_ushort_calloc (const size_t n1, const size_t n2);

gsl_matrix_ushort *
gsl_matrix_ushort_alloc_from_block (gsl_block_ushort * b,
                                   const size_t offset,
                                   const size_t n1,
                                   const size_t n2,
                                   const size_t d2);

gsl_matrix_ushort *
gsl_matrix_ushort_alloc_from_matrix (gsl_matrix_ushort * m,
                                    const size_t k1,
                                    const size_t k2,
                                    const size_t n1,
                                    const size_t n2);

gsl_vector_ushort *
gsl_vector_ushort_alloc_row_from_matrix (gsl_matrix_ushort * m,
                                        const size_t i);

gsl_vector_ushort *
gsl_vector_ushort_alloc_col_from_matrix (gsl_matrix_ushort * m,
                                        const size_t j);

void gsl_matrix_ushort_free (gsl_matrix_ushort * m);

/* Views */

_gsl_matrix_ushort_view
gsl_matrix_ushort_submatrix (gsl_matrix_ushort * m,
                            const size_t i, const size_t j,
                            const size_t n1, const size_t n2);

_gsl_vector_ushort_view
gsl_matrix_ushort_row (gsl_matrix_ushort * m, const size_t i);

_gsl_vector_ushort_view
gsl_matrix_ushort_column (gsl_matrix_ushort * m, const size_t j);

_gsl_vector_ushort_view
gsl_matrix_ushort_diagonal (gsl_matrix_ushort * m);

_gsl_vector_ushort_view
gsl_matrix_ushort_subdiagonal (gsl_matrix_ushort * m, const size_t k);

_gsl_vector_ushort_view
gsl_matrix_ushort_superdiagonal (gsl_matrix_ushort * m, const size_t k);

_gsl_matrix_ushort_view
gsl_matrix_ushort_view_array (unsigned short * base,
                             const size_t n1,
                             const size_t n2);

_gsl_matrix_ushort_view
gsl_matrix_ushort_view_array_with_tda (unsigned short * base,
                                      const size_t n1,
                                      const size_t n2,
                                      const size_t tda);

_gsl_matrix_ushort_view
gsl_matrix_ushort_view_vector (gsl_vector_ushort * v,
                              const size_t n1,
                              const size_t n2);

_gsl_matrix_ushort_view
gsl_matrix_ushort_view_vector_with_tda (gsl_vector_ushort * v,
                                       const size_t n1,
                                       const size_t n2,
                                       const size_t tda);

_gsl_matrix_ushort_const_view
gsl_matrix_ushort_const_submatrix (const gsl_matrix_ushort * m,
                                  const size_t i, const size_t j,
                                  const size_t n1, const size_t n2);

_gsl_vector_ushort_const_view
gsl_matrix_ushort_const_row (const gsl_matrix_ushort * m,
                            const size_t i);

_gsl_vector_ushort_const_view
gsl_matrix_ushort_const_column (const gsl_matrix_ushort * m,
                               const size_t j);

_gsl_vector_ushort_const_view
gsl_matrix_ushort_const_diagonal (const gsl_matrix_ushort * m);

_gsl_vector_ushort_const_view
gsl_matrix_ushort_const_subdiagonal (const gsl_matrix_ushort * m,
                                    const size_t k);

_gsl_vector_ushort_const_view
gsl_matrix_ushort_const_superdiagonal (const gsl_matrix_ushort * m,
                                      const size_t k);

_gsl_matrix_ushort_const_view
gsl_matrix_ushort_const_view_array (const unsigned short * base,
                                   const size_t n1,
                                   const size_t n2);

_gsl_matrix_ushort_const_view
gsl_matrix_ushort_const_view_array_with_tda (const unsigned short * base,
                                            const size_t n1,
                                            const size_t n2,
                                            const size_t tda);

_gsl_matrix_ushort_const_view
gsl_matrix_ushort_const_view_vector (const gsl_vector_ushort * v,
                                    const size_t n1,
                                    const size_t n2);

_gsl_matrix_ushort_const_view
gsl_matrix_ushort_const_view_vector_with_tda (const gsl_vector_ushort * v,
                                             const size_t n1,
                                             const size_t n2,
                                             const size_t tda);

/* Operations */

unsigned short   gsl_matrix_ushort_get(const gsl_matrix_ushort * m, const size_t i, const size_t j);
void    gsl_matrix_ushort_set(gsl_matrix_ushort * m, const size_t i, const size_t j, const unsigned short x);

unsigned short * gsl_matrix_ushort_ptr(gsl_matrix_ushort * m, const size_t i, const size_t j);
const unsigned short * gsl_matrix_ushort_const_ptr(const gsl_matrix_ushort * m, const size_t i, const size_t j);

void gsl_matrix_ushort_set_zero (gsl_matrix_ushort * m);
void gsl_matrix_ushort_set_identity (gsl_matrix_ushort * m);
void gsl_matrix_ushort_set_all (gsl_matrix_ushort * m, unsigned short x);

int gsl_matrix_ushort_fread (FILE * stream, gsl_matrix_ushort * m) ;
int gsl_matrix_ushort_fwrite (FILE * stream, const gsl_matrix_ushort * m) ;
int gsl_matrix_ushort_fscanf (FILE * stream, gsl_matrix_ushort * m);
int gsl_matrix_ushort_fprintf (FILE * stream, const gsl_matrix_ushort * m, const char * format);

int gsl_matrix_ushort_memcpy(gsl_matrix_ushort * dest, const gsl_matrix_ushort * src);
int gsl_matrix_ushort_swap(gsl_matrix_ushort * m1, const gsl_matrix_ushort * m2);

int gsl_matrix_ushort_swap_rows(gsl_matrix_ushort * m, const size_t i, const size_t j);
int gsl_matrix_ushort_swap_columns(gsl_matrix_ushort * m, const size_t i, const size_t j);
int gsl_matrix_ushort_swap_rowcol(gsl_matrix_ushort * m, const size_t i, const size_t j);
int gsl_matrix_ushort_transpose (gsl_matrix_ushort * m);
int gsl_matrix_ushort_transpose_memcpy (gsl_matrix_ushort * dest, const gsl_matrix_ushort * src);

unsigned short gsl_matrix_ushort_max (const gsl_matrix_ushort * m);
unsigned short gsl_matrix_ushort_min (const gsl_matrix_ushort * m);
void gsl_matrix_ushort_minmax (const gsl_matrix_ushort * m, unsigned short * min_out, unsigned short * max_out);

void gsl_matrix_ushort_max_index (const gsl_matrix_ushort * m, size_t * imax, size_t *jmax);
void gsl_matrix_ushort_min_index (const gsl_matrix_ushort * m, size_t * imin, size_t *jmin);
void gsl_matrix_ushort_minmax_index (const gsl_matrix_ushort * m, size_t * imin, size_t * jmin, size_t * imax, size_t * jmax);

int gsl_matrix_ushort_isnull (const gsl_matrix_ushort * m);

int gsl_matrix_ushort_add (gsl_matrix_ushort * a, const gsl_matrix_ushort * b);
int gsl_matrix_ushort_sub (gsl_matrix_ushort * a, const gsl_matrix_ushort * b);
int gsl_matrix_ushort_mul_elements (gsl_matrix_ushort * a, const gsl_matrix_ushort * b);
int gsl_matrix_ushort_div_elements (gsl_matrix_ushort * a, const gsl_matrix_ushort * b);
int gsl_matrix_ushort_scale (gsl_matrix_ushort * a, const double x);
int gsl_matrix_ushort_add_constant (gsl_matrix_ushort * a, const double x);

/***********************************************************************/
/* The functions below are obsolete                                    */
/***********************************************************************/
int gsl_matrix_ushort_get_row(gsl_vector_ushort * v, const gsl_matrix_ushort * m, const size_t i);
int gsl_matrix_ushort_get_col(gsl_vector_ushort * v, const gsl_matrix_ushort * m, const size_t j);
int gsl_matrix_ushort_set_row(gsl_matrix_ushort * m, const size_t i, const gsl_vector_ushort * v);
int gsl_matrix_ushort_set_col(gsl_matrix_ushort * m, const size_t j, const gsl_vector_ushort * v);

extern int gsl_check_range ;

/* inline functions if you are using GCC */

#ifdef HAVE_INLINE
extern inline
unsigned short
gsl_matrix_ushort_get(const gsl_matrix_ushort * m, const size_t i, const size_t j)
{
#ifndef GSL_RANGE_CHECK_OFF
  if (i >= m->size1)
    {
      GSL_ERROR_VAL("first index out of range", GSL_EINVAL, 0) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_VAL("second index out of range", GSL_EINVAL, 0) ;
    }
#endif
  return m->data[i * m->tda + j] ;
}

extern inline
void
gsl_matrix_ushort_set(gsl_matrix_ushort * m, const size_t i, const size_t j, const unsigned short x)
{
#ifndef GSL_RANGE_CHECK_OFF
  if (i >= m->size1)
    {
      GSL_ERROR_VOID("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_VOID("second index out of range", GSL_EINVAL) ;
    }
#endif
  m->data[i * m->tda + j] = x ;
}

extern inline
unsigned short *
gsl_matrix_ushort_ptr(gsl_matrix_ushort * m, const size_t i, const size_t j)
{
#ifndef GSL_RANGE_CHECK_OFF
  if (i >= m->size1)
    {
      GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
    }
#endif
  return (unsigned short *) (m->data + (i * m->tda + j)) ;
}

extern inline
const unsigned short *
gsl_matrix_ushort_const_ptr(const gsl_matrix_ushort * m, const size_t i, const size_t j)
{
#ifndef GSL_RANGE_CHECK_OFF
  if (i >= m->size1)
    {
      GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
    }
#endif
  return (const unsigned short *) (m->data + (i * m->tda + j)) ;
}

#endif

__END_DECLS

#endif /* __GSL_MATRIX_USHORT_H__ */