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

#include <config.h>
#include <stdlib.h>
#include <math.h>

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

/*
gsl_spmatrix_fprintf()
  Print sparse matrix to file in MatrixMarket format:

M  N  NNZ
I1 J1 A(I1,J1)
...

Note that indices start at 1 and not 0
*/

int
gsl_spmatrix_fprintf(FILE *stream, const gsl_spmatrix *m,
                     const char *format)
{
  int status;

  /* print header */
  status = fprintf(stream, "%%%%MatrixMarket matrix coordinate real general\n");
  if (status < 0)
    {
      GSL_ERROR("fprintf failed for header", GSL_EFAILED);
    }

  /* print rows,columns,nnz */
  status = fprintf(stream, "%u\t%u\t%u\n",
                   (unsigned int) m->size1,
                   (unsigned int) m->size2,
                   (unsigned int) m->nz);
  if (status < 0)
    {
      GSL_ERROR("fprintf failed for dimension header", GSL_EFAILED);
    }

  if (GSL_SPMATRIX_ISTRIPLET(m))
    {
      size_t n;

      for (n = 0; n < m->nz; ++n)
        {
          status = fprintf(stream, "%u\t%u\t", (unsigned int) m->i[n] + 1, (unsigned int) m->p[n] + 1);
          if (status < 0)
            {
              GSL_ERROR("fprintf failed", GSL_EFAILED);
            }

          status = fprintf(stream, format, m->data[n]);
          if (status < 0)
            {
              GSL_ERROR("fprintf failed", GSL_EFAILED);
            }

          status = putc('\n', stream);
          if (status == EOF)
            {
              GSL_ERROR("putc failed", GSL_EFAILED);
            }
        }
    }
  else if (GSL_SPMATRIX_ISCCS(m))
    {
      size_t j, p;

      for (j = 0; j < m->size2; ++j)
        {
          for (p = m->p[j]; p < m->p[j + 1]; ++p)
            {
              status = fprintf(stream, "%u\t%u\t", (unsigned int) m->i[p] + 1, (unsigned int) j + 1);
              if (status < 0)
                {
                  GSL_ERROR("fprintf failed", GSL_EFAILED);
                }

              status = fprintf(stream, format, m->data[p]);
              if (status < 0)
                {
                  GSL_ERROR("fprintf failed", GSL_EFAILED);
                }

              status = putc('\n', stream);
              if (status == EOF)
                {
                  GSL_ERROR("putc failed", GSL_EFAILED);
                }
            }
        }
    }
  else if (GSL_SPMATRIX_ISCRS(m))
    {
      size_t i, p;

      for (i = 0; i < m->size1; ++i)
        {
          for (p = m->p[i]; p < m->p[i + 1]; ++p)
            {
              status = fprintf(stream, "%u\t%u\t", (unsigned int) i + 1, (unsigned int) m->i[p] + 1);
              if (status < 0)
                {
                  GSL_ERROR("fprintf failed", GSL_EFAILED);
                }

              status = fprintf(stream, format, m->data[p]);
              if (status < 0)
                {
                  GSL_ERROR("fprintf failed", GSL_EFAILED);
                }

              status = putc('\n', stream);
              if (status == EOF)
                {
                  GSL_ERROR("putc failed", GSL_EFAILED);
                }
            }
        }
    }
  else
    {
      GSL_ERROR("unknown sparse matrix type", GSL_EINVAL);
    }

  return GSL_SUCCESS;
}

gsl_spmatrix *
gsl_spmatrix_fscanf(FILE *stream)
{
  gsl_spmatrix *m;
  unsigned int size1, size2, nz;
  char buf[1024];
  int found_header = 0;

  /* read file until we find rows,cols,nz header */
  while (fgets(buf, 1024, stream) != NULL)
    {
      int c;

      /* skip comments */
      if (*buf == '%')
        continue;

      c = sscanf(buf, "%u %u %u",
                 &size1, &size2, &nz);
      if (c == 3)
        {
          found_header = 1;
          break;
        }
    }

  if (!found_header)
    {
      GSL_ERROR_NULL ("fscanf failed reading header", GSL_EFAILED);
    }

  m = gsl_spmatrix_alloc_nzmax((size_t) size1, (size_t) size2, (size_t) nz, GSL_SPMATRIX_TRIPLET);
  if (!m)
    {
      GSL_ERROR_NULL ("error allocating m", GSL_ENOMEM);
    }

  {
    unsigned int i, j;
    double val;

    while (fgets(buf, 1024, stream) != NULL)
      {
        int c = sscanf(buf, "%u %u %lg", &i, &j, &val);
        if (c < 3 || (i == 0) || (j == 0))
          {
            GSL_ERROR_NULL ("error in input file format", GSL_EFAILED);
          }
        else if ((i > size1) || (j > size2))
          {
            GSL_ERROR_NULL ("element exceeds matrix dimensions", GSL_EBADLEN);
          }
        else
          {
            /* subtract 1 from (i,j) since indexing starts at 1 */
            gsl_spmatrix_set(m, i - 1, j - 1, val);
          }
      }
  }

  return m;
}

int
gsl_spmatrix_fwrite(FILE *stream, const gsl_spmatrix *m)
{
  size_t items;

  /* write header: size1, size2, nz */

  items = fwrite(&(m->size1), sizeof(size_t), 1, stream);
  if (items != 1)
    {
      GSL_ERROR("fwrite failed on size1", GSL_EFAILED);
    }

  items = fwrite(&(m->size2), sizeof(size_t), 1, stream);
  if (items != 1)
    {
      GSL_ERROR("fwrite failed on size2", GSL_EFAILED);
    }

  items = fwrite(&(m->nz), sizeof(size_t), 1, stream);
  if (items != 1)
    {
      GSL_ERROR("fwrite failed on nz", GSL_EFAILED);
    }

  /* write m->i and m->data which are size nz in all storage formats */

  items = fwrite(m->i, sizeof(size_t), m->nz, stream);
  if (items != m->nz)
    {
      GSL_ERROR("fwrite failed on row indices", GSL_EFAILED);
    }

  items = fwrite(m->data, sizeof(double), m->nz, stream);
  if (items != m->nz)
    {
      GSL_ERROR("fwrite failed on data", GSL_EFAILED);
    }

  if (GSL_SPMATRIX_ISTRIPLET(m))
    {
      items = fwrite(m->p, sizeof(size_t), m->nz, stream);
      if (items != m->nz)
        {
          GSL_ERROR("fwrite failed on column indices", GSL_EFAILED);
        }
    }
  else if (GSL_SPMATRIX_ISCCS(m))
    {
      items = fwrite(m->p, sizeof(size_t), m->size2 + 1, stream);
      if (items != m->size2 + 1)
        {
          GSL_ERROR("fwrite failed on column indices", GSL_EFAILED);
        }
    }
  else if (GSL_SPMATRIX_ISCRS(m))
    {
      items = fwrite(m->p, sizeof(size_t), m->size1 + 1, stream);
      if (items != m->size1 + 1)
        {
          GSL_ERROR("fwrite failed on column indices", GSL_EFAILED);
        }
    }

  return GSL_SUCCESS;
}

int
gsl_spmatrix_fread(FILE *stream, gsl_spmatrix *m)
{
  size_t size1, size2, nz;
  size_t items;

  /* read header: size1, size2, nz */

  items = fread(&size1, sizeof(size_t), 1, stream);
  if (items != 1)
    {
      GSL_ERROR("fread failed on size1", GSL_EFAILED);
    }

  items = fread(&size2, sizeof(size_t), 1, stream);
  if (items != 1)
    {
      GSL_ERROR("fread failed on size2", GSL_EFAILED);
    }

  items = fread(&nz, sizeof(size_t), 1, stream);
  if (items != 1)
    {
      GSL_ERROR("fread failed on nz", GSL_EFAILED);
    }

  if (m->size1 != size1)
    {
      GSL_ERROR("matrix has wrong size1", GSL_EBADLEN);
    }
  else if (m->size2 != size2)
    {
      GSL_ERROR("matrix has wrong size2", GSL_EBADLEN);
    }
  else if (nz > m->nzmax)
    {
      GSL_ERROR("matrix nzmax is too small", GSL_EBADLEN);
    }
  else
    {
      /* read m->i and m->data arrays, which are size nz for all formats */

      items = fread(m->i, sizeof(size_t), nz, stream);
      if (items != nz)
        {
          GSL_ERROR("fread failed on row indices", GSL_EFAILED);
        }

      items = fread(m->data, sizeof(double), nz, stream);
      if (items != nz)
        {
          GSL_ERROR("fread failed on data", GSL_EFAILED);
        }

      m->nz = nz;

      if (GSL_SPMATRIX_ISTRIPLET(m))
        {
          items = fread(m->p, sizeof(size_t), nz, stream);
          if (items != nz)
            {
              GSL_ERROR("fread failed on column indices", GSL_EFAILED);
            }

          /* build binary search tree for m */
          gsl_spmatrix_tree_rebuild(m);
        }
      else if (GSL_SPMATRIX_ISCCS(m))
        {
          items = fread(m->p, sizeof(size_t), size2 + 1, stream);
          if (items != size2 + 1)
            {
              GSL_ERROR("fread failed on row pointers", GSL_EFAILED);
            }
        }
      else if (GSL_SPMATRIX_ISCRS(m))
        {
          items = fread(m->p, sizeof(size_t), size1 + 1, stream);
          if (items != size1 + 1)
            {
              GSL_ERROR("fread failed on column pointers", GSL_EFAILED);
            }
        }
    }

  return GSL_SUCCESS;
}


text:
#pragma once

#include "BgfxCallback.h"
#include "FrameBuffer.h"
#include "ShaderCompiler.h"

#include <Babylon/JsRuntime.h>
#include <Babylon/JsRuntimeScheduler.h>

#include <GraphicsImpl.h>

#include <napi/napi.h>

#include <bgfx/bgfx.h>
#include <bgfx/platform.h>
#include <bimg/bimg.h>
#include <bx/allocator.h>

#include <gsl/gsl>

#include <assert.h>

#include <arcana/containers/weak_table.h>
#include <arcana/threading/cancellation.h>
#include <unordered_map>

namespace Babylon
{
    struct TextureData final
    {
        ~TextureData()
        {
            if (OwnsHandle && bgfx::isValid(Handle))
            {
                bgfx::destroy(Handle);
            }
        }

        bgfx::TextureHandle Handle{bgfx::kInvalidHandle};
        bool OwnsHandle{true};
        uint32_t Width{0};
        uint32_t Height{0};
        uint32_t Flags{0};
        uint8_t AnisotropicLevel{0};
    };

    struct UniformInfo final
    {
        uint8_t Stage{};
        bgfx::UniformHandle Handle{bgfx::kInvalidHandle};
        bool YFlip{false};
    };

    struct ProgramData final
    {
        ProgramData() = default;
        ProgramData(const ProgramData&) = delete;
        ProgramData(ProgramData&&) = delete;

        ~ProgramData()
        {
            if (bgfx::isValid(Handle))
            {
                bgfx::destroy(Handle);
            }
        }

        std::unordered_map<std::string, uint32_t> VertexAttributeLocations{};
        std::unordered_map<std::string, UniformInfo> VertexUniformInfos{};
        std::unordered_map<std::string, UniformInfo> FragmentUniformInfos{};

        bgfx::ProgramHandle Handle{bgfx::kInvalidHandle};

        struct UniformValue
        {
            std::vector<float> Data{};
            uint16_t ElementLength{};
            bool YFlip{false};
        };

        std::unordered_map<uint16_t, UniformValue> Uniforms{};

        void SetUniform(bgfx::UniformHandle handle, gsl::span<const float> data, bool YFlip, size_t elementLength = 1)
        {
            UniformValue& value = Uniforms[handle.idx];
            value.Data.assign(data.begin(), data.end());
            value.ElementLength = static_cast<uint16_t>(elementLength);
            value.YFlip = YFlip;
        }
    };

    class IndexBufferData;
    class VertexBufferData;

    struct VertexArray final
    {
        ~VertexArray()
        {
            for (auto& vertexBufferPair : VertexBuffers)
            {
                bgfx::destroy(vertexBufferPair.second.VertexLayoutHandle);
            }
        }

        struct IndexBuffer
        {
            const IndexBufferData* Data{};
        };

        IndexBuffer indexBuffer{};

        struct VertexBuffer
        {
            const VertexBufferData* Data{};
            uint32_t StartVertex{};
            bgfx::VertexLayoutHandle VertexLayoutHandle{};
        };

        std::unordered_map<uint32_t, VertexBuffer> VertexBuffers;
    };

    class NativeEngine final : public Napi::ObjectWrap<NativeEngine>
    {
        static constexpr auto JS_CLASS_NAME = "_NativeEngine";
        static constexpr auto JS_ENGINE_CONSTRUCTOR_NAME = "Engine";

    public:
        NativeEngine(const Napi::CallbackInfo& info);
        NativeEngine(const Napi::CallbackInfo& info, JsRuntime& runtime);
        ~NativeEngine();

        static void Initialize(Napi::Env env);

    private:
        void Dispose();

        void Dispose(const Napi::CallbackInfo& info);
        Napi::Value HomogeneousDepth(const Napi::CallbackInfo& info);
        void RequestAnimationFrame(const Napi::CallbackInfo& info);
        Napi::Value CreateVertexArray(const Napi::CallbackInfo& info);
        void DeleteVertexArray(const Napi::CallbackInfo& info);
        void BindVertexArray(const Napi::CallbackInfo& info);
        Napi::Value CreateIndexBuffer(const Napi::CallbackInfo& info);
        void DeleteIndexBuffer(const Napi::CallbackInfo& info);
        void RecordIndexBuffer(const Napi::CallbackInfo& info);
        void UpdateDynamicIndexBuffer(const Napi::CallbackInfo& info);
        Napi::Value CreateVertexBuffer(const Napi::CallbackInfo& info);
        void DeleteVertexBuffer(const Napi::CallbackInfo& info);
        void RecordVertexBuffer(const Napi::CallbackInfo& info);
        void UpdateDynamicVertexBuffer(const Napi::CallbackInfo& info);
        Napi::Value CreateProgram(const Napi::CallbackInfo& info);
        Napi::Value GetUniforms(const Napi::CallbackInfo& info);
        Napi::Value GetAttributes(const Napi::CallbackInfo& info);
        void SetProgram(const Napi::CallbackInfo& info);
        void SetState(const Napi::CallbackInfo& info);
        void SetZOffset(const Napi::CallbackInfo& info);
        Napi::Value GetZOffset(const Napi::CallbackInfo& info);
        void SetDepthTest(const Napi::CallbackInfo& info);
        Napi::Value GetDepthWrite(const Napi::CallbackInfo& info);
        void SetDepthWrite(const Napi::CallbackInfo& info);
        void SetColorWrite(const Napi::CallbackInfo& info);
        void SetBlendMode(const Napi::CallbackInfo& info);
        void SetMatrix(const Napi::CallbackInfo& info);
        void SetInt(const Napi::CallbackInfo& info);
        void SetIntArray(const Napi::CallbackInfo& info);
        void SetIntArray2(const Napi::CallbackInfo& info);
        void SetIntArray3(const Napi::CallbackInfo& info);
        void SetIntArray4(const Napi::CallbackInfo& info);
        void SetFloatArray(const Napi::CallbackInfo& info);
        void SetFloatArray2(const Napi::CallbackInfo& info);
        void SetFloatArray3(const Napi::CallbackInfo& info);
        void SetFloatArray4(const Napi::CallbackInfo& info);
        void SetMatrices(const Napi::CallbackInfo& info);
        void SetMatrix3x3(const Napi::CallbackInfo& info);
        void SetMatrix2x2(const Napi::CallbackInfo& info);
        void SetFloat(const Napi::CallbackInfo& info);
        void SetFloat2(const Napi::CallbackInfo& info);
        void SetFloat3(const Napi::CallbackInfo& info);
        void SetFloat4(const Napi::CallbackInfo& info);
        Napi::Value CreateTexture(const Napi::CallbackInfo& info);
        void LoadTexture(const Napi::CallbackInfo& info);
        void LoadRawTexture(const Napi::CallbackInfo& info);
        void LoadCubeTexture(const Napi::CallbackInfo& info);
        void LoadCubeTextureWithMips(const Napi::CallbackInfo& info);
        Napi::Value GetTextureWidth(const Napi::CallbackInfo& info);
        Napi::Value GetTextureHeight(const Napi::CallbackInfo& info);
        void SetTextureSampling(const Napi::CallbackInfo& info);
        void SetTextureWrapMode(const Napi::CallbackInfo& info);
        void SetTextureAnisotropicLevel(const Napi::CallbackInfo& info);
        void SetTexture(const Napi::CallbackInfo& info);
        void DeleteTexture(const Napi::CallbackInfo& info);
        Napi::Value CreateFrameBuffer(const Napi::CallbackInfo& info);
        void DeleteFrameBuffer(const Napi::CallbackInfo& info);
        void BindFrameBuffer(const Napi::CallbackInfo& info);
        void UnbindFrameBuffer(const Napi::CallbackInfo& info);
        void DrawIndexed(const Napi::CallbackInfo& info);
        void Draw(const Napi::CallbackInfo& info);
        void Clear(const Napi::CallbackInfo& info);
        Napi::Value GetRenderWidth(const Napi::CallbackInfo& info);
        Napi::Value GetRenderHeight(const Napi::CallbackInfo& info);
        void SetViewPort(const Napi::CallbackInfo& info);
        Napi::Value GetHardwareScalingLevel(const Napi::CallbackInfo& info);
        void SetHardwareScalingLevel(const Napi::CallbackInfo& info);
        Napi::Value CreateImageBitmap(const Napi::CallbackInfo& info);
        Napi::Value ResizeImageBitmap(const Napi::CallbackInfo& info);
        void GetFrameBufferData(const Napi::CallbackInfo& info);

        void Draw(bgfx::Encoder* encoder, int fillMode);

        GraphicsImpl::UpdateToken& GetUpdateToken();

        std::shared_ptr<arcana::cancellation_source> m_cancellationSource{};

        ShaderCompiler m_shaderCompiler{};

        ProgramData* m_currentProgram{nullptr};
        arcana::weak_table<std::unique_ptr<ProgramData>> m_programDataCollection{};

        JsRuntime& m_runtime;
        GraphicsImpl& m_graphicsImpl;

        JsRuntimeScheduler m_runtimeScheduler;

        std::optional<GraphicsImpl::UpdateToken> m_updateToken{};

        void ScheduleRequestAnimationFrameCallbacks();
        bool m_requestAnimationFrameCallbacksScheduled{};

        bx::DefaultAllocator m_allocator{};
        uint64_t m_engineState{BGFX_STATE_DEFAULT};

        template<int size, typename arrayType>
        void SetTypeArrayN(const Napi::CallbackInfo& info);

        template<int size>
        void SetFloatN(const Napi::CallbackInfo& info);

        template<int size>
        void SetMatrixN(const Napi::CallbackInfo& info);

        // Scratch vector used for data alignment.
        std::vector<float> m_scratch{};

        std::vector<Napi::FunctionReference> m_requestAnimationFrameCallbacks{};

        const VertexArray* m_boundVertexArray{};
        FrameBuffer* m_boundFrameBuffer{};
    };
}

ext: h


text:
#pragma once

#include <gsl/string_span>
#include <string>

#include "params.h"

namespace particles {

class ParserParams {
public:

	static PartSysParam* Parse(PartSysParamId id,
							   gsl::cstring_span<> value,
	                           float emitterLifespan,
	                           float particleLifespan,
	                           bool& success);

	static PartSysParam* Parse(gsl::cstring_span<> value,
	                           float defaultValue,
	                           float parentLifespan,
	                           bool& success);

	static PartSysParamKeyframes* ParseKeyframes(gsl::cstring_span<> value, float parentLifespan);

	static PartSysParamRandom* ParseRandom(gsl::cstring_span<> value);

	static PartSysParamSpecial* ParseSpecial(gsl::cstring_span<> value);

	static PartSysParamConstant* ParseConstant(gsl::cstring_span<> value, float defaultValue, bool& success);

};

}

max_stars_repo_name: edoipi/TemplePlus

max_issues_repo_name: edoipi/TemplePlus

max_forks_repo_name: edoipi/TemplePlus


text:
/*
 * Copyright 2021 The DAPHNE Consortium
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#pragma once

#include <runtime/local/context/DaphneContext.h>
#include <runtime/local/datastructures/DataObjectFactory.h>
#include <runtime/local/datastructures/DenseMatrix.h>

#include <cblas.h>

#include <cassert>
#include <cstddef>

// ****************************************************************************
// Struct for partial template specialization
// ****************************************************************************

template<class DTRes, class DTLhs, class DTRhs>
struct MatMul {
    static void apply(DTRes *& res, const DTLhs * lhs, const DTRhs * rhs, DCTX(ctx)) = delete;
};

// ****************************************************************************
// Convenience function
// ****************************************************************************

template<class DTRes, class DTLhs, class DTRhs>
void matMul(DTRes *& res, const DTLhs * lhs, const DTRhs * rhs, DCTX(ctx)) {
    MatMul<DTRes, DTLhs, DTRhs>::apply(res, lhs, rhs, ctx);
}

// ****************************************************************************
// (Partial) template specializations for different data/value types
// ****************************************************************************

// ----------------------------------------------------------------------------
// DenseMatrix <- DenseMatrix, DenseMatrix
// ----------------------------------------------------------------------------

template<>
struct MatMul<DenseMatrix<float>, DenseMatrix<float>, DenseMatrix<float>> {
    static void apply(DenseMatrix<float> *& res, const DenseMatrix<float> * lhs, const DenseMatrix<float> * rhs, DCTX(ctx)) {
        const auto nr1 = static_cast<int>(lhs->getNumRows());
        const auto nc1 = static_cast<int>(lhs->getNumCols());
        const auto nc2 = static_cast<int>(rhs->getNumCols());
        assert((nc1 == static_cast<int>(rhs->getNumRows())) && "#cols of lhs and #rows of rhs must be the same");

        if(res == nullptr)
            res = DataObjectFactory::create<DenseMatrix<float>>(nr1, nc2, false);

        if(nr1 == 1 && nc2 == 1) // Vector-Vector
            res->set(0, 0, cblas_sdot(nc1, lhs->getValues(), 1, rhs->getValues(),
                    static_cast<int>(rhs->getRowSkip())));
        else if(nc2 == 1)        // Matrix-Vector
            cblas_sgemv(CblasRowMajor, CblasNoTrans, nr1, nc1, 1, lhs->getValues(),
                    static_cast<int>(lhs->getRowSkip()), rhs->getValues(),
                    static_cast<int>(rhs->getRowSkip()), 0,res->getValues(),
                    static_cast<int>(res->getRowSkip()));
        else                     // Matrix-Matrix
            cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, nr1, nc2, nc1,
                    1, lhs->getValues(), static_cast<int>(lhs->getRowSkip()), rhs->getValues(),
                static_cast<int>(rhs->getRowSkip()), 0, res->getValues(), static_cast<int>(res->getRowSkip()));
    }
};

template<>
struct MatMul<DenseMatrix<double>, DenseMatrix<double>, DenseMatrix<double>> {
    static void apply(DenseMatrix<double> *& res, const DenseMatrix<double> * lhs, const DenseMatrix<double> * rhs, DCTX(ctx)) {
        const auto nr1 = static_cast<int>(lhs->getNumRows());
        const auto nc1 = static_cast<int>(lhs->getNumCols());
        const auto nc2 = static_cast<int>(rhs->getNumCols());
        assert((nc1 == static_cast<int>(rhs->getNumRows())) && "#cols of lhs and #rows of rhs must be the same");

        if(res == nullptr)
            res = DataObjectFactory::create<DenseMatrix<double>>(nr1, nc2, false);

        if(nr1 == 1 && nc2 == 1) // Vector-Vector
            res->set(0, 0, cblas_ddot(nc1, lhs->getValues(), 1, rhs->getValues(),
                static_cast<int>(rhs->getRowSkip())));
        else if(nc2 == 1)        // Matrix-Vector
            cblas_dgemv(CblasRowMajor, CblasNoTrans, nr1, nc1, 1, lhs->getValues(),
            static_cast<int>(lhs->getRowSkip()), rhs->getValues(),static_cast<int>(rhs->getRowSkip()), 0,
                res->getValues(), static_cast<int>(res->getRowSkip()));
        else                     // Matrix-Matrix
            cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, nr1, nc2, nc1,
                    1, lhs->getValues(), static_cast<int>(lhs->getRowSkip()), rhs->getValues(),
                    static_cast<int>(rhs->getRowSkip()), 0, res->getValues(), static_cast<int>(res->getRowSkip()));
    }
};

max_stars_repo_name: daphne-eu/daphne

max_issues_repo_name: daphne-eu/daphne

max_forks_repo_name: daphne-eu/daphne


text:
/*
Project 2: A Simple Ocean Model
This code models a set of equations for a 2-layer model of Ekman heat transport defined by F. Codron 2012: Ekman Heat Transport for Slab Oceans
The model is designed to take input from the Met Office Unified Model (UM) in .dat file form. As a result, this model uses the same spatial
horizontal resolution as the UM (90 latitude grid points and 144 longitude points, using spacing of 2.0 and 2.5 degrees respectively). The code
solves the heat transport using a finite difference method, and a sparse linear algebra GSL solver.
Model can be ran on an Ubuntu terminal using the following set of commands:
	gcc -Wall -I/GSL_PATH/gsl/include -c pr2_ocean_transport_v7.c
	gcc -L/GSL_PATH/gsl/lib pr2_ocean_transport_v7.o -lgsl -lgslcblas -lm
	./a.out [-v {NO_TRANSPORT, DIFUSION_ONLY, DIFUSION_EKMAN}]
where the options for -v are integers defined in enum Transport.
Created by Jake K. Eager
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <limits.h>
#include <unistd.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_spmatrix.h>
#include <gsl/gsl_splinalg.h>

#define C_D 0.0013 // drag coefficient, dimensionless
#define C_V 4000.0 // specific heat capacity of water J/kg/K
#define EMISSIVITY 0.985 // effective blackbody emissivity of ocean water
#define HEAT_TRANSFER_COEFFICIENT 3000. /* heat transfer coefficient for water W/m2/K 50-3000*/
#define RHO_AIR 1.22 // air density kg/m3
#define RHO_WATER 1027 // km/m3
#define R_PLANET (1.12384*6.371e6) // radius of planet in m (as fraction of Earth radius)
#define EPSILON (0.00001) // 0.00001 pretty good need to find appropriate value
#define OMEGA (6.501e-6) // angular frequency planet, about its central axis
// #define OMEGA (7.272e-5) // Earth value for omega
#define SIGMA (5.67e-8) // stefan boltzman constant
#define H_S 50.0 // thickness of surface layer in m
#define H_D 150.0 // thickness of deep layer below
#define D (25000.0) // horizontal diffusion coefficient m2/s
// #define D (25000.0*50.0/(H_S+H_D)) // horizontal diffusion coefficient m2/s MIGHT NEED THIS FIX, UNSURE...
#define N_LATS 90 /* number of latitude points */
#define N_LONS 144 /* number of longitude points */
#define LAT_MIN -89.0
#define DELTA_LAT 2.0
#define LON_MIN 1.25
#define DELTA_LON 2.5
#define HOURS_PER_DAY 24.
#define MINUTES_PER_HOUR 60.
#define SECONDS_PER_MINUTE 60.
#define DELTA_T (12.0*MINUTES_PER_HOUR*SECONDS_PER_MINUTE) // time step of 1 hours
#define TIME_STEPS (20000*2) // for time step of 12 hours, 20,000 time steps needed for 10,000 day run
#define TIME_OUTPUT_FREQ (25.0*HOURS_PER_DAY*MINUTES_PER_HOUR*SECONDS_PER_MINUTE) // print update frequency in seconds, but first term is number of days
#define DATA_OUTPUT_FREQ (100.0*HOURS_PER_DAY*MINUTES_PER_HOUR*SECONDS_PER_MINUTE) // output frequency in seconds, but first term is number of days
#define TIME_OUTPUT_TOL (1e-6) // days
#define T_OFFSET 0.0 // IF you want to restart run, specify start time here, so it doesn't save over stuff, and change input temp files to output ones
#define TOL (1e-6)
#define MAX_ITER 100
// input data files
// #define LAND_MASK_DATA "input_data/land_mask_day_cont_100_180.dat"
#define LAND_MASK_DATA "input_data/land_mask_no_land.dat"
#define MAX_FILE_LINE_SIZE 4000
#define SW_FLUX_NET_DOWN_DATA "input_data/ProCb/surface_net_downward_shortwave_flux.dat"
#define LW_FLUX_DOWN_DATA "input_data/ProCb/surface_downwelling_longwave_flux_in_air.dat"
// #define LW_FLUX_DOWN_DATA "input_data/ProCb/surface_net_downward_longwave_flux.dat"
#define LATENT_UP_FLUX_DATA "input_data/ProCb/surface_upward_latent_heat_flux.dat"
#define SENSIBLE_UP_FLUX_DATA "input_data/ProCb/surface_upward_sensible_heat_flux.dat"
// #define INITIAL_SURFACE_TEMP_DATA "output_data/ProCb/T_surf_1000_days.dat"
// #define INITIAL_DEEP_TEMP_DATA "output_data/ProCb/T_deep_1000_days.dat"
#define INITIAL_SURFACE_TEMP_DATA "input_data/ProCb/surface_temperature.dat"
#define INITIAL_DEEP_TEMP_DATA "input_data/ProCb/surface_temperature.dat"
#define U_WIND_DATA "input_data/ProCb/x_wind.dat"
#define V_WIND_DATA "input_data/ProCb/y_wind.dat"
#define X_STRESS_DATA "input_data/ProCb/surface_downward_eastward_stress.dat"
#define Y_STRESS_DATA "input_data/ProCb/surface_downward_northward_stress.dat"
#define LATS_FILE "input_data/lats.dat"
#define LONS_FILE "input_data/lons.dat"
// output data files
#define MAX_FNAME_CHAR 100
// #define PYTHON_EXE "python"
// #define PYTHON_SCRIPT 		"plotting_scripts/temperature_contour_plot.py"
#define OUTPUT_UPWARD_Q_FLUX "output_data/ProCb/upward_surface_Q_flux_"
#define OUTPUT_SURFACE_TEMP_DATA "output_data/ProCb/T_surf_"
#define OUTPUT_DEEP_TEMP_DATA "output_data/ProCb/T_deep_"
#define OUTPUT_X_MASS_FLUX_DATA "output_data/ProCb/eastward_flow_"
#define OUTPUT_Y_MASS_FLUX_DATA "output_data/ProCb/northward_flow_"
#define OUTPUT_VERITCAL_FLUX_DATA "output_data/ProCb/vertical_mass_flux_"
#define OUTPUT_DATA_EXT "_days.dat"

typedef enum coords { THETA, PHI, N_COORDS } Coords;
typedef enum depths { SURFACE, DEEP, N_DEPTHS } Depths;
typedef enum transport {NO_TRANSPORT, DIFUSION_ONLY, DIFUSION_EKMAN} Transport;
typedef enum errors { MALLOC_ERR, FREE_ERR, FOPEN_ERR, DEPTH_ERR, DERIVATIVE_ERR, READ_DATA_ERR, SOLVER_ERR, VERSION_ERR } Errors;

typedef struct grid_vals { // values of latitude and longitude
	double lat[N_LATS];
	double lon[N_LONS];
} Grid_vals;

static int select_version(int argc, char **argv);
static void time_stepper(int n_times, int n_lats, int n_lons, int version);

int main(int argc, char **argv)
{
	int version = select_version(argc, argv);
	time_stepper(TIME_STEPS,N_LATS,N_LONS,version);

	return 0;
}

/* Selects the whether to solve using no horizontal transport, diffusion only or diffusion and Ekman transport. */
static int select_version(int argc, char **argv)
{
	extern char *optarg;
	int c, err = 0;
	// char *case_number=NULL;
	int version=DIFUSION_EKMAN; // default case is the step function test case
	static char usage[] = "usage: %s [-v version_number] [name ...]\n";

	while ((c = getopt(argc, argv, "v:")) != -1)
		switch (c) {
		case 'v':
			version=atoi(optarg);
			if(version!=NO_TRANSPORT && version!=DIFUSION_EKMAN && version!=DIFUSION_ONLY){
				fprintf(stderr, "Version selected not valid, please enter [-v {%i, %i, %i}]\n", NO_TRANSPORT, DIFUSION_ONLY, DIFUSION_EKMAN);
				exit(VERSION_ERR);
			}
			break;
		case '?':
			err = 1;
			break;
		}
	if (err) {
		fprintf(stderr, usage, argv[0]);
		exit(1);
	}
	return version;
}

/* Checks memory allocation has been successful */
static void *xmalloc(size_t n)
{
    void *p = malloc(n);
    if(p == NULL)   {
        fprintf(stderr,"Out of memory.\n");
        exit(MALLOC_ERR);
    }
    return p;
}

/* checks pointer being freed is not NULL */
static void xfree(void *p)
{
	if(p == NULL)   {
        fprintf(stderr,"Pointer is NULL, so cannot free.\n");
        exit(FREE_ERR);
    }
	free(p);
}

/* allocates a 2d array of pointers */
static int **create_2d_int_pointer(int n, int m)
{
	int **p = (int **)xmalloc(n*sizeof(int*));
	for(int i=0;i<n;i++){
		p[i]=(int *)xmalloc(m*sizeof(int));
	}
	return p;
}

/* frees a 2d array of pointers */
static void destroy_2d_int_pointer(int **p, int n)
{
	for(int i=0;i<n;i++){
		xfree(p[i]);
	}
	xfree(p);
}

/* allocates a 2d array of pointers */
static double **create_2d_pointer(int n, int m)
{
	double **p = (double **)xmalloc(n*sizeof(double*));
	for(int i=0;i<n;i++){
		p[i]=(double *)xmalloc(m*sizeof(double));
	}
	return p;
}

/* frees a 2d array of pointers */
static void destroy_2d_pointer(double **p, int n)
{
	for(int i=0;i<n;i++){
		xfree(p[i]);
	}
	xfree(p);
}

/* allocates a 3d array of pointers */
static double ***create_3d_pointer(int n, int m, int l)
{
	double ***p = (double ***)xmalloc(n*sizeof(double**));
	for(int i=0;i<n;i++){
		p[i]=(double **)xmalloc(m*sizeof(double*));
		for(int j=0;j<m;j++){
			p[i][j]=(double *)xmalloc(l*sizeof(double));
		}
	}
	return p;
}

/* frees a 3d array of pointers */
static void destroy_3d_pointer(double ***p, int n, int m)
{
	for(int i=0;i<n;i++){
		for(int j=0;j<m;j++){
			xfree(p[i][j]);
		}
		xfree(p[i]);
	}
	xfree(p);
}

/* checks file opened correctly */
static void *xfopen(char *filename, char *command)
{
	FILE *outfile;
	outfile = fopen(filename,command);
	if(outfile==NULL){
		fprintf(stderr, "%s could not be opened.\n", filename);
		exit(FOPEN_ERR);
	}
	return outfile;
}

static void construct_fname(char* fname, size_t fname_length, const char* base, const char* time, const char* ext)
{ // construct filename from base, time of run and appropriate extention
    strcpy(fname, base);
    strcat(fname,time);
    strcat(fname,ext);
}

/* reads in a .dat files called filename, that contains data with n_lats rows and n_lons collumns, ie UM output data */
static void read_input_data(double **var, char *filename, int n_lats, int n_lons)
{
	char line[MAX_FILE_LINE_SIZE];
	char *data = line;
	FILE *input;
	input = xfopen(filename, "r" );
	int offset;
	for(int i=0;i<n_lats;i++){
		fgets(line, MAX_FILE_LINE_SIZE, input);
		for(int j=0;j<n_lons;j++){
			int found = sscanf(data, " %le%n", &var[i][j], &offset);
			if(found!=1){
				fprintf(stderr, "%s does not contain enough longitudinal points.\n", filename);
				exit(READ_DATA_ERR);
			}
			data+=offset;
		}
		data=line;
	}
  	fclose(input);
}

/* reads in a .dat files called filename, that contains integer data with n_lats rows and n_lons collumns, ie UM output data */
static void read_input_int_data(int **var, char *filename, int n_lats, int n_lons)
{
	char line[MAX_FILE_LINE_SIZE];
	char *data = line;
	FILE *input;
	input = xfopen(filename, "r" );
	int offset;
	for(int i=0;i<n_lats;i++){
		fgets(line, MAX_FILE_LINE_SIZE, input);
		for(int j=0;j<n_lons;j++){
			int found = sscanf(data, " %d%n", &var[i][j], &offset);
			if(found!=1){
				fprintf(stderr, "%s does not contain enough longitudinal points.\n", filename);
				exit(READ_DATA_ERR);
			}
			data+=offset;
		}
		data=line;
	}
  	fclose(input);
}

/* writes out put data to filename on a grid of n_lats by n_lons. */
static void write_data(double **var, char *fbasename, char *time_str, Grid_vals *grid, int n_lats, int n_lons)
{
	char *fname = (char*) xmalloc((MAX_FNAME_CHAR)*sizeof(char));
	construct_fname(fname, MAX_FNAME_CHAR, fbasename, time_str, OUTPUT_DATA_EXT);
	FILE *output;
	output = xfopen(fname, "w" );
	for(int i=0;i<n_lats;i++){
		for(int j=0;j<n_lons;j++){
			fprintf(output, "%lg\t", var[i][j]);
		}
		fprintf(output, "\n");
	}
	fclose(output);
	xfree(fname);
}

static void calc_Q_flux(double **upward_Q_flux, double **T_surf, Grid_vals *grid, int n_lats, int n_lons)
{
	double **F_net_sw_down = create_2d_pointer(n_lats,n_lons);
	double **F_lw_down = create_2d_pointer(n_lats,n_lons);
	double **F_latent_up = create_2d_pointer(n_lats,n_lons);
	double **F_sensible_up = create_2d_pointer(n_lats,n_lons);
	read_input_data(F_net_sw_down,SW_FLUX_NET_DOWN_DATA,n_lats,n_lons);
	read_input_data(F_lw_down,LW_FLUX_DOWN_DATA,n_lats,n_lons);
	read_input_data(F_latent_up,LATENT_UP_FLUX_DATA,n_lats,n_lons);
	read_input_data(F_sensible_up,SENSIBLE_UP_FLUX_DATA,n_lats,n_lons);
	for(int i=0;i<n_lats;i++){
		for(int j=0;j<n_lons;j++){
			upward_Q_flux[i][j] = F_latent_up[i][j] + F_sensible_up[i][j] + EMISSIVITY*SIGMA*pow(T_surf[i][j],4) - (F_net_sw_down[i][j] + F_lw_down[i][j]);
		}
	}
	destroy_2d_pointer(F_net_sw_down,n_lats);
	destroy_2d_pointer(F_lw_down,n_lats);
	destroy_2d_pointer(F_latent_up,n_lats);
	destroy_2d_pointer(F_sensible_up,n_lats);
}

/* convert angle in degrees to radians for evaluating trigonometric functions in calculus */
static double deg_to_rad(double theta)
{
	return theta*M_PI/180.0;
}

/* calculates the first spatial derivative of A with respect to a (which may be theta or phi), on a spherical surface. Used mainly to
calculate derivative of mass transport components M_theta and M_phi. Boundary conditions are used for a spherical geometry, using A[i][0]
to be adjacent to A[i][n_lons-1], and also A[0][j] is adjacent to the opposite point in longitude (j+-n_lats/2) at the same latitude. */
static double dA_da(double **A, int i, int j, int a, int n_lats, int n_lons)
{
	double d_theta=deg_to_rad(DELTA_LAT);
	double d_phi=deg_to_rad(DELTA_LON);
	if(a==THETA){
		if(i==n_lats-1){
			int j_temp = j+n_lons/2;
			if(j_temp>=n_lons){
				j_temp -=n_lons;
			}
			return (A[i][j_temp]-A[i-1][j])/(2*d_theta);
		}
		else if(i==0){
			int j_temp = j+n_lons/2;
			if(j_temp>=n_lons){
				j_temp -=n_lons;
			}
			return (A[i+1][j]-A[i][j_temp])/(2*d_theta);
		}
		else{
			return (A[i+1][j]-A[i-1][j])/(2*d_theta);
		}
	}
	else if(a==PHI){
		if(j==n_lons-1){
			return (A[i][0]-A[i][j-1])/(2*d_phi);
		}
		else if(j==0){
			return (A[i][j+1]-A[i][n_lons-1])/(2*d_phi);
		}
		else{
			return (A[i][j+1]-A[i][j-1])/(2*d_phi);
		}
	}
	else{
		fprintf(stderr, "Specified derivative is not respect to either theta or phi\n");
		exit(DERIVATIVE_ERR);
	}
}

/* Returns the horizontal divergence of quanitity M. determines whether there is an up or downward flow between the ocean layers */
static double calculate_div_M(double ***M, Grid_vals *grid, int i, int j, int n_lats, int n_lons)
{
	double theta=deg_to_rad(grid->lat[i]);
	double dM_theta_dtheta = dA_da(M[THETA],i,j,THETA,n_lats,n_lons);
	double dM_phi_dphi = dA_da(M[PHI],i,j,PHI,n_lats,n_lons);
	double div_M=1./R_PLANET*(dM_theta_dtheta-tan(theta)*M[THETA][i][j]+1./cos(theta)*dM_phi_dphi);
	return div_M;
}

static void calc_flow_Sv(double ***Sv_flow, double ***M, Grid_vals *grid, int n_lats, int n_lons)
{
	for(int i=0;i<n_lats;i++){
		double theta=deg_to_rad(grid->lat[i]);
		for(int j=0;j<n_lons;j++){
			Sv_flow[PHI][i][j] = M[PHI][i][j] * R_PLANET * DELTA_LAT * M_PI / 180. / RHO_WATER / 1.0e6;
			Sv_flow[THETA][i][j] = M[THETA][i][j] * R_PLANET * cos(theta) * deg_to_rad(DELTA_LON) / RHO_WATER / 1.0e6;
		}
	}
}

static void process_output_data(double ***T, double ***M, Grid_vals *grid, int n_lats, int n_lons, double time)
{ // processes output data
    double days = T_OFFSET+time/(HOURS_PER_DAY*MINUTES_PER_HOUR*SECONDS_PER_MINUTE);
    days=round(days);
    printf("Data outputted at %lg\n", days);
    int day_int = days+TIME_OUTPUT_TOL;
    char time_str[10];
    sprintf(time_str, "%d", day_int);
    // char *fname = (char*) xmalloc((MAX_FNAME_CHAR)*sizeof(char));
    double **upward_Q_flux= create_2d_pointer(n_lats,n_lons);
    calc_Q_flux(upward_Q_flux,T[SURFACE],grid,n_lats,n_lons);
    // construct_fname(fname, MAX_FNAME_CHAR, OUTPUT_UPWARD_LW_SURFACE_FLUX, time_str, OUTPUT_DATA_EXT);
    write_data(upward_Q_flux, OUTPUT_UPWARD_Q_FLUX, time_str, grid, n_lats, n_lons);
    destroy_2d_pointer(upward_Q_flux,n_lats);
    double ***Sv_flow = create_3d_pointer(N_COORDS,N_LATS,N_LONS);
    calc_flow_Sv(Sv_flow,M,grid,n_lats,n_lons);
    write_data(Sv_flow[PHI], OUTPUT_X_MASS_FLUX_DATA, time_str, grid, n_lats, n_lons);
    write_data(Sv_flow[THETA], OUTPUT_Y_MASS_FLUX_DATA, time_str, grid, n_lats, n_lons);
 	destroy_3d_pointer(Sv_flow,N_COORDS,N_LATS);
    double **div_M= create_2d_pointer(n_lats,n_lons);
    for(int i=0;i<n_lats;i++){
		for(int j=0;j<n_lons;j++){
			div_M[i][j] = calculate_div_M(M,grid,i,j,n_lats,n_lons);
		}
	}
	write_data(div_M, OUTPUT_VERITCAL_FLUX_DATA, time_str, grid, n_lats, n_lons);
	destroy_2d_pointer(div_M,n_lats);
    // construct_fname(fname, MAX_FNAME_CHAR, OUTPUT_SURFACE_TEMP_DATA, time_str, OUTPUT_DATA_EXT);
    write_data(T[SURFACE], OUTPUT_SURFACE_TEMP_DATA, time_str, grid, n_lats, n_lons);
    // construct_fname(fname, MAX_FNAME_CHAR, OUTPUT_DEEP_TEMP_DATA, time_str, OUTPUT_DATA_EXT);
    write_data(T[DEEP], OUTPUT_DEEP_TEMP_DATA, time_str, grid, n_lats, n_lons);
    // xfree(fname);
}

/* returns height of layer given index h, from enum heights. */
static double height_of_slab(int h)
{
	if(h==SURFACE){
		return H_S;
	}
	else if(h==DEEP){
		return H_D;
	}
	else{
		fprintf(stderr, "Index for depth is out of range.\n");
		exit(DEPTH_ERR);
	}
}

/* Constructs a grid with latitude and longitude points on the sphere, starting at LAT_MIN and LON_MIN, increasing in increments of DELTA_LAT
and DELTA_LON */
static void construct_grid(Grid_vals *grid, int n_lats, int n_lons, char *lats_file, char *lons_file)
{
	FILE *lats_out;
	lats_out = xfopen(lats_file, "w" );
	for(int i=0;i<n_lats;i++){
		grid->lat[i]=LAT_MIN+(i*DELTA_LAT);
		fprintf(lats_out, "%i\t%lg\n", i, grid->lat[i]);
	}
	fclose(lats_out);

	FILE *lons_out;
	lons_out = xfopen(lons_file, "w" );
	for(int j=0;j<n_lons;j++){
		grid->lon[j]=LON_MIN+(j*DELTA_LON);
		fprintf(lons_out, "%i\t%lg\n", j, grid->lon[j]);
	}
	fclose(lons_out);
}

/* Calculates the surface stress on the ocean surface as a result of the atmospheric winds near the surface.
Reads in velocity data from .dat files */
static void calculate_surface_stress(double ***tau, Grid_vals *grid, int n_lats, int n_lons)
{
	double ***v = create_3d_pointer(N_COORDS,N_LATS,N_LONS);
	read_input_data(v[THETA],V_WIND_DATA,n_lats,n_lons);
	read_input_data(v[PHI],U_WIND_DATA,n_lats,n_lons);
	for(int coord=THETA;coord<N_COORDS;coord++){
		for(int i=0;i<n_lats;i++){
			for(int j=0;j<n_lons;j++){
				tau[coord][i][j]=C_D*RHO_AIR*v[coord][i][j]*fabs(v[coord][i][j]); // retains sign of velocity
			}
		}
	}
	destroy_3d_pointer(v,N_COORDS,N_LATS);
}

/* Returns the coriolis parameter at an latitude theta */
static double calculate_coriolis_parameter(double theta)
{ // f=2*Omega*sin(theta)
	return 2*OMEGA*sin(deg_to_rad(theta));
}

/* calculates the mass flux due to Ekman transport */
static void calculate_mass_flux(double ***M, int **land_mask, Grid_vals *grid, int n_lats, int n_lons)
{
	double ***tau = create_3d_pointer(N_COORDS,N_LATS,N_LONS);
	// read_input_data(tau[THETA],Y_STRESS_DATA,n_lats,n_lons);
	// read_input_data(tau[PHI],X_STRESS_DATA,n_lats,n_lons);
	calculate_surface_stress(tau, grid, n_lats, n_lons);
	for(int i=0;i<n_lats;i++){
		double f = calculate_coriolis_parameter(grid->lat[i]);
		for(int j=0;j<n_lons;j++){
			if(land_mask[i][j]==1){
				M[THETA][i][j]=0.;
				M[PHI][i][j]=0.;
			}
			else{
				M[THETA][i][j]=(EPSILON*tau[THETA][i][j]+f*tau[PHI][i][j])/(pow(EPSILON,2)+pow(f,2));
				M[PHI][i][j]=(EPSILON*tau[PHI][i][j]-f*tau[THETA][i][j])/(pow(EPSILON,2)+pow(f,2));
			}
		}
	}
	destroy_3d_pointer(tau,N_COORDS,N_LATS);
}

/* Plots the script called filename in python */
// static void plot_data(char *filename)
// {
// 	char command[PATH_MAX];
// 	snprintf(command, sizeof(command), "%s %s", PYTHON_EXE, filename);
// 	system( command );
// }

/* Calculates the index for the matrix used in the solver, before being compressed */
static int calculate_matrix_index(int height, int lat, int lon, int n_lats, int n_lons)
{
	return height*n_lats*n_lons+lat*n_lons+lon;
}

/* At maximum or minimum latitutde, the adjacent latitude point above or bellow, respectively, is located at the same latitude,
but at the new longitudinal point (from lon to new_lon) */
static int calculate_new_lon(int lon, int n_lons)
{
	int new_lon = lon+n_lons/2;
	if(new_lon>=n_lons){
		new_lon -= n_lons;
	}
	return new_lon;
}

/* calculates and sets the matrix elements for a 2 layer ocean with 144 longitude points and 90 latitude points.
version determines whether the matrix is calculated for a no horizontal transport model, diffusion only model or the full Ekman model. */
static void calculate_matrix(gsl_spmatrix *A, double ***M, int **land_mask, Grid_vals *grid, int n_lats, int n_lons, int version)
{
	double d_theta=deg_to_rad(DELTA_LAT);
	double d_phi=deg_to_rad(DELTA_LON);
	double s_dfsn = DELTA_T*D/pow(R_PLANET,2);

	for(int height=SURFACE; height<N_DEPTHS; height++){
		double thickness = height_of_slab(height);
		double s_ekmn = DELTA_T/(RHO_WATER*thickness*R_PLANET);
		for(int lat=0; lat<n_lats; lat++){
			double theta = deg_to_rad(grid->lat[lat]);
			for(int lon=0; lon<n_lons; lon++){
				double Aii=0.0;
				int i;
				i = calculate_matrix_index(height,lat,lon,n_lats,n_lons);
				/* central point - check if land */
				if(land_mask[lat][lon]==1){
					Aii=1.;
					gsl_spmatrix_set(A, i, i, Aii);
					continue;
				}
				double dM_theta_dtheta = dA_da(M[THETA],lat,lon,THETA,n_lats,n_lons);
				double dM_phi_dphi = dA_da(M[PHI],lat,lon,PHI,n_lats,n_lons);
				double M_theta = M[THETA][lat][lon];
				double M_phi = M[PHI][lat][lon];
				double div_M = calculate_div_M(M,grid,lat,lon,n_lats,n_lons);
				double Aij=0.0;
				int j;
				/* central point */
				Aii=1.;
				if(version != NO_TRANSPORT){
					// if(height==SURFACE){
					Aii+= 2.*s_dfsn*(1./pow(d_theta,2)+1./pow(cos(theta)*d_phi,2));
					// }
					if(version==DIFUSION_EKMAN){
						if(height==SURFACE){
							Aii+= s_ekmn*(dM_theta_dtheta-M_theta*tan(theta)+1./cos(theta)*dM_phi_dphi);
							if(div_M<0.0){
								Aii+=(-DELTA_T/(RHO_WATER*thickness)*div_M);
							}
						}
						else if(height==DEEP){
							Aii+= -s_ekmn*(dM_theta_dtheta-M_theta*tan(theta)+1./cos(theta)*dM_phi_dphi);
							if(div_M>0.0){
								Aii+=(+DELTA_T/(RHO_WATER*thickness)*div_M);
							}
						}
					}
				}
				gsl_spmatrix_set(A, i, i, Aii);
				/* Neigbouring points */
				int lon_j, lat_j;
				if(version != NO_TRANSPORT){
					/* neighbouring latitude - below */
					if(lat==0){
						lon_j = calculate_new_lon(lon,n_lons);
						lat_j = lat;
						j = calculate_matrix_index(height,lat_j,lon_j,n_lats,n_lons);
					}
					else{
						lon_j = lon;
						lat_j = lat-1;
						j = calculate_matrix_index(height,lat_j,lon_j,n_lats,n_lons);
					}
					Aij=-s_dfsn*(1./pow(d_theta,2)+tan(theta)/(2*d_theta));
					if(version==DIFUSION_EKMAN){
						if(height==SURFACE){
							Aij+= s_ekmn*(-M_theta/(2.*d_theta));
						}
						else{
							Aij+= -s_ekmn*(-M_theta/(2.*d_theta));
						}
					}
					if(land_mask[lat_j][lon_j]==1){
						Aii+=Aij;
						Aij=0.;
					}
					gsl_spmatrix_set(A, i, j, Aij);
					/* neighbouring latitude - above */
					if(lat==n_lats-1){
						lon_j = calculate_new_lon(lon,n_lons);
						lat_j = lat;
						j = calculate_matrix_index(height,lat_j,lon_j,n_lats,n_lons);
					}
					else{
						lon_j = lon;
						lat_j = lat+1;
						j = calculate_matrix_index(height,lat_j,lon_j,n_lats,n_lons);
					}
					Aij=-s_dfsn*(1./pow(d_theta,2)-tan(theta)/(2*d_theta));
					if(version==DIFUSION_EKMAN){
						if(height==SURFACE){
							Aij+= s_ekmn*(M_theta/(2.*d_theta));
						}
						else{
							Aij+= -s_ekmn*(M_theta/(2.*d_theta));
						}
					}
					if(land_mask[lat_j][lon_j]==1){
						Aii+=Aij;
						Aij=0.;
					}
					gsl_spmatrix_set(A, i, j, Aij);
					/* neighbouring longitude - left */
					if(lon==0){
						lon_j = n_lons-1;
						lat_j = lat;
						j = calculate_matrix_index(height,lat_j,lon_j,n_lats,n_lons);
					}
					else{
						lon_j = lon-1;
						lat_j = lat;
						j = calculate_matrix_index(height,lat_j,lon_j,n_lats,n_lons);
					}
					Aij=-s_dfsn*1./pow(cos(theta)*d_phi,2);
					if(version==DIFUSION_EKMAN){
						if(height==SURFACE){
							Aij+= s_ekmn*(-M_phi/(2.*d_phi*cos(theta)));
						}
						else{
							Aij+= -s_ekmn*(-M_phi/(2.*d_phi*cos(theta)));
						}
					}
					if(land_mask[lat_j][lon_j]==1){
						Aii+=Aij;
						Aij=0.;
					}
					gsl_spmatrix_set(A, i, j, Aij);
					/* neighbouring longitude - right */
					if(lon==n_lons-1){
						lon_j = 0;
						lat_j = lat;
						j = calculate_matrix_index(height,lat_j,lon_j,n_lats,n_lons);
					}
					else{
						lon_j = lon+1;
						lat_j = lat;
						j = calculate_matrix_index(height,lat_j,lon_j,n_lats,n_lons);
					}
					Aij=-s_dfsn*1./pow(cos(theta)*d_phi,2);
					if(version==DIFUSION_EKMAN){
						if(height==SURFACE){
							Aij+= s_ekmn*(M_phi/(2.*d_phi*cos(theta)));
						}
						else{
							Aij+= -s_ekmn*(M_phi/(2.*d_phi*cos(theta)));
						}
					}
					if(land_mask[lat_j][lon_j]==1){
						Aii+=Aij;
						Aij=0.;
					}
					gsl_spmatrix_set(A, i, j, Aij);
					gsl_spmatrix_set(A, i, i, Aii);
					/* interaction between deep and surface layer */
					if(version==DIFUSION_EKMAN){
						if(height==SURFACE && div_M>0.0){
							Aij=(-DELTA_T/(RHO_WATER*thickness)*div_M);
							j = calculate_matrix_index(DEEP,lat,lon,n_lats,n_lons);
							gsl_spmatrix_set(A, i, j, Aij);
						}
						else if(height==DEEP && div_M<0.0){
							Aij=(+DELTA_T/(RHO_WATER*thickness)*div_M);
							j = calculate_matrix_index(SURFACE,lat,lon,n_lats,n_lons);
							gsl_spmatrix_set(A, i, j, Aij);
						}
					}
				}
			}
		}
	}
}

/* Reads in atmospheric fluxes, net stellar sw flux (subrtracts reflected sw radiation from the incident)
and incoming lw flux from atmospheric emission */
static void calculate_F_a(double **F_a, int n_lats, int n_lons)
{
	double **F_net_sw_down = create_2d_pointer(n_lats,n_lons);
	double **F_lw_down = create_2d_pointer(n_lats,n_lons);
	double **F_latent_up = create_2d_pointer(n_lats,n_lons);
	double **F_sensible_up = create_2d_pointer(n_lats,n_lons);
	read_input_data(F_net_sw_down,SW_FLUX_NET_DOWN_DATA,n_lats,n_lons);
	read_input_data(F_lw_down,LW_FLUX_DOWN_DATA,n_lats,n_lons);
	read_input_data(F_latent_up,LATENT_UP_FLUX_DATA,n_lats,n_lons);
	read_input_data(F_sensible_up,SENSIBLE_UP_FLUX_DATA,n_lats,n_lons);
	for(int i=0;i<n_lats;i++){
		for(int j=0;j<n_lons;j++){
			F_a[i][j] = (F_net_sw_down[i][j] + F_lw_down[i][j]) - (F_latent_up[i][j] + F_sensible_up[i][j]);
		}
	}
	destroy_2d_pointer(F_net_sw_down,n_lats);
	destroy_2d_pointer(F_lw_down,n_lats);
	destroy_2d_pointer(F_latent_up,n_lats);
	destroy_2d_pointer(F_sensible_up,n_lats);
}

/* Checks if convetion condition is met, then calculates the associated flux. */
static void calculate_F_c(double **F_c, double ***T, int n_lats, int n_lons)
{
	for(int i=0;i<n_lats;i++){
		for(int j=0;j<n_lons;j++){
			if(T[DEEP][i][j]>T[SURFACE][i][j]){ /* convection condition */
				F_c[i][j] = HEAT_TRANSFER_COEFFICIENT*(T[DEEP][i][j]-T[SURFACE][i][j]);
			}
			else{
				F_c[i][j] = 0.0;
			}
		}
	}
}

/* calculated and sets the vector b in matrix Ax=b. EMISSIVITY*SIGMA*T^4 is a blackbody emission from one of the layers,
to ensure an energy exchange between the layers. 1.0/(RHO_WATER*C_V*H_S) prefactor converts a flux in W/m2 to K/s */
static void calculate_vector_b(double ***T, int **land_mask, gsl_vector *b, int n_lats, int n_lons)
{ // calculates b from matrix equation Ax=b
	double **F_a = create_2d_pointer(n_lats,n_lons);
	calculate_F_a(F_a,n_lats,n_lons);
	double **F_c = create_2d_pointer(n_lats,n_lons);
	calculate_F_c(F_c,T,n_lats,n_lons);
	double b_hij = 0.0;
	for(int h=0;h<N_DEPTHS;h++){
		for(int i=0;i<n_lats;i++){
			for(int j=0;j<n_lons;j++){
				if(land_mask[i][j]==1){
					b_hij = T[h][i][j];
				}
				else{
					double depth = height_of_slab(h);
					if(h==SURFACE){
						b_hij = DELTA_T/(RHO_WATER*C_V*depth)*(F_c[i][j]+F_a[i][j]-EMISSIVITY*SIGMA*pow(T[SURFACE][i][j],4))+ T[h][i][j];
						// b_hij = DELTA_T/(RHO_WATER*C_V*depth)*(F_c[i][j]+F_a[i][j]) + T[h][i][j];
					}
					else{
						b_hij = DELTA_T/(RHO_WATER*C_V*depth)*(-F_c[i][j]) + T[h][i][j];
					}
				}
				gsl_vector_set(b, (h*n_lats*n_lons+(i*n_lons+j)), b_hij);
			}
		}
	}
	destroy_2d_pointer(F_a,n_lats);
	destroy_2d_pointer(F_c,n_lats);
}

/* Uses GSL sparse linear algebra to solve Ax=b */
static void calculate_new_T(double ***T, double ***M, int **land_mask, Grid_vals *grid, int n_lats, int n_lons, int version)
{
	int n = N_DEPTHS*n_lats*n_lons;
	gsl_spmatrix *A = gsl_spmatrix_alloc(n ,n);
	gsl_spmatrix *C;
	gsl_vector *b = gsl_vector_alloc(n);
	gsl_vector *x = gsl_vector_alloc(n);

	calculate_matrix(A,M,land_mask,grid,n_lats,n_lons,version);
	calculate_vector_b(T,land_mask,b,n_lats,n_lons);
	C = gsl_spmatrix_ccs(A);

	/* now solve the system with the GMRES iterative solver */
    const double tol = TOL;  /* solution relative tolerance */
    const size_t max_iter = MAX_ITER; /* maximum iterations */
    const gsl_splinalg_itersolve_type *S = gsl_splinalg_itersolve_gmres;
    gsl_splinalg_itersolve *work = gsl_splinalg_itersolve_alloc(S, n, 0);
    size_t iter = 0;
    double residual;
    int status;

    /* initial guess x = 0 */
    gsl_vector_set_zero(x);
 //    for(int h=0;h<N_DEPTHS;h++){
	// 	for(int i=0;i<n_lats;i++){
	// 		for(int j=0;j<n_lons;j++){
	// 		gsl_vector_set(x, (h*n_lats*n_lons+(i*n_lons+j)), T[h][i][j]);
	// 		}
	// 	}
	// }

    /* solve the system A x = b */
    do{
		status = gsl_splinalg_itersolve_iterate(C, b, tol, x, work);

		if (status != GSL_SUCCESS && iter==max_iter-1){ /* Solver exits if solution diverges */
			residual = gsl_splinalg_itersolve_normr(work);
			fprintf(stderr, "iter %zu residual = %.12e\n", iter, residual);
		  	fprintf(stderr, "Diverged\n");
		  	exit(SOLVER_ERR);
		}
	}
	while (status == GSL_CONTINUE && ++iter < max_iter);

    /* output solution */
    for(int h=0;h<N_DEPTHS;h++){ /* sets the new temperature from the gsl vector */
	    for(int i=0;i<n_lats;i++){
			for(int j=0;j<n_lons;j++){
				T[h][i][j]=gsl_vector_get(x,h*n_lats*n_lons+(i*n_lons+j));
			}
		}
	}

    gsl_splinalg_itersolve_free(work);

	gsl_spmatrix_free(A);
	gsl_spmatrix_free(C);
	gsl_vector_free(b);
	gsl_vector_free(x);
}

/* Function that steps through selected time. First creates the horizontal grid, then read in temperatures,
and steps through n_times time steps. Writes the final Temperature data to files, and plots a selected plotting script. */
static void time_stepper(int n_times, int n_lats, int n_lons, int version)
{
	Grid_vals *grid = xmalloc(sizeof(Grid_vals));
	construct_grid(grid,n_lats,n_lons,LATS_FILE,LONS_FILE);
	double ***T=create_3d_pointer(N_DEPTHS,N_LATS,N_LONS);
	read_input_data(T[SURFACE],INITIAL_SURFACE_TEMP_DATA,n_lats,n_lons);
	read_input_data(T[DEEP],INITIAL_DEEP_TEMP_DATA,n_lats,n_lons);
	int **land_mask = create_2d_int_pointer(n_lats,n_lons);
	read_input_int_data(land_mask,LAND_MASK_DATA,n_lats,n_lons);
	double ***M=create_3d_pointer(N_COORDS,N_LATS,N_LONS);
	double time = 0.;
	double days = 0.;
	for(int t=0;t<n_times;t++){
		time = ((double)t+1.)*DELTA_T;
		calculate_mass_flux(M,land_mask,grid,n_lats,n_lons);
		calculate_new_T(T,M,land_mask,grid,n_lats,n_lons,version);
		if (fmod(time,TIME_OUTPUT_FREQ)<TIME_OUTPUT_TOL){
			days = T_OFFSET+time/(HOURS_PER_DAY*MINUTES_PER_HOUR*SECONDS_PER_MINUTE);
			printf("Days passed = %lg\n", days);
		}
		if (fmod(time,DATA_OUTPUT_FREQ)<TIME_OUTPUT_TOL){
			process_output_data(T, M, grid, n_lats, n_lons, time);
		}
	}
	destroy_2d_int_pointer(land_mask,n_lats);
	process_output_data(T, M, grid, n_lats, n_lons, time);
	// plot_data(PYTHON_SCRIPT);
	destroy_3d_pointer(M,N_COORDS,N_LATS);
	destroy_3d_pointer(T,N_DEPTHS,N_LATS);
	xfree(grid);
}

ext: c


text:
#ifndef cgsl_shim_h
#define cgsl_shim_h

#include <gsl/gsl_math.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_rng.h>

#endif

max_stars_repo_name: kongzii/SGSL

max_issues_repo_name: kongzii/SGSL

max_forks_repo_name: kongzii/SGSL


text:
#pragma once

#include <map>
#include <set>

#ifdef _NOGSL
class gsl_vector;
#else
#include <gsl/gsl_vector.h>
#endif

class Interface;

class SymbolicEvaluator {
	double otherError = 1;
public:
	virtual void setInputs(Interface* inputValues_p) = 0;

    virtual void run(const gsl_vector* ctrls_p, const set<int>& nodesSubset) = 0;
    virtual void run(const gsl_vector* ctrls_p) = 0;

    virtual double getErrorOnConstraint(int nodeid, gsl_vector* grad) = 0;
    virtual double getErrorOnConstraint(int nodeid) = 0;

    virtual double getErrorForAsserts(const set<int>& nodeid, gsl_vector* grad) = 0;
	virtual double getErrorForAssert(int assertId, gsl_vector* grad)=0;

    virtual double dist(int nid) = 0;
    virtual double dist(int nid, gsl_vector* grad) = 0;

	virtual void print() = 0;
	virtual void printFull() = 0;

	void resetOtherError() {
		otherError = 1;
	}

	void setOtherError(double error) {
		otherError = error;
	}

	double getOtherError() {
		return otherError;
	}

	bool isSetOtherError() {
		return otherError != 1;
	}
};

max_stars_repo_name: natebragg/sketch-backend

max_issues_repo_name: natebragg/sketch-backend

max_forks_repo_name: natebragg/sketch-backend


text:
#pragma once

#include <gsl/gsl_complex.h>

int compare(double n1, double n2, double precision = 1e-15);

int compare(const gsl_complex& c1, const gsl_complex& c2, double precision = 1e-15);
bool operator==(const gsl_complex& c1, const gsl_complex& c2);
bool operator!=(const gsl_complex& c1, const gsl_complex& c2);

max_stars_repo_name: antoniojkim/CalcPlusPlus

max_issues_repo_name: antoniojkim/CalcPlusPlus

max_forks_repo_name: antoniojkim/CalcPlusPlus


text:
#ifndef _PKMCMC_H_
#define _PKMCMC_H_

#include <vector>
#include <string>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_matrix.h>

class pkmcmc{
    int num_data, num_pars, workspace_size;
    std::vector<double> data, model, k;
    std::vector<double> theta_0, theta_i, param_vars, min, max;
    std::vector<std::vector<double>> cov, Psi;
    std::vector<bool> limit_pars;
    gsl_spline *Pk_bao, *Pk_nw;
    gsl_interp_accel *acc_bao, *acc_nw;
    double chisq_0, chisq_i, abs_err, rel_err, k_min, k_max;

    // The function to be integrated by GSL
    double model_func(std::vector<double> &pars, int j); //done

    // Performs the integral needed for the model value at each data point
    void model_calc(std::vector<double> &pars); //done

    // Randomly selects a new parameter realization for a trial
    void get_param_real(); //done

    // Compute the chi^2 value of the current model with the data
    double calc_chi_squared(); //done

    // Perform one Metropolis-Hastings trial
    bool trial(); //done

    // Write the current accepted parameter realization to the screen
    void write_theta_screen(); //done

    // Perform a number of trials to move from the initial parameter guess to a higher likelihood region
    void burn_in(int num_burn); //done

    // Tune the acceptance ratio to be about 0.234
    void tune_vars(); //done

    public:
        // Initialization function
        pkmcmc(std::string data_file, std::string cov_file, std::string pk_bao_file, std::string pk_nw_file,
               std::vector<double> &pars, std::vector<double> &vars, bool mock_avg); //done

        // Check that the vectors have the expected sizes
        void check_init(); //done

        // Set which parameters have priors and what those priors are.
        void set_param_limits(std::vector<bool> &lim_pars, std::vector<double> &min_in,
                              std::vector<double> &max_in); //done

        // Run a number of Metropolis-Hastings trials to form the MCMC chain
        void run_chain(int num_draws, std::string reals_file, bool new_chain);

        // Free the memory with various GSL class members
        void clean_up_gsl();

};

#endif


text:
/*  Copyright 2017 International Business Machines Corporation

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.  */

#pragma once

#include "kernelpp/types.h"
#include <gsl.h>

namespace mwe
{
    using kernelpp::maybe;
    using kernelpp::status;

    /* Operations  --------------------------------------------------------- */

    maybe<int> add(int a, int b);

    status add(const gsl::span<int> a,
               const gsl::span<int> b,
                     gsl::span<int> result
              );
}

max_stars_repo_name: rayglover-ibm/cuda-bindings

max_issues_repo_name: rayglover-ibm/cuda-bindings

max_forks_repo_name: rayglover-ibm/cuda-bindings


text:
/*
 * Take the input stream and tally the distribution of 8-bit bytes.
 * Calculate the entropy.
 */

#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_sf_log.h>
#include <gsl/gsl_math.h>

int main(int argc, char **argv) {
	int tally[256];
	int ch;
	int n = 0;
	double h = 0;

	for (int i=0; i<256; i++) tally[i] = 0;
	FILE *in;

	printf ("argc = %d\n", argc);

	if (argc < 2) {
		printf ("Input from stdin.\n");
		in = stdin;
	} else {
		printf ("Input file: %s\n", argv[1]);
		in = fopen(argv[1], "r");
		if (in == NULL) {
			printf("BAD file: Usage: %s <random_file>\n", argv[0]);
			return 1;
		}
	}

	while ((ch = fgetc(in)) != EOF) {
		tally[ch] ++;
		n ++;
	}
	fclose(in);

	printf("Frequency by byte:");
	for (int i=0; i<256; i++) {
		printf(" %d", tally[i]);
	}
	printf("\n");
	printf("Input length = %d\n", n);

	for (int i=0; i<256; i++) {
		double p = (double)(tally[i])/n;
		if (p>0) {
			h = h - p*gsl_sf_log(p);
		}
	}
	// Divide by ln(2) to produce base 2 logarithms to calculate bits of entropy.
	h = h/M_LN2;
	printf("Entropy (h) = %1.10f\n", h);

	// Maximum entropy on bytes is 8 bits. Efficiency is calculated by dividing by the max.
	double eta = h/8;
	printf("Efficiency (eta) = %1.10f\n", eta);
	printf("ln(1-eta) = %1.10f\n", gsl_sf_log(1-eta));

	// This prediction is based on the regression from test5.
	printf("Predicted value = %1.10f\n", -1.035*gsl_sf_log(n) + 3.562);

	return 0;
}

max_stars_repo_name: jonkanderson/BitBalancedTableHash

max_issues_repo_name: jonkanderson/BitBalancedTableHash

max_forks_repo_name: jonkanderson/BitBalancedTableHash


text:
#ifndef GLOTTAL_PHASES_H
#define GLOTTAL_PHASES_H

#include <array>
#include <vector>
#include <cstring>

#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_vector_complex.h>
#include <gsl/gsl_fft_halfcomplex.h>
#include <gsl/gsl_fft_complex.h>

#include "constants.h"

void estimatePhases(gsl_vector *dEGG, double T0, double *Oq, double *Cq);

#endif // GLOTTAL_PHASES_H

max_stars_repo_name: ichi-rika/glottal-inverse

max_issues_repo_name: ichi-rika/glottal-inverse

max_forks_repo_name: ichi-rika/glottal-inverse


text:
#pragma once

#include <variant>

#include <gsl/gsl-lite.hpp>

#include <clrie/instruction_factory.h>

namespace appmap { namespace cil {

struct op {
    ILOrdinalOpcode op;
};

struct token_op {
    ILOrdinalOpcode op;
    mdToken token;
};

struct long_op {
    ILOrdinalOpcode op;
    uint64_t value;
};

template <ILOrdinalOpcode Op>
struct token_opcode: token_op {
    token_opcode(mdToken token): token_op{Op, token} {}
};

namespace ops {
    struct ldarg { int index; };
    struct ldloc { int index; };
    struct stloc { int index; };

    using ldfld = token_opcode<Cee_Ldfld>;
    using stfld = token_opcode<Cee_Stfld>;
    using ldftn = token_opcode<Cee_Ldftn>;
    using call = token_opcode<Cee_Call>;
    using calli = token_opcode<Cee_Calli>;
    using callvirt = token_opcode<Cee_Callvirt>;
    using newobj = token_opcode<Cee_Newobj>;

    constexpr inline auto ldnull = op{Cee_Ldnull};
    constexpr inline auto pop = op{Cee_Pop};
    constexpr inline auto dup = op{Cee_Dup};

    struct ldc: long_op {
        ldc(uint64_t val): long_op{Cee_Ldc_I8, val} {}

        template <typename Ret, typename... Args>
        ldc(Ret(*fn)(Args...)): ldc(reinterpret_cast<uint64_t>(fn)) {}
    };
}

using instruction = std::variant<ops::ldarg, ops::ldloc, ops::stloc, op, token_op, long_op>;

clrie::instruction_factory::instruction_sequence compile(const gsl::span<const instruction> code, const clrie::instruction_factory &factory);
clrie::instruction_factory::instruction_sequence compile(std::initializer_list<const instruction> code, const clrie::instruction_factory &factory);

}}

max_stars_repo_name: applandinc/appmap-dotnet

max_issues_repo_name: applandinc/appmap-dotnet

max_forks_repo_name: applandinc/appmap-dotnet


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

/* Modified by Tuomo Keskitalo to add Nelder Mead Simplex test suite */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_ieee_utils.h>

#include "test_funcs.h"

unsigned int fcount, gcount;

int
test_fdf(const char * desc, gsl_multimin_function_fdf *f,
         initpt_function initpt, const gsl_multimin_fdfminimizer_type *T);

int
test_f(const char * desc, gsl_multimin_function *f, initpt_function initpt,
       const gsl_multimin_fminimizer_type *T);

int
main (void)
{
  gsl_ieee_env_setup ();

  {
    const gsl_multimin_fdfminimizer_type *fdfminimizers[6];
    const gsl_multimin_fdfminimizer_type ** T;

    fdfminimizers[0] = gsl_multimin_fdfminimizer_steepest_descent;
    fdfminimizers[1] = gsl_multimin_fdfminimizer_conjugate_pr;
    fdfminimizers[2] = gsl_multimin_fdfminimizer_conjugate_fr;
    fdfminimizers[3] = gsl_multimin_fdfminimizer_vector_bfgs;
    fdfminimizers[4] = gsl_multimin_fdfminimizer_vector_bfgs2;
    fdfminimizers[5] = 0;

    T = fdfminimizers;

    while (*T != 0)
      {
        test_fdf("Roth", &roth, roth_initpt,*T);
        test_fdf("Wood", &wood, wood_initpt,*T);
        test_fdf("Rosenbrock", &rosenbrock, rosenbrock_initpt,*T);
        T++;
      }

    T = fdfminimizers;

    while (*T != 0)
      {
        test_fdf("NRoth", &Nroth, roth_initpt,*T);
        test_fdf("NWood", &Nwood, wood_initpt,*T);
        test_fdf("NRosenbrock", &Nrosenbrock, rosenbrock_initpt,*T);
        T++;
      }
  }

  {
    const gsl_multimin_fminimizer_type *fminimizers[4];
    const gsl_multimin_fminimizer_type ** T;

    fminimizers[0] = gsl_multimin_fminimizer_nmsimplex;
    fminimizers[1] = gsl_multimin_fminimizer_nmsimplex2;
    fminimizers[2] = gsl_multimin_fminimizer_nmsimplex2rand;
    fminimizers[3] = 0;

    T = fminimizers;

    while (*T != 0)
      {
        test_f("Roth", &roth_fmin, roth_initpt,*T);
        test_f("Wood", &wood_fmin, wood_initpt,*T);
        test_f("Rosenbrock", &rosenbrock_fmin, rosenbrock_initpt,*T);
        test_f("Spring", &spring_fmin, spring_initpt,*T);
        T++;
      }
  }

  exit (gsl_test_summary());
}

int
test_fdf(const char * desc,
         gsl_multimin_function_fdf *f,
         initpt_function initpt,
         const gsl_multimin_fdfminimizer_type *T)
{
  int status;
  size_t iter = 0;
  double step_size;

  gsl_vector *x = gsl_vector_alloc (f->n);

  gsl_multimin_fdfminimizer *s;
  fcount = 0; gcount = 0;

  (*initpt) (x);

  step_size = 0.1 * gsl_blas_dnrm2 (x);

  s = gsl_multimin_fdfminimizer_alloc(T, f->n);

  gsl_multimin_fdfminimizer_set (s, f, x, step_size, 0.1);

#ifdef DEBUG
  printf("x "); gsl_vector_fprintf (stdout, s->x, "%g");
  printf("g "); gsl_vector_fprintf (stdout, s->gradient, "%g");
#endif

  do
    {
      iter++;
      status = gsl_multimin_fdfminimizer_iterate(s);

#ifdef DEBUG
      printf("%i: \n",iter);
      printf("x "); gsl_vector_fprintf (stdout, s->x, "%g");
      printf("g "); gsl_vector_fprintf (stdout, s->gradient, "%g");
      printf("f(x) %g\n",s->f);
      printf("dx %g\n",gsl_blas_dnrm2(s->dx));
      printf("\n");
#endif

      status = gsl_multimin_test_gradient(s->gradient,1e-3);
    }
  while (iter < 5000 && status == GSL_CONTINUE);

  status |= (fabs(s->f) > 1e-5);

  gsl_test(status, "%s, on %s: %i iters (fn+g=%d+%d), f(x)=%g",
           gsl_multimin_fdfminimizer_name(s),desc, iter, fcount, gcount, s->f);

  gsl_multimin_fdfminimizer_free(s);
  gsl_vector_free(x);

  return status;
}

int
test_f(const char * desc, gsl_multimin_function *f, initpt_function initpt,
       const gsl_multimin_fminimizer_type *T)
{
  int status;
  size_t i, iter = 0;

  gsl_vector *x = gsl_vector_alloc (f->n);

  gsl_vector *step_size = gsl_vector_alloc (f->n);

  gsl_multimin_fminimizer *s;

  fcount = 0; gcount = 0;
  (*initpt) (x);

  for (i = 0; i < f->n; i++)
    gsl_vector_set (step_size, i, 1);

  s = gsl_multimin_fminimizer_alloc(T, f->n);

  gsl_multimin_fminimizer_set (s, f, x, step_size);

#ifdef DEBUG
  printf("x "); gsl_vector_fprintf (stdout, s->x, "%g");
#endif

  do
    {
      iter++;
      status = gsl_multimin_fminimizer_iterate(s);

#ifdef DEBUG
      printf("%i: \n",iter);
      printf("x "); gsl_vector_fprintf (stdout, s->x, "%g");
      printf("f(x) %g\n", gsl_multimin_fminimizer_minimum (s));
      printf("size: %g\n", gsl_multimin_fminimizer_size (s));
      printf("\n");
#endif

      status = gsl_multimin_test_size (gsl_multimin_fminimizer_size (s),
                                       1e-3);
    }
  while (iter < 5000 && status == GSL_CONTINUE);

  status |= (fabs(s->fval) > 1e-5);

  gsl_test(status, "%s, on %s: %d iter (fn=%d), f(x)=%g",
           gsl_multimin_fminimizer_name(s),desc, iter, fcount, s->fval);

  gsl_multimin_fminimizer_free(s);
  gsl_vector_free(x);
  gsl_vector_free(step_size);

  return status;
}


text:
/**
 *
 * @file qwrapper_slansy.c
 *
 *  PLASMA core_blas quark wrapper
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Julien Langou
 * @author Henricus Bouwmeester
 * @author Mathieu Faverge
 * @date 2010-11-15
 * @generated s Tue Jan  7 11:44:57 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 **/
void QUARK_CORE_slansy(Quark *quark, Quark_Task_Flags *task_flags,
                       int norm, PLASMA_enum uplo, int N,
                       const float *A, int LDA, int szeA,
                       int szeW, float *result)
{
    szeW = max(1, szeW);
    DAG_CORE_LANSY;
    QUARK_Insert_Task(quark, CORE_slansy_quark, task_flags,
        sizeof(PLASMA_enum),                &norm,  VALUE,
        sizeof(PLASMA_enum),                &uplo,  VALUE,
        sizeof(int),                        &N,     VALUE,
        sizeof(float)*szeA,     A,             INPUT,
        sizeof(int),                        &LDA,   VALUE,
        sizeof(float)*szeW,                 NULL,          SCRATCH,
        sizeof(float),                      result,        OUTPUT,
        0);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_slansy_f1(Quark *quark, Quark_Task_Flags *task_flags,
                          PLASMA_enum norm, PLASMA_enum uplo, int N,
                          const float *A, int LDA, int szeA,
                          int szeW, float *result,
                          float *fake, int szeF)
{
    szeW = max(1, szeW);
    DAG_CORE_LANSY;

    if ( result == fake ) {
        QUARK_Insert_Task(quark, CORE_slansy_quark, task_flags,
            sizeof(PLASMA_enum),                &norm,  VALUE,
            sizeof(PLASMA_enum),                &uplo,  VALUE,
            sizeof(int),                        &N,     VALUE,
            sizeof(float)*szeA,     A,             INPUT,
            sizeof(int),                        &LDA,   VALUE,
            sizeof(float)*szeW,                 NULL,          SCRATCH,
            sizeof(float)*szeF,                 result,        OUTPUT | GATHERV,
            0);
    } else {
        QUARK_Insert_Task(quark, CORE_slansy_f1_quark, task_flags,
            sizeof(PLASMA_enum),                &norm,  VALUE,
            sizeof(PLASMA_enum),                &uplo,  VALUE,
            sizeof(int),                        &N,     VALUE,
            sizeof(float)*szeA,     A,             INPUT,
            sizeof(int),                        &LDA,   VALUE,
            sizeof(float)*szeW,                 NULL,          SCRATCH,
            sizeof(float),                      result,        OUTPUT,
            sizeof(float)*szeF,                 fake,          OUTPUT | GATHERV,
            0);
    }
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_slansy_quark = PCORE_slansy_quark
#define CORE_slansy_quark PCORE_slansy_quark
#endif
void CORE_slansy_quark(Quark *quark)
{
    float *normA;
    int norm;
    PLASMA_enum uplo;
    int N;
    float *A;
    int LDA;
    float *work;

    quark_unpack_args_7(quark, norm, uplo, N, A, LDA, work, normA);
    *normA = LAPACKE_slansy_work(
        LAPACK_COL_MAJOR,
        lapack_const(norm), lapack_const(uplo),
        N, A, LDA, work);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_slansy_f1_quark = PCORE_slansy_f1_quark
#define CORE_slansy_f1_quark PCORE_slansy_f1_quark
#endif
void CORE_slansy_f1_quark(Quark *quark)
{
    float *normA;
    int norm;
    PLASMA_enum uplo;
    int N;
    float *A;
    int LDA;
    float *work;
    float *fake;

    quark_unpack_args_8(quark, norm, uplo, N, A, LDA, work, normA, fake);
    *normA = LAPACKE_slansy_work(
        LAPACK_COL_MAJOR,
        lapack_const(norm),
        lapack_const(uplo),
        N, A, LDA, work);
}


text:
/*
 -----------------------------------------------------------------
 ** Top contributors:
 **   Shiqi Wang and Suman Jana
 ** This file is part of the ReluVal project.
 ** Copyright (c) 2018-2019 by the authors listed in the file LICENSE
 ** and their institutional affiliations.
 ** All rights reserved.
 -----------------------------------------------------------------
 */

#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>
#include <time.h>
#include <cblas.h>

#ifndef MATRIX_H
#define MATRIX_H

/* Define the structure of Matrix */
struct Matrix
{
    float* data;
    int row, col;
};

/* add the constant to matrix */
void add_constant(struct Matrix* A, float alpha);

/*matrix multiplication with factors */
void matmul_with_factor(struct Matrix* A,\
                        struct Matrix* B,\
                        struct Matrix* C,\
                        float alpha,\
                        float beta);

/* matrix multiplication */
void matmul(struct Matrix* A,\
            struct Matrix* B,\
            struct Matrix* C);

/* matrix multiplication with bias */
void matmul_with_bias(struct Matrix* A,\
                      struct Matrix* B,\
                      struct Matrix* C);

/* element-wise multiplication */
void multiply(struct Matrix* A, struct Matrix* B);

/* print matrix */
void printMatrix(struct Matrix* A);

/* print matrix to the file */
void fprintMatrix(FILE *fp, struct Matrix* A);

/* takes the relu of the matrix */
void relu(struct Matrix* A);

#endif

max_stars_repo_name: cherrywoods/HorizontalCAS

max_issues_repo_name: ahmed-irfan/HorizontalCAS

max_forks_repo_name: ahmed-irfan/HorizontalCAS


text:
#pragma once

#define NOMINMAX

#ifndef _WIN32
#include <wsl/winadapter.h>
#include "directml_guids.h"
#endif

#include <unordered_map>
#include <vector>
#include <iostream>
#include <filesystem>
#include <variant>
#include <fstream>
#include <deque>
#include <optional>
#include <set>
#include <string_view>
#include <wrl/client.h>
#include <wil/result.h>
#include <gsl/gsl>
#include <DirectML.h>
#include "DirectMLX.h"
#include <rapidjson/document.h>
#include <rapidjson/istreamwrapper.h>
#include <fmt/format.h>
#include <half.hpp>

template<class... Ts> struct overload : Ts... { using Ts::operator()...; };
template<class... Ts> overload(Ts...) -> overload<Ts...>;

max_stars_repo_name: miaobin/DirectML

max_issues_repo_name: miaobin/DirectML

max_forks_repo_name: miaobin/DirectML


text:
#ifndef CPP_MICROBENCH_BENCH_OUTPUT_H
#define CPP_MICROBENCH_BENCH_OUTPUT_H

#include <string>
#include <chrono>
#include <fstream>
#include <iomanip>

#include <fmt/os.h>
#include <fmt/ostream.h>
#include <fmt/chrono.h>

#include <gsl/util>

#include "util/aligned_vector.hpp"
#include "util/alloc_array.hpp"
#include "util/trivial_aligned_array.hpp"
#include "util/plain_array.hpp"

template <typename array_type>
auto generate_array(std::size_t max) {
  using namespace std::chrono;
  auto t1 = high_resolution_clock::now();
  array_type v(max);
  auto t2 = high_resolution_clock::now();
  for (std::size_t i = 0; i < max; ++i) {
    v[i] = 1.0F / gsl::narrow_cast<float>(i+1);
  }
  auto t3 = high_resolution_clock::now();
  fmt::print("Create: {}\n", duration_cast<microseconds>(t2-t1));
  fmt::print("Fill: {}\n", duration_cast<microseconds>(t3-t2));
  return v;
}

template <typename vector_type>
auto generate_vector(std::size_t max) {
  using namespace std::chrono;
  auto t1 = high_resolution_clock::now();
  vector_type v;
  v.reserve(max);
  auto t2 = high_resolution_clock::now();
  for (std::size_t i = 0; i < max; ++i) {
    v.push_back(1.0F / gsl::narrow_cast<float>(i+1));
  }
  auto t3 = high_resolution_clock::now();
  fmt::print("Create: {}\n", duration_cast<microseconds>(t2-t1));
  fmt::print("Fill: {}\n", duration_cast<microseconds>(t3-t2));
  return v;
}

void fmt_output_vector(const auto & v) {
  auto file_name = std::string(static_cast<const char *>(__func__)) + ".txt";
  auto file = fmt::output_file(file_name);
  file.print("{}\n", v.size());

  for (const auto & x : v) {
    file.print("{:.18f}\n", x);
  }
}

void stream_output_vector(const auto & v) {
  auto file_name = std::string(static_cast<const char *>(__func__)) + ".txt";
  std::ofstream file{file_name};
  file << v.size() << '\n';

  for (const auto & x : v) {
    file << std::setprecision(18) << x << '\n';
  }
}

template <typename array_type>
auto bench_fmt_output_array(std::size_t n) {
  using namespace std::chrono;

  auto t1 = high_resolution_clock::now();
  auto v = generate_array<array_type>(n);
  auto t2 = high_resolution_clock::now();
  fmt_output_vector(v);
  auto t3 = high_resolution_clock::now();

  return std::tuple {__func__, t2 - t1, t3 - t2};
}

template <typename vector_type>
auto bench_fmt_output_vector(std::size_t n) {
  using namespace std::chrono;

  auto t1 = high_resolution_clock::now();
  auto v = generate_vector<vector_type>(n);
  auto t2 = high_resolution_clock::now();
  fmt_output_vector(v);
  auto t3 = high_resolution_clock::now();

  return std::tuple {__func__, t2 - t1, t3 - t2};
}

template <typename vector_type>
auto bench_stream_vector(std::size_t n) {
  using namespace std::chrono;

  auto t1 = high_resolution_clock::now();
  auto v = generate_vector<vector_type>(n);
  auto t2 = high_resolution_clock::now();
  stream_output_vector(v);
  auto t3 = high_resolution_clock::now();

  return std::tuple {__func__, t2 - t1, t3 - t2};
}

void print_bench(auto t) {
  using namespace std::chrono;
  auto d1 = duration_cast<microseconds>(std::get<1>(t));
  auto d2 = duration_cast<microseconds>(std::get<2>(t));
  fmt::print("\nVersion: {}\n", std::get<0>(t));
  fmt::print("Generation {} us\n", d1.count());
  fmt::print("Writing {} us\n", d2.count());
}

#endif//CPP_MICROBENCH_BENCH_OUTPUT_H


text:
#ifndef PCA_H
#define PCA_H

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>

typedef struct
{
    float** matrix;
    int l;
    int c;
}
pca_matrix;

pca_matrix empirical_mean(pca_matrix X);

pca_matrix mean_deviations(pca_matrix X, pca_matrix U);

pca_matrix matrix_multiplication(pca_matrix A, pca_matrix B);

pca_matrix transposed_matrix(pca_matrix A);

pca_matrix covariance_matrix(pca_matrix B);

void eigenvectors_and_eigenvalues(pca_matrix C, gsl_vector *eval, gsl_matrix *evec);

void pca_matrix_free(pca_matrix M);

#endif /* PCA_H */

max_stars_repo_name: gh-duarteljd/faceBIT

max_issues_repo_name: gh-duarteljd/faceBIT

max_forks_repo_name: gh-duarteljd/faceBIT


text:
/**
 *
 * @file qwrapper_dtrtri.c
 *
 *  PLASMA core_blas quark wrapper
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Julien Langou
 * @author Henricus Bouwmeester
 * @author Mathieu Faverge
 * @date 2010-11-15
 * @generated d Tue Jan  7 11:44:56 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 **/
void QUARK_CORE_dtrtri(Quark *quark, Quark_Task_Flags *task_flags,
                       PLASMA_enum uplo, PLASMA_enum diag,
                       int n, int nb,
                       double *A, int lda,
                       PLASMA_sequence *sequence, PLASMA_request *request,
                       int iinfo)
{
    QUARK_Insert_Task(
        quark, CORE_dtrtri_quark, task_flags,
        sizeof(PLASMA_enum),                &uplo,      VALUE,
        sizeof(PLASMA_enum),                &diag,      VALUE,
        sizeof(int),                        &n,         VALUE,
        sizeof(double)*nb*nb,    A,                 INOUT,
        sizeof(int),                        &lda,       VALUE,
        sizeof(PLASMA_sequence*),           &sequence,  VALUE,
        sizeof(PLASMA_request*),            &request,   VALUE,
        sizeof(int),                        &iinfo,     VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_dtrtri_quark = PCORE_dtrtri_quark
#define CORE_dtrtri_quark PCORE_dtrtri_quark
#endif
void CORE_dtrtri_quark(Quark *quark)
{
    PLASMA_enum uplo;
    PLASMA_enum diag;
    int N;
    double *A;
    int LDA;
    PLASMA_sequence *sequence;
    PLASMA_request *request;
    int iinfo;

    int info;

    quark_unpack_args_8(quark, uplo, diag, N, A, LDA, sequence, request, iinfo);
    info = LAPACKE_dtrtri_work(
        LAPACK_COL_MAJOR,
        lapack_const(uplo), lapack_const(diag),
        N, A, LDA);
    if ((sequence->status == PLASMA_SUCCESS) && (info > 0))
        plasma_sequence_flush(quark, sequence, request, iinfo + info);
}


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_blas.h>

#define m(i,j) m[i*l+j]

double* calcMys(char m[],int n,int l)
{
  int i,j,k,k1,k2,k3,k4;
  char reslist[]={ 'A', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'K',
      'L', 'M', 'N', 'P', 'Q', 'R', 'S', 'T', 'V', 'W', 'Y'};
  int resnum=sizeof(reslist)/sizeof(char);
  double *simi;
  simi=malloc(n*n*sizeof(double));
  for (i=0;i<n;i++)
  {
    simi[i*n+i]=0.0;
    for (j=i+1;j<n;j++)
    {
      double allcount=0,same=0;
      for (k=0;k<l;k++)
      {
        if (m(i,k)!='-' || m(j,k)!='-')
        {
          allcount++;
          if (m(i,k)==m(j,k))
            same++;
        }
      }
      simi[i*n+j]=same/allcount;
      simi[j*n+i]=simi[i*n+j];
    }
  }
  int *edge,edgenow=0;
  edge=malloc(5*n*2*sizeof(int));
  for (i=0;i<5*n*2;i++)
  {
    edge[i]=-1;
  }
  for (i=0;i<n;i++)
  {
    int has=0,now=0,maxn[5],max5=-1;
    while (has<5)
    {
      if (now!=i)
      {
        maxn[has]=now;
        has++;
      }
      now++;
    }
    for(j=0;j<5;j++)
      if (max5==-1||simi[i*n+maxn[max5]]>simi[i*n+maxn[j]])
        max5=j;
    for (j=now;j<n;j++)
    {
      if (simi[i*n+j]>simi[i*n+maxn[max5]])
      {
        maxn[max5]=j;
        max5=-1;
        for(k=0;k<5;k++)
          if (max5==-1||simi[i*n+maxn[max5]]>simi[i*n+maxn[k]])
            max5=k;
      }
    }
    for (j=0;j<5;j++)
    {
      edge[edgenow*2+0]=i;
      edge[edgenow*2+1]=maxn[j];
      edgenow++;
    }
  }
  free(simi);
  double *mys;
  int n1=0,n2=0;
  mys=malloc(l*l*sizeof(double));
  for (i=0;i<l*l;i++)
    mys[i]=0.0;
  for (i=0;i<l;i++)
  {
    for (j=i+1;j<l;j++)
    {
      double *prop,*pro1,*pro2,allnumber=0,add=0;
      prop=malloc(20*20*20*20*sizeof(double));
      pro1=malloc(20*20*sizeof(double));
      pro2=malloc(20*20*sizeof(double));
      int count1=0,count2=0,count3=0,count4=0;
      for (k=0;k<160000;k++)
        prop[k]=0.0;
      for (k=0;k<400;k++)
      {
        pro1[k]=0;
        pro2[k]=0;
      }

      for (k=0;k<5*n;k++)
      {
        n1=edge[k*2];
        n2=edge[k*2+1];
        if (m(n1,i)!='-' && m(n1,j)!='-' && m(n2,i)!='-' && m(n2,j)!='-')
        {
          count1=0;
          count2=0;
          count3=0;
          count4=0;
          for (k1 = 0; k1 < 20; k1++)
          {
            if (reslist[k1]==m(n1,i))
            {
              count1=k1;
              break;
            }
          }
          for (k1 = 0; k1 < 20; k1++)
          {
            if (reslist[k1]==m(n2,i))
            {
              count2=k1;
              break;
            }
          }
          for (k1 = 0; k1 < 20; k1++)
          {
            if (reslist[k1]==m(n1,j))
            {
              count3=k1;
              break;
            }
          }
          for (k1 = 0; k1 < 20; k1++)
          {
            if (reslist[k1]==m(n2,j))
            {
              count4=k1;
              break;
            }
          }
          prop[(count1*20+count3)*400+(count2*20+count4)]+=1;
          prop[(count2*20+count4)*400+(count1*20+count3)]+=1;
          allnumber+=2;
        }
      }

      gsl_matrix_view m = gsl_matrix_view_array (prop, 400, 400);
      gsl_vector *eval = gsl_vector_alloc (400);
      gsl_matrix *evec = gsl_matrix_alloc (400, 400);
      gsl_eigen_symmv_workspace * w = gsl_eigen_symmv_alloc (400);
      gsl_eigen_symmv (&m.matrix, eval, evec, w);
      gsl_eigen_symmv_free (w);
      gsl_eigen_symmv_sort (eval, evec, GSL_EIGEN_SORT_ABS_ASC);
      gsl_matrix *mid = gsl_matrix_alloc (400, 400);
      gsl_matrix *right = gsl_matrix_alloc (400, 400);
      for (k1 = 0; k1 < 400; k1++)
      {
        double ser=gsl_vector_get(eval,k1);
        for (k2=0;k2<400;k2++)
        {
          gsl_matrix_set(mid,k1,k2,0);
          gsl_matrix_set(right,k1,k2,gsl_matrix_get(evec,k2,k1));
        }
        gsl_matrix_set(mid,k1,k1,ser/(ser+1));
      }
      gsl_matrix *temp= gsl_matrix_alloc (400, 400);
      gsl_blas_dgemm (CblasNoTrans, CblasNoTrans,1.0, evec, mid,
                  0.0, temp);
      gsl_blas_dgemm (CblasNoTrans, CblasNoTrans,1.0, temp, right,
                  0.0, evec);
      gsl_matrix_free(mid);
      gsl_matrix_free(right);
      gsl_vector_free (eval);
      for (k1=0;k1<400;k1++)
        for(k2=0;k2<400;k2++)
        {
          prop[k1*400+k2]=gsl_matrix_get(evec,k1,k2);
          if (prop[k1*400+k2]<0)
            prop[k1*400+k2]=-prop[k1*400+k2];
        }
      gsl_matrix_free (evec);
      gsl_matrix_free (temp);
      for (k1=0;k1<20;k1++)
        for(k2=0;k2<20;k2++)
          for(k3=0;k3<20;k3++)
            for(k4=0;k4<20;k4++)
            {
              pro1[k1*20+k2]+=prop[(k1*20+k3)*400+(k2*20+k4)];
              pro2[k3*20+k4]+=prop[(k1*20+k3)*400+(k2*20+k4)];
            }
      if (allnumber!=0)
      {
        for (k=0;k<160000;k++)
          prop[k]=prop[k]/allnumber;
        for (k=0;k<400;k++)
        {
          pro1[k]=pro1[k]/allnumber;
          pro2[k]=pro2[k]/allnumber;
        }
      }
      allnumber=0.0;
      for (k1=0;k1<20;k1++)
        for (k2=0;k2<20;k2++)
        {
          if (pro1[k1*20+k2]==0)
          {
            continue;
          }
          else
          {
            for (k3=0;k3<20;k3++)
              for (k4=0;k4<20;k4++)
              {
                if (pro2[k3*20+k4]==0||prop[(k1*20+k3)*400+(k2*20+k4)]==0)
                {
                  continue;
                }
                else
                {
                  allnumber+=prop[(k1*20+k3)*400+(k2*20+k4)]*log(prop[(k1*20+k3)*400+(k2*20+k4)]/pro1[k1*20+k2]/pro2[k3*20+k4]);
                }
              }
          }
        }
      mys[i*l+j]=allnumber;
      mys[j*l+i]=allnumber;
      free(prop);
      free(pro1);
      free(pro2);
    }
  }
  free(edge);
  return mys;
}

max_stars_repo_name: wzmao/mbio

max_issues_repo_name: wzmao/mbio

max_forks_repo_name: wzmao/mbio


text:
#ifndef OPENMC_VOLUME_CALC_H
#define OPENMC_VOLUME_CALC_H

#include "openmc/position.h"
#include "openmc/tallies/trigger.h"

#include "pugixml.hpp"
#include "xtensor/xtensor.hpp"

#include <array>
#include <string>
#include <vector>
#include <gsl/gsl>

namespace openmc {

//==============================================================================
// Volume calculation class
//==============================================================================

class VolumeCalculation {

public:
  // Aliases, types
  struct Result {
    std::array<double, 2> volume; //!< Mean/standard deviation of volume
    std::vector<int> nuclides; //!< Index of nuclides
    std::vector<double> atoms; //!< Number of atoms for each nuclide
    std::vector<double> uncertainty; //!< Uncertainty on number of atoms
    int iterations; //!< Number of iterations needed to obtain the results
  }; // Results for a single domain

  // Constructors
  VolumeCalculation(pugi::xml_node node);

  // Methods

  //! \brief Stochastically determine the volume of a set of domains along with the
  //!   average number densities of nuclides within the domain
  //
  //! \return Vector of results for each user-specified domain
  std::vector<Result> execute() const;

  //! \brief Write volume calculation results to HDF5 file
  //
  //! \param[in] filename  Path to HDF5 file to write
  //! \param[in] results   Vector of results for each domain
  void to_hdf5(const std::string& filename, const std::vector<Result>& results) const;

  // Tally filter and map types
  enum class TallyDomain {
    UNIVERSE,
    MATERIAL,
    CELL
  };

  // Data members
  TallyDomain domain_type_; //!< Type of domain (cell, material, etc.)
  size_t n_samples_; //!< Number of samples to use
  double threshold_ {-1.0}; //!< Error threshold for domain volumes
  TriggerMetric trigger_type_ {TriggerMetric::not_active}; //!< Trigger metric for the volume calculation
  Position lower_left_; //!< Lower-left position of bounding box
  Position upper_right_; //!< Upper-right position of bounding box
  std::vector<int> domain_ids_; //!< IDs of domains to find volumes of

private:
  //! \brief Check whether a material has already been hit for a given domain.
  //! If not, add new entries to the vectors
  //
  //! \param[in] i_material Index in global materials vector
  //! \param[in,out] indices Vector of material indices
  //! \param[in,out] hits Number of hits corresponding to each material
  void check_hit(int i_material, std::vector<int>& indices,
    std::vector<int>& hits) const;

};

//==============================================================================
// Global variables
//==============================================================================

namespace model {
  extern std::vector<VolumeCalculation> volume_calcs;
}

//==============================================================================
// Non-member functions
//==============================================================================

void free_memory_volume();

} // namespace openmc

#endif // OPENMC_VOLUME_CALC_H

max_stars_repo_name: Hit-Weixg/openmc

max_issues_repo_name: mehmeturkmen/openmc

max_forks_repo_name: mehmeturkmen/openmc