text:
#include <petsc.h>
#include "flow.h"
#include "mesh.h"

#define DIM 2                   /* Geometric dimension */
#define NFIELDS 3
static const FlowFieldDescriptor PhysicsFields[] = {{"Density", "Den",1},{"Momentum", "Momentum",DIM},{"Energy", "Energy",1},{NULL,0}};
typedef enum {EULER_PAR_GAMMA,EULER_PAR_RHOR,EULER_PAR_AMACH,EULER_PAR_ITANA,EULER_PAR_SIZE} EulerParamIdx;

typedef struct {
    PetscReal rho;
    PetscReal rhoU[DIM];
    PetscReal rhoE;
} EulerNode;
typedef union {
    EulerNode eulernode;
    PetscReal vals[DIM+2];
} EulerNodeUnion;

typedef struct {
    PetscReal gamma;
    PetscReal rhoL;
    PetscReal rhoR;
    PetscReal uL;
    PetscReal uR;
    PetscReal pL;
    PetscReal pR;
    PetscReal maxTime;
    PetscReal length;
} Setup;

typedef struct {
    PetscReal pstar,ustar,rhostarL,astarL,SL,SHL,STL,rhostarR,astarR,SR,SHR,STR,gamm1,gamp1;
} StarState;

typedef struct {
    StarState starState;
    Setup setup;
    PetscReal cfl;
} ProblemSetup;

typedef struct{
    PetscReal f;
    PetscReal fprm;
} PressureFunction;

static PressureFunction f_and_fprm_rarefaction(PetscReal pstar, PetscReal pLR, PetscReal aLR, PetscReal gam, PetscReal gamm1,PetscReal gamp1) {
    // compute value of pressure function for rarefaction
    PressureFunction function;
    function.f = ((2. * aLR) / gamm1) * (pow(pstar / pLR, 0.5 * gamm1 / gam) - 1.);
    function.fprm = (aLR / pLR / gam) * pow(pstar / pLR, -0.5 * gamp1 / gam);
    return function;
}

static PressureFunction f_and_fprm_shock(PetscReal pstar, PetscReal pLR, PetscReal rhoLR, PetscReal gam, PetscReal gamm1, PetscReal gamp1){
    // compute value of pressure function for shock
    PetscReal A = 2./gamp1/rhoLR;
    PetscReal B = gamm1*pLR/gamp1;
    PetscReal sqrtterm =  PetscSqrtReal(A/(pstar+B));
    PressureFunction function;
    function.f = (pstar-pLR)*sqrtterm;
    function.fprm = sqrtterm*(1.-0.5*(pstar-pLR)/(B+pstar));
    return function;
}

#define EPS 1.e-6
#define MAXIT 100
static PetscErrorCode DetermineStarState(const Setup* setup, StarState* starState){
    // compute the speed of sound
    PetscReal aL = PetscSqrtReal(setup->gamma*setup->pL/setup->rhoL);
    PetscReal aR = PetscSqrtReal(setup->gamma*setup->pR/setup->rhoR);

    //first guess pstar based on two-rarefacation approximation
    starState->pstar = aL+aR - 0.5*(setup->gamma-1.)*(setup->uR-setup->uL);
    starState->pstar  = starState->pstar  / (aL/pow(setup->pL,0.5*(setup->gamma-1.)/setup->gamma) + aR/pow(setup->pR,0.5*(setup->gamma-1.)/setup->gamma) );
    starState->pstar = pow(starState->pstar,2.*setup->gamma/(setup->gamma-1.));
    starState->gamm1 = setup->gamma-1.;
    starState->gamp1 = setup->gamma+1.;

    PressureFunction fL;
    if (starState->pstar <= setup->pL){
        fL = f_and_fprm_rarefaction(starState->pstar, setup->pL,aL,setup->gamma,starState->gamm1,starState->gamp1);
    }else{
        fL = f_and_fprm_shock(starState->pstar,setup->pL,setup->rhoL,setup->gamma,starState->gamm1,starState->gamp1);
    }

    PressureFunction fR;
    if (starState->pstar <= setup->pR) {
        fR = f_and_fprm_rarefaction(starState->pstar, setup->pR, aR, setup->gamma, starState->gamm1, starState->gamp1);
    }else {
        fR = f_and_fprm_shock(starState->pstar, setup->pR, setup->rhoR, setup->gamma, starState->gamm1, starState->gamp1);
    }
    PetscReal delu = setup->uR-setup->uL;

    // iterate using Newton-Rapson
    if ((fL.f+fR.f+delu)> EPS) {
        // iterate using Newton-Rapson
        for(PetscInt it =0; it < MAXIT+4; it++){
            PetscReal pold = starState->pstar;
            starState->pstar = pold - (fL.f+fR.f+delu)/(fL.fprm+fR.fprm);

            if(starState->pstar < 0){
                starState->pstar = EPS;
            }

            if(2.0*PetscAbsReal(starState->pstar - pold)/(starState->pstar + pold) < EPS){
                break;
            }else{
                if(starState->pstar < setup->pL){
                    fL = f_and_fprm_rarefaction(starState->pstar, setup->pL,aL,setup->gamma,starState->gamm1,starState->gamp1);
                }else{
                    fL = f_and_fprm_shock(starState->pstar,setup->pL,setup->rhoL,setup->gamma,starState->gamm1,starState->gamp1);
                }
                if (starState->pstar<=setup->pR) {
                    fR = f_and_fprm_rarefaction(starState->pstar, setup->pR, aR, setup->gamma, starState->gamm1, starState->gamp1);
                }else {
                    fR = f_and_fprm_shock(starState->pstar, setup->pR, setup->rhoR, setup->gamma, starState->gamm1, starState->gamp1);
                }
            }

            if (it>MAXIT){
                SETERRQ(PETSC_COMM_WORLD,1,"error in Riemann.find_pstar - did not converage for pstar" );
            }
        }
    }

    // determine rest of star state
    starState->ustar = 0.5*(setup->uL+setup->uR+fR.f-fL.f);

    // left star state
    PetscReal pratio = starState->pstar/setup->pL;
    if (starState->pstar<=setup->pL) {  // rarefaction
        starState->rhostarL = setup->rhoL * PetscPowReal(pratio, 1. / setup->gamma);
        starState->astarL = aL * PetscPowReal(pratio, 0.5 * starState->gamm1 / setup->gamma);
        starState->SHL = setup->uL - aL;
        starState->STL = starState->ustar - starState->astarL;
    }else {  // #shock
        starState->rhostarL = setup->rhoL * (pratio + starState->gamm1 / starState->gamp1) / (starState->gamm1 * pratio / starState->gamp1 + 1.);
        starState->SL = setup->uL - aL * PetscSqrtReal(0.5 * starState->gamp1 / setup->gamma * pratio + 0.5 * starState->gamm1 / setup->gamma);
    }

    // right star state
    pratio = starState->pstar/setup->pR;
    if (starState->pstar<=setup->pR) {  // # rarefaction
        starState->rhostarR = setup->rhoR * PetscPowReal(pratio, 1. / setup->gamma);
        starState->astarR = aR * PetscPowReal(pratio, 0.5 * starState->gamm1 / setup->gamma);
        starState-> SHR = setup->uR + aR;
        starState->STR = starState->ustar + starState->astarR;
    }else {  // shock
        starState->rhostarR = setup->rhoR * (pratio + starState->gamm1 / starState->gamp1) / (starState->gamm1 * pratio / starState->gamp1 + 1.);
        starState->SR = setup->uR + aR * PetscSqrtReal(0.5 * starState->gamp1 / setup->gamma * pratio + 0.5 * starState->gamm1 / setup->gamma);
    }
    return 0;
}

static void SetExactSolutionAtPoint(PetscInt dim, PetscReal xDt, const Setup* setup, const StarState* starState, EulerNode* uu){
    PetscReal p;
    // compute the speed of sound
    PetscReal aL = PetscSqrtReal(setup->gamma*setup->pL/setup->rhoL);
    PetscReal aR = PetscSqrtReal(setup->gamma*setup->pR/setup->rhoR);

    for(PetscInt i =0; i < dim; i++){
        uu->rhoU[i] = 0.0;
    }

    if (xDt <= starState->ustar) {  //# left of contact surface
        if (starState->pstar <= setup->pL) {  // # rarefaction
            if (xDt <= starState->SHL) {
                uu->rho = setup->rhoL;
                p = setup->pL;
                uu->rhoU[0] = setup->uL*uu->rho;
            }else if (xDt <=starState->STL) {  //#SHL < x / t < STL
                PetscReal tmp = 2. / starState->gamp1 + (starState->gamm1 / starState->gamp1 / aL) * (setup->uL - xDt);
                uu->rho = setup->rhoL * pow(tmp, 2. / starState->gamm1);
                uu->rhoU[0] = uu->rho * (2. / starState->gamp1) * (aL + 0.5 * starState->gamm1 * setup->uL + xDt);
                p = setup->pL * pow(tmp, 2. * setup->gamma / starState->gamm1);
            }else {  //# STL < x/t < u*
                uu->rho = starState->rhostarL;
                p = starState->pstar;
                uu->rhoU[0] = uu->rho * starState->ustar;
            }
        }else{  //# shock
            if (xDt<= starState->SL) {  // # xDt < SL
                uu->rho = setup->rhoL;
                p = setup->pL;
                uu->rhoU[0] = uu->rho * setup->uL;
            }else {  //# SL < xDt < ustar
                uu->rho = starState->rhostarL;
                p = starState->pstar;
                uu->rhoU[0] = uu->rho * starState->ustar;
            }
        }
    }else{//# right of contact surface
        if (starState->pstar<=setup->pR) {  //# rarefaction
            if (xDt>= starState->SHR) {
                uu->rho = setup->rhoR;
                p = setup->pR;
                uu->rhoU[0] = uu->rho * setup->uR;
            }else if (xDt >= starState->STR) {  // # SHR < x/t < SHR
                PetscReal tmp = 2./starState->gamp1 - (starState->gamm1/starState->gamp1/aR)*(setup->uR-xDt);
                uu->rho = setup->rhoR*PetscPowReal(tmp,2./starState->gamm1);
                uu->rhoU[0] = uu->rho * (2./starState->gamp1)*(-aR + 0.5*starState->gamm1*setup->uR+xDt);
                p = setup->pR*PetscPowReal(tmp,2.*setup->gamma/starState->gamm1);
            }else{ //# u* < x/t < STR
                uu->rho = starState->rhostarR;
                p = starState->pstar;
                uu->rhoU[0] = uu->rho * starState->ustar;
            }
        }else {//# shock
            if (xDt>= starState->SR) {  // # xDt > SR
                uu->rho = setup->rhoR;
                p = setup->pR;
                uu->rhoU[0] = uu->rho * setup->uR;
            }else {//#ustar < xDt < SR
                uu->rho = starState->rhostarR;
                p = starState->pstar;
                uu->rhoU[0] = uu->rho * starState->ustar;
            }
        }
    }
    PetscReal e =  p/starState->gamm1/uu->rho;
    PetscReal E = e + 0.5*(uu->rhoU[0]/uu->rho)*(uu->rhoU[0]/uu->rho);
    uu->rhoE = uu->rho*E;
}

static PetscErrorCode SetExactSolutionRho(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx){
    ProblemSetup* prob = (ProblemSetup*)ctx;

    PetscReal xDt = (x[0]-prob->setup.length/2)/time;
    EulerNode uu;
    SetExactSolutionAtPoint(dim, xDt, &prob->setup, &prob->starState, &uu);

    u[0] = uu.rho;
    return 0;
}

static PetscErrorCode SetExactSolutionRhoU(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx){
    ProblemSetup* prob = (ProblemSetup*)ctx;

    PetscReal xDt = (x[0]-prob->setup.length/2)/time;
    EulerNode uu;
    SetExactSolutionAtPoint(dim, xDt, &prob->setup, &prob->starState, &uu);
    u[0] = uu.rhoU[0];
    u[1] = uu.rhoU[1];
    return 0;
}

static PetscErrorCode SetExactSolutionRhoE(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx){
    ProblemSetup* prob = (ProblemSetup*)ctx;

    PetscReal xDt = (x[0]-prob->setup.length/2)/time;
    EulerNode uu;
    SetExactSolutionAtPoint(dim, xDt, &prob->setup, &prob->starState, &uu);
    u[0] = uu.rhoE;
    return 0;
}

static PetscErrorCode PrintVector(DM dm, Vec v, PetscInt step, const char * fileName){
    Vec                cellgeom;
    PetscErrorCode ierr = DMPlexGetGeometryFVM(dm, NULL, &cellgeom, NULL);CHKERRQ(ierr);
    PetscInt cStart, cEnd;
        ierr = DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
    DM dmCell;
    ierr = VecGetDM(cellgeom, &dmCell);CHKERRQ(ierr);
    const PetscScalar *cgeom;
    ierr = VecGetArrayRead(cellgeom, &cgeom);CHKERRQ(ierr);
    const PetscScalar      *x;
    ierr = VecGetArrayRead(v, &x);CHKERRQ(ierr);
    // print the header for each file
    char filename[100];
    ierr = PetscSNPrintf(filename,sizeof(filename),"%s.%d.txt",fileName, step);CHKERRQ(ierr);

    PetscInt rank = 0;
    PetscInt size;
    ierr = MPI_Comm_rank(PetscObjectComm(dm), &rank);CHKERRMPI(ierr);
    ierr = MPI_Comm_size(PetscObjectComm(dm), &size);CHKERRMPI(ierr);

    for(PetscInt r =0; r < size; r++ ) {
        if(r == rank) {
            FILE *fptr;
            if(r == 0){
                fptr = fopen(filename, "w");
                fprintf(fptr, "x rho u e\n");
            }else{
                fptr = fopen(filename, "a");
            }
            for (PetscInt c = cStart; c < cEnd; ++c) {
                PetscFVCellGeom *cg;
                const EulerNode *xc;

                ierr = DMPlexPointLocalRead(dmCell, c, cgeom, &cg);
                CHKERRQ(ierr);
                ierr = DMPlexPointGlobalFieldRead(dm, c, 0, x, &xc);
                CHKERRQ(ierr);
                if(xc) {// must be real cell and not ghost
                    PetscReal u0 = xc->rhoU[0] / xc->rho;
                    fprintf(fptr, "%f %f %f %f\n", cg->centroid[0], xc->rho, u0, (xc->rhoE / xc->rho) - 0.5 * u0 * u0);
                }
            }

            fclose(fptr);
        }
        MPI_Barrier(PetscObjectComm(dm));
    }
    ierr = VecRestoreArrayRead(cellgeom, &cgeom);CHKERRQ(ierr);
    ierr = VecRestoreArrayRead(v, &x);CHKERRQ(ierr);
    return 0;
}

static PetscErrorCode MonitorError(TS ts, PetscInt step, PetscReal time, Vec u, void *ctx) {
    PetscFunctionBeginUser;
    PetscErrorCode     ierr;

    // Get the DM
    DM dm;
    ierr = TSGetDM(ts, &dm);CHKERRQ(ierr);

    // Create a copy of the u vector
    Vec e;
    ierr = DMCreateGlobalVector(dm, &e);CHKERRQ(ierr);
    ierr = PetscObjectSetName((PetscObject)e, "exact");CHKERRQ(ierr);

    // Set the values
    PetscErrorCode     (*func[3]) (PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) = {{SetExactSolutionRho, SetExactSolutionRhoU, SetExactSolutionRhoE}};
    void* ctxs[3] ={ctx, ctx, ctx};
    ierr    = DMProjectFunction(dm,time,func,ctxs,INSERT_ALL_VALUES,e);CHKERRQ(ierr);

    // just print to a file for now
    ierr = PrintVector(dm, e, step, "exact");CHKERRQ(ierr);
    ierr = PrintVector(dm, u, step, "solution");CHKERRQ(ierr);

    PetscPrintf(PETSC_COMM_WORLD, "TS %d: %f\n", step, time);

    DMRestoreGlobalVector(dm, &e);
    PetscFunctionReturn(0);
}

static PetscErrorCode PhysicsBoundary_Euler_Mirror(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *a_xI, PetscScalar *a_xG, void *ctx)
{
    const EulerNode *xI = (const EulerNode*)a_xI;
    EulerNode       *xG = (EulerNode*)a_xG;
    ProblemSetup* prob = (ProblemSetup*)ctx;
    PetscFunctionBeginUser;
    xG->rho = xI->rho;
    xG->rhoE = xI->rhoE;
    xG->rhoU[0] = xI->rhoU[0];
    xG->rhoU[1] = xI->rhoU[1];

    PetscFunctionReturn(0);
}

/* PetscReal* => EulerNode* conversion */
static PetscErrorCode PhysicsBoundary_Euler_Left(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *a_xI, PetscScalar *a_xG, void *ctx)
{
    const EulerNode *xI = (const EulerNode*)a_xI;
    EulerNode       *xG = (EulerNode*)a_xG;
    ProblemSetup* prob = (ProblemSetup*)ctx;
    PetscFunctionBeginUser;
    xG->rho = prob->setup.rhoL;
    PetscReal eT = prob->setup.rhoL*((prob->setup.pL /(prob->setup.gamma -1) / prob->setup.rhoL) + 0.5 * prob->setup.uL * prob->setup.uL);
    xG->rhoE = eT;
    xG->rhoU[0] = prob->setup.rhoL * prob->setup.uL;
    xG->rhoU[1] = 0.0;

    PetscFunctionReturn(0);
}

/* PetscReal* => EulerNode* conversion */
static PetscErrorCode PhysicsBoundary_Euler_Right(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *a_xI, PetscScalar *a_xG, void *ctx)
{
    const EulerNode *xI = (const EulerNode*)a_xI;
    EulerNode       *xG = (EulerNode*)a_xG;
    ProblemSetup*        prob = (ProblemSetup*)ctx;
    PetscFunctionBeginUser;
    xG->rho = prob->setup.rhoR;
    PetscReal eT = prob->setup.rhoR*((prob->setup.pR /(prob->setup.gamma -1)/ prob->setup.rhoR) + 0.5 * prob->setup.uR * prob->setup.uR);
    xG->rhoE = eT;
    xG->rhoU[0] = prob->setup.rhoR * prob->setup.uR;
    xG->rhoU[1] = 0.0;

    PetscFunctionReturn(0);
}

static PetscErrorCode  ComputeTimeStep(TS ts){
    DM dm;
    PetscErrorCode ierr = TSGetDM(ts, &dm);CHKERRQ(ierr);
    Vec v;
    TSGetSolution(ts, &v);
    ProblemSetup *problem;
    ierr = DMGetApplicationContext(dm, &problem);CHKERRQ(ierr);

    Vec                cellgeom;
    ierr = DMPlexGetGeometryFVM(dm, NULL, &cellgeom, NULL);CHKERRQ(ierr);
    PetscInt cStart, cEnd;
    ierr = DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
    DM dmCell;
    ierr = VecGetDM(cellgeom, &dmCell);CHKERRQ(ierr);
    const PetscScalar *cgeom;
    ierr = VecGetArrayRead(cellgeom, &cgeom);CHKERRQ(ierr);
    const PetscScalar      *x;
    ierr = VecGetArrayRead(v, &x);CHKERRQ(ierr);

    PetscReal dtMin = 1.0;

    for (PetscInt c = cStart; c < cEnd; ++c) {
        PetscFVCellGeom       *cg;
        const EulerNode           *xc;

        ierr = DMPlexPointLocalRead(dmCell, c, cgeom, &cg);CHKERRQ(ierr);
        ierr = DMPlexPointGlobalFieldRead(dm, c, 0, x, &xc);CHKERRQ(ierr);

        if(xc) {  // must be real cell and not ghost

            PetscReal rho = xc->rho;
            PetscReal u = xc->rhoU[0] / rho;
            PetscReal e = (xc->rhoE / rho) - 0.5 * u * u;
            PetscReal p = (problem->setup.gamma - 1) * rho * e;

            PetscReal a = PetscSqrtReal(problem->setup.gamma * p / rho);
            PetscReal dt = problem->cfl * cg->volume / (a + PetscAbsReal(u));
            dtMin = PetscMin(dtMin, dt);
        }
    }
    PetscInt rank;
    MPI_Comm_rank(PetscObjectComm(ts), &rank);
    printf("dtMin(%d): %f\n", rank,dtMin );

    PetscReal dtMinGlobal;
    ierr = MPI_Allreduce(&dtMin, &dtMinGlobal, 1,MPIU_REAL, MPI_MIN, PetscObjectComm(ts));

    PetscPrintf(PetscObjectComm(ts), "TimeStep: %f\n", dtMinGlobal);
//    ierr = TSSetTimeStep(ts, dtMinGlobal);CHKERRQ(ierr);

    ierr = VecRestoreArrayRead(cellgeom, &cgeom);CHKERRQ(ierr);
    ierr = VecRestoreArrayRead(v, &x);CHKERRQ(ierr);
    return 0;
}

static void ComputeFluxRho(PetscInt dim, PetscInt Nf, const PetscReal *qp, const PetscReal *n, const EulerNode *xL, const EulerNode *xR, PetscInt numConstants, const PetscScalar constants[], PetscReal *flux, void* ctx) {
//    dim	- The spatial dimension
//    Nf	- The number of fields
//    x	- The coordinates at a point on the interface
//    n	- The normal vector to the interface
//    uL	- The state vector to the left of the interface
//    uR	- The state vector to the right of the interface
//    flux	- output array of flux through the interface
//    numConstants	- number of constant parameters
//    constants	- constant parameters
//    ctx	- optional user context
    ProblemSetup* prob = (ProblemSetup*)ctx;

    // this is a hack, only add in flux from left/right
    if(PetscAbs(n[0]) > 1E-5) {
        // Setup Godunov
        Setup currentValues;

        currentValues.gamma = prob->setup.gamma;
        currentValues.length = prob->setup.length;

        if (n[0] > 0) {
            currentValues.rhoL = xL->rho;
            currentValues.uL = xL->rhoU[0] / currentValues.rhoL;
            PetscReal eL = (xL->rhoE / currentValues.rhoL) - 0.5 * currentValues.uL * currentValues.uL;
            currentValues.pL = (prob->setup.gamma - 1) * currentValues.rhoL * eL;

            currentValues.rhoR = xR->rho;
            currentValues.uR = xR->rhoU[0] / currentValues.rhoR;
            PetscReal eR = (xR->rhoE / currentValues.rhoR) - 0.5 * currentValues.uR * currentValues.uR;
            currentValues.pR = (prob->setup.gamma - 1) * currentValues.rhoR * eR;
        }else{
            currentValues.rhoR = xL->rho;
            currentValues.uR = xL->rhoU[0] / currentValues.rhoR;
            PetscReal eR = (xL->rhoE / currentValues.rhoR) - 0.5 * currentValues.uR * currentValues.uR;
            currentValues.pR = (prob->setup.gamma - 1) * currentValues.rhoR * eR;

            currentValues.rhoL = xR->rho;
            currentValues.uL = xR->rhoU[0] / currentValues.rhoL;
            PetscReal eL = (xR->rhoE / currentValues.rhoL) - 0.5 * currentValues.uL * currentValues.uL;
            currentValues.pL = (prob->setup.gamma - 1) * currentValues.rhoL * eL;
        }
        StarState result;
        DetermineStarState(&currentValues, &result);
        EulerNode exact;
        SetExactSolutionAtPoint(dim, 0.0, &currentValues, &result, &exact);

        PetscReal rho = exact.rho;
        PetscReal u = exact.rhoU[0]/rho;
        PetscReal e = (exact.rhoE / rho) - 0.5 * u * u;
        PetscReal p = (prob->setup.gamma-1) * rho * e;

        flux[0] = rho * u * PetscSignReal(n[0]);
        printf("flux qp[%f]: %f  n:%f rL:%f rR:%f\n ", qp[0], flux[0], n[0], currentValues.rhoL, currentValues.rhoR);
//        flux->rhoU[0] = (rho * u * u + p)* PetscSignReal(n[0]);
//        flux->rhoU[1] = 0.0;
//        PetscReal et = e + 0.5 * u * u;
//        flux->rhoE = (rho * u * (et + p / rho))* PetscSignReal(n[0]);

//        printf("%f,%f %f %f,%f\n", qp[0], qp[1], flux[0], n[0], n[1]);mm

    }else{
        flux[0] = 0.0;
//        flux->rhoU[0] =0.0;
//        flux->rhoU[1] = 0.0;
//        flux->rhoE = 0.0;

    }

//    F[0][i]=rho[i]*u[i]
//    F[1][i]=rho[i]*u[i]*u[i]+p[i]
//    et = e[i]+0.5*u[i]*u[i]
//    F[2][i]=rho[i]*u[i]*(et+p[i]/rho[i])

}

static void ComputeFluxU(PetscInt dim, PetscInt Nf, const PetscReal *qp, const PetscReal *n, const EulerNode *xL, const EulerNode *xR, PetscInt numConstants, const PetscScalar constants[], PetscReal *flux, void* ctx) {
//    dim	- The spatial dimension
//    Nf	- The number of fields
//    x	- The coordinates at a point on the interface
//    n	- The normal vector to the interface
//    uL	- The state vector to the left of the interface
//    uR	- The state vector to the right of the interface
//    flux	- output array of flux through the interface
//    numConstants	- number of constant parameters
//    constants	- constant parameters
//    ctx	- optional user context
    ProblemSetup* prob = (ProblemSetup*)ctx;

    // this is a hack, only add in flux from left/right
    if(PetscAbs(n[0]) > 1E-5) {
        // Setup Godunov
        Setup currentValues;

        currentValues.gamma = prob->setup.gamma;
        currentValues.length = prob->setup.length;

        if (n[0] > 0) {
            currentValues.rhoL = xL->rho;
            currentValues.uL = xL->rhoU[0] / currentValues.rhoL;
            PetscReal eL = (xL->rhoE / currentValues.rhoL) - 0.5 * currentValues.uL * currentValues.uL;
            currentValues.pL = (prob->setup.gamma - 1) * currentValues.rhoL * eL;

            currentValues.rhoR = xR->rho;
            currentValues.uR = xR->rhoU[0] / currentValues.rhoR;
            PetscReal eR = (xR->rhoE / currentValues.rhoR) - 0.5 * currentValues.uR * currentValues.uR;
            currentValues.pR = (prob->setup.gamma - 1) * currentValues.rhoR * eR;
        }else{
            currentValues.rhoR = xL->rho;
            currentValues.uR = xL->rhoU[0] / currentValues.rhoR;
            PetscReal eR = (xL->rhoE / currentValues.rhoR) - 0.5 * currentValues.uR * currentValues.uR;
            currentValues.pR = (prob->setup.gamma - 1) * currentValues.rhoR * eR;

            currentValues.rhoL = xR->rho;
            currentValues.uL = xR->rhoU[0] / currentValues.rhoL;
            PetscReal eL = (xR->rhoE / currentValues.rhoL) - 0.5 * currentValues.uL * currentValues.uL;
            currentValues.pL = (prob->setup.gamma - 1) * currentValues.rhoL * eL;
        }
        StarState result;
        DetermineStarState(&currentValues, &result);
        EulerNode exact;
        SetExactSolutionAtPoint(dim, 0.0, &currentValues, &result, &exact);

        PetscReal rho = exact.rho;
        PetscReal u = exact.rhoU[0]/rho;
        PetscReal e = (exact.rhoE / rho) - 0.5 * u * u;
        PetscReal p = (prob->setup.gamma-1) * rho * e;

//        flux[0] = (rho * u) * PetscSignReal(n[0]);
        flux[0] = (rho * u * u + p)* PetscSignReal(n[0]);
        flux[1] = 0.0;
//        PetscReal et = e + 0.5 * u * u;
//        flux->rhoE = (rho * u * (et + p / rho))* PetscSignReal(n[0]);

//        printf("%f,%f %f %f %f %f\n", qp[0], qp[1], flux->rho, flux->rhoU[0], flux->rhoE, n[0]);

    }else{
        flux[0] = 0.0;
        flux[1] = 0.0;
//        flux->rhoU[0] =0.0;
//        flux->rhoU[1] = 0.0;
//        flux->rhoE = 0.0;

    }

//    F[0][i]=rho[i]*u[i]
//    F[1][i]=rho[i]*u[i]*u[i]+p[i]
//    et = e[i]+0.5*u[i]*u[i]
//    F[2][i]=rho[i]*u[i]*(et+p[i]/rho[i])

}

static void ComputeFluxE(PetscInt dim, PetscInt Nf, const PetscReal *qp, const PetscReal *n, const EulerNode *xL, const EulerNode *xR, PetscInt numConstants, const PetscScalar constants[], PetscReal *flux, void* ctx) {
//    dim	- The spatial dimension
//    Nf	- The number of fields
//    x	- The coordinates at a point on the interface
//    n	- The normal vector to the interface
//    uL	- The state vector to the left of the interface
//    uR	- The state vector to the right of the interface
//    flux	- output array of flux through the interface
//    numConstants	- number of constant parameters
//    constants	- constant parameters
//    ctx	- optional user context
    ProblemSetup* prob = (ProblemSetup*)ctx;

    // this is a hack, only add in flux from left/right
    if(PetscAbs(n[0]) > 1E-5) {
        // Setup Godunov
        Setup currentValues;

        currentValues.gamma = prob->setup.gamma;
        currentValues.length = prob->setup.length;

        if (n[0] > 0) {
            currentValues.rhoL = xL->rho;
            currentValues.uL = xL->rhoU[0] / currentValues.rhoL;
            PetscReal eL = (xL->rhoE / currentValues.rhoL) - 0.5 * currentValues.uL * currentValues.uL;
            currentValues.pL = (prob->setup.gamma - 1) * currentValues.rhoL * eL;

            currentValues.rhoR = xR->rho;
            currentValues.uR = xR->rhoU[0] / currentValues.rhoR;
            PetscReal eR = (xR->rhoE / currentValues.rhoR) - 0.5 * currentValues.uR * currentValues.uR;
            currentValues.pR = (prob->setup.gamma - 1) * currentValues.rhoR * eR;
        }else{
            currentValues.rhoR = xL->rho;
            currentValues.uR = xL->rhoU[0] / currentValues.rhoR;
            PetscReal eR = (xL->rhoE / currentValues.rhoR) - 0.5 * currentValues.uR * currentValues.uR;
            currentValues.pR = (prob->setup.gamma - 1) * currentValues.rhoR * eR;

            currentValues.rhoL = xR->rho;
            currentValues.uL = xR->rhoU[0] / currentValues.rhoL;
            PetscReal eL = (xR->rhoE / currentValues.rhoL) - 0.5 * currentValues.uL * currentValues.uL;
            currentValues.pL = (prob->setup.gamma - 1) * currentValues.rhoL * eL;
        }
        StarState result;
        DetermineStarState(&currentValues, &result);
        EulerNode exact;
        SetExactSolutionAtPoint(dim, 0.0, &currentValues, &result, &exact);

        PetscReal rho = exact.rho;
        PetscReal u = exact.rhoU[0]/rho;
        PetscReal e = (exact.rhoE / rho) - 0.5 * u * u;
        PetscReal p = (prob->setup.gamma-1) * rho * e;

//        flux[0] = (rho * u) * PetscSignReal(n[0]);
//        flux[0] = (rho * u * u + p)* PetscSignReal(n[0]);
//        flux[1] = 0.0;
        PetscReal et = e + 0.5 * u * u;
        flux[0] = (rho * u * (et + p / rho))* PetscSignReal(n[0]);

//        printf("%f,%f %f %f %f %f\n", qp[0], qp[1], flux->rho, flux->rhoU[0], flux->rhoE, n[0]);

    }else{
        flux[0] = 0.0;
//        flux[1] = 0.0;
//        flux->rhoU[0] =0.0;
//        flux->rhoU[1] = 0.0;
//        flux->rhoE = 0.0;

    }

//    F[0][i]=rho[i]*u[i]
//    F[1][i]=rho[i]*u[i]*u[i]+p[i]
//    et = e[i]+0.5*u[i]*u[i]
//    F[2][i]=rho[i]*u[i]*(et+p[i]/rho[i])

}

int main(int argc, char **argv)
{
    PetscErrorCode ierr;
    // create the mesh
    // setup the ts
    DM dm;                 /* problem definition */
    TS ts;                 /* timestepper */

    // initialize petsc and mpi
    PetscInitialize(&argc, &argv, NULL, "HELP");

    // Setup the problem
    ProblemSetup problem;

    // case 1 - Sod problem
    problem.setup.rhoL=1.0;
    problem.setup.uL=0.0;
    problem.setup.pL=1.0;
    problem.setup.rhoR=0.125;
    problem.setup.uR=0.0;
    problem.setup.pR=0.1;
    problem.setup.maxTime = 0.25;
    problem.setup.length = 1;
    problem.setup.gamma = 1.4;
    problem.cfl = .5;

    // case 2 - 123 problem - expansion left and expansion right
//    problem.setup.rhoL=1.0;
//    problem.setup.uL=-2.0;
//    problem.setup.pL=0.4;
//    problem.setup.rhoR=1.0;
//    problem.setup.uR=2.0;
//    problem.setup.pR=0.4;
//    problem.setup.maxTime = 0.15;
//    problem.setup.length = 1;
//    problem.setup.gamma = 1.4;
//    problem.cfl = 0.4;

    ierr = TSCreate(PETSC_COMM_WORLD, &ts);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    ierr = TSSetType(ts, TSEULER);CHKERRQ(ierr);

    //PetscErrorCode DMPlexCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscBool simplex, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool interpolate, DM *dm)
    PetscReal start[] = {0.0, 0.0};
    PetscReal end[] = {problem.setup.length, 1};
    PetscInt nx[] = {100, 1};
    DMBoundaryType bcType[] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE};
    ierr = DMPlexCreateBoxMesh(PETSC_COMM_WORLD, DIM, PETSC_FALSE, nx, start, end, bcType, PETSC_TRUE, &dm);CHKERRQ(ierr);

    {
        DM dmDist;

//        ierr = DMSetBasicAdjacency(dm, PETSC_TRUE, PETSC_FALSE);CHKERRQ(ierr);
        ierr = DMPlexDistribute(dm, 1, NULL, &dmDist);CHKERRQ(ierr);
        if (dmDist) {
            ierr = DMDestroy(&dm);CHKERRQ(ierr);
            dm   = dmDist;
        }
    }

    ierr = DMSetFromOptions(dm);CHKERRQ(ierr);

    {
        DM gdm;
        ierr = DMPlexConstructGhostCells(dm, NULL, NULL, &gdm);CHKERRQ(ierr);
        ierr = DMDestroy(&dm);CHKERRQ(ierr);
        dm   = gdm;
    }

//    DMLabel label;
//    ierr = DMGetLabel(dm, "marker", &label );CHKERRQ(ierr);
//    ierr = DMLabelView(label, PETSC_VIEWER_STDOUT_WORLD);;CHKERRQ(ierr);

    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    ierr = TSSetDM(ts, dm);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    // setup the FV

    {// Setup the fields
        PetscInt f, dof;
        for (f=0,dof=0; f < 3; f++) {
            PetscInt newDof = PhysicsFields[f].components;

//            if (newDof == 1) {
                PetscFV           fvm;
                ierr = PetscFVCreate(PETSC_COMM_WORLD, &fvm);CHKERRQ(ierr);

                ierr = PetscFVSetFromOptions(fvm);CHKERRQ(ierr);
                ierr = PetscFVSetNumComponents(fvm, newDof);CHKERRQ(ierr);
                ierr = PetscFVSetSpatialDimension(fvm, DIM);CHKERRQ(ierr);
                ierr = PetscObjectSetName((PetscObject) fvm,PhysicsFields[f].fieldName);CHKERRQ(ierr);

                /* FV is now structured with one field having all physics as components */
                ierr = DMAddField(dm, NULL, (PetscObject) fvm);CHKERRQ(ierr);
//            }
//            else {
//                PetscInt j;
//
//                for (j = 0; j < newDof; j++) {
//                    char     compName[256]  = "Unknown";
//
//                    ierr = PetscSNPrintf(compName,sizeof(compName),"%s_%d",PhysicsFields[f].fieldName,j);CHKERRQ(ierr);
//                    ierr = PetscFVSetComponentName(fvm,dof+j,compName);CHKERRQ(ierr);
//                }
//            }
//            dof += newDof;
        }
    }

    // Compute the star state
    ierr = DetermineStarState(&problem.setup, &problem.starState);CHKERRQ(ierr);

    PetscDS           prob;
    ierr = DMCreateDS(dm);CHKERRQ(ierr);
    ierr = DMGetDS(dm, &prob);CHKERRQ(ierr);
    ierr = DMSetApplicationContext(dm, &problem);CHKERRQ(ierr);

    //TODO: Add flux
    ierr = PetscDSSetRiemannSolver(prob, 0,ComputeFluxRho);CHKERRQ(ierr);
    ierr = PetscDSSetContext(prob, 0, &problem);CHKERRQ(ierr);
    ierr = PetscDSSetRiemannSolver(prob, 1,ComputeFluxU);CHKERRQ(ierr);
    ierr = PetscDSSetContext(prob, 1, &problem);CHKERRQ(ierr);
    ierr = PetscDSSetRiemannSolver(prob, 2,ComputeFluxE);CHKERRQ(ierr);
    ierr = PetscDSSetContext(prob, 2, &problem);CHKERRQ(ierr);
    ierr = PetscDSSetFromOptions(prob);CHKERRQ(ierr);
    //TODO: Apply boundary

    // setup the solution vector, this olds everything
    Vec X;
    ierr = DMCreateGlobalVector(dm, &X);CHKERRQ(ierr);
    ierr = PetscObjectSetName((PetscObject) X, "solution");CHKERRQ(ierr);

//    ierr = DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, NULL);CHKERRQ(ierr);
//    ierr = DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, NULL);CHKERRQ(ierr);
//    ierr = DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, NULL);CHKERRQ(ierr);
//    ierr = DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFVM, NULL);CHKERRQ(ierr);//TODO: This is were we set the RHS function

    const PetscInt idsLeft[]= {4};
    ierr = PetscDSAddBoundary(prob, DM_BC_NATURAL_RIEMANN, "wall left", "Face Sets", 0, 0, NULL, (void (*)(void)) PhysicsBoundary_Euler_Left, NULL, 1, idsLeft, &problem);CHKERRQ(ierr);
    const PetscInt idsRight[]= {2};
    ierr = PetscDSAddBoundary(prob, DM_BC_NATURAL_RIEMANN, "wall right", "Face Sets", 0, 0, NULL, (void (*)(void)) PhysicsBoundary_Euler_Right, NULL, 1, idsRight, &problem);CHKERRQ(ierr);

    const PetscInt mirror[]= {1, 3};
    ierr = PetscDSAddBoundary(prob, DM_BC_NATURAL_RIEMANN, "mirrorWall", "Face Sets", 0, 0, NULL, (void (*)(void)) PhysicsBoundary_Euler_Mirror, NULL, 2, mirror, &problem);CHKERRQ(ierr);

//    ierr = DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, NULL);CHKERRQ(ierr);
    ierr = DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFVM, NULL);CHKERRQ(ierr);//TODO: This is were we set the RHS function
    ierr = TSSetMaxTime(ts,problem.setup.maxTime);CHKERRQ(ierr);
    ierr = TSMonitorSet(ts, MonitorError, &problem, NULL);CHKERRQ(ierr);
    ierr = TSSetTimeStep(ts, 0.0008);CHKERRQ(ierr);
    ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
    ierr = TSSetPostStep(ts, ComputeTimeStep);CHKERRQ(ierr);

    // set the initial conditions
    PetscErrorCode     (*func[3]) (PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) = {SetExactSolutionRho, SetExactSolutionRhoU, SetExactSolutionRhoE};
    void* ctxs[3] ={&problem, &problem, &problem};
    ierr    = DMProjectFunction(dm,0.0,func,ctxs,INSERT_ALL_VALUES,X);CHKERRQ(ierr);

//    PetscInt rank, size;
//    MPI_Comm_rank(PetscObjectComm(dm), &rank);
//    MPI_Comm_size(PetscObjectComm(dm), &size);
//    for(int r =0; r < size; r++) {
//        if(r == rank) {
//            printf("Rank: %d\n", r);
//            ierr = DMView(dm, PETSC_VIEWER_STDOUT_SELF);
//            CHKERRQ(ierr);
////            ierr = DMViewFromOptions(dm, NULL, "-dm_view");
////            CHKERRQ(ierr);
//        }
//        MPI_Barrier(PetscObjectComm(dm));
//    }
//    TSSetMaxSteps(ts, 1);
    ierr = TSSolve(ts,X);CHKERRQ(ierr);
//    ierr = TSGetSolveTime(ts,&ftime);CHKERRQ(ierr);
//    ierr = TSGetStepNumber(ts,&nsteps);
    //    CHKERRQ(ierr);

    return PetscFinalize();

}

max_stars_repo_name: mmcgurn/MattFlowCases

max_issues_repo_name: mmcgurn/MattFlowCases

max_forks_repo_name: mmcgurn/MattFlowCases


text:
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_multifit.h>

/* Fit
 *
 * y = X c
 *
 * where X is an n x p matrix of n observations for p variables.
 *
 * The solution includes a possible standard form Tikhonov regularization:
 *
 * c = (X^T X + lambda^2 I)^{-1} X^T y
 *
 * where lambda^2 is the Tikhonov regularization parameter.
 *
 * The function multifit_linear_svd() must first be called to
 * compute the SVD decomposition of X
 *
 * Inputs: X        - least squares matrix
 *         y        - right hand side vector
 *         tol      - singular value tolerance
 *         lambda   - Tikhonov regularization parameter lambda;
 *                    ignored if <= 0
 *         rank     - (output) effective rank
 *         c        - (output) model coefficient vector
 *         rnorm    - (output) residual norm ||y - X c||
 *         snorm    - (output) solution norm ||c||
 *         work     - workspace
 *
 * Notes:
 * 1) The dimensions of X must match work->n and work->p which are set
 *    by multifit_linear_svd()
 * 2) On input:
 *    work->A contains U
 *    work->Q contains Q
 *    work->S contains singular values
 * 3) If this function is called from gsl_multifit_wlinear(), then
 *    the input y points to work->t, which contains sqrt(W) y. Since
 *    work->t is also used as scratch workspace by this function, we
 *    do the necessary computations with y first to avoid problems.
 * 4) When lambda <= 0, singular values are truncated when:
 *    s_j <= tol * s_0
 */

static int
multifit_linear_solve (const gsl_matrix * X,
                       const gsl_vector * y,
                       const double tol,
                       const double lambda,
                       size_t * rank,
                       gsl_vector * c,
                       double *rnorm,
                       double *snorm,
                       gsl_multifit_linear_workspace * work)
{
  const size_t n = X->size1;
  const size_t p = X->size2;

  if (n != work->n || p != work->p)
    {
      GSL_ERROR("observation matrix does not match workspace", GSL_EBADLEN);
    }
  else if (n != y->size)
    {
      GSL_ERROR("number of observations in y does not match matrix",
                GSL_EBADLEN);
    }
  else if (p != c->size)
    {
      GSL_ERROR ("number of parameters c does not match matrix",
                 GSL_EBADLEN);
    }
  else if (tol <= 0)
    {
      GSL_ERROR ("tolerance must be positive", GSL_EINVAL);
    }
  else
    {
      const double lambda_sq = lambda * lambda;

      double rho_ls = 0.0;     /* contribution to rnorm from OLS */

      size_t j, p_eff;

      /* these inputs are previously computed by multifit_linear_svd() */
      gsl_matrix_view A = gsl_matrix_submatrix(work->A, 0, 0, n, p);
      gsl_matrix_view Q = gsl_matrix_submatrix(work->Q, 0, 0, p, p);
      gsl_vector_view S = gsl_vector_subvector(work->S, 0, p);

      /* workspace */
      gsl_matrix_view QSI = gsl_matrix_submatrix(work->QSI, 0, 0, p, p);
      gsl_vector_view xt = gsl_vector_subvector(work->xt, 0, p);
      gsl_vector_view D = gsl_vector_subvector(work->D, 0, p);
      gsl_vector_view t = gsl_vector_subvector(work->t, 0, n);

      /*
       * Solve y = A c for c
       * c = Q diag(s_i / (s_i^2 + lambda_i^2)) U^T y
       */

      /* compute xt = U^T y */
      gsl_blas_dgemv (CblasTrans, 1.0, &A.matrix, y, 0.0, &xt.vector);

      if (n > p)
        {
          /*
           * compute OLS residual norm = || y - U U^T y ||;
           * for n = p, U U^T = I, so no need to calculate norm
           */
          gsl_vector_memcpy(&t.vector, y);
          gsl_blas_dgemv(CblasNoTrans, -1.0, &A.matrix, &xt.vector, 1.0, &t.vector);
          rho_ls = gsl_blas_dnrm2(&t.vector);
        }

      if (lambda > 0.0)
        {
          /* xt <-- [ s(i) / (s(i)^2 + lambda^2) ] .* U^T y */
          for (j = 0; j < p; ++j)
            {
              double sj = gsl_vector_get(&S.vector, j);
              double f = (sj * sj) / (sj * sj + lambda_sq);
              double *ptr = gsl_vector_ptr(&xt.vector, j);

              /* use D as workspace for residual norm */
              gsl_vector_set(&D.vector, j, (1.0 - f) * (*ptr));

              *ptr *= sj / (sj*sj + lambda_sq);
            }

          /* compute regularized solution vector */
          gsl_blas_dgemv (CblasNoTrans, 1.0, &Q.matrix, &xt.vector, 0.0, c);

          /* compute solution norm */
          *snorm = gsl_blas_dnrm2(c);

          /* compute residual norm */
          *rnorm = gsl_blas_dnrm2(&D.vector);

          if (n > p)
            {
              /* add correction to residual norm (see eqs 6-7 of [1]) */
              *rnorm = sqrt((*rnorm) * (*rnorm) + rho_ls * rho_ls);
            }

          /* reset D vector */
          gsl_vector_set_all(&D.vector, 1.0);
        }
      else
        {
          /* Scale the matrix Q, QSI = Q S^{-1} */

          gsl_matrix_memcpy (&QSI.matrix, &Q.matrix);

          {
            double s0 = gsl_vector_get (&S.vector, 0);
            p_eff = 0;

            for (j = 0; j < p; j++)
              {
                gsl_vector_view column = gsl_matrix_column (&QSI.matrix, j);
                double sj = gsl_vector_get (&S.vector, j);
                double alpha;

                if (sj <= tol * s0)
                  {
                    alpha = 0.0;
                  }
                else
                  {
                    alpha = 1.0 / sj;
                    p_eff++;
                  }

                gsl_vector_scale (&column.vector, alpha);
              }

            *rank = p_eff;
          }

          gsl_blas_dgemv (CblasNoTrans, 1.0, &QSI.matrix, &xt.vector, 0.0, c);

          /* Unscale the balancing factors */
          gsl_vector_div (c, &D.vector);

          *snorm = gsl_blas_dnrm2(c);
          *rnorm = rho_ls;
        }

      return GSL_SUCCESS;
    }
}

ext: c

max_stars_repo_name: agroimpacts/imager

max_issues_repo_name: agroimpacts/imager

max_forks_repo_name: agroimpacts/imager


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

#define NCols 5
#define Nrows 4

float A[] = { 8, 4, 7, 3, 5, 1, 1, 3, 2, 1, 2, 3, 2, 0, 1, 1 , 2, 3, 4, 1};

float x[] = { -1, 2, -1, 1, 2 };

float y[Nrows];
float alpha = 1.0, beta = 0.0;
char tbuf[1024];
int main() {
    int i, j;

    // Print original matrix

    // y = Ax
    cblas_sgemv(CblasRowMajor, CblasNoTrans, Nrows, NCols, alpha, A, NCols, x, 1, beta, y, 1);
    // Print resulting vector
    for (j = 0; j < Nrows; j++) {
        printf(" %f\n", y[j]);
    }

    cblas_sgemv(CblasColMajor, CblasNoTrans, Nrows, NCols, alpha, A, Nrows, x, 1, beta, y, 1);
    // Print resulting vector
    for (j = 0; j < Nrows; j++) {
        printf(" %f\n", y[j]);
    }

    return 0;
}

max_stars_repo_name: ramcn/flappie-x-cpp

max_stars_repo_licenses: Naumen Condor-1.1 MS-PL

max_issues_repo_name: ramcn/flappie-x-cpp

max_issues_repo_licenses: Naumen Condor-1.1 MS-PL

max_forks_repo_name: ramcn/flappie-x-cpp

max_forks_repo_licenses: Naumen Condor-1.1 MS-PL


text:
/* specfunc/gsl_sf_log.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004 Gerard Jungman
 *
 * 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

/* Author:  G. Jungman */

#ifndef __GSL_SF_LOG_H__
#define __GSL_SF_LOG_H__

#include <gsl/gsl_sf_result.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

/* Provide a logarithm function with GSL semantics.
 *
 * exceptions: GSL_EDOM
 */
int gsl_sf_log_e(const double x, gsl_sf_result * result);
double gsl_sf_log(const double x);

/* Log(|x|)
 *
 * exceptions: GSL_EDOM
 */
int gsl_sf_log_abs_e(const double x, gsl_sf_result * result);
double gsl_sf_log_abs(const double x);

/* Complex Logarithm
 *   exp(lnr + I theta) = zr + I zi
 * Returns argument in [-pi,pi].
 *
 * exceptions: GSL_EDOM
 */
int gsl_sf_complex_log_e(const double zr, const double zi, gsl_sf_result * lnr, gsl_sf_result * theta);

/* Log(1 + x)
 *
 * exceptions: GSL_EDOM
 */
int gsl_sf_log_1plusx_e(const double x, gsl_sf_result * result);
double gsl_sf_log_1plusx(const double x);

/* Log(1 + x) - x
 *
 * exceptions: GSL_EDOM
 */
int gsl_sf_log_1plusx_mx_e(const double x, gsl_sf_result * result);
double gsl_sf_log_1plusx_mx(const double x);

#ifdef HAVE_INLINE
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>

extern inline
int
gsl_sf_log_e(const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x <= 0.0) {
    result->val = GSL_NAN;
    result->err = GSL_NAN;
    GSL_ERROR ("domain error", GSL_EDOM);
  }
  else {
    result->val = log(x);
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
}
extern inline
int
gsl_sf_log_abs_e(const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x == 0.0) {
    result->val = GSL_NAN;
    result->err = GSL_NAN;
    GSL_ERROR ("domain error", GSL_EDOM);
  }
  else {
    result->val = log(fabs(x));
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
}
#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_SF_LOG_H__ */

ext: h

max_stars_repo_name: andrewkern/segSiteHMM

max_issues_repo_name: andrewkern/segSiteHMM

max_forks_repo_name: andrewkern/segSiteHMM


text:
/* C functions for running SampleTau from Cython*/

/*System includes*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <sys/stat.h>

/*GSL includes*/
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_blas.h>

/*User includes*/
//#include "c_sample_tau.h"

static gsl_rng *ptGSLRNG;

void c_initRNG()
{
    //const gsl_rng_type * T;

    gsl_rng_env_setup();

    //T = gsl_rng_default;
    ptGSLRNG = gsl_rng_alloc (gsl_rng_mt19937);
}

void c_setRNG(unsigned long int seed)
{
    //printf("GSL RNG initialise %lu\n",seed);
    gsl_rng_set (ptGSLRNG, seed);
}

void c_freeRNG()
{
    gsl_rng_free (ptGSLRNG);
}

void normaliseLog4(double *adLogProb)
{
    double dMax = adLogProb[0], dSum = 0.0;
    int b = 0;

    for(b = 1; b < 4; b++){
        if(adLogProb[b] > dMax){
            dMax = adLogProb[b];
        }
    }

    for(b = 0; b < 4; b++){
        adLogProb[b] = adLogProb[b] - dMax;
        dSum += exp(adLogProb[b]);
    }

    for(b = 0; b < 4; b++){
        adLogProb[b] = exp(adLogProb[b])/dSum;
        //printf("b=%d,p=%f\n",b,adLogProb[b]);
    }

    return;
}

int sample4(double *adProb, double dU){
    double adCProb[4];

    adCProb[0] = adProb[0];
    adCProb[1] = adProb[1] + adCProb[0];
    adCProb[2] = adProb[2] + adCProb[1];
    adCProb[3] = 1.0;

    if(dU < adCProb[0]){
        return 0;
    }
    else if(dU < adCProb[1])
        return 1;
    else if(dU < adCProb[2]){
        return 2;
    }
    else{
        return 3;
    }
}

int c_sample_tau (long *anTau, double* adPi, double *adEta, long* anVariants, int nV, int nG, int nS)
{
    int a = 0, b = 0;
    int g = 0, h = 0, v = 0, s = 0, t = 0;
    int nchange = 0;
    double adPSB[nS][4];
    double adPSBStore[nS][4];
    double dLogProb[4];
    int** anTauIndex = NULL;
    double u = 0.0;

    anTauIndex = (int **) malloc(nV*sizeof(int*));
    if(!anTauIndex)
        goto memoryError;

    for(v = 0; v < nV; v++){
        anTauIndex[v] = (int *) malloc(nG*sizeof(int));
        if(!anTauIndex[v])
            goto memoryError;
        for(g = 0; g < nG; g++){
            for(b = 0; b < 4; b++){
                int vIndex = v*4*nG + 4*g + b;

                if(anTau[vIndex] == 1){
                    anTauIndex[v][g] = b;
                    break;
                }
            }
            //printf("%d,%d,%d\n",v,g,anTauIndex[v][g]);
        }

    }

    //loop V positions
    for(v = 0; v < nV; v++){

        //loop G strains
        for(g = 0; g < nG; g++){

            //calc contribution from all other strains
            for(s = 0; s < nS; s++){
                for(b = 0; b < 4; b++){
                    adPSBStore[s][b] = 0.0;

                    for(h = 0; h < nG; h++){
                        if(h != g){
                            int pIndex = s*nG + h;
                            int eIndex = anTauIndex[v][h]*4 + b;

                            adPSBStore[s][b] += adEta[eIndex]*adPi[pIndex];
                        }
                    }

                }
            }

            for(a = 0; a < 4; a++){
                for(s = 0; s < nS; s++){
                    for(b = 0; b < 4; b++){
                        adPSB[s][b] = adPSBStore[s][b];

                        adPSB[s][b] += adEta[a*4 + b]*adPi[s*nG + g];
                    }
                }
                dLogProb[a] = 0.0;
                for(s = 0; s < nS; s++){
                    for(b = 0; b < 4; b++){
                        int vIndex = v*4*nS + 4*s + b;
                        double temp = ((float) anVariants[vIndex])*log(adPSB[s][b]);
                        //printf("s=%d,b=%d,cont=%f\n",s,b,temp);
                        dLogProb[a] += temp;
                    }
                }
            }
  //          printf("u=%f,p1=%f,p2=%f,p3=%f,p4 =%f\n",u,dLogProb[0],dLogProb[1],dLogProb[2],dLogProb[3]);

            normaliseLog4(dLogProb);

            u = gsl_rng_uniform (ptGSLRNG);
       //     printf("u=%f,p1=%f,p2=%f,p3=%f,p4 =%f\n",u,dLogProb[0],dLogProb[1],dLogProb[2],dLogProb[3]);
            t = sample4(dLogProb, u);
         //   printf("v=%d,g=%d,tnew=%d,told=%d\n",v,g,t,anTauIndex[v][g]);
            if(t != anTauIndex[v][g]){
                int vIndex = v*4*nG + 4*g;
                anTau[vIndex + anTauIndex[v][g]] = 0;
                anTau[vIndex + t] = 1;

                nchange++;
                anTauIndex[v][g] = t;
            }
        } //finish sampling strain g

    }//finish sampling position v

    //Free up tau indices
    for(v = 0; v < nV; v++){
        free(anTauIndex[v]);
    }
    free(anTauIndex);

    return nchange;

    memoryError:
    fprintf(stderr, "Failed allocating memory in c_sample_tau\n");
    fflush(stderr);
    exit(EXIT_FAILURE);
}

max_stars_repo_name: MatthewWolff/DESMAN

max_issues_repo_name: MatthewWolff/DESMAN

max_forks_repo_name: MatthewWolff/DESMAN


text:
//           $Id: pfc_libraries.h 37984 2018-10-27 15:50:30Z p20068 $
//          $URL: https://svn01.fh-hagenberg.at/bin/cepheiden/vocational/teaching/ESD/SPS3/2018-WS/ILV/src/Snippets/bitmap-gsl/pfc_libraries.h $
//     $Revision: 37984 $
//         $Date: 2018-10-27 17:50:30 +0200 (Sa., 27 Okt 2018) $
//       $Author: p20068 $
//
//       Creator: Peter Kulczycki (peter.kulczycki<AT>fh-hagenberg.at)
// Creation Date:
//     Copyright: (c) 2018 Peter Kulczycki (peter.kulczycki<AT>fh-hagenberg.at)
//
//       License: This document contains proprietary information belonging to
//                University of Applied Sciences Upper Austria, Campus
//                Hagenberg. It is distributed under the Boost Software License,
//                Version 1.0 (see http://www.boost.org/LICENSE_1_0.txt).
//
//    Annotation: This file is part of the code snippets handed out during one
//                of my HPC lessons held at the University of Applied Sciences
//                Upper Austria, Campus Hagenberg.

#pragma once

#include "./pfc_macros.h"

// -------------------------------------------------------------------------------------------------

#if defined PFC_DETECTED_COMPILER_NVCC
   #define PFC_DO_NOT_USE_BOOST_UNITS
   #define PFC_DO_NOT_USE_GSL
   #define PFC_DO_NOT_USE_VLD
   #define PFC_DO_NOT_USE_WINDOWS
#endif

// -------------------------------------------------------------------------------------------------

#undef PFC_HAVE_VLD
#undef PFC_VLD_INCLUDED

#if __has_include (<vld.h>) && !defined PFC_DO_NOT_USE_VLD   // Visual Leak Detector (https://kinddragon.github.io/vld)
   #include <vld.h>

   #define PFC_HAVE_VLD
   #define PFC_VLD_INCLUDED

   #pragma message ("PFC: using 'Visual Leak Detector'")
#else
   #pragma message ("PFC: not using 'Visual Leak Detector'")
#endif

// -------------------------------------------------------------------------------------------------

#undef PFC_HAVE_GSL
#undef PFC_GSL_INCLUDED

#if __has_include (<gsl/gsl>) && !defined PFC_DO_NOT_USE_GSL   // Guideline Support Library (https://github.com/Microsoft/GSL)
   #include <gsl/gsl>

   #define PFC_HAVE_GSL
   #define PFC_GSL_INCLUDED

   #pragma message ("PFC: using 'Guideline Support Library'")
#else
   #pragma message ("PFC: not using 'Guideline Support Library'")
#endif

// -------------------------------------------------------------------------------------------------

#undef PFC_HAVE_BOOST_UNITS
#undef PFC_BOOST_UNITS_INCLUDED

#if __has_include (<boost/units/io.hpp>)
#if __has_include (<boost/units/systems/si/length.hpp>)
#if __has_include (<boost/units/systems/si/prefixes.hpp>) && !defined PFC_DO_NOT_USE_BOOST_UNITS
   #include <boost/units/io.hpp>                    // http://www.boost.org
   #include <boost/units/systems/si/length.hpp>     // https://sourceforge.net/projects/boost/files/boost-binaries
   #include <boost/units/systems/si/prefixes.hpp>   //

   #define PFC_HAVE_BOOST_UNITS
   #define PFC_BOOST_UNITS_INCLUDED

   #pragma message ("PFC: using 'Boost.Units'")
#else
   #pragma message ("PFC: not using 'Boost.Units'")
#endif
#endif
#endif

// -------------------------------------------------------------------------------------------------

#undef PFC_HAVE_WINDOWS
#undef PFC_WINDOWS_INCLUDED

#if __has_include (<windows.h>) && !defined PFC_DO_NOT_USE_WINDOWS
   #undef  NOMINMAX
   #define NOMINMAX

   #undef  STRICT
   #define STRICT

   #undef  VC_EXTRALEAN
   #define VC_EXTRALEAN

   #undef  WIN32_LEAN_AND_MEAN
   #define WIN32_LEAN_AND_MEAN

   #include <windows.h>

   #define PFC_HAVE_WINDOWS
   #define PFC_WINDOWS_INCLUDED

   #pragma message ("PFC: using 'windows.h'")
#else
   #pragma message ("PFC: not using 'windows.h'")
#endif


text:
/* ode-initval/rk2imp.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * 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.
 */

/* Runge-Kutta 2, Gaussian implicit */

/* Author:  G. Jungman
 */
#include <config.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv.h>

#include "odeiv_util.h"

typedef struct
{
  double *knu;
  double *ytmp;
}
rk2imp_state_t;

static void *
rk2imp_alloc (size_t dim)
{
  rk2imp_state_t *state = (rk2imp_state_t *) malloc (sizeof (rk2imp_state_t));

  if (state == 0)
    {
      GSL_ERROR_NULL ("failed to allocate space for rk2imp_state",
                      GSL_ENOMEM);
    }

  state->knu = (double *) malloc (dim * sizeof (double));

  if (state->knu == 0)
    {
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for knu", GSL_ENOMEM);
    }

  state->ytmp = (double *) malloc (dim * sizeof (double));

  if (state->ytmp == 0)
    {
      free (state->knu);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for ytmp", GSL_ENOMEM);
    }

  return state;
}

static int
rk2imp_apply (void *vstate,
              size_t dim,
              double t,
              double h,
              double y[],
              double yerr[],
              const double dydt_in[],
              double dydt_out[], const gsl_odeiv_system * sys)
{
  rk2imp_state_t *state = (rk2imp_state_t *) vstate;

  const int iter_steps = 3;
  int status = 0;
  int nu;
  size_t i;

  double *const knu = state->knu;
  double *const ytmp = state->ytmp;

  /* initialization step */
  if (dydt_in != NULL)
    {
      DBL_MEMCPY (knu, dydt_in, dim);
    }
  else
    {
      int s = GSL_ODEIV_FN_EVAL (sys, t, y, knu);
      GSL_STATUS_UPDATE (&status, s);
    }

  /* iterative solution */
  for (nu = 0; nu < iter_steps; nu++)
    {
      for (i = 0; i < dim; i++)
        {
          ytmp[i] = y[i] + 0.5 * h * knu[i];
        }
      {
        int s = GSL_ODEIV_FN_EVAL (sys, t + 0.5 * h, ytmp, knu);
        GSL_STATUS_UPDATE (&status, s);
      }
    }

  /* assignment */
  for (i = 0; i < dim; i++)
    {
      y[i] += h * knu[i];
      yerr[i] = h * h * knu[i];
      if (dydt_out != NULL)
        dydt_out[i] = knu[i];
    }

  return status;
}

static int
rk2imp_reset (void *vstate, size_t dim)
{
  rk2imp_state_t *state = (rk2imp_state_t *) vstate;

  DBL_ZERO_MEMSET (state->knu, dim);
  DBL_ZERO_MEMSET (state->ytmp, dim);

  return GSL_SUCCESS;
}

static unsigned int
rk2imp_order (void *vstate)
{
  rk2imp_state_t *state = (rk2imp_state_t *) vstate;
  state = 0; /* prevent warnings about unused parameters */
  return 2;
}

static void
rk2imp_free (void *vstate)
{
  rk2imp_state_t *state = (rk2imp_state_t *) vstate;
  free (state->knu);
  free (state->ytmp);
  free (state);
}

static const gsl_odeiv_step_type rk2imp_type = { "rk2imp",      /* name */
  1,                            /* can use dydt_in */
  0,                            /* gives exact dydt_out */
  &rk2imp_alloc,
  &rk2imp_apply,
  &rk2imp_reset,
  &rk2imp_order,
  &rk2imp_free
};

const gsl_odeiv_step_type *gsl_odeiv_step_rk2imp = &rk2imp_type;


text:
#ifndef PyGSL_ERROR_HELPER_H
#define PyGSL_ERROR_HELPER_H 1
#include <pygsl/intern.h>
#include <pygsl/utils.h>
#include <gsl/gsl_errno.h>
#include <pygsl/errorno.h>

/*
 * 22 Sep. 2009 Pierre Schnizer
 * Uncomment only if trouble with the gsl error handler (e.g. when using
 * Python with threading support (typical ufuncs). At the time of this writing
 * the error handler would call python to find the approbriate python exception
 *
 * So I used to uncomment the macro as well as the function to ensure that
 * gsl_error was not called any more within the pygsl wrappper
 */
/*
#undef GSL_ERROR
#undef GSL_ERROR_VAL
#undef GSL_ERROR_NULL
#define gsl_error()
*/
/*
 * handle gsl error flags.
 *
 * If a flag arrives check if there was already a python error. If so leave it alone.
 * We cannot return two exceptions.
 *
 * Otherwise:
 *       Should I put an exception up? E.g. some function not conforming to GSL
 *       Convention returning a flag, instead of calling gsl_error?
 *       Currently I follow that idea. But I have no more information about the reason
 *       than the flag.
 *
 * Return:
 *       GSL_SUCCESS ... No errornous call
 *       GSL_FAILURE ...    errornous call
 *
 * If you need to return the flag e.g. "int gsl_odeiv_iterate( ... " use
 * PyGSL_error_flag_to_pyint instead!
 *
 */

PyGSL_API_EXTERN int
PyGSL_error_flag(long flag);

/*
 * Handles gsl_error flags.
 * It differs from the above that it returns the integer.
 *
 * Negative values mean something like go one with the iteration. These are
 * converted to an python integer. Positive values flag a problem. These are
 * converted to python exceptions.
 */
PyGSL_API_EXTERN PyObject *
PyGSL_error_flag_to_pyint(long flag);

/*
 * Add a Python trace back frame to the python interpreter.
 * Input :
 *     module   ... the module. Pass NULL if not known.
 *     filename ... The filename to list in the stack frame. Pass NULL if not
 *                  known.
 *     funcname ... The function name to list in the stack frame. Pass NULL if
 *                  not known.
 *     lineno   ... The Linenumber where the error occurred.
 */
PyGSL_API_EXTERN void
PyGSL_add_traceback(PyObject *module, const char *filename, const char *funcname, int lineno);

PyGSL_API_EXTERN int
PyGSL_warning(const char *, const char*, int, int);

#ifndef _PyGSL_API_MODULE
/* Section for modules importing the functions */
#define PyGSL_error_flag           (*(int (*)(long))                                         PyGSL_API[PyGSL_error_flag_NUM])
#define PyGSL_error_flag_to_pyint  (*(PyObject * (*)(long))                                  PyGSL_API[PyGSL_error_flag_to_pyint_NUM])
#define PyGSL_add_traceback        (*(void (*)(PyObject *, const char *, const char *, int)) PyGSL_API[PyGSL_add_traceback_NUM])
#define PyGSL_warning              (*(int (*)(const char *, const char *, int, int))         PyGSL_API[PyGSL_warning_NUM])

#endif /* _PyGSL_API_MODULE */

#define PyGSL_ERROR_FLAG(flag)                                              \
(((long) flag == GSL_SUCCESS) && (!PyErr_Occurred())) ? GSL_SUCCESS :       \
                     PyGSL_error_flag((long) (flag))

#define PyGSL_ERROR_FLAG_TO_PYINT(flag)                                     \
(((long) flag <= 0) && (!PyErr_Occurred())) ? PyInt_FromLong((long) flag) : \
                     PyGSL_error_flag_to_pyint((long) (flag))

#endif /* PyGSL_ERROR_HELPER_H  */

max_stars_repo_name: juhnowski/FishingRod

max_issues_repo_name: juhnowski/FishingRod

max_forks_repo_name: juhnowski/FishingRod


text:
/****************************************************************************
 *                                                                          *
 *  Author : lukasz.iwaszkiewicz@gmail.com                                  *
 *  ~~~~~~~~                                                                *
 *  License : see COPYING file for details.                                 *
 *  ~~~~~~~~~                                                               *
 ****************************************************************************/

#pragma once
#include <gsl/gsl>

namespace sd {

void init ();
int lsdir (gsl::czstring path);

} // namespace sd

max_stars_repo_name: iwasz/zephyr-grbl-plotter

max_issues_repo_name: iwasz/zephyr-grbl-plotter

max_forks_repo_name: iwasz/zephyr-grbl-plotter


text:
#ifndef SPN_H
#define SPN_H
#include "init_cy_fde.h"
#include "matrix_util.h"
//#include <cblas.h>

typedef struct SemiCircularNet{
  int p;
  int d;
  DCOMPLEX *Pa_h_A;         // d
  DCOMPLEX *Pa_omega;       // 2*d

  DCOMPLEX F_A[2];         //
  DCOMPLEX h_A[2];         //

  DCOMPLEX TG_Ge_sc[4] ;   //
  DCOMPLEX DG_sc[4];       //

  DCOMPLEX temp_T_eta[4];  //
  DCOMPLEX Dh_sc[4];       //

  DCOMPLEX Psigma_G_sc[2]; //

  DCOMPLEX Psigma_h_sc[2] ;//

  DCOMPLEX DG_A[4];        //
  DCOMPLEX Dh_A[4];        //

  DCOMPLEX S[4];           //

  DCOMPLEX temp_mat[4];    //

  DCOMPLEX Psigma_omega[2];//

}SCN;

void
SCN_construct(SCN* self, int p , int d);

void
SCN_init(SCN* self);

void
SCN_init_forward(SCN* self);
void
SCN_init_backward(SCN* self);

void
SCN_destroy(SCN* self);

int
SCN_cauchy(SCN* self);

void
SCN_grad(SCN* self,  int p, int d, double  *a,  double  sigma, \
  DCOMPLEX z,DCOMPLEX *G, DCOMPLEX *omega, DCOMPLEX *omega_sc,\
  DCOMPLEX *o_grad_a, DCOMPLEX *o_grad_sigma);

/** Compute Cauchy transform of SemiCircular( returns total iterations)
* @param Z : input matrix
* @param o_G : out_put Cauchy transform
*
*/
int
cauchy_sc( int p,  int d,  double sigma, DCOMPLEX* Z, \
   int max_iter, double thres,\
   DCOMPLEX* o_G);

 /** Compute Cauchy transform of Signal-Plus-Noise model( returns total iterations)
 * @param Z : input matrix
 * @param o_G_sc : out_put Cauchy transform
 *
 */
int
cauchy_spn(int p_dim, int dim, double* a, double sigma,\
     DCOMPLEX* B,\
     int max_iter,double thres, \
     DCOMPLEX* o_G_sc, DCOMPLEX* o_omega, DCOMPLEX* o_omega_sc);

/*
Only for debug
*/
void
grad_cauchy_spn(int p, int d,  double  *a, double  sigma, \
  DCOMPLEX z, DCOMPLEX *G, DCOMPLEX *omega, DCOMPLEX *omega_sc,\
  DCOMPLEX *o_grad_a, DCOMPLEX *o_grad_sigma);

// transpose of derivation of Ge
// G : 2
// o_DGe: 2 x 2
void TG_Ge( const int p, const int d, const double sigma, \
  const DCOMPLEX *G, DCOMPLEX *o_DGe);

// transpose of derivation of cauchy_sc
// G: 2
// DG: 2 x 2
void DG(const DCOMPLEX *G,  const  DCOMPLEX *DGe,  DCOMPLEX *o_DG);

void T_eta(const int p, const int d, DCOMPLEX *o_T_eta);

void Dh(const DCOMPLEX* DG, const DCOMPLEX *T_eta, const double sigma,DCOMPLEX *o_Dh);

void Psigma_G(const int p,const int d, const double sigma, const DCOMPLEX *G, const DCOMPLEX *DGe, DCOMPLEX *o_Psigma_G);

void Psigma_h(const int p, const int d, const double sigma, const DCOMPLEX * G, const DCOMPLEX* P_sigma_G, const DCOMPLEX *T_eta,\
DCOMPLEX* o_Psigma_h);

//// Descrete

void des_DG( int p, int d, const double *a, const DCOMPLEX *W,DCOMPLEX*o_DG);

void des_Dh( const DCOMPLEX *DG, const DCOMPLEX *F,DCOMPLEX*o_Dh);

void des_Pa_h( int p, int d, const double *a, const DCOMPLEX *W, DCOMPLEX *F, DCOMPLEX *Pa_h);

/** compute gradient and loss of likelihood
* return total number of forward_iter
*
*/
int
grad_loss_cauchy_spn(  int p, int d, double  *a, double  sigma, double scale, \
  int  batch_size, double *batch, \
  double *o_grad_a, double *o_grad_sigma, double *o_loss);

#endif

max_stars_repo_name: ThayaFluss/cnl

max_issues_repo_name: ThayaFluss/cnl

max_forks_repo_name: ThayaFluss/cnl


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

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_multifit_nlin.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_blas.h>

#define SCALE 0

/*
 * This module contains an implementation of the Levenberg-Marquardt
 * algorithm for nonlinear optimization problems. This implementation
 * closely follows the following works:
 *
 * [1] H. B. Nielsen, K. Madsen, Introduction to Optimization and
 *     Data Fitting, Informatics and Mathematical Modeling,
 *     Technical University of Denmark (DTU), 2010.
 */

typedef struct
{
  gsl_matrix *A;             /* J^T J */
  gsl_matrix *A_copy;        /* copy of J^T J */
  gsl_matrix *J;             /* Jacobian J(x) */
  gsl_vector *diag;          /* D = diag(J^T J) */
  gsl_vector *rhs;           /* rhs vector = -g = -J^T f */
  gsl_vector *x_trial;       /* trial parameter vector */
  gsl_vector *f_trial;       /* trial function vector */
  gsl_vector *work;          /* workspace length p */
  long nu;                   /* nu */
  double mu;                 /* LM damping parameter mu */
  double tau;                /* initial scale factor for mu */
} lmniel_state_t;

#include "lmmisc.c"

#define LM_ONE_THIRD         (0.333333333333333)

static int lmniel_alloc (void *vstate, const size_t n, const size_t p);
static void lmniel_free(void *vstate);
static int lmniel_set(void *vstate, const gsl_vector * swts,
                      gsl_multifit_function_fdf *fdf,
                      gsl_vector *x, gsl_vector *f, gsl_vector *dx);
static int lmniel_iterate(void *vstate, const gsl_vector *swts,
                          gsl_multifit_function_fdf *fdf,
                          gsl_vector *x, gsl_vector *f, gsl_vector *dx);

static int
lmniel_alloc (void *vstate, const size_t n, const size_t p)
{
  lmniel_state_t *state = (lmniel_state_t *) vstate;

  state->A = gsl_matrix_alloc(p, p);
  if (state->A == NULL)
    {
      GSL_ERROR ("failed to allocate space for A", GSL_ENOMEM);
    }

  state->J = gsl_matrix_alloc(n, p);
  if (state->J == NULL)
    {
      GSL_ERROR ("failed to allocate space for J", GSL_ENOMEM);
    }

  state->diag = gsl_vector_alloc(p);
  if (state->diag == NULL)
    {
      GSL_ERROR ("failed to allocate space for diag", GSL_ENOMEM);
    }

  state->rhs = gsl_vector_alloc(p);
  if (state->rhs == NULL)
    {
      GSL_ERROR ("failed to allocate space for rhs", GSL_ENOMEM);
    }

  state->work = gsl_vector_alloc(p);
  if (state->work == NULL)
    {
      GSL_ERROR ("failed to allocate space for work", GSL_ENOMEM);
    }

  state->A_copy = gsl_matrix_alloc(p, p);
  if (state->A_copy == NULL)
    {
      GSL_ERROR ("failed to allocate space for A_copy", GSL_ENOMEM);
    }

  state->x_trial = gsl_vector_alloc(p);
  if (state->x_trial == NULL)
    {
      GSL_ERROR ("failed to allocate space for x_trial", GSL_ENOMEM);
    }

  state->f_trial = gsl_vector_alloc(n);
  if (state->f_trial == NULL)
    {
      GSL_ERROR ("failed to allocate space for f_trial", GSL_ENOMEM);
    }

  state->tau = 1.0e-3;

  return GSL_SUCCESS;
} /* lmniel_alloc() */

static void
lmniel_free(void *vstate)
{
  lmniel_state_t *state = (lmniel_state_t *) vstate;

  if (state->A)
    gsl_matrix_free(state->A);

  if (state->J)
    gsl_matrix_free(state->J);

  if (state->diag)
    gsl_vector_free(state->diag);

  if (state->rhs)
    gsl_vector_free(state->rhs);

  if (state->work)
    gsl_vector_free(state->work);

  if (state->A_copy)
    gsl_matrix_free(state->A_copy);

  if (state->x_trial)
    gsl_vector_free(state->x_trial);

  if (state->f_trial)
    gsl_vector_free(state->f_trial);
} /* lmniel_free() */

static int
lmniel_set(void *vstate, const gsl_vector *swts,
           gsl_multifit_function_fdf *fdf, gsl_vector *x,
           gsl_vector *f, gsl_vector *dx)
{
  int status;
  lmniel_state_t *state = (lmniel_state_t *) vstate;
  const size_t p = x->size;
  size_t i;

  /* initialize counters for function and Jacobian evaluations */
  fdf->nevalf = 0;
  fdf->nevaldf = 0;

  /* evaluate function and Jacobian at x and apply weight transform */
  status = gsl_multifit_eval_wf(fdf, x, swts, f);
  if (status)
   return status;

  if (fdf->df)
    status = gsl_multifit_eval_wdf(fdf, x, swts, state->J);
  else
    status = gsl_multifit_fdfsolver_dif_df(x, swts, fdf, f, state->J);
  if (status)
    return status;

  /* compute rhs = -J^T f */
  gsl_blas_dgemv(CblasTrans, -1.0, state->J, f, 0.0, state->rhs);

#if SCALE
  gsl_vector_set_zero(state->diag);
#else
  gsl_vector_set_all(state->diag, 1.0);
#endif

  /* set default parameters */
  state->nu = 2;

#if SCALE
  state->mu = state->tau;
#else
  /* compute mu_0 = tau * max(diag(J^T J)) */
  state->mu = -1.0;
  for (i = 0; i < p; ++i)
    {
      gsl_vector_view c = gsl_matrix_column(state->J, i);
      double result; /* (J^T J)_{ii} */

      gsl_blas_ddot(&c.vector, &c.vector, &result);
      state->mu = GSL_MAX(state->mu, result);
    }

  state->mu *= state->tau;
#endif

  return GSL_SUCCESS;
} /* lmniel_set() */

/*
lmniel_iterate()
  This function performs 1 iteration of the LM algorithm 6.18
from [1]. The algorithm is slightly modified to loop until we
find an acceptable step dx, in order to guarantee that each
function call contains a new input vector x.

Args: vstate - lm workspace
      swts   - data weights (NULL if unweighted)
      fdf    - function and Jacobian pointers
      x      - on input, current parameter vector
               on output, new parameter vector x + dx
      f      - on input, f(x)
               on output, f(x + dx)
      dx     - (output only) parameter step vector

Notes:
1) On input, the following must be initialized in state:
nu, mu, rhs, J

2) On output, the following are updated with the current iterates:
nu, mu, rhs, J

rhs needs to be set on each output, so that lmniel_gradient supplies
the correct g = J^T f
*/

static int
lmniel_iterate(void *vstate, const gsl_vector *swts,
               gsl_multifit_function_fdf *fdf, gsl_vector *x,
               gsl_vector *f, gsl_vector *dx)
{
  int status;
  lmniel_state_t *state = (lmniel_state_t *) vstate;
  gsl_matrix *J = state->J;                   /* Jacobian J(x) */
  gsl_matrix *A = state->A;                   /* J^T J */
  gsl_vector *rhs = state->rhs;               /* -g = -J^T f */
  gsl_vector *x_trial = state->x_trial;       /* trial x + dx */
  gsl_vector *f_trial = state->f_trial;       /* trial f(x + dx) */
  gsl_vector *diag = state->diag;             /* diag(D) */
  double dF;                                  /* F(x) - F(x + dx) */
  double dL;                                  /* L(0) - L(dx) */
  int foundstep = 0;                          /* found step dx */

  /* compute A = J^T J */
  status = gsl_blas_dsyrk(CblasLower, CblasTrans, 1.0, J, 0.0, A);
  if (status)
    return status;

  /* copy lower triangle to upper */
  gsl_matrix_transpose_tricpy('L', 0, A, A);

#if SCALE
  lmniel_update_diag(J, diag);
#endif

  /* loop until we find an acceptable step dx */
  while (!foundstep)
    {
      /* solve (A + mu*I) dx = g */
      status = lmniel_calc_dx(state->mu, A, rhs, dx, state);
      if (status)
        return status;

      /* compute x_trial = x + dx */
      lmniel_trial_step(x, dx, x_trial);

      /* compute f(x + dx) */
      status = gsl_multifit_eval_wf(fdf, x_trial, swts, f_trial);
      if (status)
       return status;

      /* compute dF = F(x) - F(x + dx) */
      dF = lmniel_calc_dF(f, f_trial);

      /* compute dL = L(0) - L(dx) = dx^T (mu*dx - g) */
      dL = lmniel_calc_dL(state->mu, diag, dx, rhs);

      /* check that rho = dF/dL > 0 */
      if ((dL > 0.0) && (dF >= 0.0))
        {
          /* reduction in error, step acceptable */

          double tmp;

          /* update LM parameter mu */
          tmp = 2.0 * (dF / dL) - 1.0;
          tmp = 1.0 - tmp*tmp*tmp;
          state->mu *= GSL_MAX(LM_ONE_THIRD, tmp);
          state->nu = 2;

          /* compute J <- J(x + dx) */
          if (fdf->df)
            status = gsl_multifit_eval_wdf(fdf, x_trial, swts, J);
          else
            status = gsl_multifit_fdfsolver_dif_df(x_trial, swts, fdf, f_trial, J);
          if (status)
            return status;

          /* update x <- x + dx */
          gsl_vector_memcpy(x, x_trial);

          /* update f <- f(x + dx) */
          gsl_vector_memcpy(f, f_trial);

          /* compute new rhs = -J^T f */
          gsl_blas_dgemv(CblasTrans, -1.0, J, f, 0.0, rhs);

          foundstep = 1;
        }
      else
        {
          long nu2;

          /* step did not reduce error, reject step */
          state->mu *= (double) state->nu;
          nu2 = state->nu << 1; /* 2*nu */
          if (nu2 <= state->nu)
            {
              gsl_vector_view d = gsl_matrix_diagonal(A);

              /*
               * nu has wrapped around / overflown, reset mu and nu
               * to original values and break to force another iteration
               */
              /*GSL_ERROR("nu parameter has overflown", GSL_EOVRFLW);*/
              state->nu = 2;
              state->mu = state->tau * gsl_vector_max(&d.vector);
              break;
            }
          state->nu = nu2;
        }
    } /* while (!foundstep) */

  return GSL_SUCCESS;
} /* lmniel_iterate() */

static int
lmniel_gradient(void *vstate, gsl_vector * g)
{
  lmniel_state_t *state = (lmniel_state_t *) vstate;
  gsl_vector_memcpy(g, state->rhs);
  gsl_vector_scale(g, -1.0);
  return GSL_SUCCESS;
}

static int
lmniel_jac(void *vstate, gsl_matrix * J)
{
  lmniel_state_t *state = (lmniel_state_t *) vstate;
  int s = gsl_matrix_memcpy(J, state->J);

  return s;
}

static const gsl_multifit_fdfsolver_type lmniel_type =
{
  "lmniel",
  sizeof(lmniel_state_t),
  &lmniel_alloc,
  &lmniel_set,
  &lmniel_iterate,
  &lmniel_gradient,
  &lmniel_jac,
  &lmniel_free
};

const gsl_multifit_fdfsolver_type *gsl_multifit_fdfsolver_lmniel = &lmniel_type;

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/* err/test_results.c
 *
 * 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.
 */

#include <config.h>
#include <math.h>
#include <string.h>
#include <stdio.h>
#include <stdlib.h>

#ifdef HAVE_VPRINTF
#ifdef STDC_HEADERS
#include <stdarg.h>
#else
#include <varargs.h>
#endif
#endif

#include <gsl/gsl_test.h>

static unsigned int tests = 0;
static unsigned int passed = 0;
static unsigned int failed = 0;

static unsigned int verbose = 1;

void
gsl_test (int status, const char *test_description,...)
{

  tests++;

  if (status == 0)
    {
      passed++;
      if (verbose)
	printf ("PASS: ");
    }
  else
    {
      failed++;
      if (verbose)
	printf ("FAIL: ");
    }

  if (verbose)
    {

#ifdef HAVE_VPRINTF
      va_list ap;

#ifdef STDC_HEADERS
      va_start (ap, test_description);
#else
      va_start (ap);
#endif
      vprintf (test_description, ap);
      va_end (ap);
#endif

      printf("\n");
      fflush (stdout);
    }
}

void
gsl_test_rel (double result, double expected, double relative_error,
	      const char *test_description,...)
{
  int status ;

  if (expected != 0 )
    {
      status = (fabs(result-expected)/fabs(expected) > relative_error) ;
    }
  else
    {
      status = (fabs(result) > relative_error) ;
    }

  tests++;

  if (status == 0)
    {
      passed++;
      if (verbose)
	printf ("PASS: ");
    }
  else
    {
      failed++;
      if (verbose)
	printf ("FAIL: ");

    }

  if (verbose)
    {

#ifdef HAVE_VPRINTF
      va_list ap;

#ifdef STDC_HEADERS
      va_start (ap, test_description);
#else
      va_start (ap);
#endif
      vprintf (test_description, ap);
      va_end (ap);
#endif
      if (status == 0)
	{
	  if (strlen(test_description) < 45)
	    {
	      printf(" (%g observed vs %g expected)", result, expected) ;
	    }
	  else
	    {
	      printf(" (%g obs vs %g exp)", result, expected) ;
	    }
	}
      else
	{
	  printf(" (%.18g observed vs %.18g expected)", result, expected) ;
	}

      printf ("\n") ;
      fflush (stdout);
    }
}

void
gsl_test_abs (double result, double expected, double absolute_error,
	      const char *test_description,...)
{
  int status ;

  status = fabs(result-expected) > absolute_error ;

  tests++;

  if (status == 0)
    {
      passed++;
      if (verbose)
	printf ("PASS: ");
    }
  else
    {
      failed++;
      if (verbose)
	printf ("FAIL: ");

    }

  if (verbose)
    {

#ifdef HAVE_VPRINTF
      va_list ap;

#ifdef STDC_HEADERS
      va_start (ap, test_description);
#else
      va_start (ap);
#endif
      vprintf (test_description, ap);
      va_end (ap);
#endif
      if (status == 0)
	{
	  if (strlen(test_description) < 45)
	    {
	      printf(" (%g observed vs %g expected)", result, expected) ;
	    }
	  else
	    {
	      printf(" (%g obs vs %g exp)", result, expected) ;
	    }
	}
      else
	{
	  printf(" (%.18g observed vs %.18g expected)", result, expected) ;
	}

      printf ("\n") ;
      fflush (stdout);
    }
}

void
gsl_test_factor (double result, double expected, double factor,
                 const char *test_description,...)
{
  int status;

  if (result == expected)
    {
      status = 0;
    }
  else if (expected == 0.0)
    {
      status = (result > expected || result < expected);
    }
  else
    {
      double u = result / expected;
      status = (u > factor || u < 1.0 / factor) ;
    }

  tests++;

  if (status == 0)
    {
      passed++;
      if (verbose)
	printf ("PASS: ");
    }
  else
    {
      failed++;
      if (verbose)
	printf ("FAIL: ");

    }

  if (verbose)
    {

#ifdef HAVE_VPRINTF
      va_list ap;

#ifdef STDC_HEADERS
      va_start (ap, test_description);
#else
      va_start (ap);
#endif
      vprintf (test_description, ap);
      va_end (ap);
#endif
      if (status == 0)
	{
	  if (strlen(test_description) < 45)
	    {
	      printf(" (%g observed vs %g expected)", result, expected) ;
	    }
	  else
	    {
	      printf(" (%g obs vs %g exp)", result, expected) ;
	    }
	}
      else
	{
	  printf(" (%.18g observed vs %.18g expected)", result, expected) ;
	}

      printf ("\n") ;
      fflush (stdout);
    }
}

void
gsl_test_int (int result, int expected, const char *test_description,...)
{
  int status = (result != expected) ;

  tests++;

  if (status == 0)
    {
      passed++;
      if (verbose)
	printf ("PASS: ");
    }
  else
    {
      failed++;
      if (verbose)
	printf ("FAIL: ");
    }

  if (verbose)
    {

#ifdef HAVE_VPRINTF
      va_list ap;

#ifdef STDC_HEADERS
      va_start (ap, test_description);
#else
      va_start (ap);
#endif
      vprintf (test_description, ap);
      va_end (ap);
#endif
      if (status == 0)
	{
	  printf(" (%d observed vs %d expected)", result, expected) ;
	}
      else
	{
	  printf(" (%d observed vs %d expected)", result, expected) ;
	}

      printf ("\n");
      fflush (stdout);
    }
}

void
gsl_test_str (const char * result, const char * expected,
	      const char *test_description,...)
{
  int status = strcmp(result,expected) ;

  tests++;

  if (status == 0)
    {
      passed++;
      if (verbose)
	printf ("PASS: ");
    }
  else
    {
      failed++;
      if (verbose)
	printf ("FAIL: ");
    }

  if (verbose)
    {

#ifdef HAVE_VPRINTF
      va_list ap;

#ifdef STDC_HEADERS
      va_start (ap, test_description);
#else
      va_start (ap);
#endif
      vprintf (test_description, ap);
      va_end (ap);
#endif
      if (status)
	{
	  printf(" (%s observed vs %s expected)", result, expected) ;
	}

      printf ("\n");
      fflush (stdout);
    }
}

void
gsl_test_verbose (int v)
{
  verbose = v;
}

int
gsl_test_summary (void)
{

  if (verbose && 0)		/* FIXME: turned it off, this annoys me */
    printf ("%d tests, passed %d, failed %d.\n", tests, passed, failed);

  if (failed != 0)
    {

      if (verbose && 0)		/* FIXME: turned it off, this annoys me */
	{
	  printf ("%d TEST%s FAILED.\n", failed, failed == 1 ? "" : "S");
	}
      return EXIT_FAILURE;
    }

  if (tests != passed + failed)
    {
      if (verbose)
	printf ("TEST RESULTS DO NOT ADD UP %d != %d + %d\n",
		tests, passed, failed);
      return EXIT_FAILURE;
    }

  if (passed == tests)
    {
      if (verbose && 0)		/* FIXME: turned it off, this annoys me */
	printf ("All tests passed successfully\n");
      return EXIT_SUCCESS;
    }

  return EXIT_FAILURE;
}


text:
/**
 *
 * @file core_zlag2c.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Mathieu Faverge
 * @date 2010-11-15
 * @precisions mixed zc -> ds
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 * @ingroup CORE_PLASMA_Complex64_t
 *
 *  CORE_PLASMA_zlag2c converts a PLASMA_Complex64_t matrix, A, to a
 *  PLASMA_Complex32_t matrix, B.
 *
 *******************************************************************************
 *
 * @param[in] m
 *          The number of rows of the matrices A and B. m >= 0.
 *
 * @param[in] n
 *          The number of columns of the matrices A and B. n >= 0.
 *
 * @param[in] A
 *          The PLASMA_Complex64_t m-by-n matrix to convert.
 *
 * @param[in] lda
 *          The leading dimension of the array A. lda >= max(1,m).
 *
 * @param[out] B
 *          The PLASMA_Complex32_t m-by-n matrix to convert.
 *
 * @param[in] ldb
 *          The leading dimension of the array B. ldb >= max(1,m).
 *
 * @param[out] info
 *         - 0 on successful exit.
 *         - 1 if an entry of the matrix A is greater than the SINGLE
 *         PRECISION overflow threshold, in this case, the content
 *         of B in exit is unspecified.
 *
 ******************************************************************************/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_zlag2c = PCORE_zlag2c
#define CORE_zlag2c PCORE_zlag2c
#endif
void CORE_zlag2c(int m, int n,
                 const PLASMA_Complex64_t *A, int lda,
                 PLASMA_Complex32_t *B, int ldb, int *info)
{
    *info = LAPACKE_zlag2c_work(LAPACK_COL_MAJOR, m, n, A, lda, B, ldb);
}

/***************************************************************************//**
 *
 * @ingroup CORE_PLASMA_Complex64_t
 *
 *  CORE_PLASMA_clag2z converts a PLASMA_Complex32_t matrix, A, to a
 *  PLASMA_Complex64_t matrix, B.
 *
 *******************************************************************************
 *
 * @param[in] m
 *          The number of rows of the matrices A and B. m >= 0.
 *
 * @param[in] n
 *          The number of columns of the matrices A and B. n >= 0.
 *
 * @param[in] A
 *          The PLASMA_Complex32_t m-by-n matrix to convert.
 *
 * @param[in] lda
 *          The leading dimension of the array A. lda >= max(1,m).
 *
 * @param[out] B
 *          The PLASMA_Complex64_t m-by-n matrix to convert.
 *
 * @param[in] ldb
 *          The leading dimension of the array B. ldb >= max(1,m).
 *
 ******************************************************************************/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_clag2z = PCORE_clag2z
#define CORE_clag2z PCORE_clag2z
#endif
void CORE_clag2z(int m, int n,
                 const PLASMA_Complex32_t *A, int lda,
                 PLASMA_Complex64_t *B, int ldb)
{
    LAPACKE_clag2z_work(LAPACK_COL_MAJOR, m, n, A, lda, B, ldb);
}