text:
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>

// y' = 2y + x
// https://www.wolframalpha.com/input/?i=y%27+%3D+2y+%2B+x%2C+y%280%29+%3D+1%2C+y%281%29
// https://www.wolframalpha.com/input/?i=1%2F4+*+%285+*+e%5E2+-+3%29
// 미분 방정식을 푸는 방법
// 1. 미분 방정식의 표준형을 작성한다.
//    y' <<<--- 미분항의 계수를 1로 만들어야 표준형이됨
// 2. 표준형의 y항 계수식을 통해 적분 인자를 구한다.
// 3. 적분 인자를 양변에 곱한다.
//    2y' + y = 3 -> y' + y / 2 = 3 / 2
//    여기서 적분 인자는 1 / 2가 됨
//    좀 더 정확하게는 exp^integral(1/2) 형태가 됨
// 4. 양변에 대한 계산(적분)을 진행한다.
// 5. 양변을 exp^적분인자 로 나누어 y에 대해 정리한다.
int odefunc (double x, const double y[], double f[], void *params)
{
	// 계산하고자 하는 미분 방정식의 표준형식을 여기에 기록함
    f[0] = x+2*y[0];
    return GSL_SUCCESS;
}

int * jac;

int main(void)
{
	// 1계 미분 방정식
	// 만약 2계 미분 방정식이라면 dim = 2
    int dim = 1;
	// 미방을 풀기 위한 데이터 타입
    gsl_odeiv2_system sys = {odefunc, NULL, dim, NULL};

	// 드라이버는 미방을 쉽게 풀 수 있도록 시간에 따른 변화, 제어, 스텝에 대한 래퍼
	// gsl_odeiv2_step_rkf45는 룬게쿠타 4-5차를 의미함
	// Taylor(테일러) 혹은 Maculaurin(맥클로린) 급수
	// 테일러 급수를 활용해서 미분을 4번이상 하여
	// 오차를 최소화하는 방식을 룬게쿠타 4-5차라고 표현함

	// 다음 인자는 미분방정식의 시작지점
	// 절대적인 에러 바운드와 상대적인 에러 바운드
    gsl_odeiv2_driver * d = gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rkf45, 1e-6, 1e-6, 0.0);
    int i;
    double x0 = 0.0,  xf = 10.0; /* start and end of integration interval */
    double x = x0;
    double y[1] = { 1.0  };  /* initial value */

	// 루프를 돌면서 각 케이스별 미분 방정식을 계산함
    for (i = 1; i <= xf; i++)
    {
		// 0 ~ 10 - x input
        double xi = x0 + i * (xf-x0) / xf;
		// 결국 이 코드를 다시 재해석하자면
		// x가 0일 경우, 1일 경우, ... , 10일 경우에 대해 각각을 해석한다고 보면 됩니다!
		// 앞서서 룬게쿠타 4-5차 방식으로 미분 방정식을 계산하도록
		// odefunc을 만들었고 이에 대한 설정값들로 xi를 넣고 이에 대한 y값을 계산한다.
        int status = gsl_odeiv2_driver_apply (d, &x, xi, y);

        if (status != GSL_SUCCESS)
        {
            printf ("error, return value=%d\n", status);
            break;
        }

        printf ("%.8e %.8e\n", x, y[0]);
		printf ("x = %lf, xi = %lf\n", x, xi);
    }

	// malloc 혹은 new 처럼 Heap에 메모리를 할당하기 때문에
	// free를 사용해서 메모리 해제를 해주셔야 메모리 누수를 방지할 수 있음
    gsl_odeiv2_driver_free (d);
    return 0;
}


text:
#include <malloc.h>
#include <math.h>
#include <memory.h>
#include <stdio.h>
#include <string.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_matrix_double.h>
#include "tc_mat.h"

struct tc_mat *tc_mat_ctr(uint32_t nr_, uint32_t nc_)
{
	struct tc_mat *A = (struct tc_mat *) malloc(sizeof(struct tc_mat));
	memset(A, 0, sizeof(struct tc_mat));
	tc_mat_resize(A, nr_, nc_);
	return A;
}

void tc_mat_dtr(struct tc_mat *A)
{
	if (!A)
		return;
	tc_mat_clear(A);
	free(A);
}

void tc_mat_clear(struct tc_mat *A)
{
	if (!A)
		return;
	if (A->a) {
		for(uint32_t r=0; r < A->nr; r++)
			if (A->a[r])
				free(A->a[r]);
		free(A->a);
	}
	A->a = NULL;
	A->nr = 0;
	A->nc = 0;
}

void tc_mat_print(const struct tc_mat *A)
{
	if (!A)
		return;
	for(uint32_t i=0; i < A->nr; i++) {
		for (uint32_t j=0; j < A->nc; j++)
			printf(" %13.10f", A->a[i][j]);
		printf("\n");
	}
}

void tc_mat_resize(struct tc_mat *A, uint32_t nr_, uint32_t nc_)
{
	if ((A->nr == nr_) && (A->nc == nc_))
		return;
	tc_mat_clear(A);
	A->nr = nr_;
	A->nc = nc_;
	A->a = (double **) malloc(sizeof(double *) * nr_);
	for(uint32_t r=0; r < nr_; r++)
		A->a[r] = (double *) malloc(sizeof(double) * nc_);
}

void tc_mat_copy(struct tc_mat *B, const struct tc_mat *A)
{
	if (A == B)
		return;
	if ((B->nr != A->nr) || (B->nc != A->nc))
		tc_mat_resize(B, A->nr, A->nc);
	for(uint32_t i=0; i < A->nr; i++)
		for(uint32_t j=0; j < A->nc; j++)
			B->a[i][j] = A->a[i][j];
}

void tc_mat_identity(struct tc_mat *A)
{
	for(uint32_t i=0; i < A->nr; i++)
		for(uint32_t j=0; j < A->nc; j++)
			A->a[i][j] = (i==j)? 1.0: 0.0;
}

void tc_mat_transpose(struct tc_mat *B, const struct tc_mat *A)
{
	if (A == B) {
		struct tc_mat *b = tc_mat_ctr(A->nc, A->nr);
		for(uint32_t i=0; i < A->nr; i++)
			for(uint32_t j=0; j < A->nc; j++)
				b->a[j][i] = A->a[i][j];
		tc_mat_copy(B, b);
		tc_mat_dtr(b);
		return;
	}
	if ((B->nr != A->nc) || (B->nc != A->nr))
		tc_mat_resize(B, A->nc, A->nr);
	for(uint32_t i=0; i < A->nr; i++)
		for(uint32_t j=0; j < A->nc; j++)
			B->a[j][i] = A->a[i][j];
}

int tc_mat_add(struct tc_mat *C, const struct tc_mat *A, const struct tc_mat *B)
{
	if (!C || !A || !B || (A->nr != B->nr) || (A->nc != B->nc))
		return -1;
	if ((C->nr != A->nr) || (C->nc != A->nc))
		tc_mat_resize(C, A->nr, B->nc);
	for(uint32_t i=0; i < C->nr; i++)
		for(uint32_t j=0; j < C->nc; j++)
			C->a[i][j] = A->a[i][j] + B->a[i][j];
	return 0;
}

int tc_mat_mult(struct tc_mat *C, const struct tc_mat *A, const struct tc_mat *B)
{
	if (!C || !A || !B || (A->nc != B->nr))
		return -1;
	if ((C == A) || (C == B)) {
		struct tc_mat *c = tc_mat_ctr(A->nr, B->nr);
		tc_mat_mult(c, A, B);
		tc_mat_copy(C, c);
		tc_mat_dtr(c);
		return 0;
	}
	if ((C->nc != B->nc) || (C->nr != A->nr))
		tc_mat_resize(C, A->nr, B->nc);
	for(uint32_t i=0; i < C->nr; i++) {
		for(uint32_t j=0; j < C->nc; j++) {
			C->a[i][j] = 0.0;
			for(uint32_t k=0; k < A->nc; k++)
				C->a[i][j] += A->a[i][k] * B->a[k][j];
		}
	}
	return 0;
}

int tc_mat_mult_scalar(struct tc_mat *C, double a, const struct tc_mat *B)
{
	if (!C || !B)
		return -1;
	if ((C->nr != B->nr) || (C->nc != B->nc))
		tc_mat_resize(C, B->nr, B->nc);
	for(uint32_t i=0; i < C->nr; i++)
		for(uint32_t j=0; j < C->nc; j++)
			C->a[i][j] = a * B->a[i][j];
	return 0;
}

/* Householder Reduction to Bidiagonal Form
 * Input m x n matrix A
 * Output: m x n matrix  U, products of Householder matrices
 *         m x n matrix  B, bidiagonal
 *         m x n matrix  V, products of Householder matrices
 *   such that A = U B V^T
 *
 * A Householder matrix P is a reflection transformation of the form
 * I - tau u u^T where tau = 2 / (u^T u). P is orthonormal since
 *        P^T P = (I - tau u u^T)^T (I - tau u u^T)
 *              = (I - tau u u^T)   (I - tau u u^T)
 *              =  I - 2 tau u u^T + tau^2 u (u^T u) u^T
 *              =  I - 2 tau u u^T + tau^2 u (2 / tau) u^T
 *              =  I
 * Products of orthonomal matrices are orthonormal so U and V will be
 * orthonormal.
 *
 * B is constructed via the algorithm
 *   1.  B <-- A
 *   2.  U <-- I  m x n
 *   3.  V <-- I  n x n
 *   4.  For k = 1, ..., n
 *       a. construct a Householder matrix P such that
 *             i. P is the identity on the first k-1 components, and
 *            ii. the k-th column of P B is zero below the diagonal.
 *       b.  B <--- P B
 *       c.  U <--- U P
 *       d.  If k < n - 2, increment k and repeat steps a-c for the
 *              operation from the right (the row of B P will be zero
 *              to the right of the superdiagonal).
 *
 *  Step 4a(i) is equivalent to setting the first k-1 components of u
 *  to be zero. Let v be the k-th column of P B. Since
 *      P v = (I - tau u u^T) v
 *          = (I - tau u u^T) v
 *          = v - tau (u^T v) u   (Eq. 1)
 *  we have that step 4a(ii) is equivalent to setting u to be an appropriate
 *  scalar multiple of v for all components above k.  The previous steps
 *  would have zeroed out the first k-2 components of v. Let
 *     b^2 = Sum_{i>=k} (v_i)^2
 *  and choose the sign of b to ensure s = v_k - b is nonzero. Set u = v/s
 *  on all components larger than k and set u_k = 1. Then
 *    tau = 2 / (u^T u)
 *        = 2 / (1 + (b^2 - (v_k)^2) / s^2 )
 *        = (2 s) / (s + (b^2 - (v_k)^2) / s)
 *        = (2 s) / (s - (b + v_k))
 *        = (2 s) / ((v_k - b) - (b + v_k))
 *        = (2 s) / (- 2 b)
 *        =  -s / b
 *  This leads to the zeroing out of the lower components of P v since
 *    u^T v = v_k + (b^2 - (v_k)^2) / s
 *          = v_k - (b   +  v_k)
 *          = b
 *  meaning that Eq 1 is simply
 *      P v = v - (-s/b) (b) u
 *          = v + s u
 *  and so P v is zero on all components larger than k.
 *
 *  Side note on efficiency 1 (not implemented):
 *    Step c produces U and V as products of Housholder matrices. The actual
 *    computation of these products is not the most computationally efficient
 *    method. A Housholder matrix acts on a vector v as
 *          P v = (I - tau u u^T) v
 *              =  v - tau u (v^T u)^T  (Eq. 2)
 *    Thus if we've recorded all u's and tau's, we can find U, the product of
 *    all the P's, by sequentially applying Eq 2 to the components e1, e2, etc.
 *
 *  Side note on efficiency 2 (not implemented):
 *    Since each iteration of B leads to known places where its terms have zeros,
 *    all matrix multiplications involving B can be restricted to the nonzero
 *    terms.
 *
 *  Side note on efficiency 3 (not implemented):
 *    Since the u's have more zeros in each iteration, the collection of u's
 *    can be stored in the original matrix A (and hence B) so long as the
 *    Side note on efficiency #2 was implemented.
 */
int tc_mat_bidiag_decomp(
	const struct tc_mat *A,
	struct tc_mat *U,
	struct tc_mat *B,
	struct tc_mat *V)
{
	if (!A || !U || !B || !V)
		return -1;
	if ((U->nr != A->nr) || (U->nc != A->nr))
		tc_mat_resize(U, A->nr, A->nr);
	if ((V->nr != A->nc) || (V->nc != A->nc))
		tc_mat_resize(V, A->nc, A->nc);
	if ((B->nr != A->nr) || (B->nc != A->nc))
		tc_mat_resize(B, A->nr, A->nc);
	tc_mat_copy(B, A);
	tc_mat_identity(U);
	tc_mat_identity(V);
	struct tc_mat *Ic = tc_mat_ctr(A->nc, A->nc);
	tc_mat_identity(Ic);
	struct tc_mat *Ir = tc_mat_ctr(A->nr, A->nr);
	tc_mat_identity(Ir);
	struct tc_mat *u = tc_mat_ctr(0, 0);
	struct tc_mat *uT = tc_mat_ctr(0, 0);
	struct tc_mat *u_uT = tc_mat_ctr(0, 0);
	struct tc_mat *P = tc_mat_ctr(0, 0);
	for(uint32_t k=0; (k < A->nc) && (k+1 < A->nr); k++) {
		/* from the left */
		double b_k = B->a[k][k];
		double b = 0.0;
		for(uint32_t i=k; i < A->nr; i++)
			b += B->a[i][k] * B->a[i][k];
		if (b != 0.0) { /* if b = 0, there's nothing to do */
			b = ((b_k > 0)? -1.0: 1.0) * sqrt(b);
			double s = b_k - b;
			double tau = - s / b;
			tc_mat_resize(u, A->nr, 1);
			for(uint32_t i=0; i < A->nr; i++)
				u->a[i][0] = (i > k)? B->a[i][k]/s: ((i==k)? 1.0: 0.0);
			/*  P = I - tau u u^T */
			tc_mat_transpose(uT, u);
			tc_mat_mult(u_uT, u, uT);
			tc_mat_mult_scalar(P, -tau, u_uT);
			tc_mat_add(P, Ir, P);
			/*  U B = (U P) (P B) */
			tc_mat_mult(B, P, B);
			tc_mat_mult(U, U, P);
		}
		if (k + 2 < A->nc) {
			/* from the right */
			double b_k = B->a[k][k+1];
			double b = 0.0;
			for(uint32_t j=k+1; j < A->nc; j++)
				b += B->a[k][j] * B->a[k][j];
			b = ((b_k > 0)? -1.0: 1.0) * sqrt(b);
			if (b != 0.0) { /* if b = 0, there's nothing to do */
				double s = b_k - b;
				double tau = - s / b;
				tc_mat_resize(u, A->nc, 1);
				for(uint32_t i=0; i < A->nc; i++)
					u->a[i][0] = (i > k+1)? B->a[k][i]/s: ((i==k+1)? 1.0: 0.0);
				/*  P = I - tau u u^T */
				tc_mat_transpose(uT, u);
				tc_mat_mult(u_uT, u, uT);
				tc_mat_mult_scalar(P, -tau, u_uT);
				tc_mat_add(P, Ic, P);
				/*  B V = (B P) (P V) */
				tc_mat_mult(B, B, P);
				tc_mat_mult(V, P, V);
			}
		}
	}
	tc_mat_dtr(P);
	tc_mat_dtr(u);
	tc_mat_dtr(uT);
	tc_mat_dtr(u_uT);
	tc_mat_dtr(Ic);
	tc_mat_dtr(Ir);
	return 0;
}

/* tc_mat_eigenvalues
 * only implmented for 2x2 matrices
 * Input 2 x 2 matrix A
 * Output: 2 x 1 matrix E of eigenvalues
 *
 * If we set A = [ a b ]
 *               [ c d ]
 * Then
 *  0 = det(A - x I)
 *    = (a - x)(d - x) - bc
 *    = x^2 - (a + d) x + (ad - bc)
 */
int tc_mat_eigenvalues(struct tc_mat *E, const struct tc_mat *A)
{
	if (!A || !E)
		return -1;
	if ((A->nr != 2) || (A->nc != 2))
		return -1;
	if (E->nr != A->nr)
		tc_mat_resize(E, A->nr, 1);
	double a = A->a[0][0];
	double b = A->a[1][0];
	double c = A->a[0][1];
	double d = A->a[1][1];
	double ad2 = (a + d)/2.0;
	double det = ad2*ad2 - (a*d - b*c);
	double s = (det <= 0.0)? 0.0: sqrt(det);
	E->a[0][0] = ad2 + s;
	E->a[1][0] = ad2 - s;
	return 0;
}

/* Wilkinson Shift
 * Input: m x 2 matrix A
 * Output: the Wilkinson shift
 *
 * ei= trailing 2x2 matrix of D^T D */
double tc_mat_wilkinson_shift(const struct tc_mat *A)
{
	if (!A)
		return 0.0;
	if (A->nc != 2)
		return 0.0;
	struct tc_mat *AT = tc_mat_ctr(2, A->nr);
	tc_mat_transpose(AT, A);
	struct tc_mat *AT_A = tc_mat_ctr(2, 2);
	tc_mat_mult(AT_A, AT, A);
	struct tc_mat *E = tc_mat_ctr(2, 1);
	tc_mat_eigenvalues(E, AT_A);
	double l0 = E->a[0][0];
	double l1 = E->a[1][0];
	double t22 = AT_A->a[1][1];
	tc_mat_dtr(AT);
	tc_mat_dtr(AT_A);
	tc_mat_dtr(E);
	return (fabs(l0 - t22) < fabs(l1 - t22))? l0: l1;
}

int tc_mat_svd(
	const struct tc_mat *A,
	struct tc_mat *U,
	struct tc_mat *D,
	struct tc_mat *V)
{
	if (!U || !D || !V || !A)
		return -1;
	if ((U->nr != A->nr) || (U->nc != A->nc))
		tc_mat_resize(U, A->nr, A->nc);
	if ((V->nr != A->nc) || (V->nc != A->nc))
		tc_mat_resize(V, A->nc, A->nc);
	if ((D->nr != A->nc) || (D->nc != A->nc))
		tc_mat_resize(D, A->nc, A->nc);
	if (A->nr < A->nc)
		return -2;
	/* case nc = 1 */
	if (A->nc == 1) {
		V->a[0][0] = 1.0;
		double normsq = 0.0;
		for(uint32_t i=0; i < A->nr; i++)
			normsq += A->a[i][0] * A->a[i][0];
		D->a[0][0] = sqrt(normsq);
		if (D->a[0][0] == 0.0)
			tc_mat_copy(U, A);
		else {
			double inv = 1.0 / D->a[0][0];
			for(uint32_t i=0; i < A->nr; i++)
				U->a[i][0] = A->a[i][0] * inv;
		}
		return 0;
	}
	/* Step1: A = U D V where D is bi-diagonal
	 *            A is nr x nc
	 *            U is nr x nr
	 *            D is nr x nc
	 *            V is nc x nc
	 *
	 *    If A was not a square matrix, change U and D
	 *            U  nr x nc
	 *            D  nc x nc
	 */
	tc_mat_bidiag_decomp(A, U, D, V);
	if (D->nr != D->nc) {
		struct tc_mat *U0 = tc_mat_ctr(A->nr, A->nc);
		for(uint32_t i=0; i < A->nr; i++)
			for(uint32_t j=0; j < A->nc; j++)
				U0->a[i][j] = U->a[i][j];
		tc_mat_resize(U, A->nr, A->nc);
		tc_mat_copy(U, U0);
		tc_mat_dtr(U0);

		struct tc_mat *D0 = tc_mat_ctr(A->nc, A->nc);
		for(uint32_t i=0; i < A->nc; i++)
			for(uint32_t j=0; j < A->nc; j++)
				D0->a[i][j] = D->a[i][j];
		tc_mat_resize(D, A->nc, A->nc);
		tc_mat_copy(D, D0);
		tc_mat_dtr(D0);
	}
	/* iterate until either max_iterations has been hit
	 * or the largest superdiagonal entry is below the
	 * threshold. Set threshold = 1-e8 * largest element.
	 * Any term less than 0.1 threshold will be considered
	 * to be zero.
	 */
	const uint32_t max_iterations = 100 * D->nc;
	double threshold = D->a[0][0];
	for(uint32_t i=0; i < D->nc; i++)
		if (fabs(D->a[i][i]) > threshold)
			threshold = fabs(D->a[i][i]);
	for(uint32_t i=0; i+1 < D->nc; i++)
		if (fabs(D->a[i][i+1]) > threshold)
			threshold = fabs(D->a[i][i+1]);
	threshold *= 1e-12;
	const double zero_threshold = 0.1 * threshold;
	/* Always reindex the components so that
     * D so that has the form
	 *     [D1, 0] where D1 is diagonal and
	 *     [ 0,D2]       D2 is bidiagonal
	 * Let i0 be the starting index of D2
	 */
	uint32_t i0 = 0;
	for(uint32_t I=0; I < max_iterations; I++)
	{
		/* For any zeros on the diagonal, apply a series of
		 * Givens rotations to also zero out its off-diagonal
		 * term.  Then move this component into D1.
		 */
		for(uint32_t i1=i0; i1 < D->nr; i1++) {
			if (fabs(D->a[i1][i1]) > zero_threshold)
				continue;
			if ((i1+1 == D->nr) && fabs(D->a[i1-1][i1]) > zero_threshold)
				continue;
			if ((i1+1 < D->nr) && (fabs(D->a[i1][i1+1]) > zero_threshold)) {
				for(uint32_t i=i1; i+1 < D->nr; i++) {
					/* U D = (U G) (G^T D)  */
					double alpha = D->a[i1][i+1];
					if (fabs(alpha) < zero_threshold)
						break;
					double beta = D->a[i+1][i+1];
					double gamma = sqrt(alpha*alpha + beta*beta);
					double c = beta / gamma;
					double s = alpha / gamma;
					for(uint32_t j=0; j < D->nr; j++) {
						double a = D->a[i1][j];
						double b = D->a[i+1][j];
						D->a[i1][j] = a*c - b*s;
						D->a[i+1][j] = a*s + b*c;
					}
					for(uint32_t j=0; j < U->nr; j++) {
						double a = U->a[j][i1];
						double b = U->a[j][i+1];
						U->a[j][i1] = a*c - b*s;
						U->a[j][i+1] = a*s + b*c;
					}
				}
				i1++;
			}
			/* move (i0,i1-1) down, i1 -> i0  */
			for(uint32_t j=0; j < D->nr; j++) {
				double tmp = V->a[i1][j];
				for(uint32_t k=i1; k > i0; k--)
					V->a[k][j] = V->a[k-1][j];
				V->a[i0][j] = tmp;
			}
			for(uint32_t j=0; j < D->nr; j++) {
				double tmp = U->a[j][i1];
				for(uint32_t k=i1; k > i0; k--)
					U->a[j][k] = U->a[j][k-1];
				U->a[j][i0] = tmp;
			}
			double tmp = D->a[i1][i1];
			double tmp1 = D->a[i1][i1+1];
			for(uint32_t k=i1; k > i0; k--) {
				D->a[k][k] = D->a[k-1][k-1];
				if (k+1 < D->nc)
					D->a[k][k+1] = D->a[k-1][k];
			}
			D->a[i0][i0] = tmp;
			D->a[i0][i0+1] = tmp1;
			i0++;
		}
		/* For any zeros on the superdiagonal, move the
		 * component to D1.
		 */
		for(uint32_t i=i0; i+1 < D->nr; i++) {
			if (fabs(D->a[i][i+1]) >= zero_threshold)
				continue;
			if (i == i0) {
				i0++;
				continue;
			}
			if (i+2 != D->nr)
				continue;
			/* move (i0,i) down, i+1 -> i0    */
			for(uint32_t j=0; j < D->nr; j++) {
				double tmp = V->a[i+1][j];
				for(uint32_t k=i+1; k > i0; k--)
					V->a[k][j] = V->a[k-1][j];
				V->a[i0][j] = tmp;
			}
			for(uint32_t j=0; j < D->nr; j++) {
				double tmp = U->a[j][i+1];
				for(uint32_t k=i+1; k > i0; k--)
					U->a[j][k] = U->a[j][k-1];
				U->a[j][i0] = tmp;
			}
			double tmp = D->a[i+1][i+1];
			double tmp1 = D->a[i][i+1];
			for(uint32_t k=i+1; k > i0; k--) {
				D->a[k][k] = D->a[k-1][k-1];
				D->a[k-1][k] = (k-1==i0)? tmp1: D->a[k-2][k-1];
			}
			D->a[i0][i0] = tmp;
			i0++;
			if ((i+2 == D->nr) && (i0 != i))
				i--; /* retry last term */
		}
		if (i0 >= D->nr - 1)
			break;
		/* Find largest element on super diagonal.*/
		/* Break if less than threshold.          */
		double largest_off_diag = 0.0;
		uint32_t n_zeros_on_off_diagonal = 0;
		for(uint32_t i=i0; i+1 < D->nr; i++) {
			if (fabs(D->a[i][i+1]) > largest_off_diag)
				largest_off_diag = fabs(D->a[i][i+1]);
			if (fabs(D->a[i][i+1]) < zero_threshold)
				n_zeros_on_off_diagonal++;
		}
		if (largest_off_diag < threshold)
			break;
		/* Find the next zero on the uper diagonal*/
		/* Treat [i0,i1] as a block.              */
		uint32_t i1 = i0;
		for(  ; i1+1 < D->nr; i1++)
			if (fabs(D->a[i1][i1+1]) < zero_threshold)
				break;
		/* Find Wilkerson shift.                  */
		struct tc_mat *t = tc_mat_ctr(3, 2);
		for(uint32_t i=0; i < 3; i++)
			for(uint32_t j=0; j < 2; j++)
				t->a[i][j] = (i1+i>2)? D->a[i1+i-2][i1+j-1]: 0.0;
		double mu = tc_mat_wilkinson_shift(t);
		tc_mat_dtr(t);
		double alpha = D->a[i0][i0] * D->a[i0][i0] - mu;
		double beta = D->a[i0][i0] * D->a[i0][i0+1];
		/* Apply Givens rotations G from i0 to the bottom,
		 * chasing the nonzero element until off the matrix
		 */
		for(uint32_t i=i0; i < i1; i++) {
			/* D V = (D G) (G^T V)  */
			double gamma = sqrt(alpha*alpha + beta*beta);
			double c = alpha / gamma;
			double s = -beta / gamma;
			for(uint32_t j=0; j < D->nr; j++) {
				double a = D->a[j][i+0];
				double b = D->a[j][i+1];
				D->a[j][i+0] = a*c - b*s;
				D->a[j][i+1] = a*s + b*c;
			}
			for(uint32_t j=0; j < D->nc; j++) {
				double a = V->a[i+0][j];
				double b = V->a[i+1][j];
				V->a[i+0][j] = a*c - b*s;
				V->a[i+1][j] = a*s + b*c;
			}
			/* U D = (U G) (G^T D) */
			alpha = D->a[i+0][i];
			beta = D->a[i+1][i];
			gamma = sqrt(alpha*alpha + beta*beta);
			if (fabs(gamma) > 0.0) {
				c = alpha / gamma;
				s = -beta / gamma;
				for(uint32_t j=0; j < D->nc; j++) {
					double a = D->a[i+0][j];
					double b = D->a[i+1][j];
					D->a[i+0][j] = a*c - b*s;
					D->a[i+1][j] = a*s + b*c;
				}
				for(uint32_t j=0; j < U->nr; j++) {
					double a = U->a[j][i+0];
					double b = U->a[j][i+1];
					U->a[j][i+0] = a*c - b*s;
					U->a[j][i+1] = a*s + b*c;
				}
			}
			if (i + 2 < D->nr) {
				alpha = D->a[i][i+1];
				beta = D->a[i][i+2];
			}
		}
	}
	/* now swap components to order the diagonal terms
	 * from largest to smallest
	 */
	for(uint32_t i=0; i+1 < D->nr; i++) {
		double largest = fabs(D->a[i][i]);
		uint32_t largest_i = i;
		for(uint32_t j=i+1; j < D->nr; j++) {
			if (fabs(D->a[j][j]) > largest) {
				largest_i = j;
				largest = fabs(D->a[j][j]);
			}
		}
		if (largest_i != i) {
			for(uint32_t j=0; j < D->nr; j++) {
				double tmp = V->a[i][j];
				V->a[i][j] = V->a[largest_i][j];
				V->a[largest_i][j] = tmp;
				tmp = U->a[j][i];
				U->a[j][i] = U->a[j][largest_i];
				U->a[j][largest_i] = tmp;
			}
			double tmp = D->a[i][i];
			D->a[i][i] = D->a[largest_i][largest_i];
			D->a[largest_i][largest_i] = tmp;
		}
		if (D->a[i][i] < 0) {
			D->a[i][i] = -D->a[i][i];
			for(uint32_t j=0; j < D->nr; j++)
				V->a[i][j] = -V->a[i][j];
		}
	}
	/* just to be sure, zero out all off-diagonal terms of D */
	for(uint32_t i=0; i < D->nr; i++)
		for(uint32_t j=0; j < D->nc; j++)
			if (i != j)
				D->a[i][j] = 0.0;
	/* transpose V */
	tc_mat_transpose(V, V);
	return 0;
}

ext: c

max_stars_repo_name: endolith/Truthcoin

max_issues_repo_name: endolith/Truthcoin

max_forks_repo_name: endolith/Truthcoin


text:
/*
 * Copyright (c) 2016-2021 lymastee, All rights reserved.
 * Contact: lymastee@hotmail.com
 *
 * This file is part of the gslib project.
 *
 * 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

#ifndef polygon_9cb4a5c6_8903_42fd_9232_eb81bfcd30a5_h
#define polygon_9cb4a5c6_8903_42fd_9232_eb81bfcd30a5_h

#include <gslib/std.h>
#include <gslib/math.h>
#include <ariel/config.h>

__ariel_begin__

class polygon
{
public:
    typedef gs::vector<point3> vtstream;
    typedef gs::vector<int> idstream;
    typedef vtstream::iterator vtiter;
    typedef vtstream::const_iterator const_vtiter;
    typedef idstream::iterator iditer;
    typedef idstream::const_iterator const_iditer;

public:
    polygon();
    int get_vertex_count() const { return (int)_vtstream.size(); }
    int get_index_count() const { return (int)_idstream.size(); }
    int get_face_count() const { return _faces; }
    int query_face_count();
    void write_vtstream(int start, int cnt, const point3* src);
    void append_vtstream(int cnt, const point3* src) { write_vtstream(get_vertex_count(), cnt, src); }
    void append_vtstream(const point3& src) { _vtstream.push_back(src); }
    void write_idstream(int start, int cnt, const int* src);
    void append_idstream(int cnt, const int* src) { write_idstream(get_index_count(), cnt, src); }
    void append_idstream(int idx) { _idstream.push_back(idx); }
    point3& get_vertex(int i) { return _vtstream.at(i); }
    const point3& get_vertex(int i) const { return _vtstream.at(i); }
    int& get_index(int i) { return _idstream.at(i); }
    const int& get_index(int i) const { return _idstream.at(i); }

protected:
    vtstream    _vtstream;
    idstream    _idstream;
    int         _faces;
};

__ariel_end__

#endif

ext: h

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


text:
#pragma once

#include <gsl/gsl>
#include <DirectXMath.h>
#include <vector>
#include <memory>
#include <string>
#include <cstdint>

namespace Library
{
	class Model;
	class Bone;
	class Keyframe;
	class OutputStreamHelper;
	class InputStreamHelper;

	struct BoneAnimationData final
	{
		std::uint32_t BoneIndex{ 0 };
		std::vector<std::shared_ptr<Keyframe>> Keyframes;
	};

    class BoneAnimation final
    {
    public:
		BoneAnimation(Model& model, InputStreamHelper& streamHelper);
		BoneAnimation(Model& model, const BoneAnimationData& boneAnimationData);
		BoneAnimation(Model& model, BoneAnimationData&& boneAnimationData);
		BoneAnimation(const BoneAnimation&) = default;
		BoneAnimation(BoneAnimation&& rhs) = default;
		BoneAnimation& operator=(const BoneAnimation& rhs) = default;
		BoneAnimation& operator=(BoneAnimation&& rhs) = default;
		~BoneAnimation() = default;

		Bone& GetBone();
		std::vector<std::shared_ptr<Keyframe>>& Keyframes();

		std::uint32_t GetTransform(float time, DirectX::XMFLOAT4X4& transform) const;
		void GetTransformAtKeyframe(std::uint32_t keyframeIndex, DirectX::XMFLOAT4X4& transform) const;
		void GetInteropolatedTransform(float time, DirectX::XMFLOAT4X4& transform) const;

		void Save(OutputStreamHelper& streamHelper);

    private:
		void Load(InputStreamHelper& streamHelper);
		std::uint32_t FindKeyframeIndex(float time) const;

		Model* mModel;
		std::weak_ptr<Bone> mBone;
		std::vector<std::shared_ptr<Keyframe>> mKeyframes;
    };
}


text:
#ifndef TOURNAMENTPIVOTING_H
#define TOURNAMENTPIVOTING_H

#include <stdio.h>
#include <stdlib.h>
#include "SuiteSparseQR_C.h"
#include <string.h>
#include <math.h>
#include <mpi.h>

#ifdef MKL
  #include <mkl.h>
#else
  #include <lapacke.h>
#endif

#define ASSERT(t_) if(!(t_)) printf("WARNING %s, line %d\n",__FILE__,__LINE__)
#ifndef MAX
#define MAX(a_,b_) ((a_)>(b_)?(a_):(b_))
#endif
#define DBL_EPSILON 1e-15

int preAlps_tournamentPivoting(MPI_Comm comm, int *colPtr, int *rowInd, double *a, int m,  int n,  int nnz, long col_offset, int k, long *Jc,
   double *Sval, int printSVal, int ordering);

int preAlps_tournamentPivotingQR(MPI_Comm comm, int *colPtr, int *rowInd, double *a, int m,  int n,  int nnz,
  long col_offset, int k, long *Jc, double *Sval, int printSVal, int checkFact, int printFact, int ordering);

int preAlps_tournamentPivotingCUR(MPI_Comm comm, int *colPtr, int *rowInd, double *a, int m,  int n,  int nnz,
    long col_offset, int k, long *Jr, long *Jc, double *Sval, int printSVal, int checkFact, int printFact, int ordering);

#endif

max_stars_repo_name: W-Wuxian/preAlps

max_issues_repo_name: W-Wuxian/preAlps

max_forks_repo_name: W-Wuxian/preAlps


text:
/*
 * To change this license header, choose License Headers in Project Properties.
 * To change this template file, choose Tools | Templates
 * and open the template in the editor.
 */

/*
 * File:   stdheader.h
 * Author: lars
 *
 * Created on May 19, 2020, 2:40 PM
 */

#ifndef STDHEADER_H
#define STDHEADER_H
#include <string>
#include <vector>
#include <random>
#include <chrono>
#include <limits>
#include <iostream>
#include <fstream>
#include <stdio.h>
#include <cstddef>
#include <regex>
#include <cmath>
#include <cstdlib>
#include <filesystem>
#include <tuple>
#include <map>
#include <utility>
#include <iomanip>
#include <type_traits>
#include <set>
#include <experimental/filesystem>
#include <filesystem>
#include <time.h>
#include <unordered_set>
#include <unordered_map>
//simple parallel
#include <parallel/algorithm>
//#include <experimental/algorithm>

#include <nlohmann/json.hpp>
#include <gsl/gsl_statistics_double.h>
#include "CLI11.hpp"
#include <essentia/algorithmfactory.h>
#include <essentia/streaming/algorithms/poolstorage.h>
#include <essentia/streaming/algorithms/fileoutput.h>
#include <essentia/pool.h>
#include <essentia/essentiamath.h>

#include "statInterface.h"
#include "utilities.h"
#include "featureSelektor.h"
#include "AudioFeature.h"
#include "GLCM.h"
#include "FACM.h"
#include "Haralick.h"
#include "AudioFile.h"
#include "distanceMetric.h"
#include "essentiaInterface.h"
#include "Json_Handler.h"
#include "printObjects.h"
#include "CLIHandler.h"
#include "src/hierarchCluster/fastcluster.h"

#include "simpleTests.h"
#endif /* STDHEADER_H */

max_stars_repo_name: Quodoso/SampleOrganiser

max_issues_repo_name: Quodoso/SampleOrganiser

max_forks_repo_name: Quodoso/SampleOrganiser


text:
// Copyright (c) 2021 Stig Rune Sellevag
//
// This file is distributed under the MIT License. See the accompanying file
// LICENSE.txt or http://www.opensource.org/licenses/mit-license.php for terms
// and conditions.

#ifndef SCILIB_LINALG_EIGENVALUE_H
#define SCILIB_LINALG_EIGENVALUE_H

#ifdef USE_MKL
#include <mkl.h>
#else
#include <lapacke.h>
#endif

#include <scilib/mdarray.h>
#include <scilib/linalg_impl/lapack_types.h>
#include <exception>
#include <cassert>
#include <complex>
#include <type_traits>

namespace Sci {
namespace Linalg {

// Compute eigenvalues and eigenvectors of a real symmetric matrix.
template <class Layout>
inline void eigs(Sci::Matrix_view<double, Layout> a,
                 Sci::Vector_view<double, Layout> w,
                 double abstol = -1.0 /* use default value */)
{
    static_assert(a.is_contiguous());
    static_assert(w.is_contiguous());

    assert(a.extent(0) == a.extent(1));
    assert(w.extent(0) == a.extent(0));

    const BLAS_INT n = static_cast<BLAS_INT>(a.extent(0));
    const BLAS_INT nselect = n;
    const BLAS_INT lda = n;
    const BLAS_INT ldz = nselect;

    BLAS_INT il = 1;
    BLAS_INT iu = n;
    BLAS_INT m;
    BLAS_INT info;

    double vl = 0.0;
    double vu = 0.0;

    Sci::Vector<BLAS_INT, Layout> isuppz(2 * n);
    Sci::Matrix<double, Layout> z(ldz, n);

    auto matrix_layout = LAPACK_ROW_MAJOR;
    if constexpr (std::is_same_v<Layout, stdex::layout_left>) {
        matrix_layout = LAPACK_COL_MAJOR;
    }

    info = LAPACKE_dsyevr(matrix_layout, 'V', 'A', 'U', n, a.data(), lda, vl,
                          vu, il, iu, abstol, &m, w.data(), z.data(), ldz,
                          isuppz.data());
    if (info != 0) {
        throw std::runtime_error("dsyevr failed");
    }
    Sci::copy(z.view(), a);
}

template <class Layout, class Allocator>
inline void eigs(Sci::Matrix<double, Layout, Allocator>& a,
                 Sci::Vector<double, Layout, Allocator>& w,
                 double abstol = -1.0 /* use default value */)
{
    eigs(a.view(), w.view(), abstol);
}

// Compute eigenvalues and eigenvectors of a real non-symmetric matrix.
template <class Layout>
void eig(Sci::Matrix_view<double, Layout> a,
         Sci::Matrix_view<std::complex<double>, Layout> evec,
         Sci::Vector_view<std::complex<double>, Layout> eval)
{
    using namespace Sci;

    static_assert(a.is_contiguous());
    static_assert(evec.is_contiguous());
    static_assert(eval.is_contiguous());

    assert(a.extent(0) == a.extent(1));
    assert(a.extent(0) == eval.extent(0));
    assert(a.extent(0) == evec.extent(0));
    assert(a.extent(1) == evec.extent(1));

    const BLAS_INT n = static_cast<BLAS_INT>(a.extent(1));

    Sci::Vector<double, Layout> wr(n);
    Sci::Vector<double, Layout> wi(n);
    Sci::Matrix<double, Layout> vr(n, n);
    Sci::Matrix<double, Layout> vl(n, n);

    auto matrix_layout = LAPACK_ROW_MAJOR;
    if constexpr (std::is_same_v<Layout, stdex::layout_left>) {
        matrix_layout = LAPACK_COL_MAJOR;
    }
    BLAS_INT info =
        LAPACKE_dgeev(matrix_layout, 'N', 'V', n, a.data(), n, wr.data(),
                      wi.data(), vl.data(), n, vr.data(), n);
    if (info != 0) {
        throw std::runtime_error("dgeev failed");
    }
    for (BLAS_INT i = 0; i < n; ++i) {
        std::complex<double> wii(wr(i), wi(i));
        eval(i) = wii;
        BLAS_INT j = 0;
        while (j < n) {
            if (wi(j) == 0.0) {
                evec(i, j) = std::complex<double>{vr(i, j), 0.0};
                ++j;
            }
            else {
                evec(i, j) = std::complex<double>{vr(i, j), vr(i, j + 1)};
                evec(i, j + 1) = std::complex<double>{vr(i, j), -vr(i, j + 1)};
                j += 2;
            }
        }
    }
}

template <class Layout,
          class Allocator_a,
          class Allocator_evec,
          class Allocator_eval>
void eig(Sci::Matrix<double, Layout, Allocator_a>& a,
         Sci::Matrix<std::complex<double>, Layout, Allocator_evec>& evec,
         Sci::Vector<std::complex<double>, Layout, Allocator_eval>& eval)
{
    eig(a.view(), evec.view(), eval.view());
}

} // namespace Linalg
} // namespace Sci

#endif // SCILIB_LINALG_EIGENVALUE_H

max_stars_repo_name: stigrs/scilib

max_issues_repo_name: stigrs/scilib

max_forks_repo_name: stigrs/scilib


text:
static char help[] =
"Solves a structured-grid Stokes problem with DMDA and SNES using finite elements: Q^2-P^-1.\n"
"FIXME: START OVER\n"
"\n\n";

#include <petsc.h>

typedef struct {
  double u;
  double v;
  double p;
} Field;

typedef struct {
  DM        da;
  Vec       xexact;
  double    L,     // length of domain in x direction
            H,     // length of domain in y direction
            g1,    // signed component of gravity in x-direction
            g2,    // signed component of gravity in y-direction
            mu;    // viscosity
} AppCtx;

extern PetscErrorCode FormExactSolution(AppCtx*);
extern PetscErrorCode FormFunctionLocal(DMDALocalInfo*,Field**,Field**,AppCtx*);
//extern PetscErrorCode FormJacobianLocal(DMDALocalInfo*,Field**,Mat,Mat,AppCtx*);

int main(int argc,char **argv)
{
  PetscErrorCode ierr;
  AppCtx         user;                         /* user-defined work context */
  SNES           snes;                         /* nonlinear solver */
  Vec            x;                            /* solution vector */

  PetscInitialize(&argc,&argv,(char*)0,help);

  const double rg = 1000.0 * 9.81, // = rho g; scale for body force;
                                      //     kg / (m^2 s^2);  weight of water
                  theta = PETSC_PI / 9.0; // 20 degrees
  user.L     = 1.0;
  user.H     = 1.0;
  user.g1    = rg * sin(theta);
  user.g2    = - rg * cos(theta);
  user.nu    = 0.25;  // Pa s;  typical dynamic viscosity of motor oil
  user.ppeps = 1.0;

// FIXME START OVER

  PetscBool doerror = PETSC_FALSE, exactinit = PETSC_FALSE;
  ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "stokes_", "options for stokes", ""); CHKERRQ(ierr);
  FIXME
  ierr = PetscOptionsEnd(); CHKERRQ(ierr);

  ierr = DMDACreate2d(PETSC_COMM_WORLD,
                      FIXME); CHKERRQ(ierr);
  ierr = DMSetFromOptions(user.da); CHKERRQ(ierr);
  ierr = DMSetUp(user.da); CHKERRQ(ierr);
  ierr = DMDASetUniformCoordinates(user.da, 0.0, user.H, 0.0, user.L, -1.0, -1.0); CHKERRQ(ierr);
  ierr = DMSetApplicationContext(user.da,&user); CHKERRQ(ierr);
  ierr = DMCreateGlobalVector(user.da,&x); CHKERRQ(ierr);

  ierr = DMDASetFieldName(user.da,0,"u"); CHKERRQ(ierr);
  ierr = DMDASetFieldName(user.da,1,"v"); CHKERRQ(ierr);
  ierr = DMDASetFieldName(user.da,2,"p"); CHKERRQ(ierr);

  FIXME

  SNESDestroy(&snes);  DMDestroy(&user.da);
  return PetscFinalize();
}

PetscErrorCode FormExactSolution(AppCtx* user) {
  PetscErrorCode ierr;
  DMDALocalInfo  info;
  int       i,j;
  double      hy, y;
  Field          **ax;

  ierr = DMDAGetLocalInfo(user->da,&info); CHKERRQ(ierr);
  hy = user->H / (double)(info.my - 1);
  ierr = DMDAVecGetArray(user->da,user->xexact,&ax); CHKERRQ(ierr);
  for (j = info.ys; j < info.ys+info.ym; j++) {
    y = hy * (double)j;
    for (i = info.xs; i < info.xs+info.xm; i++) {
      ax[j][i].u = (user->g1 / user->nu) * y * (user->H - y/2.0);
      ax[j][i].v = 0.0;
      ax[j][i].p = - user->g2 * (user->H - y);
    }
  }
  ierr = DMDAVecRestoreArray(user->da,user->xexact,&ax); CHKERRQ(ierr);
  return 0;
}


text:
/* cdf/exponential.c
 *
 * Copyright (C) 2003, 2007 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 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_cdf.h>

/* The exponential distribution has the form

   p(x) dx = exp(-x/mu) dx/mu

   for x = 0 ... +infty */

double
gsl_cdf_exponential_P (const double x, const double mu)
{
  if (x < 0)
    {
      return 0;
    }
  else
    {
      double P = -expm1 (-x / mu);

      return P;
    }
}

double
gsl_cdf_exponential_Q (const double x, const double mu)
{
  if (x < 0)
    {
      return 1;
    }
  else
    {
      double Q = exp (-x / mu);

      return Q;
    }
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <sys/time.h>
#include <sys/resource.h>
#include <unistd.h>

#ifdef GSL_FOUND
#include <gsl/gsl_integration.h>
#endif

#include "core_allvars.h"
#include "core_init.h"
#include "core_mymalloc.h"
#include "core_cool_func.h"
//#include "core_metal_yield.h"

/* These functions do not need to be exposed externally */
double integrand_time_to_present(const double a, void *param);
void set_units(void);
void read_snap_list(const int ThisTask);
double time_to_present(const double z);
struct table{
        double tbl[1000][17];
        int nr;
};
struct table read_table(char *fname, int ncols);
void read_metal_yield();

#ifdef HDF5
#include "io/io_save_hdf5.h"
#endif

void init(const int ThisTask)
{
    run_params.Age = mymalloc(ABSOLUTEMAXSNAPS*sizeof(*(run_params.Age)));

    set_units();

    read_snap_list(ThisTask);

    //Hack to fix deltaT for snapshot 0
    //This way, galsnapnum = -1 will not segfault.
    run_params.Age[0] = time_to_present(1000.0);//lookback time from z=1000
    run_params.Age++;

    for(int i = 0; i < run_params.Snaplistlen; i++) {
        run_params.ZZ[i] = 1 / run_params.AA[i] - 1;
        run_params.Age[i] = time_to_present(run_params.ZZ[i]);  //this is actually a lookback time
	run_params.lbtime[i] = (run_params.Age[0] - run_params.Age[i]) * run_params.UnitTime_in_s / SEC_PER_MEGAYEAR; //age of universe since bigbang
    }

    run_params.a0 = 1.0 / (1.0 + run_params.Reionization_z0);
    run_params.ar = 1.0 / (1.0 + run_params.Reionization_zr);

    read_cooling_functions();
    if(ThisTask == 0) {
        printf("cooling functions read\n");
    }

    read_metal_yield();
    if(ThisTask == 0){
	printf("metal yields read\n\n");
    }

#if 0
#ifdef HDF5
    if(HDF5Output) {
        calc_hdf5_props();
    }
#endif
#endif

}

void set_units(void)
{

    run_params.UnitTime_in_s = run_params.UnitLength_in_cm / run_params.UnitVelocity_in_cm_per_s;
    run_params.UnitTime_in_Megayears = run_params.UnitTime_in_s / SEC_PER_MEGAYEAR;
    run_params.G = GRAVITY / pow(run_params.UnitLength_in_cm, 3) * run_params.UnitMass_in_g * pow(run_params.UnitTime_in_s, 2);
    run_params.UnitDensity_in_cgs = run_params.UnitMass_in_g / pow(run_params.UnitLength_in_cm, 3);
    run_params.UnitPressure_in_cgs = run_params.UnitMass_in_g / run_params.UnitLength_in_cm / pow(run_params.UnitTime_in_s, 2);
    run_params.UnitCoolingRate_in_cgs = run_params.UnitPressure_in_cgs / run_params.UnitTime_in_s;
    run_params.UnitEnergy_in_cgs = run_params.UnitMass_in_g * pow(run_params.UnitLength_in_cm, 2) / pow(run_params.UnitTime_in_s, 2);

    run_params.EnergySNcode = run_params.EnergySN / run_params.UnitEnergy_in_cgs * run_params.Hubble_h;
    run_params.EtaSNcode = run_params.EtaSN * (run_params.UnitMass_in_g / SOLAR_MASS) / run_params.Hubble_h;

    // convert some physical input parameters to internal units
    run_params.Hubble = HUBBLE * run_params.UnitTime_in_s;

    // compute a few quantitites
    run_params.RhoCrit = 3.0 * run_params.Hubble * run_params.Hubble / (8 * M_PI * run_params.G);
}

void read_snap_list(const int ThisTask)
{
    char fname[MAX_STRING_LEN+1];

    snprintf(fname, MAX_STRING_LEN, "%s", run_params.FileWithSnapList);
    FILE *fd = fopen(fname, "r");
    if(fd == NULL) {
        printf("can't read output list in file '%s'\n", fname);
        ABORT(0);
    }

    run_params.Snaplistlen = 0;
    do {
        if(fscanf(fd, " %lg ", &(run_params.AA[run_params.Snaplistlen])) == 1) {
            run_params.Snaplistlen++;
        } else {
            break;
        }
    } while(run_params.Snaplistlen < run_params.MAXSNAPS);
    fclose(fd);

    if(ThisTask == 0) {
        printf("found %d defined times in snaplist\n", run_params.Snaplistlen);
    }
}

double time_to_present(const double z)
{
    const double end_limit = 1.0;
    const double start_limit = 1.0/(1 + z);
    double result=0.0;
#ifdef GSL_FOUND
#define WORKSIZE 1000
    gsl_function F;
    gsl_integration_workspace *workspace;
    double abserr;

    workspace = gsl_integration_workspace_alloc(WORKSIZE);
    F.function = &integrand_time_to_present;

    gsl_integration_qag(&F, start_limit, end_limit, 1.0 / run_params.Hubble,
                        1.0e-9, WORKSIZE, GSL_INTEG_GAUSS21, workspace, &result, &abserr);

    gsl_integration_workspace_free(workspace);

#undef WORKSIZE
#else
    /* Do not have GSL - let's integrate numerically ourselves */
    const double step  = 1e-7;
    const int64_t nsteps = (end_limit - start_limit)/step;
    result = 0.0;
    const double y0 = integrand_time_to_present(start_limit + 0*step, NULL);
    const double yn = integrand_time_to_present(start_limit + nsteps*step, NULL);
    for(int64_t i=1; i<nsteps; i++) {
        result  += integrand_time_to_present(start_limit + i*step, NULL);
    }

    result = (step*0.5)*(y0 + yn + 2.0*result);
#endif

    /* convert into Myrs/h (I think -> MS 23/6/2018) */
    const double time = 1.0 / run_params.Hubble * result;

    // return time to present as a function of redshift
    return time;
}

double integrand_time_to_present(const double a, void *param)
{
    (void) param;
    return 1.0 / sqrt(run_params.Omega / a + (1.0 - run_params.Omega - run_params.OmegaLambda) + run_params.OmegaLambda * a * a);
}

void read_metal_yield(void)
{
        char fname[MAX_STRING_LEN];
	int cols=0;
	int i, j, rows;

	/* READ AGB YIELD */
	if (run_params.AGBYields == 0)
	{
		double Z_std[7] = {0.0, 1e-4, 4e-4, 4e-3, 8e-3, 0.02, 0.05}; //metal grid from Karakas 10
		strcpy(fname, "src/auxdata/yields/table2d.dat");
		cols = 13; //number of columns in the file

        	struct table data;
        	data = read_table(fname, cols);
        	rows = data.nr;
		double Z[rows];

        	//count how many rows necessary for each column based on Z
        	int count = 0;
        	for (i=0; i<rows; i++)
        	{
                	Z[i] = data.tbl[i][0];
                	if (Z[i] == Z_std[0]) {
                        	count++;
			}
        	}

        	//now assign
		run_params.countagb = count;

        	for (j=0; j<7; j++){
                	int index = 0;
                	for (i=0; i<rows; i++){
                        	if (Z[i] == Z_std[j])
                        	{
                                	run_params.magb[index] = data.tbl[i][1]; //1 is the column number
                                	run_params.qCagb[index][j] = data.tbl[i][6] + data.tbl[i][7] + data.tbl[i][11];
                                	run_params.qNagb[index][j] = data.tbl[i][8] + data.tbl[i][12];
                                	run_params.qOagb[index][j] = data.tbl[i][9];
                                	run_params.Qagb[index][j] = run_params.qCagb[index][j] + run_params.qNagb[index][j] + run_params.qOagb[index][j];
                                	index++;
                        	}
                	}
        	}
	}

	/* READ SN II YIELD */
	if (run_params.SNIIYields == 0)
        {
		double Z_std[METALGRID] = {0.0, 1e-4, 4e-4, 4e-3, 8e-3, 0.02, 0.05}; //metal grid Woosley & Weaver 1995
		strcpy(fname, "src/auxdata/yields/table4a.dat");
                cols = 17; //number of columns in the file

		struct table data;
		data = read_table(fname, cols);
		rows = data.nr;
		double Z[rows];

		//count how many rows necessary for each column based on Z
		int count = 0;
        	for (i=0; i<rows; i++)
        	{
                	Z[i] = data.tbl[i][0];
                	if (Z[i] == Z_std[0]) {
                        	count++;}
        	}

		run_params.countsn = count;

        	for (j=0; j<7; j++){
                	int index = 0;
                	for (i=0; i<rows; i++){
                	        if (Z[i] == Z_std[j])
                        	{
                                	run_params.msn[index] = data.tbl[i][1];
                                	run_params.qCsn[index][j] = data.tbl[i][6] + data.tbl[i][15];
                                	run_params.qOsn[index][j] = data.tbl[i][7];
                                	run_params.qMgsn[index][j] = data.tbl[i][9];
                                	run_params.qSisn[index][j] = data.tbl[i][10];
                                	run_params.qSsn[index][j] = data.tbl[i][11];
                                	run_params.qCasn[index][j] = data.tbl[i][12];
                                	run_params.qFesn[index][j] = data.tbl[i][13];
                                	run_params.Qsn[index][j] = run_params.qCsn[index][j] + run_params.qOsn[index][j] + run_params.qMgsn[index][j] + run_params.qSisn[index][j] + run_params.qSsn[index][j] + run_params.qCasn[index][j] + run_params.qFesn[index][j];
					index++;
                        	}
                	}
        	}
	}

	else if (run_params.SNIIYields == 1) //from Nomoto et al. 2006
	{
		double Z_std[METALGRID] = {0.0, 0.001, 0.004, 0.02};
		strcpy(fname, "src/auxdata/yields/Nomoto.dat");
		cols = 13; //number of columns in the file

		struct table data;
		data = read_table(fname, cols);
		rows = data.nr;
		double Z[rows];

		//count how many rows necessary for each column based on Z
		int count = 0;
		for (i=0; i<rows; i++)
		{
			Z[i] = data.tbl[i][0];
			if (Z[i] == Z_std[0]) {
				count++;}
		}

		run_params.countsn = count;

		for (j=0; j<4; j++){
			int index = 0;
			for (i=0; i<rows; i++){
				if (Z[i] == Z_std[j])
				{
					run_params.msn[index] = data.tbl[i][1];
					run_params.qCsn[index][j] = data.tbl[i][2] + data.tbl[i][11];
					run_params.qOsn[index][j] = data.tbl[i][3];
					run_params.qMgsn[index][j] = data.tbl[i][5];
					run_params.qSisn[index][j] = data.tbl[i][6];
					run_params.qSsn[index][j] = data.tbl[i][7];
					run_params.qCasn[index][j] = data.tbl[i][8];
					run_params.qFesn[index][j] = data.tbl[i][9];
					run_params.Qsn[index][j] = run_params.qCsn[index][j] + run_params.qOsn[index][j] + run_params.qMgsn[index][j] + run_params.qSisn[index][j] + run_params.qSsn[index][j] + run_params.qCasn[index][j] + run_params.qFesn[index][j];
					index++;
				}
			}
		}
	}

        /* SNIa YIELD */
	if (run_params.SNIaYields == 0)
	{
		run_params.qCrsnia = 0.0168; //yields from Iwamoto 1999
		run_params.qFesnia = 0.587;
		run_params.qNisnia = 0.0314;
	}

	else if (run_params.SNIaYields == 1) //yields from Seitenzahl et al. 2013
	{
		run_params.qCrsnia = 0.00857;
		run_params.qFesnia = 0.622;
		run_params.qNisnia = 0.069;
	}
	printf("ends reading yieldi \n");
}

struct table read_table(char *fname, int ncols)
{
        int i, j;
        struct table dt;
        FILE *file;

        file = fopen(fname,"r");
	fscanf(file, "%*[^\n]"); //skips the first line

        //start reading the file
        i=0;
        while (!feof(file)) //while it is not at the end of the file
    {
        for (j=0; j<ncols; j++)
        {
            fscanf(file, " %lf", &dt.tbl[i][j]); //notice the space char in front of %lf to handle any number of spaces as delimiter
        }

        fscanf(file,"%*[^\n]"); //ignore the rest of the line (if ncols is less than the file's column size), and go to the next line
        i++;
    }

        fclose(file);

        //store number of rows
        dt.nr = i-1; //i is at the line AFTER the last line, so it's must be subtracted by 1

    return dt;
}

max_stars_repo_name: dptriani/sage

max_issues_repo_name: dptriani/sage

max_forks_repo_name: dptriani/sage


text:
#ifndef INTEGRATIONS_H
#define INTEGRATIONS_H
#include <armadillo>
#include <cmath>
#include <gsl/gsl_math.h>
#include <iostream>

template <typename F> class gsl_function_pp : public gsl_function {
public:
  explicit gsl_function_pp(const F &func) : gsl_function_struct(), _func(func) {
    function = &gsl_function_pp::invoke;
    params = this;
  }

private:
  const F &_func;
  static double invoke(double x, void *params) {
    return static_cast<gsl_function_pp *>(params)->_func(x);
  }
};

// wrapper functions for integration routines of the gsl
double cquad(const std::function<double(double)> &f, double a, double b,
             double relerr, double epsabs);
double qags(const std::function<double(double)> &f, double a, double b,
            double relerr, double epsabs);
double qagiu(const std::function<double(double)> &f, double a, double relerr,
             double epsabs);

#endif // INTEGRATIONS_H

max_stars_repo_name: QuaCaTeam/quaca

max_issues_repo_name: myoelmy/quaca

max_forks_repo_name: myoelmy/quaca


text:
#ifndef RADFIELD_H
#define RADFIELD_H

#include <cstdio>
#include <gsl/gsl_integration.h>
#include "types.h"
#include "sn3d.h"

namespace radfield
{
  void zero_estimators(int modelgridindex);
  void init(int my_rank, int ndo, int ndo_nonempty);
  void initialise_prev_titer_photoionestimators(void);
  void write_to_file(int modelgridindex, int timestep);
  void close_file(void);
  __host__ __device__ void update_estimators(int modelgridindex, double distance_e_cmf, double nu_cmf, const PKT *pkt_ptr, double t_current);
  void increment_lineestimator(int modelgridindex, int lineindex, double increment);
  __host__ __device__ double radfield(double nu, int modelgridindex);
  __host__ __device__ double dbb_mgi(double nu, int modelgridindex);
  void fit_parameters(int modelgridindex, int timestep);
  __host__ __device__ void set_J_normfactor(int modelgridindex, double normfactor);
  __host__ __device__ void normalise_J(int modelgridindex, double estimator_normfactor_over4pi);
  __host__ __device__ void normalise_nuJ(int modelgridindex, double estimator_normfactor_over4pi);
  __host__ __device__ double get_T_R_from_J(int modelgridindex);
  __host__ __device__ int get_Jblueindex(int lineindex);
  __host__ __device__ double get_Jb_lu(int modelgridindex, int jblueindex);
  __host__ __device__ int get_Jb_lu_contribcount(int modelgridindex, int jblueindex);
  void titer_J(int modelgridindex);
  void titer_nuJ(int modelgridindex);
  void reduce_estimators(void);
  void do_MPI_Bcast(int modelgridindex, int root, int root_node_id);
  void write_restart_data(FILE *gridsave_file);
  void read_restart_data(FILE *gridsave_file);
  __host__ __device__ void normalise_bf_estimators(int modelgridindex, double estimator_normfactor_over_H);
  double get_bfrate_estimator(int element, int lowerion, int lower, int phixstargetindex, int modelgridindex);
  void print_bfrate_contributions(int element, int lowerion, int lower, int phixstargetindex, int modelgridindex, double nnlowerlevel, double nnlowerion);
  void reset_bfrate_contributions(const int modelgridindex);
  int integrate(
    const gsl_function *f, double nu_a, double nu_b, double epsabs, double epsrel,
    size_t limit, int key, gsl_integration_workspace *workspace, double *result, double *abserr);

  __host__ __device__ inline double dbb(double nu, float T, float W)
  // returns J_nu for a dilute black body [ergs/s/sr/cm2/Hz]
  {
    return W * TWOHOVERCLIGHTSQUARED * pow(nu, 3) / expm1(HOVERKB * nu / T);
  }
}

#endif //RADFIELD_H

max_stars_repo_name: artis-mcrt/artis

max_issues_repo_name: artis-mcrt/artis

max_forks_repo_name: artis-mcrt/artis


text:
#pragma once

#include <memory>
#include <string>
#include <vector>
#include <typeinfo>
#include <type_traits>

#include <gsl/gsl>

#include <upca/config.h>

#include "arch/arch.h"

namespace upca {

template <typename BACKEND> class resolver {
  class description {

    std::string name_;
    unsigned size_;
    unsigned offset_;
    decltype(typename BACKEND::resolver_type().resolve({})) data_;

  public:
    description(std::string name, const unsigned size, const unsigned offset,
                const typename BACKEND::resolver_type &resolver)
        : name_(std::move(name)), size_(size), offset_(offset),
          data_(resolver.resolve(name_)) {}
    description(const char *name, const unsigned size, const unsigned offset,
                const typename BACKEND::resolver_type &resolver)
        : name_(name), size_(size), offset_(offset),
          data_(resolver.resolve(name_)) {}

    const std::string &name() const { return name_; }
    unsigned size() const { return size_; }
    unsigned offset() const { return offset_; }
    decltype(auto) data() const { return data_; }
  };

  std::vector<description> counters_;
  typename BACKEND::resolver_type resolver_;

public:
  using value_type = description;
  using iterator = typename std::vector<description>::iterator;
  using const_iterator = typename std::vector<description>::const_iterator;

  void add(const char *const name, const unsigned size = 8) {
    counters_.emplace_back(name, size, counters_.size(), resolver_);
  }

  void add(std::string name, const unsigned size = 8) {
    counters_.emplace_back(std::move(name), size, counters_.size(), resolver_);
  }

  size_t size() const { return counters_.size() + 1; }

  unsigned bytesize() const {
    unsigned sum = 0;
    for (const auto &pmc : counters_) {
      sum += pmc.size();
    }
    return sum;
  }

  const description &
  at(const typename std::vector<description>::size_type pos) const {
    return counters_.at(pos);
  }
  const_iterator begin() const { return counters_.cbegin(); }
  const_iterator end() const { return counters_.cend(); }
  const_iterator cbegin() const { return counters_.cbegin(); }
  const_iterator cend() const { return counters_.cend(); }

  /* Get PMU backend. Thread has to be pinned to a CPU before this is called */
  std::unique_ptr<BACKEND> configure(const unsigned additional_counters = 0) const {
    /* meh. add +1 for the "out-of-band" cycles counter :/ */
    return std::make_unique<BACKEND>(counters_, additional_counters + 1);
  }
};

/* sample pmcs */
} // namespace upca

max_stars_repo_name: hannesweisbach/ucpa

max_issues_repo_name: hannesweisbach/ucpa

max_forks_repo_name: hannesweisbach/ucpa


text:
/* -*- linux-c -*- */
/* fewbody_utils.c

   Copyright (C) 2002-2004 John M. Fregeau

   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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

#include <stdio.h>
#include <stddef.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_roots.h>
#include "fewbody.h"

/* allocate a vector */
double *fb_malloc_vector(int n)
{
	return((double *) malloc(n * sizeof(double)));
}

/* allocate a matrix */
double **fb_malloc_matrix(int nr, int nc)
{
	int i;
	double **m;

	m = (double **) malloc(nr * sizeof(double *));

	m[0] = (double *) malloc(nr * nc * sizeof(double));

	for (i=1; i<nr; i++) {
		m[i] = m[i-1] + nc;
	}

	return(m);
}

/* free a vector */
void fb_free_vector(double *v)
{
	free(v);
}

/* free a matrix */
void fb_free_matrix(double **m)
{
	free(m[0]);
	free(m);
}

/* a fast square function */
double fb_sqr(double x)
{
	return(x*x);
}

/* a fast cube function */
double fb_cub(double x)
{
	return(x*x*x);
}

/* the dot product of two vectors */
double fb_dot(double x[3], double y[3])
{
	return(x[0] * y[0] + x[1] * y[1] + x[2] * y[2]);
}

/* the modulus of a vector */
double fb_mod(double x[3])
{
	return(sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]));
}

/* the cross product of two vectors */
int fb_cross(double x[3], double y[3], double z[3])
{
	z[0] = x[1] * y[2] - x[2] * y[1];
	z[1] = x[2] * y[0] - x[0] * y[2];
	z[2] = x[0] * y[1] - x[1] * y[0];

	return(0);
}

/* something to calculate the angular momentum */
int fb_angmom(fb_obj_t *star, int nstar, double L[3])
{
	int i, k;
	double **l;

	l = fb_malloc_matrix(nstar, 3);

	for (k=0; k<3; k++) {
		L[k] = 0.0;
	}

	for (i=0; i<nstar; i++) {
		fb_cross(star[i].x, star[i].v, l[i]);
		for (k=0; k<3; k++) {
			L[k] += star[i].m * l[i][k];
		}
	}

	fb_free_matrix(l);

	return(0);
}

/* something to calculate the internal angular momentum */
void fb_angmomint(fb_obj_t *star, int nstar, double L[3])
{
	int i, k;

	for (k=0; k<3; k++) {
		L[k] = 0.0;
	}

	for (i=0; i<nstar; i++) {
		for (k=0; k<3; k++) {
			L[k] += star[i].Lint[k];
		}
	}
}

/* something to calculate the internal energy */
double fb_einttot(fb_obj_t *star, int nstar)
{
	int i;
	double eint=0.0;

	for (i=0; i<nstar; i++) {
		eint += star[i].Eint;
	}

	return(eint);
}

/* calculates the total potential energy of the system */
double fb_petot(fb_obj_t *star, int nstar)
{
	int i, j, k;
	double pe=0.0, r[3];

	for (i=0; i<nstar; i++) {
		for (j=i+1; j<nstar; j++) {
			for (k=0; k<3; k++) {
				r[k] = star[j].x[k] - star[i].x[k];
			}
			pe += -star[i].m * star[j].m / fb_mod(r);
		}
	}

	return(pe);
}

/* calculates the total kinetic energy of the system */
double fb_ketot(fb_obj_t *star, int nstar)
{
	int i;
	double ke=0.0;

	for (i=0; i<nstar; i++) {
		ke += 0.5 * star[i].m * fb_dot(star[i].v, star[i].v);
	}

	return(ke);
}

/* calculates the potential energy of the bound members of the system */
double fb_outerpetot(fb_obj_t **obj, int nobj)
{
	int i, j, k;
	double pe=0.0, r[3];

	for (i=0; i<nobj; i++) {
		for (j=i+1; j<nobj; j++) {
			for (k=0; k<3; k++) {
				r[k] = obj[j]->x[k] - obj[i]->x[k];
			}
			pe += -obj[i]->m * obj[j]->m / fb_mod(r);
		}
	}

	return(pe);
}

/* calculates the kinetic energy of the bound members of the system */
double fb_outerketot(fb_obj_t **obj, int nobj)
{
	int i;
	double ke=0.0;

	for (i=0; i<nobj; i++) {
		ke += 0.5 * obj[i]->m * fb_dot(obj[i]->v, obj[i]->v);
	}

	return(ke);
}

/* solve the Kepler equation for the eccentric anomaly, given the mean anomaly and eccentricity */
double fb_kepler(double e, double mean_anom)
{
	int status, iter;
	double ecc_anom, params[2];
	gsl_function F;
	const gsl_root_fsolver_type *T;
	gsl_root_fsolver *s;

	/* set up the root solver */
	F.function = &fb_keplerfunc;
	F.params = &params;

	/* set the parameters */
	params[0] = e;
	params[1] = mean_anom;

	T = gsl_root_fsolver_brent;
	s = gsl_root_fsolver_alloc(T);
	gsl_root_fsolver_set(s, &F, 0.0, 2.0*FB_CONST_PI);

	/* get eccentric anomaly by root-finding */
	iter = 0;
	do {
		iter++;
		gsl_root_fsolver_iterate(s);
		status = gsl_root_test_interval(gsl_root_fsolver_x_lower(s), gsl_root_fsolver_x_upper(s), \
						FB_ROOTSOLVER_ABS_ACC, FB_ROOTSOLVER_REL_ACC);
	} while (status == GSL_CONTINUE && iter < FB_ROOTSOLVER_MAX_ITER);

	if (iter >= FB_ROOTSOLVER_MAX_ITER) {
		fprintf(stderr, "Root finder failed to converge.\n");
		exit(1);
	}

	/* we've got the root */
	ecc_anom = gsl_root_fsolver_root(s);

	/* free memory associated with root solver */
	gsl_root_fsolver_free(s);

	return(ecc_anom);
}

/* the Kepler function for the root finder */
double fb_keplerfunc(double ecc_anom, void *params)
{
	double e, mean_anom;

	e = ((double *)params)[0];
	mean_anom = ((double *)params)[1];

	return(ecc_anom - e * sin(ecc_anom) - mean_anom);
}

/* calculate the relative tidal acceleration */
double fb_reltide(fb_obj_t *bin, fb_obj_t *single, double r)
{
	double arel, atid;

	arel = bin->m / fb_sqr(bin->a * (1.0 + bin->e));
	/* Factor in numerator is (single->m + bin->m) instead of single->m for the case of
	   small single->m, for which we want to at some point resolve "bin" to get the
	   motion of "single".  This will eventually be fixed when we we use full collapsed
	   binary member positions in the integration, and not just the CM of the binary. */
	atid = 2.0 * (single->m + bin->m) / fb_cub(r) * bin->a * (1.0 + bin->e);

	return(atid/arel);
}

max_stars_repo_name: gnodvi/cosmos

max_issues_repo_name: gnodvi/cosmos

max_forks_repo_name: gnodvi/cosmos


text:
/*
 * ParticleFilter.h
 *
 *  Created on: Mar 21, 2012
 *      Author: nair
 */

#ifndef PARTICLEFILTER_H_
#define PARTICLEFILTER_H_

#include <QSettings>
#include <QFileInfo>

#include <opencv2/opencv.hpp>
#include <Eigen/Dense>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

namespace tracking
{

class ParticleFilterParams
{
public:
    ParticleFilterParams();
    ~ParticleFilterParams();

    void parse(std::string filename);

public:

    int NUMBER_OF_PARTICLES;

    /* standard deviations for gaussian sampling in transition model */
    float TRANS_X_STD;
    float TRANS_Y_STD;
    float TRANS_Z_STD;

    /* autoregressive dynamics parameters for transition model */
    float DYNAMICS_A1;
    float DYNAMICS_A2;
    float DYNAMICS_B0;

    bool DYNAMIC_MODEL_CV;

    float LAMBDA;
};

struct Particle
{
    float x;          /**< current x coordinate */
    float y;          /**< current y coordinate */
    float z;          /**< current z coordinate */
    float xp;         /**< previous x coordinate */
    float yp;         /**< previous y coordinate */
    float zp;         /**< previous z coordinate */
    float x0;         /**< original x coordinate */
    float y0;         /**< original y coordinate */
    float z0;			/**< original z coordinate */
    float x_velocity; /**< current x_velocity coordinate */
    float y_velocity; /**< current y_velocity coordinate */
    float z_velocity; /**< current z_velocity coordinate */

    float w;          /**< weight */
};

class ParticleFilter
{
public:

    ParticleFilter(uint& filter_id);
    ~ParticleFilter();

    void initialize(const cv::Vec3d& initPose);

    void predict();

    void predict_with_simple_model();

    void predict_with_cv_model();

    void correct();

    void resizeParticleSet(int n);

    void setParticleWeightFromDistanceMeasure(int id, float dist);

    void setParticleWeight(int id, float w);

    int nOfParticles()
    {
        return mParams.NUMBER_OF_PARTICLES;
    }

    static int particle_cmp(const void* p1,const void* p2 )
    {
        Particle* _p1 = (Particle*)p1;
        Particle* _p2 = (Particle*)p2;

        if( _p1->w > _p2->w )
            return -1;
        if( _p1->w < _p2->w )
            return 1;
        return 0;
    }

    void getParticlePose(int id, cv::Vec3d& pose);

    void getOutputPose(cv::Vec3d& pose);

    void getOutputVelocity(cv::Vec3d& vel);

    uint getFilterId();

    float getLikelihood();

private:

    void normalize();

    void resample();

private:

    ParticleFilterParams mParams;

    Particle* mParticles;

    cv::Vec3d mOutputPose;

    Eigen::MatrixXd A_Matrix;
    Eigen::MatrixXd G_Matrix;
    Eigen::MatrixXd W_Matrix;

    Eigen::MatrixXd State_Matrix;
    Eigen::MatrixXd State_Prev_Matrix;

    gsl_rng* rng;

    uint filterId;
    float mMaxLikelihood;
};

}
#endif /* PARTICLEFILTER_H_ */

max_stars_repo_name: cognitivesystems/smartcamera

max_issues_repo_name: cognitivesystems/smartcamera

max_forks_repo_name: cognitivesystems/smartcamera


text:
/**
 *
 * @file testing_cgebrd.c
 *
 *  PLASMA testing 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 c Tue Jan  7 11:45:17 2014
 *
 **/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

#include <plasma.h>
#include <cblas.h>
#include <lapacke.h>
#include <core_blas.h>
#include "testing_cmain.h"

#undef REAL
#define COMPLEX

static int check_orthogonality(int, int, int, PLASMA_Complex32_t*, int, float);
static int check_reduction(int, int, float*, PLASMA_Complex32_t*, int, PLASMA_Complex32_t*, int, PLASMA_Complex32_t*, int, float);
static int check_solution(int, float*, float*, float);

int testing_cgebrd(int argc, char **argv)
{
    int tree = 0;

    if ( argc < 1 ){
        goto usage;
    } else {
        tree = atoi(argv[0]);
    }

    /* Check for number of arguments*/
    if ( ((tree == 0) && (argc != 4)) ||
         ((tree != 0) && (argc != 5)) ){
      usage:
        USAGE("GEBRD", "MODE M N LDA [RH]",
              "   - MODE : 0: flat, 1: tree (RH needed)\n"
              "   - M    : number of rows of the matrix A\n"
              "   - N    : number of columns of the matrix A\n"
              "   - LDA  : leading dimension of the matrix A\n"
              "   - RH   : Size of each subdomains\n");
        return -1;
    }

    int M   = atoi(argv[1]);
    int N   = atoi(argv[2]);
    int LDA = atoi(argv[3]);
    int rh;
    if ( tree ) {
        rh = atoi(argv[4]);
        PLASMA_Set(PLASMA_HOUSEHOLDER_MODE, PLASMA_TREE_HOUSEHOLDER);
        PLASMA_Set(PLASMA_HOUSEHOLDER_SIZE, rh);
    }

    if (LDA < M){
        printf("LDA should be >= M !\n");
        return -1;
    }

    float eps  = LAPACKE_slamch_work('e');
    float dmax = 1.0;
    PLASMA_enum jobu  = PlasmaVec;
    PLASMA_enum jobvt = PlasmaVec;
    int info_orthou    = 0;
    int info_orthovt   = 0;
    int info_solution  = 0;
    int info_reduction = 0;
    int minMN = min(M, N);
    int mode  = 4;
    float rcond = (float) minMN;
    int INFO=-1;

    PLASMA_Complex32_t *A1   = (PLASMA_Complex32_t *)malloc(LDA*N*sizeof(PLASMA_Complex32_t));
    float *S1               = (float *)            malloc(minMN*sizeof(float));
    float *S2               = (float *)            malloc(minMN*sizeof(float));
    float *E2               = (float *)            malloc(minMN*sizeof(float));
    PLASMA_Complex32_t *work = (PLASMA_Complex32_t *)malloc(3*max(M, N)* sizeof(PLASMA_Complex32_t));
    PLASMA_Complex32_t *A2 = NULL;
    PLASMA_Complex32_t *U  = NULL;
    PLASMA_Complex32_t *VT = NULL;
    PLASMA_desc *T;

    /* Check if unable to allocate memory */
    if ( (!A1) || (!S1) || (!S2) || (!work) ) {
        printf("Out of Memory \n ");
        return -2;
    }

    PLASMA_Disable(PLASMA_AUTOTUNING);
    PLASMA_Set(PLASMA_TILE_SIZE, 120);
    PLASMA_Set(PLASMA_INNER_BLOCK_SIZE, 20);

    PLASMA_Enable(PLASMA_WARNINGS);
    PLASMA_Enable(PLASMA_ERRORS);
    PLASMA_Alloc_Workspace_cgebrd(M, N, &T);

    /*----------------------------------------------------------
    *  TESTING CGEBRD
    */
    /* Initialize A1 */
    LAPACKE_clatms_work( LAPACK_COL_MAJOR, M, N,
                         lapack_const(PlasmaDistUniform), ISEED,
                         lapack_const(PlasmaNonsymPosv), S1, mode, rcond,
                         dmax, M, N,
                         lapack_const(PlasmaNoPacking), A1, LDA, work );
    free(work);

    /* Copy A1 for check */
    if ( (jobu == PlasmaVec) && (jobvt == PlasmaVec) ) {
        A2 = (PLASMA_Complex32_t *)malloc(LDA*N*sizeof(PLASMA_Complex32_t));
        LAPACKE_clacpy_work(LAPACK_COL_MAJOR, 'A', M, N, A1, LDA, A2, LDA);
    }
    if ( jobu == PlasmaVec ) {
        U = (PLASMA_Complex32_t *)malloc(M*M*sizeof(PLASMA_Complex32_t));
        LAPACKE_claset_work(LAPACK_COL_MAJOR, 'A', M, M, 0., 1., U, M);
    }
    if ( jobvt == PlasmaVec ) {
        VT = (PLASMA_Complex32_t *)malloc(N*N*sizeof(PLASMA_Complex32_t));
        LAPACKE_claset_work(LAPACK_COL_MAJOR, 'A', N, N, 0., 1., VT, N);
    }

    /* PLASMA CGEBRD */
    INFO = PLASMA_cgebrd(jobu, jobvt, M, N, A1, LDA, S2, E2, T, U, M, VT, N);
    if(INFO!=0){
            printf(" ERROR OCCURED INFO %d\n",INFO);
            goto fin;
    }

    printf("\n");
    printf("------ TESTS FOR PLASMA CGEBRD ROUTINE -------  \n");
    printf("        Size of the Matrix %d by %d\n", M, N);
    printf("\n");
    printf(" The matrix A is randomly generated for each test.\n");
    printf("============\n");
    printf(" The relative machine precision (eps) is to be %e \n",eps);
    printf(" Computational tests pass if scaled residuals are less than 60.\n");
    /* call SVD solver using lapack routine for our resulting bidiag [S E] */
    int NCVT  = 0;
    int NRU   = 0;
    int NCC   = 0;
    if (jobu == PlasmaVec)
        NRU = M;
    if (jobvt == PlasmaVec)
        NCVT = N;
    PLASMA_enum uplo = M >=N ? PlasmaUpper : PlasmaLower;
    INFO = LAPACKE_cbdsqr(LAPACK_COL_MAJOR, lapack_const(uplo),
                            minMN, NCVT, NRU, NCC,
                            S2, E2,
                            VT, N, U, M, NULL, 1 );
    /* Check the orthogonality, reduction and the singular values */
    if ( jobu == PlasmaVec )
        info_orthou = check_orthogonality(PlasmaLeft, M, M, U, M, eps);

    if ( jobvt == PlasmaVec )
        info_orthovt = check_orthogonality(PlasmaRight, N, N, VT, N, eps);

    /*
     * WARNING: For now, Q is associated to Band tridiagonal reduction and
     * not to the final tridiagonal reduction, so we can not call the check
     * Need to accumulate the orthogonal transformations
     * during the bulge chasing to be able to perform the next test!
     */
    if ( (jobu == PlasmaVec) && (jobvt == PlasmaVec) )
        info_reduction = check_reduction(M, N, S2, A2, LDA, U, M, VT, N, eps);

    info_solution = check_solution(minMN, S1, S2, eps);

    if ( (info_solution == 0) & (info_orthou == 0) &
         (info_orthovt == 0) & (info_reduction == 0) ) {
        if (M >= N) {
           printf("***************************************************\n");
           printf(" ---- TESTING CGEBRD .. M >= N ........... PASSED !\n");
           printf("***************************************************\n");
        }
        else {
           printf("***************************************************\n");
           printf(" ---- TESTING CGEBRD .. M < N ............ PASSED !\n");
           printf("***************************************************\n");
        }
    }
    else {
        if (M >= N) {
           printf("************************************************\n");
           printf(" - TESTING CGEBRD .. M >= N .. FAILED !\n");
           printf("************************************************\n");
        }
        else {
           printf("************************************************\n");
           printf(" - TESTING CGEBRD .. M < N .. FAILED !\n");
           printf("************************************************\n");
        }
    }

fin:
    if ( A2 != NULL ) free(A2);
    if ( U  != NULL ) free(U);
    if ( VT != NULL ) free(VT);
    free(A1); free(S1); free(S2);free(E2);
    PLASMA_Dealloc_Handle_Tile(&T);

    return 0;
}

/*-------------------------------------------------------------------
 * Check the orthogonality of U VT
 */
static int check_orthogonality(int side, int M, int N, PLASMA_Complex32_t *Q, int LDQ, float eps)
{
    float  done =  1.0;
    float  mdone  = -1.0;
    float  normQ, result;
    int     info_ortho;
    int     minMN = min(M, N);
    float *work = (float *)malloc(minMN*sizeof(float));

    /* Build the idendity matrix */
    PLASMA_Complex32_t *Id = (PLASMA_Complex32_t *) malloc(minMN*minMN*sizeof(PLASMA_Complex32_t));
    LAPACKE_claset_work(LAPACK_COL_MAJOR, 'A', minMN, minMN, 0., 1., Id, minMN);

    /* Perform Id - Q'Q */
    if (M >= N)
        cblas_cherk(CblasColMajor, CblasUpper, CblasConjTrans, N, M, done, Q, LDQ, mdone, Id, N);
    else
        cblas_cherk(CblasColMajor, CblasUpper, CblasNoTrans,   M, N, done, Q, LDQ, mdone, Id, M);

    normQ = LAPACKE_clansy_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), 'U', minMN, Id, minMN, work);

    if (getenv("PLASMA_TESTING_VERBOSE"))
        printf( "||Q||_oo=%f\n", normQ );

    result = normQ / (M * eps);
    if (side == PlasmaLeft)
    {
        printf(" ======================================================\n");
        printf(" ||Id-U'*U||_oo / (M*eps)          : %15.3E \n",  result );
        printf(" ======================================================\n");
    }
    else
    {
        printf(" ======================================================\n");
        printf(" ||Id-VT'*VT||_oo / (M*eps)          : %15.3E \n",  result );
        printf(" ======================================================\n");
    }

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

    free(work); free(Id);

    return info_ortho;
}

/*------------------------------------------------------------
 *  Check the bidiagonal reduction
 */
static int check_reduction(int M, int N, float *S, PLASMA_Complex32_t *A, int LDA,
                           PLASMA_Complex32_t *U, int LDU, PLASMA_Complex32_t *VT, int LDVT, float eps )
{
    PLASMA_Complex32_t zone = 1.0;
    PLASMA_Complex32_t mzone  = -1.0;
    float Anorm, Rnorm, result;
    int info_reduction;
    int i;
    int maxMN = max(M, N);
    int minMN = min(M, N);

    PLASMA_Complex32_t *TEMP     = (PLASMA_Complex32_t *)malloc(maxMN*maxMN*sizeof(PLASMA_Complex32_t));
    PLASMA_Complex32_t *Residual = (PLASMA_Complex32_t *)malloc(maxMN*maxMN*sizeof(PLASMA_Complex32_t));
    float *work = (float *)malloc(maxMN*sizeof(float));

    memset((void*)TEMP, 0, maxMN*maxMN*sizeof(PLASMA_Complex32_t));

    /* Compute TEMP =  U * SIGMA */
    /*
     * U is of size (M,M),
     * SIGMA is of size (M,N) but contains only MINMN singular values.
     * Thus TEMP is of size (M,N)
     */
    LAPACKE_clacpy_work(LAPACK_COL_MAJOR, lapack_const(PlasmaUpperLower), M, M, U, LDU, TEMP, M);
    for (i = 0; i < minMN; i++){
            cblas_csscal(M, S[i], &(TEMP[i*M]),1);
    }
    /* Compute Residual = A - U * SIGMA * VT */
    /*
     * VT is of size (N,N)
     */
    LAPACKE_clacpy_work(LAPACK_COL_MAJOR, lapack_const(PlasmaUpperLower), M, N, A, LDA, Residual, M);
    cblas_cgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, M, N, N, CBLAS_SADDR(mzone), TEMP, M,  VT, LDVT, CBLAS_SADDR(zone), Residual,     M);

    /* Compute the norms */
    Rnorm = LAPACKE_clange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), M, N, Residual, M,   work);
    Anorm = LAPACKE_clange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), M, N, A,       LDA, work);

    if (getenv("PLASMA_TESTING_VERBOSE"))
        printf( "||A||_oo=%f\n||A - U*B*VT||_oo=%e\n", Anorm, Rnorm );

    result = Rnorm / ( Anorm * maxMN * eps);
    printf(" ======================================================\n");
    printf(" ||A-U*SIGMA*VT'||_oo/(||A||_oo.N.eps) : %15.3E \n",  result );
    printf(" ======================================================\n");

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

    free(TEMP);
    free(Residual);
    free(work);

    return info_reduction;
}

/*------------------------------------------------------------
 *  Check the eigenvalues
 */
static int check_solution(int N, float *E1, float *E2, float eps)
{
    int info_solution, i;
    float resid;
    float maxtmp;
    float maxel = fabs( fabs(E1[0]) - fabs(E2[0]) );
    float maxeig = max( fabs(E1[0]), fabs(E2[0]) );
    for (i = 1; i < N; i++){
        resid   = fabs(fabs(E1[i])-fabs(E2[i]));
        maxtmp  = max(fabs(E1[i]), fabs(E2[i]));

        /* Update */
        maxeig = max(maxtmp, maxeig);
        maxel  = max(resid,  maxel );
    }

    maxel = maxel / (maxeig * N * eps);
    printf(" ======================================================\n");
    printf(" | S - singularcomputed | / (|S| * N * eps) : %15.3E \n",  maxel );
    printf(" ======================================================\n");

    if ( isnan(maxel) || isinf(maxel) || (maxel > 100) ) {
        printf("-- The singular values are suspicious ! \n");
        info_solution = 1;
    }
    else{
        printf("-- The singular values are CORRECT ! \n");
        info_solution = 0;
    }
    return info_solution;
}


text:
/*
 * Copyright (c) 2011-2020 The University of Tennessee and The University
 *                         of Tennessee Research Foundation.  All rights
 *                         reserved.
 * Copyright (c) 2013      Inria. All rights reserved.
 *
 * @precisions normal z -> c d s
 *
 */

#include "dplasma.h"
#include "dplasma/types.h"
#include "dplasmaaux.h"
#include "parsec/data_dist/matrix/two_dim_rectangle_cyclic.h"
#include "parsec/data_dist/matrix/vector_two_dim_cyclic.h"
#include "cores/core_blas.h"

#include "zpltmg_chebvand.h"
#include "zpltmg_fiedler.h"
#include "zpltmg_hankel.h"
#include "zpltmg_toeppd.h"

#include <string.h>
#include <lapacke.h>

/**
 *******************************************************************************
 *
 *  Generic case
 *
 *******************************************************************************
 */
struct zpltmg_args_s {
    dplasma_enum_t            mtxtype;
    unsigned long long int seed;
    dplasma_complex64_t     *W;
};
typedef struct zpltmg_args_s zpltmg_args_t;

static int
dplasma_zpltmg_generic_operator( parsec_execution_stream_t *es,
                                 const parsec_tiled_matrix_t *descA,
                                 void *_A,
                                 dplasma_enum_t uplo, int m, int n,
                                 void *op_data )
{
    int tempmm, tempnn, ldam;
    zpltmg_args_t     *args = (zpltmg_args_t*)op_data;
    dplasma_complex64_t *A    = (dplasma_complex64_t*)_A;
    (void)es;
    (void)uplo;

    tempmm = (m == (descA->mt-1)) ? (descA->m - m * descA->mb) : descA->mb;
    tempnn = (n == (descA->nt-1)) ? (descA->n - n * descA->nb) : descA->nb;
    ldam   = BLKLDD( descA, m );

    if ( args->mtxtype == dplasmaMatrixCircul ) {
        return CORE_zpltmg_circul(
            tempmm, tempnn, A, ldam,
            descA->m, m*descA->mb, n*descA->nb, args->W );
    } else {
        return CORE_zpltmg(
            args->mtxtype, tempmm, tempnn, A, ldam,
            descA->m, descA->n, m*descA->mb, n*descA->nb, args->seed );
    }
}

/**
 *******************************************************************************
 *
 * @ingroup dplasma_complex64
 *
 * dplasma_zpltmg_generic - Generic wrapper for cases that are based on the map
 * function. This is the default for many test matrices generation.
 *
 * See parsec_apply() for further information.
 *
 *******************************************************************************
 *
 * @param[in,out] parsec
 *          The parsec context of the application that will run the operation.
 *
 * @param[in] mtxtype
 *          Type of matrix to be generated.
 *
 * @param[in,out] A
 *          Descriptor of the distributed matrix A to generate. Any tiled matrix
 *          descriptor can be used.
 *
 * @param[out] W
 *          Workspace required by some generators.
 *
 * @param[in] seed
 *          The seed used in the random generation.
 *
 *******************************************************************************
 *
 * @return
 *          \retval -i if the ith parameters is incorrect.
 *          \retval 0 on success.
 *
 *******************************************************************************
 *
 * @sa dplasma_cpltmg_genvect
 * @sa dplasma_dpltmg_genvect
 * @sa dplasma_spltmg_genvect
 *
 ******************************************************************************/
static inline int
dplasma_zpltmg_generic( parsec_context_t *parsec,
                        dplasma_enum_t mtxtype,
                        parsec_tiled_matrix_t *A,
                        dplasma_complex64_t *W,
                        unsigned long long int seed)
{
    parsec_taskpool_t *parsec_zpltmg = NULL;
    zpltmg_args_t *params = (zpltmg_args_t*)malloc(sizeof(zpltmg_args_t));

    params->mtxtype = mtxtype;
    params->seed    = seed;
    params->W       = W;

    parsec_zpltmg = parsec_apply_New( dplasmaUpperLower, A, dplasma_zpltmg_generic_operator, params );
    if ( parsec_zpltmg != NULL )
    {
        parsec_context_add_taskpool(parsec, (parsec_taskpool_t*)parsec_zpltmg);
        dplasma_wait_until_completion(parsec);
        parsec_apply_Destruct( parsec_zpltmg );
        return 0;
    }
    return -101;
}

/**
 *******************************************************************************
 *
 * @ingroup dplasma_complex64
 *
 * dplasma_zpltmg_genvect - Generic wrapper for cases that are using two
 * datatypes: the default one, and one describing a vector.
 *
 *******************************************************************************
 *
 * @param[in,out] parsec
 *          The parsec context of the application that will run the operation.
 *
 * @param[in] mtxtype
 *          Type of matrix to be generated.
 *
 * @param[in,out] A
 *          Descriptor of the distributed matrix A to generate. Any tiled matrix
 *          descriptor can be used.
 *
 * @param[in] seed
 *          The seed used in the random generation.
 *
 *******************************************************************************
 *
 * @return
 *          \retval -i if the ith parameters is incorrect.
 *          \retval 0 on success.
 *
 *******************************************************************************
 *
 * @sa dplasma_cpltmg_genvect
 * @sa dplasma_dpltmg_genvect
 * @sa dplasma_spltmg_genvect
 *
 ******************************************************************************/
static inline int
dplasma_zpltmg_genvect( parsec_context_t *parsec,
                        dplasma_enum_t mtxtype,
                        parsec_tiled_matrix_t *A,
                        unsigned long long int seed )
{
    size_t vectorsize = 0;
    parsec_taskpool_t* tp;

    switch( mtxtype ) {
    case dplasmaMatrixChebvand:
        tp = (parsec_taskpool_t*)parsec_zpltmg_chebvand_new( seed,
                                                            A );
        vectorsize = 2 * A->nb * sizeof(dplasma_complex64_t);
        break;

    case dplasmaMatrixFiedler:
        tp = (parsec_taskpool_t*)parsec_zpltmg_fiedler_new( seed,
                                                            A );
        vectorsize = A->mb * sizeof(dplasma_complex64_t);
        break;

    case dplasmaMatrixHankel:
        tp = (parsec_taskpool_t*)parsec_zpltmg_hankel_new( seed,
                                                           A );
        vectorsize = A->mb * sizeof(dplasma_complex64_t);
        break;

    case dplasmaMatrixToeppd:
        tp = (parsec_taskpool_t*)parsec_zpltmg_toeppd_new( seed,
                                                           A );
        vectorsize = 2 * A->mb * sizeof(dplasma_complex64_t);
        break;

    default:
        return -2;
    }

    if (tp != NULL) {
        parsec_zpltmg_hankel_taskpool_t *zpltmg_tp = (parsec_zpltmg_hankel_taskpool_t*)tp;

        /* Default type */
        dplasma_add2arena_tile( &zpltmg_tp->arenas_datatypes[PARSEC_zpltmg_hankel_DEFAULT_ADT_IDX],
                                A->mb*A->nb*sizeof(dplasma_complex64_t),
                                PARSEC_ARENA_ALIGNMENT_SSE,
                                parsec_datatype_double_complex_t, A->mb );

        /* Vector type */
        dplasma_add2arena_tile( &zpltmg_tp->arenas_datatypes[PARSEC_zpltmg_hankel_VECTOR_ADT_IDX],
                                vectorsize,
                                PARSEC_ARENA_ALIGNMENT_SSE,
                                parsec_datatype_double_complex_t, A->mb );

        parsec_context_add_taskpool(parsec, tp);
        dplasma_wait_until_completion(parsec);

        dplasma_matrix_del2arena( &zpltmg_tp->arenas_datatypes[PARSEC_zpltmg_hankel_DEFAULT_ADT_IDX] );
        dplasma_matrix_del2arena( &zpltmg_tp->arenas_datatypes[PARSEC_zpltmg_hankel_VECTOR_ADT_IDX ] );
        parsec_taskpool_free(tp);
        return 0;
    }
    return -101;
}

/**
 *******************************************************************************
 *
 * @ingroup dplasma_complex64
 *
 * dplasma_zpltmg_circul - Generates a Circulant test matrix by tiles.
 *
 *******************************************************************************
 *
 * @param[in,out] parsec
 *          The parsec context of the application that will run the operation.
 *
 * @param[in,out] A
 *          Descriptor of the distributed matrix A to generate. Any tiled matrix
 *          descriptor can be used.
 *
 * @param[in] seed
 *          The seed used in the random generation.
 *
 *******************************************************************************
 *
 * @return
 *          \retval -i if the ith parameters is incorrect.
 *          \retval 0 on success.
 *
 *******************************************************************************
 *
 * @sa dplasma_cpltmg_circul
 * @sa dplasma_dpltmg_circul
 * @sa dplasma_spltmg_circul
 *
 ******************************************************************************/
static inline int
dplasma_zpltmg_circul( parsec_context_t *parsec,
                       parsec_tiled_matrix_t *A,
                       unsigned long long int seed )
{
    int info;
    dplasma_complex64_t *V = (dplasma_complex64_t*) malloc( A->m * sizeof(dplasma_complex64_t) );

    CORE_zplrnt( A->m, 1, V, A->m, A->m, 0, 0, seed );

    info = dplasma_zpltmg_generic(parsec, dplasmaMatrixCircul, A, V, seed);

    free(V);
    return info;
}

/**
 *******************************************************************************
 *
 * @ingroup dplasma_internal
 *
 * dplasma_zpltmg_condex - Generates a Condex test matrix by tiles.
 *
 *******************************************************************************
 *
 * @param[in,out] parsec
 *          The parsec context of the application that will run the operation.
 *
 * @param[in,out] A
 *          Descriptor of the distributed matrix A to generate. Any tiled matrix
 *          descriptor can be used.
 *
 *******************************************************************************
 *
 * @return
 *          \retval -i if the ith parameters is incorrect.
 *          \retval 0 on success.
 *
 *******************************************************************************
 *
 * @sa dplasma_cpltmg_condex
 * @sa dplasma_dpltmg_condex
 * @sa dplasma_spltmg_condex
 *
 ******************************************************************************/
static inline int
dplasma_zpltmg_condex( parsec_context_t *parsec,
                       parsec_tiled_matrix_t *A )
{
    /* gallery('condex', A->m, 4, 100.) */
    dplasma_complex64_t theta = 100.;
    parsec_matrix_block_cyclic_t *twodA = (parsec_matrix_block_cyclic_t *)A;
    parsec_matrix_block_cyclic_t Q;
    parsec_matrix_block_cyclic_init( &Q, PARSEC_MATRIX_COMPLEX_DOUBLE, PARSEC_MATRIX_TILE,
                               A->super.myrank,
                               A->mb, A->nb, A->mb*A->mt, 3, 0, 0, A->m, 3, 1, 1, twodA->grid.krows, twodA->grid.kcols, twodA->grid.ip, twodA->grid.jq );
    Q.mat = parsec_data_allocate((size_t)Q.super.nb_local_tiles *
                                (size_t)Q.super.bsiz *
                                (size_t)parsec_datadist_getsizeoftype(Q.super.mtype));
    parsec_data_collection_set_key((parsec_data_collection_t*)&Q, "Q");

    if (A->super.myrank == 0) {
        dplasma_complex64_t *Qmat;

        Qmat = (dplasma_complex64_t*)(Q.mat);

        /* Initialize the Q matrix */
        CORE_zpltmg_condexq( A->m, A->n, Qmat, Q.super.lm );

        /*
         * Conversion to tile layout
         */
        {
            dplasma_complex64_t *W = (dplasma_complex64_t*) malloc (A->mb * sizeof(dplasma_complex64_t) );
            int *leaders = NULL;
            int i, nleaders;

            /* Get all the cycles leaders and length
             * They are the same for each independent problem (each panel) */
            GKK_getLeaderNbr( Q.super.lmt, A->nb, &nleaders, &leaders );

            /* shift cycles. */
            for(i=0; i<nleaders; i++) {

                /* cycle #i belongs to this thread, so shift it */
                memcpy(W, Qmat + leaders[i*3] * A->mb, A->mb * sizeof(dplasma_complex64_t) );
                CORE_zshiftw(leaders[i*3], leaders[i*3+1], A->mt, A->nb, A->mb, Qmat, W);
            }

            free(leaders); free(W);
        }
    }

    dplasma_zlaset( parsec, dplasmaUpperLower, 0., 1. + theta, A );
    dplasma_zgemm( parsec, dplasmaNoTrans, dplasmaConjTrans,
                   -theta, (parsec_tiled_matrix_t*)&Q,
                           (parsec_tiled_matrix_t*)&Q,
                   1.,     A );

    parsec_data_free(Q.mat);
    parsec_tiled_matrix_destroy((parsec_tiled_matrix_t*)&Q);
    return 0;
}

/**
 *******************************************************************************
 *
 * @ingroup dplasma_complex64
 *
 * dplasma_zpltmg_house - Generates a Householder test matrix by tiles.
 *
 *******************************************************************************
 *
 * @param[in,out] parsec
 *          The parsec context of the application that will run the operation.
 *
 * @param[in,out] A
 *          Descriptor of the distributed matrix A to generate. Any tiled matrix
 *          descriptor can be used.
 *
 * @param[in] seed
 *          The seed used in the random generation.
 *
 *******************************************************************************
 *
 * @return
 *          \retval -i if the ith parameters is incorrect.
 *          \retval 0 on success.
 *
 *******************************************************************************
 *
 * @sa dplasma_cpltmg_house
 * @sa dplasma_dpltmg_house
 * @sa dplasma_spltmg_house
 *
 ******************************************************************************/
static inline int
dplasma_zpltmg_house( parsec_context_t *parsec,
                      parsec_tiled_matrix_t *A,
                      unsigned long long int seed )
{
    /* gallery('house', random, 0 ) */
    parsec_vector_two_dim_cyclic_t V;
    dplasma_complex64_t *Vmat, tau;

    parsec_vector_two_dim_cyclic_init( &V, PARSEC_MATRIX_COMPLEX_DOUBLE, PARSEC_VECTOR_DISTRIB_DIAG,
                                A->super.myrank,
                                A->mb, A->m, 0, A->m, 1, 1 );
    V.mat = parsec_data_allocate((size_t)V.super.nb_local_tiles *
                                (size_t)V.super.bsiz *
                                (size_t)parsec_datadist_getsizeoftype(V.super.mtype));
    parsec_data_collection_set_key((parsec_data_collection_t*)&V, "V");
    Vmat = (dplasma_complex64_t*)(V.mat);

    /* Initialize Householder vector */
    if (A->super.myrank == 0) {
        CORE_zplrnt( A->m, 1, Vmat, A->m, A->m, 0, 0, seed );
        LAPACKE_zlarfg_work( A->m, Vmat, Vmat+1, 1, &tau );
        Vmat[0] = 1.;
    }

#if defined(PARSEC_HAVE_MPI)
    /* If we don't need to broadcast, don't do it, this way we don't require MPI to be initialized */
    if( A->super.nodes > 1 )
        MPI_Bcast( &tau, 1, parsec_datatype_double_complex_t, 0, *(MPI_Comm*)dplasma_pcomm );
#endif

    /* Compute the Householder matrix I - tau v * v' */
    dplasma_zlaset( parsec, dplasmaUpperLower, 0., 1., A);
    dplasma_zgerc( parsec, -tau,
                   (parsec_tiled_matrix_t*)&V,
                   (parsec_tiled_matrix_t*)&V,
                   A );

    parsec_data_free(V.mat);
    parsec_tiled_matrix_destroy((parsec_tiled_matrix_t*)&V);

    return 0;
}

/**
 *******************************************************************************
 *
 * @ingroup dplasma_complex64
 *
 * dplasma_zpltmg - Generates a special test matrix by tiles.
 *
 *******************************************************************************
 *
 * @param[in,out] parsec
 *          The parsec context of the application that will run the operation.
 *
 * @param[in] mtxtype
 *          See PLASMA_zpltmg() for possible values and information on generated
 *          matrices.
 *
 * @param[in,out] A
 *          Descriptor of the distributed matrix A to generate. Any tiled matrix
 *          descriptor can be used.
 *          On exit, the matrix A generated.
 *
 * @param[in] seed
 *          The seed used in the random generation.
 *
 *******************************************************************************
 *
 * @return
 *          \retval -i if the ith parameters is incorrect.
 *          \retval 0 on success.
 *
 *******************************************************************************
 *
 * @sa dplasma_cpltmg
 * @sa dplasma_dpltmg
 * @sa dplasma_spltmg
 *
 ******************************************************************************/
int
dplasma_zpltmg( parsec_context_t *parsec,
                dplasma_enum_t mtxtype,
                parsec_tiled_matrix_t *A,
                unsigned long long int seed)
{

    switch( mtxtype ) {
    case dplasmaMatrixCircul:
        return dplasma_zpltmg_circul(parsec, A, seed);
        break;

    case dplasmaMatrixChebvand:
        return dplasma_zpltmg_genvect(parsec, mtxtype,
                                      A, seed);
        break;

    case dplasmaMatrixCondex:
        return dplasma_zpltmg_condex(parsec, A);
        break;

    case dplasmaMatrixFiedler:
        return dplasma_zpltmg_genvect(parsec, mtxtype,
                                      A, seed);
        break;

    case dplasmaMatrixHankel:
        return dplasma_zpltmg_genvect(parsec, mtxtype,
                                      A, seed);
        break;

    case dplasmaMatrixHouse:
        return dplasma_zpltmg_house(parsec, A, seed);
        break;

    case dplasmaMatrixToeppd:
        return dplasma_zpltmg_genvect(parsec, mtxtype,
                                      A, seed);
        break;

    case dplasmaMatrixCauchy:
    case dplasmaMatrixCompan:
    case dplasmaMatrixDemmel:
    case dplasmaMatrixDorr:
    case dplasmaMatrixFoster:
    case dplasmaMatrixHadamard:
    case dplasmaMatrixHilb:
    case dplasmaMatrixInvhess:
    case dplasmaMatrixKms:
    case dplasmaMatrixLangou:
    case dplasmaMatrixLehmer:
    case dplasmaMatrixLotkin:
    case dplasmaMatrixMinij:
    case dplasmaMatrixMoler:
    case dplasmaMatrixOrthog:
    case dplasmaMatrixParter:
    case dplasmaMatrixRandom:
    case dplasmaMatrixRiemann:
    case dplasmaMatrixRis:
    case dplasmaMatrixWilkinson:
    case dplasmaMatrixWright:
        return dplasma_zpltmg_generic(parsec, mtxtype, A, NULL, seed);
        break;
    default:
        return -2;
    }

    return 0;
}

max_stars_repo_name: therault/dplasma

max_issues_repo_name: therault/dplasma

max_forks_repo_name: therault/dplasma