text:
/* File description : Implementation of neural networks
 * Author : Bruno Faria
 * Year : 2014
 * Work place : University of Aveiro
*/
#include <time.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "neural_network.h"

//#include "mkl.h"
//#include <atlas_enum.h>
//#include <lapacke.h>
//#include <lapacke.h>
const char * nn_error_strings[19] = { "API succeed",
								  	  "API could not allocate memory for network object",
								  	  "API could not evaluate network object",
								  	  "API could not rescale data",
								  	  "neural network not initialized",
								  	  "API bad input",
								  	  "API weight initialization",
 								  	  "API could not duplicate network",
								  	  "API could not copy data from src to dest, networks are different",
								  	  "API could not copy data :s",
								  	  "API could not solve system of equations (lapack least squares)",
								  	  "API lapack error",
								  	  "API could not free memory (what?)",
								  	  "API train",
								  	  "API train 2",
								  	  "API train 3",
								  	  "API could not duplicate network",
								  	  "API could not get inverse",
								  	  "API could not compute J, G, and H"
								  };

#define log2(x) log(x)/0.30102999566
/* Random number generator vars */
unsigned long mt_state[MT_N];
unsigned long mt_constants[2];
int mt_state_idx;
bool rand_state = false;

/* global stack error */
int32_t NN_API_ERROR = 0;

/* =======================================
 * = Allocate and destroy neural network =
 * =======================================
*/
net_t * create_nn(int32_t num_layers, int32_t * layers_conf, int32_t * layers_func )
{
   int32_t i = 0, j = 0;
   net_t * nn = NULL;

   C_CHECK_CONDITION( num_layers == 0 || num_layers > 200 || num_layers < 1, NN_API_BAD_INPUT);

   /* allocate a neural network object */
   C_SAFE_CALL( nn = (net_t *)mem_alloc(sizeof(net_t), true ) );

   /* allocate the layers field pointers */
   C_SAFE_CALL( nn->layer = (net_layer_t **)mem_alloc(num_layers * sizeof(net_layer_t *), true) );

   /* set num layers */
   nn -> num_layers = num_layers;
   nn -> init_state = false;

   /* allocate layers vector */
   C_SAFE_CALL( nn->layer[0] = (net_layer_t *)mem_alloc(num_layers * sizeof(net_layer_t), true) );

   /* generate layers pointers */
   for(i = 1; i < num_layers; i++) { nn->layer[i] = nn->layer[0] + i; }

   /* allocate space for the layers */
   for(i = 0; i < num_layers; i++)
   {

   		C_CHECK_CONDITION( layers_conf[i] <= 0 || layers_conf[i] > 400, NN_API_BAD_INPUT );

   		/* insert number of inputs and outputs of the layer */
   		if(i == 0)
   		{
   			C_CHECK_CONDITION( layers_func[i] !=  NN_NONE  || layers_func[i] > net_func_nelem, NN_API_BAD_INPUT );
   			/* input layer */
   			nn->layer[i]->num_inputs = 0;
   			nn->layer[i]->num_outputs = layers_conf[i];
            nn->layer[i]->eval_funct = NN_NONE;
            nn->layer[i]->weights = NULL;
            nn->layer[i]->bias = NULL;
			C_SAFE_CALL( nn->layer[i]->results = (double*)mem_alloc( layers_conf[i] * sizeof(double), true ));

		    /* allocate memory for input and output mapping */
   			C_SAFE_CALL( nn->xmin = (double *)mem_alloc(layers_conf[i] * sizeof(double),true ));
   			C_SAFE_CALL( nn->xmax = (double *)mem_alloc(layers_conf[i] * sizeof(double),true ));

   		}else{

   			C_CHECK_CONDITION( layers_func[i] ==  NN_NONE  || layers_func[i] > net_func_nelem, NN_API_BAD_INPUT );

   			/* allocate memory for input and output mapping */
   			C_SAFE_CALL( nn->ymin = (double *)mem_alloc(layers_conf[i] * sizeof(double),true ));
   			C_SAFE_CALL( nn->ymax = (double *)mem_alloc(layers_conf[i] * sizeof(double),true ));

   			/* remaining layers */
   			nn->layer[i]->num_inputs = layers_conf[i-1];
   			nn->layer[i]->num_outputs = layers_conf[i];
        	nn->layer[i]->eval_funct = layers_func[i];

            C_SAFE_CALL( nn->layer[i]->weights 	 	 = (double**)mem_alloc( layers_conf[i] 				  	* sizeof(double 	*), true ) );
            C_SAFE_CALL( nn->layer[i]->results 	 	 = (double *)mem_alloc( layers_conf[i] 					* sizeof(double		 ), true ) );
            C_SAFE_CALL( nn->layer[i]->bias 		 = (double *)mem_alloc( layers_conf[i] 					* sizeof(double		 ), true ) );
            C_SAFE_CALL( nn->layer[i]->weights[0] 	 = (double *)mem_alloc( layers_conf[i-1] * layers_conf[i] * sizeof(double		 ), true ) );

            for(j = 1; j < layers_conf[i]; j++)
            {
            	nn->layer[i]->weights[j] = nn->layer[i]->weights[0] + j * layers_conf[i-1];
            }
       }
   }
   nn -> init_state = true;
   nn -> use_mapping = true;
   return nn;

error:
    /* something went wrong so clean memory and exit */
    destroy_nn( nn );

    return NULL;
}

int32_t destroy_nn(net_t * nn)
{
    int32_t i = 0;
    net_layer_t * layer;

    C_CHECK_CONDITION( nn == NULL, NN_API_BAD_INPUT );
    C_CHECK_CONDITION( nn->num_layers == 0, NN_API_BAD_INPUT );

    /* first deallocate the layers */
    if(nn->layer != NULL)
    {
    	if(nn->layer[0] != NULL)
    	{
			layer = nn->layer[0];
			mem_free(layer->results);
    		for(i = 1; i < nn->num_layers; i++ )
    		{
    			layer = nn->layer[i];
    			if(layer->weights != NULL)
    			{
    				mem_free( layer->weights[0] );
        			mem_free( layer->weights 	);
        		}
        		mem_free( layer->results 	 );
       	 		mem_free( layer->bias    	 );
    		}
    		/* now free the main object */
   			mem_free( nn->layer[0] );
   		}
   		mem_free( nn->layer );
    }

    mem_free( nn->xmin );
    mem_free( nn->xmax );
    mem_free( nn->ymin );
    mem_free( nn->ymax );

    mem_free( nn );

    return 0;
error:
    /* what? cannot deallocate memory? why? */
    return -1;
}

void * mem_alloc(int32_t size, bool zero_set)
{
     void * p = malloc(size);

     if( p == NULL){ NN_API_ERROR = NN_API_NETWORK_ALLOC_ERROR; return NULL;}

     if(zero_set == true)
     {
        memset(p, 0, size);
     }
     NN_API_ERROR = NN_API_SUCCESS;
     return p;
}

void mem_free(void * ptr)
{
	C_CHECK_CONDITION(ptr == NULL, NN_API_BAD_INPUT);
    free( ptr );
    ptr = NULL;
    NN_API_ERROR = NN_API_SUCCESS;
	return;
error:
	NN_API_ERROR = NN_API_FREE;
	return;
}

/* get the neural network response for a given input */
int32_t get_nn_response(net_t * nn)
{
	int32_t i = 0, n = 0, l = 0;
	double w_sum = 0.0f;

	net_layer_t * layer = NULL, * prev_layer = NULL;

	C_CHECK_CONDITION( nn == NULL, 				  NN_API_BAD_INPUT 		);
	C_CHECK_CONDITION( nn -> init_state == false, NN_API_BAD_INPUT 		);
	C_CHECK_CONDITION( nn -> layer == NULL,		  NN_API_MALFORMED_NN	);

	memcpy(nn->layer[0]->results, nn->in_layer,nn->layer[0]->num_outputs * sizeof(double));

	/* first layer hols the data to be evaluated */
	if(nn->use_mapping){
		C_SAFE_CALL( transform_in_out_data(nn, nn->layer[0]->results, nn->layer[0]->num_outputs, false) );
	}

	for(l = 1; l < nn->num_layers; l++)
	{
		C_CHECK_CONDITION( nn -> layer[l] 	== NULL,	NN_API_MALFORMED_NN	);
		C_CHECK_CONDITION( nn -> layer[l-1] == NULL,	NN_API_MALFORMED_NN	);

		/* remap layer (to increase readibility) */
		layer = nn->layer[l];
		prev_layer = nn->layer[l-1];

		C_CHECK_CONDITION( layer -> bias 		== NULL,	NN_API_MALFORMED_NN	);
		C_CHECK_CONDITION( layer -> weights 	== NULL,	NN_API_MALFORMED_NN	);
		C_CHECK_CONDITION( layer -> weights[0] 	== NULL,	NN_API_MALFORMED_NN	);

		for(n = 0; n < layer->num_outputs; n++)
		{

			w_sum = layer->bias[n];

			for( i = 0; i < layer->num_inputs; i++)
				w_sum += prev_layer->results[i] * layer->weights[n][i];

			// apply layer function
			switch( layer->eval_funct )
			{
				case NN_LOGSIG: { w_sum =  1 / (1 + exp(-w_sum)); 			break; }
				case NN_TANSIG: { w_sum =  2 / (1 + exp(-2 * w_sum)) - 1;	break; }
				default: 		{ 											break; }
			}

			layer->results[n] = w_sum;
		}
	}

	/* transform the data and output it to the outer layer */
	layer = nn->layer[nn->num_layers-1];
	memcpy(nn->out_layer, layer->results, layer->num_outputs * sizeof(double));
	if(nn->use_mapping)
	{
		C_SAFE_CALL( transform_in_out_data(nn, nn->out_layer, layer->num_outputs, true) );
	}

	NN_API_ERROR = NN_API_SUCCESS;
	return 0;
error:
	NN_API_ERROR = NN_API_RESCALE_DATA;
	return -1;
}

inline double map_min_max( double x, double xmin, double xmax, double ymin, double ymax, bool reverse )
{
	if(reverse)
	{
		return (ymax-ymin)*(x+1.0f)/2.0f + ymin;
	}else{
		return 2.0f * (x - xmin)/(xmax-xmin) - 1.0f;
	}
}

int32_t transform_in_out_data(net_t * nn, double * in_out, int32_t size, bool Mode)
{
	int32_t i = 0;

	C_CHECK_CONDITION( size == 0, NN_API_BAD_INPUT);

	for(i = 0; i < size; i++){
		if(Mode){
			in_out[i] = map_min_max(in_out[i], 0, 0, nn->ymin[i], nn->ymax[i], Mode);
		}else{
			in_out[i] = map_min_max(in_out[i], nn->xmin[i], nn->xmax[i], 0, 0, Mode);
		}
	}
	NN_API_ERROR = NN_API_SUCCESS;
	return 0;

error:
	NN_API_ERROR = NN_API_RESCALE_DATA;
	return -1;
}

/*
 * Compute stuff
 */

double compute_JGH(net_t * nn, double **J, double *G, double **H, double ** samples, double ** targets, double * errors, double * nn_out, double * nn_temp_out,
 int32_t * perm, int32_t num_samples, int32_t num_outputs, int32_t G_length)
{
	double sse = 0;
	int32_t i = 0, j = 0, k = 0, l = 0, s = 0, ct = 0;
	net_layer_t * layer = NULL;
	double weight_s = 0;

	/* verify conditions */
	C_CHECK_CONDITION( nn 			== NULL, 				  NN_API_BAD_INPUT 	);
	C_CHECK_CONDITION( J  			== NULL, 				  NN_API_BAD_INPUT 	);
	C_CHECK_CONDITION( G  			== NULL,				  NN_API_BAD_INPUT	);
	C_CHECK_CONDITION( samples  	== NULL, 				  NN_API_BAD_INPUT 	);
	C_CHECK_CONDITION( targets  	== NULL,				  NN_API_BAD_INPUT	);
	C_CHECK_CONDITION( errors  		== NULL, 				  NN_API_BAD_INPUT 	);
	C_CHECK_CONDITION( nn_out  		== NULL,				  NN_API_BAD_INPUT	);
	C_CHECK_CONDITION( nn_temp_out  == NULL, 				  NN_API_BAD_INPUT 	);
	C_CHECK_CONDITION( num_samples  == 0, 					  NN_API_BAD_INPUT 	);
	C_CHECK_CONDITION( num_outputs  == 0,					  NN_API_BAD_INPUT	);
	C_CHECK_CONDITION( G_length  	== 0,					  NN_API_BAD_INPUT	);

	/* compute sample errors and jacobian */
	for(s = 0; s < num_samples; s++)
	{
		ct = 0;

		/* short circuit the pointers */
		nn->in_layer = samples[perm[s]]; nn->out_layer = nn_out;

		/* first evaluate the network for this sample without any perturbation*/
		C_SAFE_CALL( get_nn_response( nn ) );

		/* compute errors and sse */
		for(i = 0; i < num_outputs; i++)
		{
			errors[i*num_samples + s] = targets[perm[s]][i] - nn_out[i];
			sse += errors[i*num_samples + s]*errors[i*num_samples + s];
		}

		nn->out_layer = nn_temp_out;

		/* the lines of the jacobian for this sample are computed here */
		for(l = 1; l < nn->num_layers; l++)
		{
			layer = nn->layer[l];

			for(i = 0; i < layer->num_outputs; i++)
			{
				/* weights perturbation */
				for(j = 0; j < layer->num_inputs; j++)
				{

					weight_s = layer->weights[i][j];

					/* apply perturbation parameter */
					layer->weights[i][j] += (PERTURBATION * ( 1 + abs(weight_s)));

					/* get neural network response */
					C_SAFE_CALL( get_nn_response( nn ) );

					/* restore weight */
					layer->weights[i][j] = weight_s;

					/* compute derivative */
					for(k = 0; k < num_outputs; k++){
						J[k*num_samples+s][ct] = (nn_out[k] - nn_temp_out[k]) / (PERTURBATION * ( 1 + abs(weight_s)));
					}
					ct++;
				}

				/* bias perturbation */
				weight_s = layer->bias[i];

				/* apply perturbation parameter */
				layer->bias[i] += (PERTURBATION * ( 1 + abs(weight_s)));

				/* get neural network response */
				C_SAFE_CALL( get_nn_response( nn ) );

				/* restore weight */
				layer->bias[i] = weight_s;

				/* compute derivative */
				for(k = 0; k < num_outputs; k++)
					J[k*num_samples+s][ct] = (nn_out[k] - nn_temp_out[k]) / (PERTURBATION * ( 1 + abs(weight_s)));
				ct++;
			}

		}
	}

	/* compute gradient and aproximation to the hessian matrix */
	memset(G, 0, G_length * sizeof(double));
	memset(H[0], 0, G_length * G_length * sizeof(double));

	for(i = 0; i < num_samples * num_outputs; i++ )
	{
		for(j = 0; j < G_length; j++)
		{
			G[j] += J[i][j] * errors[i];
		}
	}

	/* symmetric matrix ??? Do we have singularities? */
	for(i = 0; i < G_length; i++)
	{
		for(j = 0; j < G_length; j++)
		{
			for(k = 0; k < num_samples * num_outputs; k++)
				H[i][j] += J[k][i]*J[k][j];
		}
	}
	NN_API_ERROR = NN_API_SUCCESS;;
	return sse;
error:
	NN_API_ERROR = NN_API_COMPUTE_JGH;
	return 0;
}
/* =====================================
 * = Neural network training functions =
 * =====================================
 */
int32_t nn_train(net_t * nn, double **samples, double ** targets, int32_t num_samples, int32_t max_epochs, double min_error )
{
	double ** J = NULL, ** A = NULL, * A_temp = NULL, * G = NULL, * errors = NULL, * nn_out = NULL, * nn_t_out = NULL, * work = NULL, * S = NULL;
	double niu = 0.01, current_error = 1e20, new_error = 0;
	double **r_samples = NULL, **r_targets = NULL;
	int32_t * iwork = NULL;
	int32_t l = 0, J_ydim = 0, nlvl = 0, epoch = 0, *perm = NULL, lwork = 0, s = 0, num_outputs = 0, num_inputs = 0;
	int32_t i = 0;
	net_layer_t *layer = NULL;
	net_t * test_nn = NULL;

	/* memory for rescaled samples */

	/* conditions check */
	C_CHECK_CONDITION( nn == NULL, 							 NN_API_NOT_INIT  );
	C_CHECK_CONDITION( nn->init_state == false, 			 NN_API_NOT_INIT  );
	C_CHECK_CONDITION( num_samples <= 0, 					 NN_API_BAD_INPUT );
	C_CHECK_CONDITION( max_epochs < 0 && max_epochs > 10000, NN_API_BAD_INPUT );
	C_CHECK_CONDITION( min_error < 0 && min_error > 1000,	 NN_API_BAD_INPUT );
	C_CHECK_CONDITION( samples == NULL, 					 NN_API_BAD_INPUT );
	C_CHECK_CONDITION( targets == NULL, 					 NN_API_BAD_INPUT );

	/* init weights */
	C_SAFE_CALL( init_nn_weights(nn) );

	/* first and last layer : required to know the number of inputs and outputs and so allocate memory */
	layer = nn->layer[nn->num_layers-1];
	num_outputs = layer->num_outputs;

	layer = nn->layer[0];
	num_inputs = layer->num_outputs;

	/* compute jacobian number of lines */
	for(l = 1; l < nn->num_layers; l++)
	{
		layer = nn->layer[l];
		J_ydim += (layer->num_inputs + 1) * layer-> num_outputs;
	}
	printf("J_ydim %d\n",J_ydim);

	/* initialize required memory */
	C_SAFE_CALL( J 		  = (double **)mem_alloc(num_outputs * num_samples	* sizeof(double  *), true) );
	C_SAFE_CALL( A 		  = (double **)mem_alloc(J_ydim     				* sizeof(double  *), true) );
	C_SAFE_CALL( G 	 	  = (double  *)mem_alloc(J_ydim		  				* sizeof(double   ), true) );
	C_SAFE_CALL( errors   = (double  *)mem_alloc(num_outputs * num_samples 	* sizeof(double   ), true) );
	C_SAFE_CALL( perm     = (int32_t *)mem_alloc(num_samples 				* sizeof(int32_t  ), true) );
	C_SAFE_CALL( nn_out   = (double  *)mem_alloc(num_outputs		   		* sizeof(double   ), true) );
	C_SAFE_CALL( nn_t_out = (double  *)mem_alloc(num_outputs		   		* sizeof(double   ), true) );

	C_SAFE_CALL( J[0]     	  = (double  *)mem_alloc(num_outputs * num_samples *  J_ydim * sizeof(double ), true) );
	C_SAFE_CALL( A[0]      	  = (double  *)mem_alloc(J_ydim * J_ydim  				 	 * sizeof(double ), true) );
	C_SAFE_CALL( A_temp   	  = (double  *)mem_alloc(J_ydim * J_ydim  					 * sizeof(double ), true) );
	C_SAFE_CALL( r_samples    = (double **)mem_alloc(num_samples 						 * sizeof(double*), true) );
	C_SAFE_CALL( r_samples[0] = (double  *)mem_alloc(num_samples * num_inputs 			 * sizeof(double*), true) );
	C_SAFE_CALL( r_targets    = (double **)mem_alloc(num_samples 						 * sizeof(double*), true) );
	C_SAFE_CALL( r_targets[0] = (double  *)mem_alloc(num_samples * num_outputs 			 * sizeof(double*), true) );

	/* calculate pointers address */
	for(l = 1; l < J_ydim; l++)
		A[l] = A[0] + l * J_ydim;

	for(l = 1; l < num_outputs * num_samples; l++)
		J[l] = J[0] + l * J_ydim;

	for(l=1; l < num_samples; l++){
		r_samples[l] = r_samples[0] + l * num_inputs;
		r_targets[l] = r_targets[0] + l * num_outputs;
	}

	/* query lapack for the best size of the work area */
	nlvl = (int32_t)(floor(log2((double)J_ydim/(25.0f+1.0f)))+1.0f);

	C_SAFE_CALL( S 		= (double  *) mem_alloc( J_ydim 							 * sizeof(double ), true) );
	C_SAFE_CALL( iwork 	= (int32_t *) mem_alloc( ( 3 * J_ydim * nlvl + 11 * J_ydim ) * sizeof(int32_t), true) );
	lwork = lpck_solve_ls_svd_pivot(J_ydim, -1, NULL, 0, S, iwork, A[0], G);
	C_SAFE_CALL( work = (double *)mem_alloc(lwork * sizeof(double), true) );

	/* grab scalling parameters */
	nn_grab_rescale_params( nn, samples, targets, num_samples );

	/* rescale samples and targets to -1, 1 range */
	for(l = 0; l < num_samples; l++){
		memcpy(r_samples[l], samples[l], num_inputs * sizeof(double));
		C_SAFE_CALL( transform_in_out_data(nn, r_samples[l], num_inputs, false) );

		memcpy(r_targets[l], targets[l], num_outputs * sizeof(double));
		for(i = 0; i < num_outputs; i++){
				r_targets[l][i] = map_min_max(r_targets[l][i], nn->ymin[i], nn->ymax[i], 0, 0, false);
		}

	}

	/* generate a duplicate network */
	C_SAFE_CALL( test_nn = duplicate_network( nn ) );

	/* disable rescaling */
	nn->use_mapping = false;

	/* generate random presentation */
	random_perm(num_samples, perm);
	C_SAFE_CALL( copy_nn_data( test_nn, nn ) );

	/* compute jacobian, hessian and stuff */
	current_error = compute_JGH(test_nn, J, G, A, r_samples, r_targets, errors, nn_out, nn_t_out, perm, num_samples, num_outputs, J_ydim);

	/* lets train the neural network */
	while(sqrt(current_error) > min_error && epoch < max_epochs && niu < 1e20)
	{

		/* copy network to a safe place */
		C_SAFE_CALL( copy_nn_data( test_nn, nn ) );

		memcpy(A_temp, A[0], J_ydim * J_ydim * sizeof(double));

		/* apply niu */
		for(i = 0; i < J_ydim; i++)
			 A_temp[i*J_ydim+i] += niu ;

		/* call lapack to solve our system of equations */
		C_SAFE_CALL( lpck_solve_ls_svd_pivot(J_ydim, 0, work, lwork, S, iwork, A_temp, G) );

		/* apply changes to neural network weights (lapack gives solution in G) */
		nn_change_weights_bias(test_nn, G);

		/* compute new error */
		new_error = 0.0f;
		for(s = 0; s < num_samples; s++)
		{
			/* short circuit the pointers */
			test_nn->in_layer = r_samples[perm[s]]; test_nn->out_layer = nn_out;

			/* compute neural network response for this input */
			C_SAFE_CALL( get_nn_response( test_nn ) );

			/* compute sample errors */
			for(i = 0; i < num_outputs; i++)
				new_error += (r_targets[perm[s]][i] - nn_out[i]) * (r_targets[perm[s]][i] - nn_out[i]);
		}

		/* check wether or not we should accept this network */
		if(new_error < current_error )
		{
			/* copy network */
			C_SAFE_CALL( copy_nn_data( nn, test_nn ) );
			current_error = new_error;
			niu *= 0.1;
			random_perm(num_samples, perm);

			/* compute new hessian, jacobian and gradient */
			current_error = compute_JGH(test_nn, J, G, A, r_samples, r_targets, errors, nn_out, nn_t_out, perm, num_samples, num_outputs, J_ydim);

		}else{
			niu *= 10;
		}

		printf("epoch: %d error : %f niu : %f\n",epoch,sqrt(current_error),niu);

		epoch++;
	}

	/* clean memory and exit */
	destroy_nn( test_nn );

	/* enable rescaling */
	nn->use_mapping = true;

	if(r_samples != NULL)
	{
		mem_free(r_samples[0]);
		mem_free(r_samples);
	}

	if(r_targets != NULL)
	{
		mem_free(r_targets[0]);
		mem_free(r_targets);
	}

	if(J != NULL)
	{
		mem_free( J[0] 	);
		mem_free( J 	);
	}
	if(A != NULL)
	{
		mem_free( A[0] 	);
		mem_free( A 	);
	}
	mem_free( A_temp	);
	mem_free( G 		);
	mem_free( errors 	);
	mem_free( perm 		);
	mem_free( nn_out	);
	mem_free( nn_t_out 	);
	mem_free( work 		);
	mem_free( iwork		);
	mem_free( S			);
	NN_API_ERROR = NN_API_SUCCESS;

	return 0;
error:
	destroy_nn( test_nn );

	/* clean memory and exit */
	if(r_samples != NULL)
	{
		mem_free(r_samples[0]);
		mem_free(r_samples);
	}

	if(r_targets != NULL)
	{
		mem_free(r_targets[0]);
		mem_free(r_targets);
	}

	if(J != NULL)
	{
		mem_free( J[0] 	);
		mem_free( J 	);
	}
	if(A != NULL)
	{
		mem_free( A[0] 	);
		mem_free( A 	);
	}
	mem_free( G 		);
	mem_free( errors 	);
	mem_free( perm 		);
	mem_free( nn_out	);
	mem_free( nn_t_out 	);
	mem_free( work 		);
	mem_free( iwork		);
	mem_free( S			);

	NN_API_ERROR = NN_API_TRAIN;
	return -1;
}

/* =====================================
 * = Neural network training functions =
 * =====================================
 */
int32_t nn_train2(net_t * nn, double **samples, double ** targets, int32_t num_samples, int32_t max_epochs, double min_error, double train_perc )
{
	double ** J = NULL, ** A = NULL, * A_temp = NULL, * G = NULL, * errors = NULL, * nn_out = NULL, * nn_t_out = NULL, * work = NULL, * S = NULL;
	double niu = 0.001, current_error = 1e20, new_error = 0, error_perc = 1- train_perc, grb = 0;
	double **r_samples = NULL, **r_targets = NULL;
	int32_t * iwork = NULL;
	int32_t l = 0, J_ydim = 0, nlvl = 0, epoch = 0, *perm = NULL, lwork = 0, s = 0, num_outputs = 0, num_inputs = 0;
	int32_t i = 0, train_s_index = 0, error_e_index = 0;
	net_layer_t *layer = NULL;
	net_t * test_nn = NULL;

	/* conditions check */
	C_CHECK_CONDITION( nn == NULL, 							 NN_API_NOT_INIT  );
	C_CHECK_CONDITION( nn->init_state == false, 			 NN_API_NOT_INIT  );
	C_CHECK_CONDITION( num_samples <= 0, 					 NN_API_BAD_INPUT );
	C_CHECK_CONDITION( max_epochs < 0 && max_epochs > 10000, NN_API_BAD_INPUT );
	C_CHECK_CONDITION( min_error < 0 && min_error > 1000,	 NN_API_BAD_INPUT );
	C_CHECK_CONDITION( samples == NULL, 					 NN_API_BAD_INPUT );
	C_CHECK_CONDITION( targets == NULL, 					 NN_API_BAD_INPUT );
	C_CHECK_CONDITION( train_perc <= 0 || train_perc >= 1, 	 NN_API_BAD_INPUT );
	C_CHECK_CONDITION( train_perc*num_samples < 10,		 	 NN_API_BAD_INPUT );
	C_CHECK_CONDITION( error_perc*num_samples < 10,		 	 NN_API_BAD_INPUT );

	/* init weights */
	C_SAFE_CALL( init_nn_weights(nn) );

	/* compute jacobian number of lines */
	for(l = 1; l < nn->num_layers; l++){
		layer = nn->layer[l];
		J_ydim += (layer->num_inputs + 1) * layer-> num_outputs;
	}

	/* first and last layer : required to know the number of inputs and outputs and so allocate memory */
	layer = nn->layer[nn->num_layers-1];
	num_outputs = layer->num_outputs;

	layer = nn->layer[0];
	num_inputs = layer->num_outputs;

	/* initialize required memory */
	C_SAFE_CALL( J 		  = (double **)mem_alloc(num_outputs * num_samples	* sizeof(double  *), true) );
	C_SAFE_CALL( A 		  = (double **)mem_alloc(J_ydim     				* sizeof(double  *), true) );
	C_SAFE_CALL( G 	 	  = (double  *)mem_alloc(J_ydim		  				* sizeof(double   ), true) );
	C_SAFE_CALL( errors   = (double  *)mem_alloc(num_outputs * num_samples 	* sizeof(double   ), true) );
	C_SAFE_CALL( perm     = (int32_t *)mem_alloc(num_samples 				* sizeof(int32_t  ), true) );
	C_SAFE_CALL( nn_out   = (double  *)mem_alloc(num_outputs		   		* sizeof(double   ), true) );
	C_SAFE_CALL( nn_t_out = (double  *)mem_alloc(num_outputs		   		* sizeof(double   ), true) );

	C_SAFE_CALL( J[0]     	  = (double  *)mem_alloc(num_outputs * num_samples *  J_ydim * sizeof(double ), true) );
	C_SAFE_CALL( A[0]      	  = (double  *)mem_alloc(J_ydim * J_ydim  				 	 * sizeof(double ), true) );
	C_SAFE_CALL( A_temp   	  = (double  *)mem_alloc(J_ydim * J_ydim  					 * sizeof(double ), true) );
	C_SAFE_CALL( r_samples    = (double **)mem_alloc(num_samples 						 * sizeof(double*), true) );
	C_SAFE_CALL( r_samples[0] = (double  *)mem_alloc(num_samples * num_inputs 			 * sizeof(double*), true) );
	C_SAFE_CALL( r_targets    = (double **)mem_alloc(num_samples 						 * sizeof(double*), true) );
	C_SAFE_CALL( r_targets[0] = (double  *)mem_alloc(num_samples * num_outputs 			 * sizeof(double*), true) );

	/* calculate pointers address */
	for(l = 1; l < J_ydim; l++)
		A[l] = A[0] + l * J_ydim;

	for(l = 1; l < num_outputs * num_samples; l++)
		J[l] = J[0] + l * J_ydim;

	for(l=1; l < num_samples; l++){
		r_samples[l] = r_samples[0] + l * num_inputs;
		r_targets[l] = r_targets[0] + l * num_outputs;
	}

	/* query lapack for the best size of the work area */
	nlvl = (int32_t)(floor(log2((double)J_ydim/(25.0f+1.0f)))+1.0f);

	C_SAFE_CALL( S 		= (double  *) mem_alloc( J_ydim 							 * sizeof(double ), true) );
	C_SAFE_CALL( iwork 	= (int32_t *) mem_alloc( ( 3 * J_ydim * nlvl + 11 * J_ydim ) * sizeof(int32_t), true) );
	lwork = lpck_solve_ls_svd_pivot(J_ydim, -1, NULL, 0, S, iwork, A[0], G);
	C_SAFE_CALL( work = (double *)mem_alloc(lwork * sizeof(double), true) );

	/* grab scalling parameters */
	nn_grab_rescale_params( nn, samples, targets, num_samples );

	/* rescale samples and targets to -1, 1 range */
	for(l = 0; l < num_samples; l++){
		memcpy(r_samples[l], samples[l], num_inputs * sizeof(double));
		C_SAFE_CALL( transform_in_out_data(nn, r_samples[l], num_inputs, false) );

		memcpy(r_targets[l], targets[l], num_outputs * sizeof(double));
		for(i = 0; i < num_outputs; i++){
			r_targets[l][i] = map_min_max(r_targets[l][i], nn->ymin[i], nn->ymax[i], 0, 0, false);
		}
	}

	/* generate a duplicate network */
	C_SAFE_CALL( test_nn = duplicate_network( nn ) );

	/* copy network to a safe place */
	C_SAFE_CALL( copy_nn_data( test_nn, nn ) );

	/* generate a random permutation */
	random_perm(num_samples, perm);

	/* training start index and error end index */
	train_s_index = num_samples * error_perc + 1;
	error_e_index = num_samples * error_perc;

	/* disable rescaling */
	nn->use_mapping = false;

	/* ASSUMPTION: FIRST SAMPLES ON PERM ARE FOR ERROR CALCULATION*/
	C_SAFE_CALL( current_error = compute_JGH(test_nn, J, G, A, r_samples, r_targets, errors, nn_out, nn_t_out, perm, error_e_index, num_outputs, J_ydim) );

	/* compute new hessian, jacobian and gradient */
	C_SAFE_CALL(	grb = compute_JGH(test_nn, J, G, A, r_samples, r_targets,
							&errors[train_s_index*num_outputs], nn_out, nn_t_out,
							&perm[train_s_index],
							num_samples - train_s_index,
							num_outputs, J_ydim);
				);

	/* lets train the neural network */
	while(sqrt(current_error) > min_error && epoch < max_epochs && niu < 1e20)
	{

		/* copy network to a safe place */
		C_SAFE_CALL( copy_nn_data( test_nn, nn ) );

		memcpy(A_temp, A[0], J_ydim * J_ydim * sizeof(double));

		/* apply niu */
		for(i = 0; i < J_ydim; i++)
			 A_temp[i*J_ydim+i] += niu ;

		/* call lapack to solve our system of equations */
		C_SAFE_CALL( lpck_solve_ls_svd_pivot(J_ydim, 0, work, lwork, S, iwork, A_temp, G) );

		/* apply changes to neural network weights (lapack gives solution in G) */
		nn_change_weights_bias(test_nn, G);

		/* compute new error */
		new_error = 0.0f;
		for(s = 0; s < error_e_index; s++)
		{
			/* short circuit the pointers */
			test_nn->in_layer = r_samples[perm[s]]; test_nn->out_layer = nn_out;

			/* compute neural network response for this input */
			C_SAFE_CALL( get_nn_response( test_nn ) );

			/* compute sample errors */
			for(i = 0; i < num_outputs; i++)
			{
				new_error += (r_targets[perm[s]][i] - nn_out[i]) * (r_targets[perm[s]][i] - nn_out[i]);
			}
		}

		/* check wether or not we should accept this network */
		if(new_error < current_error )
		{
			/* copy network */
			C_SAFE_CALL( copy_nn_data( nn, test_nn ) );
			current_error = new_error;
			/* compute new hessian, jacobian and gradient */
			C_SAFE_CALL(
				grb = compute_JGH(test_nn, J, G, A, r_samples, r_targets,
								&errors[train_s_index*num_outputs], nn_out, nn_t_out,
								&perm[train_s_index],
								num_samples - train_s_index,
								num_outputs, J_ydim);
			);
			niu *= 0.1;
		}else{

			niu *= 10;
		}

		printf("epoch: %d error : %f niu : %f\n",epoch,sqrt(current_error),niu);

		epoch++;
	}

	/* clean memory and exit */
	destroy_nn( test_nn );

	/* enable rescaling */
	nn->use_mapping = true;

	if(r_samples != NULL)
	{
		mem_free(r_samples[0]);
		mem_free(r_samples);
	}

	if(r_targets != NULL)
	{
		mem_free(r_targets[0]);
		mem_free(r_targets);
	}

	if(J != NULL)
	{
		mem_free( J[0] 	);
		mem_free( J 	);
	}
	if(A != NULL)
	{
		mem_free( A[0] 	);
		mem_free( A 	);
	}
	mem_free( A_temp 	);
	mem_free( G 		);
	mem_free( errors 	);
	mem_free( perm 		);
	mem_free( nn_out	);
	mem_free( nn_t_out 	);
	mem_free( work 		);
	mem_free( iwork		);
	mem_free( S			);
	NN_API_ERROR = NN_API_SUCCESS;

	return 0;
error:
	destroy_nn( test_nn );

	/* clean memory and exit */
	if(r_samples != NULL)
	{
		mem_free(r_samples[0]);
		mem_free(r_samples);
	}

	if(r_targets != NULL)
	{
		mem_free(r_targets[0]);
		mem_free(r_targets);
	}

	if(J != NULL)
	{
		mem_free( J[0] 	);
		mem_free( J 	);
	}
	if(A != NULL)
	{
		mem_free( A[0] 	);
		mem_free( A 	);
	}
	mem_free( A_temp 	);
	mem_free( G 		);
	mem_free( errors 	);
	mem_free( perm 		);
	mem_free( nn_out	);
	mem_free( nn_t_out 	);
	mem_free( work 		);
	mem_free( iwork		);
	mem_free( S			);

	NN_API_ERROR = NN_API_TRAIN2;
	return -1;
}

int32_t nn_train3(net_t * nn, double **samples, double ** targets, int32_t num_samples, int32_t max_epochs, double min_error )
{
	double ** J = NULL, ** A = NULL, *A_temp = NULL, * G = NULL, * errors = NULL, * nn_out = NULL, * nn_t_out = NULL, * work = NULL, * S = NULL, * inv_ipiv = NULL, * inv_work = NULL;
	double niu = 0.01, current_error = 1e20, new_error = 0,  alpha = 0, beta = 1, gamma = 0, Ew = 0, Ed = 0, trace = 0;
	double **r_samples = NULL, **r_targets = NULL;
	int32_t * iwork = NULL;
	int32_t l = 0, J_ydim = 0, nlvl = 0, epoch = 0, *perm = NULL, lwork = 0, s = 0, num_outputs = 0, num_inputs  =0;
	int32_t i = 0;
	net_layer_t *layer = NULL;
	net_t * test_nn = NULL;

	/* conditions check */
	C_CHECK_CONDITION( nn == NULL, 							 NN_API_NOT_INIT  );
	C_CHECK_CONDITION( nn->init_state == false, 			 NN_API_NOT_INIT  );
	C_CHECK_CONDITION( num_samples <= 0, 					 NN_API_BAD_INPUT );
	C_CHECK_CONDITION( max_epochs < 0 && max_epochs > 10000, NN_API_BAD_INPUT );
	C_CHECK_CONDITION( min_error < 0 && min_error > 1000,	 NN_API_BAD_INPUT );
	C_CHECK_CONDITION( samples == NULL, 					 NN_API_BAD_INPUT );
	C_CHECK_CONDITION( targets == NULL, 					 NN_API_BAD_INPUT );

	/* init weights */
	C_SAFE_CALL( init_nn_weights(nn) );

	/* compute jacobian number of lines */
	for(l = 1; l < nn->num_layers; l++){
		layer = nn->layer[l];
		J_ydim += (layer->num_inputs + 1) * layer-> num_outputs;
	}

	/* first and last layer : required to know the number of inputs and outputs and so allocate memory */
	layer = nn->layer[nn->num_layers-1];
	num_outputs = layer->num_outputs;

	layer = nn->layer[0];
	num_inputs = layer->num_outputs;

	/* initialize required memory */
	C_SAFE_CALL( J 		  = (double **)mem_alloc(num_outputs * num_samples	* sizeof(double  *), true) );
	C_SAFE_CALL( A 		  = (double **)mem_alloc(J_ydim     				* sizeof(double  *), true) );
	C_SAFE_CALL( G 	 	  = (double  *)mem_alloc(J_ydim		  				* sizeof(double   ), true) );
	C_SAFE_CALL( errors   = (double  *)mem_alloc(num_outputs * num_samples 	* sizeof(double   ), true) );
	C_SAFE_CALL( perm     = (int32_t *)mem_alloc(num_samples 				* sizeof(int32_t  ), true) );
	C_SAFE_CALL( nn_out   = (double  *)mem_alloc(num_outputs		   		* sizeof(double   ), true) );
	C_SAFE_CALL( nn_t_out = (double  *)mem_alloc(num_outputs		   		* sizeof(double   ), true) );

	C_SAFE_CALL( J[0]     	  = (double  *)mem_alloc(num_outputs * num_samples *  J_ydim * sizeof(double ), true) );
	C_SAFE_CALL( A[0]      	  = (double  *)mem_alloc(J_ydim * J_ydim  				 	 * sizeof(double ), true) );
	C_SAFE_CALL( A_temp   	  = (double  *)mem_alloc(J_ydim * J_ydim  					 * sizeof(double ), true) );
	C_SAFE_CALL( r_samples    = (double **)mem_alloc(num_samples 						 * sizeof(double*), true) );
	C_SAFE_CALL( r_samples[0] = (double  *)mem_alloc(num_samples * num_inputs 			 * sizeof(double*), true) );
	C_SAFE_CALL( r_targets    = (double **)mem_alloc(num_samples 						 * sizeof(double*), true) );
	C_SAFE_CALL( r_targets[0] = (double  *)mem_alloc(num_samples * num_outputs 			 * sizeof(double*), true) );

	/* calculate pointers address */
	for(l = 1; l < J_ydim; l++)
		A[l] = A[0] + l * J_ydim;

	for(l = 1; l < num_outputs * num_samples; l++)
		J[l] = J[0] + l * J_ydim;

	for(l=1; l < num_samples; l++){
		r_samples[l] = r_samples[0] + l * num_inputs;
		r_targets[l] = r_targets[0] + l * num_outputs;
	}

	/* query lapack for the best size of the work area */
	nlvl = (int32_t)(floor(log2((double)J_ydim/(25.0f+1.0f)))+1.0f);

	C_SAFE_CALL( S 		= (double  *) mem_alloc( J_ydim 							 * sizeof(double ), true) );
	C_SAFE_CALL( iwork 	= (int32_t *) mem_alloc( ( 3 * J_ydim * nlvl + 11 * J_ydim ) * sizeof(int32_t), true) );
	lwork = lpck_solve_ls_svd_pivot(J_ydim, -1, NULL, 0, S, iwork, A[0], G);
	C_SAFE_CALL( work = (double *)mem_alloc(lwork * sizeof(double), true) );

	/* allocate space for matrix inversion */
	C_SAFE_CALL( inv_ipiv = (double *) mem_alloc( (J_ydim + 1) * sizeof(double), true ) );
	C_SAFE_CALL( inv_work = (double *) mem_alloc( J_ydim * J_ydim * sizeof(double), true ) );

	/* grab scalling parameters */
	nn_grab_rescale_params( nn, samples, targets, num_samples );

	/* rescale samples and targets to -1, 1 range */
	for(l = 0; l < num_samples; l++){
		memcpy(r_samples[l], samples[l], num_inputs * sizeof(double));
		C_SAFE_CALL( transform_in_out_data(nn, r_samples[l], num_inputs, false) );

		memcpy(r_targets[l], targets[l], num_outputs * sizeof(double));
		for(i = 0; i < num_outputs; i++){
			r_targets[l][i] = map_min_max(r_targets[l][i], nn->ymin[i], nn->ymax[i], 0, 0, false);
		}
	}

	/* generate a duplicate network */
	C_SAFE_CALL( test_nn = duplicate_network( nn ) );

	/* copy network to a safe place */
	C_SAFE_CALL( copy_nn_data( test_nn, nn ) );

	/* generate a random permutation */
	random_perm(num_samples, perm);

	C_SAFE_CALL( Ed = compute_JGH(test_nn, J, G, A, r_samples, r_targets, errors, nn_out, nn_t_out, perm, num_samples, num_outputs, J_ydim) );

	/* compute Ew */
	Ew = 0.0;
	for( l = 1; l < test_nn -> num_layers; l++){
		layer = test_nn -> layer[l];
		for(i = 0; i < layer -> num_inputs * layer -> num_outputs; i++)
			Ew += layer -> weights[0][i] * layer -> weights[0][i];
		for(i = 0; i < layer->num_outputs; i++)
			Ew += layer->bias[i];
	}
	current_error = alpha * Ew + beta *Ed;

	/* disable rescaling */
	nn->use_mapping = false;

	/* lets train the neural network */
	while(sqrt(current_error) > min_error && epoch < max_epochs && niu < 1e10)
	{

		/* copy network to a safe place */
		C_SAFE_CALL( copy_nn_data( test_nn, nn ) );

		memcpy(A_temp, A[0], J_ydim * J_ydim * sizeof(double) );
		for(i=0; i < J_ydim * J_ydim; i++)
			A_temp[i] *= 2.0f * beta;

		/* apply niu and alpha regularization parameter */
		for(i = 0; i < J_ydim; i++)
			A_temp[i*J_ydim+i] += 2.0f*niu + 2.0f*alpha;

		/* call lapack to solve our system of equations */
		C_SAFE_CALL( lpck_solve_ls_svd_pivot(J_ydim, 0, work, lwork, S, iwork, A_temp, G) );

		/* apply changes to neural network weights (lapack gives solution in G) */
		nn_change_weights_bias(test_nn, G);

		/* compute new error */
		Ed = 0.0f;
		for(s = 0; s < num_samples; s++)
		{
			/* short circuit the pointers */
			test_nn->in_layer = r_samples[s]; test_nn->out_layer = nn_out;

			/* compute neural network response for this input */
			C_SAFE_CALL( get_nn_response( test_nn ) );

			/* compute sample errors */
			for(i = 0; i < num_outputs; i++)
				Ed += (r_targets[s][i] - nn_out[i])*(r_targets[s][i] - nn_out[i]);
		}

		/* compute Ew */
		Ew = 0.0;
		for( l = 1; l < test_nn -> num_layers; l++){
			layer = test_nn -> layer[l];
			for(i = 0; i < layer -> num_inputs * layer -> num_outputs; i++)
				Ew += layer -> weights[0][i] * layer -> weights[0][i];
			for(i = 0; i < layer->num_outputs; i++)
				Ew += layer->bias[i];
		}

		new_error = alpha * Ew + beta * Ed;

		/* check wether or not we should accept this network */
		if(new_error < current_error )
		{
			/* copy network */
			C_SAFE_CALL( copy_nn_data( nn, test_nn ) );
			current_error = new_error;
			niu *= 0.1;

			/* update parameters */
			memcpy(A_temp, A[0], J_ydim * J_ydim * sizeof(double) );
			//for(i=0; i < J_ydim * J_ydim; i++)
			//	A_temp[i] *= 2.0f * beta;

			/* apply niu and alpha regularization parameter */
			//for(i = 0; i < J_ydim; i++)
			//	A_temp[i*J_ydim+i] += 2.0f*alpha;

			C_SAFE_CALL( lpck_get_inv(J_ydim, inv_work, J_ydim * J_ydim, inv_ipiv, A_temp) );

			/* compute trace of the inverse */
			trace = 0.0;
			for(i = 0; i < J_ydim; i++)
				trace += A_temp[i*J_ydim+i];

			gamma = (double)(J_ydim) - (alpha * trace);
			beta =fabs(((double)(num_samples) - gamma)/(2.0f*Ed));
			alpha = (double)J_ydim /(2.0f*Ew+trace);

			/* compute new hessian, jacobian and gradient */
			C_SAFE_CALL( current_error = compute_JGH(test_nn, J, G, A, r_samples, r_targets,
								errors, nn_out, nn_t_out, perm, num_samples, num_outputs, J_ydim);
						);

		}else{
			niu *= 10;
		}

		printf("epoch: %d  alpha : %f beta : %f gamma : (%f ,%d, %d) trace : %f error : %f (%f, %f) niu : %f\n",epoch,alpha, beta, gamma, J_ydim,num_samples, trace,sqrt(current_error), Ed, Ew,niu);

		epoch++;
	}

	/* clean memory and exit */
	destroy_nn( test_nn );

	/* enable rescaling */
	nn->use_mapping = true;

	if(r_samples != NULL)
	{
		mem_free(r_samples[0]);
		mem_free(r_samples);
	}

	if(r_targets != NULL)
	{
		mem_free(r_targets[0]);
		mem_free(r_targets);
	}

	if(J != NULL)
	{
		mem_free( J[0] 	);
		mem_free( J 	);
	}
	if(A != NULL)
	{
		mem_free( A[0] 	);
		mem_free( A 	);
	}
	mem_free( A_temp    );
	mem_free( G 		);
	mem_free( errors 	);
	mem_free( perm 		);
	mem_free( nn_out	);
	mem_free( nn_t_out 	);
	mem_free( work 		);
	mem_free( iwork		);
	mem_free( S			);
	mem_free( inv_ipiv  );
	mem_free( inv_work  );

	NN_API_ERROR = NN_API_SUCCESS;

	return 0;
error:
	destroy_nn( test_nn );

	/* clean memory and exit */
	if(r_samples != NULL)
	{
		mem_free(r_samples[0]);
		mem_free(r_samples);
	}

	if(r_targets != NULL)
	{
		mem_free(r_targets[0]);
		mem_free(r_targets);
	}

	if(J != NULL)
	{
		mem_free( J[0] 	);
		mem_free( J 	);
	}
	if(A != NULL)
	{
		mem_free( A[0] 	);
		mem_free( A 	);
	}
	mem_free( G 		);
	mem_free( errors 	);
	mem_free( perm 		);
	mem_free( nn_out	);
	mem_free( nn_t_out 	);
	mem_free( work 		);
	mem_free( iwork		);
	mem_free( S			);

	NN_API_ERROR = NN_API_TRAIN3;
	return -1;
}
/* create a network with the same characteristics of the one that has been given */
net_t * duplicate_network(net_t * nn )
{
	int32_t * layers_func = NULL, * layers_conf = NULL;
	int32_t i = 0;
	net_t * new_nn = NULL;
	net_layer_t * layer;
	C_CHECK_CONDITION( nn == NULL, 					NN_API_BAD_INPUT );
	C_CHECK_CONDITION( nn -> layer == NULL, 		NN_API_BAD_INPUT );
	C_CHECK_CONDITION( nn -> init_state == false, 	NN_API_BAD_INPUT );
	C_CHECK_CONDITION( nn -> num_layers == 0, 		NN_API_BAD_INPUT );

	/* first allocate temporary memory */
	C_SAFE_CALL( layers_func = (int32_t *) mem_alloc( nn->num_layers * sizeof(int32_t), true ) );
	C_SAFE_CALL( layers_conf = (int32_t *) mem_alloc( nn->num_layers * sizeof(int32_t), true ) );

	/* get the parameters from this network */
	for(i = 0; i < nn->num_layers; i++)
	{
		layer = nn -> layer[i];
		C_CHECK_CONDITION( layer == NULL, NN_API_MALFORMED_NN );
		layers_conf[i] = layer->num_outputs;
		layers_func[i] = layer->eval_funct;

	}

	/* call network initializer */
	C_SAFE_CALL( new_nn = create_nn(nn->num_layers, layers_conf, layers_func ) );

	/* should we also copy the information? */

	/* input layer rescale parameters */
	for(i = 0; i < (nn->layer[0])->num_outputs; i++){
		new_nn -> xmin[i] = nn -> xmin[i];
		new_nn -> xmax[i] = nn -> xmax[i];
	}
	for(i = 0; i < (nn->layer[nn->num_layers-1])->num_outputs; i++){
		new_nn -> ymin[i] = nn -> ymin[i];
		new_nn -> ymax[i] = nn -> ymax[i];
	}

	/* free auxiliary memory areas */
	mem_free(layers_func);
	mem_free(layers_conf);

	NN_API_ERROR = NN_API_SUCCESS;
	return new_nn;
error:

	/* free auxiliary memory areas */
	mem_free(layers_func);
	mem_free(layers_conf);

	NN_API_ERROR = NN_API_DUPLICATE_NN;
	return NULL;
}

/* copy data from one network into another */
int32_t copy_nn_data(net_t *dest, net_t *src)
{
	net_layer_t * l_src = NULL, * l_dest = NULL;
	int32_t l = 0;

	C_CHECK_CONDITION( dest -> num_layers != src -> num_layers, NN_API_COPY_DIFF );
	C_CHECK_CONDITION( src -> num_layers == 0, 					NN_API_COPY_DIFF );

	for( l = 1; l < src -> num_layers; l++)
	{
		l_src = src->layer[l];
		l_dest = dest->layer[l];

		C_CHECK_CONDITION( l_src->num_inputs  != l_dest->num_inputs,  NN_API_COPY_DIFF );
		C_CHECK_CONDITION( l_src->num_outputs != l_dest->num_outputs, NN_API_COPY_DIFF );

		/* proceed with the copy */
		memcpy(l_dest->weights[0], l_src->weights[0], l_src->num_inputs * l_src->num_outputs * sizeof(double));
		memcpy(l_dest->bias, l_src->bias, l_src->num_outputs * sizeof(double));

	}
	dest->use_mapping = src->use_mapping;

	NN_API_ERROR = NN_API_SUCCESS;
	return 0;
error:
	NN_API_ERROR = NN_API_COPY;
	return -1;
}

void nn_grab_rescale_params(net_t * nn, double ** samples, double ** targets, int32_t num_samples )
{
	int32_t i = 0, j = 0;
	net_layer_t * layer = NULL;
	double min = 0, max = 0;
	layer = nn->layer[0];

	/* first find maximum and minimum of samples */
	min = samples[0][0]; max = min;

	for(i = 0; i < layer->num_outputs; i++)
	{
		for(j = 0; j < num_samples; j++ )
		{
			if(max < samples[j][i])
				max = samples[j][i];
			if(min > samples[j][i])
				min = samples[j][i];
		}
		/* store input parameters */
		nn->xmin[i] = min;
		nn->xmax[i] = max;
	}

	/* do it for targets*/
	layer = nn->layer[nn->num_layers-1];
	min = targets[0][0]; max = min;
	for(i = 0; i < layer->num_outputs; i++)
	{
		for(j = 0; j < num_samples; j++ )
		{
			if(max < targets[j][i])
				max = targets[j][i];
			if(min > targets[j][i])
				min = targets[j][i];
		}
		/* store output parameters */
		nn->ymin[i] = min;
		nn->ymax[i] = max;
	}
}

/*
 * Change weights and bias based on the values provided in G
 */
void nn_change_weights_bias(net_t * nn, double * G)
{
	net_layer_t * layer = NULL;
	int32_t l = 0, i = 0, j = 0, ct = 0;

	/* for each layer */
	for(l = 1; l < nn->num_layers; l++)
	{
		layer = nn->layer[l];

		for(i = 0; i < layer->num_outputs; i++)
		{
			/* weights */
			for(j = 0; j < layer->num_inputs; j++)
			{
				 layer->weights[i][j] -= G[ct++];
			}

			/* bias */
			layer->bias[i] -= G[ct++];
		}
	}

}

int32_t init_nn_weights(net_t * nn)
{
	int32_t i = 0, l = 0, n = 0;
	net_layer_t * layer = NULL;
	double norm_weights = 0.0f, norm_bias = 0.0f, beta = 0.0f;
	double f_range_up = 0.0f, f_range_down = 0.0f;

	C_CHECK_CONDITION( nn->num_layers == 0, NN_API_BAD_INPUT );

	for(l = 1; l < nn->num_layers; l++)
	{
		layer = nn->layer[l];

		/* scale random values according to the function used */
		switch(layer->eval_funct)
		{
			case      NN_LOGSIG: { f_range_up =  1.0f; f_range_down =  0.0f; break; }
			case      NN_TANSIG: { f_range_up =  1.0f; f_range_down = -1.0f; break; }
			default:  			 { f_range_up =  1.0f; f_range_down = -1.0f; break; }
		}

		beta = 0.7f*pow((double)layer->num_outputs, 1/(double)layer->num_inputs);

		/* for each neuron in the layer */
		norm_weights = 0.0f; norm_bias = 0.0f;
		for( n = 0; n < layer->num_outputs; n++)
		{
			for(i = 0; i < layer->num_inputs; i++ )
			{
				layer->weights[n][i] = (double)(mt_random()) * RANDMAXV * (f_range_up - f_range_down) - f_range_down;
				norm_weights += layer->weights[n][i] * layer->weights[n][i];
			}

			/* randomize bias */
			norm_weights = sqrt(norm_weights);
			norm_bias = sqrt(norm_bias);
			for(i = 0; i < layer->num_inputs; i++ )
				layer->weights[n][i] *=  (beta/norm_weights);

			// generate bias
			if(layer->weights[n][0]>0)
				layer->bias[n] = beta * (-1.0f + (double)n*2/(double)(layer->num_outputs));
			else
				layer->bias[n] = -1.0*beta * (-1.0f + (double)n*2/(double)(layer->num_outputs));
		}

	}

	NN_API_ERROR = NN_API_SUCCESS;
	return 0;

error:

	NN_API_ERROR = NN_API_WEIGHT_INIT;
	return -1;
}

int32_t lpck_solve_nxn_ls(int32_t J_ydim, int32_t query, int32_t * work, int32_t lwork, double *eq_matrix, double * eq_sol)
{
	int32_t col_size = J_ydim;
	int32_t car = 1;
	int32_t info = 0;
	if(query == -1)
	{
		return J_ydim;
	}else{
		dgesv_( &col_size, &car, eq_matrix, &col_size, work, eq_sol, &col_size, &info );

		/* Check for the exact singularity */
    	C_CHECK_CONDITION(info != 0, NN_API_LS_SOLVE);
    }
    NN_API_ERROR = NN_API_SUCCESS;
    return info;

error:
	if( info > 0 ) {
                printf( "\t\tThe diagonal element of the triangular factor of A,\n" );
                printf( "\t\tU(%i,%i) is zero, so that A is singular;\n", info, info );
                printf( "\t\tthe solution could not be computed.\n" );
                exit( 1 );
        }
    NN_API_ERROR = NN_API_LS;
	return -1;
}

int32_t lpck_solve_ls_svd_pivot(int32_t J_ydim, int32_t query, double * work, int32_t lwork, double * S, int32_t * iwork,  double *eq_matrix, double * eq_sol)
{
	int32_t col_size = J_ydim;
	int32_t nrhs = 1;
	int32_t info = 0;
	int32_t rank = 0;
	double work_size = 0;
	double rcond = -1;			/* default machine precision */

	if(query == -1)
	{
		dgelsd_( &col_size,		/* m dimension */
				&col_size,		/* n dimension */
				&nrhs, 			/* number of right hand sided */
				eq_matrix, 		/* A matrix */
				&col_size, 		/* leading dimension of A */
				eq_sol,			/* B vector */
				&col_size,		/* leading dimension of B */
				S, 				/* array with M columns */
				&rcond, 		/* which precision? */
				&rank, 			/* the matrix rank */
				&work_size,		/* query workspace size */
				&query,			/* the query */
				iwork, 			/* workspace array */
				&info );
		return (int32_t)work_size;
	}else{

		dgelsd_( &col_size,		/* m dimension */
				&col_size,		/* n dimension */
				&nrhs, 			/* number of right hand sided */
				eq_matrix, 		/* A matrix */
				&col_size, 		/* leading dimension of A */
				eq_sol,			/* B vector */
				&col_size,		/* leading dimension of B */
				S, 				/* array with M columns */
				&rcond, 		/* which precision? */
				&rank, 			/* the matrix rank */
				work,			/* query workspace size */
				&lwork,			/* the query */
				iwork, 			/* workspace array */
				&info );
		/* Check for the exact singularity */
    	C_CHECK_CONDITION(info != 0, NN_API_LS_SOLVE);
    }
    NN_API_ERROR = NN_API_SUCCESS;
    return info;

error:
	if( info > 0 ) {
        printf( "\t\tThe diagonal element of the triangular factor of A,\n" );
        printf( "\t\tU(%i,%i) is zero, so that A is singular;\n", info, info );
        printf( "\t\tthe solution could not be computed.\n" );
    }

    NN_API_ERROR = NN_API_LS;
	return -1;
}

int32_t lpck_solve_ls(int32_t J_ydim, int32_t query, double * work, int32_t lwork, double * eq_matrix,  double * eq_sol)
{
	char  N = 'N';
	int32_t col_size = J_ydim;
	int32_t car = 1;
	int32_t info = 0;
	double work_size = 7000;

	if (query==-1)
	{

		dgels_(  &N,			         /* use matrix M or its transpose */
			    &col_size,           /* number of lines of matrix M   */
				&col_size,           /* number of rows of matrix M    */
				&car,                /* number of colums of B and x   */
				eq_matrix,           /* matrix M                      */
				&col_size,           /* size -> max(1,vector_size)    */
				eq_sol,              /* sol vector / out vector       */
				&col_size,           /* LDB                           */
				&work_size,          /* work array                    */
				&query,              /* LWORK                         */
				&info);
		return (int32_t)work_size;
	}else{

		dgels_(  &N,				     /* use matrix M or its transpose */
			    &col_size,           /* number of lines of matrix M   */
				&col_size,           /* number of rows of matrix M    */
				&car,                /* number of colums of B and x   */
				eq_matrix,           /* matrix M                      */
				&col_size,           /* size -> max(1,vector_size)    */
				eq_sol,        		 /* sol vector / out vector       */
				&col_size,           /* LDB                           */
				work,				 /* work array                    */
				&lwork,              /* LWORK                         */
				&info);

		C_CHECK_CONDITION(info != 0, NN_API_LS_SOLVE);
	}

	NN_API_ERROR = NN_API_SUCCESS;
	return info;
error:
	if( info > 0 )
	{
        printf( "\t\tThe diagonal element %i of the triangular factor ", info );
        printf( "\t\tof A is zero, so that A does not have full rank;\n" );
        printf( "\t\tthe least squares solution could not be computed.\n" );
    }else{
    	printf("Lwork value %d\n", lwork);

    }
    NN_API_ERROR = NN_API_LS;
	return -1;
}

int32_t lpck_get_inv(int32_t J_ydim, double * work, int32_t lwork, double * ipiv, double * eq_matrix)
{
	int32_t col_size = J_ydim;
	int32_t info = 0;
	int32_t i = 0;

	dgetrf_(&col_size,&col_size, eq_matrix, &col_size, ipiv, &info);

	if(info!= 0){
		goto error;
	}
	dgetri_(&col_size, eq_matrix, &col_size, ipiv, work, &lwork, &info);

	if(info!= 0){
		goto error;
	}

	if(info!= 0)
		goto error;

	NN_API_ERROR = NN_API_SUCCESS;
	return info;
error:
	if(info > 0)
		printf("decomposition successful but U is singular, u_%d%d=0\n", i,i);
	else
		printf("decomposition not successful %dth parameter has an illegal value\n", i);

	NN_API_ERROR = NN_API_GET_INV;
	return -1;
}

void random_perm(int32_t num_p, int32_t * ranpat)
{
	unsigned int p, np, op;
	unsigned int r;
	for( p = 0 ; p < (unsigned int)num_p; p++ ) { ranpat[p] = (unsigned int)p; }

	for( p = 0 ; p < (unsigned int)num_p; p++)
	{
		r =  (unsigned int)( mt_random() * 1.0 * ( RANDMAXV  * (num_p - p )));
		np = p + r ;
		op = ranpat[p];
		ranpat[p] = ranpat[np];
		ranpat[np] = op;
	}
}

void init_mt(void)
{
    int32_t i = 0;
    srand((int32_t)time(NULL));
    /* init mersenne twister */
	mt_state[0]= rand() & 0xffffffffUL;
    for (i = 1; i < MT_N; i++)
	{
            mt_state[i] = (1812433253UL * (mt_state[i - 1] ^ (mt_state[i-1] >> 30)) + i);
        	/* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
       		/* In the previous versions, MSBs of the seed affect   */
        	/* only MSBs of the array mt[].                        */
        	/* 2002/01/09 modified by Makoto Matsumoto             */
        	mt_state[i] &= 0xffffffffUL;
        	/* for >32 bit machines */
    }

	mt_constants[0] = 0x0L;
	mt_constants[1] = 0x9908b0dfUL;
	mt_state_idx = MT_N+1;
}

unsigned long mt_random(void){
	unsigned long y;
	if(rand_state == false){ init_mt(); rand_state = true; }
	/* Updating state */
	if( mt_state_idx >= MT_N){
		int k;
		for(k = 0; k < MT_N - MT_M; k++)
		{
			y = ( mt_state[ k ] & MT_UPPER_MASK ) | ( mt_state[ k + 1 ] & MT_LOWER_MASK );
			mt_state[ k ] = mt_state[ k + MT_M ] ^ ( y>> 1) ^ mt_constants[ y & 0x1UL ];
		}
		for(; k < MT_N - 1; k++){
			y = ( mt_state[ k ] & MT_UPPER_MASK ) | ( mt_state[ k + 1 ] & MT_LOWER_MASK );
			mt_state[ k ] = mt_state[ k + ( MT_M - MT_N ) ] ^ ( y>> 1) ^ mt_constants[ y & 0x1UL ];
		}
		y = (mt_state[ MT_N - 1 ] & MT_UPPER_MASK ) | ( mt_state[0] & MT_LOWER_MASK );
		mt_state[ MT_N - 1 ] = mt_state[ MT_M - 1 ] ^ (y >> 1) ^ mt_constants[ y & 0x1UL ];
		mt_state_idx = 0;
	}
	/* Generating */
	y = mt_state[ mt_state_idx++ ];
	/* tempering */
	y ^= (y >> 11);
	y ^= (y << 7 )  & 0x9d2c5680UL;
	y ^= (y << 15 ) & 0xefc60000UL;
	y ^= (y >> 18);
	return y;
}

max_stars_repo_name: BrunoFFaria/neural-networks

max_issues_repo_name: BrunoFFaria/neural-networks

max_forks_repo_name: BrunoFFaria/neural-networks


text:
#ifndef QDM_CONFIG_H
#define QDM_CONFIG_H 1

#include <gsl/gsl_vector.h>

typedef struct {
  unsigned long int rng_seed;

  int acc_check;

  int burn_discovery;
  int iter_discovery;

  int burn_analysis;
  int iter_analysis;

  int knot_try;
  int spline_df;
  int thin;

  gsl_vector *months;
  gsl_vector *knots;

  gsl_vector *tau_highs;
  gsl_vector *tau_lows;

  double theta_min;
  double theta_tune_sd;

  double xi_high;
  double xi_low;
  double xi_prior_mean;
  double xi_prior_var;
  double xi_tune_sd;

  double bound;

  int debug;
  int tau_table;
} qdm_config;

int
qdm_config_read(
    qdm_config **config,
    const char *file_path,
    const char *group_path
);

void
qdm_config_fwrite(
    FILE *f,
    const qdm_config *cfg
);

#endif /* QDM_CONFIG_H */

ext: h

max_stars_repo_name: calebcase/qdm

max_issues_repo_name: calebcase/qdm

max_forks_repo_name: calebcase/qdm


text:
#ifndef SOURCE_PRINTTOFILE_H_
#define SOURCE_PRINTTOFILE_H_

// *** release/print_to_file.h ***
// Author: Kevin Wolz, date: 08/2018
//
// Utilities to write the results obtained by robustness_main into files.

#include <iostream>
#include <iomanip>
#include <fstream>
#include <string>
#include <sstream>
#include <vector>
#include <stdlib.h>
#include <sys/stat.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

#include "utils.h"
#include "md5hash.h"

// Appends the robustness (and, if analytical output, pdf) value to a file.
void WriteData(std::ofstream& datafile, double& robustness, double& pdf);

// Appends the initial conditions of a run (e.g. dimensions of data, covariance
// matrix etc.) to a file.
void WriteMetadata(std::ofstream& metafile, const unsigned& datdim,
                   gsl_matrix* covar_matrix, const unsigned& pardim,
                   gsl_vector* prior_param_values, gsl_matrix* prior_matrix,
                   gsl_matrix* g_matrix, const unsigned& num_samplings);

// Appends the metadata about a run (e.g. data hash ID, filepaths, timestamp,
// etc) to a file.
void WriteLogs(std::ofstream& logfile, std::string& meta_filepath,
               const int& op_mode, std::string& data_filepath,
               const int& success);

#endif // SOURCE_PRINTTOFILE_H_

max_stars_repo_name: kevguitar/Robustness

max_issues_repo_name: kevguitar/Robustness

max_forks_repo_name: kevguitar/Robustness


text:
/***************************************************************************/
/*                                                                         */
/* LU_decomp.h - LU Decomposition class for mruby                          */
/* Copyright (C) 2015 Paolo Bosetti                                        */
/* paolo[dot]bosetti[at]unitn.it                                           */
/* Department of Industrial Engineering, University of Trento              */
/*                                                                         */
/* This library is free software.  You can redistribute it and/or          */
/* modify it under the terms of the GNU GENERAL PUBLIC LICENSE 2.0.        */
/*                                                                         */
/* This library 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           */
/* Artistic License 2.0 for more details.                                  */
/*                                                                         */
/* See the file LICENSE                                                    */
/*                                                                         */
/***************************************************************************/

#ifndef LU_DECOMP_H
#define LU_DECOMP_H

#include <stdlib.h>
#include <string.h>
#include <errno.h>
#include <unistd.h>
#include <stdio.h>
#include <stdlib.h>
#include <sys/param.h>
#include <time.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_permutation.h>

#include "mruby.h"
#include "mruby/variable.h"
#include "mruby/string.h"
#include "mruby/data.h"
#include "mruby/class.h"
#include "mruby/value.h"
#include "mruby/array.h"
#include "mruby/numeric.h"
#include "mruby/compile.h"

#define E_LU_DECOMP_ERROR (mrb_class_get(mrb, "LUDecompError"))

/***********************************************\
 LU Decomposition
\***********************************************/

typedef struct {
  gsl_matrix *mat;
  gsl_permutation *p;
  int sgn;
  size_t size;
} lu_decomp_data_s;

// Garbage collector handler, for play_data struct
// if play_data contains other dynamic data, free it too!
// Check it with GC.start
void lu_decomp_destructor(mrb_state *mrb, void *p_);

// Utility function for getting the struct out of the wrapping IV @data
void mrb_lu_decomp_get_data(mrb_state *mrb, mrb_value self, lu_decomp_data_s **data);

void mrb_gsl_lu_decomp_init(mrb_state *mrb);

#endif // LU_DECOMP_H

max_stars_repo_name: UniTN-Mechatronics/mruby-gsl

max_issues_repo_name: UniTN-Mechatronics/mruby-gsl

max_forks_repo_name: UniTN-Mechatronics/mruby-gsl


text:
/**
 *
 * @file sgebrd.c
 *
 *  PLASMA computational routines
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Azzam Haidar
 * @date 2010-11-15
 * @generated s Tue Jan  7 11:45:10 2014
 *
 **/
#include <lapacke.h>
#include "common.h"
/***************************************************************************//**
 *
 * @ingroup float
 *
 *  PLASMA_sgebrd - reduces a general complex M-by-N matrix A to upper or lower
 *  bidiagonal form B using a two-stage approach
 *  First stage: reduction to band bidiagonal form (orthogonal matrices Q1 and P1);
 *  Second stage: reduction from band to bidiagonal form (orthogonal matrices
 *  Q2 and P2).
 *  Let Q = Q1 * Q2 be the global left unitary transformation;
 *  Let P = P1 * P2 be the global right unitary transformation;
 *  Q**T * A * P = B.
 *  If M >= N, B is upper bidiagonal; if M < N, B is lower bidiagonal.
 *  Not LAPACK Compliant for now!
 *  Note: T is incomplete and contains only the block reflectors of the first stage.
 *  Therefore, Q and P can not be built completely.
 *
 *******************************************************************************
 *
 * @param[in] jobq
 *          Specifies options for computing all or part of the matrix Q.
 *          Intended usage:
 *          = PlasmaVec: all M columns of Q are returned in array Q;
 *          = PlasmaNoVec: not referenced.
 *
 * @param[in] jobp
 *          Specifies options for computing all or part of the matrix V**T.
 *          Intended usage:
 *          = PlasmaVec: all N columns of P are returned in array P;
 *          = PlasmaNoVec: not referenced.
 *
 * @param[in] M
 *          The number of rows of the matrix A. M >= 0.
 *
 * @param[in] N
 *          The number of columns of the matrix A. N >= 0.
 *
 * @param[in,out] A
 *          On entry, the M-by-N matrix A.
 *          On exit,
 *          if M >= N, the diagonal and the first superdiagonal are
 *            overwritten with the upper bidiagonal matrix B; the
 *            elements below the diagonal, with the array T, represent
 *            the unitary matrix Q as a product of elementary
 *            reflectors, and the elements above the first superdiagonal,
 *            with the array T, represent the unitary matrix P as
 *            a product of elementary reflectors;
 *          if M < N, the diagonal and the first subdiagonal are
 *            overwritten with the lower bidiagonal matrix B; the
 *            elements below the first subdiagonal, with the array T,
 *            represent the unitary matrix Q as a product of
 *            elementary reflectors, and the elements above the diagonal,
 *            with the array T, represent the unitary matrix P as
 *            a product of elementary reflectors.
 *
 * @param[in] LDA
 *          The leading dimension of the array A. LDA >= max(1,M).
 *
 * @param[out] D
 *          On exit, the diagonal elements of the bidiagonal matrix:
 *          D(i) = A(i,i).
 *          Dimension (min(M,N)).
 *
 * @param[out] E
 *          On exit, the off-diagonal elements of the bidiagonal matrix:
 *          if M >= N, E(i) = A(i,i+1) for i = 1,2,...,N-1;
 *          if M < N, E(i) = A(i+1,i) for i = 1,2,...,M-1.
 *          Dimension (min(M,N)-1).
 *
 * @param[out] descT
 *          On entry, descriptor as return by PLASMA_Alloc_Workspace_sgebrd
 *          On exit, contains auxiliary factorization data.
 *
 * @param[out] Q
 *          On exit, if jobz = PlasmaVec, then if return value = 0,
 *          Q contains the M-by-M unitary matrix Q.
 *          If jobz = PlasmaNoVec, then it is not referenced.
 *
 * @param[in] LDQ
 *          The leading dimension of the array Q. LDQ >= M.
 *
 * @param[out] P
 *          On exit, if jobz = PlasmaVec, then if return value = 0,
 *          P contains the N-by-N unitary matrix P.
 *          If jobz = PlasmaNoVec, then it is not referenced.
 *
 * @param[in] LDP
 *          The leading dimension of the array P. LDP >= N.
 *
 *******************************************************************************
 *
 * @return
 *          \retval PLASMA_SUCCESS successful exit
 *          \retval <0 if -i, the i-th argument had an illegal value
 *
 *******************************************************************************
 *
 * @sa PLASMA_sgebrd_Tile
 * @sa PLASMA_sgebrd_Tile_Async
 * @sa PLASMA_cgebrd
 * @sa PLASMA_dgebrd
 * @sa PLASMA_sgebrd
 *
 ******************************************************************************/
int PLASMA_sgebrd(PLASMA_enum jobq, PLASMA_enum jobp,
                  int M, int N,
                  float *A, int LDA,
                  float *D,
                  float *E,
                  PLASMA_desc *descT,
                  float *Q, int LDQ,
                  float *P, int LDP)
{
    int NB;
    int status;
    plasma_context_t *plasma;
    PLASMA_sequence *sequence = NULL;
    PLASMA_request request = PLASMA_REQUEST_INITIALIZER;
    PLASMA_desc descA;

    plasma = plasma_context_self();
    if (plasma == NULL) {
        plasma_fatal_error("PLASMA_sgebrd", "PLASMA not initialized");
        return PLASMA_ERR_NOT_INITIALIZED;
    }

    /* Check input arguments */
    if (jobq != PlasmaNoVec  && jobq !=PlasmaVec) {
        plasma_error("PLASMA_sgebrd", "illegal value of jobq");
        return -1;
    }
    if (jobp != PlasmaNoVec && jobp != PlasmaVec) {
        plasma_error("PLASMA_sgebrd", "illegal value of jobp");
        return -2;
    }
    if (M < 0) {
        plasma_error("PLASMA_sgebrd", "illegal value of M");
        return -3;
    }
    if (N < 0) {
        plasma_error("PLASMA_sgebrd", "illegal value of N");
        return -4;
    }
    if (LDA < max(1, M)) {
        plasma_error("PLASMA_sgebrd", "illegal value of LDA");
        return -6;
    }
    if (LDQ < 1) {
        plasma_error("PLASMA_sgebrd", "illegal value of LDQ");
        return -9;
    }
    if (LDP < 1) {
        plasma_error("PLASMA_sgebrd", "illegal value of LDP");
        return -11;
    }
    /* Quick return */
    if (min(M, N) == 0) {
        return PLASMA_SUCCESS;
    }

    /* Tune NB & IB depending on M & N; Set NBNB */
    status = plasma_tune(PLASMA_FUNC_SGEBRD, M, N, 0);
    if (status != PLASMA_SUCCESS) {
        plasma_error("PLASMA_sgebrd", "plasma_tune() failed");
        return status;
    }

    /* Set NT */
    NB = PLASMA_NB;

    plasma_sequence_create(plasma, &sequence);

    if ( PLASMA_TRANSLATION == PLASMA_OUTOFPLACE ) {
        plasma_sooplap2tile( descA, A, NB, NB, LDA, N, 0, 0, M, N, sequence, &request,
                             plasma_desc_mat_free(&(descA)) );
    } else {
        plasma_siplap2tile( descA, A, NB, NB, LDA, N, 0, 0, M, N,
                            sequence, &request);
    }

    /* Call the tile interface */
    PLASMA_sgebrd_Tile_Async(jobq, jobp, &descA, D, E, descT, Q, LDQ, P, LDP, sequence, &request);

    if ( PLASMA_TRANSLATION == PLASMA_OUTOFPLACE ) {
        plasma_sooptile2lap( descA, A, NB, NB, LDA, N,  sequence, &request);
        plasma_dynamic_sync();
        plasma_desc_mat_free(&descA);
    } else {
        plasma_siptile2lap( descA, A, NB, NB, LDA, N,  sequence, &request);
        plasma_dynamic_sync();
    }

    status = sequence->status;
    plasma_sequence_destroy(plasma, sequence);
    return status;
}

/***************************************************************************//**
 *
 * @ingroup float_Tile
 *
 *  PLASMA_sgebrd_Tile - reduces a general complex M-by-N matrix A to upper or lower
 *  bidiagonal form B using a two-stage approach
 *  First stage: reduction to band bidiagonal form (orthogonal matrices Q1 and P1);
 *  Second stage: reduction from band to bidiagonal form (orthogonal matrices
 *  Q2 and P2).
 *  Let Q = Q1 * Q2 be the global left unitary transformation;
 *  Let P = P1 * P2 be the global right unitary transformation;
 *  Q**T * A * P = B.
 *  If M >= N, B is upper bidiagonal; if M < N, B is lower bidiagonal.
 *  Note: T is incomplete and contains only the block reflectors of the first stage.
 *  Therefore, Q and P can not be built completely.
 *  Tile equivalent of PLASMA_sgebrd().
 *  Operates on matrices stored by tiles.
 *  All matrices are passed through descriptors.
 *  All dimensions are taken from the descriptors.
 *
 *******************************************************************************
 *
 * @param[in] jobq
 *          Specifies options for computing all or part of the matrix Q.
 *          Intended usage:
 *          = PlasmaVec: all M columns of Q are returned in array Q;
 *          = PlasmaNoVec: not referenced.
 *
 * @param[in] jobp
 *          Specifies options for computing all or part of the matrix V**T.
 *          Intended usage:
 *          = PlasmaVec: all M columns of Q are returned in array Q;
 *          = PlasmaNoVec: not referenced.
 *
 * @param[in,out] A
 *          On entry, the M-by-N matrix A.
 *          On exit,
 *          if M >= N, the diagonal and the first superdiagonal are
 *            overwritten with the upper bidiagonal matrix B; the
 *            elements below the diagonal, with the array T, represent
 *            the unitary matrix Q as a product of elementary
 *            reflectors, and the elements above the first superdiagonal,
 *            with the array T, represent the unitary matrix P as
 *            a product of elementary reflectors;
 *          if M < N, the diagonal and the first subdiagonal are
 *            overwritten with the lower bidiagonal matrix B; the
 *            elements below the first subdiagonal, with the array T,
 *            represent the unitary matrix Q as a product of
 *            elementary reflectors, and the elements above the diagonal,
 *            with the array T, represent the unitary matrix P as
 *            a product of elementary reflectors.
 *
 * @param[out] D
 *          The real array containing the diagonal elements
 *          of the bidiagonal matrix B:
 *          D(i) = A(i,i).
 *          Dimension (min(M,N)).
 *
 * @param[out] E
 *          The real array containing the off-diagonal elements
 *          of the bidiagonal matrix B:
 *          if M >= N, E(i) = A(i,i+1) for i = 1,2,...,N-1;
 *          if M < N, E(i) = A(i+1,i) for i = 1,2,...,M-1.
 *          Dimension (min(M,N)-1).
 *
 * @param[out] T
 *          On exit, contains auxiliary factorization data.
 *
 * @param[out] Q
 *          On exit, if jobz = PlasmaVec, then if return value = 0,
 *          Q contains the M-by-M unitary matrix Q.
 *          If jobz = PlasmaNoVec, then it is not referenced.
 *
 * @param[in] LDQ
 *          The leading dimension of the array Q. LDQ >= M.
 *
 * @param[out] P
 *          On exit, if jobz = PlasmaVec, then if return value = 0,
 *          P contains the N-by-N unitary matrix P.
 *          If jobz = PlasmaNoVec, then it is not referenced.
 *
 * @param[in] LDP
 *          The leading dimension of the array P. LDP >= N.
 *
 *******************************************************************************
 *
 * @return
 *          \return PLASMA_SUCCESS successful exit
 *
 *******************************************************************************
 *
 * @sa PLASMA_sgebrd
 * @sa PLASMA_sgebrd_Tile_Async
 * @sa PLASMA_cgebrd_Tile
 * @sa PLASMA_dgebrd_Tile
 * @sa PLASMA_sgebrd_Tile
 *
 ******************************************************************************/
int PLASMA_sgebrd_Tile(PLASMA_enum jobq, PLASMA_enum jobp,
                       PLASMA_desc *A,
                       float *D, float *E,
                       PLASMA_desc *T,
                       float *Q, int LDQ,
                       float *P, int LDP)
{
    plasma_context_t *plasma;
    PLASMA_sequence *sequence = NULL;
    PLASMA_request request = PLASMA_REQUEST_INITIALIZER;
    int status;

    plasma = plasma_context_self();
    if (plasma == NULL) {
        plasma_fatal_error("PLASMA_sgebrd_Tile", "PLASMA not initialized");
        return PLASMA_ERR_NOT_INITIALIZED;
    }
    plasma_sequence_create(plasma, &sequence);
    PLASMA_sgebrd_Tile_Async(jobq, jobp, A, D, E, T, Q, LDQ, P, LDP, sequence, &request);
    plasma_dynamic_sync();
    status = sequence->status;
    plasma_sequence_destroy(plasma, sequence);
    return status;
}

/***************************************************************************//**
 *
 * @ingroup float_Tile_Async
 *
 *  PLASMA_sgebrd_Tile_Async - reduces a general complex M-by-N matrix A to upper or lower
 *  bidiagonal form B using a two-stage approach
 *  First stage: reduction to band bidiagonal form (orthogonal matrices Q1 and P1);
 *  Second stage: reduction from band to bidiagonal form (orthogonal matrices
 *  Q2 and P2).
 *  Let Q = Q1 * Q2 be the global left unitary transformation;
 *  Let P = P1 * P2 be the global right unitary transformation;
 *  Q**T * A * P = B.
 *  If M >= N, B is upper bidiagonal; if M < N, B is lower bidiagonal.
 *  Note: T is incomplete and contains only the block reflectors of the first stage.
 *  Therefore, Q and P can not be built completely.
 *  Non-blocking equivalent of PLASMA_sgebrd_Tile().
 *  May return before the computation is finished.
 *  Allows for pipelining of operations at runtime.
 *
 *******************************************************************************
 *
 * @param[in] sequence
 *          Identifies the sequence of function calls that this call belongs to
 *          (for completion checks and exception handling purposes).
 *
 * @param[out] request
 *          Identifies this function call (for exception handling purposes).
 *
 *******************************************************************************
 *
 * @sa PLASMA_sgebrd
 * @sa PLASMA_sgebrd_Tile
 * @sa PLASMA_cgebrd_Tile_Async
 * @sa PLASMA_dgebrd_Tile_Async
 * @sa PLASMA_sgebrd_Tile_Async
 *
 ******************************************************************************/
int PLASMA_sgebrd_Tile_Async(PLASMA_enum jobu, PLASMA_enum jobvt,
                             PLASMA_desc *A,
                             float *S, float *E,
                             PLASMA_desc *T,
                             float *U, int LDU,
                             float *VT, int LDVT,
                             PLASMA_sequence *sequence, PLASMA_request *request)
{
    PLASMA_desc descA ;
    PLASMA_desc descT ;
    PLASMA_desc descU, descVT;
    float *AB;
    int M     ;
    int N     ;
    int MINMN ;
    int NB    ;
    int LDAB  ;
    int i;
    int status;

    plasma_context_t *plasma;
    plasma = plasma_context_self();

    if (plasma == NULL) {
        plasma_fatal_error("PLASMA_sgebrd_Tile_Async", "PLASMA not initialized");
        return PLASMA_ERR_NOT_INITIALIZED;
    }
    if (sequence == NULL) {
        plasma_fatal_error("PLASMA_sgebrd_Tile_Async", "NULL sequence");
        return PLASMA_ERR_UNALLOCATED;
    }
    if (request == NULL) {
        plasma_fatal_error("PLASMA_sgebrd_Tile_Async", "NULL request");
        return PLASMA_ERR_UNALLOCATED;
    }
    /* Check sequence status */
    if (sequence->status == PLASMA_SUCCESS)
        request->status = PLASMA_SUCCESS;
    else
        return plasma_request_fail(sequence, request, PLASMA_ERR_SEQUENCE_FLUSHED);

    /* Check descriptors for correctness */
    if (plasma_desc_check(A) != PLASMA_SUCCESS) {
        plasma_error("PLASMA_sgebrd_Tile_Async", "invalid first descriptor");
        return plasma_request_fail(sequence, request, PLASMA_ERR_ILLEGAL_VALUE);
    } else {
        descA = *A;
    }
    if (plasma_desc_check(T) != PLASMA_SUCCESS) {
        plasma_error("PLASMA_sgebrd_Tile_Async", "invalid fourth descriptor");
        return plasma_request_fail(sequence, request, PLASMA_ERR_ILLEGAL_VALUE);
    } else {
        descT = *T;
    }
    /* Check input arguments */
    if (jobu != PlasmaNoVec  && jobu != PlasmaVec) {
        plasma_error("PLASMA_sgebrd_Tile_Async", "illegal value of jobu");
        return PLASMA_ERR_NOT_SUPPORTED;
    }
    if (jobvt != PlasmaNoVec && jobvt != PlasmaVec) {
        plasma_error("PLASMA_sgebrd_Tile_Async", "illegal value of jobvt");
        return PLASMA_ERR_NOT_SUPPORTED;
    }
    if (descA.nb != descA.mb) {
        plasma_error("PLASMA_sgebrd_Tile_Async", "only square tiles supported");
        return plasma_request_fail(sequence, request, PLASMA_ERR_ILLEGAL_VALUE);
    }

    #if defined(ENABLE_TIMER)
    PLASMA_Double_t timelpk=0.0,timeaplQ2=0.0, timeT=0.0;
    PLASMA_Double_t timeB=0.0,timeblg=0.0,timeaplQ1=0.0,timeconv1=0.0,timeconv2=0.0,timeall=0.0;
    timeall = PLASMA_Wtime();
    #endif
    PLASMA_enum uplo = descA.m >= descA.n ? PlasmaUpper : PlasmaLower;
    M     = descA.m;
    N     = descA.n;
    MINMN = min(M,N);
    NB    = min(descA.mb,MINMN);
    LDAB  = 3*NB+1;
    /*=======================================
     *  case M<NB or N<NB call lapack
     *=======================================*/
    if( ( M<= NB) || ( N <=NB ) ){
        /* convert the tile descA to lapack A */
        if ( PLASMA_TRANSLATION == PLASMA_OUTOFPLACE ) {
            AB = (float *) plasma_shared_alloc(plasma, M*N, PlasmaRealFloat);
            if (AB == NULL) {
                plasma_error("PLASMA_sgesvd_Tile_Async", "plasma_shared_alloc(AB-0-) failed");
                plasma_shared_free(plasma, AB);
                return PLASMA_ERR_OUT_OF_RESOURCES;
            }
            plasma_sooptile2lap( descA, AB, NB, NB, M, N,  sequence, request);
        } else {
            AB = descA.mat;
            plasma_siptile2lap( descA, AB, NB, NB, M, N,  sequence, request);
        }
        plasma_dynamic_sync();
        /*=======================================
         *  calling LAPACK SGESVD
         *=======================================*/
        #if defined(ENABLE_TIMER)
        plasma_dynamic_sync();
        timelpk   = PLASMA_Wtime();
        #endif
        plasma_setlapack_multithreads(plasma->world_size);
        /* call SVD solver using lapack routine  */
        float *TAUQ = (float *) plasma_shared_alloc(plasma, MINMN, PlasmaRealFloat);
        float *TAUP = (float *) plasma_shared_alloc(plasma, MINMN, PlasmaRealFloat);
        status = LAPACKE_sgebrd(LAPACK_COL_MAJOR, M, N, AB, M, S, E, TAUQ, TAUP);
        if(jobu == PlasmaVec){
            LAPACKE_slacpy(LAPACK_COL_MAJOR, lapack_const(PlasmaLower), M, MINMN, AB, M, U, LDU );
            LAPACKE_sorgbr(LAPACK_COL_MAJOR, 'Q', M, M, N, U, LDU, TAUQ );
        }
        if(jobvt == PlasmaVec){
            LAPACKE_slacpy(LAPACK_COL_MAJOR, lapack_const(PlasmaUpper), MINMN, N, AB, M, VT, LDVT );
            LAPACKE_sorgbr(LAPACK_COL_MAJOR, 'P', N, N, MINMN, VT, LDVT, TAUP );
        }
        if(status != 0){
            plasma_error("PLASMA_sgesvd","SGESVD");
        }
        sequence->status = status;
        plasma_setlapack_sequential(plasma);
        #if defined(ENABLE_TIMER)
        timelpk = PLASMA_Wtime()-timelpk;
        printf("  Finish Eigensolver-lpkonly timing= %lf  threads %d\n" ,timelpk, plasma->world_size);
        #endif
        /*=======================================
         *  END of calling SVD solver
         *=======================================*/
        /* convert the lapack to the tile descA */
        if ( PLASMA_TRANSLATION == PLASMA_OUTOFPLACE ) {
            //plasma_sooplap2tile_noalloc( descA, AB, NB, NB,  M, N, 0, 0, M, N, sequence, request);
            plasma_parallel_call_5( plasma_pslapack_to_tile,
            float*, AB,
            int,                 M,
            PLASMA_desc,         descA,
            PLASMA_sequence*,    sequence,
            PLASMA_request*,     request);
            plasma_dynamic_sync();
            free(AB);
        } else {
            plasma_siplap2tile( descA, AB, NB, NB,  M, N, 0, 0, M, N,
                                sequence, request);
        }
        plasma_dynamic_sync();
        free(TAUQ);
        free(TAUP);
        return PLASMA_SUCCESS;
    }
    /*=======================================
     *  END OF case M<NB or N<NB
     *=======================================*/
    /*
     * Allocate workspace for band storage of the band matrix A
     * AB looks like:
     *       __________________________________
     * NB   |               zero               |
     *       ----------------------------------
     * NB+1 |               band A             |
     *       ----------------------------------
     * NB   |_______________zero_______________|
     *
     * */
    AB = (float *)plasma_shared_alloc(plasma, LDAB*MINMN, PlasmaRealFloat);
    memset( AB, 0, LDAB * MINMN * sizeof(float) );
    if (AB == NULL) {
        plasma_error("PLASMA_sgebrd_Tile_Async", "plasma_shared_alloc(AB) failed");
        plasma_shared_free(plasma, AB);
        return PLASMA_ERR_OUT_OF_RESOURCES;
    }
    /*=======================================
     *  calling Reduction from DENSE to BAND
     *  then convert matrix to band form
     *=======================================*/
    //plasma_dynamic_sync();
    #if defined(ENABLE_TIMER)
    timeB   = PLASMA_Wtime();
    #endif
    /*
     * Reduction to BAND bidiagonal form
     * May be further optimized using the algo described in Trefethen
     */
    /* if (plasma->householder == PLASMA_FLAT_HOUSEHOLDER) { */
        plasma_dynamic_call_4(plasma_psgebrd_ge2gb,
            PLASMA_desc, descA,
            PLASMA_desc, descT,
            PLASMA_sequence*, sequence,
            PLASMA_request*, request);
    /* } */
    /* else { */
    /*     plasma_dynamic_call_4(plasma_psgebrd_ge2gb_rh, */
    /*         PLASMA_desc, descA, */
    /*         PLASMA_desc, descT, */
    /*         PLASMA_sequence*, sequence, */
    /*         PLASMA_request*, request); */
    /* } */
    //plasma_dynamic_sync();
    plasma_dynamic_call_6( plasma_psgbcpy_t2bl,
        PLASMA_enum, descA.m >= descA.n ? PlasmaUpper : PlasmaLower,
        PLASMA_desc, descA,
        float*, &(AB[NB]),
        int, LDAB,
        PLASMA_sequence*, sequence,
        PLASMA_request*, request);
    plasma_dynamic_sync();
    status = sequence->status;
    if (status != PLASMA_SUCCESS) {
        plasma_error("PLASMA_sgebrd","psgebrd_ge2gb+pzcopy");
        return status;
    }
    #if defined(ENABLE_TIMER)
    timeB   = PLASMA_Wtime()-timeB;
    printf("\n  Finish Band red    timing= %lf \n",timeB);
    #endif
    /*=======================================
     *  END of calling Reduction to BAND
     *=======================================*/
    /*=======================================
     *  calling Reduction from BAND to bidiag
     *=======================================*/
    float *VQ2   = NULL;
    float *VP2   = NULL;
    float *TAUQ2 = NULL;
    float *TAUP2 = NULL;
    float *TQ2   = NULL;
    float *TP2   = NULL;
    int Vblksiz, blkcnt, LDT, LDV;
    int WANTZ   = 0;

    if( jobu == PlasmaNoVec )
        WANTZ=0;
    else
        WANTZ=1;

    /* Vblksiz correspond to the blocking used when applying V2 to the matrix U
     * it is similar to IB in LAPACK SORMQR.
     * blkcnt is the number of diamond or tile of Vs */
    /* Note that in case PlamaVec requested, the V2 and T2 are stored by the
     * bulgechasing function in a special format:
     * for V2s: it store the V2(LDV,Vblksiz) of each diamond in a tile storage meaning
     * that V2_1 is stored then V2_2,..., V2_blkcnt.
     * blkcnt is the number of diamond.
     * */
    Vblksiz = min(NB,48);
    LDT     = Vblksiz;
    /* data for U */
    if( jobu == PlasmaVec ) {
        findVTsiz(MINMN, NB, Vblksiz, &blkcnt, &LDV);
        TAUQ2   = (float *) plasma_shared_alloc(plasma,     blkcnt*Vblksiz, PlasmaRealFloat);
        VQ2     = (float *) plasma_shared_alloc(plasma, LDV*blkcnt*Vblksiz, PlasmaRealFloat);
        TQ2     = (float *) plasma_shared_alloc(plasma, LDT*blkcnt*Vblksiz, PlasmaRealFloat);
        if ( (TAUQ2 == NULL) || (VQ2 == NULL) || (TQ2 == NULL) ) {
            plasma_error("PLASMA_sgebrd", "plasma_shared_alloc() failed");
            plasma_shared_free(plasma, TAUQ2);
            plasma_shared_free(plasma, VQ2);
            plasma_shared_free(plasma, TQ2);
            return PLASMA_ERR_OUT_OF_RESOURCES;
        }
        memset(TAUQ2, 0,     blkcnt*Vblksiz*sizeof(float));
        memset(VQ2,   0, LDV*blkcnt*Vblksiz*sizeof(float));
        memset(TQ2,   0, LDT*blkcnt*Vblksiz*sizeof(float));
    }
    else {
        TAUQ2   = (float *) plasma_shared_alloc(plasma, 2*MINMN, PlasmaRealFloat);
        VQ2     = (float *) plasma_shared_alloc(plasma, 2*MINMN, PlasmaRealFloat);
        if ( (TAUQ2 == NULL) || (VQ2 == NULL) ) {
            plasma_error("PLASMA_sgebrd", "plasma_shared_alloc() failed");
            plasma_shared_free(plasma, TAUQ2);
            plasma_shared_free(plasma, VQ2);
            return PLASMA_ERR_OUT_OF_RESOURCES;
        }
        memset(TAUQ2, 0, 2*MINMN*sizeof(float));
        memset(VQ2,   0, 2*MINMN*sizeof(float));
    }
    /* data for VT */
    if( jobvt == PlasmaVec ) {
        findVTsiz(MINMN, NB, Vblksiz, &blkcnt, &LDV);
        TAUP2   = (float *) plasma_shared_alloc(plasma,     blkcnt*Vblksiz, PlasmaRealFloat);
        VP2     = (float *) plasma_shared_alloc(plasma, LDV*blkcnt*Vblksiz, PlasmaRealFloat);
        TP2     = (float *) plasma_shared_alloc(plasma, LDT*blkcnt*Vblksiz, PlasmaRealFloat);
        if ( (TAUP2 == NULL) || (VP2 == NULL) || (TP2 == NULL) ) {
            plasma_error("PLASMA_sgebrd", "plasma_shared_alloc() failed");
            plasma_shared_free(plasma, TAUP2);
            plasma_shared_free(plasma, VP2);
            plasma_shared_free(plasma, TP2);
            return PLASMA_ERR_OUT_OF_RESOURCES;
        }
        memset(TAUP2, 0,     blkcnt*Vblksiz*sizeof(float));
        memset(VP2,   0, LDV*blkcnt*Vblksiz*sizeof(float));
        memset(TP2,   0, LDT*blkcnt*Vblksiz*sizeof(float));
    }
    else {
        TAUP2   = (float *) plasma_shared_alloc(plasma, 2*MINMN, PlasmaRealFloat);
        VP2     = (float *) plasma_shared_alloc(plasma, 2*MINMN, PlasmaRealFloat);
        if ( (TAUP2 == NULL) || (VP2 == NULL) ) {
            plasma_error("PLASMA_sgebrd", "plasma_shared_alloc() failed");
            plasma_shared_free(plasma, TAUQ2);
            plasma_shared_free(plasma, VQ2);
            return PLASMA_ERR_OUT_OF_RESOURCES;
        }
        memset(TAUP2, 0, 2*MINMN*sizeof(float));
        memset(VP2,   0, 2*MINMN*sizeof(float));
    }
    /*=======================================
     *  calling bulge chasing
     *=======================================*/
    #if defined(ENABLE_TIMER)
    timeblg  = PLASMA_Wtime();
    #endif
    plasma_parallel_call_16(plasma_psgebrd_gb2bd_v1,
        PLASMA_enum,         uplo,
        int,                 MINMN,
        int,                 NB,
        int,                 Vblksiz,
        float*, AB,
        int,                 LDAB,
        float*, VQ2,
        float*, TAUQ2,
        float*, VP2,
        float*, TAUP2,
        float*,             S,
        float*,             E,
        int,                 WANTZ,
        int,                 WANTZ,
        PLASMA_sequence*,    sequence,
        PLASMA_request*,     request);
    /* WARNING: If plasma_psgebrd_gb2bd is implemented through a dynamic call, don't
     * forget to synchronize */
    plasma_dynamic_sync();
    #if defined(ENABLE_TIMER)
    timeblg   = PLASMA_Wtime()-timeblg;
    printf("  Finish Bulge       timing= %lf \n" ,timeblg);
    #endif
    /*=======================================
     *  END of calling bulge chasing
     *=======================================*/
    /*=======================================
     *  generate U from the bulge
     *=======================================*/
    if (jobu == PlasmaVec){
        memset(U,   0, M*LDU*sizeof(float));
        /* Initialize U to Identity */
        for(i=0; i<M; i++){
            U[i+i*LDU] = 1.0;
        }
        #if defined(ENABLE_TIMER)
        timeT  = PLASMA_Wtime();
        #endif
        /* compute T2 */
        plasma_static_call_8(plasma_pslarft_blgtrd,
            int,                 MINMN,
            int,                 NB,
            int,                 Vblksiz,
            float*, VQ2,
            float*, TQ2,
            float*, TAUQ2,
            PLASMA_sequence*,    sequence,
            PLASMA_request*,     request);
        #if defined(ENABLE_TIMER)
        plasma_dynamic_sync();
        timeT  = PLASMA_Wtime()-timeT;
        printf("  Finish compute TU2  timing= %lf \n" ,timeT);
        timeaplQ2  = PLASMA_Wtime();
        #endif
        /* apply Q2 from Left */
        plasma_static_call_14(plasma_psormqr_blgtrd,
            PLASMA_enum,         PlasmaLeft,
            PLASMA_enum,         PlasmaNoTrans,
            int,                 MINMN,
            int,                 NB,
            int,                 MINMN,
            int,                 Vblksiz,
            int,                 WANTZ,
            float*, VQ2,
            float*, TQ2,
            float*, TAUQ2,
            float*, U,
            int,                 LDU,
            PLASMA_sequence*,    sequence,
            PLASMA_request*,     request);
        #if defined(ENABLE_TIMER)
        plasma_dynamic_sync();
        timeaplQ2  = PLASMA_Wtime()-timeaplQ2;
        printf("  Finish compute U2  timing= %lf \n" ,timeaplQ2);
        #endif
        /*=======================================
         *  apply Q1 from the reduction to band
         *=======================================*/
        /* CASE NB>N, Q1 doesn't need to be applied, only bulge chasing has been done */
        if( NB < N ){
            #if defined(ENABLE_TIMER)
            plasma_dynamic_sync();
            timeconv1   = PLASMA_Wtime();
            #endif
            if ( PLASMA_TRANSLATION == PLASMA_OUTOFPLACE ) {
                plasma_sooplap2tile( descU, U, NB, NB, LDU, M, 0, 0, M, M, sequence, request, plasma_desc_mat_free(&(descU)) );
            } else {
                plasma_siplap2tile( descU, U, NB, NB, LDU, M, 0, 0, M, M, sequence, request);
            }
            #if defined(ENABLE_TIMER)
            timeconv1    = PLASMA_Wtime()-timeconv1;
            timeaplQ1   = PLASMA_Wtime();
            #endif
            /* Accumulate the transformations from the first stage */
            if(M<N){
                plasma_dynamic_call_7(plasma_psormqr,
                    PLASMA_enum, PlasmaLeft,
                    PLASMA_enum, PlasmaNoTrans,
                    PLASMA_desc, plasma_desc_submatrix(descA, descA.mb, 0, descA.m-descA.mb, descA.n-descA.nb),
                    PLASMA_desc, plasma_desc_submatrix(descU, descU.mb, 0, descU.m-descU.mb, descU.n),
                    PLASMA_desc, plasma_desc_submatrix(descT, descT.mb, 0, descT.m-descT.mb, descT.n-descT.nb),
                    PLASMA_sequence*, sequence,
                    PLASMA_request*, request);
            }
            else {
                plasma_dynamic_call_7(plasma_psormqr,
                    PLASMA_enum, PlasmaLeft,
                    PLASMA_enum, PlasmaNoTrans,
                    PLASMA_desc, descA,
                    PLASMA_desc, descU,
                    PLASMA_desc, descT,
                    PLASMA_sequence*, sequence,
                    PLASMA_request*, request);
            }
            #if defined(ENABLE_TIMER)
            plasma_dynamic_sync();
            timeaplQ1   = PLASMA_Wtime()-timeaplQ1;
            printf("  Finish compute U1  timing= %lf \n" ,timeaplQ1);
            timeconv2    = PLASMA_Wtime();
            #endif
            if ( PLASMA_TRANSLATION == PLASMA_OUTOFPLACE ) {
                plasma_sooptile2lap( descU, U, NB, NB, LDU, M, sequence, request );
                plasma_dynamic_sync();
                plasma_desc_mat_free(&descU);
            } else {
                plasma_siptile2lap( descU, U, NB, NB, LDU, M, sequence, request );
                plasma_dynamic_sync();
            }
            #if defined(ENABLE_TIMER)
            plasma_dynamic_sync();
            timeconv2    = PLASMA_Wtime()-timeconv2;
            printf("  Finish convert  U  timing= %lf \n" ,timeconv1+timeconv2);
            #endif
        } /* END of ( NB < N ) */
    }
    /*=======================================
     *  END of calling computing U
     *=======================================*/
    /*=======================================
     *  generate VT from the bulge
     *=======================================*/
    if (jobvt == PlasmaVec){
        memset(VT,   0, N*LDVT*sizeof(float));
        /* Initialize VT to Identity */
        for(i=0; i<N; i++){
            VT[i+i*LDVT] = 1.0;
        }
        #if defined(ENABLE_TIMER)
        timeT  = PLASMA_Wtime();
        #endif
        /* compute T2 */
        plasma_static_call_8(plasma_pslarft_blgtrd,
            int,                 MINMN,
            int,                 NB,
            int,                 Vblksiz,
            float*, VP2,
            float*, TP2,
            float*, TAUP2,
            PLASMA_sequence*,    sequence,
            PLASMA_request*,     request);
        #if defined(ENABLE_TIMER)
        plasma_dynamic_sync();
        timeT  = PLASMA_Wtime()-timeT;
        printf("  Finish compute TV2  timing= %lf \n" ,timeT);
        timeaplQ2  = PLASMA_Wtime();
        #endif
        /* apply Q2 from Left */
        plasma_static_call_14(plasma_psormqr_blgtrd,
            PLASMA_enum,         PlasmaRight,
            PLASMA_enum,         PlasmaTrans,
            int,                 MINMN,
            int,                 NB,
            int,                 MINMN,
            int,                 Vblksiz,
            int,                 WANTZ,
            float*, VP2,
            float*, TP2,
            float*, TAUP2,
            float*, VT,
            int,                 LDVT,
            PLASMA_sequence*,    sequence,
            PLASMA_request*,     request);
        #if defined(ENABLE_TIMER)
        plasma_dynamic_sync();
        timeaplQ2  = PLASMA_Wtime()-timeaplQ2;
        printf("  Finish compute V2 timing= %lf \n" ,timeaplQ2);
        #endif
        /*=======================================
         *  apply Q1 from the reduction to band
         *=======================================*/
        /* CASE NB>N, Q1 doesn't need to be applied, only bulge chasing has been done */
        if( NB < N ){
            #if defined(ENABLE_TIMER)
            plasma_dynamic_sync();
            timeconv1   = PLASMA_Wtime();
            #endif
            if ( PLASMA_TRANSLATION == PLASMA_OUTOFPLACE ) {
                plasma_sooplap2tile( descVT, VT, NB, NB, LDVT, N, 0, 0, N, N, sequence, request, plasma_desc_mat_free(&(descVT)) );
            } else {
                plasma_siplap2tile( descVT, VT, NB, NB, LDVT, N, 0, 0, N, N, sequence, request);
            }
            #if defined(ENABLE_TIMER)
            timeconv1    = PLASMA_Wtime()-timeconv1;
            timeaplQ1   = PLASMA_Wtime();
            #endif
            /* Accumulate the transformations from the first stage */
            if(M<N){
                plasma_dynamic_call_7(plasma_psormlq,
                    PLASMA_enum, PlasmaRight,
                    PLASMA_enum, PlasmaNoTrans,
                    PLASMA_desc, descA,
                    PLASMA_desc, descVT,
                    PLASMA_desc, descT,
                    PLASMA_sequence*, sequence,
                    PLASMA_request*, request);
            }
            else {
                plasma_dynamic_call_7(plasma_psormlq,
                    PLASMA_enum, PlasmaRight,
                    PLASMA_enum, PlasmaNoTrans,
                    PLASMA_desc, plasma_desc_submatrix(descA, 0, descA.nb, descA.m-descA.mb, descA.n-descA.nb),
                    PLASMA_desc, plasma_desc_submatrix(descVT,0, descVT.nb, descVT.m, descVT.n-descVT.nb),
                    PLASMA_desc, plasma_desc_submatrix(descT, 0, descT.nb, descT.m-descT.mb, descT.n-descT.nb),
                    PLASMA_sequence*, sequence,
                    PLASMA_request*, request);
            }
            #if defined(ENABLE_TIMER)
            plasma_dynamic_sync();
            timeaplQ1   = PLASMA_Wtime()-timeaplQ1;
            printf("  Finish compute V1  timing= %lf \n" ,timeaplQ1);
            timeconv2    = PLASMA_Wtime();
            #endif

            if ( PLASMA_TRANSLATION == PLASMA_OUTOFPLACE ) {
                plasma_sooptile2lap( descVT, VT, NB, NB, LDVT, N, sequence, request );
                plasma_dynamic_sync();
                plasma_desc_mat_free(&descVT);
            } else {
                plasma_siptile2lap( descVT, VT, NB, NB, LDVT, N, sequence, request );
                plasma_dynamic_sync();
            }
            #if defined(ENABLE_TIMER)
            plasma_dynamic_sync();
            timeconv2    = PLASMA_Wtime()-timeconv2;
            printf("  Finish convert VT  timing= %lf \n" ,timeconv1+timeconv2);
            #endif
        } /* END of ( NB < N ) */
    }
    /*=======================================
     *  END of calling computing VT
     *=======================================*/
    #if defined(ENABLE_TIMER)
    timeall = PLASMA_Wtime()-timeall;
    printf("  Finish full Bidiagonalisation threads %d  N %d  timeall= %lf \n", plasma->world_size, N, timeall);
    #endif
    if( jobu  == PlasmaVec )
        plasma_shared_free(plasma, TQ2);
    if( jobvt == PlasmaVec )
        plasma_shared_free(plasma, TP2);
    plasma_shared_free(plasma, VQ2);
    plasma_shared_free(plasma, TAUQ2);
    plasma_shared_free(plasma, VP2);
    plasma_shared_free(plasma, TAUP2);
    plasma_shared_free(plasma, AB);
    return PLASMA_SUCCESS;
}

ext: c