text:
#ifndef KABSCH_RMSD
#define KABSCH_RMSD

#include <cblas.h>
#include <lapacke.h>

#include <iterator>

#include <iostream>
#include <fstream>

#include <tuple>

namespace kabsch
{

/**

	Calculates the RMSD

	@param P - coordinates of molecule P
	@param Q - coordinates of molecule Q
	@param N - number of atoms

	@return double - the RMSD

**/
template <class M>
double rmsd(
	M *P,
	M *Q,
	const unsigned int N)
{
    double rmsd {0.0};
    const unsigned int D {3};
	const unsigned int size = N*D;

    for(unsigned int i = 0; i < size; ++i) {
        rmsd += (P[i] - Q[i])*(P[i] - Q[i]);
    }

    return sqrt(rmsd/N);
}

template <class M>
M* centroid(
	M *coordinates,
	unsigned int n_atoms)
{
    double x {0};
    double y {0};
    double z {0};
    // unsigned int size = sizeof(coordinates);
    // unsigned int n_atoms = size / 3;

	const unsigned int size = n_atoms*3;

    unsigned int i = 0;
    while(i<size)
    {
        x += coordinates[i++];
        y += coordinates[i++];
        z += coordinates[i++];
    }

    x /= n_atoms;
    y /= n_atoms;
    z /= n_atoms;

	M *centroid = new M[3];
	centroid[0] = x;
	centroid[1] = y;
	centroid[2] = z;

    return centroid;
}

template <class Matrix>
Matrix* multiply(Matrix *A, Matrix *B,
    const int M,
    const int N,
    const int K)
{
    double one = 1.0;

	// TODO Matrix C[M*N] = {0}
	// Compilers doesnt like the following syntax
	Matrix *C = new Matrix[M*N] {0};

    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
        M, N, K, one,
        A, K,
        B, N, one,
        C, N);

    return C;
}

template <class Matrix>
Matrix* transpose_multiply(
	Matrix *A,
	Matrix *B,
    const int M,
    const int N,
    const int K)
{
    double one = 1.0;

    Matrix *C = new Matrix[M*N] {0};

    cblas_dgemm(CblasRowMajor, CblasTrans, CblasNoTrans,
        M, N, K, one,
        A, M,
        B, N, one,
        C, N);

    return C;
}

template <class Matrix>
std::tuple<Matrix*, Matrix*, Matrix*> matrix_svd(
	Matrix *A,
	int rows,
	int cols)
{
    // lapack_int LAPACKE_dgesvd( int matrix_layout, char jobu, char jobvt,
    //     lapack_int m, lapack_int n,
    //     double* a, lapack_int lda,
    //     double* s, double* u, lapack_int ldu,
    //     double* vt, lapack_int ldvt,
    //     double* superb );

    Matrix *U = new Matrix[cols*rows] {0};
    Matrix *S = new Matrix[rows] {0};
    Matrix *VT = new Matrix[cols*rows] {0};
    Matrix *superb = new Matrix[cols*rows] {0};
    int info;

    info = LAPACKE_dgesvd(LAPACK_ROW_MAJOR, 'A', 'A',
        rows, cols,
        A, rows,
        S,
        U, rows,
        VT, rows,
        superb);

    if(info > 0)
    {
        // TODO What if failed
    }

	delete[] superb;

    return std::make_tuple(U, S, VT);
}

template <class M>
double determinant3x3(M A)
{
    // determinant of a square 3x3 matrix
    double det = A[0]*A[4]*A[8]
        +A[1]*A[5]*A[6]
        +A[2]*A[3]*A[7]
        -A[2]*A[4]*A[6]
        -A[1]*A[3]*A[8]
        -A[0]*A[5]*A[7];

    return det;
}

template <class M, class T>
M* kabsch(
	M *P,
	M *Q,
	const T n_atoms)
{
	// M *U = new M[3*3] {0};
	// U[0] = 1.0;
	// U[4] = 1.0;
	// U[8] = 1.0;
	// return U;

    M *U;
    M *S;
    M *V;

    M *C = transpose_multiply(P, Q, 3, 3, n_atoms);

	std::tie(U, S, V) = matrix_svd(C, 3, 3);

    // Getting the sign of the det(U)*(V) to decide whether we need to correct
    // our rotation matrix to ensure a right-handed coordinate system.
    if(determinant3x3(U)*determinant3x3(V) < 0.0)
    {
		U[3*0+2] = -U[3*0+2];
		U[3*2+2] = -U[3*1+2];
		U[3*1+2] = -U[3*2+2];
    }

    M *rotation = multiply(U, V, 3, 3, 3);

	delete[] C;
	delete[] U;
	delete[] S;
	delete[] V;

    return rotation;
}

template <class M, class T>
M* kabsch_rotate(
	M *P,
	M *Q,
	T n_atoms)
{
    M *U = kabsch(P, Q, n_atoms);
    M *product = multiply(P, U, n_atoms, 3, 3);
	delete[] U;
    return product;
}

template <class M, class T>
double kabsch_rmsd(
	M *P,
	M *Q,
	T n_atoms)
{
    M *P_rotated = kabsch_rotate(P, Q, n_atoms);
	double rmsdval = rmsd(P_rotated, Q, n_atoms);
	delete[] P_rotated;
	return rmsdval;
}

} // namespace rmsd

#endif

ext: h

max_stars_repo_name: chemspacelab/clockwork

max_issues_repo_name: chemspacelab/clockwork

max_forks_repo_name: chemspacelab/clockwork


text:
#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_math.h>

#include "allvars.h"
#include "proto.h"

#ifdef RADTRANSFER
#define sigma 6.3e-18/All.UnitLength_in_cm/All.UnitLength_in_cm	// in code units
#endif

#ifdef CS_MODEL
#include "cs_metals.h"
#endif

/*! \file accel.c
 *  \brief driver routines to carry out force computation
 */

/*! This routine computes the accelerations for all active particles.  First, the gravitational forces are
 * computed. This also reconstructs the tree, if needed, otherwise the drift/kick operations have updated the
 * tree to make it fullu usable at the current time.
 *
 * If gas particles are presented, the `interior' of the local domain is determined. This region is guaranteed
 * to contain only particles local to the processor. This information will be used to reduce communication in
 * the hydro part.  The density for active SPH particles is computed next. If the number of neighbours should
 * be outside the allowed bounds, it will be readjusted by the function ensure_neighbours(), and for those
 * particle, the densities are recomputed accordingly. Finally, the hydrodynamical forces are added.
 */
void compute_accelerations(int mode)
{
#ifdef RADTRANSFER
  double timeeach = 0, timeall = 0, tstart = 0, tend = 0;
#endif
#if defined(BUBBLES) || defined(MULTI_BUBBLES)
  double hubble_a;
#endif

  if(ThisTask == 0)
    {
      printf("Start force computation...\n");
      fflush(stdout);
    }

#ifdef REIONIZATION
  heating();
#endif

  CPU_Step[CPU_MISC] += measure_time();

#ifdef PMGRID
  if(All.PM_Ti_endstep == All.Ti_Current)
    {
      long_range_force();

      CPU_Step[CPU_MESH] += measure_time();

    }
#endif

#ifndef ONLY_PM

  gravity_tree();		/* computes gravity accel. */

  if(All.TypeOfOpeningCriterion == 1 && All.Ti_Current == 0)
    gravity_tree();		/* For the first timestep, we redo it
				 * to allow usage of relative opening
				 * criterion for consistent accuracy.
				 */
#endif

  if(All.Ti_Current == 0 && RestartFlag == 0 && header.flag_ic_info == FLAG_SECOND_ORDER_ICS)
    second_order_ics();		/* produces the actual ICs from the special second order IC file */

#ifdef FORCETEST
  gravity_forcetest();
#endif

  if(All.TotN_gas > 0)
    {
      /***** density *****/
      if(ThisTask == 0)
	{
	  printf("Start density computation...\n");
	  fflush(stdout);
	}

#ifdef CS_MODEL
      CPU_Step[CPU_MISC] += measure_time();

#if defined(CS_SNI) || defined(CS_SNII)
      cs_flag_SN_starparticles();	/* mark SNI star particles */
#endif
      cs_copy_densities();
      cs_find_low_density_tail();

      CPU_Step[CPU_CSMISC] += measure_time();
#endif

#ifndef VORONOI
      density();		/* computes density, and pressure */
#else
      voronoi_mesh();
      voronoi_density();
#endif

#if (defined(SMOOTH_PHI) || defined(SMOOTH_ROTB) || defined(BSMOOTH))
      smoothed_values();
#endif

#if defined(SNIA_HEATING)
      snIa_heating();
#endif

#if defined(CS_MODEL) && defined(CS_FEEDBACK)

      CPU_Step[CPU_CSMISC] += measure_time();

      cs_find_hot_neighbours();

      cs_promotion();
      cs_copy_densities();
      CPU_Step[CPU_CSMISC] += measure_time();
      density();		/* recalculate densities again */
      CPU_Step[CPU_CSMISC] += measure_time();
#endif

      /***** update smoothing lengths in tree *****/
#ifndef VORONOI
      force_update_hmax();
#endif

      /***** hydro forces *****/
      if(ThisTask == 0)
	{
	  printf("Start hydro-force computation...\n");
	  fflush(stdout);
	}

#ifndef VORONOI
      hydro_force();		/* adds hydrodynamical accelerations  and computes du/dt  */
#else
      voronoi_hydro_force();
#endif

#ifdef CONDUCTION
      if(All.Conduction_Ti_endstep == All.Ti_Current)
	conduction();
#endif

#ifdef CR_DIFFUSION
      if(All.CR_Diffusion_Ti_endstep == All.Ti_Current)
	cosmic_ray_diffusion();
#endif

#ifdef RADTRANSFER
      if(Flag_FullStep)		/* only do it for full timesteps */
	{
	  All.Radiation_Ti_endstep = All.Ti_Current;

#ifndef NOGRAVITY
	  /***** compute eddington tensor *****/
	  if(ThisTask == 0)
	    {
	      printf("Start Eddington tensor computation...\n");
	      fflush(stdout);
	    }

	  eddington();

	  if(ThisTask == 0)
	    {
	      printf("done Eddington tensor! \n");
	      fflush(stdout);
	    }

	  star_density();
#else
	  if(All.Time == All.TimeBegin)
	    {
	      if(ThisTask == 0)
		{
		  printf("Start Eddington tensor computation...\n");
		  fflush(stdout);
		}

	      eddington();

	      if(ThisTask == 0)
		{
		  printf("done Eddington tensor! \n");
		  fflush(stdout);
		}
	    }
#endif

	  star_density();

	  /***** set simple initial conditions *****/
	  if(All.Time == All.TimeBegin)
	    {
	      if(ThisTask == 0)
		{
		  printf("Setting simple inits...\n");
		  fflush(stdout);
		}

	      radtransfer_set_simple_inits();

	      if(ThisTask == 0)
		{
		  printf("done with simple inits! \n");
		  fflush(stdout);
		}
	    }

	  /***** evolve the transport of radiation *****/
	  if(ThisTask == 0)
	    {
	      printf("start radtransfer...\n");
	      fflush(stdout);
	    }

	  tstart = second();

	  radtransfer();

	  tend = second();
	  timeeach = timediff(tstart, tend);
	  MPI_Allreduce(&timeeach, &timeall, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);

	  if(ThisTask == 0)
	    {
	      printf("time consumed is %g \n", timeall);
	      printf("done with radtransfer! \n");
	      fflush(stdout);
	    }

	  All.Radiation_Ti_begstep = All.Radiation_Ti_endstep;
	}
#endif

#ifdef MHM
      /***** kinetic feedback *****/
      kinetic_feedback_mhm();
#endif

#ifdef BLACK_HOLES
      /***** black hole accretion and feedback *****/
      blackhole_accretion();
#ifdef FOF
      /* this will find new black hole seed halos */
      if(All.Time >= All.TimeNextBlackHoleCheck)
	{
	  fof_fof(-1);

	  if(All.ComovingIntegrationOn)
	    All.TimeNextBlackHoleCheck *= All.TimeBetBlackHoleSearch;
	  else
	    All.TimeNextBlackHoleCheck += All.TimeBetBlackHoleSearch;
	}
#endif
#endif

#ifdef COOLING	/**** radiative cooling and star formation *****/

#ifdef CS_MODEL
      cs_cooling_and_starformation();
#else

#ifdef SFR
      cooling_and_starformation();
#else
      cooling_only();
#endif

#endif
      CPU_Step[CPU_COOLINGSFR] += measure_time();
#endif /*ends COOLING */

#if defined(CS_MODEL) && defined(CS_ENRICH)
#ifndef CS_FEEDBACK
      Flag_phase = 0;		/* no destinction between phases */

      cs_update_weights();
      CPU_Step[CPU_WEIGHTS_HOT] += measure_time();
      cs_enrichment();
      CPU_Step[CPU_ENRICH_HOT] += measure_time();
#else

      Flag_phase = 1;		/* COLD phase Flag_phase  = 1 */

      cs_update_weights();
      CPU_Step[CPU_WEIGHTS_HOT] += measure_time();
      cs_enrichment();
      CPU_Step[CPU_ENRICH_HOT] += measure_time();

      Flag_phase = 2;		/* HOT phase Flag_phase = 2 */

      cs_update_weights();
      CPU_Step[CPU_WEIGHTS_COLD] += measure_time();
      cs_enrichment();
      CPU_Step[CPU_ENRICH_COLD] += measure_time();

      Flag_phase = 0;
#endif
#endif

#ifdef CS_TESTS
      cs_energy_test();
#endif

#ifndef BH_BUBBLES
#ifdef BUBBLES
      /**** bubble feedback *****/
      if(All.Time >= All.TimeOfNextBubble)
	{
#ifdef FOFs
	  fof_fof(-1);
	  bubble();
#else
	  bubble();
#endif
	  if(All.ComovingIntegrationOn)
	    {
	      hubble_a = hubble_function(All.Time);
	      All.TimeOfNextBubble *= (1.0 + All.BubbleTimeInterval * hubble_a);
	    }
	  else
	    All.TimeOfNextBubble += All.BubbleTimeInterval / All.UnitTime_in_Megayears;

	  if(ThisTask == 0)
	    printf("Time of the bubble generation: %g\n", 1. / All.TimeOfNextBubble - 1.);
	}
#endif
#endif

#if defined(MULTI_BUBBLES) && defined(FOF)
      if(All.Time >= All.TimeOfNextBubble)
	{
	  fof_fof(-1);

	  if(All.ComovingIntegrationOn)
	    {
	      hubble_a = hubble_func(All.Time);
	      All.TimeOfNextBubble *= (1.0 + All.BubbleTimeInterval * hubble_a);
	    }
	  else
	    All.TimeOfNextBubble += All.BubbleTimeInterval / All.UnitTime_in_Megayears;

	  if(ThisTask == 0)
	    printf("Time of the bubble generation: %g\n", 1. / All.TimeOfNextBubble - 1.);
	}
#endif

    }

  if(ThisTask == 0)
    {
      printf("force computation done.\n");
      fflush(stdout);
    }
}

ext: c

max_stars_repo_name: egpbos/egp

max_issues_repo_name: egpbos/egp

max_forks_repo_name: egpbos/egp


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

/* falsepos.c -- falsepos root finding algorithm

   The false position algorithm uses bracketing by linear interpolation.

   If a linear interpolation step would decrease the size of the
   bracket by less than a bisection step would then the algorithm
   takes a bisection step instead.

   The last linear interpolation estimate of the root is used. If a
   bisection step causes it to fall outside the brackets then it is
   replaced by the bisection estimate (x_upper + x_lower)/2.

*/

#include <config.h>

#include <stddef.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <float.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_roots.h>

#include "roots.h"

typedef struct
  {
    double f_lower, f_upper;
  }
falsepos_state_t;

static int falsepos_init (void * vstate, gsl_function * f, double * root, double x_lower, double x_upper);
static int falsepos_iterate (void * vstate, gsl_function * f, double * root, double * x_lower, double * x_upper);

static int
falsepos_init (void * vstate, gsl_function * f, double * root, double x_lower, double x_upper)
{
  falsepos_state_t * state = (falsepos_state_t *) vstate;

  double f_lower, f_upper ;

  *root = 0.5 * (x_lower + x_upper);

  SAFE_FUNC_CALL (f, x_lower, &f_lower);
  SAFE_FUNC_CALL (f, x_upper, &f_upper);

  state->f_lower = f_lower;
  state->f_upper = f_upper;

  if ((f_lower < 0.0 && f_upper < 0.0) || (f_lower > 0.0 && f_upper > 0.0))
    {
      GSL_ERROR ("endpoints do not straddle y=0", GSL_EINVAL);
    }

  return GSL_SUCCESS;

}

static int
falsepos_iterate (void * vstate, gsl_function * f, double * root, double * x_lower, double * x_upper)
{
  falsepos_state_t * state = (falsepos_state_t *) vstate;

  double x_linear, f_linear;
  double x_bisect, f_bisect;

  double x_left = *x_lower ;
  double x_right = *x_upper ;

  double f_lower = state->f_lower;
  double f_upper = state->f_upper;

  double w ;

  if (f_lower == 0.0)
    {
      *root = x_left ;
      *x_upper = x_left;
      return GSL_SUCCESS;
    }

  if (f_upper == 0.0)
    {
      *root = x_right ;
      *x_lower = x_right;
      return GSL_SUCCESS;
    }

  /* Draw a line between f(*lower_bound) and f(*upper_bound) and
     note where it crosses the X axis; that's where we will
     split the interval. */

  x_linear = x_right - (f_upper * (x_left - x_right) / (f_lower - f_upper));

  SAFE_FUNC_CALL (f, x_linear, &f_linear);

  if (f_linear == 0.0)
    {
      *root = x_linear;
      *x_lower = x_linear;
      *x_upper = x_linear;
      return GSL_SUCCESS;
    }

  /* Discard the half of the interval which doesn't contain the root. */

  if ((f_lower > 0.0 && f_linear < 0.0) || (f_lower < 0.0 && f_linear > 0.0))
    {
      *root = x_linear ;
      *x_upper = x_linear;
      state->f_upper = f_linear;
      w = x_linear - x_left ;
    }
  else
    {
      *root = x_linear ;
      *x_lower = x_linear;
      state->f_lower = f_linear;
      w = x_right - x_linear;
    }

  if (w < 0.5 * (x_right - x_left))
    {
      return GSL_SUCCESS ;
    }

  x_bisect = 0.5 * (x_left + x_right);

  SAFE_FUNC_CALL (f, x_bisect, &f_bisect);

  if ((f_lower > 0.0 && f_bisect < 0.0) || (f_lower < 0.0 && f_bisect > 0.0))
    {
      *x_upper = x_bisect;
      state->f_upper = f_bisect;
      if (*root > x_bisect)
        *root = 0.5 * (x_left + x_bisect) ;
    }
  else
    {
      *x_lower = x_bisect;
      state->f_lower = f_bisect;
      if (*root < x_bisect)
        *root = 0.5 * (x_bisect + x_right) ;
    }

  return GSL_SUCCESS;
}

static const gsl_root_fsolver_type falsepos_type =
{"falsepos",                            /* name */
 sizeof (falsepos_state_t),
 &falsepos_init,
 &falsepos_iterate};

const gsl_root_fsolver_type  * gsl_root_fsolver_falsepos = &falsepos_type;

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
// Copyright (c) Microsoft Corporation. All rights reserved.
// Licensed under the MIT License.

#pragma once

#include "core/common/inlined_containers_fwd.h"
#include "core/framework/kernel_registry.h"
#include "core/graph/graph_viewer.h"
#include <gsl/gsl>

namespace onnxruntime {

/**
  Returns a list of nodes that are preferred on CPU.
  They are commonly shape-related computation subgraphs.
  @param graph Graph viewer
  @param provider_type The target execution provider type
  @param kernel_registries Kernel registries for the target EP
  @param tentative_nodes Nodes that are tentative to be placed on on target EP
  */
InlinedHashSet<NodeIndex> GetCpuPreferredNodes(const GraphViewer& graph,
                                               const std::string& provider_type,
                                               gsl::span<const KernelRegistry* const> kernel_registries,
                                               gsl::span<const NodeIndex> tentative_nodes);

}  // namespace onnxruntime

max_issues_repo_name: Grimig/onnxruntime

max_forks_repo_name: Grimig/onnxruntime


text:
/**
  * Copyright 2017 José Manuel Abuín Mosquera <josemanuel.abuin@usc.es>
  *
  * This file is part of Matrix Market Suite.
  *
  * Matrix Market Suite 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.
  *
  * Matrix Market Suite 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 Matrix Market Suite. If not, see <http://www.gnu.org/licenses/>.
  */

#include <cblas.h>
#include "JacobiSolver.h"

int JacobiSolver(unsigned long *II, unsigned long *J, double *A, unsigned long M, unsigned long N, unsigned long long nz, double *b, unsigned long M_Vector, unsigned long N_Vector, unsigned long long nz_vector,double *x2, int iterationNumber) {

	//Jacobi Method as shown in the example from https://en.wikipedia.org/wiki/Jacobi_method

	unsigned long 	k 	= 0;
	int 		stop 	= 0;

	unsigned int	i	= 0;

	double		result	= 0.0;

	//Initial solution
	double *x1=(double *) calloc(nz_vector,sizeof(double));
	//double *x2=(double *) calloc(nz_vector,sizeof(double));

	double *res=(double *) calloc(nz_vector,sizeof(double));

	double *LU=(double *) calloc(nz,sizeof(double));
	double *Dinv=(double *) calloc(nz,sizeof(double));

	double *T=(double *) calloc(nz,sizeof(double));
	double *C=(double *) calloc(nz_vector,sizeof(double));

	if(!isDiagonallyDominant(A,M,N,nz)){
		fprintf(stderr, "[%s] The matrix is not diagonally dominant\n",__func__);
		//return 0;
	}

	getLUValues(LU, A,M,N,nz);

	getDInvValues(Dinv, A,M,N,nz);

	for(i = 0; i< N; i++){
		x1[i] = 1.0;
	}

	//T=-D^{-1}(L+U)
	/*
	void cblas_dgemm (const CBLAS_LAYOUT Layout, const CBLAS_TRANSPOSE transa, const CBLAS_TRANSPOSE transb, const MKL_INT m, const MKL_INT n, const MKL_INT k, const double alpha, const double *a, const MKL_INT lda, const double *b, const MKL_INT ldb, const double beta, double *c, const MKL_INT ldc);
	*/
	cblas_dgemm(CblasRowMajor,CblasNoTrans,CblasNoTrans, M, N, N, -1.0, Dinv, N, LU, N, 0.0, T, N);
	//writeDenseCoordinateMatrixRowLine("stdout", T,M,N,nz);

	//C=-D^{-1}b
	cblas_dgemv(CblasRowMajor, CblasNoTrans, M,N , 1.0, Dinv, N, b, 1, 0.0, C, 1);

	//writeDenseVector("stdout", C,nz_vector,1, nz_vector);
	unsigned long maxIterations = M*2;

	if(iterationNumber != 0 ){
		maxIterations = iterationNumber;
	}

	while(!stop) {

		// x^{(1)}= Tx^{(0)}+C
		//x2 = T*x1
		cblas_dgemv(CblasRowMajor, CblasNoTrans, M,N , 1.0, T, N, x1, 1, 0.0, x2, 1);

		//x2 = x2+C
		cblas_daxpy(nz_vector,1.0,C, 1, x2, 1);

		//res = A*x - b
		memcpy(res, b, nz_vector*sizeof(double));
		cblas_dgemv(CblasRowMajor, CblasNoTrans, M,N , 1.0, A, N, x2, 1, -1.0, res, 1);

		result = vectorSumElements(res,N);
		if((fabs(result)<=EPSILON)||(k == maxIterations)){
			//fprintf(stderr,"Sum vector res is %lg\n",result);
			stop = 1;
		}

		memcpy(x1, x2, nz_vector*sizeof(double));
		//writeDenseVector("stdout", x1,nz_vector,1, nz_vector);
		k++;
	}

	//memcpy(b, x2, N*sizeof(double));

	free(x1);
	//free(x2);
	free(res);

	fprintf(stderr, "[%s] Number of iterations %lu\n",__func__,k);

	return 1;

}


text:
#pragma once

#include <string>
#include <vector>

#include <gsl/gsl_complex.h>
#include <MathEngine.h>

extern MathEngine engine;

bool printDifference(const std::string& input, const expression expr, const double output, const double expected);
void requireIsEqual(const std::string& input, const double expected);

bool printDifference(const std::string& input, const expression expr, const std::string& output, const std::string& expected);
void requireIsEqual(const std::string& input, const std::string& expected, bool evaluate = false);

bool printDifference(const std::string& input, const expression expr, const expression output, const std::string& expected);
void requireExprIsEqual(const std::string& input, const std::string& expected);

void requireErrorIsEqual(const std::string& input, const std::string& expected);

max_stars_repo_name: antoniojkim/CalcPlusPlus

max_issues_repo_name: antoniojkim/CalcPlusPlus

max_forks_repo_name: antoniojkim/CalcPlusPlus


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

#include <config.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_permute_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_complex_math.h>

#include <gsl/gsl_linalg.h>

/* Factorise a general N x N complex matrix A into,
 *
 *   P A = L U
 *
 * where P is a permutation matrix, L is unit lower triangular and U
 * is upper triangular.
 *
 * L is stored in the strict lower triangular part of the input
 * matrix. The diagonal elements of L are unity and are not stored.
 *
 * U is stored in the diagonal and upper triangular part of the
 * input matrix.
 *
 * P is stored in the permutation p. Column j of P is column k of the
 * identity matrix, where k = permutation->data[j]
 *
 * signum gives the sign of the permutation, (-1)^n, where n is the
 * number of interchanges in the permutation.
 *
 * See Golub & Van Loan, Matrix Computations, Algorithm 3.4.1 (Gauss
 * Elimination with Partial Pivoting).
 */

int
gsl_linalg_complex_LU_decomp (gsl_matrix_complex * A, gsl_permutation * p, int *signum)
{
  if (A->size1 != A->size2)
    {
      GSL_ERROR ("LU decomposition requires square matrix", GSL_ENOTSQR);
    }
  else if (p->size != A->size1)
    {
      GSL_ERROR ("permutation length must match matrix size", GSL_EBADLEN);
    }
  else
    {
      const size_t N = A->size1;
      size_t i, j, k;

      *signum = 1;
      gsl_permutation_init (p);

      for (j = 0; j < N - 1; j++)
        {
          /* Find maximum in the j-th column */

          gsl_complex ajj = gsl_matrix_complex_get (A, j, j);
          double max = gsl_complex_abs (ajj);
          size_t i_pivot = j;

          for (i = j + 1; i < N; i++)
            {
              gsl_complex aij = gsl_matrix_complex_get (A, i, j);
              double ai = gsl_complex_abs (aij);

              if (ai > max)
                {
                  max = ai;
                  i_pivot = i;
                }
            }

          if (i_pivot != j)
            {
              gsl_matrix_complex_swap_rows (A, j, i_pivot);
              gsl_permutation_swap (p, j, i_pivot);
              *signum = -(*signum);
            }

          ajj = gsl_matrix_complex_get (A, j, j);

          if (!(GSL_REAL(ajj) == 0.0 && GSL_IMAG(ajj) == 0.0))
            {
              for (i = j + 1; i < N; i++)
                {
                  gsl_complex aij_orig = gsl_matrix_complex_get (A, i, j);
                  gsl_complex aij = gsl_complex_div (aij_orig, ajj);
                  gsl_matrix_complex_set (A, i, j, aij);

                  for (k = j + 1; k < N; k++)
                    {
                      gsl_complex aik = gsl_matrix_complex_get (A, i, k);
                      gsl_complex ajk = gsl_matrix_complex_get (A, j, k);

                      /* aik = aik - aij * ajk */

                      gsl_complex aijajk = gsl_complex_mul (aij, ajk);
                      gsl_complex aik_new = gsl_complex_sub (aik, aijajk);

                      gsl_matrix_complex_set (A, i, k, aik_new);
                    }
                }
            }
        }

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_complex_LU_solve (const gsl_matrix_complex * LU, const gsl_permutation * p, const gsl_vector_complex * b, gsl_vector_complex * x)
{
  if (LU->size1 != LU->size2)
    {
      GSL_ERROR ("LU matrix must be square", GSL_ENOTSQR);
    }
  else if (LU->size1 != p->size)
    {
      GSL_ERROR ("permutation length must match matrix size", GSL_EBADLEN);
    }
  else if (LU->size1 != b->size)
    {
      GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
    }
  else if (LU->size2 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      /* Copy x <- b */

      gsl_vector_complex_memcpy (x, b);

      /* Solve for x */

      gsl_linalg_complex_LU_svx (LU, p, x);

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_complex_LU_svx (const gsl_matrix_complex * LU, const gsl_permutation * p, gsl_vector_complex * x)
{
  if (LU->size1 != LU->size2)
    {
      GSL_ERROR ("LU matrix must be square", GSL_ENOTSQR);
    }
  else if (LU->size1 != p->size)
    {
      GSL_ERROR ("permutation length must match matrix size", GSL_EBADLEN);
    }
  else if (LU->size1 != x->size)
    {
      GSL_ERROR ("matrix size must match solution/rhs size", GSL_EBADLEN);
    }
  else
    {
      /* Apply permutation to RHS */

      gsl_permute_vector_complex (p, x);

      /* Solve for c using forward-substitution, L c = P b */

      gsl_blas_ztrsv (CblasLower, CblasNoTrans, CblasUnit, LU, x);

      /* Perform back-substitution, U x = c */

      gsl_blas_ztrsv (CblasUpper, CblasNoTrans, CblasNonUnit, LU, x);

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_complex_LU_refine (const gsl_matrix_complex * A, const gsl_matrix_complex * LU, const gsl_permutation * p, const gsl_vector_complex * b, gsl_vector_complex * x, gsl_vector_complex * residual)
{
  if (A->size1 != A->size2)
    {
      GSL_ERROR ("matrix a must be square", GSL_ENOTSQR);
    }
  if (LU->size1 != LU->size2)
    {
      GSL_ERROR ("LU matrix must be square", GSL_ENOTSQR);
    }
  else if (A->size1 != LU->size2)
    {
      GSL_ERROR ("LU matrix must be decomposition of a", GSL_ENOTSQR);
    }
  else if (LU->size1 != p->size)
    {
      GSL_ERROR ("permutation length must match matrix size", GSL_EBADLEN);
    }
  else if (LU->size1 != b->size)
    {
      GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
    }
  else if (LU->size1 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      /* Compute residual, residual = (A * x  - b) */

      gsl_vector_complex_memcpy (residual, b);

      {
        gsl_complex one = GSL_COMPLEX_ONE;
        gsl_complex negone = GSL_COMPLEX_NEGONE;
        gsl_blas_zgemv (CblasNoTrans, one, A, x, negone, residual);
      }

      /* Find correction, delta = - (A^-1) * residual, and apply it */

      gsl_linalg_complex_LU_svx (LU, p, residual);

      {
        gsl_complex negone= GSL_COMPLEX_NEGONE;
        gsl_blas_zaxpy (negone, residual, x);
      }

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_complex_LU_invert (const gsl_matrix_complex * LU, const gsl_permutation * p, gsl_matrix_complex * inverse)
{
  size_t i, n = LU->size1;

  int status = GSL_SUCCESS;

  gsl_matrix_complex_set_identity (inverse);

  for (i = 0; i < n; i++)
    {
      gsl_vector_complex_view c = gsl_matrix_complex_column (inverse, i);
      int status_i = gsl_linalg_complex_LU_svx (LU, p, &(c.vector));

      if (status_i)
        status = status_i;
    }

  return status;
}

gsl_complex
gsl_linalg_complex_LU_det (gsl_matrix_complex * LU, int signum)
{
  size_t i, n = LU->size1;

  gsl_complex det = gsl_complex_rect((double) signum, 0.0);

  for (i = 0; i < n; i++)
    {
      gsl_complex zi = gsl_matrix_complex_get (LU, i, i);
      det = gsl_complex_mul (det, zi);
    }

  return det;
}

double
gsl_linalg_complex_LU_lndet (gsl_matrix_complex * LU)
{
  size_t i, n = LU->size1;

  double lndet = 0.0;

  for (i = 0; i < n; i++)
    {
      gsl_complex z = gsl_matrix_complex_get (LU, i, i);
      lndet += log (gsl_complex_abs (z));
    }

  return lndet;
}

gsl_complex
gsl_linalg_complex_LU_sgndet (gsl_matrix_complex * LU, int signum)
{
  size_t i, n = LU->size1;

  gsl_complex phase = gsl_complex_rect((double) signum, 0.0);

  for (i = 0; i < n; i++)
    {
      gsl_complex z = gsl_matrix_complex_get (LU, i, i);

      double r = gsl_complex_abs(z);

      if (r == 0)
        {
          phase = gsl_complex_rect(0.0, 0.0);
          break;
        }
      else
        {
          z = gsl_complex_div_real(z, r);
          phase = gsl_complex_mul(phase, z);
        }
    }

  return phase;
}