void preprocess1(int panelSz, int D, float *XX, float*X, float*Y, float lamda, float*Z, float*B){
    // panelSz: the number of points to process in each round
    int i,j;
    // step 1: compute all X[i]'*X[j] (i>j)
    for (i=panelSz-1;i>0;i--)
        for (j=i-1;j>=0;j--){
            XX[i*panelSz+j] = cblas_sdot(D, &(X[i*D]), 1, &(X[j*D]), 1);
         // printf("XX[%d]=%8.4f, X[%d]=%8.4f, X[%d]=%8.4f\n", i*panelSz+j, XX[i*panelSz+j], i*D, X[i*D], j*D, X[j*D]);
        }

	
    // step 2: compute all Z vectors
    // Z0=lamda*X[0], B=lamda*X[0]*Y[0]
    cblas_scopy(D, X, 1, Z, 1);  

    cblas_sscal(D, lamda, Z, 1);
    float alpha=lamda*Y[0];
    cblas_scopy(D, X, 1, B, 1);
    cblas_sscal(D, alpha, B, 1);
    for (i=1; i<panelSz;i++){
        cblas_scopy(D, &(X[i*D]), 1, &(Z[i*D]),1);
        // Z[i] = lamda*(X[i] - sum_{j<i} XX[i,j]*Z[j]);
        for (j=0;j<i;j++){
            cblas_saxpy(D, -1*XX[i*panelSz+j], &(Z[j*D]), 1, &(Z[i*D]), 1);
        }
        cblas_sscal(D, lamda, &(Z[i*D]), 1);
        // B = lamda*(Y[i] - X[i]*B) X[i] + B;
        float alpha = lamda*(Y[i]-cblas_sdot(D, &(X[i*D]), 1, B, 1));
        cblas_saxpy(D, alpha, &(X[i*D]), 1, B, 1);
    }
}

/* Trains a network by presenting an example and 
 * adjusts the weights by stochastic gradient 
 * descent to reduce a squared hinge loss
 */
void train(nnet_t* n, sparse_t* v, int target){
    int i;
    /* Forward pass */
    cblas_scopy(n->hidden,n->b1,1,n->a1,1);
    for(i=0; i<v->nz; i++){
        cblas_saxpy(n->hidden, v->x[i], n->W1[v->idx[i]], 1, n->a1, 1);
    }
    activation(n->a1,n->x1,n->g1,n->hidden);
    n->a2 = n->b2 + cblas_sdot(n->hidden, n->W2, 1, n->x1, 1);
    activation(&n->a2,&n->x2,&n->g2,1);
    if(target*n->x2 > 1)
        /* Hinge loss, no error -> no need to backpropagate */
        return;
    /* Backward pass */
    n->d2 = (target-n->x2)*n->g2;
    cblas_scopy(n->hidden,n->W2,1,n->d1,1);
    for(i=0; i<n->hidden; i++)
        n->d1[i] *= n->d2*n->g1[i];
    n->b2 += n->eta*n->d2;
    cblas_saxpy(n->hidden, n->eta*n->d2, n->x1, 1, n->W2, 1);
    cblas_saxpy(n->hidden, n->eta, n->d1, 1, n->b1, 1);
    /* Sparse inputs imply sparse gradients.
     * This update saves a lot of computation
     * compared to general purpose neural net
     * implementations.
     */
    for(i=0; i<v->nz; i++){
        cblas_saxpy(n->hidden, n->eta*v->x[i], n->d1, 1, n->W1[v->idx[i]], 1);
    }
}

    void DecoderBinaural::process(const double* inputs, double* outputs)
	{
        m_decoder->process(inputs, m_channels_vector_double);

        --m_index;
        m_channels_inputs_left[0][m_index] = m_channels_vector_double[0];
        outputs[1] = outputs[0] = cblas_sdot(m_impulses_size, m_channels_inputs_left[0]+m_index, 1, m_impulses_vector[0], 1);
        for(int i = 1; i < m_number_of_virtual_channels; i++)
        {
            m_channels_inputs_left[i][m_index] = m_channels_vector_double[i];
            m_channels_inputs_right[i][m_index] = m_channels_vector_double[i];
            outputs[0] += cblas_sdot(m_impulses_size, m_channels_inputs_left[i]+m_index, 1, m_impulses_vector[m_number_of_virtual_channels - i], 1);
            outputs[1] += cblas_sdot(m_impulses_size, m_channels_inputs_right[i]+m_index, 1, m_impulses_vector[i], 1);
        }
        if(m_index <= 0)
        {
            m_index = m_impulses_size;
            cblas_scopy(m_impulses_size, m_channels_inputs_left[0], 1, m_channels_inputs_left[0]+m_impulses_size, 1);
            for(int i = 1; i < m_number_of_virtual_channels; i++)
            {
                cblas_scopy(m_impulses_size, m_channels_inputs_left[i], 1, m_channels_inputs_left[i]+m_impulses_size, 1);
                cblas_scopy(m_impulses_size, m_channels_inputs_right[i], 1, m_channels_inputs_right[i]+m_impulses_size, 1);
            }
        }
    }

void PolarLines::process(float* vector)
{
    cblas_saxpy(m_number_of_sources * 2, 1., m_values_step, 1, m_values_old, 1);
    if(m_counter++ >= m_ramp)
    {
        cblas_scopy(m_number_of_sources * 2, m_values_new, 1, m_values_old, 1);
        memset(m_values_step, 0, sizeof(float) * m_number_of_sources * 2);
        m_counter    = 0;
    }
    cblas_scopy(m_number_of_sources * 2, m_values_old, 1, vector, 1);
}

void hoa_optim_perform(t_hoa_optim *x, t_object *dsp64, float **ins, long numins, float **outs, long numouts, long sampleframes, long flags, void *userparam)
{
	for(int i = 0; i < numins; i++)
    {
        cblas_scopy(sampleframes, ins[i], 1, x->f_ins+i, numins);
    }
	for(int i = 0; i < sampleframes; i++)
    {
        x->f_optim->process(x->f_ins + numins * i, x->f_outs + numouts * i);
    }
    for(int i = 0; i < numouts; i++)
    {
        cblas_scopy(sampleframes, x->f_outs+i, numouts, outs[i], 1);
    }
}

void hoa_wider_3D_perform_offset(t_hoa_wider_3D *x, t_object *dsp, float **ins, long numins, float **outs, long numouts, long sampleframes, long f,void *up)
{
	for(int i = 0; i < numins - 1; i++)
    {
        cblas_scopy(sampleframes, ins[i], 1, x->f_ins+i, numins - 1);
    }
	for(int i = 0; i < sampleframes; i++)
    {
        x->f_wider->process(x->f_ins + (numins - 1) * i, x->f_outs + numouts * i);
    }
    for(int i = 0; i < numouts; i++)
    {
        cblas_scopy(sampleframes, x->f_outs+i, numouts, outs[i], 1);
    }
}

void SENNA_nn_temporal_convolution(float *output, int output_frame_size,
                                   float *weights, float *biases, float *input,
                                   int input_frame_size, int n_frames,
                                   int k_w) {
#ifdef USE_BLAS
  if (k_w == 1) {
    if (biases) {
      int t;
      for (t = 0; t < n_frames; t++)
        cblas_scopy(output_frame_size, biases, 1,
                    output + t * output_frame_size, 1);
    }
    cblas_sgemm(CblasColMajor, CblasTrans, CblasNoTrans, output_frame_size,
                n_frames, input_frame_size, 1.0, weights, input_frame_size,
                input, input_frame_size, (biases ? 1.0 : 0.0), output,
                output_frame_size);
  } else
#endif
  {
    int t;

    for (t = 0; t < n_frames - k_w + 1; t++)
      SENNA_nn_linear(output + t * output_frame_size, output_frame_size,
                      weights, biases, input + t * input_frame_size,
                      input_frame_size * k_w);
  }
}

void THBlas_copy(long size, const real *x, long xStride, real *y, long yStride)
{
  if(size == 1)
  {
    xStride = 1;
    yStride = 1;
  }

#if USE_CBLAS
  if( (size < INT_MAX) && (xStride < INT_MAX) && (yStride < INT_MAX) )
  {
#ifdef USE_DOUBLE
    cblas_dcopy(size, x, xStride, y, yStride);
#else
    cblas_scopy(size, x, xStride, y, yStride);
#endif
    return;
  }
#endif
  {
    long i;
    for(i = 0; i < size; i++)
      y[i*yStride] = x[i*xStride];
  }
}

void THBlas_(copy)(long n, real *x, long incx, real *y, long incy)
{
  if(n == 1)
  {
    incx = 1;
    incy = 1;
  }

#if defined(USE_BLAS) && (defined(TH_REAL_IS_DOUBLE) || defined(TH_REAL_IS_FLOAT))
  if( (n <= INT_MAX) && (incx <= INT_MAX) && (incy <= INT_MAX) )
  {
    int i_n = (int)n;
    int i_incx = (int)incx;
    int i_incy = (int)incy;

#if defined(TH_REAL_IS_DOUBLE)
    cblas_dcopy(i_n, x, i_incx, y, i_incy);
#else
    cblas_scopy(i_n, x, i_incx, y, i_incy);
#endif
    return;
  }
#endif
  {
    long i;
    for(i = 0; i < n; i++)
      y[i*incy] = x[i*incx];
  }
}

JNIEXPORT void JNICALL Java_uncomplicate_neanderthal_CBLAS_scopy
(JNIEnv *env, jclass clazz, jint N,
 jobject X, jint offsetX, jint incX,
 jobject Y, jint offsetY, jint incY) {

    float *cX = (float *) (*env)->GetDirectBufferAddress(env, X);
    float *cY = (float *) (*env)->GetDirectBufferAddress(env, Y);
    cblas_scopy(N, cX + offsetX, incX, cY + offsetY, incY);
};

  void channel_to_img(const float* const src, const int height, const int width,
                      const int channel, float* const dst)
  {
    const int sz = height * width;

    for (int c = 0; c < channel; ++c)
    {
      cblas_scopy(sz, src + c * sz, 1, dst + c, channel);
    }
  }

void CsoundObject_readCallback(void *inRefCon, UInt32 inNumberFrames, Float32 **audioData)
{
    CsoundObject *self = (CsoundObject *) inRefCon;
    MYFLT *csoundOut = csoundGetSpout(self->csound);
    csoundPerformKsmps(self->csound);
    
    for (size_t channel = 0; channel < self->channelsCount; ++channel) {
        
        cblas_scopy((SInt32)inNumberFrames, &csoundOut[channel], 2, audioData[channel], 1);
    }
}

void HoaToolsAudioProcessor::processBlock(AudioSampleBuffer& buffer, MidiBuffer& midiMessages)
{
    int i;
    int numins = getNumInputChannels();
    int numouts = getNumOutputChannels();
    int nharmo = NHARMO;
    int vectorsize = buffer.getNumSamples();

    for(i = 0; i < numins; i++)
    {
        cblas_scopy(vectorsize, buffer.getReadPointer(i), 1, m_input_vector+i, numins);
        m_lines->setRadius(i, m_sources->sourceGetRadius(i));
        m_lines->setAzimuth(i, m_sources->sourceGetAzimuth(i));
        if(m_sources->sourceGetExistence(i))
            m_map->setMute(i, 0);
        else
            m_map->setMute(i, 1);
    }
    for(; i < 16; i++)
    {
        m_map->setMute(i, 1);
    }
    for(i = 0; i < vectorsize; i++)
    {
        m_lines->process(m_lines_vector);
        for(int j = 0; j < numins; j++)
            m_map->setRadius(j, m_lines_vector[j]);
        for(int j = 0; j < numins; j++)
            m_map->setAzimuth(j, m_lines_vector[j+numins]);

        m_map->process(m_input_vector+ numins * i, m_harmo_vector + nharmo * i);
        m_optim->process(m_harmo_vector + nharmo * i, m_harmo_vector + nharmo * i);
        m_decoder->process(m_harmo_vector + nharmo * i, m_output_vector + numouts * i);
        m_meter->process(m_output_vector + numouts * i);
    }

    for(i = 0; i < numouts; i++)
    {
        cblas_scopy(vectorsize, m_output_vector+i, numouts, buffer.getWritePointer(i), 1);
    }
}

/* Given an input vector v, compute the output of the network. */
float value(nnet_t* n, sparse_t* v){
    int i;
    cblas_scopy(n->hidden,n->b1,1,n->a1,1);
    for(i=0; i<v->nz; i++){
        cblas_saxpy(n->hidden, v->x[i], n->W1[v->idx[i]], 1, n->a1, 1);
    }
    activation(n->a1,n->x1,n->g1,n->hidden);
    n->a2 = n->b2;
    n->a2 += cblas_sdot(n->hidden, n->W2, 1, n->x1, 1);
    activation(&n->a2,&n->x2,&n->g2,1);
    return n->x2;
}

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

int main(void)
{
	const int p = 10;
	float *x = malloc(p * sizeof *x);
	float *y = malloc(p * sizeof *y);
	float *z = malloc(p * sizeof *y);
	float *orig = malloc(p * sizeof *orig);
	int i;

	if (!(x && y && orig && z)) {
		perror("malloc");
		exit(1);
	}

	for (i = 0; i < p; ++i) {
		y[i] = i;
		x[i] = 2 * i;
		z[i] = 0;
	}
	cblas_scopy(p, y, 1, orig, 1);
	cblas_saxpy(p, -1, x, 1, y, 1);

	/* There's GEMV, but you have to put one vector into the
	 * diagonal of a matrix because... uh...
	 * Somethin'?
	 *
	 * But WAIT!  It's not stupid because there are compact
	 * ways to represent banded matrices, and a diagonal matrix
	 * can be represented by a single row, kind of like it
	 * was a ... VECTOR!?!?!?
	 *
	 * <neo>Whoa.</neo>
	 */
	cblas_sgbmv(CblasRowMajor,
		    CblasNoTrans, /* Don't transpose A */
		    p, p,
		    0,		/* bands below the diagonal */
		    0,		/* bands above the diagonal */
		    1,		/* alpha */
		    y,		/* some data, at last! */
		    1,		/* LDA is 1st dim of A */
		    y,
		    1,
		    1,
		    z,
		    1);
	for (i = 0; i < p; ++i) {
		printf("y(%f) - x(%f) = %f\n", orig[i], x[i], y[i]);
		printf("(y-x)^2 = %f\n", z[i]);
	}
	return 0;
}

extern "C" void
magma_strdtype2cbHLsym_withQ_v2(magma_int_t n, magma_int_t nb, 
                                float *A, magma_int_t lda, 
                                float *V, magma_int_t ldv, 
                                float *TAU,
                                magma_int_t st, magma_int_t ed, 
                                magma_int_t sweep, magma_int_t Vblksiz, 
                                float *work) 
{

    /*
     WORK (workspace) float real array, dimension NB
    */

    magma_int_t ione = 1;
    magma_int_t vpos, taupos;

    float conjtmp;

    float c_one = MAGMA_S_ONE;

    magma_int_t ldx = lda-1;
    magma_int_t len = ed - st + 1;
    magma_int_t lem = min(ed+nb, n) - ed;

    if(lem>0){
        magma_bulge_findVTAUpos(n, nb, Vblksiz, sweep-1, st-1, ldv, &vpos, &taupos);
        /* apply remaining right coming from the top block */
        lapackf77_slarfx("R", &lem, &len, V(vpos), TAU(taupos), A(ed+1, st), &ldx, work);
    }
    if(lem>1){
        magma_bulge_findVTAUpos(n, nb, Vblksiz, sweep-1, ed, ldv, &vpos, &taupos);

        /* remove the first column of the created bulge */
        *V(vpos)  = c_one;
        //memcpy(V(vpos+1), A(ed+2, st), (lem-1)*sizeof(float));
        cblas_scopy(lem-1, A(ed+2, st), ione, V(vpos+1), ione);
        memset(A(ed+2, st),0,(lem-1)*sizeof(float));

        /* Eliminate the col at st */
        lapackf77_slarfg( &lem, A(ed+1, st), V(vpos+1), &ione, TAU(taupos) );

        /* apply left on A(J1:J2,st+1:ed) */
        len = len-1; /* because we start at col st+1 instead of st. col st is the col that has been removed;*/
        conjtmp = MAGMA_S_CNJG(*TAU(taupos));
        lapackf77_slarfx("L", &lem, &len, V(vpos),  &conjtmp, A(ed+1, st+1), &ldx, work);
    }

}

    void DecoderBinaural::process(const float* const* inputs, float** outputs)
	{
        unsigned int i;
        for(i = 0; i < m_number_of_harmonics; i++)
        {
            cblas_scopy(m_vector_size, inputs[i], 1,  m_input_matrix+i*m_vector_size, 1);
        }
        
        cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, (m_impulses_size * 2), m_vector_size, m_number_of_harmonics, 1.,
                    m_impulses_matrix, m_number_of_harmonics,
                    m_input_matrix,  m_vector_size,
                    0., m_result_matrix,  m_vector_size);
        
        for(i = 0; i < m_vector_size; i++)
        {
            cblas_saxpy(m_impulses_size, 1.f, m_result_matrix + i, m_vector_size, m_linear_vector_left + i, 1);
            outputs[0][i] = m_linear_vector_left[i];
        }
        
        for(i = 0; i < m_vector_size; i++)
        {
            cblas_saxpy(m_impulses_size, 1.f, m_result_matrix + i + m_vector_size * m_impulses_size, m_vector_size, m_linear_vector_right + i, 1);
            outputs[1][i] = m_linear_vector_right[i];
        }
        
        cblas_scopy(m_impulses_size-1, m_linear_vector_left+m_vector_size, 1, m_linear_vector_left, 1);
        cblas_scopy(m_impulses_size-1, m_linear_vector_right+m_vector_size, 1, m_linear_vector_right, 1);
        
#ifdef __APPLE__
        vDSP_vclr(m_linear_vector_left + m_impulses_size - 1, 1, m_vector_size);
        vDSP_vclr(m_linear_vector_right + m_impulses_size - 1, 1, m_vector_size);
#else
        memset(m_linear_vector_left + m_impulses_size - 1, 0, m_vector_size * sizeof(float));
        memset(m_linear_vector_right + m_impulses_size - 1, 0, m_vector_size * sizeof(float));
#endif
	}

void hoa_map_3D_tilde_perform(t_hoa_map_3D_tilde *x, t_object *dsp64, float **ins, long numins, float **outs, long numouts, long sampleframes, long flags, void *userparam)
{
    for(int i = 0; i < sampleframes; i++)
    {
		x->f_lines->process(x->f_lines_vector);
		x->f_map->setRadius(0, x->f_lines_vector[0]);
		x->f_map->setAzimuth(0, x->f_lines_vector[1]);
        x->f_map->setElevation(0, x->f_lines_vector[2]);
        x->f_map->process(&ins[0][i], x->f_sig_outs + numouts * i);
    }
    for(int i = 0; i < numouts; i++)
    {
        cblas_scopy(sampleframes, x->f_sig_outs+i, numouts, outs[i], 1);
    }
}

void hoa_map_3D_tilde_perform_multisources(t_hoa_map_3D_tilde *x, t_object *dsp64, float **ins, long numins, float **outs, long numouts, long sampleframes, long flags, void *userparam)
{
	int nsources = x->f_map->getNumberOfSources();
    for(int i = 0; i < numins; i++)
    {
        cblas_scopy(sampleframes, ins[i], 1, x->f_sig_ins+i, numins);
    }
    for(int i = 0; i < sampleframes; i++)
    {
        x->f_lines->process(x->f_lines_vector);
		for(int j = 0; j < nsources; j++)
			x->f_map->setRadius(j, x->f_lines_vector[j]);
        for(int j = 0; j < nsources; j++)
			x->f_map->setAzimuth(j, x->f_lines_vector[j + nsources]);
        for(int j = 0; j < nsources; j++)
			x->f_map->setElevation(j, x->f_lines_vector[j + nsources * 2]);
        
        x->f_map->process(x->f_sig_ins + numins * i, x->f_sig_outs + numouts * i);
    }
    for(int i = 0; i < numouts; i++)
    {
        cblas_scopy(sampleframes, x->f_sig_outs+i, numouts, outs[i], 1);
    }
}