text:
/* spswap.c
 *
 * Copyright (C) 2014 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_matrix.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spmatrix.h>

#include "avl.c"

/*
gsl_spmatrix_transpose()
  Replace the sparse matrix src by its transpose,
keeping the matrix in the same storage format

Inputs: A - (input/output) sparse matrix to transpose.
*/

int
gsl_spmatrix_transpose(gsl_spmatrix * m)
{
  if (GSL_SPMATRIX_ISTRIPLET(m))
    {
      size_t n;

      /* swap row/column indices */
      for (n = 0; n < m->nz; ++n)
        {
          size_t tmp = m->p[n];
          m->p[n] = m->i[n];
          m->i[n] = tmp;
        }

      /* need to rebuild AVL tree, or element searches won't
       * work correctly with transposed indices */
      gsl_spmatrix_tree_rebuild(m);
    }
  else
    {
      GSL_ERROR("unknown sparse matrix type", GSL_EINVAL);
    }

  /* swap dimensions */
  if (m->size1 != m->size2)
    {
      size_t tmp = m->size1;
      m->size1 = m->size2;
      m->size2 = tmp;
    }

  return GSL_SUCCESS;
}

/*
gsl_spmatrix_transpose2()
  Replace the sparse matrix src by its transpose either by
  swapping its row and column indices if it is in triplet storage,
  or by switching its major if it is in compressed storage.

Inputs: A - (input/output) sparse matrix to transpose.
*/

int
gsl_spmatrix_transpose2(gsl_spmatrix * m)
{
  if (GSL_SPMATRIX_ISTRIPLET(m))
    {
      return gsl_spmatrix_transpose(m);
    }
  else if (GSL_SPMATRIX_ISCCS(m))
    {
      m->sptype = GSL_SPMATRIX_CRS;
    }
  else if (GSL_SPMATRIX_ISCRS(m))
    {
      m->sptype = GSL_SPMATRIX_CCS;
    }
  else
    {
      GSL_ERROR("unknown sparse matrix type", GSL_EINVAL);
    }

  /* swap dimensions */
  if (m->size1 != m->size2)
    {
      size_t tmp = m->size1;
      m->size1 = m->size2;
      m->size2 = tmp;
    }

  return GSL_SUCCESS;
}

int
gsl_spmatrix_transpose_memcpy(gsl_spmatrix *dest, const gsl_spmatrix *src)
{
  const size_t M = src->size1;
  const size_t N = src->size2;

  if (M != dest->size2 || N != dest->size1)
    {
      GSL_ERROR("dimensions of dest must be transpose of src matrix",
                GSL_EBADLEN);
    }
  else if (dest->sptype != src->sptype)
    {
      GSL_ERROR("cannot copy matrices of different storage formats",
                GSL_EINVAL);
    }
  else
    {
      int s = GSL_SUCCESS;
      const size_t nz = src->nz;

      if (dest->nzmax < src->nz)
        {
          s = gsl_spmatrix_realloc(src->nz, dest);
          if (s)
            return s;
        }

      if (GSL_SPMATRIX_ISTRIPLET(src))
        {
          size_t n;
          void *ptr;

          for (n = 0; n < nz; ++n)
            {
              dest->i[n] = src->p[n];
              dest->p[n] = src->i[n];
              dest->data[n] = src->data[n];

              /* copy binary tree data */
              ptr = avl_insert(dest->tree_data->tree, &dest->data[n]);
              if (ptr != NULL)
                {
                  GSL_ERROR("detected duplicate entry", GSL_EINVAL);
                }
            }
        }
      else if (GSL_SPMATRIX_ISCCS(src))
        {
          size_t *Ai = src->i;
          size_t *Ap = src->p;
          double *Ad = src->data;
          size_t *ATi = dest->i;
          size_t *ATp = dest->p;
          double *ATd = dest->data;
          size_t *w = (size_t *) dest->work;
          size_t p, j;

          /* initialize to 0 */
          for (p = 0; p < M + 1; ++p)
            ATp[p] = 0;

          /* compute row counts of A (= column counts for A^T) */
          for (p = 0; p < nz; ++p)
            ATp[Ai[p]]++;

          /* compute row pointers for A (= column pointers for A^T) */
          gsl_spmatrix_cumsum(M, ATp);

          /* make copy of row pointers */
          for (j = 0; j < M; ++j)
            w[j] = ATp[j];

          for (j = 0; j < N; ++j)
            {
              for (p = Ap[j]; p < Ap[j + 1]; ++p)
                {
                  size_t k = w[Ai[p]]++;
                  ATi[k] = j;
                  ATd[k] = Ad[p];
                }
            }
        }
      else if (GSL_SPMATRIX_ISCRS(src))
        {
          size_t *Aj = src->i;
          size_t *Ap = src->p;
          double *Ad = src->data;
          size_t *ATj = dest->i;
          size_t *ATp = dest->p;
          double *ATd = dest->data;
          size_t *w = (size_t *) dest->work;
          size_t p, i;

          /* initialize to 0 */
          for (p = 0; p < N + 1; ++p)
            ATp[p] = 0;

          /* compute column counts of A (= row counts for A^T) */
          for (p = 0; p < nz; ++p)
            ATp[Aj[p]]++;

          /* compute column pointers for A (= row pointers for A^T) */
          gsl_spmatrix_cumsum(N, ATp);

          /* make copy of column pointers */
          for (i = 0; i < N; ++i)
            w[i] = ATp[i];

          for (i = 0; i < M; ++i)
            {
              for (p = Ap[i]; p < Ap[i + 1]; ++p)
                {
                  size_t k = w[Aj[p]]++;
                  ATj[k] = i;
                  ATd[k] = Ad[p];
                }
            }
        }
      else
        {
          GSL_ERROR("unknown sparse matrix type", GSL_EINVAL);
        }

      dest->nz = nz;

      return s;
    }
} /* gsl_spmatrix_transpose_memcpy() */


text:
/*
 *  Copyright (C) 2021 FISCO BCOS.
 *  SPDX-License-Identifier: Apache-2.0
 *  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.
 *
 * @brief interface for BlockHeader
 * @file BlockHeader.h
 * @author: yujiechen
 * @date: 2021-03-22
 */
#pragma once
#include "Exceptions.h"
#include "ProtocolTypeDef.h"
#include <bcos-crypto/interfaces/crypto/CryptoSuite.h>
#include <bcos-utilities/DataConvertUtility.h>
#include <gsl/span>

namespace bcos
{
namespace protocol
{
class BlockHeader
{
public:
    using Ptr = std::shared_ptr<BlockHeader>;
    using ConstPtr = std::shared_ptr<const BlockHeader>;
    using BlockHeadersPtr = std::shared_ptr<std::vector<BlockHeader::Ptr> >;
    explicit BlockHeader(bcos::crypto::CryptoSuite::Ptr _cryptoSuite) : m_cryptoSuite(_cryptoSuite)
    {}

    virtual ~BlockHeader() {}

    virtual void decode(bytesConstRef _data) = 0;
    virtual void encode(bytes& _encodeData) const = 0;

    virtual bcos::crypto::HashType hash() const { return {}; }

    virtual void populateFromParents(BlockHeadersPtr _parents, BlockNumber _number)
    {
        // set parentInfo
        ParentInfoList parentInfoList;
        for (auto parentHeader : *_parents)
        {
            ParentInfo parentInfo;
            parentInfo.blockNumber = parentHeader->number();
            parentInfo.blockHash = parentHeader->hash();
            parentInfoList.emplace_back(parentInfo);
        }
        setParentInfo(std::move(parentInfoList));
        setNumber(_number);
    }

    virtual void clear() = 0;

    // verifySignatureList verifys the signatureList
    virtual void verifySignatureList() const
    {
        auto signatures = signatureList();
        auto sealers = sealerList();
        if (signatures.size() < sealers.size())
        {
            BOOST_THROW_EXCEPTION(InvalidBlockHeader() << errinfo_comment(
                                      "Invalid blockHeader for the size of sealerList "
                                      "is smaller than the size of signatureList"));
        }
        for (auto signature : signatures)
        {
            auto sealerIndex = signature.index;
            auto signatureData = signature.signature;
            if (!m_cryptoSuite->signatureImpl()->verify(
                    std::shared_ptr<const bytes>(&((sealers)[sealerIndex]), [](const bytes*) {}),
                    hash(), bytesConstRef(signatureData.data(), signatureData.size())))
            {
                BOOST_THROW_EXCEPTION(
                    InvalidSignatureList()
                    << errinfo_comment("Invalid signatureList for verify failed, signatureData:" +
                                       *toHexString(signatureData)));
            }
        }
    }
    virtual void populateEmptyBlock(
        BlockNumber _number, int64_t _sealerId, int64_t _timestamp = utcTime())
    {
        setNumber(_number);
        setSealer(_sealerId);
        setTimestamp(_timestamp);
    }

    // version returns the version of the blockHeader
    virtual int32_t version() const = 0;
    // parentInfo returns the parent information, including (parentBlockNumber, parentHash)
    virtual gsl::span<const ParentInfo> parentInfo() const = 0;
    // txsRoot returns the txsRoot of the current block
    virtual bcos::crypto::HashType txsRoot() const = 0;
    // receiptsRoot returns the receiptsRoot of the current block
    virtual bcos::crypto::HashType receiptsRoot() const = 0;
    // stateRoot returns the stateRoot of the current block
    virtual bcos::crypto::HashType stateRoot() const = 0;
    // number returns the number of the current block
    virtual BlockNumber number() const = 0;
    virtual u256 gasUsed() const = 0;
    virtual int64_t timestamp() const = 0;
    // sealer returns the sealer that generate this block
    virtual int64_t sealer() const = 0;
    // sealerList returns the current sealer list
    virtual gsl::span<const bytes> sealerList() const = 0;
    virtual bytesConstRef extraData() const = 0;
    virtual gsl::span<const Signature> signatureList() const = 0;
    virtual gsl::span<const uint64_t> consensusWeights() const = 0;

    virtual void setVersion(int32_t _version) = 0;
    virtual void setParentInfo(gsl::span<const ParentInfo> const& _parentInfo) = 0;
    virtual void setParentInfo(ParentInfoList&& _parentInfo) = 0;

    virtual void setTxsRoot(bcos::crypto::HashType _txsRoot) = 0;
    virtual void setReceiptsRoot(bcos::crypto::HashType _receiptsRoot) = 0;
    virtual void setStateRoot(bcos::crypto::HashType _stateRoot) = 0;
    virtual void setNumber(BlockNumber _blockNumber) = 0;
    virtual void setGasUsed(u256 _gasUsed) = 0;
    virtual void setTimestamp(int64_t _timestamp) = 0;
    virtual void setSealer(int64_t _sealerId) = 0;

    virtual void setSealerList(gsl::span<const bytes> const& _sealerList) = 0;
    virtual void setSealerList(std::vector<bytes>&& _sealerList) = 0;

    virtual void setConsensusWeights(gsl::span<const uint64_t> const& _weightList) = 0;
    virtual void setConsensusWeights(std::vector<uint64_t>&& _weightList) = 0;

    virtual void setExtraData(bytes const& _extraData) = 0;
    virtual void setExtraData(bytes&& _extraData) = 0;

    virtual void setSignatureList(gsl::span<const Signature> const& _signatureList) = 0;
    virtual void setSignatureList(SignatureList&& _signatureList) = 0;
    virtual bcos::crypto::CryptoSuite::Ptr cryptoSuite() { return m_cryptoSuite; }

protected:
    bcos::crypto::CryptoSuite::Ptr m_cryptoSuite;
};
}  // namespace protocol
}  // namespace bcos

ext: h

max_stars_repo_name: contropist/FISCO-BCOS

max_issues_repo_name: contropist/FISCO-BCOS

max_forks_repo_name: contropist/FISCO-BCOS


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

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

double
gsl_histogram_max (const gsl_histogram * h)
{
  const int n = h->n;

  return h->range[n];
}

double
gsl_histogram_min (const gsl_histogram * h)
{
  return h->range[0];
}

size_t
gsl_histogram_bins (const gsl_histogram * h)
{
  return h->n;
}

ext: c

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


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

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

#ifdef MPISENDRECV_CHECKSUM

#undef MPI_Sendrecv

int MPI_Check_Sendrecv(void *sendbuf, int sendcount, MPI_Datatype sendtype,
		       int dest, int sendtag, void *recvbufreal, int recvcount,
		       MPI_Datatype recvtype, int source, int recvtag, MPI_Comm comm, MPI_Status * status)
{
  int checksumtag = 1000, errtag = 2000;
  int i, iter = 0, err_flag, err_flag_imported, size_sendtype, size_recvtype;
  long long sendCheckSum, recvCheckSum, importedCheckSum;
  unsigned char *p, *buf, *recvbuf;

  if(dest != source)
    endrun(3);

  MPI_Type_size(sendtype, &size_sendtype);
  MPI_Type_size(recvtype, &size_recvtype);

  if(dest == ThisTask)
    {
      memcpy(recvbufreal, sendbuf, recvcount * size_recvtype);
      return 0;
    }

  if(!(buf = mymalloc(recvcount * size_recvtype + 1024)))
    endrun(6);

  for(i = 0, p = buf; i < recvcount * size_recvtype + 1024; i++)
    *p++ = 255;

  recvbuf = buf + 512;

  MPI_Sendrecv(sendbuf, sendcount, sendtype, dest, sendtag,
	       recvbuf, recvcount, recvtype, source, recvtag, comm, status);

  for(i = 0, p = buf; i < 512; i++, p++)
    {
      if(*p != 255)
	{
	  printf
	    ("MPI-ERROR: Task=%d/%s: Recv occured before recv buffer. message-size=%d from %d, i=%d c=%d\n",
	     ThisTask, getenv("HOST"), recvcount, dest, i, *p);
	  fflush(stdout);
	  endrun(6);
	}
    }

  for(i = 0, p = recvbuf + recvcount * size_recvtype; i < 512; i++, p++)
    {
      if(*p != 255)
	{
	  printf
	    ("MPI-ERROR: Task=%d/%s: Recv occured after recv buffer. message-size=%d from %d, i=%d c=%d\n",
	     ThisTask, getenv("HOST"), recvcount, dest, i, *p);
	  fflush(stdout);
	  endrun(6);
	}
    }

  for(i = 0, p = sendbuf, sendCheckSum = 0; i < sendcount * size_sendtype; i++, p++)
    sendCheckSum += *p;

  importedCheckSum = 0;

  if(dest > ThisTask)
    {
      if(sendcount > 0)
	MPI_Ssend(&sendCheckSum, sizeof(sendCheckSum), MPI_BYTE, dest, checksumtag, MPI_COMM_WORLD);
      if(recvcount > 0)
	MPI_Recv(&importedCheckSum, sizeof(importedCheckSum), MPI_BYTE, dest, checksumtag, MPI_COMM_WORLD,
		 status);
    }
  else
    {
      if(recvcount > 0)
	MPI_Recv(&importedCheckSum, sizeof(importedCheckSum), MPI_BYTE, dest, checksumtag, MPI_COMM_WORLD,
		 status);
      if(sendcount > 0)
	MPI_Ssend(&sendCheckSum, sizeof(sendCheckSum), MPI_BYTE, dest, checksumtag, MPI_COMM_WORLD);
    }

  checksumtag++;

  for(i = 0, p = recvbuf, recvCheckSum = 0; i < recvcount * size_recvtype; i++, p++)
    recvCheckSum += *p;

  err_flag = err_flag_imported = 0;

  if(recvCheckSum != importedCheckSum)
    {
      printf
	("MPI-ERROR: Receive error on task=%d/%s from task=%d, message size=%d, sendcount=%d checksums= %d %d  %d %d. Try to fix it...\n",
	 ThisTask, getenv("HOST"), source, recvcount, sendcount, (int) (recvCheckSum >> 32),
	 (int) recvCheckSum, (int) (importedCheckSum >> 32), (int) importedCheckSum);
      fflush(stdout);

      err_flag = 1;
    }

  if(dest > ThisTask)
    {
      MPI_Ssend(&err_flag, 1, MPI_INT, dest, errtag, MPI_COMM_WORLD);
      MPI_Recv(&err_flag_imported, 1, MPI_INT, dest, errtag, MPI_COMM_WORLD, status);
    }
  else
    {
      MPI_Recv(&err_flag_imported, 1, MPI_INT, dest, errtag, MPI_COMM_WORLD, status);
      MPI_Ssend(&err_flag, 1, MPI_INT, dest, errtag, MPI_COMM_WORLD);
    }
  errtag++;

  if(err_flag > 0 || err_flag_imported > 0)
    {
      printf("Task=%d is on %s, wants to send %d and has checksum=%d %d of send data\n",
	     ThisTask, getenv("HOST"), sendcount, (int) (sendCheckSum >> 32), (int) sendCheckSum);
      fflush(stdout);

      do
	{
	  sendtag++;
	  recvtag++;

	  for(i = 0, p = recvbuf; i < recvcount * size_recvtype; i++, p++)
	    *p = 0;

	  if((iter & 1) == 0)
	    {
	      if(dest > ThisTask)
		{
		  if(sendcount > 0)
		    MPI_Ssend(sendbuf, sendcount, sendtype, dest, sendtag, MPI_COMM_WORLD);
		  if(recvcount > 0)
		    MPI_Recv(recvbuf, recvcount, recvtype, dest, recvtag, MPI_COMM_WORLD, status);
		}
	      else
		{
		  if(recvcount > 0)
		    MPI_Recv(recvbuf, recvcount, recvtype, dest, recvtag, MPI_COMM_WORLD, status);
		  if(sendcount > 0)
		    MPI_Ssend(sendbuf, sendcount, sendtype, dest, sendtag, MPI_COMM_WORLD);
		}
	    }
	  else
	    {
	      if(iter > 5)
		{
		  printf("we're trying to send each byte now on task=%d (iter=%d)\n", ThisTask, iter);
		  fflush(stdout);
		  if(dest > ThisTask)
		    {
		      for(i = 0, p = sendbuf; i < sendcount * size_sendtype; i++, p++)
			MPI_Ssend(p, 1, MPI_BYTE, dest, i, MPI_COMM_WORLD);
		      for(i = 0, p = recvbuf; i < recvcount * size_recvtype; i++, p++)
			MPI_Recv(p, 1, MPI_BYTE, dest, i, MPI_COMM_WORLD, status);
		    }
		  else
		    {
		      for(i = 0, p = recvbuf; i < recvcount * size_recvtype; i++, p++)
			MPI_Recv(p, 1, MPI_BYTE, dest, i, MPI_COMM_WORLD, status);
		      for(i = 0, p = sendbuf; i < sendcount * size_sendtype; i++, p++)
			MPI_Ssend(p, 1, MPI_BYTE, dest, i, MPI_COMM_WORLD);
		    }
		}
	      else
		{
		  MPI_Sendrecv(sendbuf, sendcount, sendtype, dest, sendtag,
			       recvbuf, recvcount, recvtype, source, recvtag, comm, status);
		}
	    }

	  importedCheckSum = 0;

	  for(i = 0, p = sendbuf, sendCheckSum = 0; i < sendcount * size_sendtype; i++, p++)
	    sendCheckSum += *p;

	  printf("Task=%d gas send_checksum=%d %d\n", ThisTask, (int) (sendCheckSum >> 32),
		 (int) sendCheckSum);
	  fflush(stdout);

	  if(dest > ThisTask)
	    {
	      if(sendcount > 0)
		MPI_Ssend(&sendCheckSum, sizeof(sendCheckSum), MPI_BYTE, dest, checksumtag, MPI_COMM_WORLD);
	      if(recvcount > 0)
		MPI_Recv(&importedCheckSum, sizeof(importedCheckSum), MPI_BYTE, dest, checksumtag,
			 MPI_COMM_WORLD, status);
	    }
	  else
	    {
	      if(recvcount > 0)
		MPI_Recv(&importedCheckSum, sizeof(importedCheckSum), MPI_BYTE, dest, checksumtag,
			 MPI_COMM_WORLD, status);
	      if(sendcount > 0)
		MPI_Ssend(&sendCheckSum, sizeof(sendCheckSum), MPI_BYTE, dest, checksumtag, MPI_COMM_WORLD);
	    }

	  for(i = 0, p = recvbuf, recvCheckSum = 0; i < recvcount; i++, p++)
	    recvCheckSum += *p;

	  err_flag = err_flag_imported = 0;

	  if(recvCheckSum != importedCheckSum)
	    {
	      printf
		("MPI-ERROR: Again (iter=%d) a receive error on task=%d/%s from task=%d, message size=%d, checksums= %d %d  %d %d. Try to fix it...\n",
		 iter, ThisTask, getenv("HOST"), source, recvcount, (int) (recvCheckSum >> 32),
		 (int) recvCheckSum, (int) (importedCheckSum >> 32), (int) importedCheckSum);
	      fflush(stdout);
	      err_flag = 1;
	    }

	  if(dest > ThisTask)
	    {
	      MPI_Ssend(&err_flag, 1, MPI_INT, dest, errtag, MPI_COMM_WORLD);
	      MPI_Recv(&err_flag_imported, 1, MPI_INT, dest, errtag, MPI_COMM_WORLD, status);
	    }
	  else
	    {
	      MPI_Recv(&err_flag_imported, 1, MPI_INT, dest, errtag, MPI_COMM_WORLD, status);
	      MPI_Ssend(&err_flag, 1, MPI_INT, dest, errtag, MPI_COMM_WORLD);
	    }

	  if(err_flag == 0 && err_flag_imported == 0)
	    break;

	  errtag++;
	  checksumtag++;
	  iter++;
	}
      while(iter < 10);

      if(iter >= 10)
	{
	  char buf[1000];
	  int length;
	  FILE *fd;

	  sprintf(buf, "send_data_%d.dat", ThisTask);
	  fd = fopen(buf, "w");
	  length = sendcount * size_sendtype;
	  fwrite(&length, 1, sizeof(int), fd);
	  fwrite(sendbuf, sendcount, size_sendtype, fd);
	  fclose(fd);

	  sprintf(buf, "recv_data_%d.dat", ThisTask);
	  fd = fopen(buf, "w");
	  length = recvcount * size_recvtype;
	  fwrite(&length, 1, sizeof(int), fd);
	  fwrite(recvbuf, recvcount, size_recvtype, fd);
	  fclose(fd);

	  printf("MPI-ERROR: Even 10 trials proved to be insufficient on task=%d/%s. Stopping\n", ThisTask,
		 getenv("HOST"));
	  fflush(stdout);
	  endrun(10);
	}
    }

  memcpy(recvbufreal, recvbuf, recvcount * size_recvtype);

  myfree(buf);

  return 0;
}

#endif

max_stars_repo_name: egpbos/egp

max_issues_repo_name: egpbos/egp

max_forks_repo_name: egpbos/egp


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

#ifndef PICTOLEV_TILES_H
#define PICTOLEV_TILES_H

#include <algorithm>
#include <cstddef>
#include <iterator>
#include <type_traits>
#include <vector>
#include <gsl/gsl_assert>
#include <gsl/gsl_util>
#include <gsl/span>
#include "grid_size.h"

template<class T>
using image_fragment = std::vector<gsl::span<T>>;

template<class T>
auto buffer_to_image(gsl::span<T> buffer, grid_size size) {
	Expects(gsl::narrow_cast<std::size_t>(buffer.size()) == grid_area(size));
	image_fragment<T> image(size.height);
	std::generate(image.begin(), image.end(), [buffer, width = size.width]() mutable {
		gsl::span<T> row = buffer.subspan(0, width);
		buffer = buffer.subspan(width);
		return row;
	});
	return image;
}

template<class T>
using tile_vector = std::vector<image_fragment<T>>;

template<class ForwardIt>
auto image_to_tile_list(ForwardIt first, ForwardIt last, std::size_t tile_size) {
	Expects(std::distance(first, last) % tile_size == 0);
	Expects(first == last || first->size() % tile_size == 0);
	using container_type = std::remove_reference_t<typename ForwardIt::reference>;
	using value_type = std::remove_reference_t<typename container_type::reference>;
	const grid_size size {
		first != last ? first->size() / tile_size : 0,
		std::distance(first, last) / tile_size,
	};
	tile_vector<value_type> tiles(grid_area(size), image_fragment<value_type>(tile_size));
	auto out = tiles.begin();
	while (first != last) {
		ForwardIt second = first;
		std::advance(second, tile_size);
		for (gsl::index x = 0; x != first->size(); x += tile_size) {
			std::transform(first, second, out->begin(), [=](const auto& row) {
				return gsl::span<value_type>(row).subspan(x, tile_size);
			});
			++out;
		}
		first = second;
	}
	return tiles;
}

#endif

max_stars_repo_name: SirEmentaler/PicToLev

max_issues_repo_name: SirEmentaler/PicToLev

max_forks_repo_name: SirEmentaler/PicToLev


text:
/* Implementation for Gillespie's direct method */

#include <math.h>
#include <stdint.h>
#include <time.h>
#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

#include "ssa.h"
#include "dm.h"

const gsl_rng_type *SSARNGT;
gsl_rng *SSARNG;

INDEX nextreaction(double *propensities, INDEX m)
{
	double r1 = gsl_ran_flat(SSARNG, 0.0, 1.0) * propensities[m - 1];
	INDEX i = 0;
	for(i = 0; i < m; i++)
		if(r1 <= propensities[i])
			return i;
	return -1;
}

/*
Perform one step of the Gillespie direct method, selecting a
reaction to perform and altering molecule counts and time
accordingly.
*/
void ssa_dmstep(INDEX n, INDEX m, double *propensities, INDEX **R, INDEX **P, double *k, COUNT *x, double *t)
{
	propensities[0] = ssa_h(R[0], x, n) * k[0];
	INDEX i;
	for (i = 1; i < m; i++)
		propensities[i] = ssa_h(R[i], x, n) * k[i] + propensities[i - 1];

	if (propensities[m - 1] > 0.0) {
		*t += gsl_ran_exponential(SSARNG, 1 / propensities[m - 1]);
		INDEX next = nextreaction(propensities, m);
		ssa_doreaction(R[next], P[next], x, n);
	}
}

void ssa_dm(INDEX **R, INDEX **P, INDEX n, INDEX m, double *k, COUNT *x, double T)
{
	gsl_rng_env_setup();
	SSARNGT = gsl_rng_default;
	SSARNG = gsl_rng_alloc(SSARNGT);
	gsl_rng_set(SSARNG, time(NULL));

        double t = 0.0;
        double *propensities = malloc(sizeof(double) * m);

	do {
		ssa_dmstep(n, m, propensities, R, P, k, x, &t);
		ssa_printstate(t, x, n);
	} while(t < T && propensities[m - 1] > 0.0);

	free(propensities);
}

max_stars_repo_name: lgrozinger/ssapy

max_issues_repo_name: lgrozinger/ssapy

max_forks_repo_name: lgrozinger/ssapy


text:
/* movstat/test_minmax.c
 *
 * Copyright (C) 2018 Patrick Alken
 *
 * 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 <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_movstat.h>

/* compute filtered data by explicitely constructing window and finding min/max */
int
slow_minmax(const gsl_movstat_end_t etype, const gsl_vector * x, gsl_vector * y_min, gsl_vector * y_max,
            const int H, const int J)
{
  const size_t n = x->size;
  const int K = H + J + 1;
  double *window = malloc(K * sizeof(double));
  size_t i;

  for (i = 0; i < n; ++i)
    {
      size_t wsize = gsl_movstat_fill(etype, x, i, H, J, window);
      gsl_vector_view v = gsl_vector_view_array(window, wsize);
      double min, max;

      gsl_vector_minmax(&v.vector, &min, &max);
      gsl_vector_set(y_min, i, min);
      gsl_vector_set(y_max, i, max);
    }

  free(window);

  return GSL_SUCCESS;
}

static double
func_min(const size_t n, double x[], void * params)
{
  gsl_vector_view v = gsl_vector_view_array(x, n);
  (void) params;
  return gsl_vector_min(&v.vector);
}

static double
func_max(const size_t n, double x[], void * params)
{
  gsl_vector_view v = gsl_vector_view_array(x, n);
  (void) params;
  return gsl_vector_max(&v.vector);
}

static void
test_minmax_x(const double tol, const gsl_vector * x, const int H, const int J,
              const gsl_movstat_end_t endtype, const char * desc)
{
  const size_t n = x->size;
  gsl_vector * u_min = gsl_vector_alloc(n);
  gsl_vector * y_min = gsl_vector_alloc(n);
  gsl_vector * z_min = gsl_vector_alloc(n);
  gsl_vector * u_max = gsl_vector_alloc(n);
  gsl_vector * y_max = gsl_vector_alloc(n);
  gsl_vector * z_max = gsl_vector_alloc(n);
  gsl_movstat_workspace * w = gsl_movstat_alloc2(H, J);
  gsl_movstat_function F1, F2;
  char buf[2048];

  F1.function = func_min;
  F1.params = NULL;

  F2.function = func_max;
  F2.params = NULL;

  /* compute moving min/max */
  gsl_movstat_min(endtype, x, u_min, w);
  gsl_movstat_max(endtype, x, u_max, w);
  gsl_movstat_minmax(endtype, x, y_min, y_max, w);

  /* compute moving min/max with slow brute force method */
  slow_minmax(endtype, x, z_min, z_max, H, J);

  sprintf(buf, "test_minmax: %s min endtype=%d n=%zu H=%d J=%d", desc, endtype, n, H, J);
  compare_vectors(tol, u_min, z_min, buf);

  sprintf(buf, "test_minmax: %s max endtype=%d n=%zu H=%d J=%d", desc, endtype, n, H, J);
  compare_vectors(tol, u_max, z_max, buf);

  sprintf(buf, "test_minmax: %s minmax(minimum) endtype=%d n=%zu H=%d J=%d", desc, endtype, n, H, J);
  compare_vectors(tol, y_min, z_min, buf);

  sprintf(buf, "test_minmax: %s minmax(maximum) endtype=%d n=%zu H=%d J=%d", desc, endtype, n, H, J);
  compare_vectors(tol, y_max, z_max, buf);

  /* in-place tests */

  gsl_vector_memcpy(u_min, x);
  gsl_vector_memcpy(u_max, x);

  gsl_movstat_min(endtype, u_min, u_min, w);
  gsl_movstat_max(endtype, u_max, u_max, w);

  sprintf(buf, "test_minmax: %s in-place min endtype=%d n=%zu H=%d J=%d", desc, endtype, n, H, J);
  compare_vectors(tol, u_min, z_min, buf);

  sprintf(buf, "test_minmax: %s in-place max endtype=%d n=%zu H=%d J=%d", desc, endtype, n, H, J);
  compare_vectors(tol, u_max, z_max, buf);

  /* user-defined function tests */

  gsl_movstat_apply(endtype, &F1, x, z_min, w);
  sprintf(buf, "n=%zu H=%d J=%d endtype=%u min user", n, H, J, endtype);
  compare_vectors(tol, z_min, y_min, buf);

  gsl_movstat_apply(endtype, &F2, x, z_max, w);
  sprintf(buf, "n=%zu H=%d J=%d endtype=%u max user", n, H, J, endtype);
  compare_vectors(tol, z_max, y_max, buf);

  gsl_vector_free(u_min);
  gsl_vector_free(y_min);
  gsl_vector_free(z_min);
  gsl_vector_free(u_max);
  gsl_vector_free(y_max);
  gsl_vector_free(z_max);
  gsl_movstat_free(w);
}

/* test alternating sequence [a,b,a,b,...] input */
static void
test_minmax_alt(const double tol, const size_t n, const int H, const int J,
                const gsl_movstat_end_t endtype)
{
  const double a = 5.0;
  const double b = -5.0;
  gsl_vector * x = gsl_vector_alloc(n);
  size_t i;

  for (i = 0; i < n; ++i)
    {
      if (i % 2 == 0)
        gsl_vector_set(x, i, a);
      else
        gsl_vector_set(x, i, b);
    }

  test_minmax_x(tol, x, H, J, endtype, "alternating");

  gsl_vector_free(x);
}

/* test noisy sine wave input */
static void
test_minmax_sine(const double tol, const size_t n, const int H, const int J,
                 const gsl_movstat_end_t endtype, gsl_rng * rng_p)
{
  gsl_vector * x = gsl_vector_alloc(n);

  /* construct noisy sine signal */
  test_noisy_sine(0.5, x, rng_p);

  test_minmax_x(tol, x, H, J, endtype, "noisy_sine");

  gsl_vector_free(x);
}

/* test random input */
static void
test_minmax_random(const double tol, const size_t n, const int H, const int J,
                   const gsl_movstat_end_t endtype, gsl_rng * rng_p)
{
  gsl_vector * x = gsl_vector_alloc(n);

  /* construct random input signal */
  random_vector(x, rng_p);

  test_minmax_x(tol, x, H, J, endtype, "random");

  gsl_vector_free(x);
}

static void
test_minmax(gsl_rng * rng_p)
{
  /* alternating input */

  test_minmax_alt(GSL_DBL_EPSILON, 1000, 7, 7, GSL_MOVSTAT_END_PADZERO);
  test_minmax_alt(GSL_DBL_EPSILON, 1000, 5, 2, GSL_MOVSTAT_END_PADZERO);
  test_minmax_alt(GSL_DBL_EPSILON, 500, 1, 3, GSL_MOVSTAT_END_PADZERO);
  test_minmax_alt(GSL_DBL_EPSILON, 20, 50, 10, GSL_MOVSTAT_END_PADZERO);
  test_minmax_alt(GSL_DBL_EPSILON, 20, 10, 50, GSL_MOVSTAT_END_PADZERO);

  /* noisy sine wave input */

  test_minmax_sine(GSL_DBL_EPSILON, 1000, 5, 7, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_minmax_sine(GSL_DBL_EPSILON, 2000, 0, 2, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_minmax_sine(GSL_DBL_EPSILON, 500, 3, 0, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_minmax_sine(GSL_DBL_EPSILON, 20, 50, 50, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_minmax_sine(GSL_DBL_EPSILON, 20, 10, 50, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_minmax_sine(GSL_DBL_EPSILON, 20, 50, 10, GSL_MOVSTAT_END_PADZERO, rng_p);

  test_minmax_sine(GSL_DBL_EPSILON, 500, 5, 7, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_minmax_sine(GSL_DBL_EPSILON, 1000, 10, 20, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_minmax_sine(GSL_DBL_EPSILON, 1000, 3, 3, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_minmax_sine(GSL_DBL_EPSILON, 20, 50, 50, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_minmax_sine(GSL_DBL_EPSILON, 20, 10, 50, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_minmax_sine(GSL_DBL_EPSILON, 20, 50, 10, GSL_MOVSTAT_END_PADVALUE, rng_p);

  test_minmax_sine(GSL_DBL_EPSILON, 500, 5, 7, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_minmax_sine(GSL_DBL_EPSILON, 1000, 10, 20, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_minmax_sine(GSL_DBL_EPSILON, 1000, 3, 3, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_minmax_sine(GSL_DBL_EPSILON, 1000, 30, 5, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_minmax_sine(GSL_DBL_EPSILON, 1000, 5, 30, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_minmax_sine(GSL_DBL_EPSILON, 20, 50, 50, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_minmax_sine(GSL_DBL_EPSILON, 20, 10, 50, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_minmax_sine(GSL_DBL_EPSILON, 20, 50, 10, GSL_MOVSTAT_END_TRUNCATE, rng_p);

  /* random input */

  test_minmax_random(GSL_DBL_EPSILON, 1000, 0, 0, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_minmax_random(GSL_DBL_EPSILON, 1000, 5, 7, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_minmax_random(GSL_DBL_EPSILON, 2000, 0, 2, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_minmax_random(GSL_DBL_EPSILON, 500, 3, 0, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_minmax_random(GSL_DBL_EPSILON, 500, 10, 5, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_minmax_random(GSL_DBL_EPSILON, 500, 5, 10, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_minmax_random(GSL_DBL_EPSILON, 20, 50, 50, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_minmax_random(GSL_DBL_EPSILON, 20, 10, 50, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_minmax_random(GSL_DBL_EPSILON, 20, 50, 10, GSL_MOVSTAT_END_PADZERO, rng_p);

  test_minmax_random(GSL_DBL_EPSILON, 1000, 0, 0, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_minmax_random(GSL_DBL_EPSILON, 1000, 5, 7, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_minmax_random(GSL_DBL_EPSILON, 2000, 0, 2, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_minmax_random(GSL_DBL_EPSILON, 500, 3, 0, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_minmax_random(GSL_DBL_EPSILON, 500, 10, 5, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_minmax_random(GSL_DBL_EPSILON, 500, 5, 10, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_minmax_random(GSL_DBL_EPSILON, 20, 50, 50, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_minmax_random(GSL_DBL_EPSILON, 20, 10, 50, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_minmax_random(GSL_DBL_EPSILON, 20, 50, 10, GSL_MOVSTAT_END_PADVALUE, rng_p);

  test_minmax_random(GSL_DBL_EPSILON, 1000, 0, 0, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_minmax_random(GSL_DBL_EPSILON, 1000, 5, 7, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_minmax_random(GSL_DBL_EPSILON, 2000, 0, 2, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_minmax_random(GSL_DBL_EPSILON, 500, 3, 0, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_minmax_random(GSL_DBL_EPSILON, 500, 10, 5, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_minmax_random(GSL_DBL_EPSILON, 500, 5, 10, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_minmax_random(GSL_DBL_EPSILON, 20, 50, 50, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_minmax_random(GSL_DBL_EPSILON, 20, 10, 50, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_minmax_random(GSL_DBL_EPSILON, 20, 50, 10, GSL_MOVSTAT_END_TRUNCATE, rng_p);
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
#include <stdio.h>
#include <gsl/gsl_cdf.h>

int
main (void)
{
  double P, Q;
  double x = 2.0;

  P = gsl_cdf_ugaussian_P (x);
  printf ("prob(x < %f) = %f\n", x, P);

  Q = gsl_cdf_ugaussian_Q (x);
  printf ("prob(x > %f) = %f\n", x, Q);

  x = gsl_cdf_ugaussian_Pinv (P);
  printf ("Pinv(%f) = %f\n", P, x);

  x = gsl_cdf_ugaussian_Qinv (Q);
  printf ("Qinv(%f) = %f\n", Q, x);

  return 0;
}

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/* eigen/genv.c
 *
 * Copyright (C) 2007 Patrick Alken
 *
 * 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 <stdlib.h>
#include <math.h>

#include <config.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_vector_complex.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_errno.h>

/*
 * This module computes the eigenvalues and eigenvectors of a
 * real generalized eigensystem A x = \lambda B x. Left and right
 * Schur vectors are optionally computed as well.
 *
 * This file contains routines based on original code from LAPACK
 * which is distributed under the modified BSD license.
 */

static int genv_get_right_eigenvectors(const gsl_matrix *S,
                                       const gsl_matrix *T,
                                       gsl_matrix *Z,
                                       gsl_matrix_complex *evec,
                                       gsl_eigen_genv_workspace *w);
static void genv_normalize_eigenvectors(gsl_vector_complex *alpha,
                                        gsl_matrix_complex *evec);

/*
gsl_eigen_genv_alloc()
  Allocate a workspace for solving the generalized eigenvalue problem.
The size of this workspace is O(7n).

Inputs: n - size of matrices

Return: pointer to workspace
*/

gsl_eigen_genv_workspace *
gsl_eigen_genv_alloc(const size_t n)
{
  gsl_eigen_genv_workspace *w;

  if (n == 0)
    {
      GSL_ERROR_NULL ("matrix dimension must be positive integer",
                      GSL_EINVAL);
    }

  w = (gsl_eigen_genv_workspace *) calloc (1, sizeof (gsl_eigen_genv_workspace));

  if (w == 0)
    {
      GSL_ERROR_NULL ("failed to allocate space for workspace", GSL_ENOMEM);
    }

  w->size = n;
  w->Q = NULL;
  w->Z = NULL;

  w->gen_workspace_p = gsl_eigen_gen_alloc(n);

  if (w->gen_workspace_p == 0)
    {
      gsl_eigen_genv_free(w);
      GSL_ERROR_NULL ("failed to allocate space for gen workspace", GSL_ENOMEM);
    }

  /* compute the full Schur forms */
  gsl_eigen_gen_params(1, 1, 1, w->gen_workspace_p);

  w->work1 = gsl_vector_alloc(n);
  w->work2 = gsl_vector_alloc(n);
  w->work3 = gsl_vector_alloc(n);
  w->work4 = gsl_vector_alloc(n);
  w->work5 = gsl_vector_alloc(n);
  w->work6 = gsl_vector_alloc(n);

  if (w->work1 == 0 || w->work2 == 0 || w->work3 == 0 ||
      w->work4 == 0 || w->work5 == 0 || w->work6 == 0)
    {
      gsl_eigen_genv_free(w);
      GSL_ERROR_NULL ("failed to allocate space for additional workspace", GSL_ENOMEM);
    }

  return (w);
} /* gsl_eigen_genv_alloc() */

/*
gsl_eigen_genv_free()
  Free workspace w
*/

void
gsl_eigen_genv_free(gsl_eigen_genv_workspace *w)
{
  RETURN_IF_NULL (w);

  if (w->gen_workspace_p)
    gsl_eigen_gen_free(w->gen_workspace_p);

  if (w->work1)
    gsl_vector_free(w->work1);

  if (w->work2)
    gsl_vector_free(w->work2);

  if (w->work3)
    gsl_vector_free(w->work3);

  if (w->work4)
    gsl_vector_free(w->work4);

  if (w->work5)
    gsl_vector_free(w->work5);

  if (w->work6)
    gsl_vector_free(w->work6);

  free(w);
} /* gsl_eigen_genv_free() */

/*
gsl_eigen_genv()

Solve the generalized eigenvalue problem

A x = \lambda B x

for the eigenvalues \lambda and right eigenvectors x.

Inputs: A     - general real matrix
        B     - general real matrix
        alpha - (output) where to store eigenvalue numerators
        beta  - (output) where to store eigenvalue denominators
        evec  - (output) where to store eigenvectors
        w     - workspace

Return: success or error
*/

int
gsl_eigen_genv (gsl_matrix * A, gsl_matrix * B, gsl_vector_complex * alpha,
                gsl_vector * beta, gsl_matrix_complex *evec,
                gsl_eigen_genv_workspace * w)
{
  const size_t N = A->size1;

  /* check matrix and vector sizes */

  if (N != A->size2)
    {
      GSL_ERROR ("matrix must be square to compute eigenvalues", GSL_ENOTSQR);
    }
  else if ((N != B->size1) || (N != B->size2))
    {
      GSL_ERROR ("B matrix dimensions must match A", GSL_EBADLEN);
    }
  else if (alpha->size != N || beta->size != N)
    {
      GSL_ERROR ("eigenvalue vector must match matrix size", GSL_EBADLEN);
    }
  else if (w->size != N)
    {
      GSL_ERROR ("matrix size does not match workspace", GSL_EBADLEN);
    }
  else if (evec->size1 != N)
    {
      GSL_ERROR ("eigenvector matrix has wrong size", GSL_EBADLEN);
    }
  else
    {
      int s;
      gsl_matrix Z;

      /*
       * We need a place to store the right Schur vectors, so we will
       * treat evec as a real matrix and store them in the left
       * half - the factor of 2 in the tda corresponds to the
       * complex multiplicity
       */
      Z.size1 = N;
      Z.size2 = N;
      Z.tda = 2 * N;
      Z.data = evec->data;
      Z.block = 0;
      Z.owner = 0;

      s = gsl_eigen_gen_QZ(A, B, alpha, beta, w->Q, &Z, w->gen_workspace_p);

      if (w->Z)
        {
          /* save right Schur vectors */
          gsl_matrix_memcpy(w->Z, &Z);
        }

      /* only compute eigenvectors if we found all eigenvalues */
      if (s == GSL_SUCCESS)
        {
          /* compute eigenvectors */
          s = genv_get_right_eigenvectors(A, B, &Z, evec, w);

          if (s == GSL_SUCCESS)
            genv_normalize_eigenvectors(alpha, evec);
        }

      return s;
    }
} /* gsl_eigen_genv() */

/*
gsl_eigen_genv_QZ()

Solve the generalized eigenvalue problem

A x = \lambda B x

for the eigenvalues \lambda and right eigenvectors x. Optionally
compute left and/or right Schur vectors Q and Z which satisfy:

A = Q S Z^t
B = Q T Z^t

where (S, T) is the generalized Schur form of (A, B)

Inputs: A     - general real matrix
        B     - general real matrix
        alpha - (output) where to store eigenvalue numerators
        beta  - (output) where to store eigenvalue denominators
        evec  - (output) where to store eigenvectors
        Q     - (output) if non-null, where to store left Schur vectors
        Z     - (output) if non-null, where to store right Schur vectors
        w     - workspace

Return: success or error
*/

int
gsl_eigen_genv_QZ (gsl_matrix * A, gsl_matrix * B,
                   gsl_vector_complex * alpha, gsl_vector * beta,
                   gsl_matrix_complex * evec,
                   gsl_matrix * Q, gsl_matrix * Z,
                   gsl_eigen_genv_workspace * w)
{
  if (Q && (A->size1 != Q->size1 || A->size1 != Q->size2))
    {
      GSL_ERROR("Q matrix has wrong dimensions", GSL_EBADLEN);
    }
  else if (Z && (A->size1 != Z->size1 || A->size1 != Z->size2))
    {
      GSL_ERROR("Z matrix has wrong dimensions", GSL_EBADLEN);
    }
  else
    {
      int s;

      w->Q = Q;
      w->Z = Z;

      s = gsl_eigen_genv(A, B, alpha, beta, evec, w);

      w->Q = NULL;
      w->Z = NULL;

      return s;
    }
} /* gsl_eigen_genv_QZ() */

/********************************************
 *           INTERNAL ROUTINES              *
 ********************************************/

/*
genv_get_right_eigenvectors()
  Compute right eigenvectors of the Schur form (S, T) and then
backtransform them using the right Schur vectors to get right
eigenvectors of the original system.

Inputs: S     - upper quasi-triangular Schur form of A
        T     - upper triangular Schur form of B
        Z     - right Schur vectors
        evec  - (output) where to store eigenvectors
        w     - workspace

Return: success or error

Notes: 1) based on LAPACK routine DTGEVC
       2) eigenvectors are stored in the order that their
          eigenvalues appear in the Schur form
*/

static int
genv_get_right_eigenvectors(const gsl_matrix *S, const gsl_matrix *T,
                            gsl_matrix *Z,
                            gsl_matrix_complex *evec,
                            gsl_eigen_genv_workspace *w)
{
  const size_t N = w->size;
  const double small = GSL_DBL_MIN * N / GSL_DBL_EPSILON;
  const double big = 1.0 / small;
  const double bignum = 1.0 / (GSL_DBL_MIN * N);
  size_t i, j, k, end;
  int is;
  double anorm, bnorm;
  double temp, temp2, temp2r, temp2i;
  double ascale, bscale;
  double salfar, sbeta;
  double acoef, bcoefr, bcoefi, acoefa, bcoefa;
  double creala, cimaga, crealb, cimagb, cre2a, cim2a, cre2b, cim2b;
  double dmin, xmax;
  double scale;
  size_t nw, na;
  int lsa, lsb;
  int complex_pair;
  gsl_complex z_zero, z_one;
  double bdiag[2] = { 0.0, 0.0 };
  double sum[4];
  int il2by2;
  size_t jr, jc, ja;
  double xscale;
  gsl_vector_complex_view ecol;
  gsl_vector_view re, im, re2, im2;

  GSL_SET_COMPLEX(&z_zero, 0.0, 0.0);
  GSL_SET_COMPLEX(&z_one, 1.0, 0.0);

  /*
   * Compute the 1-norm of each column of (S, T) excluding elements
   * belonging to the diagonal blocks to check for possible overflow
   * in the triangular solver
   */

  anorm = fabs(gsl_matrix_get(S, 0, 0));
  if (N > 1)
    anorm += fabs(gsl_matrix_get(S, 1, 0));
  bnorm = fabs(gsl_matrix_get(T, 0, 0));

  gsl_vector_set(w->work1, 0, 0.0);
  gsl_vector_set(w->work2, 0, 0.0);

  for (j = 1; j < N; ++j)
    {
      temp = temp2 = 0.0;
      if (gsl_matrix_get(S, j, j - 1) == 0.0)
        end = j;
      else
        end = j - 1;

      for (i = 0; i < end; ++i)
        {
          temp += fabs(gsl_matrix_get(S, i, j));
          temp2 += fabs(gsl_matrix_get(T, i, j));
        }

      gsl_vector_set(w->work1, j, temp);
      gsl_vector_set(w->work2, j, temp2);

      for (i = end; i < GSL_MIN(j + 2, N); ++i)
        {
          temp += fabs(gsl_matrix_get(S, i, j));
          temp2 += fabs(gsl_matrix_get(T, i, j));
        }

      anorm = GSL_MAX(anorm, temp);
      bnorm = GSL_MAX(bnorm, temp2);
    }

  ascale = 1.0 / GSL_MAX(anorm, GSL_DBL_MIN);
  bscale = 1.0 / GSL_MAX(bnorm, GSL_DBL_MIN);

  complex_pair = 0;
  for (k = 0; k < N; ++k)
    {
      size_t je = N - 1 - k;

      if (complex_pair)
        {
          complex_pair = 0;
          continue;
        }

      nw = 1;
      if (je > 0)
        {
          if (gsl_matrix_get(S, je, je - 1) != 0.0)
            {
              complex_pair = 1;
              nw = 2;
            }
        }

      if (!complex_pair)
        {
          if (fabs(gsl_matrix_get(S, je, je)) <= GSL_DBL_MIN &&
              fabs(gsl_matrix_get(T, je, je)) <= GSL_DBL_MIN)
            {
              /* singular matrix pencil - unit eigenvector */
              for (i = 0; i < N; ++i)
                gsl_matrix_complex_set(evec, i, je, z_zero);

              gsl_matrix_complex_set(evec, je, je, z_one);

              continue;
            }

          /* clear vector */
          for (i = 0; i < N; ++i)
            gsl_vector_set(w->work3, i, 0.0);
        }
      else
        {
          /* clear vectors */
          for (i = 0; i < N; ++i)
            {
              gsl_vector_set(w->work3, i, 0.0);
              gsl_vector_set(w->work4, i, 0.0);
            }
        }

      if (!complex_pair)
        {
          /* real eigenvalue */

          temp = 1.0 / GSL_MAX(GSL_DBL_MIN,
                               GSL_MAX(fabs(gsl_matrix_get(S, je, je)) * ascale,
                                       fabs(gsl_matrix_get(T, je, je)) * bscale));
          salfar = (temp * gsl_matrix_get(S, je, je)) * ascale;
          sbeta = (temp * gsl_matrix_get(T, je, je)) * bscale;
          acoef = sbeta * ascale;
          bcoefr = salfar * bscale;
          bcoefi = 0.0;

          /* scale to avoid underflow */
          scale = 1.0;
          lsa = fabs(sbeta) >= GSL_DBL_MIN && fabs(acoef) < small;
          lsb = fabs(salfar) >= GSL_DBL_MIN && fabs(bcoefr) < small;
          if (lsa)
            scale = (small / fabs(sbeta)) * GSL_MIN(anorm, big);
          if (lsb)
            scale = GSL_MAX(scale, (small / fabs(salfar)) * GSL_MIN(bnorm, big));

          if (lsa || lsb)
            {
              scale = GSL_MIN(scale,
                        1.0 / (GSL_DBL_MIN *
                               GSL_MAX(1.0,
                                 GSL_MAX(fabs(acoef), fabs(bcoefr)))));
              if (lsa)
                acoef = ascale * (scale * sbeta);
              else
                acoef *= scale;

              if (lsb)
                bcoefr = bscale * (scale * salfar);
              else
                bcoefr *= scale;
            }

          acoefa = fabs(acoef);
          bcoefa = fabs(bcoefr);

          /* first component is 1 */
          gsl_vector_set(w->work3, je, 1.0);
          xmax = 1.0;

          /* compute contribution from column je of A and B to sum */

          for (i = 0; i < je; ++i)
            {
              gsl_vector_set(w->work3, i,
                bcoefr*gsl_matrix_get(T, i, je) -
                acoef * gsl_matrix_get(S, i, je));
            }
        }
      else
        {
          gsl_matrix_const_view vs =
            gsl_matrix_const_submatrix(S, je - 1, je - 1, 2, 2);
          gsl_matrix_const_view vt =
            gsl_matrix_const_submatrix(T, je - 1, je - 1, 2, 2);

          /* complex eigenvalue */

          gsl_schur_gen_eigvals(&vs.matrix,
                                &vt.matrix,
                                &bcoefr,
                                &temp2,
                                &bcoefi,
                                &acoef,
                                &temp);
          if (bcoefi == 0.0)
            {
              GSL_ERROR("gsl_schur_gen_eigvals failed on complex block", GSL_FAILURE);
            }

          /* scale to avoid over/underflow */
          acoefa = fabs(acoef);
          bcoefa = fabs(bcoefr) + fabs(bcoefi);
          scale = 1.0;

          if (acoefa*GSL_DBL_EPSILON < GSL_DBL_MIN && acoefa >= GSL_DBL_MIN)
            scale = (GSL_DBL_MIN / GSL_DBL_EPSILON) / acoefa;
          if (bcoefa*GSL_DBL_EPSILON < GSL_DBL_MIN && bcoefa >= GSL_DBL_MIN)
            scale = GSL_MAX(scale, (GSL_DBL_MIN/GSL_DBL_EPSILON) / bcoefa);
          if (GSL_DBL_MIN*acoefa > ascale)
            scale = ascale / (GSL_DBL_MIN * acoefa);
          if (GSL_DBL_MIN*bcoefa > bscale)
            scale = GSL_MIN(scale, bscale / (GSL_DBL_MIN*bcoefa));
          if (scale != 1.0)
            {
              acoef *= scale;
              acoefa = fabs(acoef);
              bcoefr *= scale;
              bcoefi *= scale;
              bcoefa = fabs(bcoefr) + fabs(bcoefi);
            }

          /* compute first two components of eigenvector */

          temp = acoef * gsl_matrix_get(S, je, je - 1);
          temp2r = acoef * gsl_matrix_get(S, je, je) -
                   bcoefr * gsl_matrix_get(T, je, je);
          temp2i = -bcoefi * gsl_matrix_get(T, je, je);

          if (fabs(temp) >= fabs(temp2r) + fabs(temp2i))
            {
              gsl_vector_set(w->work3, je, 1.0);
              gsl_vector_set(w->work4, je, 0.0);
              gsl_vector_set(w->work3, je - 1, -temp2r / temp);
              gsl_vector_set(w->work4, je - 1, -temp2i / temp);
            }
          else
            {
              gsl_vector_set(w->work3, je - 1, 1.0);
              gsl_vector_set(w->work4, je - 1, 0.0);
              temp = acoef * gsl_matrix_get(S, je - 1, je);
              gsl_vector_set(w->work3, je,
                (bcoefr*gsl_matrix_get(T, je - 1, je - 1) -
                 acoef*gsl_matrix_get(S, je - 1, je - 1)) / temp);
              gsl_vector_set(w->work4, je,
                bcoefi*gsl_matrix_get(T, je - 1, je - 1) / temp);
            }

          xmax = GSL_MAX(fabs(gsl_vector_get(w->work3, je)) +
                         fabs(gsl_vector_get(w->work4, je)),
                         fabs(gsl_vector_get(w->work3, je - 1)) +
                         fabs(gsl_vector_get(w->work4, je - 1)));

          /* compute contribution from column je and je - 1 */

          creala = acoef * gsl_vector_get(w->work3, je - 1);
          cimaga = acoef * gsl_vector_get(w->work4, je - 1);
          crealb = bcoefr * gsl_vector_get(w->work3, je - 1) -
                   bcoefi * gsl_vector_get(w->work4, je - 1);
          cimagb = bcoefi * gsl_vector_get(w->work3, je - 1) +
                   bcoefr * gsl_vector_get(w->work4, je - 1);
          cre2a = acoef * gsl_vector_get(w->work3, je);
          cim2a = acoef * gsl_vector_get(w->work4, je);
          cre2b = bcoefr * gsl_vector_get(w->work3, je) -
                  bcoefi * gsl_vector_get(w->work4, je);
          cim2b = bcoefi * gsl_vector_get(w->work3, je) +
                  bcoefr * gsl_vector_get(w->work4, je);

          for (i = 0; i < je - 1; ++i)
            {
              gsl_vector_set(w->work3, i,
                -creala * gsl_matrix_get(S, i, je - 1) +
                crealb * gsl_matrix_get(T, i, je - 1) -
                cre2a * gsl_matrix_get(S, i, je) +
                cre2b * gsl_matrix_get(T, i, je));
              gsl_vector_set(w->work4, i,
                -cimaga * gsl_matrix_get(S, i, je - 1) +
                cimagb * gsl_matrix_get(T, i, je - 1) -
                cim2a * gsl_matrix_get(S, i, je) +
                cim2b * gsl_matrix_get(T, i, je));
            }
        }

      dmin = GSL_MAX(GSL_DBL_MIN,
               GSL_MAX(GSL_DBL_EPSILON*acoefa*anorm,
                       GSL_DBL_EPSILON*bcoefa*bnorm));

      /* triangular solve of (a A - b B) x = 0 */

      il2by2 = 0;
      for (is = (int) je - (int) nw; is >= 0; --is)
        {
          j = (size_t) is;

          if (!il2by2 && j > 0)
            {
              if (gsl_matrix_get(S, j, j - 1) != 0.0)
                {
                  il2by2 = 1;
                  continue;
                }
            }

          bdiag[0] = gsl_matrix_get(T, j, j);
          if (il2by2)
            {
              na = 2;
              bdiag[1] = gsl_matrix_get(T, j + 1, j + 1);
            }
          else
            na = 1;

          if (nw == 1)
            {
              gsl_matrix_const_view sv =
                gsl_matrix_const_submatrix(S, j, j, na, na);
              gsl_vector_view xv, bv;

              bv = gsl_vector_subvector(w->work3, j, na);

              /*
               * the loop below expects the solution in the first column
               * of sum, so set stride to 2
               */
              xv = gsl_vector_view_array_with_stride(sum, 2, na);

              gsl_schur_solve_equation(acoef,
                                       &sv.matrix,
                                       bcoefr,
                                       bdiag[0],
                                       bdiag[1],
                                       &bv.vector,
                                       &xv.vector,
                                       &scale,
                                       &temp,
                                       dmin);
            }
          else
            {
              double bdat[4];
              gsl_matrix_const_view sv =
                gsl_matrix_const_submatrix(S, j, j, na, na);
              gsl_vector_complex_view xv =
                gsl_vector_complex_view_array(sum, na);
              gsl_vector_complex_view bv =
                gsl_vector_complex_view_array(bdat, na);
              gsl_complex z;

              bdat[0] = gsl_vector_get(w->work3, j);
              bdat[1] = gsl_vector_get(w->work4, j);
              if (na == 2)
                {
                  bdat[2] = gsl_vector_get(w->work3, j + 1);
                  bdat[3] = gsl_vector_get(w->work4, j + 1);
                }

              GSL_SET_COMPLEX(&z, bcoefr, bcoefi);

              gsl_schur_solve_equation_z(acoef,
                                         &sv.matrix,
                                         &z,
                                         bdiag[0],
                                         bdiag[1],
                                         &bv.vector,
                                         &xv.vector,
                                         &scale,
                                         &temp,
                                         dmin);
            }

          if (scale < 1.0)
            {
              for (jr = 0; jr <= je; ++jr)
                {
                  gsl_vector_set(w->work3, jr,
                    scale * gsl_vector_get(w->work3, jr));
                  if (nw == 2)
                    {
                      gsl_vector_set(w->work4, jr,
                        scale * gsl_vector_get(w->work4, jr));
                    }
                }
            }

          xmax = GSL_MAX(scale * xmax, temp);

          for (jr = 0; jr < na; ++jr)
            {
              gsl_vector_set(w->work3, j + jr, sum[jr*na]);
              if (nw == 2)
                gsl_vector_set(w->work4, j + jr, sum[jr*na + 1]);
            }

          if (j > 0)
            {
              xscale = 1.0 / GSL_MAX(1.0, xmax);
              temp = acoefa * gsl_vector_get(w->work1, j) +
                     bcoefa * gsl_vector_get(w->work2, j);
              if (il2by2)
                {
                  temp = GSL_MAX(temp,
                           acoefa * gsl_vector_get(w->work1, j + 1) +
                           bcoefa * gsl_vector_get(w->work2, j + 1));
                }

              temp = GSL_MAX(temp, GSL_MAX(acoefa, bcoefa));
              if (temp > bignum * xscale)
                {
                  for (jr = 0; jr <= je; ++jr)
                    {
                      gsl_vector_set(w->work3, jr,
                        xscale * gsl_vector_get(w->work3, jr));
                      if (nw == 2)
                        {
                          gsl_vector_set(w->work4, jr,
                            xscale * gsl_vector_get(w->work4, jr));
                        }
                    }
                  xmax *= xscale;
                }

              for (ja = 0; ja < na; ++ja)
                {
                  if (complex_pair)
                    {
                      creala = acoef * gsl_vector_get(w->work3, j + ja);
                      cimaga = acoef * gsl_vector_get(w->work4, j + ja);
                      crealb = bcoefr * gsl_vector_get(w->work3, j + ja) -
                               bcoefi * gsl_vector_get(w->work4, j + ja);
                      cimagb = bcoefi * gsl_vector_get(w->work3, j + ja) +
                               bcoefr * gsl_vector_get(w->work4, j + ja);
                      for (jr = 0; jr <= j - 1; ++jr)
                        {
                          gsl_vector_set(w->work3, jr,
                            gsl_vector_get(w->work3, jr) -
                            creala * gsl_matrix_get(S, jr, j + ja) +
                            crealb * gsl_matrix_get(T, jr, j + ja));
                          gsl_vector_set(w->work4, jr,
                            gsl_vector_get(w->work4, jr) -
                            cimaga * gsl_matrix_get(S, jr, j + ja) +
                            cimagb * gsl_matrix_get(T, jr, j + ja));
                        }
                    }
                  else
                    {
                      creala = acoef * gsl_vector_get(w->work3, j + ja);
                      crealb = bcoefr * gsl_vector_get(w->work3, j + ja);
                      for (jr = 0; jr <= j - 1; ++jr)
                        {
                          gsl_vector_set(w->work3, jr,
                            gsl_vector_get(w->work3, jr) -
                            creala * gsl_matrix_get(S, jr, j + ja) +
                            crealb * gsl_matrix_get(T, jr, j + ja));
                        }
                    } /* if (!complex_pair) */
                } /* for (ja = 0; ja < na; ++ja) */
            } /* if (j > 0) */

          il2by2 = 0;
        } /* for (i = 0; i < je - nw; ++i) */

      for (jr = 0; jr < N; ++jr)
        {
          gsl_vector_set(w->work5, jr,
            gsl_vector_get(w->work3, 0) * gsl_matrix_get(Z, jr, 0));
          if (nw == 2)
            {
              gsl_vector_set(w->work6, jr,
                gsl_vector_get(w->work4, 0) * gsl_matrix_get(Z, jr, 0));
            }
        }

      for (jc = 1; jc <= je; ++jc)
        {
          for (jr = 0; jr < N; ++jr)
            {
              gsl_vector_set(w->work5, jr,
                gsl_vector_get(w->work5, jr) +
                gsl_vector_get(w->work3, jc) * gsl_matrix_get(Z, jr, jc));
              if (nw == 2)
                {
                  gsl_vector_set(w->work6, jr,
                    gsl_vector_get(w->work6, jr) +
                    gsl_vector_get(w->work4, jc) * gsl_matrix_get(Z, jr, jc));
                }
            }
        }

      /* store the eigenvector */

      if (complex_pair)
        {
          ecol = gsl_matrix_complex_column(evec, je - 1);
          re = gsl_vector_complex_real(&ecol.vector);
          im = gsl_vector_complex_imag(&ecol.vector);

          ecol = gsl_matrix_complex_column(evec, je);
          re2 = gsl_vector_complex_real(&ecol.vector);
          im2 = gsl_vector_complex_imag(&ecol.vector);
        }
      else
        {
          ecol = gsl_matrix_complex_column(evec, je);
          re = gsl_vector_complex_real(&ecol.vector);
          im = gsl_vector_complex_imag(&ecol.vector);
        }

      for (jr = 0; jr < N; ++jr)
        {
          gsl_vector_set(&re.vector, jr, gsl_vector_get(w->work5, jr));
          if (complex_pair)
            {
              gsl_vector_set(&im.vector, jr, gsl_vector_get(w->work6, jr));
              gsl_vector_set(&re2.vector, jr, gsl_vector_get(w->work5, jr));
              gsl_vector_set(&im2.vector, jr, -gsl_vector_get(w->work6, jr));
            }
          else
            {
              gsl_vector_set(&im.vector, jr, 0.0);
            }
        }

      /* scale eigenvector */
      xmax = 0.0;
      if (complex_pair)
        {
          for (j = 0; j < N; ++j)
            {
              xmax = GSL_MAX(xmax,
                             fabs(gsl_vector_get(&re.vector, j)) +
                             fabs(gsl_vector_get(&im.vector, j)));
            }
        }
      else
        {
          for (j = 0; j < N; ++j)
            {
              xmax = GSL_MAX(xmax, fabs(gsl_vector_get(&re.vector, j)));
            }
        }

      if (xmax > GSL_DBL_MIN)
        {
          xscale = 1.0 / xmax;
          for (j = 0; j < N; ++j)
            {
              gsl_vector_set(&re.vector, j,
                             gsl_vector_get(&re.vector, j) * xscale);
              if (complex_pair)
                {
                  gsl_vector_set(&im.vector, j,
                                 gsl_vector_get(&im.vector, j) * xscale);
                  gsl_vector_set(&re2.vector, j,
                                 gsl_vector_get(&re2.vector, j) * xscale);
                  gsl_vector_set(&im2.vector, j,
                                 gsl_vector_get(&im2.vector, j) * xscale);
                }
            }
        }
    } /* for (k = 0; k < N; ++k) */

  return GSL_SUCCESS;
} /* genv_get_right_eigenvectors() */

/*
genv_normalize_eigenvectors()
  Normalize eigenvectors so that their Euclidean norm is 1

Inputs: alpha - eigenvalue numerators
        evec  - eigenvectors
*/

static void
genv_normalize_eigenvectors(gsl_vector_complex *alpha,
                            gsl_matrix_complex *evec)
{
  const size_t N = evec->size1;
  size_t i;     /* looping */
  gsl_complex ai;
  gsl_vector_complex_view vi;
  gsl_vector_view re, im;
  double scale; /* scaling factor */

  for (i = 0; i < N; ++i)
    {
      ai = gsl_vector_complex_get(alpha, i);
      vi = gsl_matrix_complex_column(evec, i);

      re = gsl_vector_complex_real(&vi.vector);

      if (GSL_IMAG(ai) == 0.0)
        {
          scale = 1.0 / gsl_blas_dnrm2(&re.vector);
          gsl_blas_dscal(scale, &re.vector);
        }
      else if (GSL_IMAG(ai) > 0.0)
        {
          im = gsl_vector_complex_imag(&vi.vector);

          scale = 1.0 / gsl_hypot(gsl_blas_dnrm2(&re.vector),
                                  gsl_blas_dnrm2(&im.vector));
          gsl_blas_zdscal(scale, &vi.vector);

          vi = gsl_matrix_complex_column(evec, i + 1);
          gsl_blas_zdscal(scale, &vi.vector);
        }
    }
} /* genv_normalize_eigenvectors() */


text:
#include <cblas.h>
#include <stdio.h>

/* compilar: gcc fonte.c -o exec -lblas */

int main( void ) {
  int lda = 3;
  float A[] = { 1, 0, -1,
                2, 3, 1 };

  int ldb = 2;
  float B[] = { 5, 1,
                6, -1,
                7, 2 };

  int ldc = 2;
  float C[] = { 0.00, 0.00,
                0.00, 0.00 };

  /* Compute C = A B */

  cblas_sgemm( CblasRowMajor, CblasNoTrans, CblasNoTrans, 2, 2, 3, 1.0, A, lda, B, ldb, 0.0, C, ldc );

  printf( "[ %g, %g\n", C[ 0 ], C[ 1 ] );
  printf( "  %g, %g ]\n", C[ 2 ], C[ 3 ] );

  return( 0 );
}

max_stars_repo_name: lellisls/PAD

max_issues_repo_name: lellisls/PAD

max_forks_repo_name: lellisls/PAD


text:
#include<stdio.h>
#include<math.h>
#include"recur.c"
#include <gsl/gsl_integration.h>
#include"gsl_integration.c"
#define NMAX 100
#include"timer.h"
#include"timer.c"
#include"adjust.c"

/* The integral is evaluated for n=0 to n=100 using both a derived recursion formula
 * and the GSL general integrator. The values of both are stored in their own
 * respective arrays. The values and their absolute error are printed in order to
 * compare the accuracy of the two methods. Then, both functions are timed to
 * determine which one is faster at generating the results. Finally, the timings of
 * the two functions are divided to provide a ratio describing the difference in
 * efficiency of the two methods.
*/
int main (void)
{
    double vals1[NMAX + 1], vals2[NMAX + 1], time1a, time1b, time2a, time2b, tmin = 1., tmax = 2.;      // Initializing all the "double" variables for use in the program//
    int nmin = 0, nmax = NMAX, count1 = 1000, count2 = 800;

    integral_recur (nmin, nmax, vals1); // Recursion formula stores I(0) to I(100) in the array vals1//
    integral_gen (nmin, nmax, vals2);   // GSL integration algorithm stores I(0) to I(100) in the array vals2//

    printf
        ("      --n--       --Recur--        --GSL Integrator--   --Absolute Error*10^9--\n");
    int k;

    for (k = 0; k < 101; k++)
    {
        printf ("   %6d     %15.15f     %15.15f      %15.15f \n", k, vals1[k],
            vals2[k], fabs (1.e9 * (vals1[k] - vals2[k])));
    }                           // Prints the values for each term in the sequence as generated by both of the algorithms //
    int p;

    printf
        ("\n--Timing for recurrence function--\nTime per call (usec) Total Time (sec)      Count\n");
    do
    {
        timer_start ();         // The timer is started.
        for (p = 0; p < count1; p++)    // The function is called "count1" times.
        {
            integral_recur (nmin, nmax, vals1); // The recursion formula stores each of the sequence values in the array vals1//
        }
        time1a = timer_stop (); // The timer is stopped.
        time1b = time1a / count1;       // The time per function call is calculated as the total time divided by the number of function calls.
        printf (" %10.2f usec     %10.6f sec    %10d\n", time1b * 1.e6,
            time1a, count1);
        count1 = adjust_rep_count (count1, time1a, tmin, tmax); // The number of function calls "count1" is adjusted.
    }
    while ((time1a > tmax) || (time1a < tmin)); // This is repeated until the total time is within the desired range.
    printf ("Time per recur function call = %f usec \n", time1b * 1.e6);
    printf
        ("\n--Timing for GSL integration function--\nTime per call (usec) Total Time (sec)      Count\n");
    do
    {
        timer_start ();         // The timer is started.
        for (p = 0; p < count2; p++)    // The function is called "count2" times.
        {
            integral_gen (nmin, nmax, vals2);   // The GSL integration formula stores each of the sequence values in the array vals1//
        }
        time2a = timer_stop (); // The timer is stopped.
        time2b = time2a / count2;       // The time per function call is calculated as the total time divided by the number of function calls.
        printf (" %10.2f usec     %10.6f sec    %10d\n", time2b * 1.e6,
            time2a, count2);
        count2 = adjust_rep_count (count2, time2a, tmin, tmax); // The number of function calls "count2" is adjusted.
    }
    while ((time2a > tmax) || (time2a < tmin)); // This is repeated until the total time is within the desired range.
    printf ("Time per GSL function call = %f usec \n", time2b * 1.e6);

    printf ("Ratio of times Tgsl / Trecur = %f\n", time2b / time1b);    // The ratio of times per function call of each function is printed//
    return 0;
}


text:
/**
 *
 * @generated d Tue Jan  7 11:45:26 2014
 *
 **/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <cblas.h>
#include <lapacke.h>
#include <plasma.h>
#include <core_blas.h>
#include "auxiliary.h"

/*-------------------------------------------------------------------
 * Check the orthogonality of Q
 */

int d_check_orthogonality(int M, int N, int LDQ, double *Q)
{
    double alpha, beta;
    double normQ;
    int info_ortho;
    int i;
    int minMN = min(M, N);
    double eps;
    double *work = (double *)malloc(minMN*sizeof(double));

    eps = LAPACKE_dlamch_work('e');
    alpha = 1.0;
    beta  = -1.0;

    /* Build the idendity matrix USE DLASET?*/
    double *Id = (double *) malloc(minMN*minMN*sizeof(double));
    memset((void*)Id, 0, minMN*minMN*sizeof(double));
    for (i = 0; i < minMN; i++)
        Id[i*minMN+i] = (double)1.0;

    /* Perform Id - Q'Q */
    if (M >= N)
        cblas_dsyrk(CblasColMajor, CblasUpper, CblasTrans, N, M, alpha, Q, LDQ, beta, Id, N);
    else
        cblas_dsyrk(CblasColMajor, CblasUpper, CblasNoTrans, M, N, alpha, Q, LDQ, beta, Id, M);

    normQ = LAPACKE_dlansy_work(LAPACK_COL_MAJOR, 'i', 'u', minMN, Id, minMN, work);

    printf("============\n");
    printf("Checking the orthogonality of Q \n");
    printf("||Id-Q'*Q||_oo / (N*eps) = %e \n",normQ/(minMN*eps));

    if ( isnan(normQ / (minMN * eps)) || (normQ / (minMN * eps) > 10.0) ) {
        printf("-- Orthogonality is suspicious ! \n");
        info_ortho=1;
    }
    else {
        printf("-- Orthogonality is CORRECT ! \n");
        info_ortho=0;
    }

    free(work); free(Id);

    return info_ortho;
}

/*------------------------------------------------------------
 *  Check the factorization QR
 */

int d_check_QRfactorization(int M, int N, double *A1, double *A2, int LDA, double *Q)
{
    double Anorm, Rnorm;
    double alpha, beta;
    int info_factorization;
    int i,j;
    double eps;

    eps = LAPACKE_dlamch_work('e');

    double *Ql       = (double *)malloc(M*N*sizeof(double));
    double *Residual = (double *)malloc(M*N*sizeof(double));
    double *work              = (double *)malloc(max(M,N)*sizeof(double));

    alpha=1.0;
    beta=0.0;

    if (M >= N) {
        /* Extract the R */
        double *R = (double *)malloc(N*N*sizeof(double));
        memset((void*)R, 0, N*N*sizeof(double));
        LAPACKE_dlacpy_work(LAPACK_COL_MAJOR,'u', M, N, A2, LDA, R, N);

        /* Perform Ql=Q*R */
        memset((void*)Ql, 0, M*N*sizeof(double));
        cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, M, N, N, (alpha), Q, LDA, R, N, (beta), Ql, M);
        free(R);
    }
    else {
        /* Extract the L */
        double *L = (double *)malloc(M*M*sizeof(double));
        memset((void*)L, 0, M*M*sizeof(double));
        LAPACKE_dlacpy_work(LAPACK_COL_MAJOR,'l', M, N, A2, LDA, L, M);

    /* Perform Ql=LQ */
        memset((void*)Ql, 0, M*N*sizeof(double));
        cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, M, N, M, (alpha), L, M, Q, LDA, (beta), Ql, M);
        free(L);
    }

    /* Compute the Residual */
    for (i = 0; i < M; i++)
        for (j = 0 ; j < N; j++)
            Residual[j*M+i] = A1[j*LDA+i]-Ql[j*M+i];

    Rnorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'i', M, N, Residual, M, work);
    Anorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'i', M, N, A2, LDA, work);

    if (M >= N) {
        printf("============\n");
        printf("Checking the QR Factorization \n");
        printf("-- ||A-QR||_oo/(||A||_oo.N.eps) = %e \n",Rnorm/(Anorm*N*eps));
    }
    else {
        printf("============\n");
        printf("Checking the LQ Factorization \n");
        printf("-- ||A-LQ||_oo/(||A||_oo.N.eps) = %e \n",Rnorm/(Anorm*N*eps));
    }

    if (isnan(Rnorm / (Anorm * N *eps)) || (Rnorm / (Anorm * N * eps) > 10.0) ) {
        printf("-- Factorization is suspicious ! \n");
        info_factorization = 1;
    }
    else {
        printf("-- Factorization is CORRECT ! \n");
        info_factorization = 0;
    }

    free(work); free(Ql); free(Residual);

    return info_factorization;
}

/*------------------------------------------------------------------------
 *  Check the factorization of the matrix A2
 */

int d_check_LLTfactorization(int N, double *A1, double *A2, int LDA, int uplo)
{
    double Anorm, Rnorm;
    double alpha;
    int info_factorization;
    int i,j;
    double eps;

    eps = LAPACKE_dlamch_work('e');

    double *Residual = (double *)malloc(N*N*sizeof(double));
    double *L1       = (double *)malloc(N*N*sizeof(double));
    double *L2       = (double *)malloc(N*N*sizeof(double));
    double *work              = (double *)malloc(N*sizeof(double));

    memset((void*)L1, 0, N*N*sizeof(double));
    memset((void*)L2, 0, N*N*sizeof(double));

    alpha= 1.0;

    LAPACKE_dlacpy_work(LAPACK_COL_MAJOR,' ', N, N, A1, LDA, Residual, N);

    /* Dealing with L'L or U'U  */
    if (uplo == PlasmaUpper){
        LAPACKE_dlacpy_work(LAPACK_COL_MAJOR,'u', N, N, A2, LDA, L1, N);
        LAPACKE_dlacpy_work(LAPACK_COL_MAJOR,'u', N, N, A2, LDA, L2, N);
        cblas_dtrmm(CblasColMajor, CblasLeft, CblasUpper, CblasTrans, CblasNonUnit, N, N, (alpha), L1, N, L2, N);
    }
    else{
        LAPACKE_dlacpy_work(LAPACK_COL_MAJOR,'l', N, N, A2, LDA, L1, N);
        LAPACKE_dlacpy_work(LAPACK_COL_MAJOR,'l', N, N, A2, LDA, L2, N);
        cblas_dtrmm(CblasColMajor, CblasRight, CblasLower, CblasTrans, CblasNonUnit, N, N, (alpha), L1, N, L2, N);
    }

    /* Compute the Residual || A -L'L|| */
    for (i = 0; i < N; i++)
        for (j = 0; j < N; j++)
           Residual[j*N+i] = L2[j*N+i] - Residual[j*N+i];

    Rnorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'i', N, N, Residual, N, work);
    Anorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'i', N, N, A1, LDA, work);

    printf("============\n");
    printf("Checking the Cholesky Factorization \n");
    printf("-- ||L'L-A||_oo/(||A||_oo.N.eps) = %e \n",Rnorm/(Anorm*N*eps));

    if ( isnan(Rnorm/(Anorm*N*eps)) || (Rnorm/(Anorm*N*eps) > 10.0) ){
        printf("-- Factorization is suspicious ! \n");
        info_factorization = 1;
    }
    else{
        printf("-- Factorization is CORRECT ! \n");
        info_factorization = 0;
    }

    free(Residual); free(L1); free(L2); free(work);

    return info_factorization;
}

/*--------------------------------------------------------------
 * Check the gemm
 */
double d_check_gemm(PLASMA_enum transA, PLASMA_enum transB, int M, int N, int K,
                   double alpha, double *A, int LDA,
                   double *B, int LDB,
                   double beta, double *Cplasma,
                   double *Cref, int LDC,
                   double *Cinitnorm, double *Cplasmanorm, double *Clapacknorm )
{
    double beta_const = -1.0;
    double Rnorm;
    double *work = (double *)malloc(max(K,max(M, N))* sizeof(double));

    *Cinitnorm   = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'i', M, N, Cref,    LDC, work);
    *Cplasmanorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'i', M, N, Cplasma, LDC, work);

    cblas_dgemm(CblasColMajor, (CBLAS_TRANSPOSE)transA, (CBLAS_TRANSPOSE)transB, M, N, K,
                (alpha), A, LDA, B, LDB, (beta), Cref, LDC);

    *Clapacknorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'i', M, N, Cref, LDC, work);

    cblas_daxpy(LDC * N, (beta_const), Cplasma, 1, Cref, 1);

    Rnorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'i', M, N, Cref, LDC, work);

    free(work);

    return Rnorm;
}

/*--------------------------------------------------------------
 * Check the trsm
 */
double d_check_trsm(PLASMA_enum side, PLASMA_enum uplo, PLASMA_enum trans, PLASMA_enum diag,
                   int M, int NRHS, double alpha,
                   double *A, int LDA,
                   double *Bplasma, double *Bref, int LDB,
                   double *Binitnorm, double *Bplasmanorm, double *Blapacknorm )
{
    double beta_const = -1.0;
    double Rnorm;
    double *work = (double *)malloc(max(M, NRHS)* sizeof(double));
    /*double eps = LAPACKE_dlamch_work('e');*/

    *Binitnorm   = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'i', M, NRHS, Bref,    LDB, work);
    *Bplasmanorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'm', M, NRHS, Bplasma, LDB, work);

    cblas_dtrsm(CblasColMajor, (CBLAS_SIDE)side, (CBLAS_UPLO)uplo,
                (CBLAS_TRANSPOSE)trans, (CBLAS_DIAG)diag, M, NRHS,
                (alpha), A, LDA, Bref, LDB);

    *Blapacknorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'm', M, NRHS, Bref, LDB, work);

    cblas_daxpy(LDB * NRHS, (beta_const), Bplasma, 1, Bref, 1);

    Rnorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'm', M, NRHS, Bref, LDB, work);
    Rnorm = Rnorm / *Blapacknorm;
    /* max(M,NRHS) * eps);*/

    free(work);

    return Rnorm;
}

/*--------------------------------------------------------------
 * Check the solution
 */

double d_check_solution(int M, int N, int NRHS, double *A, int LDA,
                      double *B,  double *X, int LDB,
                      double *anorm, double *bnorm, double *xnorm )
{
/*     int info_solution; */
    double Rnorm = -1.00;
    double zone  =  1.0;
    double mzone = -1.0;
    double *work = (double *)malloc(max(M, N)* sizeof(double));

    *anorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'i', M, N,    A, LDA, work);
    *xnorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'i', M, NRHS, X, LDB, work);
    *bnorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'i', N, NRHS, B, LDB, work);

    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, M, NRHS, N, (zone), A, LDA, X, LDB, (mzone), B, LDB);

    Rnorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'i', N, NRHS, B, LDB, work);

    free(work);

    return Rnorm;
}


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

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

#include <time.h>

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

#include <gsl/gsl_errno.h>

#include "complex_internal.h"
#include "bitreverse.h"
#include "urand.c"

void my_error_handler (const char *reason, const char *file,
                       int line, int err);

int
main (int argc, char *argv[])
{
  double *data, *fft_data;
  gsl_fft_complex_wavetable * cw;
  unsigned int i, logn;
  int result;
  int status;
  clock_t start, end;
  int resolution = CLOCKS_PER_SEC;
  unsigned int n = 1;

  gsl_set_error_handler (&my_error_handler);

  if (argc == 2)
    {
      n = strtol (argv[1], NULL, 0);
    }
  else
    {
      printf ("test n\n");
      exit (EXIT_FAILURE);
    }

  cw = gsl_fft_complex_wavetable_alloc (n);
  cwork = gsl_fft_complex_wavetable_workspace_alloc (n);
  status = gsl_fft_complex_init (n, cw, cwork);
  status = gsl_fft_complex_generate (n, cw, cwork);

  data = (double *) malloc (n * 2 * sizeof (double));
  fft_data = (double *) malloc (n * 2 * sizeof (double));

  for (i = 0; i < n; i++)
    {
      REAL(data,1,i) = urand ();
      IMAG(data,1,i) = urand ();
    }

  /* compute the fft */

  memcpy (fft_data, data, n * 2 * sizeof (double));

  start = clock ();
  i = 0;
  do
    {
      status = gsl_fft_complex_forward (fft_data, 1, n, cw, cwork);
      i++;
      end = clock ();
    }
  while (end < start + resolution && status == 0);

  if (status == 0)
    {
      printf ("n = %d gsl_fft_complex_forward %f seconds\n", n, (end - start) / ((double) i) / ((double) CLOCKS_PER_SEC));
    }
  else
    {
      printf ("MR fft failed\n");
    }

  /* compute the fft with radix2 */
  memcpy (fft_data, data, n * 2 * sizeof (double));

  result = fft_binary_logn(n) ;

  if (result == -1) {
    exit(EXIT_SUCCESS) ;
  } else {
    logn = result ;
  }

  start = clock ();
  i = 0;
  do
    {
      status = gsl_fft_complex_radix2_forward (fft_data, 1, n);
      i++;
      end = clock ();
    }
  while (end < start + resolution && status == 0);

  if (status == 0)
    {
      printf ("n = %d gsl_fft_complex_radix2_forward %f seconds\n", n, (end - start) / ((double) i) / ((double) CLOCKS_PER_SEC));
    }
  else
    {
      printf ("fft_radix2: not a power of 2\n");
    }

  start = clock ();
  i = 0;
  do
    {
      status = fft_complex_bitreverse_order (fft_data, 1, n, logn);
      i++;
      end = clock ();
    }
  while (end < start + resolution && status == 0);

  printf ("n = %d gsl_fft_complex_bitreverse_order %f seconds\n", n, (end - start) / ((double) i) / ((double) CLOCKS_PER_SEC));

  return 0;
}

void
my_error_handler (const char *reason, const char *file, int line, int err)
{
  printf ("error: %s in %s at %d (gsl_errno=%d)\n", reason, file, line, err);
}

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
#ifndef _INC_REGRESSION_
#define _INC_REGRESSION_
/**
 * contains the functions needed to handle the linear regression
 * this entails, creating hvectors (for given x vectors in the parameter space)
 * estimatingBeta, this needs the inverse covariance matrix and the training data
 * creating hmatrix, a done once kind of job, slap together hvectors for each of the design points
 */

#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_sf.h>

/**
 * \brief creates a vector of regression coefficients at x_location
 *
 * create a h_vector evaluated at given x_location.,
 * this currently only creates a set of linear basis fns, you could add another loop which
 * would add squares etc
 *
 * on calling h_vector should be initliased to a vector of length nregression_fns which
 * in the default case nregression_fns = nparams + 1
 *
 * the value of options->nregression_fns is going to be very important for this
 *
 * this pointer is set to one of the fns in this file by optstruct.c:setup_regression
 */

void (*makeHVector)(gsl_vector *h_vector, gsl_vector *x_location, int nparams);

//void makeHVector(gsl_vector *h_vector, gsl_vector *x_location, int nparams);

void makeHMatrix(gsl_matrix *h_matrix, gsl_matrix *xmodel, int nmodel_points, int nparams, int nregresion_fns);

void makeHMatrix_fnptr(gsl_matrix *h_matrix, gsl_matrix *xmodel, int nmodel_points, int nparams, int nregression_fns,
											 void (*makeHVector_ptr)(gsl_vector *h_vector, gsl_vector *x_location, int nparams));

void estimateBeta(gsl_vector *beta_vector, gsl_matrix *h_matrix, gsl_matrix* cinverse, gsl_vector *trainingvector, int nmodel_points, int nregression_fns);

void makeHVector_linear(	gsl_vector *h_vector, gsl_vector *x_location, int nparams);
void makeHVector_quadratic( gsl_vector *h_vector, gsl_vector *x_location, int nparams);
void makeHVector_trivial(gsl_vector *h_vector, gsl_vector *x_location, int nparams);
void makeHVector_cubic( gsl_vector *h_vector, gsl_vector *x_location, int nparams);
#endif

max_stars_repo_name: jackdawjackdaw/emulator

max_issues_repo_name: MADAI/MADAIEmulator

max_forks_repo_name: MADAI/MADAIEmulator