text:
/*
 * This module must not depend on any other as it is needed during the startup
 * of pygsl.
 */
#include <Python.h>
#include <gsl/gsl_errno.h>
#include <pygsl/errorno.h>

static int
add_errno(PyObject * dict, int num, char * name)
{
     PyObject * item;

     item = PyInt_FromLong(num);
     if(item == NULL){
	  fprintf(stderr, "Failed to generate PyInt with value %d for errno %s\n",
		  num, name);
	  return -1;
     }
     if(PyDict_SetItemString(dict, name, item) == -1){
	  fprintf(stderr, "Failed to add PyInt %p with value %d to dict %p for errno %s\n",
		  (void *) item, num, (void *) dict, name);
	  return -1;
     }
     return 0;
}
#define ADD_ERRNO(ERRNO, ERRNOSTR) if(add_errno(dict, ERRNO, ERRNOSTR) != 0) goto fail

static PyMethodDef errornoMethods[] = {
     {NULL, NULL, 0, NULL}
};

DL_EXPORT(void) initerrno(void)
{
     PyObject *dict=NULL, *m=NULL;

     m = Py_InitModule("errno", errornoMethods);
     assert(m);

     dict = PyModule_GetDict(m);
     if(!dict)
	  goto fail;

     ADD_ERRNO(GSL_SUCCESS , "GSL_SUCCESS" );
     ADD_ERRNO(GSL_FAILURE , "GSL_FAILURE" );
     ADD_ERRNO(GSL_CONTINUE, "GSL_CONTINUE");
     ADD_ERRNO(GSL_EDOM    , "GSL_EDOM"    );
     ADD_ERRNO(GSL_ERANGE  , "GSL_ERANGE"  );
     ADD_ERRNO(GSL_EFAULT  , "GSL_EFAULT"  );
     ADD_ERRNO(GSL_EINVAL  , "GSL_EINVAL"  );
     ADD_ERRNO(GSL_EFAILED , "GSL_EFAILED" );
     ADD_ERRNO(GSL_EFACTOR , "GSL_EFACTOR" );
     ADD_ERRNO(GSL_ESANITY , "GSL_ESANITY" );
     ADD_ERRNO(GSL_ENOMEM  , "GSL_ENOMEM"  );
     ADD_ERRNO(GSL_EBADFUNC, "GSL_EBADFUNC");
     ADD_ERRNO(GSL_ERUNAWAY, "GSL_ERUNAWAY");
     ADD_ERRNO(GSL_EMAXITER, "GSL_EMAXITER");
     ADD_ERRNO(GSL_EZERODIV, "GSL_EZERODIV");
     ADD_ERRNO(GSL_EBADTOL , "GSL_EBADTOL" );
     ADD_ERRNO(GSL_ETOL    , "GSL_ETOL"    );
     ADD_ERRNO(GSL_EUNDRFLW, "GSL_EUNDRFLW");
     ADD_ERRNO(GSL_EOVRFLW , "GSL_EOVRFLW" );
     ADD_ERRNO(GSL_ELOSS   , "GSL_ELOSS"   );
     ADD_ERRNO(GSL_EROUND  , "GSL_EROUND"  );
     ADD_ERRNO(GSL_EBADLEN , "GSL_EBADLEN" );
     ADD_ERRNO(GSL_ENOTSQR , "GSL_ENOTSQR" );
     ADD_ERRNO(GSL_ESING   , "GSL_ESING"   );
     ADD_ERRNO(GSL_EDIVERGE, "GSL_EDIVERGE");
     ADD_ERRNO(GSL_EUNSUP  , "GSL_EUNSUP"  );
     ADD_ERRNO(GSL_EUNIMPL , "GSL_EUNIMPL" );
     ADD_ERRNO(GSL_ECACHE  , "GSL_ECACHE"  );
     ADD_ERRNO(GSL_ETABLE  , "GSL_ETABLE"  );
     ADD_ERRNO(GSL_ENOPROG , "GSL_ENOPROG" );
     ADD_ERRNO(GSL_ENOPROGJ, "GSL_ENOPROGJ");
     ADD_ERRNO(GSL_ETOLF   , "GSL_ETOLF"   );
     ADD_ERRNO(GSL_ETOLX   , "GSL_ETOLX"   );
     ADD_ERRNO(GSL_ETOLG   , "GSL_ETOLG"   );
     ADD_ERRNO(GSL_EOF     , "GSL_EOF"     );
     ADD_ERRNO(PyGSL_ESTRIDE, "PyGSL_ESTRIDE");

     return;
 fail:
     fprintf(stderr, "Initialisation of module errorno failed!\n");
     return;
}


text:
#pragma once

#include <halley/utils/utils.h>
#include <vector>
#include <gsl/gsl>
#include "halley/data_structures/flat_map.h"
#include "halley/maths/vector2.h"
#include "halley/maths/rect.h"
#include "halley/file/path.h"
#include <map>
#include <unordered_map>
#include <cstdint>
#include <utility>
#include <set>
#include "halley/data_structures/maybe.h"
#include "halley/maths/colour.h"
#include "halley/maths/vector4.h"

namespace Halley {
	class String;

	class SerializerOptions {
	public:
		constexpr static int maxVersion = 1;

		int version = 0;
		bool exhaustiveDictionary = false;
		std::function<std::optional<size_t>(const String& string)> stringToIndex;
		std::function<const String&(size_t index)> indexToString;

		SerializerOptions() = default;
		SerializerOptions(int version)
			: version(version)
		{}
	};

	class SerializerState {};

	class ByteSerializationBase {
	public:
		ByteSerializationBase(SerializerOptions options)
			: options(std::move(options))
		{}

		SerializerState* setState(SerializerState* state);

		template <typename T>
		T* getState() const
		{
			return static_cast<T*>(state);
		}

		int getVersion() const { return version; }
		void setVersion(int v) { version = v; }

	protected:
		SerializerOptions options;

	private:
		SerializerState* state = nullptr;
		int version = 0;
	};

	class Serializer : public ByteSerializationBase {
	public:
		Serializer(SerializerOptions options);
		explicit Serializer(gsl::span<gsl::byte> dst, SerializerOptions options);

		template <typename T, typename std::enable_if<std::is_convertible<T, std::function<void(Serializer&)>>::value, int>::type = 0>
		static Bytes toBytes(const T& f, SerializerOptions options = {})
		{
			auto dry = Serializer(options);
			f(dry);
			Bytes result(dry.getSize());
			auto s = Serializer(gsl::as_writable_bytes(gsl::span<Halley::Byte>(result)), options);
			f(s);
			return result;
		}

		template <typename T, typename std::enable_if<!std::is_convertible<T, std::function<void(Serializer&)>>::value, int>::type = 0>
		static Bytes toBytes(const T& value, SerializerOptions options = {})
		{
			return toBytes([&value](Serializer& s) { s << value; }, options);
		}

		size_t getSize() const { return size; }

		Serializer& operator<<(bool val) { return serializePod(val); }
		Serializer& operator<<(int8_t val) { return serializeInteger(val); }
		Serializer& operator<<(uint8_t val) { return serializeInteger(val); }
		Serializer& operator<<(int16_t val) { return serializeInteger(val); }
		Serializer& operator<<(uint16_t val) { return serializeInteger(val); }
		Serializer& operator<<(int32_t val) { return serializeInteger(val); }
		Serializer& operator<<(uint32_t val) { return serializeInteger(val); }
		Serializer& operator<<(int64_t val) { return serializeInteger(val); }
		Serializer& operator<<(uint64_t val) { return serializeInteger(val); }
		Serializer& operator<<(float val) { return serializePod(val); }
		Serializer& operator<<(double val) { return serializePod(val); }

		Serializer& operator<<(const std::string& str);
		Serializer& operator<<(const String& str);
		Serializer& operator<<(const StringUTF32& str);
		Serializer& operator<<(const Path& path);
		Serializer& operator<<(gsl::span<const gsl::byte> span);
		Serializer& operator<<(const Bytes& bytes);

		template <typename T>
		Serializer& operator<<(const std::vector<T>& val)
		{
			unsigned int sz = static_cast<unsigned int>(val.size());
			*this << sz;
			for (unsigned int i = 0; i < sz; i++) {
				*this << val[i];
			}
			return *this;
		}

		template <typename T>
		Serializer& operator<<(const std::list<T>& val)
		{
			unsigned int sz = static_cast<unsigned int>(val.size());
			*this << sz;
			for (const auto& v: val) {
				*this << v;
			}
			return *this;
		}

		template <typename T, typename U>
		Serializer& operator<<(const FlatMap<T, U>& val)
		{
			*this << static_cast<unsigned int>(val.size());
			for (auto& kv : val) {
				*this << kv.first << kv.second;
			}
			return *this;
		}

		template <typename K, typename V, typename Cmp, typename Allocator>
		Serializer& operator<<(const std::map<K, V, Cmp, Allocator>& val)
		{
			*this << static_cast<unsigned int>(val.size());
			for (auto& kv : val) {
				*this << kv.first << kv.second;
			}
			return *this;
		}

		template <typename T, typename U>
		Serializer& operator<<(const std::unordered_map<T, U>& val)
		{
			std::map<T, U> m;
			for (auto& kv: val) {
				m[kv.first] = kv.second;
			}
			return (*this << m);
		}

		template <typename T>
		Serializer& operator<<(const std::set<T>& val)
		{
			unsigned int sz = static_cast<unsigned int>(val.size());
			*this << sz;
			for (auto& v: val) {
				*this << v;
			}
			return *this;
		}

		template <typename T>
		Serializer& operator<<(const Vector2D<T>& val)
		{
			return *this << val.x << val.y;
		}

		template <typename T>
		Serializer& operator<<(const Vector4D<T>& val)
		{
			return *this << val.x << val.y << val.z << val.w;
		}

		template <typename T>
		Serializer& operator<<(const Colour4<T>& val)
		{
			return *this << val.r << val.g << val.b << val.a;
		}

		template <typename T>
		Serializer& operator<<(const Rect2D<T>& val)
		{
			return *this << val.getTopLeft() << val.getBottomRight();
		}

		template <typename T, typename U>
		Serializer& operator<<(const std::pair<T, U>& p)
		{
			return *this << p.first << p.second;
		}

		template <typename T>
		Serializer& operator<<(const std::optional<T>& p)
		{
			if (p) {
				return *this << true << p.value();
			} else {
				return *this << false;
			}
		}

		template <typename T>
		Serializer& operator<<(const Range<T>& p)
		{
			return *this << p.start << p.end;
		}

		template <typename T, std::enable_if_t<std::is_enum<T>::value == true, int> = 0>
		Serializer& operator<<(const T& val)
		{
			using B = typename std::underlying_type<T>::type;
			return *this << B(val);
		}

		template <typename T, std::enable_if_t<std::is_enum<T>::value == false, int> = 0>
		Serializer& operator<<(const T& val)
		{
			val.serialize(*this);
			return *this;
		}

		size_t getPosition() const { return size; }

	private:
		size_t size = 0;
		gsl::span<gsl::byte> dst;
		bool dryRun;

		template <typename T>
		Serializer& serializePod(T val)
		{
			if (!dryRun) {
				memcpy(dst.data() + size, &val, sizeof(T));
			}
			size += sizeof(T);
			return *this;
		}

		template <typename T>
		Serializer& serializeInteger(T val)
		{
			if (options.version >= 1) {
				// Variable-length
				if constexpr (std::is_signed_v<T>) {
					serializeVariableInteger(static_cast<uint64_t>(val >= 0 ? val : -(val + 1)), val < 0);
				} else {
					serializeVariableInteger(val, {});
				}
				return *this;
			} else {
				// Fixed length
				return serializePod(val);
			}
		}

		void serializeVariableInteger(uint64_t val, std::optional<bool> sign);
	};

	class Deserializer : public ByteSerializationBase {
	public:
		Deserializer(gsl::span<const gsl::byte> src, SerializerOptions options = {});
		Deserializer(const Bytes& src, SerializerOptions options = {});

		template <typename T>
		static T fromBytes(const Bytes& src, SerializerOptions options = {})
		{
			T result;
			Deserializer s(src, std::move(options));
			s >> result;
			return result;
		}

		template <typename T>
		static T fromBytes(gsl::span<const gsl::byte> src, SerializerOptions options = {})
		{
			T result;
			Deserializer s(src, std::move(options));
			s >> result;
			return result;
		}

		template <typename T>
		static void fromBytes(T& target, const Bytes& src, SerializerOptions options = {})
		{
			Deserializer s(src, std::move(options));
			s >> target;
		}

		template <typename T>
		static void fromBytes(T& target, gsl::span<const gsl::byte> src, SerializerOptions options = {})
		{
			Deserializer s(src, std::move(options));
			s >> target;
		}

		Deserializer& operator>>(bool& val) { return deserializePod(val); }
		Deserializer& operator>>(int8_t& val) { return deserializeInteger(val); }
		Deserializer& operator>>(uint8_t& val) { return deserializeInteger(val); }
		Deserializer& operator>>(int16_t& val) { return deserializeInteger(val); }
		Deserializer& operator>>(uint16_t& val) { return deserializeInteger(val); }
		Deserializer& operator>>(int32_t& val) { return deserializeInteger(val); }
		Deserializer& operator>>(uint32_t& val) { return deserializeInteger(val); }
		Deserializer& operator>>(int64_t& val) { return deserializeInteger(val); }
		Deserializer& operator>>(uint64_t& val) { return deserializeInteger(val); }
		Deserializer& operator>>(float& val) { return deserializePod(val); }
		Deserializer& operator>>(double& val) { return deserializePod(val); }

		Deserializer& operator>>(std::string& str);
		Deserializer& operator>>(String& str);
		Deserializer& operator>>(StringUTF32& str);
		Deserializer& operator>>(Path& p);
		Deserializer& operator>>(gsl::span<gsl::byte> span);
		Deserializer& operator>>(Bytes& bytes);

		template <typename T>
		Deserializer& operator>>(std::vector<T>& val)
		{
			unsigned int sz;
			*this >> sz;
			ensureSufficientBytesRemaining(sz); // Expect at least one byte per vector entry

			val.clear();
			val.reserve(sz);
			for (unsigned int i = 0; i < sz; i++) {
				val.push_back(T());
				*this >> val[i];
			}
			return *this;
		}

		template <typename T>
		Deserializer& operator>>(std::list<T>& val)
		{
			unsigned int sz;
			*this >> sz;
			ensureSufficientBytesRemaining(sz); // Expect at least one byte per vector entry

			val.clear();
			for (unsigned int i = 0; i < sz; i++) {
				T v;
				*this >> v;
				val.push_back(std::move(v));
			}
			return *this;
		}

		template <typename T, typename U>
		Deserializer& operator>>(FlatMap<T, U>& val)
		{
			unsigned int sz;
			*this >> sz;
			ensureSufficientBytesRemaining(sz * 2); // Expect at least two bytes per map entry

			std::vector<std::pair<T, U>> tmpData(sz);
			for (unsigned int i = 0; i < sz; i++) {
				*this >> tmpData[i].first >> tmpData[i].second;
			}
			val = FlatMap<T, U>(boost::container::ordered_unique_range_t(), tmpData.begin(), tmpData.end());
			return *this;
		}

		template <typename K, typename V, typename Cmp, typename Allocator>
		Deserializer& operator >> (std::map<K, V, Cmp, Allocator>& val)
		{
			unsigned int sz;
			*this >> sz;
			ensureSufficientBytesRemaining(size_t(sz) * 2); // Expect at least two bytes per map entry

			for (unsigned int i = 0; i < sz; i++) {
				K key;
				V value;
				*this >> key >> value;
				val[key] = std::move(value);
			}
			return *this;
		}

		template <typename T, typename U>
		Deserializer& operator >> (std::unordered_map<T, U>& val)
		{
			unsigned int sz;
			*this >> sz;
			ensureSufficientBytesRemaining(sz * 2); // Expect at least two bytes per map entry

			for (unsigned int i = 0; i < sz; i++) {
				T key;
				U value;
				*this >> key >> value;
				val[key] = std::move(value);
			}
			return *this;
		}

		template <typename T>
		Deserializer& operator>>(std::set<T>& val)
		{
			unsigned int sz;
			*this >> sz;
			ensureSufficientBytesRemaining(sz); // Expect at least one byte per set entry

			val.clear();
			for (unsigned int i = 0; i < sz; i++) {
				T v;
				*this >> v;
				val.insert(std::move(v));
			}
			return *this;
		}

		template <typename T>
		Deserializer& operator>>(Vector2D<T>& val)
		{
			*this >> val.x;
			*this >> val.y;
			return *this;
		}

		template <typename T>
		Deserializer& operator>>(Vector4D<T>& val)
		{
			*this >> val.x;
			*this >> val.y;
			*this >> val.z;
			*this >> val.w;
			return *this;
		}

		template <typename T>
		Deserializer& operator>>(Colour4<T>& val)
		{
			*this >> val.r;
			*this >> val.g;
			*this >> val.b;
			*this >> val.a;
			return *this;
		}

		template <typename T>
		Deserializer& operator>>(Rect2D<T>& val)
		{
			Vector2D<T> p1, p2;
			*this >> p1;
			*this >> p2;
			val = Rect2D<T>(p1, p2);
			return *this;
		}

		template <typename T, typename U>
		Deserializer& operator>>(std::pair<T, U>& p)
		{
			return *this >> p.first >> p.second;
		}

		template <typename T>
		Deserializer& operator>>(Range<T>& p)
		{
			return *this >> p.start >> p.end;
		}

		template <typename T>
		Deserializer& operator>>(std::optional<T>& p)
		{
			bool present;
			*this >> present;
			if (present) {
				T tmp;
				*this >> tmp;
				p = tmp;
			} else {
				p = std::optional<T>();
			}
			return *this;
		}

		template <typename T, std::enable_if_t<std::is_enum<T>::value == true, int> = 0>
		Deserializer& operator>>(T& val)
		{
			typename std::underlying_type<T>::type tmp;
			*this >> tmp;
			val = T(tmp);
			return *this;
		}

		template <typename T, std::enable_if_t<std::is_enum<T>::value == false, int> = 0>
		Deserializer& operator>>(T& val)
		{
			val.deserialize(*this);
			return *this;
		}

		template <typename T>
		void peek(T& val)
		{
			const auto oldPos = pos;
			*this >> val;
			pos = oldPos;
		}

		size_t getPosition() const { return pos; }

	private:
		size_t pos = 0;
		gsl::span<const gsl::byte> src;

		template <typename T>
		Deserializer& deserializePod(T& val)
		{
			ensureSufficientBytesRemaining(sizeof(T));
			memcpy(&val, src.data() + pos, sizeof(T));
			pos += sizeof(T);
			return *this;
		}

		template <typename T>
		Deserializer& deserializeInteger(T& val)
		{
			if (options.version >= 1) {
				// Variable-length
				bool sign;
				uint64_t temp;
				deserializeVariableInteger(temp, sign, std::is_signed_v<T>);
				if (sign) {
					int64_t signedTemp = -int64_t(temp) - 1;
					val = static_cast<T>(signedTemp);
				} else {
					val = static_cast<T>(temp);
				}
				return *this;
			} else {
				// Fixed length
				return deserializePod(val);
			}
		}

		void deserializeVariableInteger(uint64_t& val, bool& sign, bool isSigned);

		void ensureSufficientBytesRemaining(size_t bytes);
		size_t getBytesRemaining() const;
	};
}

ext: h

max_stars_repo_name: JLouhela/halley

max_issues_repo_name: JLouhela/halley

max_forks_repo_name: JLouhela/halley


text:
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <gsl/gsl_statistics_double.h>
#include "tz_error.h"
#include "tz_constant.h"
#include "tz_iarray.h"
#include "tz_image_lib.h"
#include "tz_objdetect.h"
#include "tz_imatrix.h"
#include "tz_stack_math.h"
#include "tz_stack_bwdist.h"
#include "tz_pixel_array.h"
#include "tz_stack_bwmorph.h"
#include "tz_dimage_lib.h"
#include "tz_darray.h"
#include "tz_u8array.h"
#include "tz_stack_objlabel.h"
#include "tz_voxel_linked_list.h"
#include "tz_voxel_graphics.h"
#include "tz_dmatrix.h"
#include "tz_stack_sampling.h"
#include "tz_stack_threshold.h"

INIT_EXCEPTION_MAIN(e)

int main(int argc, char* argv[])
{
  char neuron_name[100];
  if (argc == 1) {
    strcpy(neuron_name, "fly_neuron");
  } else {
    strcpy(neuron_name, argv[1]);
  }

  char file_path[100];
  sprintf(file_path, "../data/%s/mask.tif", neuron_name);

  Stack *mask = Read_Stack(file_path);

  sprintf(file_path, "../data/%s/seeds.pa", neuron_name);
  Pixel_Array *pa = Pixel_Array_Read(file_path);
  double *pa_array = (double *) pa->array;

  sprintf(file_path, "../data/%s/seed_offset.ar", neuron_name);
  int seed_number = pa->size;
  int *seed_offset = (int *) malloc(sizeof(int) * seed_number);
  iarray_read(file_path, seed_offset, &seed_number);

  int max_idx;
  double max_dist = darray_max(pa_array, pa->size, &max_idx);
  max_idx = seed_offset[max_idx];

  Stack_Label_Object_Dist_N(mask, NULL, max_idx, 1, 2, max_dist * 2.0, 26);

  Stack_Threshold_Binarize(mask, 1);

  sprintf(file_path, "../data/%s/blobmask.tif", neuron_name);
  Write_Stack(file_path, mask);

  free(seed_offset);
  Kill_Stack(mask);

  return 0;
}

ext: c


text:
#include "mpi.h"
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>
#include <string.h>
#include <stdbool.h>

#ifdef __APPLE__

#include <cblas.h>
#include <libiomp/omp.h>

#define set_num_threads(x) openblas_set_num_threads(x)
#define get_num_threads() openblas_get_num_threads()
#else
#include <mkl.h>
#include <omp.h>
#define set_num_threads(x) mkl_set_num_threads(x)
#define get_num_threads() mkl_get_num_threads()
#endif

double norm_inf_mat_Y11(double *X, int n);

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

    char *files[] = {
            //"1200_m1_t2_hi_c363417000.7.mat.bin",
            //"1200_m1_t7_hi_c536234414.4.mat.bin",
            //"2340_m2_a1e15_to_c3291.0.mat.bin",
            //"5900_m1_t1_to_c3205771604.1.mat.bin",
            //"1200_m1_t3_hi_c70487569.2.mat.bin",
            //"2340_m1_t1_to_c1446918556.5.mat.bin",
            //"2340_m2_a1e5_to_c8.0.mat.bin",
            //"5900_m2_a1e15_to_c4198.3.mat.bin",
            "494_bus.mtx.bin",
            //"bcsstk13.mtx.bin",
            //"bcsstk27.mtx.bin",
            //"ex9.mtx.bin",
            //"msc01050.mtx.bin",
            //"bcsstk15.mtx.bin",
            //"cage9.mtx.bin",
            //"gyro.mtx.bin",
            "tomography.mtx.bin"
    };
    int num_of_files = sizeof(files) / sizeof(char *);

    char *file_dir = NULL;
    int option = 0, omp_threads = 1;
    while ((option = getopt(argc, argv, "f:t:")) != -1) {
        switch (option) {
            case 'f':
                file_dir = optarg;
                break;
            case 't':
                omp_threads = atoi(optarg);
                break;
            default:
                printf("Usage: mpi_mexp -f string -t num_threads \n");
                return 0;
        }
    }
    if (file_dir == NULL) {
        printf("Usage: mpi_mexp -f string -t num_threads \n");
        return 0;
    }

    set_num_threads(omp_threads);
    omp_set_num_threads(omp_threads);

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

    int qProcs = (int) sqrt(numTasks);
    if (qProcs - sqrt(numTasks) != 0) {
        if (rank == 0)
            printf("np must be a square number\n");
        MPI_Finalize();
        return 0;
    }

    int q_root = qProcs - 1;
    int world_root = numTasks - 1;
    int n, dims[2] = {qProcs, qProcs}, periods[2] = {0, 0}, reorder = 0, coords[2], rc, i;
    MPI_Comm cartComm;
    MPI_Cart_create(MPI_COMM_WORLD, 2, dims, periods, reorder, &cartComm);
    MPI_Comm_rank(cartComm, &rank);
    MPI_Cart_coords(cartComm, rank, 2, coords);

    MPI_Comm rowComm, colComm;
    MPI_Comm_split(MPI_COMM_WORLD, coords[0], rank, &rowComm);
    MPI_Comm_split(MPI_COMM_WORLD, coords[1], rank, &colComm);

    int row_rank, col_rank;
    MPI_Comm_rank(colComm, &col_rank);
    MPI_Comm_rank(rowComm, &row_rank);

    char filename[1024];
    while (num_of_files-- > 0) {

        memset(filename, 0, strlen(filename));
        strcpy(filename, file_dir);
        strcat(filename, files[num_of_files]);

        MPI_Barrier(MPI_COMM_WORLD);

        if (rank == 0)
            printf("\n%s\n", filename);
        MPI_File matFile;
        rc = MPI_File_open(cartComm, filename, MPI_MODE_RDONLY, MPI_INFO_NULL, &matFile);
        if (rc && (rank == 0)) {
            printf("Unable to open file %s\n", filename);
            fflush(stdout);
        }

        MPI_Status status;
        MPI_File_read(matFile, &n, 1, MPI_INT, &status);

        int blockSize = 0, lastBlock = 0;
        blockSize = (n + qProcs) / qProcs;  /* number of rows in _block_ */
        lastBlock = n + 1 - (qProcs - 1) * blockSize;

        if (rank == 0) {
            printf("size: %d, %d ,%d\n", n, blockSize, lastBlock);
        }

        double *LocalA = (double *) malloc(blockSize * blockSize * sizeof(double));

        MPI_Datatype readFileType;
        MPI_Type_vector(blockSize, blockSize, n, MPI_DOUBLE, &readFileType);
        MPI_Type_commit(&readFileType);
        MPI_Offset offset = (MPI_Offset) (1 * sizeof(int) +
                                          sizeof(double) * (blockSize * coords[1] + blockSize * n * coords[0]));
        MPI_File_set_view(matFile, offset, MPI_DOUBLE, readFileType,
                          "native", MPI_INFO_NULL);

        MPI_File_read_at(matFile, 0, LocalA, blockSize * blockSize, MPI_DOUBLE, &status);
        MPI_Type_free(&readFileType);
        if (coords[0] == q_root) {
            memset(LocalA + (lastBlock - 1) * blockSize, 0,
                   (blockSize + 1 - lastBlock) * blockSize * sizeof(double));
        }
        if (coords[1] == q_root) {
            int kk = 0;
            for (kk = 0; kk < blockSize; ++kk) {
                memset(LocalA + kk * blockSize + lastBlock - 1, 0, (blockSize + 1 - lastBlock) * sizeof(double));
            }
        }

        double loc_d_var;
        double *loc_vec_var = (double *) malloc(blockSize * sizeof(double));
        memset(loc_vec_var, 0, blockSize * sizeof(double));
        offset = (MPI_Offset) (1 * sizeof(int) + sizeof(double) * n * n);
        MPI_File_set_view(matFile, offset, MPI_DOUBLE, MPI_DOUBLE,
                          "native", MPI_INFO_NULL);
        offset = coords[0] * blockSize;
        MPI_File_read_at_all(matFile, offset, loc_vec_var, coords[0] == qProcs - 1 ? lastBlock : blockSize, MPI_DOUBLE,
                             &status);

        MPI_File_close(&matFile);

        if (coords[1] == qProcs - 1) {
            cblas_dcopy(blockSize, loc_vec_var, 1, LocalA + lastBlock - 1, blockSize);
        }

        double mpi_start = MPI_Wtime();
        //compute inf norm of matA
        for (i = 0; i < (col_rank == q_root ? lastBlock - 1 : blockSize); ++i) {
            loc_vec_var[i] = cblas_dasum((row_rank == q_root ? lastBlock - 1 : blockSize), LocalA + i * blockSize,
                                         1);
        }

        if (row_rank == q_root) {
            MPI_Reduce(MPI_IN_PLACE, loc_vec_var, blockSize, MPI_DOUBLE, MPI_SUM, q_root, rowComm);
        } else {
            MPI_Reduce(loc_vec_var, loc_vec_var, blockSize, MPI_DOUBLE, MPI_SUM, q_root, rowComm);
        }

        if (row_rank == q_root) {
            loc_d_var = 0;
            for (i = 0; i < (col_rank == q_root ? lastBlock - 1 : blockSize); ++i) {
                if (loc_vec_var[i] > loc_d_var)
                    loc_d_var = loc_vec_var[i];
            }
            if (col_rank == q_root) {
                MPI_Reduce(MPI_IN_PLACE, &loc_d_var, 1, MPI_DOUBLE, MPI_MAX, q_root, colComm);
            } else {
                MPI_Reduce(&loc_d_var, &loc_d_var, 1, MPI_DOUBLE, MPI_MAX, q_root, colComm);
            }
        }

        MPI_Bcast(&loc_d_var, 1, MPI_DOUBLE, world_root, MPI_COMM_WORLD);

        //init Y martix
        for (i = 0; i < (col_rank == q_root ? lastBlock - 1 : blockSize); ++i) {
            cblas_dscal((row_rank == q_root ? lastBlock - 1 : blockSize), -1 / loc_d_var, LocalA + i * blockSize,
                        1);
        }
        if (row_rank == q_root) {
            cblas_dscal((col_rank == q_root ? lastBlock - 1 : blockSize), 1 / loc_d_var, LocalA + lastBlock - 1,
                        blockSize);
        }
        if (row_rank == col_rank) {
            for (i = 0; i < (row_rank == q_root ? lastBlock - 1 : blockSize); ++i) {
                LocalA[i * blockSize + i] += 1;
            }
        }
        if (rank == world_root) {
            LocalA[(lastBlock - 1) * (blockSize + 1)] = 1;
        }

        //loop
        double tolerance = 1e-8;
        int loopBreak = 0, iteration = 0;
        double *col_mat = (double *) malloc(qProcs * blockSize * blockSize * sizeof(double));

        double *row_mat = (double *) malloc(qProcs * blockSize * blockSize * sizeof(double));

        while (true) {
            //norm mat
            //compute inf norm of matA
            for (i = 0; i < (col_rank == q_root ? lastBlock - 1 : blockSize); ++i) {
                loc_vec_var[i] = cblas_dasum((row_rank == q_root ? lastBlock - 1 : blockSize),
                                             LocalA + i * blockSize, 1);
            }

            if (row_rank == q_root) {
                MPI_Reduce(MPI_IN_PLACE, loc_vec_var, blockSize, MPI_DOUBLE, MPI_SUM, q_root, rowComm);
            } else {
                MPI_Reduce(loc_vec_var, loc_vec_var, blockSize, MPI_DOUBLE, MPI_SUM, q_root, rowComm);
            }

            if (row_rank == q_root) {
                loc_d_var = 0;
                for (i = 0; i < (col_rank == q_root ? lastBlock - 1 : blockSize); ++i) {
                    if (loc_vec_var[i] > loc_d_var)
                        loc_d_var = loc_vec_var[i];
                }
                if (col_rank == q_root) {
                    MPI_Reduce(MPI_IN_PLACE, &loc_d_var, 1, MPI_DOUBLE, MPI_MAX, q_root, colComm);
                } else {
                    MPI_Reduce(&loc_d_var, &loc_d_var, 1, MPI_DOUBLE, MPI_MAX, q_root, colComm);
                }

                size_t index = cblas_idamax((col_rank == q_root ? lastBlock - 1 : blockSize),
                                            LocalA + lastBlock - 1,
                                            blockSize);
                double norm2b = 0;
                MPI_Reduce(&LocalA[index * blockSize + lastBlock - 1], &norm2b, 1, MPI_DOUBLE, MPI_MAX, q_root,
                           colComm);

                if (col_rank == q_root) {
                    loopBreak = (loc_d_var / norm2b < tolerance) ? 1 : 0;
                }

            }
            MPI_Bcast(&loopBreak, 1, MPI_INT, world_root, MPI_COMM_WORLD);
            if (loopBreak > 0)
                break;

            iteration++;

            //gather loc mat
            MPI_Allgather(LocalA, blockSize * blockSize, MPI_DOUBLE, col_mat, blockSize * blockSize, MPI_DOUBLE,
                          colComm);
            MPI_Allgather(LocalA, blockSize * blockSize, MPI_DOUBLE, row_mat, blockSize * blockSize, MPI_DOUBLE,
                          rowComm);
            cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, blockSize, blockSize, blockSize,
                        1, row_mat, blockSize, col_mat, blockSize, 0.0, LocalA, blockSize);
            for (i = 1; i < qProcs; ++i) {
                cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, blockSize, blockSize, blockSize,
                            1, row_mat + i * (blockSize * blockSize), blockSize,
                            col_mat + i * (blockSize * blockSize), blockSize, 1, LocalA, blockSize);
            }
        }

        free(row_mat);
        free(col_mat);
        free(loc_vec_var);

        if (rank == world_root) {
            printf("NumProc: %d\n", numTasks);
            printf("mpi_mexp_iter: %d\t", iteration);
            printf("mpi_mexp_time: %f\n", MPI_Wtime() - mpi_start);

        }

        if (row_rank == q_root) {
            for (i = 0; i < (row_rank == qProcs - 1 ? lastBlock : blockSize); ++i) {
                if (fabs(LocalA[i * blockSize + lastBlock - 1] - 1) > 1e-3) {
                    printf("MPI_MEXP_ERR %f \n", LocalA[i * blockSize + lastBlock - 1]);
                    break;
                }
            }
        }
        free(LocalA);

        /*
        * MPI_CG
        */
        rc = MPI_File_open(cartComm, filename, MPI_MODE_RDONLY, MPI_INFO_NULL, &matFile);
        if (rc && (rank == 0)) {
            printf("Unable to open file %s\n", filename);
            fflush(stdout);
        }

        MPI_File_read(matFile, &n, 1, MPI_INT, &status);
        blockSize = (n + qProcs - 1) / qProcs;  /* number of rows in _block_ */
        lastBlock = n - (qProcs - 1) * blockSize;

        LocalA = (double *) malloc(blockSize * blockSize * sizeof(double));

        MPI_Type_vector(blockSize, blockSize, n, MPI_DOUBLE, &readFileType);
        MPI_Type_commit(&readFileType);
        offset = (MPI_Offset) (1 * sizeof(int) +
                               sizeof(double) * (blockSize * coords[1] + blockSize * n * coords[0]));
        MPI_File_set_view(matFile, offset, MPI_DOUBLE, readFileType,
                          "native", MPI_INFO_NULL);

        MPI_File_read_at(matFile, 0, LocalA, blockSize * blockSize, MPI_DOUBLE, &status);
        MPI_Type_free(&readFileType);
        if (coords[0] == qProcs - 1) {
            memset(LocalA + lastBlock * blockSize, 0, (blockSize - lastBlock) * blockSize * sizeof(double));
        }
        if (coords[1] == qProcs - 1) {
            int kk = 0;
            for (kk = 0; kk < blockSize; ++kk) {
                memset(LocalA + kk * blockSize + lastBlock, 0, (blockSize - lastBlock) * sizeof(double));
            }
        }

        double *LocalB = (double *) malloc(blockSize * sizeof(double));
        memset(LocalB, 0, blockSize * sizeof(double));
        offset = (MPI_Offset) (1 * sizeof(int) + sizeof(double) * n * n);
        MPI_File_set_view(matFile, offset, MPI_DOUBLE, MPI_DOUBLE,
                          "native", MPI_INFO_NULL);
        offset = coords[1] * blockSize;
        MPI_File_read_at_all(matFile, offset, LocalB, coords[1] == qProcs - 1 ? lastBlock : blockSize, MPI_DOUBLE,
                             &status);

        MPI_File_close(&matFile);
        mpi_start = MPI_Wtime();

        double alpha, beta, rho_new = 0.0, rho_old = 0.0;
        iteration = 0;
        double *loc_Vr = LocalB; //(double *) malloc(blockSize * sizeof(double));
        double *loc_Vp = (double *) malloc(blockSize * sizeof(double));
        double *loc_Vw = (double *) malloc(blockSize * sizeof(double));
        double *loc_Vx = (double *) malloc(blockSize * sizeof(double));

        //localr <- localB
        memset(loc_Vx, 0, blockSize * sizeof(double));
        //cblas_dcopy(blockSize, LocalB, 1, loc_Vr, 1);
        double b_norm2 = 0;
        if (col_rank == qProcs - 1) {
            double loc_b_dot = cblas_ddot((row_rank == qProcs - 1 ? lastBlock : blockSize), loc_Vr, 1, loc_Vr, 1);
            MPI_Reduce(&loc_b_dot, &rho_new, 1, MPI_DOUBLE, MPI_SUM, qProcs - 1, rowComm);
            b_norm2 = sqrt(rho_new);
        }

        loopBreak = 0;
        while (true) {
            if (rank == world_root) {
                loopBreak = (sqrt(rho_new) / b_norm2 < tolerance) ? 1 : 0;
            }
            MPI_Bcast(&loopBreak, 1, MPI_INT, world_root, MPI_COMM_WORLD);
            if (loopBreak > 0)
                break;

            iteration++;
            if (iteration == 1) {
                cblas_dcopy(blockSize, loc_Vr, 1, loc_Vp, 1);
            } else {
                if (rank == world_root) {
                    beta = rho_new / rho_old;
                }
                MPI_Bcast(&beta, 1, MPI_DOUBLE, world_root, MPI_COMM_WORLD);
                cblas_dscal((row_rank == qProcs - 1 ? lastBlock : blockSize), beta, loc_Vp, 1);
                cblas_daxpy((row_rank == qProcs - 1 ? lastBlock : blockSize), 1.0, loc_Vr, 1, loc_Vp, 1);
            }

            cblas_dgemv(CblasRowMajor, CblasNoTrans, (col_rank == qProcs - 1 ? lastBlock : blockSize),
                        (row_rank == qProcs - 1 ? lastBlock : blockSize),
                        1.0, LocalA, blockSize, loc_Vp, 1, 0.0, loc_Vw, 1);

            if (col_rank == row_rank) {
                MPI_Reduce(MPI_IN_PLACE, loc_Vw, blockSize, MPI_DOUBLE, MPI_SUM, col_rank, rowComm);
            } else {
                MPI_Reduce(loc_Vw, loc_Vw, blockSize, MPI_DOUBLE, MPI_SUM, col_rank, rowComm);
            }

            MPI_Bcast(loc_Vw, (row_rank == qProcs - 1 ? lastBlock : blockSize), MPI_DOUBLE, row_rank, colComm);

            if (col_rank == qProcs - 1) {
                double loc_pw_dot = cblas_ddot((row_rank == qProcs - 1 ? lastBlock : blockSize), loc_Vp, 1, loc_Vw,
                                               1), glo_pw_dot;
                MPI_Reduce(&loc_pw_dot, &glo_pw_dot, 1, MPI_DOUBLE, MPI_SUM, qProcs - 1, rowComm);
                if (rank == world_root) {
                    alpha = rho_new / glo_pw_dot;
                }
            }
            MPI_Bcast(&alpha, 1, MPI_DOUBLE, world_root, MPI_COMM_WORLD);

            cblas_daxpy((row_rank == qProcs - 1 ? lastBlock : blockSize), -alpha, loc_Vw, 1, loc_Vr, 1);

            if (coords[0] == qProcs - 1) {
                rho_old = rho_new;
                cblas_daxpy((row_rank == qProcs - 1 ? lastBlock : blockSize), alpha, loc_Vp, 1, loc_Vx, 1);
                double loc_r_dot = cblas_ddot(blockSize, loc_Vr, 1, loc_Vr, 1);

                MPI_Reduce(&loc_r_dot, &rho_new, 1, MPI_DOUBLE, MPI_SUM, qProcs - 1, rowComm);
            }
        }

        if (rank == world_root) {
            printf("mpi_cg_it: %d\t", iteration);
            printf("mpi_cg_time: %f\n", MPI_Wtime() - mpi_start);

        }
        free(LocalA);
        free(loc_Vp);
        free(loc_Vr);
        free(loc_Vw);

        if (col_rank == qProcs - 1) {
            for (rc = 0; rc < (row_rank == qProcs - 1 ? lastBlock : blockSize); ++rc) {
                if (fabs(loc_Vx[rc] - 1) > 1e-3) {
                    printf("MPI_CG_ERR %f \n", loc_Vx[rc]);
                    break;
                }
            }
        }
        free(loc_Vx);
    }

    MPI_Comm_free(&cartComm);
    MPI_Comm_free(&colComm);
    MPI_Comm_free(&rowComm);

    MPI_Finalize();
}

max_stars_repo_name: baishuai/MEXP

max_issues_repo_name: baishuai/MEXP

max_forks_repo_name: baishuai/MEXP


text:
#pragma once

#include "RTTI.h"
#include <d3d11.h>
#include <stack>
#include <gsl\gsl>
#include <vector>

namespace Library
{
	class RenderTarget : public RTTI
	{
		RTTI_DECLARATIONS(RenderTarget, RTTI)

	public:
		RenderTarget() = default;
		RenderTarget(const RenderTarget&) = delete;
		RenderTarget(RenderTarget&&) = default;
		RenderTarget& operator=(const RenderTarget&) = delete;
		RenderTarget& operator=(RenderTarget&&) = default;
		virtual ~RenderTarget() = default;

		virtual void Begin() = 0;
		virtual void End() = 0;

	protected:
		struct RenderTargetData
		{
			uint32_t ViewCount() const { return gsl::narrow_cast<uint32_t>(RenderTargetViews.size()); }

			std::vector<ID3D11RenderTargetView*> RenderTargetViews;
			gsl::not_null<ID3D11DepthStencilView*> DepthStencilView;
			D3D11_VIEWPORT Viewport;

			RenderTargetData(const gsl::span<ID3D11RenderTargetView*>& renderTargetViews, gsl::not_null<ID3D11DepthStencilView*> depthStencilView, const D3D11_VIEWPORT& viewport) :
				RenderTargetViews(renderTargetViews.begin(), renderTargetViews.end()), DepthStencilView(depthStencilView), Viewport(viewport) { }
		};

		void Begin(gsl::not_null<ID3D11DeviceContext*> deviceContext, const gsl::span<ID3D11RenderTargetView*>& renderTargetViews, gsl::not_null<ID3D11DepthStencilView*> depthStencilView, const D3D11_VIEWPORT& viewport);
		void End(gsl::not_null<ID3D11DeviceContext*> deviceContext);
		void RebindCurrentRenderTargets(gsl::not_null<ID3D11DeviceContext*> deviceContext);

	private:
		static std::stack<RenderTargetData> sRenderTargetStack;
	};
}


text:
/*
* -----------------------------------------------------------------
*  In Situ Adaptive Tabulation Library --- isat_lib.c
*  Version: 1.6180
*  Date: Nov 15, 2010
* -----------------------------------------------------------------
*  Programmer: Americo Barbosa da Cunha Junior
*              americo.cunhajr@gmail.com
* -----------------------------------------------------------------
*  Copyright (c) 2010 by Americo Barbosa da Cunha Junior
*
*  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.
*
*  A copy of the GNU General Public License is available in
*  LICENSE.txt or http://www.gnu.org/licenses/.
* -----------------------------------------------------------------
*  This is the implementation file of a library
*  to work with ISAT algorithm.
* -----------------------------------------------------------------
*/

#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>

#include "../include/thrm_lib.h"
#include "../include/ell_lib.h"
#include "../include/ode_lib.h"
#include "../include/isat_lib.h"

/*
*------------------------------------------------------------
*   isat_alloc
*
*   This function alocates memory for a struct isat.
*
*   Output:
*   isat - pointer to a struct isat
*
*   last update: Jul 13, 2010
*------------------------------------------------------------
*/

isat_wrk *isat_alloc()
{
    /* memory allocation for isat_wrk */
    isat_wrk *isat = NULL;
    isat = (isat_wrk *) malloc(sizeof(isat_wrk));
    if ( isat == NULL )
        return NULL;

    /* setting isat_wrk elements equal NULL and 0.0 */
    isat->root     = NULL;
    isat->lf       = 0;
    isat->nd       = 0;
    isat->add      = 0;
    isat->grw      = 0;
    isat->rtv      = 0;
    isat->dev      = 0;
    isat->hgt      = 0;
    isat->max_lf   = 0;
    isat->time_add = 0.0;
    isat->time_grw = 0.0;
    isat->time_rtv = 0.0;
    isat->time_dev = 0.0;

    return isat;
}
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   isat_free
*
*   This function frees the memory used by a struct isat.
*
*   Input:
*   isat - pointer to a struct isat
*
*   last update: May 10, 2009
*------------------------------------------------------------
*/

void isat_free(void **isat_bl)
{
    isat_wrk *isat = NULL;

    /* checking if isat_bl is NULL */
    if ( *isat_bl == NULL )
        return;

    isat = (isat_wrk *) (*isat_bl);

    /* releasing allocated memory by isat_wrk elements */
    if( isat->root != NULL )
    {
        bst_node_free((void **)&(isat->root));
        isat->root = NULL;
    }

    /* releasing allocated memory by isat_wrk */
    free(*isat_bl);
    *isat_bl = NULL;

    return;
}
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   isat_set
*
*   This function initiates ISAT binary search tree.
*
*   Input:
*   isat - pointer to a struct isat
*
*   last update: May 13, 2009
*------------------------------------------------------------
*/

int isat_set(isat_wrk *isat)
{
    /* checking if isat is NULL */
    if ( isat == NULL )
        return GSL_EINVAL;

    /* memory allocation for binary search tree root */
    if ( isat->root == NULL )
    {
        isat->root = bst_node_alloc();
        if ( isat->root == NULL )
        {
            free(isat);
            isat = NULL;
            return GSL_EINVAL;
        }
    }

    return GSL_SUCCESS;
}
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   isat_statistics
*
*   This function prints on the screen informations
*   about the isat workspace.
*
*   Input:
*   isat - pointer to a struct isat
*
*   last update: Jul 13, 2010
*------------------------------------------------------------
*/

void isat_statistics(isat_wrk *isat)
{
    double time_add;
    double time_grw;
    double time_rtv;
    double time_dev;

    time_add = (double) (isat->time_add / CLOCKS_PER_SEC) / isat->add;
    time_grw = (double) (isat->time_grw / CLOCKS_PER_SEC) / isat->grw;
    time_rtv = (double) (isat->time_rtv / CLOCKS_PER_SEC) / isat->rtv;
    time_dev = (double) (isat->time_dev / CLOCKS_PER_SEC) / isat->dev;

    printf("\n ISAT statistics:");
    printf("\n # of adds         = %d",   isat->add);
    printf("\n # of grows        = %d",   isat->grw);
    printf("\n # of retrieves    = %d",   isat->rtv);
    printf("\n # of dir. eval.   = %d",   isat->dev);
    printf("\n # of leaves       = %d",   isat->lf);
    printf("\n # of nodes        = %d",   isat->nd);
    printf("\n tree height       = %d\n", isat->hgt);

    printf("\n average values for CPU time (s):");
    printf("\n add: %+.6e",time_add);
    printf("\n grw: %+.6e",time_grw);
    printf("\n rtv: %+.6e",time_rtv);
    printf("\n dev: %+.6e",time_dev);

    return;
}
/*------------------------------------------------------------*/

/*
* -----------------------------------------------------------------
*   isat_input
*
*   This function receives ISAT parameters.
*
*   Input:
*   max_lf - maximum value of tree leaves
*   etol   - ISAT error tolerance
*   n0     - factor to multiply by unit roundoff
*
*   Output:
*   success or error
*
*   last update: Feb 25, 2010
* -----------------------------------------------------------------
*/

int isat_input(unsigned int *max_lf,double *etol,double *n0)
{
    printf("\n Input ISAT parameters:\n");

    printf("\n maximum of tree leaves:");
    scanf("%d", max_lf);
    printf("\n %d\n", *max_lf);
    if( *max_lf <= 0 )
	GSL_ERROR(" max_lf must be a positive integer (max_lf > 0)",GSL_EINVAL);

    printf("\n ISAT error tolerance:");
    scanf("%lf", etol);
    printf("\n %+.1e\n", *etol);
    if( *etol < 0.0 )
	GSL_ERROR(" etol must be grather than zero (etol > 0.0)",GSL_EINVAL);

    printf("\n Unit roundoff multiple:");
    scanf("%lf", n0);
    printf("\n %+.1e\n", *n0);
    if( *n0 < 0.0 )
	GSL_ERROR(" n0 must be grather than zero (n0 > 0.0)",GSL_EINVAL);

    return GSL_SUCCESS;
}
/*----------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   isat_eoa_matrix
*
*   This function computes the EOA matrix in Cholesky form.
*
*   Input:
*   A    - gradient matrix
*   n0   - factor to multiply by unit roundoff
*   etol - error tolerance
*
*   Output:
*   L    - EOA Cholesky matrix
*
*   last update: Feb 19, 2010
*------------------------------------------------------------
*/

void isat_eoa_mtrx(gsl_matrix *A,double etol,double n0,
						gsl_matrix *L)
{
    unsigned int i;
    double  eps_max   = 0.5;
    double  eps_min   = etol/(n0*DBL_EPSILON);
    gsl_matrix *Aetol = NULL;
    gsl_matrix *V     = NULL;
    gsl_vector *sig   = NULL;

    /* memory allocation */
    Aetol = gsl_matrix_calloc(A->size2,A->size2);
    V     = gsl_matrix_calloc(A->size2,A->size2);
    sig   = gsl_vector_calloc(A->size2);

    /* Aetol := A */
    gsl_matrix_memcpy(Aetol,A);

    /* Aetol := (1/etol).A */
    gsl_matrix_scale (Aetol,1.0/etol);

    /* Aetol = U*sig*V^T */
    ell_psd2eig(Aetol,V,sig);

    /* eliminating small and large singular values */
    for ( i = 0; i < sig->size; i++ )
        sig->data[i] = GSL_MIN( GSL_MAX(sig->data[i],eps_max), eps_min);

    /* V*sig^2*V^T = L*L^T */
    ell_eig2chol(V,sig,L);

    /* releasing allocated memory */
    gsl_vector_free(sig);
    gsl_matrix_free(V);
    gsl_matrix_free(Aetol);
    sig   = NULL;
    V     = NULL;
    Aetol = NULL;

    return;
}
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   isat_lerror
*
*   This function computes the local error defined as
*
*   eps = 2-norm(Rl(phi)- R(phi)) where
*
*   Rl(phi) = R(phi0) + A*(phi-phi0).
*
*   Input:
*   Rphi  - reaction mapping of phi
*   Rphi0 - reaction mapping of phi0
*   A     - mapping gradient matrix
*   phi   - query composition
*   phi0  - initial composition
*
*   Output:
*   eps   - local error
*
*   last update: Feb 19, 2010
*------------------------------------------------------------
*/

double isat_lerror(gsl_vector *Rphi,gsl_vector *Rphi0,
		gsl_matrix *A,gsl_vector *phi,gsl_vector *phi0)
{
    double eps;
    gsl_vector *Rlphi = NULL;

    /* memory allocation for Rlphi */
    Rlphi = gsl_vector_calloc(phi->size);

    /* Rlphi := R(phi0) + A*(phi-phi0) */
    linear_approx(phi,phi0,Rphi0,A,Rlphi);

    /* Rlphi := Rl(phi) - R(phi) */
    gsl_vector_sub(Rlphi,Rphi);

    /* eps := 2-norm(Rl(phi) - R(phi)) */
    eps = gsl_blas_dnrm2(Rlphi);

    /* releasing allocated memory */
    gsl_vector_free(Rlphi);
    Rlphi = NULL;

    return eps;
}
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   isat4
*
*   This function executes the 4-th version of the
*   in situ adaptive tabulation algorithm.
*
*
*   Input:
*   isat      - isat workspace
*   thrm_data - thermochemistry workspace
*   etol      - error tolerance
*   n0        - factor to multiply by unit roundoff
*   t0        - initial time
*   delta_t   - time step
*   phi       - query composition
*   A         - mapping gradient matrix
*   L         - EOA Cholesky matrix
*
*   Output:
*   Rphi      - reaction mapping
*   success or error
*
*   last update: Nov 15, 2010
*------------------------------------------------------------
*/

int isat4(isat_wrk *isat,void *thrm_data,void *cvode_mem,
	    double etol,double n0,double t0,double delta_t,
    gsl_vector *phi,gsl_matrix *A,gsl_matrix *L,gsl_vector *Rphi)
{
    clock_t cpu_start = clock();

    /* ISAT first step */
    if( isat->lf == 0 )
    {
        int flag;
        bst_leaf *first_leaf = NULL;

	/* memory allocation for first_leaf */
        first_leaf = bst_leaf_alloc();
	if ( first_leaf == NULL )
                return GSL_ENOMEM;

        /* performing direct integration */
        flag = odesolver_reinit(thrm_eqs,thrm_data,t0,
                                        ATOL,RTOL,phi,cvode_mem);
        if ( flag != GSL_SUCCESS )
            return flag;

        flag = odesolver(cvode_mem,delta_t,Rphi);
        if ( flag != GSL_SUCCESS )
            return flag;

        /* computing the mapping gradient matrix */
        flag = gradient(thrm_eqs,thrm_data,cvode_mem,t0,delta_t,
                                                ATOL,RTOL,phi,Rphi,A);
        if ( flag != GSL_SUCCESS )
            return flag;

        /* computing EOA Cholesky matrix */
        isat_eoa_mtrx(A,etol,n0,L);

        /* setting first leaf elements */
        bst_leaf_set(phi,Rphi,A,L,first_leaf);
        isat->root->r_leaf = first_leaf;

        /* updating leaves and height counters */
        isat->lf++;
        isat->hgt = bst_height(isat->root);

        return GSL_SUCCESS;
    }

    /* ISAT second and following steps */
    int side, flag;
    double dot         = 0.0;
    bst_leaf *end_leaf = NULL;
    bst_node *end_node = NULL;

    /* searching for the near composition in binary search tree */
    if( isat->lf > 1 )
	side = bst_search(isat->root,phi,&end_node,&end_leaf);
    else
    {
        end_node = isat->root;
        end_leaf = isat->root->r_leaf;
        side     = RIGHT;
    }

    /* checking if the near composition is inside EOA */
    dot = ell_pt_in(phi,end_leaf->phi,end_leaf->L);
    if( gsl_fcmp(dot,1.0,ATOL) < 1 )
    {
        /* computing linear approximation */
        linear_approx(phi,end_leaf->phi,end_leaf->Rphi,end_leaf->A,Rphi);

        /* updating counters */
        isat->rtv++;
	isat->time_rtv += clock() - cpu_start;

        return GSL_SUCCESS;
    }
    else
    {
        double lerror = 0.0;

        /* performing direct integration */
        flag = odesolver_reinit(thrm_eqs,thrm_data,t0,
                                        ATOL,RTOL,phi,cvode_mem);
        if ( flag != GSL_SUCCESS )
            return flag;

        flag = odesolver(cvode_mem,delta_t,Rphi);
        if ( flag != GSL_SUCCESS )
            return flag;

        /* computing ISAT local error */
        lerror = isat_lerror(Rphi,end_leaf->Rphi,
                            end_leaf->A,phi,end_leaf->phi);

        /* checking if lerror is greater than etol */
        if( gsl_fcmp(lerror,etol,ATOL) < 1 )
        {
            /* growing the EOA */
            ell_pt_modify(phi,end_leaf->phi,end_leaf->L);

            /* updating counters */
            isat->grw++;
	    isat->time_grw += clock() - cpu_start;

            return GSL_SUCCESS;
        }
        else
        {
	    /* checking if the maximum # of leaves was excced */
	    if( isat->lf > isat->max_lf )
	    {
	        /*  updating counter */
	        isat->dev++;
		isat->time_dev +=  clock() - cpu_start;

		return GSL_SUCCESS;
	    }

            bst_leaf *new_leaf = NULL;

	    /* memory allocation */
            new_leaf = bst_leaf_alloc();
            if ( new_leaf == NULL )
                return GSL_ENOMEM;

            /* computing the mapping gradient matrix */
            flag = gradient(thrm_eqs,thrm_data,cvode_mem,t0,delta_t,
                                                ATOL,RTOL,phi,Rphi,A);
            if ( flag != GSL_SUCCESS )
		return flag;

            /* computing EOA Cholesky matrix */
            isat_eoa_mtrx(A,etol,n0,L);

            /* setting new_leaf elements */
            bst_leaf_set(phi,Rphi,A,L,new_leaf);

            /* cheking if binary search tree has more than one leaf */
            if( isat->lf > 1 )
            {
                bst_node *new_node = NULL;

		/* memory allocation for new node */
                new_node = bst_node_alloc();
                if ( new_node == NULL )
                    return GSL_ENOMEM;

		/* setting the leaves of the new node */
                bst_node_set(end_leaf,new_leaf,new_node);

		/* adding the new node to the tree */
                bst_node_add(side,end_node,new_node);
            }
	    else
                bst_node_set(end_leaf,new_leaf,end_node);

            /* updating counters */
            isat->add++;
            isat->lf++;
            isat->nd++;
            isat->hgt = bst_height(isat->root);
	    isat->time_add += clock() - cpu_start;

            return GSL_SUCCESS;
        }
    }
}
/*------------------------------------------------------------*/

max_stars_repo_name: americocunhajr/CRFlowLib

max_issues_repo_name: americocunhajr/CRFlowLib

max_forks_repo_name: americocunhajr/CRFlowLib


text:
/*
* -----------------------------------------------------------------
*  util_lib.h
*  Utilities Library
*  Version: 2.0
*  Last Update: Nov 3, 2019
*
*  Programmer: Americo Barbosa da Cunha Junior
*              americo.cunhajr@gmail.com
* -----------------------------------------------------------------
*  Copyright (c) 2010-2019, Americo Barbosa da Cunha Junior
*  All rights reserved.
* -----------------------------------------------------------------
*  This is the header file for UTIL_LIB module, a computational
*  library with routines for miscelaneous applications.
* -----------------------------------------------------------------
*/

#ifndef __UTIL_LIB_H__
#define __UTIL_LIB_H__

#include <time.h>

#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_rng.h>

/*
*------------------------------------------------------------
* function prototypes
*------------------------------------------------------------
*/

int gsl_fprint_vector(FILE* file,
                      gsl_vector *v);

int gsl_print_vector(gsl_vector *v);

int gsl_print_matrix_line(unsigned int i,
                          gsl_matrix *M);

int gsl_print_matrix_column(unsigned int j,
                            gsl_matrix *M);

int gsl_print_matrix(gsl_matrix *M);

int gsl_rand_vector(gsl_rng *r,
                    gsl_vector *v);

int gsl_rand_matrix(gsl_rng *r,
                    gsl_matrix *M);

int gsl_diagonal_matrix(gsl_vector *v,
                        gsl_matrix *D);

void util_time_used(clock_t cpu_start,
                    clock_t cpu_end,
                    time_t wall_start,
                    time_t wall_end);

#endif /* __UTIL_LIB_H__ */


text:
#define _BSD_SOURCE
#include <bcm2835.h>
#include <stdlib.h>
#include <unistd.h>
#include <stdio.h>
#include <math.h>
#include "csvparser.h"
#include <time.h>
#include "mailbox.h"
#include "gpu_fft.h"
#include "spi.h"
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_complex_math.h>

gsl_matrix_complex *compute_fft(int *points, int length) { //length should be specified as length = log2(FFT_length)
    unsigned t[2];
    struct GPU_FFT_COMPLEX *dataIn,*dataOut;
    struct GPU_FFT *fftinfo;

    int mb = mbox_open();
    gpu_fft_prepare(mb,length,GPU_FFT_FWD,1,&fftinfo);

    dataIn = fftinfo->in;
    dataOut = fftinfo->out;

    // Fill the fft data input array with signal
    int i;
    int FFT_length = 1 << length;
    for (i = 0; i < FFT_length; i++) {
	dataIn[i].re = points[i];
	dataIn[i].im = 0;
    }

    usleep(1); // Yield to OS
    clock_t t1 = clock();
    gpu_fft_execute(fftinfo);
    clock_t t2 = clock();
    float diff = (1000000*(float)(t2 - t1) / CLOCKS_PER_SEC );
    printf("Computed FFT of length %d in %f microseconds\n",FFT_length,diff);

    gsl_matrix_complex *matrix = gsl_matrix_complex_alloc(FFT_length,1); //Initialize complex fft vector
    // Fill it with complex values
    for (i = 0; i < FFT_length; i++) {
       gsl_complex complex = gsl_complex_rect(dataOut[i].re,dataOut[i].im);
       gsl_matrix_complex_set(matrix, i, 0,complex); //Fill matrix rows
    }
    gpu_fft_release(fftinfo);
    return matrix;
}

/*
scale_matrix function

matrix - comlex fft matrix to scale

This function scales matrix by it's size.

I use it to normalize fft vector.
*/
void scale_matrix(gsl_matrix_complex *matrix) {
     size_t size = matrix->size1;
     double scale = 1.0/size;
     gsl_complex complex_scale = gsl_complex_rect(scale,0); // This is just a scale represented by a complex number
     gsl_matrix_complex_scale(matrix,complex_scale); //Scale matrix by the length of the signal
}

/*
abs_matrix function

matrix - comlex fft matrix

This function takes an absolute valus of the complex matrix and returns a real matrix as a result.

*/
gsl_matrix *abs_matrix(gsl_matrix_complex *matrix) {
     size_t size = matrix->size1;
     gsl_matrix *absolute = gsl_matrix_alloc(size,1);
     for (int i = 0; i < size; i++) {
         gsl_complex complex = gsl_matrix_complex_get(matrix,i,0);
	 double abs = gsl_complex_abs(complex);
	 gsl_matrix_set(absolute,i,0,abs);
     }
     return absolute;
}

/*
half_matrix function

matrix - real values matrix

This function "chops" matrix in half to get the positive frequency spectrum.
*/
gsl_matrix *half_matrix(gsl_matrix *matrix) {
     size_t size = matrix->size1;
     size_t half_size = (size / 2)+1;
     gsl_matrix *half = gsl_matrix_alloc(half_size,1); //Take only half of the absolute of ff
     for (int i = 0; i < half_size; i++) {
         double value = gsl_matrix_get(matrix,i,0);
	 gsl_matrix_set(half,i,0,value);
     }
     return half;
}

/*
save_waveform_to_file function

data - signal
length - signal length
filename - name of the file to be saved to

Saves the signal to the file.
*/
void save_waveform_to_file(int *data, int length, char *filename) {
    FILE *f = fopen(filename,"w");
    printf("Length: %d\n",length);
    for (int i = 0; i < length; i++) {
	fprintf(f,"%d,\n",data[i]);
    }
    fclose(f);
}

/*
print_result_to_file function

result - vector
filename - name of the file to be saved to

Saves the vector to the file.
*/
void print_result_to_file(gsl_matrix *result, char *filename) {
        FILE *f = fopen(filename,"w");
	int size = result->size1;
	for (int i = 0; i < size; i++) {
            double value = gsl_matrix_get(result,i,0);
	    fprintf(f,"%d,\n",value);
	}
	fclose(f);
}

/*
print_fft_to_file function

result - complex fft vector

Saves the vector to the file.
*/
void print_fft_to_file(gsl_matrix_complex *result) {
        FILE *real = fopen("real.csv","w");
        FILE *imag = fopen("imag.csv","w");

	int size = result->size1;
	for (int i = 0; i < size; i++) {
            gsl_complex value = gsl_matrix_complex_get(result,i,0);
	    fprintf(real,"%f,\n",value.dat[0]);
	    fprintf(imag,"%f,\n",value.dat[1]);
	}
	fclose(real);
	fclose(imag);
}

void display_harmonics(gsl_matrix *fft_vector) {
	double value_60 = gsl_matrix_get(fft_vector,60,0);
	double amplitude = value_60 * (3.72/4096);
	printf("Amplitude of 60Hz is: %f\n",amplitude);
}

/*
extract_50_harmonics function

fft_vector - complex fft vector

Exracts 50 fundamental hormonics
60,120,180, etc

Returns the harmonics array of size 50.
*/
int *extract_50_harmonics(gsl_matrix *fft_vector) {
	int *harmonics = (int *)malloc(sizeof(int)*50);;
	for (int i = 1; i < 51; i++) {
		double value = gsl_matrix_get(fft_vector,i*60,0);
		int amplitude = (int)(value * (1000.0/4096) / 1.14);
		harmonics[i-1] = amplitude;
	}
	return harmonics;
}

/*
perform_fft_analysis function

points - signal

performs fft analysis
*/
void perform_fft_analysis(int *points) {
    clock_t t1 = clock();
    int length = 1 << WINDOW_SIZE; // Sampling frequency
    gsl_matrix_complex *fft_matrix = compute_fft(points,WINDOW_SIZE); // performs fft transformation of the signal
    scale_matrix(fft_matrix); // normalizes the resulting fft vector
    gsl_matrix *absolute = abs_matrix(fft_matrix); // takes an absolute value of the fft vector
    gsl_matrix *half = half_matrix(absolute); // and chops of the negative part.
    gsl_matrix_scale(half,2); // Scales by two
    clock_t t2 = clock();
    float diff = ((float)(t2 - t1) / CLOCKS_PER_SEC );
    printf("FFT analysis took: %f seconds\n",diff);

    int *harmonics = extract_50_harmonics(half); // extracts harmonics
    for (int i = 1; i < 51; i++) { // diplays harmonics aplitudes.
	   printf("Amplitude at %d Hz is %d mV\n",i*60,harmonics[i-1]);
    }

    // frees the resources.
    gsl_matrix_free(half);
    gsl_matrix_free(absolute);
    gsl_matrix_complex_free(fft_matrix);
    free(harmonics);
}

/*
rms function

signal - signal
frequency - signal frequency

Computes rms value of the given signal.
*/
int rms(int *signal, int frequency) {
    clock_t t1 = clock();
    int sum = 0;
    for (int i = 0; i < frequency; i++) {
       int data = signal[i];
	sum += data*data;
    }
    sum /= frequency;
    double root = sqrt((double)sum);
    clock_t t2 = clock();
    float diff = ((float)(t2 - t1) / CLOCKS_PER_SEC );
    printf("RMS analysis took: %f seconds\n",diff);
    printf("RMS value: %f\n",root);
    return (int)root;
}

int main(int argc, char *argv[]) {
    if (argc < 2) {
        printf("Usage %s channel\n",argv[0]);
	    return 1;
    }

    int channels[8] = {0};

     // Extract which channels to sample from the command arguments.
    for (int arg_index = 0; arg_index < argc - 1; arg_index++) {
         int channel = atoi(argv[arg_index + 1]);
	     printf("Sampling channel %d\n",channel);
         channels[channel] = 1;
    }

    if (init()) {
        printf("Init failed.\n");
	    return 1;
    } else {
	    printf("Succesfull init of bcm2835.\n");
    }
    printf("Starting bcm2835 setup...\n");
    setup();
    printf("bcm2835 setup finished.\n");

    for (int i = 0; i < 1; i++) { // Specify sampling duration here! If we want to sample indefinetely, just write while(1) {}
	    int frequency = 1 << WINDOW_SIZE;
        sample *s = sample_data(channels,frequency);
	    perform_fft_analysis(s->channel_data[1]);
	    sample_deallocate(s);
    }

    printf("Ending spi connection ...\n");
    end();
    printf("bcm2835 has ended.\n");
}

max_stars_repo_name: seadsystem/FFT-Harmonic-Extraction

max_issues_repo_name: seadsystem/FFT-Harmonic-Extraction

max_forks_repo_name: seadsystem/FFT-Harmonic-Extraction


text:
#ifndef FUSED_LASSO_E_H
#define FUSED_LASSO_E_H

#include "Matrix.h"
#include "LADMM.h"
#include <string>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <stdio.h>      /* printf */
#include <time.h>
#include <fstream>
#include <algorithm>
#include <iomanip>
#include <ctime>
#include <math.h>

//This class solves problem: min_x f(x)+ g(x)+ h(Mx) by LADMM;
// where g(x)=\frac{lambda2}{2}\|x\|_2^2 +lambda1\|x\|_1;
// f(x)= 1/2*\|Ax- b\|_2^2= \sum_i phi_i(A_ix);
// f_i(x)= 0.5*(x- b_i)^2, h_i(x)= \|x\|, Mx= (x^1- x^2, ..., x_n- x_1).

template<typename L, typename D>
class Fused_lasso_e: public LADMM<L, D>
{
private:

D lambda1;

D lambda2;

D lambda3;

Matrix<L,D> my_M;

Matrix<L,D> my_A;

protected:

public:

  Fused_lasso_e(const char* Matrix_file,const char* Matrix_file2, D val_lambda1, D val_lambda2, D val_lambda3)
  :LADMM<L,D>(),my_A(Matrix_file), my_M(Matrix_file2)
  {
    lambda1=val_lambda1;
    lambda2=val_lambda2;
    lambda3=val_lambda3;
    this->lambda_f= 1;
  }

Fused_lasso_e(const char* Matrix_file,D val_lambda1, D val_lambda2, D val_lambda3)
  :LADMM<L,D>(),my_A(Matrix_file)
  {
    lambda1=val_lambda1;
    lambda2=val_lambda2;
    lambda3=val_lambda3;
    this->lambda_f= 1;
    my_M.construct_fused_matrix(my_A);
  }

  L get_n(){return my_M.nfeatures;}
  L get_m(){return my_M.nsamples;}

  inline D gradient_of_f_j(D x1, L i){
      return x1- my_A.b[i];
  }

  inline D value_of_g_j(D x, L i){
        return lambda2*x*x/2+lambda1*fabs(x);
      }

  inline D value_of_f_j(D x1, L i){
  		return 0.5*(x1- my_A.b[i])*(x1- my_A.b[i]);
	}

  inline D value_of_h_j(D x, L i){
  	return lambda3*fabs(x);
  }

  inline D prox_of_h_j(D x1,D x2, L i){
  	D newx;
  	if (x1- lambda3/x2> my_M.b[i]){
  		newx= x1- lambda3/x2;
	  }
	  else if(x1+ lambda3/x2< my_M.b[i]){
	  	newx= x1+ lambda3/x2;
	  }
	  else{
	  	newx= my_M.b[i];
	  }
	  return newx;
  }

  inline D prox_of_g_j(D x1,D x2, L i){
        D new_x;
        if(x1*x2> lambda1)
         new_x=(x1*x2- lambda1)/(x2+lambda2);
        else if(x1*x2< -lambda1)
         new_x=(x1*x2+ lambda1)/(x2+ lambda2);
        else
         new_x=0;
        return new_x;
      }

  inline void set_matrix_M(){
    this->data_M=my_M;
  }

  inline void set_matrix_A(){
    this->data_A=my_A;
  }
/*
  inline void set_matrix_A(){
    this->data_A.nsamples=0;
    this->data_A.nfeatures=this->data_M.nfeatures;
    this->data_A.nnz=0;
    this->data_A.ptr.resize(1,0);
    this->data_A.ptr_t.resize(this->data_M.nfeatures+1,0);
  }
*/

  void LADMM_solver(D beta_0,  D rho,vector<D> & x0,vector<D> & y0, vector<D> & lambda0, L max_nb_outer, L p_N_1, string filename1, D time){
    this->ADMM_solve_with_Linear(beta_0, rho,x0,y0, lambda0, max_nb_outer,  p_N_1, filename1, time);
}

};

#endif


text:
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <float.h>
#include <math.h>
#include <time.h>
#include <sys/time.h>
#include <cblas.h>

/* Your function must have the following signature: */

void sgemm( int m, int n, float *A, float *C );

/* The benchmarking program */

int main( int argc, char **argv )
{
	srand(time(NULL));

	if(argc != 3)
	{
		printf("Incorrect number of arguments!\n");
		printf("Usage: ./bench-test M N\n");
		exit(0);
	}

	int m  = atoi(argv[1]);
	int n  = atoi(argv[2]);

	// Change these values to modify range checking
	if(m < 1000 || m > 10000 || n < 32 || n > 10000)
	{
		printf("Dimensions provided are out of the range of matrices we will be testing!  Change tester.c to modify range checking.\n");
		exit(0);
	}

	/* Allocate and fill 2 random matrices A, C */
	float *A = (float*) malloc( m * n * sizeof(float) );
	float *C = (float*) malloc( m * m * sizeof(float) );

	for( int i = 0; i < m*n; i++ ) A[i] = 2 * drand48() - 1;
	for( int i = 0; i < m*m; i++ ) C[i] = 2 * drand48() - 1;

	/* measure Gflop/s rate; time a sufficiently long sequence of calls to eliminate noise */
	double Gflop_s, seconds = -1.0;
	for( int n_iterations = 1; seconds < 0.1; n_iterations *= 2 )
	{
		/* warm-up */
		sgemm( m, n, A, C );

		/* measure time */
		struct timeval start, end;
		gettimeofday( &start, NULL );
		for( int i = 0; i < n_iterations; i++ )
			sgemm( m,n, A, C );
		gettimeofday( &end, NULL );
		seconds = (end.tv_sec - start.tv_sec) + 1.0e-6 * (end.tv_usec - start.tv_usec);

		/* compute Gflop/s rate */
		Gflop_s = 2e-9 * n_iterations * m * m * n / seconds;
	}

	printf( "%d by %d matrix \t %g Gflop/s\n", m, n, Gflop_s );

	/* Ensure that error does not exceed the theoretical error bound */

	/* Set initial C to 0 and do matrix multiply of A*B */
	memset( C, 0, sizeof( float ) * m * m );
	sgemm( m,n, A, C );

	/* Subtract A*B from C using standard sgemm (note that this should be 0 to within machine roundoff) */
	cblas_sgemm( CblasColMajor,CblasNoTrans,CblasTrans, m,m,n, -1, A,m, A,m, 1, C,m );

	/* Subtract the maximum allowed roundoff from each element of C */
	for( int i = 0; i < m*n; i++ ) A[i] = fabs( A[i] );
	for( int i = 0; i < m*m; i++ ) C[i] = fabs( C[i] );
	cblas_sgemm( CblasColMajor,CblasNoTrans,CblasTrans, m,m,n, -3.0*FLT_EPSILON*n, A,m, A,m, 1, C,m );

	/* After this test if any element in C is still positive something went wrong in square_sgemm */
	for( int i = 0; i < m * m; i++ )
	  	if( C[i] > 0 ) {
			printf( "FAILURE: error in matrix multiply exceeds an acceptable margin\n" );
		return -1;
	}

	/* release memory */
	free( C );
	free( A );
	return 0;
}

max_stars_repo_name: jasonaibrahim/AcademicWorks

max_issues_repo_name: jasonaibrahim/AcademicWorks

max_forks_repo_name: jasonaibrahim/AcademicWorks


text:
#include <math.h>
#include <stdio.h>
#include <mpfr.h>
#include <petsc.h>
#include <stdlib.h>
#include <string.h>

#include "ellipsoid/ellipsoid.h"

int testCoordinateTransform()
{
  //Verify that the ellipsoidal to Cartesian transformation is a bijection
  //the tolerance is 10 significant figures
  EllipsoidalSystem e;
  initEllipsoidalSystem(&e, 3.0, 2.0, 1.0, 32);

  double relErrorX, relErrorY, relErrorZ;

  //loop over octants
  for(int sx=-1; sx<=1; sx+=2) { //sx = [-1,1]
    for(int sy=-1; sy<=1; sy+=2) { //sy = [-1,1]
      for(int sz=-1; sz<=1; sz+=2) { //sz = [-1,1]
	//loop over brick in cartesian space
	for(double x=4; x<5+tol; x+=1.5) {
	  for(double y=4; y<5+tol; y+=1.5) {
	    for(double z=4; z<5+tol; z+=1.5) {
	      Point temp = { .x1 = sx*x, .x2 = sy*y, .x3 = sz*z };
	      cartesianToEllipsoidal(&e, &temp);
	      ellipsoidToCartesian(&e, &temp);
	      relErrorX = fabs(sx*x - temp.x1)/x;
	      relErrorY = fabs(sy*y - temp.x2)/y;
	      relErrorZ = fabs(sz*z - temp.x3)/z;
	      //check error
	      if( relErrorX > tol ||
		  relErrorY > tol ||
		  relErrorZ > tol) {
		printf("testCoordinateTransform failed\n");
		printf("relative error(x): %16.16f\n", relErrorX);
		printf("relative error(y): %16.16f\n", relErrorY);
		printf("relative error(z): %16.16f\n", relErrorZ);
		return 0;
	      }
	    }
	  }
	}
      }
    }
  }
  return 1;
}

int testLameType()
{

  static char formula1[18] = "K^0_1 L^0_1 M^0_1";
  static char formula2[30] = "K^0_2 K^1_2 L^0_2 M^0_2 N^0_2";
  static char formula3[66] = "K^0_5 K^1_5 K^2_5 L^0_5 L^1_5 L^2_5 M^0_5 M^1_5 M^2_5 N^0_5 N^1_5";

  char **checks = (char**) malloc(sizeof(char*)*3);
  checks[0] = (char*) malloc(sizeof(char)*17);
  checks[1] = (char*) malloc(sizeof(char)*29);
  checks[2] = (char*) malloc(sizeof(char)*65);

  EllipsoidalSystem e;
  initEllipsoidalSystem(&e, 3.0, 2.0, 1.0, 32);
  int index[3] = { 1, 2, 5 };
  int k;
  int n;
  for(int i=0; i<3; ++i) {
    n = index[i];
    for(int p=0; p < 2*n + 1; ++p) {
      k = p*6;
      checks[i][k]   = getLameTypeT(n, p);
      checks[i][k+1] = '^';
      checks[i][k+2] = getLameTypeTp(n, p) + '0';
      checks[i][k+3] = '_';
      checks[i][k+4] = n + '0';
      if( p != 2*n )
	checks[i][k+5] = ' ';
    }
  }
  //check first string
  for(int i=0; i < 17; ++i) {
    if(checks[0][i] != formula1[i]) {
      return 0;
    }
  }
  //check second string
  for(int i=0; i < 29; ++i) {
    if(checks[1][i] != formula2[i]) {
      return 0;
    }
  }
  //check first string
  for(int i=0; i < 65; ++i) {
    if(checks[2][i] != formula3[i]) {
      return 0;
    }
  }
  free(checks[0]);
  free(checks[1]);
  free(checks[2]);
  free(checks);
  return 1;
}

int testI()
{
  //Using results from Dassios, test the computation of I^p_n for n = 0, 1, 2
  EllipsoidalSystem e;
  initEllipsoidalSystem(&e, 3.0, 2.0, 1.0, 32);
  double a = e.a;
  double b = e.b;
  double c = e.c;
  double h = e.h;
  double k = e.k;
  double l = a;
  double rootSum = sqrt((a*a*a*a - b*b*c*c) + (b*b*b*b - a*a*c*c) + (c*c*c*c - a*a*b*b));
  double DassiosLambda = (1/3.)*(a*a + b*b + c*c) + (1/3.)*rootSum;
  double DassiosLambdaPrime = (1/3.)*(a*a + b*b + c*c) - (1/3.)*rootSum;

  double errors[18];
  int derror[18] = { 1, 2, 3, 3, 4, 4, 5, 5, 6, 6, 6, 7, 8, 9, 10, 11, 12, 13 };
  int ecount = 0;
  //Dassios D1 error
  double lhs = 3 * (DassiosLambda + DassiosLambdaPrime);
  double rhs = 2 * (a*a + b*b + c*c);
  errors[ecount] = fabs(lhs - rhs); ecount++;
  //Dassios D2 error
  lhs = 3 * DassiosLambda * DassiosLambdaPrime;
  rhs = a*a*b*b + a*a*c*c + b*b*c*c;
  errors[ecount] = fabs(lhs - rhs); ecount ++;
  //Dassios D3 error
  lhs = ((-1)*(b*b - c*c)*(DassiosLambda - a*a)) + (k*k*(DassiosLambda - b*b)) + (-h*h*(DassiosLambda - c*c));
  rhs = ((-1)*(b*b - c*c)*(DassiosLambdaPrime - a*a)) + (k*k*(DassiosLambdaPrime - b*b)) + (-h*h*(DassiosLambdaPrime - c*c));
  errors[ecount] = fabs(lhs-0); ecount++;
  errors[ecount ] = fabs(rhs-0); ecount++;
  //Dassios D4 error
  lhs = (-a*a*(b*b - c*c)*(DassiosLambda - a*a)) + (b*b*k*k*(DassiosLambda - b*b)) + (-c*c*h*h*(DassiosLambda - c*c));
  rhs = (-a*a*(b*b - c*c)*(DassiosLambdaPrime - a*a)) + (b*b*k*k*(DassiosLambdaPrime - b*b)) + (-c*c*h*h*(DassiosLambdaPrime - c*c));
  double analytical = h*h*k*k*(b*b - c*c);
  errors[ecount] = fabs(lhs - analytical); ecount++;
  errors[ecount] = fabs(rhs - analytical); ecount++;
  //Dassios D5 error
  double alpha[3] = { a, b, c };
  lhs = 0;
  rhs = 0;
  for(int i=0; i<3; ++i) {
    lhs += alpha[i]*alpha[i]/(alpha[i]*alpha[i] - DassiosLambda);
    rhs += alpha[i]*alpha[i]/(alpha[i]*alpha[i] - DassiosLambdaPrime);
  }
  errors[ecount] = fabs(lhs - 3); ecount++;
  errors[ecount] = fabs(rhs - 3); ecount++;
  //Dassios D6 error
  double anal1 = h*h*k*k*(b*b - c*c);
  double num1  = 3*(b*b - c*c)*(DassiosLambda - a*a)*(DassiosLambdaPrime - a*a);
  double anal2 = -h*h*k*k*(b*b - c*c);
  double num2  = 3*k*k*(DassiosLambda - b*b)*(DassiosLambdaPrime - b*b);
  double anal3 = h*h*k*k*(b*b - c*c);
  double num3  = 3*h*h*(DassiosLambda - c*c)*(DassiosLambdaPrime - c*c);
  errors[ecount] = fabs(num1 - anal1); ecount++;
  errors[ecount] = fabs(num2 - anal2); ecount++;
  errors[ecount] = fabs(num3 - anal3); ecount++;

  //Calculate I^p_n
  int nmax = 2;
  int Isize = 0;
  for(int n=0; n < nmax+1; ++n) {
    for(int p=0; p < 2*n+1; ++p) {
      Isize++;
    }
  }
  double *Ival = (double*) malloc(sizeof(double)*Isize);
  int i=0;
  for(int n=0; n < nmax+1; ++n) {
    for(int p=0; p < 2*n+1; ++p) {
      calcI(&e, n, p, l, 1, 1, Ival+i); //Ival[i] = calcI(&e, n, p, l, 1, 1);
      i++;
    }
  }

  //Dassios D7: their h_3 = our h and their h_2 = our k
  double analyticalSumTest1 = 1.0/(l*sqrt(l*l - h*h)*sqrt(l*l - k*k));
  double numericalSumTest1 = 0;
  for(int k=1; k<4; ++k) {
    numericalSumTest1 += Ival[k];
  }
  errors[ecount] = fabs(analyticalSumTest1 - numericalSumTest1); ecount++;

  //Dassios D8: their alpha1 = our a, their alpha2 = our b, their alpha3 = our c
  double analyticalSumTest2 = Ival[0] - (l*l - a*a)/(l*sqrt(l*l - h*h)*sqrt(l*l - k*k));
  double numericalSumTest2 = a*a*Ival[1] + b*b*Ival[2] + c*c*Ival[3];
  errors[ecount] = fabs(analyticalSumTest2 - numericalSumTest2); ecount++;
  //Dassios D9: their Lambda = our DassiosLambda, their LambdaPrime = our DassiosLambdaPrime
  double analyticalSumTest3 = 1.0/(2.0*(DassiosLambda - a*a + l*l)*l*sqrt(l*l - h*h)*sqrt(l*l - k*k)) - (1.0/2.0)*(Ival[1]/(DassiosLambda - a*a) + Ival[2]/(DassiosLambda - b*b) + Ival[3]/(DassiosLambda - c*c));
  double numericalSumTest3 = Ival[4];
  errors[ecount] = fabs(analyticalSumTest3 - numericalSumTest3); ecount++;
  //Dassios D10: their Lambda = our DassiosLambda, their LambdaPrime = our DassiosLambdaPrime
  double analyticalSumTest4 = 1.0/(2.0*(DassiosLambdaPrime - a*a + l*l)*l*sqrt(l*l - h*h)*sqrt(l*l - k*k)) - (1.0/2.0)*(Ival[1]/(DassiosLambdaPrime - a*a) + Ival[2]/(DassiosLambdaPrime - b*b) + Ival[3]/(DassiosLambdaPrime - c*c));
  double numericalSumTest4 = Ival[5];
  errors[ecount] = fabs(analyticalSumTest4 - numericalSumTest4); ecount++;
  //Dassios D11: their h1*h1 = our (b*b - c*c)
  double analyticalSumTest5 = (1.0/(h*h)) * (Ival[2] - Ival[1]);
  double numericalSumTest5 = Ival[6];
  errors[ecount] = fabs(analyticalSumTest5 - numericalSumTest5); ecount++;
  //Dassios D12: their h1*h1 = our (b*b - c*c)
  double analyticalSumTest6 = (1.0/(k*k)) * (Ival[3] - Ival[1]);
  double numericalSumTest6 = Ival[7];
  errors[ecount] = fabs(analyticalSumTest6 - numericalSumTest6); ecount++;
  //Dassios D13: their h1*h1 = our (b*b - c*c)
  double analyticalSumTest7 = (1.0/(b*b - c*c)) * (Ival[3] - Ival[2]);
  double numericalSumTest7 = Ival[8];
  errors[ecount] = fabs(analyticalSumTest7 - numericalSumTest7); ecount++;

  //checking errors
  for(int i=0; i < 18; ++i) {
    if(errors[i] > tol) {
      printf("Dassios D%d failed with error %8.8e\n", derror[i], errors[i]);
      free(Ival);
      return 0;
    }
  }

  free(Ival);
  return 1;
}

int testNormalization()
{
  //Check calculated \gamma^p_n against analytic results for n = 0, 1, 2
  double a = 3.0;
  double b = 2.0;
  double c = 1.0;
  EllipsoidalSystem e;
  initEllipsoidalSystem(&e, a, b, c, 32);
  //Dassios (B14)
  double firstTerm = (a*a + b*b + c*c)/3.0;
  double secondTerm = sqrt((a*a*a*a - b*b*c*c) + (b*b*b*b - a*a*c*c) + (c*c*c*c - a*a*b*b))/3.0;
  double LambdaD = firstTerm + secondTerm;
  double LambdaDprime = firstTerm - secondTerm;
  //Dassios (B16)-(B20)
  double hx = sqrt(b*b - c*c);
  double hy = e.k;
  double hz = e.h;
  double analytic[9] = { 4*PETSC_PI,
			 4*PETSC_PI/3 * hy*hy*hz*hz,
			 4*PETSC_PI/3 * hx*hx*hz*hz,
			 4*PETSC_PI/3 * hx*hx*hy*hy,
			 -8*PETSC_PI/5 * (LambdaD - LambdaDprime)*(LambdaD - a*a)*(LambdaD - b*b)*(LambdaD - c*c),
			 8*PETSC_PI/5 * (LambdaD - LambdaDprime)*(LambdaDprime - a*a)*(LambdaDprime - b*b)*(LambdaDprime - c*c),
			 4*PETSC_PI/15 * hx*hx*hy*hy*hz*hz*hz*hz,
			 4*PETSC_PI/15 * hx*hx*hy*hy*hy*hy*hz*hz,
			 4*PETSC_PI/15 * hx*hx*hx*hx*hy*hy*hz*hz };
  printf("analytic[0] = %4.4e\n", analytic[0]);
  printf("analytic[1] = %4.4e\n", analytic[1]);
  double value, estimate, error;
  int counter = 0;
  for(int n=0; n<3; ++n) {
    for(int p=0; p<2*n+1; ++p) {
      value = analytic[counter];
      calcNormalization(&e, n, p, &estimate);
      error = estimate - value;
      if(fabs(error) > tol) {
	printf("testnormalization failed iteration %d at (n,p) = (%d,%d) with error %8.8e\n", counter, n, p, fabs(error));
	return 0;
      }
      counter++;
    }
  }
  return 1;

}

void testfunc(mpfr_t *x, mpfr_t *val, FuncInfo4 *ctx)
{
  double a2 = (*ctx).a2;
  double b2 = (*ctx).b2;
  double c2 = (*ctx).c2;
  double botVar = (*ctx).botVar;

  mpfr_t temp, s;
  mpfr_init2(temp, mpfr_get_prec(*x));
  mpfr_init2(s   , mpfr_get_prec(*x));
  //change of variables
  //s = (1-x)/x
  mpfr_d_sub(s, 1.0, *x, MPFR_RNDZ);
  mpfr_div(s, s, *x, MPFR_RNDN);

  mpfr_add_d(temp, s, a2, MPFR_RNDN);
  mpfr_add_d(*val, s, b2, MPFR_RNDN);
  mpfr_mul(*val, *val, temp, MPFR_RNDN);
  mpfr_add_d(temp, s, c2, MPFR_RNDN);
  mpfr_mul(*val, *val, temp, MPFR_RNDN);
  mpfr_sqrt(*val, *val, MPFR_RNDN);

  mpfr_add_d(temp, s, botVar, MPFR_RNDN);
  mpfr_mul(*val, *val, temp, MPFR_RNDN);

  mpfr_d_div(*val, 1.0, *val, MPFR_RNDN);

  //multiply by dt = ds * (1+t)^2
  mpfr_add_d(temp, s, 1.0, MPFR_RNDN);
  mpfr_mul(temp, temp, temp, MPFR_RNDN);
  mpfr_mul(*val, *val, temp, MPFR_RNDN);
  //mpfr_clears(temp, s, NULL);
}

int testSurfaceOperatorEigenvalues()
{
  EllipsoidalSystem e;
  double a = 3.0;
  double b = 2.0;
  double c = 1.0;
  initEllipsoidalSystem(&e, a, b, c, 32);
  double l = a;
  int n = 1;
  double analytic[3];
  mpfr_t integral;

  mpfr_t mpfrzero, mpfrone;
  mpfr_inits(mpfrzero, mpfrone, integral, NULL);
  mpfr_set_d(mpfrzero, 0.0, MPFR_RNDN);
  mpfr_set_d(mpfrone, 1.0, MPFR_RNDN);

  //Integrals are from Ritter, need to transform integrals to finite interval
  FuncInfo4 ctx = { &e, a*a, b*b, c*c, a*a };
  integrateMPFR((PetscErrorCode (*)(mpfr_t*, mpfr_t*, void*)) testfunc, &e, mpfrzero, mpfrone, 14, &integral, &ctx);
  analytic[0] = (a*b*c * mpfr_get_d(integral, MPFR_RNDN) - 1.0)/2.0;
  ctx.botVar = b*b;
  integrateMPFR((PetscErrorCode (*)(mpfr_t*, mpfr_t*, void*)) testfunc, &e, mpfrzero, mpfrone, 14, &integral, &ctx);
  analytic[1] = (a*b*c * mpfr_get_d(integral, MPFR_RNDN) - 1.0)/2.0;
  ctx.botVar = c*c;
  integrateMPFR((PetscErrorCode (*)(mpfr_t*, mpfr_t*, void*)) testfunc, &e, mpfrzero, mpfrone, 14, &integral, &ctx);
  analytic[2] = (a*b*c * mpfr_get_d(integral, MPFR_RNDN) - 1.0)/2.0;

  double *ev = (double*) malloc(sizeof(double)*(3));
  for(int p=0; p < 3; ++p)
    calcSurfaceOperatorEigenvalues(&e, n, p, l, 1, 1, ev+p); //ev[p] = calcSurfaceOperatorEigenvalues(&e, n, p, l, 1, 1);
  for(int p=0; p < 3; ++p) {
    if(fabs(ev[p] - analytic[p]) > tol) {
      printf("analytic (%8.8e) doesn't match numerical (%8.8e), error = (%8.8e) for integral %d\n", analytic[p], ev[p], fabs(ev[p] - analytic[p]), p+1);
      return 0;
    }
  }
  mpfr_clears(mpfrzero, mpfrone, integral, NULL);
  free(ev);
  return 1;

}

int compareIntegration()
{

  //int n, p given in calcNormalization
  int n = 1;
  int p = 1;

  int digits = 32;

  //set up ellipsoidal system
  EllipsoidalSystem e;
  double xA = 3.0;
  double yB = 2.0;
  double zC = 1.0;
  initEllipsoidalSystem(&e, xA, yB, zC, digits);

  mpfr_t bound_a, bound_b;
  mpfr_inits2(4*digits, bound_a, bound_b, NULL);

  FuncInfo2 ctx1 = { .e = &e, .n = n, .p = p, .numeratorType = 0, .denomSign = 1 };

  mpfr_t a, b, c;
  mpfr_inits(a, b, c, NULL);
  double integralMPFR, integralMidpoint;

  integrateMPFR((PetscErrorCode (*)(mpfr_t*, mpfr_t*, void*)) normFunction1, &e, e.hp_h, e.hp_k, 14, &integralMPFR, &ctx1);
  integrateMidpoint((PetscErrorCode (*)(mpfr_t*, mpfr_t*, void*)) normFunction1, e.hp_h, e.hp_k, 10, &integralMidpoint, &ctx1);

  printf("MPFR integral is: %15.15f\n", integralMPFR);
  printf("Midpoint integral is: %15.15f\n", integralMidpoint);
  printf("oh wow\n");

  mpfr_clears(a, b, c, NULL);

  return 0;
}

int main(int argc, char **argv)
{
  mpfr_set_default_prec(4*64);
  /*
  EllipsoidalSystem e;
  initEllipsoidalSystem(&e, 3.0, 2.0, 1.0);
  printf("h: %15.15f\nh2: %15.15f\nk: %15.15f\nk2: %15.15f\n", e.h, e.h2, e.k, e.k2);
  int n = 1;
  int p = 1;
  double l = 1.0;

  //printf("calcLame(1, 1, 1.0) = %15.15f\n", calcLame(&e, 1, 1, 1.0));

  printf("calcNormalization(%d, %d, %15.15f) = %15.15f\n", n, p, l, calcNormalization(&e, n, p));
  */

  if(testCoordinateTransform())
    printf("testCoordinateTransform         passed\n");
  else
    printf("testCoordinateTransform         failed\n");
  if(testLameType())
    printf("testLameType                    passed\n");
  else
    printf("testLameType                    failed\n");
  if(testI())
    printf("testI                           passed\n");
  else
    printf("testI                           failed\n");
  if(testNormalization())
    printf("testNormalization               passed\n");
  else
    printf("testNormalization               failed\n");
  if(testSurfaceOperatorEigenvalues())
    printf("testSurfaceOperatorEigenvalues  passed\n");
  else
    printf("testSurfaceOperatorEigenvalues  failed\n");

  //compareIntegration();

  /*
  double h = .1;
  mpfr_t x, val;
  EllipsoidalSystem e;
  initEllipsoidalSystem(&e, 3.0, 2.0, 1.0);
  FuncInfo2 ctx1 = { .e = &e, .n = 1, .p = 2, .numeratorType = 0, .denomSign = -1 };
  mpfr_inits2(4*16, x, val, NULL);
  for(int k=0; k < 10; ++k) {
    mpfr_set_d(x, k*h, MPFR_RNDN);
    normFunction1(&x, &val, &ctx1);
    printf("normfunction1(%8.8f) = %15.15f\n", k*h, mpfr_get_d(val, MPFR_RNDN));
  }
  */
  return 0;
}

max_stars_repo_name: tom-klotz/ellipsoid-solvation

max_issues_repo_name: tom-klotz/ellipsoid-solvation

max_forks_repo_name: tom-klotz/ellipsoid-solvation


text:
//
//  herk.h
//  Linear Algebra Template Library
//
//  Created by Rodney James on 1/6/12.
//  Copyright (c) 2012 University of Colorado Denver. All rights reserved.
//

#ifndef _herk_h
#define _herk_h

/// @file herk.h Performs complex matrix-matrix multiplication operations resulting in Hermitian matrices.

#include <cctype>
#include "latl.h"

namespace LATL
{
   /// @brief Performs multiplcation of a complex matrix with its conjugate transpose.
   ///
   /// For complex matrix A, complex Hermitian matrix C, and real scalars alpha and beta
   ///
   ///        C := alpha*A*A'+beta*C  or  C := alpha*A'*A+beta*C
   /// is computed.
   /// @return 0 if success.
   /// @return -i if the ith argument is invalid.
   /// @tparam real_t Floating point type.
   /// @param uplo Specifies whether the upper or lower triangular part of the Hermitian matrix C
   /// is to be referenced:
   ///
   ///        if uplo = 'U' or 'u' then C is upper triangular,
   ///        if uplo = 'L' or 'l' then C is lower triangular.
   /// @param trans Specifies the operation to be perfomed as follows:
   ///
   ///        if trans = 'N' or 'n' then C := alpha*A*A'+beta*C
   ///        if trans = 'C' or 'c' then C := alpha*A'*A+beta*C
   /// @param n Specifies the order of the complex Hermitian matrix C.  n>=0
   /// @param k Specifies the other dimension of the complex matrix A (see below). k>=0
   /// @param alpha Real scalar.
   /// @param A Pointer to complex matrix.
   ///
   ///        if trans = 'N' or 'n' then A is n-by-k
   ///        if trans = 'C' or 'c' then A is k-by-n
   /// @param ldA Column length of the matrix A. If trans = 'N' or 'n' ldA>=n, otherwise ldA>=k.
   /// @param beta Real scalar.
   /// @param C Pointer to complex Hermitian n-by-n matrix C.
   /// Only the upper or lower triangular part of C is referenced, depending on the value of uplo above.
   /// @param ldC Column length of the matrix C.  ldC>=n
   /// @ingroup BLAS

   template <typename real_t>
   int HERK(char uplo, char trans, int_t n, int_t k, real_t alpha, complex<real_t> *A, int_t ldA, real_t beta, complex<real_t> *C, int_t ldC)
   {
      using std::toupper;
      uplo=toupper(uplo);
      trans=toupper(trans);
      using std::real;
      using std::conj;
      const complex<real_t> zero(0.0,0.0);
      const real_t rzero=0.0;
      const real_t one=1.0;
      int_t i,j,l;
      complex<real_t> *a,*c,*at;
      complex<real_t> t;
      real_t s;

      if((uplo!='U')&&(uplo!='L'))
         return -1;
      else if((trans!='N')&&(trans!='C'))
         return -2;
      else if(n<0)
         return -3;
      else if(k<0)
         return -4;
      else if(ldA<((trans=='N')?n:k))
         return -7;
      else if(ldC<n)
         return -10;
      else if((n==0)||(((alpha==rzero)||(k==0))&&(beta==one)))
         return 0;

      if(alpha==zero)
      {
         if(uplo=='U')
         {
            c=C;
            for(j=0;j<n;j++)
            {
               for(i=0;i<j;i++)
                  c[i]*=beta;
               c[j]=beta*real(c[j]);
               c+=ldC;
            }

         }
         else
         {
            c=C;
            for(j=0;j<n;j++)
            {
               c[j]=beta*real(c[j]);
               for(i=j+1;i<n;i++)
                  c[i]*=beta;
               c+=ldC;
            }
         }
      }
      else if(trans=='N')
      {
         if(uplo=='U')
         {
            c=C;
            for(j=0;j<n;j++)
            {
               for(i=0;i<j;i++)
                  c[i]*=beta;
               c[j]=beta*real(c[j]);
               a=A;
               for(l=0;l<k;l++)
               {
                  t=alpha*conj(a[j]);
                  for(i=0;i<j;i++)
                     c[i]+=t*a[i];
                  c[j]=real(c[j])+real(t*a[i]);
                  a+=ldA;
               }
               c+=ldC;
            }
         }
         else
         {
            c=C;
            for(j=0;j<n;j++)
            {
               c[j]=beta*real(c[j]);
               for(i=j+1;i<n;i++)
                  c[i]*=beta;
               a=A;
               for(l=0;l<k;l++)
               {
                  t=alpha*conj(a[j]);
                  c[j]=real(c[j])+real(t*a[j]);
                  for(i=j+1;i<n;i++)
                     c[i]+=t*a[i];
                  a+=ldA;
               }
               c+=ldC;
            }
         }
      }
      else
      {
         if(uplo=='U')
         {
            c=C;
            at=A;
            for(j=0;j<n;j++)
            {
               a=A;
               for(i=0;i<j;i++)
               {
                  t=zero;
                  for(l=0;l<k;l++)
                     t+=conj(a[l])*at[l];
                  c[i]=alpha*t+beta*c[i];
                  a+=ldA;
               }
               s=rzero;
               for(l=0;l<k;l++)
                  s+=real(conj(at[l])*at[l]);
               c[j]=alpha*s+beta*real(c[j]);
               at+=ldA;
               c+=ldC;
            }
         }
         else
         {
            at=A;
            c=C;
            for(j=0;j<n;j++)
            {
               s=rzero;
               for(l=0;l<k;l++)
                  s+=real(conj(at[l])*at[l]);
               c[j]=alpha*s+beta*real(c[j]);
               a=A+(j+1)*ldA;
               for(i=j+1;i<n;i++)
               {
                  t=zero;
                  for(l=0;l<k;l++)
                     t+=conj(a[l])*at[l];
                  c[i]=alpha*t+beta*c[i];
                  a+=ldA;
               }
               at+=ldA;
               c+=ldC;
            }
         }
      }
      return 0;
   }
#ifdef __latl_cblas
#include <cblas.h>

   template <> int HERK<float>(char uplo, char trans, int_t n, int_t k, float alpha, complex<float> *A, int_t ldA, float beta, complex<float> *C, int_t ldC)
   {
      using std::toupper;
      uplo=toupper(uplo);
      trans=toupper(trans);
      if((uplo!='U')&&(uplo!='L'))
         return -1;
      else if((trans!='N')&&(trans!='C'))
         return -2;
      else if(n<0)
         return -3;
      else if(k<0)
         return -4;
      else if(ldA<((trans=='N')?n:k))
         return -7;
      else if(ldC<n)
         return -10;

      const CBLAS_UPLO Uplo=(uplo=='U')?CblasUpper:CblasLower;
      const CBLAS_TRANSPOSE Trans=(trans=='N')?CblasNoTrans:((trans=='T')?CblasTrans:CblasConjTrans);

      cblas_cherk(CblasColMajor,Uplo,Trans,n,k,alpha,A,ldA,beta,C,ldC);

      return 0;
   }

   template <> int HERK<double>(char uplo, char trans, int_t n, int_t k, double alpha, complex<double> *A, int_t ldA, double beta, complex<double> *C, int_t ldC)
   {
      using std::toupper;
      uplo=toupper(uplo);
      trans=toupper(trans);
      if((uplo!='U')&&(uplo!='L'))
         return -1;
      else if((trans!='N')&&(trans!='C'))
         return -2;
      else if(n<0)
         return -3;
      else if(k<0)
         return -4;
      else if(ldA<((trans=='N')?n:k))
         return -7;
      else if(ldC<n)
         return -10;

      const CBLAS_UPLO Uplo=(uplo=='U')?CblasUpper:CblasLower;
      const CBLAS_TRANSPOSE Trans=(trans=='N')?CblasNoTrans:((trans=='T')?CblasTrans:CblasConjTrans);

      cblas_zherk(CblasColMajor,Uplo,Trans,n,k,alpha,A,ldA,beta,C,ldC);

      return 0;
   }

#endif
}
#endif

max_stars_repo_name: langou/latl

max_issues_repo_name: langou/latl

max_forks_repo_name: langou/latl


text:
#ifndef Photometry_h
#define Photometry_h
/** This is free and unencumbered software released into the public domain.
The authors of ISIS do not claim copyright on the contents of this file.
For more details about the LICENSE terms and the AUTHORS, you will
find files of those names at the top level of this repository. **/

/* SPDX-License-Identifier: CC0-1.0 */

#include <string>
#include <vector>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_roots.h>

namespace Isis {
  class Pvl;
  class PhotoModel;
  class AtmosModel;
  class NormModel;
  /**
   * @author ????-??-?? Unknown
   *
   * @internal
   *  @history 2007-08-02 Steven Lambright - Fixed memory leak
   *  @history 2008-03-07 Janet Barrett - Added SetPhotomWl method to allow
   *                      the application to set the p_normWavelength variable for
   *                      use by MoonAlbedo normalization.
   *  @history 2008-06-18 Steven Koechle - Fixed Documentation Errors
   *  @history 2008-07-09 Steven Lambright - Fixed unit test
   *  @history 2011-08-19 Sharmila Prasad - Implemented brentminimizer using GSL
   *  @history 2011-09-15 Sharmila Prasad - Implemented brent's root solver using GSL
   */
  class Photometry {
    public:
      Photometry(Pvl &pvl);
      Photometry() {};
      virtual ~Photometry();

      //! Calculate the surface brightness
      void Compute(double pha, double inc, double ema, double dn,
                   double &albedo, double &mult, double &base);
      void Compute(double pha, double inc, double ema, double deminc,
                   double demema, double dn, double &albedo,
                   double &mult, double &base);

      //! Set the wavelength
      virtual void SetPhotomWl(double wl);

      //! Double precision version of bracketing algorithm ported from Python.
      //! Solution bracketing for 1-D minimization routine.
      static void minbracket(double &xa, double &xb, double &xc, double &fa,
          double &fb, double &fc, double Func(double par, void *params),
          void *params);

      //! Brent's method 1-D minimization routine using GSL's r8Brent minimization Algorithm
      static int brentminimizer(double x_lower, double x_upper, gsl_function *Func,
          double & x_minimum, double tolerance);

      //! GSL's the Brent-Dekker method (Brent's method) combines an interpolation strategy
      //! with the bisection algorithm to estimate the root of the quadratic function.
      static int brentsolver(double x_lo, double x_hi, gsl_function *Func, double tolerance, double &root);

      PhotoModel *GetPhotoModel() const {
        return p_phtPmodel;
      }

      AtmosModel *GetAtmosModel() const {
        return p_phtAmodel;
      }

      NormModel *GetNormModel() const {
        return p_phtNmodel;
      }

    protected:
      AtmosModel *p_phtAmodel;
      PhotoModel *p_phtPmodel;
      NormModel *p_phtNmodel;
  };
};

#endif