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

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_combination.h>

gsl_combination *
gsl_combination_alloc (const size_t n, const size_t k)
{
  gsl_combination * c;

  if (n == 0)
    {
      GSL_ERROR_VAL ("combination parameter n must be positive integer",
                        GSL_EDOM, 0);
    }
  if (k > n)
    {
      GSL_ERROR_VAL ("combination length k must be an integer less than or equal to n",
                        GSL_EDOM, 0);
    }
  c = (gsl_combination *) malloc (sizeof (gsl_combination));

  if (c == 0)
    {
      GSL_ERROR_VAL ("failed to allocate space for combination struct",
                        GSL_ENOMEM, 0);
    }

  if (k > 0)
    {
      c->data = (size_t *) malloc (k * sizeof (size_t));

      if (c->data == 0)
        {
          free (c);             /* exception in constructor, avoid memory leak */

          GSL_ERROR_VAL ("failed to allocate space for combination data",
                         GSL_ENOMEM, 0);
        }
    }
  else
    {
      c->data = 0;
    }

  c->n = n;
  c->k = k;

  return c;
}

gsl_combination *
gsl_combination_calloc (const size_t n, const size_t k)
{
  size_t i;

  gsl_combination * c =  gsl_combination_alloc (n, k);

  if (c == 0)
    return 0;

  /* initialize combination to identity */

  for (i = 0; i < k; i++)
    {
      c->data[i] = i;
    }

  return c;
}

void
gsl_combination_init_first (gsl_combination * c)
{
  const size_t k = c->k ;
  size_t i;

  /* initialize combination to identity */

  for (i = 0; i < k; i++)
    {
      c->data[i] = i;
    }
}

void
gsl_combination_init_last (gsl_combination * c)
{
  const size_t k = c->k ;
  size_t i;
  size_t n = c->n;

  /* initialize combination to identity */

  for (i = 0; i < k; i++)
    {
      c->data[i] = n - k + i;
    }
}

void
gsl_combination_free (gsl_combination * c)
{
  if (c->k > 0) free (c->data);
  free (c);
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
/* Copyright 2016. The Regents of the University of California.
 * Copyright 2016. Martin Uecker.
 * All rights reserved. Use of this source code is governed by
 * a BSD-style license which can be found in the LICENSE file.
 *
 * Authors:
 * 2016 Jonathan Tamir <jtamir@eecs.berkeley.edu>
 * 2016 Martin Uecker <martin.uecker@med.uni-goettingen.de>
 */

#include <assert.h>
#include <complex.h>

#include "misc/misc.h"

#ifdef USE_MACPORTS
#include <cblas_openblas.h>
#else
#include <cblas.h>
#endif

#ifdef USE_CUDA
#include "num/gpuops.h"

#include <cublas.h>
#endif

#include "blas.h"

void blas_cgemm(char transa, char transb, long M, long N,  long K, const complex float alpha, long lda, const complex float A[K][lda], long ldb, const complex float B[N][ldb], const complex float beta, long ldc, complex float C[N][ldc])
{
#ifdef USE_CUDA
#define CUCOMPLEX(x) (((union { cuComplex cu; complex float std; }){ .std = (x) }).cu)
        if (cuda_ondevice(A)) {

                cublasCgemm(transa, transb, M, N, K,  CUCOMPLEX(alpha),
                                (const cuComplex*)A, lda,
                                (const cuComplex*)B, ldb, CUCOMPLEX(beta),
                                (cuComplex*)C, ldc);
        } else
#endif
        cblas_cgemm(CblasColMajor, transa, transb, M, N, K, (void*)&alpha, (void*)A, lda, (void*)B, ldb, (void*)&beta, (void*)C, ldc);
}

void (blas_matrix_multiply)(long M, long N, long K, complex float C[N][M], const complex float A[K][M], const complex float B[N][K])
{
	blas_cgemm(CblasNoTrans, CblasNoTrans, M, N, K, 1., M, A, K, B, 0., M, C);
}

void (blas_csyrk)(char uplo, char trans, long N, long K, const complex float alpha, long lda, const complex float A[][lda], complex float beta, long ldc, complex float C[][ldc])
{
	assert('U' == uplo);
	assert(('T' == trans) || ('N' == trans));

	cblas_csyrk(CblasColMajor, CblasUpper, ('T' == trans) ? CblasTrans : CblasNoTrans, N, K, (void*)&alpha, (void*)A, lda, (void*)&beta, (void*)C, ldc);
}

ext: c

max_stars_repo_name: youngwook-mri/bart

max_issues_repo_name: youngwook-mri/bart

max_forks_repo_name: youngwook-mri/bart


text:
#include <sys/types.h> // per utilizzare variabili di tipo pthread_t
#include <stdlib.h> // utilità generale (allocazione, controllo processi, etc.)
#include <pthread.h> // per utilizzo dei thread
#include <string.h> // per elaborazione su stringhe
#include <math.h>

#include <gsl/gsl_vector.h> // per usare i vettori
#include <gsl/gsl_matrix.h> // per usare le matrici

#include <gsl/gsl_eigen.h> // per calcolare autovalori e autovettori

#include <gsl/gsl_statistics.h> // per usare funzioni di probabilità

#include <gsl/gsl_linalg.h> // per usare funzioni su matrici
#include <gsl/gsl_permutation.h> // per usare le permutazioni

#define debug 0 // variabile per settare le stampe di controllo
#define troncamento 100 // per approssimazione alla 3° cifra significativa dopo la virgola

int controllo = 1; // variabile per verificare la chiamata di funzioni

double **data_set_origine, **data_set, **data_compression, **autovettori;
double *autovalori;
int dimensione = 0;

struct registro{ // puntatori alle celle per eseguire il prodotto matriciale A*B
	double *mtx_1, *mtx_2; //mtx_1: riga di A, mtx_2: colonna di B
	int *indici; // posizione in cui salvare il singolo valore dell'operazione
} catalogo;

struct gsl{ // puntatori agli autovalori e autovettori di tipo gsl
	gsl_vector *vector; // autovalori
	gsl_matrix *matrix; // autovettori
} puntatori_gsl;

void stampa(int tipo, double **mt1, double *mt2, char nome[], int nr, int nc){ // il tipo indica che tipo di stampa effettuo
	int i,j = 0;
	switch(tipo){
		case 0: //stampa dimensione
			printf("%s: %d \n", nome, nr); //Dimensioni
		break;
		case 1: //stampa la matrice con il proprio nome
			printf("Stampa matrice %s: \n", nome);
			for(i = 0; i < nr; i++){
				for(j = 0; j < nc; j++){
					if(j == 0){
						printf("|   ");
					}
					printf("%2.2f\t", mt1[i][j]);
					if(j == nc-1){
						printf("|");
					}
				}
				printf("\n");
			}
		break;
		case 2://stampa media in forma estesa
			printf("Stampa valor medio: \n");
			for(i = 0; i < nr; i++){
				printf("media[%d]=%2.2f \n", i, mt2[i]);
			}
		break;
		case 3: //stampa gli autovalori in forma matriciale
			printf("Stampa autovalori(covarianza): \n");
			for(i = 0; i < dimensione; i++){
				if(i == 0){
					printf("|   ");
				}
				printf("%2.4g\t", autovalori[i]);
				if(i == dimensione-1){
					printf("| \n");
				}
			}
		break;
		case 5: //stampa errori verifica
			printf("Errore la verifica non è andata a buon fine\n");
			printf("La matrice Verifica è diversa dalla matrice di partenza \n");
		break;
		case 6://stampa matrici monodimensionali in maniera estesa
			printf("Stampa Matrice Estesa %s: \n", nome);
			for(i = 0; i < dimensione; i++){
				printf("%s[%d]=%2.2f \n", nome, i, mt2[i]);
			}
		break;
		case 7: //stampa matrici bidimensionali in maniera estesa
			printf("Stampa Matrice Estesa %s : \n", nome);
			for(i = 0; i < dimensione; i++){
				for(j = 0; j < dimensione; j++){
					printf("%s[%d][%d]=%2.2f \n", nome, i, j, mt1[i][j]);
				}
			}
		break;
	}
	printf("\n");
}

void *thread_function(void *candidato){ // singola iterazione del prodotto matriciale
	int i = 0;
    int *indici;
	double *mtx_1, *mtx_2;
	struct registro *catalogo;

	// casting degli argomenti della thread_function
	catalogo=(struct registro*)candidato;
	mtx_1=catalogo->mtx_1;
	mtx_2=catalogo->mtx_2;
	indici=catalogo->indici;

	/*
		esecuzione del prodotto matriciale (a singoli elementi): si esegue il singolo prodotto, il risultato
		lo si salva sulla posizione corrente del'array tmp[0]; infine si somma il precedente risultato al valore presente
		nella cella 0 di tmp (tmp[0][0] svolgerà il ruolo di sommatore, senza dover creare una variabile temporanea).
	*/
	for(i=0; i<dimensione; i++){
		if(debug==1){printf("mtx_1[%d] = %2.2f  ,  mtx_2[%d] = %2.2f\n", i, mtx_1[i], i, mtx_2[i]);}
		mtx_1[i]=mtx_1[i]*mtx_2[i];
		if(debug==1){printf("singolo prodotto: %2.2f\n", mtx_1[i]);}
		if(i>0){
			mtx_1[0]=mtx_1[0]+mtx_1[i];
			if(debug==1){printf("somma: %2.2f\n", mtx_1[0]);}
		}
	} i=0;

	// inserimento della somma nell'apposita cella della matrice globale: data_compression
	data_compression[indici[0]][indici[1]]=mtx_1[0];
	if(debug==1){
		printf("data_compression[%d][%d]=%2.2f\n", indici[0], indici[1], data_compression[indici[0]][indici[1]]);
	}

	return NULL;
}

double *media(int nrighe){ // media sulle righe del data set
	int i, j = 0;
	double somma = 0;
	double *media;

	media=malloc(nrighe*sizeof(double)); // alloccamento dinamico
	if(media==NULL){printf("ERRORE: media non alloccata.\n");}

	for(i=0; i<nrighe; i++){
		for(j=0; j<dimensione; j++){
			somma = somma+data_set[i][j]; // somma su riga i-esima del data_set
			if(debug==1){printf("somma = %2.2f\n", somma);}
		}
		media[i]=somma/dimensione; // valor medio della i-esima riga del data_set
		somma=0; // reset per riga successiva
		if(debug){printf("Media: %2.2f \n", media[i]);}
	}

	return media;
}

double **covarianza(double *media, int nrighe){ // matrice di covarianza del data set
	int i, j, k = 0;
	double **covarianza;
	double	*tmp1, *tmp2;

	// alloccamento dinamico
	tmp1=malloc(nrighe*sizeof(double));
	if(tmp1==NULL){printf("ERRORE - covarianza: tmp1 non alloccato.\n"); exit(0);}

	tmp2=malloc(nrighe*sizeof(double));
	if(tmp2==NULL){printf("ERRORE - covarianza: tmp2 non alloccato.\n"); exit(0);}

	covarianza=malloc(dimensione*sizeof(double*));
	if(covarianza==NULL){printf("ERRORE - covarianza: covarianza non alloccato.\n"); exit(0);}

	for(i=0; i<dimensione; i++){
		covarianza[i]=malloc(dimensione*sizeof(double));
		if(covarianza[i]==NULL){printf("ERRORE - covarianza: covarianza[%d] non alloccato.\n", i); exit(0);}
	}

	/*
		dei 4 for, definiamo "principali" quelli riferiti alla variabile 'i' e 'j' in quanto vanno ad identificare la cella  (i, j) della matrice di
		covarianza; i for riferiti alla variabile 'k' sono for di settaggio dei due array temporanei, che andranno a considerare le colonne
		del data set.
		Ad ogni ciclo dell' i-for si fisserà una riga (partendo dalla riga 0), e si andrà a "combinare" con tutte le colonne del data set (a
		partire dalla riga 0).
	*/

	for(i=0; i<dimensione; i++){
		for(k=0; k<nrighe; k++){ // fissiamo colonna i-esima
			tmp1[k]=data_set[k][i];
			if(debug==1){
				stampa(6, NULL, tmp1, "tmp1", 0, 0);
			}
		}
		for(j=0; j<dimensione; j++){
			for(k=0; k<nrighe; k++){ // consideriamo colonna j-esima
				tmp2[k]=data_set[k][j];
				if(debug==1){
					stampa(6, NULL, tmp2, "tmp2", 0, 0);
				}
			}
			covarianza[i][j]=gsl_stats_covariance(tmp1, 1, tmp2, 1, nrighe); // calcolo della covarianza
		}
		if(debug==1){
			stampa(7, covarianza, NULL, "Covarianza", 0, 0);
		}
	}

	free(tmp1);
	free(tmp2);

	return covarianza;
}

struct gsl eigen(double **covarianza){ // calcolo autovalori e autovettori della matrice di covarianza
	int i=0, j=0;
	gsl_matrix *matrice, *evec;
	gsl_vector *evel;
	gsl_eigen_symmv_workspace *campo;
	struct gsl pointer_gsl;

	// allocazione
	matrice=gsl_matrix_alloc(dimensione, dimensione); // alloccamento matrice bid: "dimensione" righe e "dimensione" colonne
	if(matrice==NULL){printf("ERRORE - eigen: matrice non alloccato"); exit(0);}

	evec=gsl_matrix_alloc(dimensione, dimensione); // alloccamento matrice bid: "dimensione" righe e "dimensione" colonne
	if(evec==NULL){printf("ERRORE - eigen: evec non alloccato"); exit(0);}

	evel=gsl_vector_calloc(dimensione); // alloccamento vettore: "dimensione" posizioni
	if(evel==NULL){printf("ERRORE - eigen: evel non alloccato"); exit(0);}

	campo=gsl_eigen_symmv_alloc(dimensione); // alloccamento del workspace, dimensione pari a: "dimensione"
	if(campo==NULL){printf("ERRORE - eigen: campo non alloccato"); exit(0);}

	// copiamo "covarianza" nella "matrice"
	for(i=0; i<dimensione; i++){
		for(j=0; j<dimensione; j++){
			if(debug==1){printf("covarianza[%d][%d] = %2.2f\n", i, j, covarianza[i][j]);}
			gsl_matrix_set(matrice, i, j, covarianza[i][j]); // inserimento del valore nella "matrice"
			if(debug==1){printf("matrice[%d][%d] = %2.2f\n", i, j, gsl_matrix_get(matrice, i, j));}
		}
	}

	controllo=gsl_eigen_symmv(matrice, evel, evec, campo); // calcolo degli autovalori e autovettori
	if(controllo!=0){printf("ERRORE - eigen: gsl_eigen_symmv fallito.\n"); exit(0);}

	// salvataggio dei puntatori agli autovalori  e autovettori nella struttura "gsl"
	pointer_gsl.vector=evel;
	pointer_gsl.matrix=evec;

	gsl_matrix_free(matrice);
	gsl_eigen_symmv_free(campo);

	return pointer_gsl;
}

void ordinamento(gsl_matrix *vettori){ // ordinamento decrescente degli autovalori
	int i,j = 0;
	double tmp = 0;

	for(i=0; i<dimensione; i++){
		if(autovalori[i]<0){
			autovalori[i]=autovalori[i]*(-1);
		}
	}

	/*
		dato che c'è una corrispondenza biunivoca fra l'array degli autovalori e la matrice degli autovettori, allora
		ordinando gli autovalori si ordinano anche gli autovettori, per mantenere la precedente proprietà.
	*/
	for(i = 0; i < dimensione; i++){
		for(j = 0; j < dimensione; j++){
			if(debug==1){printf("autovalori[%d] = %2.2f  ,  autovalori[%d] = %2.2f\n", j, autovalori[j], i, autovalori[i]);}
			if(autovalori[j] < autovalori[i]){
				tmp=autovalori[i];
				if(debug==1){printf("tmp = %2.2f\n", tmp);}
				autovalori[i]=autovalori[j];
				autovalori[j]=tmp;
				if(debug==1){printf("autovalori[%d] = %2.2f  ,  autovalori[%d] = %2.2f\n", j, autovalori[j], i, autovalori[i]);}
				gsl_matrix_swap_columns(vettori, i, j);
				if(controllo!=0){printf("ERRORE - ordinamento: gsl_matrix_swap_columns fallito.\n"); exit(0);}
			}
		}
		tmp=0;
	}
}

double **sottrazione_matriciale(int nrighe, double *media){ // data set "differenziato"
	int i, j = 0;
	double **sottrazione;

	// alloccamento dinamico
	sottrazione=malloc(nrighe*sizeof(double*));
	if(sottrazione==NULL){printf("ERRORE - sottrazione_matriciale: sottrazione non alloccato.\n"); exit(0);}
	for(i=0; i<nrighe; i++){
		sottrazione[i]=malloc(dimensione*sizeof(double));
		if(sottrazione[i]==NULL){printf("ERRORE - sottrazione_matriciale: sottrazione[%d] non alloccato.\n", i); exit(0);}
	}

	// ogni riga del data set viene differenziato per la media della riga stessa
	for(i=0; i<nrighe; i++){
		for(j=0; j<dimensione; j++){
			sottrazione[i][j]=data_set[i][j]-media[i];
			if(debug==1){printf("data_set[%d][%d] = %2.2f  ,  media[%d] = %2.2f  ,  sottrazione[%d][%d] = %2.2f\n",
					i, j, data_set[i][j], i, media[i], i, j, sottrazione[i][j]);
			}
		}
	}

	return sottrazione;
}

int prodotto_matriciale(int righe_1, int colonne_1, int righe_2, int colonne_2, double **matrice_1, double **matrice_2){ //Funzione per fare il prodotto matriciale
	int i, i_1, j_1, contatore, fine = 0;
	int *indici;
	double *mtx_1, *mtx_2;
	pthread_t *thread;

	struct registro *enciclopedia; // array di strutture registro

	if(colonne_1 != righe_2){
		printf("ERRORE - prodotto_matriciale: non è soddisfatta la condizione del prodotto matriciale.\n");
		return -1;
	}

	enciclopedia=malloc(righe_1*colonne_2*sizeof(struct registro)); //creazione struttura
	if(enciclopedia==NULL){printf("ERRORE - prodotto_matriciale: enciclopedia non alloccata.\n"); exit(0);}

	/*
		Definiamo come "prodotto", il prodotto dei valori (in stessa posizione) fra due array, e successiva somma dei singoli
		prodotti.
		Definiamo come "singola iterazione" o "singola procedura" il prodotto fra la riga i-esima della matrice A per la colonna
		j-esima della matrice B, al variare di i e j.
	*/

		// settaggio dei parametri da inviare ai thread (obbligatori, altrimenti parte il segmentation fault)
		contatore=0;
		i_1=0;
		j_1=0;

		/*
			Il seguente while serve a riempire l'array di strutture che contengono: la riga di A, la colonna di B e gli indici di posizione
			in cui salvare il risultato della singola procedura del prodotto matriciale.
		*/
		while(fine==0){ // alloccamento dinamico di due array "temporanei" e salvataggio indirizzi nell'array tmp

			mtx_1=malloc(colonne_1*sizeof(double)); // riga di A
			if(mtx_1==NULL){printf("ERRORE - prodotto_matriciale: mtx_1 non alloccato.\n"); exit(0);}

			mtx_2=malloc(righe_2*sizeof(double*)); // colonna di B
			if(mtx_2==NULL){printf("ERRORE - prodotto_matriciale: mtx_2 non alloccato.\n"); exit(0);}

			indici=malloc(2*sizeof(int)); // array degli indici di posizione del data_compression
			if(indici==NULL){printf("ERRORE - prodotto_matriciale: indici non alloccato.\n"); exit(0);}

			if(debug==1){
				printf("contatore = %d, colonne_2=%d, righe_1=%d, i_1=%d, j_1=%d\n",
																	contatore, colonne_2, righe_1, i_1, j_1);
				}

			for(i = 0; i < colonne_1; i++){
				mtx_1[i]=matrice_1[i_1][i]; // copia della riga i_1 di A in mtx_1
				mtx_2[i]=matrice_2[i][j_1]; // copia della colonna j_1 di B in mtx_2
			}

			if(debug==1){
				stampa(6, NULL, mtx_1, "MTX_1", 0, 0);
				stampa(6, NULL, mtx_2, "MTX_2", 0, 0);
			}

			// salvataggio dimensioni ed indici di posizione nella singola struttura
			indici[0]=i_1;
			indici[1]=j_1;

			// salvataggio dei parametri in una struttura "registro"
			enciclopedia[contatore].mtx_1 = mtx_1;
			enciclopedia[contatore].mtx_2 = mtx_2;
			enciclopedia[contatore].indici = indici;

			contatore++; // rappresenta l'indice dell'array di struttura
			if(debug==1){
				printf("contatore = %d\n", contatore);
			}

			if(debug==1){
				printf("j_1+1 = %d  ,  colonne_2 = %d  ,  i_1+1 = %d  ,  righe_1 = %d\n", j_1+1, colonne_2, i_1+1, righe_1);
				}
			if(j_1+1<colonne_2){ // andiamo alla colonna successiva di B
				j_1++;
				if(debug==1){
					printf("j_1 = %d\n", j_1);
				}
			}else if(i_1+1<righe_1){ // andiamo alla riga successiva di A
				i_1++;
				j_1=0; // la riga parte dalla colonna 0
				if(debug==1){
					printf("i_1 = %d\n", i_1);
				}
			}else{
				fine=1; // sono stati passati tutti i valori, il while è finito
			}
		}
		contatore=0;
		i=0;
		i_1=0;
		j_1=0;
		fine=0;

		// alloccamento array di thread
		thread=malloc(righe_1*colonne_2*sizeof(pthread_t));
		if(thread==NULL){
			printf("ERRORE - prodotto_matriciale: thread non alloccato.\n");
			exit(0);
		}

		// verifica del corretto contenuto delle strutture
		if(debug==1){
			int j = 0;
			for(i=0; i<righe_1*colonne_2; i++){
				printf("Stampa del contenuto delle strutture.\n");
				printf("------------\n");
				for(j=0; j<colonne_1; j++){
					printf("enciclopedia[%d].mtx_1[%d]= %2.2f\n", i, j, enciclopedia[i].mtx_1[j]);
				}
				for(j=0; j<colonne_1; j++){
					printf("enciclopedia[%d].mtx_2[%d]= %2.2f\n", i, j, enciclopedia[i].mtx_2[j]);
				}
				for(j=0; j<2; j++){
					printf("enciclopedia[%d].indici[%d]=%d\n", i, j, enciclopedia[i].indici[j]);
				}
			}
			printf("------------\n");
		}

		for(i=0; i<righe_1*colonne_2; i++){

			/*
			chiamata del singolo thread, gli argomenti sono: riga della matrice 1, colonna della matrice 2, indice della riga considerata e
			indice della  colonna considerata.
			*/
			//printf("------------\n");
			controllo=pthread_create(&thread[i], NULL, (void*)thread_function, (void*)&enciclopedia[i]);
			if(controllo!=0){printf("ERRORE - prodotto_matriciale: pthread_create fallito.\n"); exit(0);}
			//printf("------------\n");
		}
		controllo=pthread_join(thread[(righe_1*colonne_2)-1], NULL); // il main attende l'ultimo thread prima di proseguire
		if(controllo!=0){printf("ERRORE - prodotto_matriciale: pthread_join fallito.\n"); exit(0);}

		free(mtx_1);
		free(mtx_2);
		free(indici);
		free(enciclopedia);

		return 0;
}

double **inversione_matriciale(gsl_matrix *evec){ // matrice degli autovettori invertita
	int i, j = 0;
	int *segno;
	double **inversione;
	gsl_matrix *inverse;
	gsl_permutation *permutazione;

	// alloccamento dinamico
	inversione=malloc(dimensione*sizeof(double*));
	if(inversione==NULL){printf("ERRORE - inversione_matriciale: inversione non alloccato.\n"); exit(0);}
	for(i=0; i<dimensione; i++){
		inversione[i]=malloc(dimensione*sizeof(double));
		if(inversione[i]==NULL){printf("ERRORE - inversione_matriciale: inversione[%d] non alloccato.\n", i); exit(0);}
	}

	inverse=gsl_matrix_alloc(dimensione, dimensione); //alloccamento della matrice inversa
	if(inverse==NULL){printf("ERRORE - inversione_matriciale: inverse non alloccato.\n"); exit(0);}

	permutazione=gsl_permutation_calloc(dimensione); // alloccamento dell'array di permutazione
	if(permutazione==NULL){printf("ERRORE - inversione_matriciale: permutazione non alloccato.\n"); exit(0);}

	segno=malloc(1*sizeof(int)); // inizializzazione della variabile
	if(segno==NULL){printf("ERRORE - inversione_matriciale: segno non alloccato.\n"); exit(0);}

	controllo=gsl_linalg_LU_decomp(evec, permutazione, segno); //decompoposizione (obbligatorio)
	if(controllo!=0){printf("ERRORE - inversione_matriciale: gsl_linalg_LU_decomp fallito.\n");}

	controllo=gsl_linalg_LU_invert(evec, permutazione, inverse); //inversione matriciale
	if(controllo!=0){printf("ERRORE - inversione_matriciale: gsl_linalg_LU_invert fallito.\n");}

	// matrice inversa inserita in una matrice "non di libreria"
	for(i=0; i<dimensione; i++){
		for(j=0; j<dimensione; j++){
			inversione[i][j]=gsl_matrix_get(inverse, i, j);
		}
	}
	gsl_matrix_free(inverse);
	gsl_permutation_free(permutazione);
	free(segno);

	return inversione;
}

double **addizione_matriciale(int nrighe, double *media){ // data compression "addizionato"
	int i, j = 0;
	double **addizione;

	// alloccamento dinamico
	addizione=malloc(nrighe*sizeof(double*));
	if(addizione==NULL){printf("ERRORE - addizione_matriciale: addizione non alloccato.\n"); exit(0);}
	for(i=0; i<nrighe; i++){
		addizione[i]=malloc(dimensione*sizeof(double));
		if(addizione[i]==NULL){printf("ERRORE - addizione_matriciale: addizione[%d] non alloccato.\n", i); exit(0);}
	}

	// ogni riga del data compression viene addizzionato per la media (della stessa riga) del data set
	for(i=0; i<nrighe; i++){
		for(j=0; j<dimensione; j++){
			addizione[i][j]=data_compression[i][j]+media[i];
			if(debug==1){printf("data_compression[%d][%d] = %2.2f  ,  media[%d] = %2.2f  ,  addizione[%d][%d] = %2.2f\n",
					i, j, data_compression[i][j], i, media[i], i, j, addizione[i][j]);
			}
		}
	}

	return addizione;
}

double **inversione_operazione(int nrighe, double *media, gsl_matrix *evec){ // invertiamo l'algoritmo
	double **verifica;

	// verifica = data_compression*autovettori(invertiti) + media_righe_data_set = data_set

	// operazione di inversione della matrice
	verifica=inversione_matriciale(evec);
	//stampa(1, verifica, NULL, "Inversa", dimensione, 0);

	// operazione di prodotto matriciale
	if(0 != prodotto_matriciale(nrighe, dimensione, dimensione, dimensione, data_compression, verifica)){
		printf("ERRORE - inversione_operazione: funzione prodotto_matriciale ha riportato un errore, esecuzione arrestata.\n");
		exit(0);
	}

	//stampa(1, data_compression, NULL, "Prodotto Matriciale", nrighe, 0);

	verifica=addizione_matriciale(nrighe, media); // operazione di somma matriciale

	return verifica;
}

void trasponi(int nr, int nc){ // esegue la trasposizione della matrice
	int i, j;
	j = 0;

	for(j = 0; j < nc; j++){
		for(i = 0; i < nr; i++){
				data_set[j][i] = data_set_origine[i][j];
			}
	}
}

void ottieni(char stringa[], int riga){ // estrazione del data set
	int set=0, last_set=0;
	int i, j=0, k;
	char buf[100];
	double a;

	k=0;
	for(i = 0; i < strlen(stringa); i++){ //per leggere tutti i caratteri
		if( (((int)stringa[i]) == 43) || (((int)stringa[i]) == 45) || (((int)stringa[i]) == 46) || ((((int)stringa[i]) > 47) && (((int)stringa[i]) < 58))){
			set=1; //dico che erano caratteri
			buf[j]=stringa[i]; //salvo il numero
			j++; //incremento il puntatore
		}else{
			set=0;
		}
		if((set==0 && last_set==1) || (set==1 && i == (strlen(stringa)-1))){ //controllo se ho letto tutto il numero
			a=atof(buf); //converto la stringa in un double
			data_set_origine[riga][k]=a; //salvo il numero
			k++; //incremento la colonna
			strcpy(buf, ""); //svuoto il buf
			j=0;//pongo il puntatore del buf a 0 così da poter leggere il numero seguente
		}
		last_set=set; //dico qual'era il valore precedente
	}
}

int colonne(char stringa[]){ //conto il numero di colonne
	int set=0, last_set=0; //variabili per puntare il carattere e il penultimo carattere
	int i, j=0; //variabili locali, j è il numero di colonne

	for(i = 0; i < strlen(stringa); i++){
		/*controllo se il carattere in questione è un numero, il segno + o -, oppure il .*/
		if( (((int)stringa[i]) == 43) || (((int)stringa[i]) == 45) || (((int)stringa[i]) == 46) || ((((int)stringa[i]) > 47) && (((int)stringa[i]) < 58))){
			set=1;
		}else{
			set=0;
		}
		if((last_set == 1 && set == 0)||(set == 1 && i == (strlen(stringa) -1))){ //controllo se ho finito di leggere il numero
			j++;
		}
		last_set=set;
	}
	return j; //ritorno il numero di colonne
}

int lettura(){ // funzione per leggere il file in input
	FILE *fp; // file per la lettura
	char pathname[1000]; //percorso file
	char buf[2000];
	char *s;
	int i=0, nrighe=0, ncolonne=1;
	int j=0;

	/*prendo in ingresso il nome del file*/
	printf("Inserisci l'indirizzo del file: \n");
	scanf("%s", pathname);
	if( (fp = fopen(pathname, "r")) == NULL ){ // apro il file e controllo se è avvenuto con successo
		printf("ERRORE - lettura: apertura file %s fallito.\n", pathname);
		exit(0);
	}

	/*inizio la lettura del file*/
	while(1){ // conto il numero di righe e colonne
		stpcpy(buf, "");
		s = fgets(buf, sizeof(buf), fp);
		if(s == NULL) break; // se la linea è vuota ho letto tutto il file
		nrighe++; // conto le righe
		if(j == 0){
			ncolonne=colonne(buf);
			j++;
		}
	}

	data_set_origine=malloc(nrighe*sizeof(double*));//allocco la memoria del data set
	data_set=malloc(ncolonne*sizeof(double*));

	if(data_set == NULL){
		printf("ERRORE - lettura: data_set non alloccato.\n");
		exit(0);
	}

	if(data_set_origine == NULL){
		printf("ERRORE - lettura: data_set non alloccato.\n");
		exit(0);
	}

	for(i = 0; i < nrighe; i++){//alloco tutta la memoria
		data_set_origine[i]=malloc(ncolonne*sizeof(double));
		data_set[i]=malloc(nrighe*sizeof(double));

		if(data_set[i] == NULL){
			printf("ERRORE - lettura: data_set[%d] non alloccato.\n", i);
			exit(0);
		}

		if(data_set_origine[i] == NULL){
			printf("ERRORE - lettura: data_set_origine[%d] non alloccato.\n", i);
			exit(0);
		}
	}

	rewind(fp);

	for(i = 0; i < nrighe; i++){
		stpcpy(buf, "");
		s = fgets(buf, sizeof(buf), fp);
		ottieni(buf, i);
	}

	fclose(fp);

	trasponi(nrighe, ncolonne);

	dimensione = nrighe;
	return ncolonne;
}

int main(int argc, char **argv){
	int nrighe, i, j = 0; //variabili di gestione
	double *valor_medio, **matrice_covarianza, **sottrazione, **verifica; //array per il salvataggio dei dati
	struct gsl pointer_gsl; //inizializzazione struttura

	if(debug){printf("MODALITA' DEBUG ATTIVA\n");}

	nrighe=lettura(); // apertura e lettura del data set, nrighe=ncolonne perchè abbiamo trasposto la matrice

	stampa(0, NULL, NULL, "nrighe", nrighe, 0);
	stampa(0, NULL, NULL, "Dimensione del Data_set", dimensione, 0); // stampa dimensioni del data set

	//stampa(1, data_set_origine, NULL, "Data Set", dimensione, nrighe);
	stampa(1, data_set, NULL, "Data Set Trasposto", nrighe, dimensione);

	// alloccamento dinamico
	data_compression=malloc(nrighe*sizeof(double*));
	if(data_compression==NULL){printf("ERRORE - main: data_compression non alloccato.\n"); exit(0);}
	for(i=0; i<nrighe; i++){
		data_compression[i]=malloc(dimensione*sizeof(double));
		if(data_compression[i]==NULL){printf("ERRORE - main: data_compression[%d] non alloccato.\n", i); exit(0);}
	}

	autovalori=malloc(dimensione*sizeof(double));
	if(autovalori==NULL){printf("ERRORE - main: autovalori non alloccato.\n"); exit(0);}

	autovettori=malloc(dimensione*sizeof(double*));
	if(autovettori==NULL){printf("ERRORE - main: autovettori non alloccato.\n"); exit(0);}
	for(i=0; i<dimensione; i++){
		autovettori[i]=malloc(dimensione*sizeof(double));
		if(autovettori==NULL){printf("ERRORE - main: autovettori[%d] non alloccato.\n", i); exit(0);}
	}

	valor_medio=media(nrighe); // calcolo dei valor medi
	//stampa(2, NULL, valor_medio, NULL, nrighe, 0);

	matrice_covarianza=covarianza(valor_medio, nrighe); // calcolo della matrice di covarianza
	//stampa(1, matrice_covarianza, NULL, "Covarianza", dimensione, dimensione);

	pointer_gsl=eigen(matrice_covarianza); // calcolo degli autovalori e autovettori

	for(i=0; i<dimensione; i++){ // settaggio degli autovalori nell'apposito array
		autovalori[i]=gsl_vector_get(pointer_gsl.vector, i);
	}

	//stampa(3, NULL, NULL, NULL, 0, 0); // stampa autovalori (non ordinanti)

	gsl_vector_free(pointer_gsl.vector);

	ordinamento(pointer_gsl.matrix); // ordinamento decrescente di autovalori e autovettori

	for(i=0; i<dimensione; i++){ // settaggio degli autovettori negll'apposito array
		for(j=0; j<dimensione; j++){
			autovettori[i][j]=gsl_matrix_get(pointer_gsl.matrix, i, j); // estrazione del valore da "evec"
			if(autovettori[i]==NULL){printf("ERRORE - main: gsl_matrix_get fallito.\n"); exit(0);}
		}
	}

	//stampa(3, NULL, NULL, NULL, 0, 0); // stampa autovalori
	//stampa(1, autovettori, NULL, "Autovettori", dimensione, dimensione); // stampa autovettori

	sottrazione=sottrazione_matriciale(nrighe, valor_medio); // operazione di sottrazione matriciale
	//stampa(1, sottrazione, NULL, "Sottrazione", nrighe, 0);

	if(0 != prodotto_matriciale(nrighe, dimensione, dimensione, dimensione, sottrazione, autovettori)){ // operazione di prodotto matriciale
		printf("ERRORE - main: funzione prodotto_matriciale ha riportato un errore, esecuzione arrestata.\n");
		return 0;
	}

	stampa(1, data_compression, NULL, "Data Compression", nrighe, dimensione);

	free(sottrazione);

	verifica=inversione_operazione(nrighe, valor_medio, pointer_gsl.matrix); // operazione inversa (ricavo il data set)
	stampa(1,  verifica, NULL, "Verifica", nrighe, dimensione);

	gsl_matrix_free(pointer_gsl.matrix);

	stampa(1, data_set, NULL, "Dataset", nrighe, dimensione);

	for(i=0; i<nrighe; i++){ // verifica della buona riuscita dell'operazione d'inversione
		for(j=0; j<dimensione; j++){
			//tronco alla 3° cifra decimale per trascurare piccole variazioni, dovute ad errori di approssimazione
			if((floor(data_set[i][j]*troncamento)/troncamento) != ((floor(verifica[i][j]*troncamento))/troncamento)){
				stampa(5, NULL, NULL, NULL, 0, 0);
			}
		}
	}
	printf("programma terminato.\n");

	return 0;
}


text:
#ifndef __ARNOLDI_Implicit_restart_Arnoldi_H__
#define __ARNOLDI_Implicit_restart_Arnoldi_H__

#include "Macros.h"
#include "timer.h"
#include "cuda_supp.h"
#include "memory_operations.h"
#include "Arnoldi_Driver.h"
#include "Select_Shifts.h"
#include "QR_Shifts.h"
#include "Matrix_Vector_emulator.h"
#include "file_operations.h"
//to set number of threads in openblas
#include <cblas.h>

extern "C" int openblas_get_num_threads(void);
extern "C" void openblas_set_num_threads(int);

real Implicit_restart_Arnoldi_GPU_data(cublasHandle_t handle, bool verbose, int N, user_map_vector Axb, void *user_struct, char which[2], int k, int m, complex real* eigenvaluesA, real tol, int max_iter, real *eigenvectors_real=NULL, real *eigenvectors_imag=NULL, int BLASTreads=1);

real Implicit_restart_Arnoldi_GPU_data_Matrix_Exponent(cublasHandle_t handle, bool verbose, int N,  user_map_vector Axb_exponent_invert, void *user_struct_exponent_invert, user_map_vector Axb, void *user_struct, char which[2], char which_exponent[2], int k, int m, complex real* eigenvaluesA, real tol, int max_iter, bool is_rotated, real *eigenvectors_real_d=NULL, real *eigenvectors_imag_d=NULL, int BLASThreads=1);

#endif

ext: h


text:
// Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at
// the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights
// reserved. See files LICENSE and NOTICE for details.
//
// This file is part of CEED, a collection of benchmarks, miniapps, software
// libraries and APIs for efficient high-order finite element and spectral
// element discretizations for exascale applications. For more information and
// source code availability see http://github.com/ceed.
//
// The CEED research is supported by the Exascale Computing Project 17-SC-20-SC,
// a collaborative effort of two U.S. Department of Energy organizations (Office
// of Science and the National Nuclear Security Administration) responsible for
// the planning and preparation of a capable exascale ecosystem, including
// software, applications, hardware, advanced system engineering and early
// testbed platforms, in support of the nation's exascale computing imperative.

//                        libCEED + PETSc Example: CEED BPs 3-6 with Multigrid
//
// This example demonstrates a simple usage of libCEED with PETSc to solve the
// CEED BP benchmark problems, see http://ceed.exascaleproject.org/bps.
//
// The code uses higher level communication protocols in DMPlex.
//
// Build with:
//
//     make multigrid [PETSC_DIR=</path/to/petsc>] [CEED_DIR=</path/to/libceed>]
//
// Sample runs:
//
//     multigrid -problem bp3
//     multigrid -problem bp4
//     multigrid -problem bp5 -ceed /cpu/self
//     multigrid -problem bp6 -ceed /gpu/cuda
//
//TESTARGS -ceed {ceed_resource} -test -problem bp3 -degree 3

/// @file
/// CEED BPs 1-6 multigrid example using PETSc
const char help[] = "Solve CEED BPs using p-multigrid with PETSc and DMPlex\n";

#include <stdbool.h>
#include <string.h>
#include <ceed.h>
#include <petsc.h>
#include <petscdmplex.h>
#include <petscksp.h>
#include <petscsys.h>

#include "bps.h"
#include "include/bpsproblemdata.h"
#include "include/petscmacros.h"
#include "include/petscutils.h"
#include "include/matops.h"
#include "include/structs.h"
#include "include/libceedsetup.h"

#if PETSC_VERSION_LT(3,12,0)
#ifdef PETSC_HAVE_CUDA
#include <petsccuda.h>
// Note: With PETSc prior to version 3.12.0, providing the source path to
//       include 'cublas_v2.h' will be needed to use 'petsccuda.h'.
#endif
#endif

int main(int argc, char **argv) {
  PetscInt ierr;
  MPI_Comm comm;
  char filename[PETSC_MAX_PATH_LEN],
       ceed_resource[PETSC_MAX_PATH_LEN] = "/cpu/self";
  double my_rt_start, my_rt, rt_min, rt_max;
  PetscInt degree = 3, q_extra, *l_size, *xl_size, *g_size, dim = 3, fine_level,
           mesh_elem[3] = {3, 3, 3}, num_comp_u = 1, num_levels = degree, *level_degrees;
  PetscScalar *r;
  PetscScalar eps = 1.0;
  PetscBool test_mode, benchmark_mode, read_mesh, write_solution;
  PetscLogStage solve_stage;
  DM  *dm, dm_orig;
  SNES snes_dummy;
  KSP ksp;
  PC pc;
  Mat *mat_O, *mat_pr, mat_coarse;
  Vec *X, *X_loc, *mult, rhs, rhs_loc;
  PetscMemType mem_type;
  UserO *user_O;
  UserProlongRestr *user_pr;
  Ceed ceed;
  CeedData *ceed_data;
  CeedVector rhs_ceed, target;
  CeedQFunction qf_error, qf_restrict, qf_prolong;
  CeedOperator op_error;
  BPType bp_choice;
  CoarsenType coarsen;

  ierr = PetscInitialize(&argc, &argv, NULL, help);
  if (ierr) return ierr;
  comm = PETSC_COMM_WORLD;

  // Parse command line options
  ierr = PetscOptionsBegin(comm, NULL, "CEED BPs in PETSc", NULL); CHKERRQ(ierr);
  bp_choice = CEED_BP3;
  ierr = PetscOptionsEnum("-problem",
                          "CEED benchmark problem to solve", NULL,
                          bp_types, (PetscEnum)bp_choice, (PetscEnum *)&bp_choice,
                          NULL); CHKERRQ(ierr);
  num_comp_u = bp_options[bp_choice].num_comp_u;
  test_mode = PETSC_FALSE;
  ierr = PetscOptionsBool("-test",
                          "Testing mode (do not print unless error is large)",
                          NULL, test_mode, &test_mode, NULL); CHKERRQ(ierr);
  benchmark_mode = PETSC_FALSE;
  ierr = PetscOptionsBool("-benchmark",
                          "Benchmarking mode (prints benchmark statistics)",
                          NULL, benchmark_mode, &benchmark_mode, NULL);
  CHKERRQ(ierr);
  write_solution = PETSC_FALSE;
  ierr = PetscOptionsBool("-write_solution",
                          "Write solution for visualization",
                          NULL, write_solution, &write_solution, NULL);
  CHKERRQ(ierr);
  ierr = PetscOptionsScalar("-eps",
                            "Epsilon parameter for Kershaw mesh transformation",
                            NULL, eps, &eps, NULL);
  if (eps > 1 || eps <= 0) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE,
                                      "-eps %D must be (0,1]", eps);
  degree = test_mode ? 3 : 2;
  ierr = PetscOptionsInt("-degree", "Polynomial degree of tensor product basis",
                         NULL, degree, &degree, NULL); CHKERRQ(ierr);
  if (degree < 1) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE,
                             "-degree %D must be at least 1", degree);
  q_extra = bp_options[bp_choice].q_extra;
  ierr = PetscOptionsInt("-q_extra", "Number of extra quadrature points",
                         NULL, q_extra, &q_extra, NULL); CHKERRQ(ierr);
  ierr = PetscOptionsString("-ceed", "CEED resource specifier",
                            NULL, ceed_resource, ceed_resource,
                            sizeof(ceed_resource), NULL); CHKERRQ(ierr);
  coarsen = COARSEN_UNIFORM;
  ierr = PetscOptionsEnum("-coarsen",
                          "Coarsening strategy to use", NULL,
                          coarsen_types, (PetscEnum)coarsen,
                          (PetscEnum *)&coarsen, NULL); CHKERRQ(ierr);
  read_mesh = PETSC_FALSE;
  ierr = PetscOptionsString("-mesh", "Read mesh from file", NULL,
                            filename, filename, sizeof(filename), &read_mesh);
  CHKERRQ(ierr);
  if (!read_mesh) {
    PetscInt tmp = dim;
    ierr = PetscOptionsIntArray("-cells","Number of cells per dimension", NULL,
                                mesh_elem, &tmp, NULL); CHKERRQ(ierr);
  }
  ierr = PetscOptionsEnd(); CHKERRQ(ierr);

  // Set up libCEED
  CeedInit(ceed_resource, &ceed);
  CeedMemType mem_type_backend;
  CeedGetPreferredMemType(ceed, &mem_type_backend);

  // Setup DM
  if (read_mesh) {
    ierr = DMPlexCreateFromFile(PETSC_COMM_WORLD, filename, PETSC_TRUE, &dm_orig);
    CHKERRQ(ierr);
  } else {
    ierr = DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, PETSC_FALSE, mesh_elem, NULL,
                               NULL, NULL, PETSC_TRUE, &dm_orig); CHKERRQ(ierr);
  }

  {
    DM dm_dist = NULL;
    PetscPartitioner part;

    ierr = DMPlexGetPartitioner(dm_orig, &part); CHKERRQ(ierr);
    ierr = PetscPartitionerSetFromOptions(part); CHKERRQ(ierr);
    ierr = DMPlexDistribute(dm_orig, 0, NULL, &dm_dist); CHKERRQ(ierr);
    if (dm_dist) {
      ierr = DMDestroy(&dm_orig); CHKERRQ(ierr);
      dm_orig = dm_dist;
    }
  }

  // Apply Kershaw mesh transformation
  ierr = Kershaw(dm_orig, eps); CHKERRQ(ierr);

  VecType vec_type;
  switch (mem_type_backend) {
  case CEED_MEM_HOST: vec_type = VECSTANDARD; break;
  case CEED_MEM_DEVICE: {
    const char *resolved;
    CeedGetResource(ceed, &resolved);
    if (strstr(resolved, "/gpu/cuda")) vec_type = VECCUDA;
    else if (strstr(resolved, "/gpu/hip/occa"))
      vec_type = VECSTANDARD; // https://github.com/CEED/libCEED/issues/678
    else if (strstr(resolved, "/gpu/hip")) vec_type = VECHIP;
    else vec_type = VECSTANDARD;
  }
  }
  ierr = DMSetVecType(dm_orig, vec_type); CHKERRQ(ierr);
  ierr = DMSetFromOptions(dm_orig); CHKERRQ(ierr);

  // Allocate arrays for PETSc objects for each level
  switch (coarsen) {
  case COARSEN_UNIFORM:
    num_levels = degree;
    break;
  case COARSEN_LOGARITHMIC:
    num_levels = ceil(log(degree)/log(2)) + 1;
    break;
  }
  ierr = PetscMalloc1(num_levels, &level_degrees); CHKERRQ(ierr);
  fine_level = num_levels - 1;

  switch (coarsen) {
  case COARSEN_UNIFORM:
    for (int i=0; i<num_levels; i++) level_degrees[i] = i + 1;
    break;
  case COARSEN_LOGARITHMIC:
    for (int i=0; i<num_levels - 1; i++) level_degrees[i] = pow(2,i);
    level_degrees[fine_level] = degree;
    break;
  }
  ierr = PetscMalloc1(num_levels, &dm); CHKERRQ(ierr);
  ierr = PetscMalloc1(num_levels, &X); CHKERRQ(ierr);
  ierr = PetscMalloc1(num_levels, &X_loc); CHKERRQ(ierr);
  ierr = PetscMalloc1(num_levels, &mult); CHKERRQ(ierr);
  ierr = PetscMalloc1(num_levels, &user_O); CHKERRQ(ierr);
  ierr = PetscMalloc1(num_levels, &user_pr); CHKERRQ(ierr);
  ierr = PetscMalloc1(num_levels, &mat_O); CHKERRQ(ierr);
  ierr = PetscMalloc1(num_levels, &mat_pr); CHKERRQ(ierr);
  ierr = PetscMalloc1(num_levels, &l_size); CHKERRQ(ierr);
  ierr = PetscMalloc1(num_levels, &xl_size); CHKERRQ(ierr);
  ierr = PetscMalloc1(num_levels, &g_size); CHKERRQ(ierr);

  // Setup DM and Operator Mat Shells for each level
  for (CeedInt i=0; i<num_levels; i++) {
    // Create DM
    ierr = DMClone(dm_orig, &dm[i]); CHKERRQ(ierr);
    ierr = DMGetVecType(dm_orig, &vec_type); CHKERRQ(ierr);
    ierr = DMSetVecType(dm[i], vec_type); CHKERRQ(ierr);
    PetscInt dim;
    ierr = DMGetDimension(dm[i], &dim); CHKERRQ(ierr);
    ierr = SetupDMByDegree(dm[i], level_degrees[i], num_comp_u, dim,
                           bp_options[bp_choice].enforce_bc, bp_options[bp_choice].bc_func);
    CHKERRQ(ierr);

    // Create vectors
    ierr = DMCreateGlobalVector(dm[i], &X[i]); CHKERRQ(ierr);
    ierr = VecGetLocalSize(X[i], &l_size[i]); CHKERRQ(ierr);
    ierr = VecGetSize(X[i], &g_size[i]); CHKERRQ(ierr);
    ierr = DMCreateLocalVector(dm[i], &X_loc[i]); CHKERRQ(ierr);
    ierr = VecGetSize(X_loc[i], &xl_size[i]); CHKERRQ(ierr);

    // Operator
    ierr = PetscMalloc1(1, &user_O[i]); CHKERRQ(ierr);
    ierr = MatCreateShell(comm, l_size[i], l_size[i], g_size[i], g_size[i],
                          user_O[i], &mat_O[i]); CHKERRQ(ierr);
    ierr = MatShellSetOperation(mat_O[i], MATOP_MULT,
                                (void(*)(void))MatMult_Ceed); CHKERRQ(ierr);
    ierr = MatShellSetOperation(mat_O[i], MATOP_GET_DIAGONAL,
                                (void(*)(void))MatGetDiag); CHKERRQ(ierr);
    ierr = MatShellSetVecType(mat_O[i], vec_type); CHKERRQ(ierr);

    // Level transfers
    if (i > 0) {
      // Interp
      ierr = PetscMalloc1(1, &user_pr[i]); CHKERRQ(ierr);
      ierr = MatCreateShell(comm, l_size[i], l_size[i-1], g_size[i], g_size[i-1],
                            user_pr[i], &mat_pr[i]); CHKERRQ(ierr);
      ierr = MatShellSetOperation(mat_pr[i], MATOP_MULT,
                                  (void(*)(void))MatMult_Prolong);
      CHKERRQ(ierr);
      ierr = MatShellSetOperation(mat_pr[i], MATOP_MULT_TRANSPOSE,
                                  (void(*)(void))MatMult_Restrict);
      CHKERRQ(ierr);
      ierr = MatShellSetVecType(mat_pr[i], vec_type); CHKERRQ(ierr);
    }
  }
  ierr = VecDuplicate(X[fine_level], &rhs); CHKERRQ(ierr);

  // Print global grid information
  if (!test_mode) {
    PetscInt P = degree + 1, Q = P + q_extra;

    const char *used_resource;
    CeedGetResource(ceed, &used_resource);

    ierr = VecGetType(X[0], &vec_type); CHKERRQ(ierr);

    ierr = PetscPrintf(comm,
                       "\n-- CEED Benchmark Problem %d -- libCEED + PETSc + PCMG --\n"
                       "  PETSc:\n"
                       "    PETSc Vec Type                     : %s\n"
                       "  libCEED:\n"
                       "    libCEED Backend                    : %s\n"
                       "    libCEED Backend MemType            : %s\n"
                       "  Mesh:\n"
                       "    Number of 1D Basis Nodes (p)       : %d\n"
                       "    Number of 1D Quadrature Points (q) : %d\n"
                       "    Global Nodes                       : %D\n"
                       "    Owned Nodes                        : %D\n"
                       "    DoF per node                       : %D\n"
                       "  Multigrid:\n"
                       "    Number of Levels                   : %d\n",
                       bp_choice+1, vec_type, used_resource,
                       CeedMemTypes[mem_type_backend],
                       P, Q, g_size[fine_level]/num_comp_u, l_size[fine_level]/num_comp_u,
                       num_comp_u, num_levels); CHKERRQ(ierr);
  }

  // Create RHS vector
  ierr = VecDuplicate(X_loc[fine_level], &rhs_loc); CHKERRQ(ierr);
  ierr = VecZeroEntries(rhs_loc); CHKERRQ(ierr);
  ierr = VecGetArrayAndMemType(rhs_loc, &r, &mem_type); CHKERRQ(ierr);
  CeedVectorCreate(ceed, xl_size[fine_level], &rhs_ceed);
  CeedVectorSetArray(rhs_ceed, MemTypeP2C(mem_type), CEED_USE_POINTER, r);

  // Set up libCEED operators on each level
  ierr = PetscMalloc1(num_levels, &ceed_data); CHKERRQ(ierr);
  for (int i=0; i<num_levels; i++) {
    // Print level information
    if (!test_mode && (i == 0 || i == fine_level)) {
      ierr = PetscPrintf(comm,"    Level %D (%s):\n"
                         "      Number of 1D Basis Nodes (p)     : %d\n"
                         "      Global Nodes                     : %D\n"
                         "      Owned Nodes                      : %D\n",
                         i, (i? "fine" : "coarse"), level_degrees[i] + 1,
                         g_size[i]/num_comp_u, l_size[i]/num_comp_u); CHKERRQ(ierr);
    }
    ierr = PetscMalloc1(1, &ceed_data[i]); CHKERRQ(ierr);
    ierr = SetupLibceedByDegree(dm[i], ceed, level_degrees[i], dim, q_extra,
                                dim, num_comp_u, g_size[i], xl_size[i], bp_options[bp_choice],
                                ceed_data[i], i==(fine_level), rhs_ceed, &target);
    CHKERRQ(ierr);
  }

  // Gather RHS
  CeedVectorTakeArray(rhs_ceed, MemTypeP2C(mem_type), NULL);
  ierr = VecRestoreArrayAndMemType(rhs_loc, &r); CHKERRQ(ierr);
  ierr = VecZeroEntries(rhs); CHKERRQ(ierr);
  ierr = DMLocalToGlobal(dm[fine_level], rhs_loc, ADD_VALUES, rhs); CHKERRQ(ierr);
  CeedVectorDestroy(&rhs_ceed);

  // Create the restriction/interpolation QFunction
  CeedQFunctionCreateIdentity(ceed, num_comp_u, CEED_EVAL_NONE, CEED_EVAL_INTERP,
                              &qf_restrict);
  CeedQFunctionCreateIdentity(ceed, num_comp_u, CEED_EVAL_INTERP, CEED_EVAL_NONE,
                              &qf_prolong);

  // Set up libCEED level transfer operators
  ierr = CeedLevelTransferSetup(ceed, num_levels, num_comp_u, ceed_data,
                                level_degrees,
                                qf_restrict, qf_prolong); CHKERRQ(ierr);

  // Create the error QFunction
  CeedQFunctionCreateInterior(ceed, 1, bp_options[bp_choice].error,
                              bp_options[bp_choice].error_loc, &qf_error);
  CeedQFunctionAddInput(qf_error, "u", num_comp_u, CEED_EVAL_INTERP);
  CeedQFunctionAddInput(qf_error, "true_soln", num_comp_u, CEED_EVAL_NONE);
  CeedQFunctionAddOutput(qf_error, "error", num_comp_u, CEED_EVAL_NONE);

  // Create the error operator
  CeedOperatorCreate(ceed, qf_error, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE,
                     &op_error);
  CeedOperatorSetField(op_error, "u", ceed_data[fine_level]->elem_restr_u,
                       ceed_data[fine_level]->basis_u, CEED_VECTOR_ACTIVE);
  CeedOperatorSetField(op_error, "true_soln",
                       ceed_data[fine_level]->elem_restr_u_i,
                       CEED_BASIS_COLLOCATED, target);
  CeedOperatorSetField(op_error, "error", ceed_data[fine_level]->elem_restr_u_i,
                       CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE);

  // Calculate multiplicity
  for (int i=0; i<num_levels; i++) {
    PetscScalar *x;

    // CEED vector
    ierr = VecZeroEntries(X_loc[i]); CHKERRQ(ierr);
    ierr = VecGetArray(X_loc[i], &x); CHKERRQ(ierr);
    CeedVectorSetArray(ceed_data[i]->x_ceed, CEED_MEM_HOST, CEED_USE_POINTER, x);

    // Multiplicity
    CeedElemRestrictionGetMultiplicity(ceed_data[i]->elem_restr_u,
                                       ceed_data[i]->x_ceed);
    CeedVectorSyncArray(ceed_data[i]->x_ceed, CEED_MEM_HOST);

    // Restore vector
    ierr = VecRestoreArray(X_loc[i], &x); CHKERRQ(ierr);

    // Creat mult vector
    ierr = VecDuplicate(X_loc[i], &mult[i]); CHKERRQ(ierr);

    // Local-to-global
    ierr = VecZeroEntries(X[i]); CHKERRQ(ierr);
    ierr = DMLocalToGlobal(dm[i], X_loc[i], ADD_VALUES, X[i]);
    CHKERRQ(ierr);
    ierr = VecZeroEntries(X_loc[i]); CHKERRQ(ierr);

    // Global-to-local
    ierr = DMGlobalToLocal(dm[i], X[i], INSERT_VALUES, mult[i]);
    CHKERRQ(ierr);
    ierr = VecZeroEntries(X[i]); CHKERRQ(ierr);

    // Multiplicity scaling
    ierr = VecReciprocal(mult[i]);
  }

  // Set up Mat
  for (int i=0; i<num_levels; i++) {
    // User Operator
    user_O[i]->comm = comm;
    user_O[i]->dm = dm[i];
    user_O[i]->X_loc = X_loc[i];
    ierr = VecDuplicate(X_loc[i], &user_O[i]->Y_loc); CHKERRQ(ierr);
    user_O[i]->x_ceed = ceed_data[i]->x_ceed;
    user_O[i]->y_ceed = ceed_data[i]->y_ceed;
    user_O[i]->op = ceed_data[i]->op_apply;
    user_O[i]->ceed = ceed;

    if (i > 0) {
      // Prolongation/Restriction Operator
      user_pr[i]->comm = comm;
      user_pr[i]->dmf = dm[i];
      user_pr[i]->dmc = dm[i-1];
      user_pr[i]->loc_vec_c = X_loc[i-1];
      user_pr[i]->loc_vec_f = user_O[i]->Y_loc;
      user_pr[i]->mult_vec = mult[i];
      user_pr[i]->ceed_vec_c = user_O[i-1]->x_ceed;
      user_pr[i]->ceed_vec_f = user_O[i]->y_ceed;
      user_pr[i]->op_prolong = ceed_data[i]->op_prolong;
      user_pr[i]->op_restrict = ceed_data[i]->op_restrict;
      user_pr[i]->ceed = ceed;
    }
  }

  // Setup dummy SNES for AMG coarse solve
  ierr = SNESCreate(comm, &snes_dummy); CHKERRQ(ierr);
  ierr = SNESSetDM(snes_dummy, dm[0]); CHKERRQ(ierr);
  ierr = SNESSetSolution(snes_dummy, X[0]); CHKERRQ(ierr);

  // -- Jacobian matrix
  ierr = DMSetMatType(dm[0], MATAIJ); CHKERRQ(ierr);
  ierr = DMCreateMatrix(dm[0], &mat_coarse); CHKERRQ(ierr);
  ierr = SNESSetJacobian(snes_dummy, mat_coarse, mat_coarse, NULL,
                         NULL); CHKERRQ(ierr);

  // -- Residual evaluation function
  ierr = SNESSetFunction(snes_dummy, X[0], FormResidual_Ceed,
                         user_O[0]); CHKERRQ(ierr);

  // -- Form Jacobian
  ierr = SNESComputeJacobianDefaultColor(snes_dummy, X[0], mat_O[0],
                                         mat_coarse, NULL); CHKERRQ(ierr);

  // Set up KSP
  ierr = KSPCreate(comm, &ksp); CHKERRQ(ierr);
  {
    ierr = KSPSetType(ksp, KSPCG); CHKERRQ(ierr);
    ierr = KSPSetNormType(ksp, KSP_NORM_NATURAL); CHKERRQ(ierr);
    ierr = KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT,
                            PETSC_DEFAULT); CHKERRQ(ierr);
  }
  ierr = KSPSetFromOptions(ksp); CHKERRQ(ierr);
  ierr = KSPSetOperators(ksp, mat_O[fine_level], mat_O[fine_level]);
  CHKERRQ(ierr);

  // Set up PCMG
  ierr = KSPGetPC(ksp, &pc); CHKERRQ(ierr);
  PCMGCycleType pcmg_cycle_type = PC_MG_CYCLE_V;
  {
    ierr = PCSetType(pc, PCMG); CHKERRQ(ierr);

    // PCMG levels
    ierr = PCMGSetLevels(pc, num_levels, NULL); CHKERRQ(ierr);
    for (int i=0; i<num_levels; i++) {
      // Smoother
      KSP smoother;
      PC smoother_pc;
      ierr = PCMGGetSmoother(pc, i, &smoother); CHKERRQ(ierr);
      ierr = KSPSetType(smoother, KSPCHEBYSHEV); CHKERRQ(ierr);
      ierr = KSPChebyshevEstEigSet(smoother, 0, 0.1, 0, 1.1); CHKERRQ(ierr);
      ierr = KSPChebyshevEstEigSetUseNoisy(smoother, PETSC_TRUE); CHKERRQ(ierr);
      ierr = KSPSetOperators(smoother, mat_O[i], mat_O[i]); CHKERRQ(ierr);
      ierr = KSPGetPC(smoother, &smoother_pc); CHKERRQ(ierr);
      ierr = PCSetType(smoother_pc, PCJACOBI); CHKERRQ(ierr);
      ierr = PCJacobiSetType(smoother_pc, PC_JACOBI_DIAGONAL); CHKERRQ(ierr);

      // Work vector
      if (i < num_levels - 1) {
        ierr = PCMGSetX(pc, i, X[i]); CHKERRQ(ierr);
      }

      // Level transfers
      if (i > 0) {
        // Interpolation
        ierr = PCMGSetInterpolation(pc, i, mat_pr[i]); CHKERRQ(ierr);
      }

      // Coarse solve
      KSP coarse;
      PC coarse_pc;
      ierr = PCMGGetCoarseSolve(pc, &coarse); CHKERRQ(ierr);
      ierr = KSPSetType(coarse, KSPPREONLY); CHKERRQ(ierr);
      ierr = KSPSetOperators(coarse, mat_coarse, mat_coarse); CHKERRQ(ierr);

      ierr = KSPGetPC(coarse, &coarse_pc); CHKERRQ(ierr);
      ierr = PCSetType(coarse_pc, PCGAMG); CHKERRQ(ierr);

      ierr = KSPSetOptionsPrefix(coarse, "coarse_"); CHKERRQ(ierr);
      ierr = PCSetOptionsPrefix(coarse_pc, "coarse_"); CHKERRQ(ierr);
      ierr = KSPSetFromOptions(coarse); CHKERRQ(ierr);
      ierr = PCSetFromOptions(coarse_pc); CHKERRQ(ierr);
    }

    // PCMG options
    ierr = PCMGSetType(pc, PC_MG_MULTIPLICATIVE); CHKERRQ(ierr);
    ierr = PCMGSetNumberSmooth(pc, 3); CHKERRQ(ierr);
    ierr = PCMGSetCycleType(pc, pcmg_cycle_type); CHKERRQ(ierr);
  }

  // First run, if benchmarking
  if (benchmark_mode) {
    ierr = KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, 1);
    CHKERRQ(ierr);
    ierr = VecZeroEntries(X[fine_level]); CHKERRQ(ierr);
    my_rt_start = MPI_Wtime();
    ierr = KSPSolve(ksp, rhs, X[fine_level]); CHKERRQ(ierr);
    my_rt = MPI_Wtime() - my_rt_start;
    ierr = MPI_Allreduce(MPI_IN_PLACE, &my_rt, 1, MPI_DOUBLE, MPI_MIN, comm);
    CHKERRQ(ierr);
    // Set maxits based on first iteration timing
    if (my_rt > 0.02) {
      ierr = KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, 5);
      CHKERRQ(ierr);
    } else {
      ierr = KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, 20);
      CHKERRQ(ierr);
    }
  }

  // Timed solve
  ierr = VecZeroEntries(X[fine_level]); CHKERRQ(ierr);
  ierr = PetscBarrier((PetscObject)ksp); CHKERRQ(ierr);

  // -- Performance logging
  ierr = PetscLogStageRegister("Solve Stage", &solve_stage); CHKERRQ(ierr);
  ierr = PetscLogStagePush(solve_stage); CHKERRQ(ierr);

  // -- Solve
  my_rt_start = MPI_Wtime();
  ierr = KSPSolve(ksp, rhs, X[fine_level]); CHKERRQ(ierr);
  my_rt = MPI_Wtime() - my_rt_start;

  // -- Performance logging
  ierr = PetscLogStagePop();

  // Output results
  {
    KSPType ksp_type;
    PCMGType pcmg_type;
    KSPConvergedReason reason;
    PetscReal rnorm;
    PetscInt its;
    ierr = KSPGetType(ksp, &ksp_type); CHKERRQ(ierr);
    ierr = KSPGetConvergedReason(ksp, &reason); CHKERRQ(ierr);
    ierr = KSPGetIterationNumber(ksp, &its); CHKERRQ(ierr);
    ierr = KSPGetResidualNorm(ksp, &rnorm); CHKERRQ(ierr);
    ierr = PCMGGetType(pc, &pcmg_type); CHKERRQ(ierr);
    if (!test_mode || reason < 0 || rnorm > 1e-8) {
      ierr = PetscPrintf(comm,
                         "  KSP:\n"
                         "    KSP Type                           : %s\n"
                         "    KSP Convergence                    : %s\n"
                         "    Total KSP Iterations               : %D\n"
                         "    Final rnorm                        : %e\n",
                         ksp_type, KSPConvergedReasons[reason], its,
                         (double)rnorm); CHKERRQ(ierr);
      ierr = PetscPrintf(comm,
                         "  PCMG:\n"
                         "    PCMG Type                          : %s\n"
                         "    PCMG Cycle Type                    : %s\n",
                         PCMGTypes[pcmg_type],
                         PCMGCycleTypes[pcmg_cycle_type]); CHKERRQ(ierr);
    }
    if (!test_mode) {
      ierr = PetscPrintf(comm,"  Performance:\n"); CHKERRQ(ierr);
    }
    {
      PetscReal max_error;
      ierr = ComputeErrorMax(user_O[fine_level], op_error, X[fine_level], target,
                             &max_error); CHKERRQ(ierr);
      PetscReal tol = 5e-2;
      if (!test_mode || max_error > tol) {
        ierr = MPI_Allreduce(&my_rt, &rt_min, 1, MPI_DOUBLE, MPI_MIN, comm);
        CHKERRQ(ierr);
        ierr = MPI_Allreduce(&my_rt, &rt_max, 1, MPI_DOUBLE, MPI_MAX, comm);
        CHKERRQ(ierr);
        ierr = PetscPrintf(comm,
                           "    Pointwise Error (max)              : %e\n"
                           "    CG Solve Time                      : %g (%g) sec\n",
                           (double)max_error, rt_max, rt_min); CHKERRQ(ierr);
      }
    }
    if (benchmark_mode && (!test_mode)) {
      ierr = PetscPrintf(comm,
                         "    DoFs/Sec in CG                     : %g (%g) million\n",
                         1e-6*g_size[fine_level]*its/rt_max,
                         1e-6*g_size[fine_level]*its/rt_min);
      CHKERRQ(ierr);
    }
  }

  if (write_solution) {
    PetscViewer vtk_viewer_soln;

    ierr = PetscViewerCreate(comm, &vtk_viewer_soln); CHKERRQ(ierr);
    ierr = PetscViewerSetType(vtk_viewer_soln, PETSCVIEWERVTK); CHKERRQ(ierr);
    ierr = PetscViewerFileSetName(vtk_viewer_soln, "solution.vtu"); CHKERRQ(ierr);
    ierr = VecView(X[fine_level], vtk_viewer_soln); CHKERRQ(ierr);
    ierr = PetscViewerDestroy(&vtk_viewer_soln); CHKERRQ(ierr);
  }

  // Cleanup
  for (int i=0; i<num_levels; i++) {
    ierr = VecDestroy(&X[i]); CHKERRQ(ierr);
    ierr = VecDestroy(&X_loc[i]); CHKERRQ(ierr);
    ierr = VecDestroy(&mult[i]); CHKERRQ(ierr);
    ierr = VecDestroy(&user_O[i]->Y_loc); CHKERRQ(ierr);
    ierr = MatDestroy(&mat_O[i]); CHKERRQ(ierr);
    ierr = PetscFree(user_O[i]); CHKERRQ(ierr);
    if (i > 0) {
      ierr = MatDestroy(&mat_pr[i]); CHKERRQ(ierr);
      ierr = PetscFree(user_pr[i]); CHKERRQ(ierr);
    }
    ierr = CeedDataDestroy(i, ceed_data[i]); CHKERRQ(ierr);
    ierr = DMDestroy(&dm[i]); CHKERRQ(ierr);
  }
  ierr = PetscFree(level_degrees); CHKERRQ(ierr);
  ierr = PetscFree(dm); CHKERRQ(ierr);
  ierr = PetscFree(X); CHKERRQ(ierr);
  ierr = PetscFree(X_loc); CHKERRQ(ierr);
  ierr = PetscFree(mult); CHKERRQ(ierr);
  ierr = PetscFree(mat_O); CHKERRQ(ierr);
  ierr = PetscFree(mat_pr); CHKERRQ(ierr);
  ierr = PetscFree(ceed_data); CHKERRQ(ierr);
  ierr = PetscFree(user_O); CHKERRQ(ierr);
  ierr = PetscFree(user_pr); CHKERRQ(ierr);
  ierr = PetscFree(l_size); CHKERRQ(ierr);
  ierr = PetscFree(xl_size); CHKERRQ(ierr);
  ierr = PetscFree(g_size); CHKERRQ(ierr);
  ierr = VecDestroy(&rhs); CHKERRQ(ierr);
  ierr = VecDestroy(&rhs_loc); CHKERRQ(ierr);
  ierr = MatDestroy(&mat_coarse); CHKERRQ(ierr);
  ierr = KSPDestroy(&ksp); CHKERRQ(ierr);
  ierr = SNESDestroy(&snes_dummy); CHKERRQ(ierr);
  ierr = DMDestroy(&dm_orig); CHKERRQ(ierr);
  CeedVectorDestroy(&target);
  CeedQFunctionDestroy(&qf_error);
  CeedQFunctionDestroy(&qf_restrict);
  CeedQFunctionDestroy(&qf_prolong);
  CeedOperatorDestroy(&op_error);
  CeedDestroy(&ceed);
  return PetscFinalize();
}

max_stars_repo_name: AdelekeBankole/libCEED

max_issues_repo_name: AdelekeBankole/libCEED

max_forks_repo_name: AdelekeBankole/libCEED


text:
/*
    learn.c
*/
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_multifit.h>
#include "learn.h"
#include "feature.h"
#include "imatrix.h"
#include "dmatrix.h"
#include "util.h"
#include "likelihood.h"
#include "hyper.h"

void mvslda_learn(document *data, double **resp, double *alpha, double beta, double nu2, double sigma2,
    int nclass, int nlex, int dlenmax, int nresp, int maxiter, double **phi, double **theta, double **eta,
    int **n_mz, int **n_zw, FILE *likp, FILE *hyperp, unsigned long int random_seed){
    document *dp;
    int ndocs;
    int *n_m;
    int *n_z;
    int ***topics;
    int word_index;
    int word_num;
    double sum_alpha;
    double *left;
    double *center;
    double *log_right;
    double *p_z;
    double *log_p_z;
    double *cum_sum_p_z;
    double log_Z, sum_p_r, sum_empirical_z;
    double temp_prediction;
    double lik;
    double **temp_phi;
    double **temp_theta;
    double **empirical_z;
    int z;
    int it;
    int m, w, t, i, j, k;
    const gsl_rng_type *T;
    gsl_rng *r;

    gsl_rng_env_setup();
    T = gsl_rng_default;
    r = gsl_rng_alloc(T);

    // count data length
    for(dp = data, ndocs = 0;(dp->len) != -1;dp++, ndocs++)
        ;

    // initialize buffers
    if((n_m = calloc(ndocs,sizeof(int))) == NULL){
        fprintf(stderr,"mvslda_learn:: cannot allocate n_m.\n");
        return;
    }
    if((n_z = calloc(nclass,sizeof(int))) == NULL){
        fprintf(stderr,"mvslda_learn:: cannot allocate n_z.\n");
        return;
    }
    if((left = calloc(nclass,sizeof(double))) == NULL){
        fprintf(stderr,"mvslda_learn:: cannot allocate left.\n");
        return;
    }
    if((center = calloc(nclass,sizeof(double))) == NULL){
        fprintf(stderr,"mvslda_learn:: cannot allocate center.\n");
        return;
    }
    if((log_right = calloc(nclass,sizeof(double))) == NULL){
        fprintf(stderr,"mvslda_learn:: cannot allocate log_cright.\n");
        return;
    }
    if((p_z = calloc(nclass,sizeof(double))) == NULL){
        fprintf(stderr,"mvslda_learn:: cannot allocate p_z.\n");
        return;
    }
    if((log_p_z = calloc(nclass,sizeof(double))) == NULL){
        fprintf(stderr,"mvslda_learn:: cannot allocate log_p_z.\n");
        return;
    }
    if((cum_sum_p_z = calloc((nclass+1),sizeof(double))) == NULL){
        fprintf(stderr,"mvslda_learn:: cannot allocate cum_sum_p_z.\n");
        return;
    }
    if((topics = calloc(ndocs,sizeof(int **))) == NULL){
        fprintf(stderr,"mvslda_learn:: cannot allocate topics.\n");
        return;
    }
    if((empirical_z = dmatrix(ndocs, nclass)) == NULL){
        fprintf(stderr,"mvslda_learn:: cannot allocate empirical_z.\n");
        return;
    }
    if((temp_phi = dmatrix(nlex, nclass)) == NULL){
        fprintf(stderr,"mvslda_learn:: cannot allocate temp_phi.\n");
        exit(1);
    }
    if((temp_theta = dmatrix(ndocs, nclass)) == NULL){
        fprintf(stderr,"mvslda_learn:: cannot allocate temp_theta.\n");
        exit(1);
    }

    printf("Number of documents          = %d\n",ndocs);
    printf("Number of unique words       = %d\n",nlex);
    printf("Number of latent classes     = %d\n",nclass);
    printf("Number of responses          = %d\n",nresp);
    printf("Number of iteration          = %d\n",maxiter);

    // choose an arbitrary topic as first topic for word
    gsl_rng_set(r, random_seed);
    for(dp = data, m = 0;(dp->len) != -1;dp++, m++){
        if((topics[m] = calloc((dp->len), sizeof(int *))) == NULL){
            fprintf(stderr,"mvslda_learn:: cannot allocate topics[m].\n");
            return;
        }
        for(w = 0;w < (dp->len);w++){
            if((topics[m][w] = calloc((dp->cnt[w]), sizeof(int))) == NULL){
                fprintf(stderr,"mvslda_learn:: cannot allocate topics[m][w].\n");
                return;
            }
            word_index = dp->id[w];
            word_num = dp->cnt[w];
            for(i = 0;i < word_num;i++){
                z = (int)gsl_rng_uniform_int(r, nclass);
                n_mz[m][z] += 1;
                n_m[m] += 1;
                n_zw[z][word_index] += 1;
                n_z[z] += 1;
                topics[m][w][i] = z;
            }
        }
    }
    // initialize eta ([nclass, nresp])
    for(k = 0;k < nclass;k++)
        for(t = 0;t < nresp;t++)
            eta[k][t] = 0.0;

    // learning main
    for(it = 0;it < maxiter;it++){
        printf("iteration %2d/%3d..\n", it + 1, maxiter);
        fflush(stdout);
        sum_alpha = 0.0;
        for(k = 0;k < nclass;k++)
            sum_alpha += alpha[k];
        for (dp = data, m = 0; (dp->len) != -1; dp++, m++){
            // for words
            for(w = 0;w < (dp->len);w++){
                word_index = dp->id[w];
                word_num = dp->cnt[w];
                for(i = 0;i < word_num;i++){
                    z = topics[m][w][i];
                    n_mz[m][z] -= 1;
                    n_m[m] -= 1;
                    n_zw[z][word_index] -= 1;
                    n_z[z] -= 1;

                    // compute conditional distribution log_p_z
                    // log_p_z left ... theta term
                    for(k = 0;k < nclass;k++){
                        left[k] = (double)n_mz[m][k] + alpha[k];
                        left[k] /= ((double)n_m[m] + sum_alpha);
                    }
                    // log_p_z center ... phi term
                    for(k = 0;k < nclass;k++){
                        center[k] = (double)n_zw[k][word_index] + beta;
                        center[k] /= ((double)n_z[k] + (double)nlex * beta);
                    }
                    // temporal log_p_z (left and center)
                    for(k = 0; k < nclass;k++){
                        log_p_z[k] = log(left[k]) + log(center[k]);
                    }
                    // p_z right ... eta term
                    sum_empirical_z = 0.0;
                    for(k = 0; k < nclass;k++){
                        empirical_z[m][k] = (double)n_mz[m][k];
                        sum_empirical_z += (double)n_mz[m][k];
                    }
                    for(k = 0; k < nclass;k++){
                        empirical_z[m][k] = empirical_z[m][k] / sum_empirical_z;
                    }
                    for(t = 0; t < nresp; t++){
                        temp_prediction = 0.0;
                        for(k = 0;k < nclass;k++){
                            temp_prediction += eta[k][t] * empirical_z[m][k]; // dot(eta, z_d)
                        }
                        for(k = 0;k < nclass;k++){
                            log_right[k] = 1.0;
                            log_right[k] *= 1.0 / (2 * sigma2);
                            log_right[k] *= (eta[k][t] / (double)n_m[m]);
                            log_right[k] *= (2 * (resp[m][t] - temp_prediction) - (eta[k][t] / (double)n_m[m]));
                            log_p_z[k] += log_right[k];
                        }
                    }
                    // conditional distribution log_p_z
                    // log_Z = logsumexp(logP_k1 + logP_k2 + ... logP_kK)
                    log_Z = logsumexp(log_p_z, nclass);
                    for(k = 0;k < nclass;k++){
                        p_z[k] = exp(log_p_z[k] - log_Z); // normalize to obtain probabilities
                    }
                    // random sampling from p_z
                    z = sampling_multinomial(r, p_z, cum_sum_p_z, nclass);
                    // update buffers
                    n_mz[m][z] += 1;
                    n_m[m] += 1;
                    n_zw[z][word_index] += 1;
                    n_z[z] += 1;
                    topics[m][w][i] = z;
                }
            }
        }

        // for eta update
        //least squares for dot(Z, eta) = resp
        for(m = 0; m < ndocs; m++){
            sum_empirical_z = 0.0;
            for(k = 0; k < nclass;k++){
                empirical_z[m][k] = (double)n_mz[m][k];
                sum_empirical_z += (double)n_mz[m][k];
            }
            for(k = 0; k < nclass;k++){
                empirical_z[m][k] = empirical_z[m][k] / sum_empirical_z;
            }
        }
        for(t = 0; t < nresp; t++){
            double chisq;
            gsl_matrix *Z, *cov;
            gsl_vector *y, *c;
            gsl_multifit_linear_workspace * work = gsl_multifit_linear_alloc (ndocs, nclass);
            Z = gsl_matrix_alloc(ndocs, nclass);
            cov = gsl_matrix_alloc(nclass, nclass);
            y = gsl_vector_alloc(ndocs);
            c = gsl_vector_alloc(nclass);
            for(m = 0;m < ndocs;m++){
                gsl_vector_set(y, m, resp[m][t]);
                for(k = 0;k < nclass;k++){
                    gsl_matrix_set(Z, m, k, empirical_z[m][k]);
                }
            }
            gsl_multifit_linear(Z, y, c, cov, &chisq, work);
            for(k = 0;k < nclass;k++){
                eta[k][t] = gsl_vector_get(c, k);
            }
            gsl_multifit_linear_free(work);
            gsl_matrix_free(Z);
            gsl_matrix_free(cov);
            gsl_vector_free(y);
            gsl_vector_free(c);
        }

        // update hyperparameters.
        update_alpha(alpha, n_m, n_mz, ndocs, nclass);
        beta = update_beta(beta, n_z, n_zw, nclass, nlex);

        // compute likelihood.
        lik = loglikelihood(n_mz, n_zw, n_m, nclass, nlex, ndocs, nresp, resp, alpha, beta, eta, empirical_z, nu2, sigma2);
        printf("\tlikelihood ... %.8f\n",lik);
        printf("\talpha = \n\t");
        for(k = 0;k < nclass;k++)
            printf("%.8f ",alpha[k]);
        printf("\n\tbeta ... %.2f\n",beta);
        fprintf(likp,"%.8f\n",lik);
        for(k = 0;k < nclass;k++)
            fprintf(hyperp,"%.8f,",alpha[k]);
        fprintf(hyperp,"%.8f\n",beta);
    }

    // compute matrix phi ([nlex, nclass] matrix)
    for(w = 0;w < nlex;w++)
        for(k = 0;k < nclass;k++)
            temp_phi[w][k] = (double)n_zw[k][w] + beta;
    normalize_matrix_col(phi, temp_phi, nlex, nclass);

    // compute matrix theta ([ndocs, nclass])
    for(m = 0;m < ndocs;m++)
        for(k = 0;k < nclass;k++)
            temp_theta[m][k] = (double)n_mz[m][k] + alpha[k];
    normalize_matrix_row(theta, temp_theta, ndocs, nclass);

    free(n_m);
    free(n_z);
    free(left);
    free(center);
    free(log_right);
    free(p_z);
    free(log_p_z);
    free(cum_sum_p_z);

    for(dp = data, m = 0;(dp->len) != -1;dp++, m++){
        for(w = 0;w < (dp->len);w++){
            free(topics[m][w]);
        }
        free(topics[m]);
    }
    free(topics);
    free_dmatrix(temp_phi, nlex);
    free_dmatrix(temp_theta, ndocs);
    free_dmatrix(empirical_z, ndocs);

    return;
}

int sampling_multinomial(gsl_rng *r, double *p, double *cum_sum_p, int len_p){
    int k, z;
    double sampling;

    cum_sum_p[0] = 0.0;
    for(k = 0;k < len_p;k++){
        cum_sum_p[k+1] = cum_sum_p[k] + p[k];
    }
    sampling = gsl_rng_uniform(r);
    for(k = 0;k < len_p;k++){
        if((sampling >= cum_sum_p[k]) && (sampling < cum_sum_p[k+1])){
            z = k;
            break;
        }
    }
    return z;
}


text:
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_rng.h>
#include <math.h>
#include <time.h>
#include <stdbool.h>
#include "bb_smc.h"

/* ==========================================================
Last update: 1.5.20 -- v2.0

ABC-SMC (parallel version)
Author: Marco Esposito
=========================================================== */

// =========== Define some constants =============
#define size_data 50

#define a0 0.5
#define b0 0.5
#define std_ker 0.3
#define striding_size 1
#define n 20
// =====================================

// ARGUMENT NUMBER 1: time at which the simulation starts
// ARGUMENT NUMBER 2: seed
// ARGUMENT NUMBER 3: number of particles to consider
// for the simulation (lesser for shorter runs)

int main(int argc, char *argv[]) {

	if (argc != 4) {

		printf("argc = %d\n", argc);
		printf("Usage: %s\n", argv[0]);
		printf("      %s [parameters] \n", argv[0]);
		exit(1);

	}

	/* ==================== VARIABLES DECLARATION ======================== */
	int time = atoi(argv[1]);
	int seed = atoi(argv[2]);
	int n_particles = atoi(argv[3]);

	double theta_star;
	double theta;
	double weight;
	double epsilon_t;
	int status;

	double* weights;
	weights = malloc((n_particles) * sizeof(double));

	int sample_accepted;
	/* ==================================================================== */

	/* ============== IMPORT DATA =========== */
	int y[size_data];
	status = import_data_BB(y, size_data);
	if (status == 1) {
		printf("Data not found\n");
		exit(1);
	}
	double mean_y = mean_int(y, size_data);
	// printf("(MAIN) -- mean_y: %f\n", mean_y);
	/* ======================================== */

	/* =========== READ THRESHOLD ============ */
	status = read_threshold(&epsilon_t, time);
	if (status == 1) {
		printf("Threshold not found\n");
		exit(1);
	}
	/* ======================================== */

	/* ================== set up GSL RNG =================== */
	gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
	gsl_rng_set (r, seed+101); // set seed
	/* ===================================================== */

	sample_accepted = 1;

	/* =========================================================== */
	/* ================= START THE ABC ROUND ===================== */
	/* =========================================================== */
	while (sample_accepted == 1) {

		// Select particle
		status = select_particle(weights, &theta, time, n_particles, r);
		if (status == 1) {
			printf("File with particles at time %d, %d particles not present\n", time, n_particles);
			exit(1);
		}

		// Perturb the parameter
		status = ker_perturb(&theta_star, r, theta, a0, b0, std_ker);
		if (status == 1) {
			continue;
		}

		// - If the distance function is too large, start over
		// - If all the conditions are fulfilled, the sample is accepted
		sample_accepted = verify_par_BB(theta_star, mean_y, r, size_data, n, epsilon_t, time, seed);
		if (sample_accepted == 1) {
			continue;
		}

	} //END OF ABC ROUND (WHILE LOOP) -- PARTICLE WAS ACCEPTED

	calculate_weight_BB(&weight, weights, theta, theta_star, n_particles, a0, b0, std_ker);

	write_particles_BB(theta_star, time, seed, striding_size, weight);

	free(weights);

	return 0; //return of main

} //End of the program


text:
/*
 * This file is part of the Visual Computing Library (VCL) release under the
 * license.
 *
 * Copyright (c) 2017 Basil Fierz
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */
#pragma once

// VCL configuration
#include <vcl/config/global.h>

// C++ Standard library
#include <memory>
#include <vector>

// GSL
#include <gsl/gsl>

// VCL
#include <vcl/hid/windows/hid.h>
#include <vcl/hid/spacenavigator.h>

namespace Vcl { namespace HID { namespace Windows
{
	//! Implementation of a 3Dconnexion space mouse
	class SpaceNavigatorHID : public AbstractHID, public SpaceNavigator
	{
	private:
		struct InputData
		{
			//! Current time to live for telling
			//! if the device was unplugged while sending data
			int timeToLive;

			//! Indicate if the data is dirty
			bool isDirty;

			//! Axis data
			std::array<float,6> axes;

			//! Check if the data is zero
			bool isZero()
			{
				return (0 == axes[0] && 0 == axes[1] && 0 == axes[2] &&
						0 == axes[3] && 0 == axes[4] && 0 == axes[5]    );
			}

			//! Maximum time to live
			static const int MaxTimeToLive = 5;
		};

	public:
		SpaceNavigatorHID(std::unique_ptr<GenericHID> dev, bool poll_3d_mouse = false);

		//! Reset device when activating the program
		void onActivateApp(BOOL active, DWORD dwThreadID);

		//! Handle device input
		bool processInput(HWND window_handle, UINT input_code, PRAWINPUT raw_input) override;

	private:

		/*!
		 * \brief Does all the preprocessing of the rawinput device data before
		 *        finally calling the move3D method.
		 *
		 * If polling is enabled (_poll3DMouse == true) this method is called
		 * from the windows timer message handler
		 * If polling is not enabled (_poll3DMouse == false) this method is
		 * called directly from the WM_INPUT handler
		 */
		void on3DMouseInput();

		/*!
		 * \brief onSpaceMouseMove is invoked when new 3d mouse data is
		 *        available.
		 *
		 * \param motion_data Contains the displacement data, using a
		 *                    right-handed coordinate system with z down.
		 *                    See 'Programing for the 3dmouse' document
		 *                    available at www.3dconnexion.com.
		 *                    Entries 0, 1, 2 is the incremental pan zoom
		 *                    displacement vector (x,y,z).
		 *                    Entries 3, 4, 5 is the incremental rotation vector
		 *                    (NOT Euler angles).
		 */
		void onSpaceMouseMove(std::array<float, 6> motion_data);

		/*!
		 * \brief onSpaceMouseKeyDown processes the 3d mouse key presses
		 *
		 * \param virtual_key 3d mouse key code
		 */
		void onSpaceMouseKeyDown(UINT virtual_key);

		/*!
		 * \brief onSpaceMouseKeyUp processes the 3d mouse key releases
		 *
		 * \param virtual_key 3d mouse key code
		 */
		void onSpaceMouseKeyUp(UINT virtual_key);

	private:
		//! Process a raw input message
		bool translateRawInputData(UINT input_code, PRAWINPUT raw_input);

		//! Axis input data
		InputData _deviceData;

		//! Button input data
		uint32_t _keystate{ 0 };

		//! Last time the data was updated.
		//! Use to calculate distance traveled since last event
		DWORD _last3DMouseInputTime{ 0 };

	private: // Polling support methods

		//! Start the timer
		void startTimer(HWND hwnd);

		//! Timer callback
		void onTimer(UINT_PTR event_id);

		//! Kill the currently running timer
		void killPollingTimer();

	private: // Polling support
		//! 3D mouse is in polling mode
		bool _poll3DMouse{ false };

		//! Polling period. Default is 50 Hz)
		UINT _pollingPeriod3DMouse{ 20 };

		//! 3DMouse data polling timer
		//! Only used if _poll3DMouse == true
		UINT_PTR _timer3DMouse{ 0 };
	};
}}}


text:
#include <stdio.h>
#include <gsl/gsl_statistics.h>
#include <stdlib.h>
#include <float.h>
#include <assert.h>

#ifndef N
#define N 64
#endif

#define K N/2

#ifndef S
#define S K/2
#endif

#ifndef SEED
#define SEED 1
#endif

#include <klee/klee.h>
int
main(int argc, char* argv[])
{
  double data[K], w[K];
  int i, j;

  srand(SEED);

  int index[K];
  for (i=0; i<K; i++)
    index[i]=i;
  for (i=0; i<K-1; i++){
    j = rand()%(K-1);
    int tmp = index[i];
    index[i] = index[j];
    index[j] = tmp;
  }

  // read concrete input file
  char *inname = argv[1];
  FILE *infile = fopen(inname, "r");
  float cinput[N]; i=0;
  while(i<N && fscanf(infile, "%f", &cinput[i])>0){
      i++;
  }
  assert(i>=N);
  fclose(infile);
  // read concrete inputs
  for (i = 0; i < K-S; i++){
    w[index[i]]=cinput[index[i]];
  }
  for (i = 0; i<K; i++){
    data[i] = cinput[i+K];
  }

  // read symbolic inputs
  float sinput[K];
  klee_make_symbolic(&sinput, sizeof(sinput), "sinput");

  for (i = 0 ; i <S ; i++){
    w[index[K-1-i]] = sinput[i] / FLT_MAX * 100;
  }

  double wsd, wmean;

  wmean   = gsl_stats_wmean(w, 1, data, 1, K);
  wsd     = gsl_stats_wsd_m(w, 1, data, 1, K, wmean);

  return 0;
}

max_stars_repo_name: snipekill/FPGen

max_issues_repo_name: snipekill/FPGen

max_forks_repo_name: snipekill/FPGen


text:
/*
 * BaseMaster.h
 *
 *  Created on: Aug 7, 2014
 *      Author: thardin
 */

#ifndef BASEMASTER_H_
#define BASEMASTER_H_

#include "FMIClient.h"
#include <zmq.hpp>
#ifdef USE_GPL
#include <gsl/gsl_multiroots.h>
#include "common/fmigo_storage.h"
using namespace fmigo_storage;
#endif
#include <chrono>

namespace fmitcp_master {
    class BaseMaster {
        int rendezvous;
        //t1 = time at which to perform next step
        std::chrono::high_resolution_clock::time_point t1;

#ifdef USE_MPI
        std::string m_mpi_str;
#endif
    protected:
        std::vector<FMIClient*> m_clients;
        std::vector<WeakConnection> m_weakConnections;
        OutputRefsType clientWeakRefs;

        InputRefsValuesType initialNonReals;  //for loop solver

#ifdef USE_GPL
        class FmigoStorage m_fmigoStorage;//(std::vector<size_t>(0));
#endif

    public:

        zmq::socket_t rep_socket;
        bool initing, paused, running;
        bool zmqControl;
        int m_pendingRequests;
        double t; //current time

        explicit BaseMaster(zmq::context_t &context, std::vector<FMIClient*> clients, std::vector<WeakConnection> weakConnections);
        virtual ~BaseMaster() {
          info("%i rendezvous\n", rendezvous);
          int messages = 0;
          for(auto client: m_clients)
            messages += client->messages;
          info("%i messages\n", messages);
        }

#ifdef USE_GPL
        static int loop_residual_f(const gsl_vector *x, void *params, gsl_vector *f);

        inline FmigoStorage & get_storage(){return m_fmigoStorage;}
        void storage_alloc(const std::vector<FMIClient*> &clients){
            vector<size_t> states({});
            vector<size_t> indicators({});
            for(auto client: clients) {
                states.push_back(client->getNumContinuousStates());
                indicators.push_back(client->getNumEventIndicators());
            }
            get_storage().allocate_storage(states,indicators);
        }
#endif
      virtual std::string getFieldNames() const {return "";}
      virtual void writeFields(bool last, FILE *outfile) {}
        void solveLoops();
        virtual void prepare() {};
        virtual void runIteration(double t, double dt) = 0;

#define on(name) void name(FMIClient* slave) {}
        on(onSlaveInstantiated)
        on(onSlaveInitialized)
        on(onSlaveTerminated)
        on(onSlaveFreed)
        on(onSlaveStepped)
        on(onSlaveGotVersion)
        on(onSlaveSetReal)
        on(onSlaveGotState)
        on(onSlaveSetState)
        on(onSlaveFreedState)
        on(onSlaveDirectionalDerivative)

        //T is needed because func maybe of a function in Client (from which FMIClient is derived)
        void queueMessage(std::vector<FMIClient*> fmus, std::string str) {
            for (auto it = fmus.begin(); it != fmus.end(); it++) {
                (*it)->queueMessage(str);
            }
        }

        //like queueMessage() but only for one FMU
        void queueMessage(FMIClient *fmu, std::string str) {
            fmu->queueMessage(str);
        }

        //queueMessage() followed by wait() (blocking)
        void sendWait(std::vector<FMIClient*> fmus, std::string str) {
            queueMessage(fmus, str);
            wait();
        }

        //like sendWait() but only for one FMU (blocking)
        void sendWait(FMIClient *fmu, std::string str) {
            queueMessage(fmu, str);
            wait();
        }

        void wait();

        void resetT1();
        //wait for t1, advance t1 by timeStep
        void waitupT1(double timeStep);
        void handleZmqControl();

        void queueValueRequests();
        //if cset_set == true then only delete values for FMUs whose ID is in cset
        void deleteCachedValues(bool cset_set = false, const fmitcp::int_set& cset = fmitcp::int_set());
    };
};

#endif /* BASEMASTER_H_ */

max_stars_repo_name: Tjoppen/fmigo

max_issues_repo_name: Tjoppen/fmigo

max_forks_repo_name: Tjoppen/fmigo


text:
#pragma once

#include <gsl/pointers>  // contains gsl::not_null

class Consultant;

class PaymentCalculator {
 public:
  double calculate() const { return {}; }

  void setConsultant(const Consultant &c) { consultant_ = &consultant_; }

  void setTaxPercentage(double tax) { taxPercentage_ = tax; }

 private:
  gsl::not_null<const Consultant *> consultant_;
  double taxPercentage_;
};


text:
/* histogram/file.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

#include <config.h>
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_block.h>
#include <gsl/gsl_histogram.h>

int
gsl_histogram_fread (FILE * stream, gsl_histogram * h)
{
  int status = gsl_block_raw_fread (stream, h->range, h->n + 1, 1);

  if (status)
    return status;

  status = gsl_block_raw_fread (stream, h->bin, h->n, 1);
  return status;
}

int
gsl_histogram_fwrite (FILE * stream, const gsl_histogram * h)
{
  int status = gsl_block_raw_fwrite (stream, h->range, h->n + 1, 1);

  if (status)
    return status;

  status = gsl_block_raw_fwrite (stream, h->bin, h->n, 1);
  return status;
}

int
gsl_histogram_fprintf (FILE * stream, const gsl_histogram * h,
		       const char *range_format, const char *bin_format)
{
  size_t i;
  const size_t n = h->n;

  for (i = 0; i < n; i++)
    {
      int status = fprintf (stream, range_format, h->range[i]);

      if (status < 0)
	{
	  GSL_ERROR ("fprintf failed", GSL_EFAILED);
	}

      status = putc (' ', stream);

      if (status == EOF)
	{
	  GSL_ERROR ("putc failed", GSL_EFAILED);
	}

      status = fprintf (stream, range_format, h->range[i + 1]);

      if (status < 0)
	{
	  GSL_ERROR ("fprintf failed", GSL_EFAILED);
	}

      status = putc (' ', stream);

      if (status == EOF)
	{
	  GSL_ERROR ("putc failed", GSL_EFAILED);
	}

      status = fprintf (stream, bin_format, h->bin[i]);

      if (status < 0)
	{
	  GSL_ERROR ("fprintf failed", GSL_EFAILED);
	}

      status = putc ('\n', stream);

      if (status == EOF)
	{
	  GSL_ERROR ("putc failed", GSL_EFAILED);
	}
    }

  return GSL_SUCCESS;
}

int
gsl_histogram_fscanf (FILE * stream, gsl_histogram * h)
{
  size_t i;
  const size_t n = h->n;
  double upper;

  for (i = 0; i < n; i++)
    {
      int status = fscanf (stream,
			   "%lg %lg %lg", h->range + i, &upper,
			   h->bin + i);

      if (status != 3)
	{
	  GSL_ERROR ("fscanf failed", GSL_EFAILED);
	}
    }

  h->range[n] = upper;

  return GSL_SUCCESS;
}


text:
#include <gsl/gsl_math.h>
#include "gsl_cblas.h"
#include "cblas.h"

void
cblas_zgemv (const enum CBLAS_ORDER order, const enum CBLAS_TRANSPOSE TransA,
	     const int M, const int N, const void *alpha, const void *A,
	     const int lda, const void *X, const int incX, const void *beta,
	     void *Y, const int incY)
{
#define BASE double
#include "source_gemv_c.h"
#undef BASE
}


text:
/*
 *   dotProduct.c
 *
 *   An MPI program that computes the inner product of two vectors
 *
 *   Written by cetinsamet -*- cetin.samet@metu.edu.tr
 *   April, 2019
 *
 */
#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <cblas.h>

int main(int argc, char *argv[]) {

    int i;
    int rank, size;

    double minStart, maxEnd;
    double dotP, locDotP;

    MPI_Init(&argc, &argv);
    MPI_Comm_size(MPI_COMM_WORLD, &size);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);

    int ORDER	= atoi(argv[1]);
    int locSize = ORDER / size;

    double* locU;	locU = (double*) malloc(locSize * sizeof(double));
    double* locV;	locV = (double*) malloc(locSize * sizeof(double));

    for (i = 0; i < locSize; ++i) {
        locU[i] = i;
        locV[i] = i;
    }

    double pStart = MPI_Wtime();    // <-- Start Time
    locDotP = cblas_ddot(locSize, locU, 1, locV, 1);
    MPI_Reduce(&locDotP, &dotP, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
    double pEnd = MPI_Wtime();      // <-- Stop Time

    MPI_Reduce(&pStart,  &minStart, 1, MPI_DOUBLE, MPI_MIN, 0, MPI_COMM_WORLD);
    MPI_Reduce(&pEnd, &maxEnd, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD);

    if (rank == 0)
    	printf("Dot product: %f\tVector size: %d\tElapsed time: %f\n", dotP, ORDER, maxEnd - minStart);

    MPI_Finalize();
    return 0;
}

max_stars_repo_name: cetinsamet/dot-product

max_issues_repo_name: cetinsamet/dot-product

max_forks_repo_name: cetinsamet/dot-product


text:
#ifndef CTETRA_DOS_H
#define CTETRA_DOS_H

#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include "input.h"
#include "ecache.h"
#include "tetra.h"
#include "submesh.h"
#include "fermi.h"

double* Tetra_AllDosList(InputFn Efn, int na, int nb, int nc, int num_bands, gsl_matrix *R, double *Es, int num_dos);
double* Tetra_DosList(InputFn Efn, int na, int nb, int nc, int num_bands, gsl_matrix *R, double *Es, int num_dos);
double* Tetra_DosEnergyDerivList(InputFn Efn, int na, int nb, int nc, int num_bands, gsl_matrix *R, double *Es, int num_dos, double num_electrons, double *fermi, double *dos_fermi, double *dos_deriv_fermi);
double Tetra_TotalDos(double E, EnergyCache *Ecache);
double Tetra_TotalDosEnergyDeriv(double E, EnergyCache *Ecache);
double DosContrib(double E, double E1, double E2, double E3, double E4, double num_tetra);
double DosEnergyDerivContrib(double E, double E1, double E2, double E3, double E4, double num_tetra);

#endif // CTETRA_DOS_H

max_stars_repo_name: tflovorn/ctetra

max_issues_repo_name: tflovorn/ctetra

max_forks_repo_name: tflovorn/ctetra