text:
#include <stdio.h>
#include <stdlib.h>
#include <cblas.h>
#include <math.h>
#include "../inc/knnring.h"

// Application Entry Point
knnresult kNN(double * X, double * Y, int n, int m, int d, int k)
{

	// Calculate distances matrix D - D is row-major and nxm
	double* tempD = calculateD(X, Y, n, m, d, k);

	double* D = calloc(m*n, sizeof(double));

	// Transpose D to mxn
	//cblas_dimatcopy(CblasRowMajor, CblasTrans, n, m, 1.0, D, m, n);

	double* identityMat = calloc(n*n, sizeof(double));
	for(int f=0; f<n; f++)
		identityMat[f*n + f] = 1;

	cblas_dgemm(CblasRowMajor, CblasTrans, CblasNoTrans, m, n, n, 1, tempD, m, identityMat, n, 0, D, n);

	// Free memory
	free(identityMat);
	free(tempD);

	// Create results struct
	knnresult results;
	results.k = k;
	results.m = m;
	results.ndist = calloc(m*k, sizeof(double));
	results.nidx = calloc(m*k, sizeof(int));

	// Create ids array for the X array
	int* ids = calloc(n, sizeof(int));

	// K-Select using quickselect to find k smallest distances for each point of Y
	for(int j = 0; j < m; j++){

		// Re-set ids vector before executing quickselect each time
		for(int i = 0; i < n; i++)
				ids[i] = i;

		// Call quickselect on the row of D to sort it up to its k-th element
		kselect(D + j * n, ids, n, k);

		// Quicksort the results
		quickSort(D + j * n, ids, 0, n-1);

		// Write results (point id and distance)
		for(int l = 0; l < k; l++){

			results.ndist[j * k + l] = D[j * n + l];
			results.nidx[j * k + l] = ids[l];
		}
	}

	// Free memory and return results
	free(ids);
	return results;

}

// Function used to calculate D = sqrt(sum(X.^2,2) -2* X*Y.' + sum(Y.^2,2).');
double* calculateD(double * X, double * Y, int n, int m, int d, int k){

		// Temporary variable
		double temp = 0;

		// Distance matrix for results
		double* D = calloc(n*m, sizeof(double));

		// Euclidean norm vectors
		double* normX = calloc(n, sizeof(double));
		double* normY = calloc(m, sizeof(double));

		// Matrice to store -2*X*Y'
		double* XY = calloc(n*m, sizeof(double));

		// Calculate -2*XY (https://software.intel.com/en-us/mkl-tutorial-c-multiplying-matrices-using-dgemm)
		cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans, n, m, d, -2, X, d, Y, d, 0, XY, m);

		// Calculate sum(X.^2,2), sum(Y.^2,2)
		for(int i = 0; i < n; i++){
			temp = cblas_dnrm2(d, X+i*d, 1);
			normX[i] = temp * temp;
		}
		for(int i = 0; i < m; i++){
			temp = cblas_dnrm2(d, Y+i*d, 1);
			normY[i] = temp * temp;
		}

		// XY = sum(X.^2,2) -2* X*Y.'
	        for (int i=0; i<n; i++)
	    		for(int j=0; j<m; j++)
				XY[i*m+j] += normX[i] + normY[j];

		// D = sqrt(sum(X.^2,2) -2* X*Y.' + sum(Y.^2,2).');
		for(int i = 0; i < n*m; i++)
			D[i] = sqrt(fabs(XY[i]));

		// Free memory
		free(normX);
		free(normY);
		free(XY);

		return D;
}

// A utility function to swap two elements
void swap_double(double *a, double *b){
    double t = *a;
    *a = *b;
    *b = t;
}

// A utility function to swap two elements
void swap_int(int *a, int *b){
    int t = *a;
    *a = *b;
    *b = t;
}

// QuickSort Partition function. low and high are the range of indexes in arr where partition should work
int partition (double arr[], int *ids, int low, int high){

    // Select a pivot and initialize flag to position of smallest element before pivot
    double pivot = arr[high];
    int i = (low - 1);

    // Go through the array examining each element
    for (int j = low; j <= high - 1; j++)
    {
        // If current element is smaller than the pivot, increment i and swap it out with the one currently pointed by i
        if (arr[j] < pivot)
        {
            i++;
            // Swap distances and corresponding point ids
            swap_double(&arr[i], &arr[j]);
            swap_int(&ids[i], &ids[j]);
        }
    }

    // Finally place pivot in its correct position in the array and return the position as the middle point
    swap_double(&arr[i + 1], &arr[high]);
    swap_int(&ids[i + 1], &ids[high]);
    return (i + 1);
}

// Returns the median using the QuickSelect algorithm
double kselect(double arr[], int *ids, int length, int k){
	return quickSelect(arr, ids, length, k);
}

// Returns the idx-th element of arr when arr is sorted. idx is the index (starting from 1) of the point we want to find when the array is sorted.
double quickSelect(double arr[], int *ids, int length, int idx){

    // Check to end recursion
    if (length == 1){
        return arr[0];
    }

    // Select last array element as pivot
    double pivot = arr[length - 1];
    // Get index of pivot after we partition the array
    int pivotIndex = partition(arr, ids, 0, length - 1);

    // Create the higher and lower arrays that occur after partitioning in QuickSort fashion
    int lowerLength = pivotIndex;
    pivotIndex++;
    int higherLength = (length - (lowerLength + 1));
    // At this point pivotIndex, lowerLength and higherLength all start from 1 not 0

    // Variable to store result of following recursive calls
    double result = 0;

    // This means that the point we're looking (median in our case) is in the lower partition
    if (idx <= lowerLength){
        result = quickSelect(arr, ids, lowerLength, idx);
    }
    // This means that idx-th element is our pivot point
    else if(idx == pivotIndex){
        result = pivot;
    }
    // This means that the median is in the higher partition
    else{
        result = quickSelect(arr + pivotIndex, ids + pivotIndex, higherLength, idx - pivotIndex);
    }

    // Return result
    return result;
}

// Simple quicksort implementation that also sorts the ids array according to the way the main array was sorted
void quickSort(double arr[], int *ids, int low, int high)
{
    if (low < high)
    {
        // idx is the partition point
        int idx = partition(arr, ids, low, high);

        // Divide and conquer
        quickSort(arr, ids, low, idx - 1);
        quickSort(arr, ids, idx + 1, high);
    }
}


text:
#ifndef __GSL_PERMUTE_MATRIX_H__
#define __GSL_PERMUTE_MATRIX_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <gsl/gsl_permute_matrix_complex_long_double.h>
#include <gsl/gsl_permute_matrix_complex_double.h>
#include <gsl/gsl_permute_matrix_complex_float.h>

#include <gsl/gsl_permute_matrix_long_double.h>
#include <gsl/gsl_permute_matrix_double.h>
#include <gsl/gsl_permute_matrix_float.h>

#include <gsl/gsl_permute_matrix_ulong.h>
#include <gsl/gsl_permute_matrix_long.h>

#include <gsl/gsl_permute_matrix_uint.h>
#include <gsl/gsl_permute_matrix_int.h>

#include <gsl/gsl_permute_matrix_ushort.h>
#include <gsl/gsl_permute_matrix_short.h>

#include <gsl/gsl_permute_matrix_uchar.h>
#include <gsl/gsl_permute_matrix_char.h>

#endif /* __GSL_PERMUTE_MATRIX_H__ */

ext: h

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:
/*insert GNU extensions*/
#define _GNU_SOURCE
/*in particular, use of NAN extension*/

/* used by <errno.h> */
extern int errno;

/* GSL ----------- */
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>
/* --------------- */

#define GET(x,i) gsl_vector_get(x,i)
#define SET(x,i,y) gsl_vector_set(x,i,y)

struct multimin_params {
  double step_size;
  double tol;
  unsigned maxiter;
  double epsabs;
  double maxsize;
  unsigned method;
  unsigned verbosity;
};

void
multimin(size_t,double *,double *,
	 const unsigned *,const double *,const double *,
	 void (*) (const size_t,const double *,void *,double *),
	 void (*) (const size_t,const double *, void *,double *),
	 void (*) (const size_t,const double *, void *,double *,double *),
	 void *,
	 const struct multimin_params);

max_stars_repo_name: dortegadelv/IMaDNA

max_issues_repo_name: dortegadelv/IMaDNA


text:
#ifndef SOURCE_LINGAUSS_H_
#define SOURCE_LINGAUSS_H_

// *** source/lin_gauss.h ***
// Author: Kevin Wolz, date: 06/2018
//
// This program defines the class "LinGauss" describing linear Gaussian models.

#include <iostream>
#include <fstream>
#include <string>
#include <sstream>
#include <vector>
#include <stdlib.h>
#include <sys/stat.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_integration.h>
#include "utils.h"

class LinGauss {
 public:
	unsigned int datdim_;
	unsigned int pardim_;
	// Ensures that g_matrix_, covar_matrix_, prior_matrix_, prior_param_values_,
	// data_mean_values_ are initialized and != 0.
	bool initialized_;
	// Ensures that cover_matrix_ is inverted and stored in inv_cover_matrix_.
	bool inverted_;
	// Ensures that this instance of LinGauss is initialized and fisher_matrix_
	// is calculated.
	bool fisherized_;
	gsl_matrix* g_matrix_; 						// design matrix of the model
	gsl_matrix* covar_matrix_;				// data covariance matrix
	gsl_matrix* inv_covar_matrix_; 		// inverse data covariance matrix
	gsl_matrix* fisher_matrix_; 			// Fisher matrix
	gsl_matrix* prior_matrix_; 			 	// prior precision matrix
	gsl_vector* prior_param_values_; 	// prior fiducial parameter values
	gsl_vector* param_values_; 				// ML estimator values
	gsl_vector* data_mean_values_; 		// mean of the data vector

	LinGauss(unsigned int num_data, unsigned int num_params);
	void DisallocateMembers();
	void CalculateFisherMatrix();
	void CalculateMlEstimator(gsl_vector* data);
	double CalculateChisquared(gsl_vector* data, gsl_vector* params);
};
#endif // SOURCE_LINGAUSS_H_

max_stars_repo_name: kevguitar/Robustness

max_issues_repo_name: kevguitar/Robustness

max_forks_repo_name: kevguitar/Robustness


text:
#pragma once

#include <renderer/Shading.h>
#include <renderer/VertexSpecification.h>
#include <gsl/span>

#include <cassert>

namespace ad {
namespace graphics {

struct GenericDrawer
{
    VertexArrayObject mVertexArray;
    Program mProgram;

    std::size_t mVertexCount{0};
    std::size_t mInstanceCount{1};

    template <class T_vertex>
    VertexBufferObject addVertexBuffer(
        const std::initializer_list<AttributeDescription> & aAttributes,
        gsl::span<T_vertex> aVertices)
    {
        if (mVertexCount)
        {
            assert(static_cast<std::size_t>(aVertices.size()) == mVertexCount);
        }
        else
        {
            mVertexCount = aVertices.size();
        }

        return loadVertexBuffer(mVertexArray,
                                aAttributes,
                                sizeof(T_vertex),
                                sizeof(T_vertex)*aVertices.size(),
                                aVertices.data());
    }

    void render() const
    {
        glBindVertexArray(mVertexArray);
        glUseProgram(mProgram);

        glDrawArraysInstanced(GL_TRIANGLE_STRIP,
                              0,
                              mVertexCount,
                              mInstanceCount);
    }
};

} // namespace graphics
} // namespace ad

max_stars_repo_name: Adnn/Graphics

max_issues_repo_name: Adnn/Graphics

max_forks_repo_name: ShredEagle/graphics


text:
/**
 * @file coroutine/linux.h
 * @author github.com/luncliff (luncliff@gmail.com)
 * @copyright CC BY 4.0
 */
#ifndef COROUTINE_SYSTEM_WRAPPER_H
#define COROUTINE_SYSTEM_WRAPPER_H
#if !(defined(__linux__))
#error "expect Linux platform for this file"
#endif
#include <sys/epoll.h> // for Linux epoll

#include <coroutine/return.h>
#include <gsl/gsl>

/**
 * @defgroup Linux
 */

namespace coro {

/**
 * @brief RAII wrapping for epoll file descriptor
 * @ingroup Linux
 */
class epoll_owner final {
    int64_t epfd;

  public:
    /**
     * @brief create a fd with `epoll`. Throw if the function fails.
     * @see kqeueue
     * @throw system_error
     */
    epoll_owner() noexcept(false);
    /**
     * @brief close the current epoll file descriptor
     */
    ~epoll_owner() noexcept;
    epoll_owner(const epoll_owner&) = delete;
    epoll_owner(epoll_owner&&) = delete;
    epoll_owner& operator=(const epoll_owner&) = delete;
    epoll_owner& operator=(epoll_owner&&) = delete;

  public:
    /**
     * @brief bind the fd to epoll
     * @param fd
     * @param req
     * @see epoll_ctl
     * @throw system_error
     */
    void try_add(uint64_t fd, epoll_event& req) noexcept(false);

    /**
     * @brief unbind the fd to epoll
     * @param fd
     * @see epoll_ctl
     */
    void remove(uint64_t fd);

    /**
     * @brief fetch all events for the given kqeueue descriptor
     * @param wait_ms millisecond to wait
     * @param list
     * @return ptrdiff_t
     * @see epoll_wait
     * @throw system_error
     *
     * Timeout is not an error for this function
     */
    ptrdiff_t wait(uint32_t wait_ms,
                   gsl::span<epoll_event> list) noexcept(false);

  public:
    /**
     * @brief return temporary awaitable object for given event
     * @param req input for `change` operation
     * @see change
     *
     * There is no guarantee of reusage of returned awaiter object
     * When it is awaited, and `req.udata` is null(0),
     * the value is set to `coroutine_handle<void>`
     *
     * ```cpp
     * auto edge_in_async(epoll_owner& ep, int64_t fd) -> frame_t {
     *     epoll_event req{};
     *     req.events = EPOLLET | EPOLLIN | EPOLLONESHOT;
     *     req.data.ptr = nullptr;
     *     co_await ep.submit(fd, req);
     * }
     * ```
     */
    [[nodiscard]] auto submit(int64_t fd, epoll_event& req) noexcept {
        class awaiter final : public suspend_always {
            epoll_owner& ep;
            int64_t fd;
            epoll_event& req;

          public:
            constexpr awaiter(epoll_owner& _ep, int64_t _fd, epoll_event& _req)
                : ep{_ep}, fd{_fd}, req{_req} {
            }

          public:
            void await_suspend(coroutine_handle<void> coro) noexcept(false) {
                if (req.data.ptr == nullptr)
                    req.data.ptr = coro.address();
                return ep.try_add(fd, req);
            }
        };
        return awaiter{*this, fd, req};
    }
};

/**
 * @brief RAII + stateful `eventfd`
 * @see https://github.com/grpc/grpc/blob/master/src/core/lib/iomgr/is_epollexclusive_available.cc
 * @ingroup Linux
 *
 * If the object is signaled(`set`),
 * the bound `epoll_owner` will yield suspended coroutine through `epoll_event`'s user data.
 *
 * Its object can be `co_await`ed multiple times
 */
class event final {
    uint64_t state;

  public:
    event() noexcept(false);
    ~event() noexcept;
    event(const event&) = delete;
    event(event&&) = delete;
    event& operator=(const event&) = delete;
    event& operator=(event&&) = delete;

    uint64_t fd() const noexcept;
    bool is_set() const noexcept;
    void set() noexcept(false);
    void reset() noexcept(false);
};

/**
 * @brief Bind the given `event`(`eventfd`) to `epoll_owner`(Epoll)
 *
 * @param ep  epoll_owner
 * @param efd event
 * @see event
 * @return awaitable struct for the binding
 * @ingroup Linux
 */
auto wait_in(epoll_owner& ep, event& efd) {
    class awaiter : epoll_event {
        epoll_owner& ep;
        event& efd;

      public:
        /**
         * @brief Prepares one-time registration
         */
        awaiter(epoll_owner& _ep, event& _efd) noexcept
            : epoll_event{}, ep{_ep}, efd{_efd} {
            this->events = EPOLLET | EPOLLIN | EPOLLONESHOT;
        }

        bool await_ready() const noexcept {
            return efd.is_set();
        }
        /**
         * @brief Wait for `write` to given `eventfd`
         */
        void await_suspend(coroutine_handle<void> coro) noexcept(false) {
            this->data.ptr = coro.address();
            return ep.try_add(efd.fd(), *this);
        }
        /**
         * @brief Reset the given event object when resumed
         */
        void await_resume() noexcept {
            return efd.reset();
        }
    };
    return awaiter{ep, efd};
}

} // namespace coro

#endif // COROUTINE_SYSTEM_WRAPPER_H

max_stars_repo_name: lanza/coroutine

max_issues_repo_name: lanza/coroutine

max_forks_repo_name: lanza/coroutine


text:
/////////////////////////////////////////////////////////////////////
// = NMatrix
//
// A linear algebra library for scientific computation in Ruby.
// NMatrix is part of SciRuby.
//
// NMatrix was originally inspired by and derived from NArray, by
// Masahiro Tanaka: http://narray.rubyforge.org
//
// == Copyright Information
//
// SciRuby is Copyright (c) 2010 - 2014, Ruby Science Foundation
// NMatrix is Copyright (c) 2012 - 2014, John Woods and the Ruby Science Foundation
//
// Please see LICENSE.txt for additional copyright notices.
//
// == Contributing
//
// By contributing source code to SciRuby, you agree to be bound by
// our Contributor Agreement:
//
// * https://github.com/SciRuby/sciruby/wiki/Contributor-Agreement
//
// == gemv.h
//
// Header file for interface with ATLAS's CBLAS gemv functions and
// native templated version of LAPACK's gemv function.
//

#ifndef GEMV_H
# define GEMV_H

extern "C" { // These need to be in an extern "C" block or you'll get all kinds of undefined symbol errors.
#if defined HAVE_CBLAS_H
  #include <cblas.h>
#elif defined HAVE_ATLAS_CBLAS_H
  #include <atlas/cblas.h>
#endif
}

namespace nm { namespace math {

/*
 * GEneral Matrix-Vector multiplication: based on dgemv.f from Netlib.
 *
 * This is an extremely inefficient algorithm. Recommend using ATLAS' version instead.
 *
 * Template parameters: LT -- long version of type T. Type T is the matrix dtype.
 */
template <typename DType>
inline bool gemv(const enum CBLAS_TRANSPOSE Trans, const int M, const int N, const DType* alpha, const DType* A, const int lda,
          const DType* X, const int incX, const DType* beta, DType* Y, const int incY) {
  int lenX, lenY, i, j;
  int kx, ky, iy, jx, jy, ix;

  typename LongDType<DType>::type temp;

  // Test the input parameters
  if (Trans < 111 || Trans > 113) {
    rb_raise(rb_eArgError, "GEMV: TransA must be CblasNoTrans, CblasTrans, or CblasConjTrans");
    return false;
  } else if (lda < std::max(1, N)) {
    fprintf(stderr, "GEMV: N = %d; got lda=%d", N, lda);
    rb_raise(rb_eArgError, "GEMV: Expected lda >= max(1, N)");
    return false;
  } else if (incX == 0) {
    rb_raise(rb_eArgError, "GEMV: Expected incX != 0\n");
    return false;
  } else if (incY == 0) {
    rb_raise(rb_eArgError, "GEMV: Expected incY != 0\n");
    return false;
  }

  // Quick return if possible
  if (!M or !N or (*alpha == 0 and *beta == 1)) return true;

  if (Trans == CblasNoTrans) {
    lenX = N;
    lenY = M;
  } else {
    lenX = M;
    lenY = N;
  }

  if (incX > 0) kx = 0;
  else          kx = (lenX - 1) * -incX;

  if (incY > 0) ky = 0;
  else          ky =  (lenY - 1) * -incY;

  // Start the operations. In this version, the elements of A are accessed sequentially with one pass through A.
  if (*beta != 1) {
    if (incY == 1) {
      if (*beta == 0) {
        for (i = 0; i < lenY; ++i) {
          Y[i] = 0;
        }
      } else {
        for (i = 0; i < lenY; ++i) {
          Y[i] *= *beta;
        }
      }
    } else {
      iy = ky;
      if (*beta == 0) {
        for (i = 0; i < lenY; ++i) {
          Y[iy] = 0;
          iy += incY;
        }
      } else {
        for (i = 0; i < lenY; ++i) {
          Y[iy] *= *beta;
          iy += incY;
        }
      }
    }
  }

  if (*alpha == 0) return false;

  if (Trans == CblasNoTrans) {

    // Form  y := alpha*A*x + y.
    jx = kx;
    if (incY == 1) {
      for (j = 0; j < N; ++j) {
        if (X[jx] != 0) {
          temp = *alpha * X[jx];
          for (i = 0; i < M; ++i) {
            Y[i] += A[j+i*lda] * temp;
          }
        }
        jx += incX;
      }
    } else {
      for (j = 0; j < N; ++j) {
        if (X[jx] != 0) {
          temp = *alpha * X[jx];
          iy = ky;
          for (i = 0; i < M; ++i) {
            Y[iy] += A[j+i*lda] * temp;
            iy += incY;
          }
        }
        jx += incX;
      }
    }

  } else { // TODO: Check that indices are correct! They're switched for C.

    // Form  y := alpha*A**DType*x + y.
    jy = ky;

    if (incX == 1) {
      for (j = 0; j < N; ++j) {
        temp = 0;
        for (i = 0; i < M; ++i) {
          temp += A[j+i*lda]*X[j];
        }
        Y[jy] += *alpha * temp;
        jy += incY;
      }
    } else {
      for (j = 0; j < N; ++j) {
        temp = 0;
        ix = kx;
        for (i = 0; i < M; ++i) {
          temp += A[j+i*lda] * X[ix];
          ix += incX;
        }

        Y[jy] += *alpha * temp;
        jy += incY;
      }
    }
  }

  return true;
}  // end of GEMV

#if (defined(HAVE_CBLAS_H) || defined(HAVE_ATLAS_CBLAS_H))
template <>
inline bool gemv(const enum CBLAS_TRANSPOSE Trans, const int M, const int N, const float* alpha, const float* A, const int lda,
          const float* X, const int incX, const float* beta, float* Y, const int incY) {
  cblas_sgemv(CblasRowMajor, Trans, M, N, *alpha, A, lda, X, incX, *beta, Y, incY);
  return true;
}

template <>
inline bool gemv(const enum CBLAS_TRANSPOSE Trans, const int M, const int N, const double* alpha, const double* A, const int lda,
          const double* X, const int incX, const double* beta, double* Y, const int incY) {
  cblas_dgemv(CblasRowMajor, Trans, M, N, *alpha, A, lda, X, incX, *beta, Y, incY);
  return true;
}

template <>
inline bool gemv(const enum CBLAS_TRANSPOSE Trans, const int M, const int N, const Complex64* alpha, const Complex64* A, const int lda,
          const Complex64* X, const int incX, const Complex64* beta, Complex64* Y, const int incY) {
  cblas_cgemv(CblasRowMajor, Trans, M, N, alpha, A, lda, X, incX, beta, Y, incY);
  return true;
}

template <>
inline bool gemv(const enum CBLAS_TRANSPOSE Trans, const int M, const int N, const Complex128* alpha, const Complex128* A, const int lda,
          const Complex128* X, const int incX, const Complex128* beta, Complex128* Y, const int incY) {
  cblas_zgemv(CblasRowMajor, Trans, M, N, alpha, A, lda, X, incX, beta, Y, incY);
  return true;
}
#endif

}} // end of namespace nm::math

#endif // GEMM_H

max_stars_repo_name: cjfuller/nmatrix

max_issues_repo_name: cjfuller/nmatrix

max_forks_repo_name: cjfuller/nmatrix


text:
#pragma once

#include <memory>
#include <iterator>
#include <cstddef>
#include <gsl/gsl>

namespace dr {

/// \brief round `s` up to the nearest multiple of n
template<typename T>
T round_up(T s, unsigned int n) { return ((s + n - 1) / n) * n; }

template<typename T, typename Allocator = std::allocator<T>>
struct gap_buffer {
  using value_type      = T;
  using allocator_type  = Allocator;
  using size_type       = std::size_t;
  using difference_type = ptrdiff_t;
  using reference       = value_type&;
  using const_reference = const value_type&;
  using pointer         = typename std::allocator_traits<Allocator>::pointer;
  using const_pointer   = typename std::allocator_traits<Allocator>::const_pointer;

  struct const_iterator;

  struct iterator {
    using self_type         = iterator;
    using container_type    = gap_buffer;

    using value_type        = container_type::value_type;
    using difference_type   = container_type::difference_type;
    using reference         = container_type::reference;
    using pointer           = container_type::pointer;
    using iterator_category = std::random_access_iterator_tag;

    explicit iterator(gap_buffer* container = nullptr, difference_type offset = 0)
        : container(container), offset(offset) { }

    operator const_iterator() const {
      return const_iterator(container, offset);
    }

    reference operator [](difference_type i) const {
      return (*container)[offset + i];
    }

    reference operator *() const {
      return (*container)[offset];
    }

    pointer operator ->() const {
      return &(*container)[offset];
    }

    self_type& operator ++() {
      offset++;
      return *this;
    }

    self_type operator ++(int) {
      self_type retval = *this;
      this->operator ++();
      return retval;
    }

    self_type& operator --() {
      offset--;
      return *this;
    }

    self_type operator --(int) {
      self_type retval = *this;
      this->operator --();
      return retval;
    }

    bool operator ==(const self_type& other) const {
      return container == other.container && offset == other.offset;
    }

    bool operator !=(const self_type& other) const {
      return !(*this == other);
    }

    bool operator <(const self_type& other) const {
      Expects(container == other.container);
      return offset < other.offset;
    }

    bool operator >(const self_type& other) const {
      return other < *this;
    }

    bool operator <=(const self_type& other) const {
      return !(other < *this);
    }

    bool operator >=(const self_type& other) const {
      return !(*this < other);
    }

    self_type& operator +=(difference_type n) {
      offset += n;
      return *this;
    }

    friend
    self_type operator +(self_type it, difference_type n) {
      it.offset += n;
      return it;
    }

    friend
    self_type operator +(difference_type n, self_type it) {
      it.offset += n;
      return it;
    }

    self_type& operator -=(difference_type n) {
      offset -= n;
      return *this;
    }

    self_type operator -(difference_type n) const {
      return self_type(container, offset - n);
    }

    difference_type operator -(const self_type& other) const {
      Expects(container == other.container);
      return offset - other.offset;
    }

    friend class gap_buffer;

  private:
    gap_buffer* container;
    difference_type offset;
  };

  // Here we repeat ourselves, that is, DRY principle is violated.
  // We can get rid of the repetition with CRTP. But let's keep it
  // as is.
  struct const_iterator {
    using self_type         = const_iterator;
    using container_type    = gap_buffer;

    using value_type        = container_type::value_type;
    using difference_type   = container_type::difference_type;
    using reference         = container_type::const_reference;
    using pointer           = container_type::const_pointer;
    using iterator_category = std::random_access_iterator_tag;

    explicit const_iterator(const gap_buffer* container = nullptr, difference_type offset = 0)
        : container(container),
          offset(offset) { }

    reference operator [](difference_type i) const {
      return (*container)[offset + i];
    }

    reference operator *() const {
      return (*container)[offset];
    }

    pointer operator ->() const {
      return &(*container)[offset];
    }

    self_type& operator ++() {
      offset++;
      return *this;
    }

    self_type operator ++(int) {
      self_type retval = *this;
      this->operator ++();
      return retval;
    }

    self_type& operator --() {
      offset--;
      return *this;
    }

    self_type operator --(int) {
      self_type retval = *this;
      this->operator --();
      return retval;
    }

    bool operator ==(const self_type& other) const {
      return container == other.container && offset == other.offset;
    }

    bool operator !=(const self_type& other) const {
      return !(*this == other);
    }

    bool operator <(const self_type& other) const {
      Expects(container == other.container);
      return offset < other.offset;
    }

    bool operator >(const self_type& other) const {
      return other < *this;
    }

    bool operator <=(const self_type& other) const {
      return !(other < *this);
    }

    bool operator >=(const self_type& other) const {
      return !(*this < other);
    }

    self_type& operator +=(difference_type n) {
      offset += n;
      return *this;
    }

    friend
    self_type operator +(self_type it, difference_type n) {
      it.offset += n;
      return it;
    }

    friend
    self_type operator +(difference_type n, self_type it) {
      it.offset += n;
      return it;
    }

    self_type& operator -=(difference_type n) {
      offset -= n;
      return *this;
    }

    self_type operator -(difference_type n) const {
      return self_type(container, offset - n);
    }

    difference_type operator -(const self_type& other) const {
      Expects(container == other.container);
      return offset - other.offset;
    }

    friend class gap_buffer;

  private:
    const gap_buffer* container;
    difference_type offset;
  };

  using reverse_iterator       = std::reverse_iterator<iterator>;
  using const_reverse_iterator = std::reverse_iterator<const_iterator>;

private:
  static constexpr float incremental_factor = 0.2;
  static constexpr size_type default_size = 8;
  static constexpr size_type alignment = 8;

public:
  explicit gap_buffer(size_type count = default_size) {
    if (count == 0) {
      start = finish = gap_start = nullptr;
      gap_size = 0;
    }
    else {
      count = round_up(count, alignment);

      start = allocate_and_construct(count);
      finish = start + count;
      gap_start = start;
      gap_size = count;
    }
  }

  gap_buffer(size_type count, const T& value)
      : gap_buffer(count) {
    for (size_type i = 0; i < count; ++i) push_back(value);
  }

  template<typename InputIt>
  gap_buffer(InputIt first, InputIt last) {
    difference_type n = std::distance(first, last);
    size_type len = round_up(std::max(default_size, size_type(n)), alignment);

    start = data_allocator.allocate(len);
    finish = start + len;

    int except_flag = 0;
    try {
      gap_start = std::uninitialized_copy(first, last, start);
      except_flag = 1;
      std::uninitialized_default_construct(gap_start, finish);
    }
    catch (...) {
      switch (except_flag) {
      case 1:
        std::destroy(start, start + n);
        // FALL THROUGH
      case 0:
        data_allocator.deallocate(start, finish - start);
      default:;// DO NOTHING
      }
      throw;
    }

    gap_size = len - n;
  }

  gap_buffer(const gap_buffer& rhs)
      : gap_buffer(rhs.begin(), rhs.end()) { }

  gap_buffer(gap_buffer&& rhs) noexcept
      : gap_buffer(0) { swap(rhs); }

  gap_buffer(std::initializer_list<T> ilist)
      : gap_buffer(ilist.begin(), ilist.end()) { }

  gap_buffer& operator =(const gap_buffer& rhs) {
    gap_buffer temp(rhs);
    swap(temp);
    return *this;
  }

  gap_buffer& operator =(gap_buffer&& rhs) noexcept {
    swap(rhs);
    return *this;
  }

  ~gap_buffer() {
    destroy_and_deallocate(start, finish);
    start = finish = gap_start = nullptr;
    gap_size = 0;
  }

  void swap(gap_buffer& rhs) {
    using std::swap;
    swap(start, rhs.start);
    swap(finish, rhs.finish);
    swap(gap_start, rhs.gap_start);
    swap(gap_size, rhs.gap_size);
  }

  void assign(size_type count, const T& value) { *this = gap_buffer(count, value); }

  template<typename InputIt>
  void assign(InputIt first, InputIt last) { *this = gap_buffer(first, last); }

  void assign(std::initializer_list<T> ilist) { *this = gap_buffer(ilist); }

  allocator_type get_allocator() const { return data_allocator; }

  // ------ basis START HERE ------

  const_reference operator [](size_type pos) const {
    if (start + pos < gap_start) return *(start + pos);
    else return *(start + gap_size + pos);
  }

  iterator erase(const_iterator first, const_iterator last) {
    Expects(first.container == this && last.container == this);
    difference_type num_to_erase = std::distance(first, last);
    relocate_gap(first.offset);
    std::fill_n(gap_start + gap_size, num_to_erase, T{});
    gap_size += num_to_erase;
    return iterator(this, first.offset);
  }

  /// \return iterator to the first inserted element
  template<class InputIt>
  iterator insert(const_iterator pos, InputIt first, InputIt last) {
    Expects(this == pos.container && pos <= end());

    difference_type num_to_insert = std::distance(first, last);
    if (gap_size >= num_to_insert) {
      relocate_gap(pos.offset);

      std::copy(first, last, gap_start);
      gap_start += num_to_insert;
      gap_size -= num_to_insert;
      return iterator(this, pos.offset);
    }
    else {
      size_type old_size = size();
      size_type old_capacity = capacity();
      auto default_delta = size_type(old_capacity * incremental_factor);
      size_type delta = round_up(std::max(default_delta, num_to_insert - gap_size), alignment);
      size_type new_capacity = std::max(old_capacity + delta, default_size);

      gap_buffer temp(new_capacity);

      relocate_gap(pos.offset);
      auto cursor = std::copy(start, gap_start, temp.start);
      cursor = std::copy(first, last, cursor);
      std::copy(gap_start + gap_size, finish, cursor);

      swap(temp);

      gap_start = start + old_size + num_to_insert;
      gap_size = finish - gap_start;

      return iterator(this, pos.offset);
    }
  }

  void reserve(size_type new_cap = 0) {
    if (capacity() >= new_cap) return;
    if (new_cap > max_size()) throw std::length_error("new_cap should be less than max_size()");

    size_type old_size = size();
    size_type new_capacity = round_up(new_cap, alignment);

    gap_buffer temp(new_capacity);
    auto cursor = std::copy(start, gap_start, temp.start);
    std::copy(gap_start + gap_size, finish, cursor);
    swap(temp);

    gap_start = start + old_size;
    gap_size = finish - gap_start;
  }

  size_type size() const noexcept { return finish - start - gap_size; }
  size_type max_size() const noexcept { return size_type(1 << 31); }
  size_type capacity() const noexcept { return finish - start; }

  // ------ basis END HERE ------

  reference operator [](size_type pos) {
    return const_cast<reference>(
        static_cast<const gap_buffer&>(*this)[pos]
    );
  }

  const_reference at(size_type pos) const {
    if (pos >= size()) throw std::out_of_range("index out of range");
    return (*this)[pos];
  }

  reference at(size_type pos) {
    return const_cast<reference>(
        static_cast<const gap_buffer&>(*this).at(pos)
    );
  }

  const_reference front() const { return (*this)[0]; }
  reference front() {
    return const_cast<reference>(
        static_cast<const gap_buffer&>(*this).front()
    );
  }

  const_reference back() const { return (*this)[size() - 1]; }
  reference back() {
    return const_cast<reference>(
        static_cast<const gap_buffer&>(*this).back()
    );
  }

  T* data() noexcept  = delete;
  const T* data() const noexcept = delete;

  [[nodiscard]] bool empty() const noexcept { return size() == 0; }

  void shrink_to_fit() {
    gap_buffer temp(begin(), end());
    swap(temp);
  }

  void clear() { erase(begin(), end()); }

  void resize(size_type count, const value_type& value = value_type{}) {
    if (count < size())
      erase(begin() + count, end());
    else
      insert(end(), count - size(), value);
  }

  iterator insert(const_iterator pos, const T& value) {
    return insert(pos, &value, &value + 1);
  }

  iterator insert(const_iterator pos, T&& value) {
    return insert(pos, std::make_move_iterator(&value), std::make_move_iterator(&value + 1));
  }

  iterator insert(const_iterator pos, size_type count, const T& value) {
    for (size_type i = 0; i < count; ++i)
      pos = insert(pos, value);
    return iterator(this, pos.offset);
  }

  iterator insert(const_iterator pos, std::initializer_list<T> ilist) {
    return insert(pos, ilist.begin(), ilist.end());
  }

  void erase(const_iterator pos) { erase(pos, pos + 1); }

  template<typename ... Args>
  iterator emplace(const_iterator pos, Args&& ... args) {
    return insert(pos, T(std::forward<Args>(args)...));
  }

  template<typename ... Args>
  reference emplace_back(Args&& ... args) {
    return *emplace(end(), std::forward<Args>(args)...);
  }

  void push_back(const T& value) { insert(end(), &value, &value + 1); }
  void push_back(T&& value) { insert(end(), std::make_move_iterator(&value), std::make_move_iterator(&value + 1)); }
  void pop_back() { erase(end() - 1); }

  const_iterator begin() const noexcept { return const_iterator(this); }
  const_iterator cbegin() const noexcept { return begin(); }

  const_iterator end() const noexcept { return const_iterator(this, size()); }
  const_iterator cend() const noexcept { return end(); }

  iterator begin() noexcept { return iterator(this); }
  iterator end() noexcept { return iterator(this, size()); }

  const_reverse_iterator rbegin() const noexcept { return const_reverse_iterator(end()); }
  const_reverse_iterator crbegin() const noexcept { return rbegin(); }

  const_reverse_iterator rend() const noexcept { return const_reverse_iterator(begin()); }
  const_reverse_iterator crend() const noexcept { return rend(); }

  reverse_iterator rbegin() noexcept { return reverse_iterator(end()); }
  reverse_iterator rend() noexcept { return reverse_iterator(begin()); }

  friend
  bool operator ==(const gap_buffer& lhs, const gap_buffer& rhs) {
    return std::equal(lhs.begin(), lhs.end(), rhs.begin(), rhs.end());
  }

  friend
  bool operator !=(const gap_buffer& lhs, const gap_buffer& rhs) {
    return !(lhs == rhs);
  }

  friend
  bool operator <(const gap_buffer& lhs, const gap_buffer& rhs) {
    auto[left, right] = std::mismatch(lhs.begin(), lhs.end(), rhs.begin(), rhs.end());

    if (right == rhs.end()) return false;
    else if (left == lhs.end()) return true;
    else return *left < *right;
  }

  friend
  bool operator >(const gap_buffer& lhs, const gap_buffer& rhs) {
    return rhs < lhs;
  }

  friend
  bool operator <=(const gap_buffer& lhs, const gap_buffer& rhs) {
    return !(rhs < lhs);
  }

  friend
  bool operator >=(const gap_buffer& lhs, const gap_buffer& rhs) {
    return !(lhs < rhs);
  }

  // additional
  template<typename InputIt>
  void append(InputIt first, InputIt last) {
    insert(end(), first, last);
  }

  void append(const T& value) { push_back(value); }

  void append(T&& value) { push_back(std::move(value)); }

  template<typename InputIt>
  void replace(const_iterator f1, const_iterator l1, InputIt f2, InputIt l2) {
    auto cursor = erase(f1, l1);
    insert(cursor, f2, l2);
  }

  template<typename InputIt>
  void replace(const_iterator pos, InputIt first, InputIt last) {
    replace(pos, pos + 1, first, last);
  }

  gap_buffer substr(const_iterator first, const_iterator last) const {
    return substr_impl<gap_buffer>(first, last);
  }

protected:

  template<typename U>
  U substr_impl(const_iterator first, const_iterator last) const {
    return U(first, last);
  }

  void relocate_gap(difference_type offset) {
    if (gap_start != start + offset) {
      if (gap_start < start + offset)
        std::move(gap_start /**/ + gap_size,
                  start + offset + gap_size,
                  gap_start);
      else
        std::move_backward(start + offset,
                           gap_start,
                           gap_start + gap_size);

      gap_start = start + offset;
    }
  }

  pointer allocate_and_construct(size_type n) {
    pointer result = data_allocator.allocate(n);
    std::uninitialized_default_construct_n(result, n);
    return result;
  }

  void destroy_and_deallocate(pointer start, pointer finish) {
    std::destroy(start, finish);
    if (start)
      data_allocator.deallocate(start, finish - start);
  }

private:
  Allocator data_allocator;

  pointer start;
  pointer finish;
  pointer gap_start;
  size_type gap_size;
};

template<typename T, typename Allocator>
void swap(gap_buffer<T, Allocator>& lhs, gap_buffer<T, Allocator>& rhs) {
  lhs.swap(rhs);
}

}

max_stars_repo_name: lie-yan/gapbuffer

max_issues_repo_name: lie-yan/gapbuffer

max_forks_repo_name: lie-yan/gapbuffer


text:
#include <jni.h>
#include <assert.h>
#include <lapacke.h>

/* Cholesky */

JNIEXPORT jint Java_JAMAJni_CholeskyDecomposition_dpotrf (JNIEnv *env, jclass klass, jint matrix_layout, jchar uplo, jint n, jdoubleArray a, jint lda){

    double *aElems;
    int info;

    aElems = (*env)-> GetDoubleArrayElements (env, a, NULL);

    assert(aElems);

    info = LAPACKE_dpotrf((int) matrix_layout, (char) uplo, (lapack_int) n, aElems, (lapack_int) lda);

    (*env)-> ReleaseDoubleArrayElements (env, a, aElems, 0);

    return info;

}

JNIEXPORT jint Java_JAMAJni_CholeskyDecomposition_dpotri(JNIEnv *env, jclass klass, jint matrix_layout, jchar uplo, jint n, jdoubleArray a, jint lda){

    double *aElems;
    int info;

    aElems = (*env)-> GetDoubleArrayElements (env, a, NULL);

    assert(aElems);

    info = LAPACKE_dpotri((int) matrix_layout, (char) uplo, (lapack_int) n, aElems, (lapack_int) lda);

    (*env)-> ReleaseDoubleArrayElements (env, a, aElems, 0);

    return info;

}

JNIEXPORT jint Java_JAMAJni_CholeskyDecomposition_dpotrs (JNIEnv *env, jclass klass, jint matrix_layout, jchar uplo, jint n, jint nrhs, jdoubleArray a, jint lda, jdoubleArray b, jint ldb){

    double *aElems, *bElems;
    int info;

    aElems = (*env)-> GetDoubleArrayElements (env, a, NULL);
    bElems = (*env)-> GetDoubleArrayElements (env, b, NULL);

    assert(aElems && bElems);

    info = LAPACKE_dpotrs ((int) matrix_layout, (char) uplo, (lapack_int) n, (lapack_int) nrhs, aElems, (lapack_int) lda, bElems, (lapack_int) ldb);

    (*env)-> ReleaseDoubleArrayElements (env, a, aElems, JNI_ABORT);
    (*env)-> ReleaseDoubleArrayElements (env, b, bElems, 0);

    return info;

}

ext: c


text:
static char help[] =
"Solves time-dependent heat equation in 2D using TS.  Option prefix -ht_.\n"
"Equation is  u_t = D_0 laplacian u + f.  Domain is (0,1) x (0,1).\n"
"Boundary conditions are non-homogeneous Neumann in x and periodic in y.\n"
"Energy is conserved (for these particular conditions/source) and an extra\n"
"'monitor' is demonstrated.  Discretization is by centered finite differences.\n"
"Converts the PDE into a system  X_t = G(t,X) (PETSc type 'nonlinear') by\n"
"method of lines.  Uses backward Euler time-stepping by default.\n";

#include <petsc.h>

typedef struct {
  PetscReal D0;    // conductivity
} HeatCtx;

static PetscReal f_source(PetscReal x, PetscReal y) {
    return 3.0 * PetscExpReal(-25.0 * (x-0.6) * (x-0.6))
               * PetscSinReal(2.0*PETSC_PI*y);
}

static PetscReal gamma_neumann(PetscReal y) {
    return PetscSinReal(6.0 * PETSC_PI * y);
}

extern PetscErrorCode Spacings(DMDALocalInfo*, PetscReal*, PetscReal*);
extern PetscErrorCode EnergyMonitor(TS, PetscInt, PetscReal, Vec, void*);
extern PetscErrorCode FormRHSFunctionLocal(DMDALocalInfo*, PetscReal, PetscReal**,
                                           PetscReal**, HeatCtx*);
extern PetscErrorCode FormRHSJacobianLocal(DMDALocalInfo*, PetscReal, PetscReal**,
                                           Mat, Mat, HeatCtx*);

int main(int argc,char **argv) {
  PetscErrorCode ierr;
  HeatCtx        user;
  TS             ts;
  Vec            u;
  DM             da;
  DMDALocalInfo  info;
  PetscReal      t0, tf;
  PetscBool      monitorenergy = PETSC_FALSE;

  ierr = PetscInitialize(&argc,&argv,NULL,help); if (ierr) return ierr;

  user.D0  = 1.0;
  ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "ht_", "options for heat", ""); CHKERRQ(ierr);
  ierr = PetscOptionsReal("-D0","constant thermal diffusivity",
           "heat.c",user.D0,&user.D0,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsBool("-monitor","also display total heat energy at each step",
           "heat.c",monitorenergy,&monitorenergy,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsEnd(); CHKERRQ(ierr);

//STARTDMDASETUP
  ierr = DMDACreate2d(PETSC_COMM_WORLD,
      DM_BOUNDARY_NONE, DM_BOUNDARY_PERIODIC, DMDA_STENCIL_STAR,
      5,4,PETSC_DECIDE,PETSC_DECIDE,  // default to hx=hx=0.25 grid
      1,1,                            // degrees of freedom, stencil width
      NULL,NULL,&da); CHKERRQ(ierr);
  ierr = DMSetFromOptions(da); CHKERRQ(ierr);
  ierr = DMSetUp(da); CHKERRQ(ierr);
  ierr = DMCreateGlobalVector(da,&u); CHKERRQ(ierr);
//ENDDMDASETUP

//STARTTSSETUP
  ierr = TSCreate(PETSC_COMM_WORLD,&ts); CHKERRQ(ierr);
  ierr = TSSetProblemType(ts,TS_NONLINEAR); CHKERRQ(ierr);
  ierr = TSSetDM(ts,da); CHKERRQ(ierr);
  ierr = TSSetApplicationContext(ts,&user); CHKERRQ(ierr);
  ierr = DMDATSSetRHSFunctionLocal(da,INSERT_VALUES,
           (DMDATSRHSFunctionLocal)FormRHSFunctionLocal,&user); CHKERRQ(ierr);
  ierr = DMDATSSetRHSJacobianLocal(da,
           (DMDATSRHSJacobianLocal)FormRHSJacobianLocal,&user); CHKERRQ(ierr);
  if (monitorenergy) {
      ierr = TSMonitorSet(ts,EnergyMonitor,&user,NULL); CHKERRQ(ierr);
  }
  ierr = TSSetType(ts,TSBDF); CHKERRQ(ierr);
  ierr = TSSetTime(ts,0.0); CHKERRQ(ierr);
  ierr = TSSetMaxTime(ts,0.1); CHKERRQ(ierr);
  ierr = TSSetTimeStep(ts,0.001); CHKERRQ(ierr);
  ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP); CHKERRQ(ierr);
  ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
//ENDTSSETUP

  // report on set up
  ierr = TSGetTime(ts,&t0); CHKERRQ(ierr);
  ierr = TSGetMaxTime(ts,&tf); CHKERRQ(ierr);
  ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,
           "solving on %d x %d grid for t0=%g to tf=%g ...\n",
           info.mx,info.my,t0,tf); CHKERRQ(ierr);

  // solve
  ierr = VecSet(u,0.0); CHKERRQ(ierr);   // initial condition
  ierr = TSSolve(ts,u); CHKERRQ(ierr);

  VecDestroy(&u);  TSDestroy(&ts);  DMDestroy(&da);
  return PetscFinalize();
}

PetscErrorCode Spacings(DMDALocalInfo *info, PetscReal *hx, PetscReal *hy) {
    if (hx)  *hx = 1.0 / (PetscReal)(info->mx-1);
    if (hy)  *hy = 1.0 / (PetscReal)(info->my);   // periodic direction
    return 0;
}

//STARTMONITOR
PetscErrorCode EnergyMonitor(TS ts, PetscInt step, PetscReal time, Vec u,
                             void *ctx) {
    PetscErrorCode ierr;
    HeatCtx        *user = (HeatCtx*)ctx;
    PetscReal      lenergy = 0.0, energy, dt, hx, hy, **au;
    PetscInt       i,j;
    MPI_Comm       com;
    DM             da;
    DMDALocalInfo  info;

    ierr = TSGetDM(ts,&da); CHKERRQ(ierr);
    ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(da,u,&au); CHKERRQ(ierr);
    for (j = info.ys; j < info.ys + info.ym; j++) {
        for (i = info.xs; i < info.xs + info.xm; i++) {
            if ((i == 0) || (i == info.mx-1))
                lenergy += 0.5 * au[j][i];
            else
                lenergy += au[j][i];
        }
    }
    ierr = DMDAVecRestoreArrayRead(da,u,&au); CHKERRQ(ierr);
    ierr = Spacings(&info,&hx,&hy); CHKERRQ(ierr);
    lenergy *= hx * hy;
    ierr = PetscObjectGetComm((PetscObject)(da),&com); CHKERRQ(ierr);
    ierr = MPI_Allreduce(&lenergy,&energy,1,MPIU_REAL,MPIU_SUM,com); CHKERRQ(ierr);
    ierr = TSGetTimeStep(ts,&dt); CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"  energy = %9.2e     nu = %8.4f\n",
                energy,user->D0*dt/(hx*hy)); CHKERRQ(ierr);
    return 0;
}
//ENDMONITOR

//STARTRHSFUNCTION
PetscErrorCode FormRHSFunctionLocal(DMDALocalInfo *info,
                                    PetscReal t, PetscReal **au,
                                    PetscReal **aG, HeatCtx *user) {
  PetscErrorCode ierr;
  PetscInt   i, j, mx = info->mx;
  PetscReal  hx, hy, x, y, ul, ur, uxx, uyy;

  ierr = Spacings(info,&hx,&hy); CHKERRQ(ierr);
  for (j = info->ys; j < info->ys + info->ym; j++) {
      y = hy * j;
      for (i = info->xs; i < info->xs + info->xm; i++) {
          x = hx * i;
          // apply Neumann b.c.s
          ul = (i == 0) ? au[j][i+1] + 2.0 * hx * gamma_neumann(y)
                        : au[j][i-1];
          ur = (i == mx-1) ? au[j][i-1] : au[j][i+1];
          uxx = (ul - 2.0 * au[j][i]+ ur) / (hx*hx);
          // DMDA is periodic in y
          uyy = (au[j-1][i] - 2.0 * au[j][i]+ au[j+1][i]) / (hy*hy);
          aG[j][i] = user->D0 * (uxx + uyy) + f_source(x,y);
      }
  }
  return 0;
}
//ENDRHSFUNCTION

//STARTRHSJACOBIAN
PetscErrorCode FormRHSJacobianLocal(DMDALocalInfo *info,
                                    PetscReal t, PetscReal **au,
                                    Mat J, Mat P, HeatCtx *user) {
    PetscErrorCode ierr;
    PetscInt         i, j, ncols;
    const PetscReal  D = user->D0;
    PetscReal        hx, hy, hx2, hy2, v[5];
    MatStencil       col[5],row;

    ierr = Spacings(info,&hx,&hy); CHKERRQ(ierr);
    hx2 = hx * hx;  hy2 = hy * hy;
    for (j = info->ys; j < info->ys+info->ym; j++) {
        row.j = j;  col[0].j = j;
        for (i = info->xs; i < info->xs+info->xm; i++) {
            // set up a standard 5-point stencil for the row
            row.i = i;
            col[0].i = i;
            v[0] = - 2.0 * D * (1.0 / hx2 + 1.0 / hy2);
            col[1].j = j-1;  col[1].i = i;    v[1] = D / hy2;
            col[2].j = j+1;  col[2].i = i;    v[2] = D / hy2;
            col[3].j = j;    col[3].i = i-1;  v[3] = D / hx2;
            col[4].j = j;    col[4].i = i+1;  v[4] = D / hx2;
            ncols = 5;
            // if at the boundary, edit the row back to 4 nonzeros
            if (i == 0) {
                ncols = 4;
                col[3].j = j;  col[3].i = i+1;  v[3] = 2.0 * D / hx2;
            } else if (i == info->mx-1) {
                ncols = 4;
                col[3].j = j;  col[3].i = i-1;  v[3] = 2.0 * D / hx2;
            }
            ierr = MatSetValuesStencil(P,1,&row,ncols,col,v,INSERT_VALUES); CHKERRQ(ierr);
        }
    }

    ierr = MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    ierr = MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    if (J != P) {
        ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
        ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    }
    return 0;
}
//ENDRHSJACOBIAN


text:
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multiroots.h>

struct rparams
  {
    double a;
    double b;
  };

int
rosenbrock_f (const gsl_vector * x, void *params, gsl_vector * f)
{
  double a = ((struct rparams *) params)->a;
  double b = ((struct rparams *) params)->b;

  double x0 = gsl_vector_get (x, 0);
  double x1 = gsl_vector_get (x, 1);

  double y0 = a * (1 - x0);
  double y1 = b * (x1 - x0 * x0);

  gsl_vector_set (f, 0, y0);
  gsl_vector_set (f, 1, y1);

  return GSL_SUCCESS;
}

#ifdef DERIV
int
rosenbrock_df (const gsl_vector * x, void *params, gsl_matrix * df)
{
  double a = ((struct rparams *) params)->a;
  double b = ((struct rparams *) params)->b;

  double x0 = gsl_vector_get (x, 0);

  double df00 = -a;
  double df01 = 0;
  double df10 = -2 * b  * x0;
  double df11 = b;

  gsl_matrix_set (df, 0, 0, df00);
  gsl_matrix_set (df, 0, 1, df01);
  gsl_matrix_set (df, 1, 0, df10);
  gsl_matrix_set (df, 1, 1, df11);

  return GSL_SUCCESS;
}

int
rosenbrock_fdf (const gsl_vector * x, void *params,
                gsl_vector * f, gsl_matrix * df)
{
  rosenbrock_f (x, params, f);
  rosenbrock_df (x, params, df);

  return GSL_SUCCESS;
}
#endif

#ifdef DERIV
#define SOLVER gsl_multiroot_fdfsolver
#define SOLVER_TYPE gsl_multiroot_fdfsolver_type
#else
#define SOLVER gsl_multiroot_fsolver
#define SOLVER_TYPE gsl_multiroot_fsolver_type
#endif

int
main (void)
{
  const SOLVER_TYPE *T;
  SOLVER *s;

  int status;
  size_t i, iter = 0;

  const size_t n = 2;
  struct rparams p = {1.0, 10.0};
#ifdef DERIV
  gsl_multiroot_function_fdf f = {&rosenbrock_f, &rosenbrock_df, &rosenbrock_fdf, n, &p};
#else
  gsl_multiroot_function f = {&rosenbrock_f, n, &p};
#endif

  double x_init[2] = {-10.0, -5.0};
  gsl_vector *x = gsl_vector_alloc (n);

  gsl_vector_set (x, 0, x_init[0]);
  gsl_vector_set (x, 1, x_init[1]);

#ifdef DERIV
  T = gsl_multiroot_fdfsolver_gnewton;
  s = gsl_multiroot_fdfsolver_alloc (T, &f, x);
#else
  T = gsl_multiroot_fsolver_hybrids;
  s = gsl_multiroot_fsolver_alloc (T, &f, x);
#endif

  print_state (iter, s);

  do
    {
      iter++;
#ifdef DERIV
      status = gsl_multiroot_fdfsolver_iterate (s);
#else
      status = gsl_multiroot_fsolver_iterate (s);
#endif

      print_state (iter, s);

      if (status)
        break;

      status = gsl_multiroot_test_residual (s->f, 0.0000001);
    }
  while (status == GSL_CONTINUE && iter < 1000);

  printf ("status = %s\n", gsl_strerror (status));

#ifdef DERIV
  gsl_multiroot_fdfsolver_free (s);
#else
  gsl_multiroot_fsolver_free (s);
#endif

  gsl_vector_free (x);
}

int
print_state (size_t iter, SOLVER * s)
{
  printf ("iter = %3u x = % 15.8f % 15.8f  f(x) = % .3e % .3e\n",
          iter,
          gsl_vector_get (s->x, 0), gsl_vector_get (s->x, 1),
          gsl_vector_get (s->f, 0), gsl_vector_get (s->f, 1));
}


text:
static char help[] =
"Simple ODE system to test implicit and IMEX methods when *both* RHSFunction()\n"
"and IFunction() are supplied.  Jacobians are supplied as well.  The problem\n"
"has form\n"
"   F(t,u,dudt) = G(t,u)\n"
"where  u(t) = (x(t),y(t))  and dF/d(dudt) is invertible, so the problem is a\n"
"2D ODE IVP.  Functions F (IFunction) and G (RHSFunction) are actually linear,\n"
"but the ProblemType is set to TS_NONLINEAR type.  By default the system is\n"
"   x' + y' +  x + 2y =  x + 8y\n"
"        y' + 3x + 4y = 4x + 4y\n"
"The implicit methods like BEULER, THETA, CN, BDF should be able to handle\n"
"this system.  Apparently ROSW can also handle it.  However, ARKIMEX and EIMEX\n"
"require dF/d(dudt) to be the identity.  (But consider option\n"
"-ts_arkimex_fully_implicit.)  If option -identity_in_F is given then an\n"
"equivalent system is formed,\n"
"   x'      +  x + 2y = 8y\n"
"        y' + 3x + 4y = 4x + 4y\n"
"Both of these systems are trivial transformations of the system\n"
"x'=6y-x, y'=x.  With the initial conditions chosen below, x(0)=1 and y(0)=3,\n"
"the exact solution to the above systems is x(t) = -3 e^{-3t} + 4 e^{2t},\n"
"y(t) = e^{-3t} + 2 e^{2t}, and we compute the final numerical error\n"
"accordingly.  Default type is BDF.  See the manual pages for the various\n"
"methods (e.g. TSROSW, TSARKIMEX, TSEIMEX, ...) for further information.\n\n";

#include <petsc.h>

typedef struct {
  PetscBool  identity_in_F;
} Ctx;

extern PetscErrorCode FormIFunction(TS,PetscReal,Vec,Vec,Vec F,void*);
extern PetscErrorCode FormIJacobian(TS,PetscReal,Vec,Vec,PetscReal,Mat,Mat,void*);
extern PetscErrorCode FormRHSFunction(TS,PetscReal,Vec,Vec,void*);
extern PetscErrorCode FormRHSJacobian(TS,PetscReal,Vec,Mat,Mat,void*);

int main(int argc,char **argv)
{
  PetscErrorCode ierr;
  TS             ts;
  Vec            u, uexact;
  Mat            JF,JG;
  Ctx            user;
  PetscReal      tf,xf,yf,errnorm;
  PetscInt       steps;

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

  user.identity_in_F = PETSC_FALSE;
  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"",
           "Simple F(t,u,u')=G(t,u) ODE system.","TS"); CHKERRQ(ierr);
  ierr = PetscOptionsBool("-identity_in_F","set up system so the dF/d(dudt) = I",
           "ex54.c",user.identity_in_F,&(user.identity_in_F),NULL);CHKERRQ(ierr);
  ierr = PetscOptionsEnd(); CHKERRQ(ierr);

  ierr = VecCreate(PETSC_COMM_WORLD,&u); CHKERRQ(ierr);
  ierr = VecSetSizes(u,PETSC_DECIDE,2); CHKERRQ(ierr);
  ierr = VecSetFromOptions(u); CHKERRQ(ierr);

  ierr = MatCreate(PETSC_COMM_WORLD,&JF); CHKERRQ(ierr);
  ierr = MatSetSizes(JF,PETSC_DECIDE,PETSC_DECIDE,2,2); CHKERRQ(ierr);
  ierr = MatSetFromOptions(JF); CHKERRQ(ierr);
  ierr = MatSetUp(JF); CHKERRQ(ierr);
  ierr = MatDuplicate(JF,MAT_DO_NOT_COPY_VALUES,&JG); CHKERRQ(ierr);

  ierr = TSCreate(PETSC_COMM_WORLD,&ts); CHKERRQ(ierr);
  ierr = TSSetApplicationContext(ts,&user); CHKERRQ(ierr);

  ierr = TSSetIFunction(ts,NULL,FormIFunction,&user); CHKERRQ(ierr);
  ierr = TSSetIJacobian(ts,JF,JF,FormIJacobian,&user); CHKERRQ(ierr);

  ierr = TSSetRHSFunction(ts,NULL,FormRHSFunction,&user); CHKERRQ(ierr);
  ierr = TSSetRHSJacobian(ts,JG,JG,FormRHSJacobian,&user); CHKERRQ(ierr);

  ierr = TSSetProblemType(ts,TS_NONLINEAR); CHKERRQ(ierr);
  // helps with ARKIMEX if you also set -ts_arkimex_fully_implicit:
  // ierr = TSSetEquationType(ts,TS_EQ_IMPLICIT); CHKERRQ(ierr);
  ierr = TSSetType(ts,TSBDF); CHKERRQ(ierr);
  ierr = TSSetTime(ts,0.0); CHKERRQ(ierr);
  ierr = TSSetMaxTime(ts,1.0); CHKERRQ(ierr);
  ierr = TSSetTimeStep(ts,0.01); CHKERRQ(ierr);
  ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP); CHKERRQ(ierr);
  ierr = TSSetFromOptions(ts);CHKERRQ(ierr);

  ierr = VecSetValue(u,0,1.0,INSERT_VALUES); CHKERRQ(ierr);
  ierr = VecSetValue(u,1,3.0,INSERT_VALUES); CHKERRQ(ierr);
  ierr = VecAssemblyBegin(u); CHKERRQ(ierr);
  ierr = VecAssemblyEnd(u); CHKERRQ(ierr);
  ierr = TSSolve(ts,u); CHKERRQ(ierr);

  ierr = TSGetStepNumber(ts,&steps); CHKERRQ(ierr);
  ierr = TSGetTime(ts,&tf); CHKERRQ(ierr);
  xf = -3.0 * PetscExpReal(-3.0*tf) + 4.0 * PetscExpReal(2.0*tf);
  yf = PetscExpReal(-3.0*tf) + 2.0 * PetscExpReal(2.0*tf);
  //ierr = PetscPrintf(PETSC_COMM_WORLD,
  //         "xf = %.6f, yf = %.6f\n",xf,yf); CHKERRQ(ierr);
  ierr = VecDuplicate(u,&uexact); CHKERRQ(ierr);
  ierr = VecSetValue(uexact,0,xf,INSERT_VALUES); CHKERRQ(ierr);
  ierr = VecSetValue(uexact,1,yf,INSERT_VALUES); CHKERRQ(ierr);
  ierr = VecAssemblyBegin(uexact); CHKERRQ(ierr);
  ierr = VecAssemblyEnd(uexact); CHKERRQ(ierr);
  ierr = VecAXPY(u,-1.0,uexact); CHKERRQ(ierr); // u <- u + (-1.0) uexact
  ierr = VecNorm(u,NORM_2,&errnorm); CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,
           "error norm at tf = %.6f from %d steps:  |u-u_exact| =  %.5e\n",
           tf,steps,errnorm); CHKERRQ(ierr);

  VecDestroy(&u);  VecDestroy(&uexact);
  MatDestroy(&JF);  MatDestroy(&JG);
  TSDestroy(&ts);
  return PetscFinalize();
}

PetscErrorCode FormIFunction(TS ts, PetscReal t, Vec u, Vec dudt,
                             Vec F, void *user) {
    PetscErrorCode  ierr;
    const PetscReal *au, *adudt;
    PetscReal       *aF;
    PetscBool       flag = ((Ctx*)user)->identity_in_F;

    ierr = VecGetArrayRead(u,&au); CHKERRQ(ierr);
    ierr = VecGetArrayRead(dudt,&adudt); CHKERRQ(ierr);
    ierr = VecGetArray(F,&aF);
    if (flag) {
        aF[0] = adudt[0] + au[0] + 2.0 * au[1];
        aF[1] = adudt[1] + 3.0 * au[0] + 4.0 * au[1];
    } else {
        aF[0] = adudt[0] + adudt[1] + au[0] + 2.0 * au[1];
        aF[1] = adudt[1] + 3.0 * au[0] + 4.0 * au[1];
    }
    ierr = VecRestoreArrayRead(u,&au); CHKERRQ(ierr);
    ierr = VecRestoreArrayRead(dudt,&adudt); CHKERRQ(ierr);
    ierr = VecRestoreArray(F,&aF); CHKERRQ(ierr);
    return 0;
}

// computes  J = dF/du + a dF/d(dudt)
PetscErrorCode FormIJacobian(TS ts, PetscReal t, Vec u, Vec dudt,
                             PetscReal a, Mat J, Mat P, void *user) {
    PetscErrorCode ierr;
    PetscInt       row[2] = {0, 1},  col[2] = {0, 1};
    PetscReal      v[4];
    PetscBool      flag = ((Ctx*)user)->identity_in_F;

    if (flag) {
        v[0] = a + 1.0;    v[1] = 2.0;
        v[2] = 3.0;        v[3] = a + 4.0;
    } else {
        v[0] = a + 1.0;    v[1] = a + 2.0;
        v[2] = 3.0;        v[3] = a + 4.0;
    }
    ierr = MatSetValues(P,2,row,2,col,v,INSERT_VALUES); CHKERRQ(ierr);
    ierr = MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    ierr = MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    if (J != P) {
        ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
        ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    }
    return 0;
}

PetscErrorCode FormRHSFunction(TS ts, PetscReal t, Vec u,
                               Vec G, void *user) {
    PetscErrorCode  ierr;
    const PetscReal *au;
    PetscReal       *aG;
    PetscBool       flag = ((Ctx*)user)->identity_in_F;

    ierr = VecGetArrayRead(u,&au); CHKERRQ(ierr);
    ierr = VecGetArray(G,&aG); CHKERRQ(ierr);
    if (flag) {
        aG[0] = 8.0 * au[1];
        aG[1] = 4.0 * au[0] + 4.0 * au[1];
    } else {
        aG[0] = au[0] + 8.0 * au[1];
        aG[1] = 4.0 * au[0] + 4.0 * au[1];
    }
    ierr = VecRestoreArrayRead(u,&au); CHKERRQ(ierr);
    ierr = VecRestoreArray(G,&aG); CHKERRQ(ierr);
    return 0;
}

PetscErrorCode FormRHSJacobian(TS ts, PetscReal t, Vec u, Mat J, Mat P,
                               void *user) {
    PetscErrorCode ierr;
    PetscInt       row[2] = {0, 1},  col[2] = {0, 1};
    PetscReal      v[4];
    PetscBool      flag = ((Ctx*)user)->identity_in_F;

    if (flag) {
        v[0] = 0.0;    v[1] = 8.0;
        v[2] = 4.0;    v[3] = 4.0;
    } else {
        v[0] = 1.0;    v[1] = 8.0;
        v[2] = 4.0;    v[3] = 4.0;
    }
    ierr = MatSetValues(P,2,row,2,col,v,INSERT_VALUES); CHKERRQ(ierr);
    ierr = MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    ierr = MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    if (J != P) {
        ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
        ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    }
    return 0;
}


text:
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_cblas.h>

#include "tests.h"

void
test_spr (void) {
const double flteps = 1e-4, dbleps = 1e-6;
  {
   int order = 101;
   int uplo = 121;
   int N = 2;
   float alpha = -0.3f;
   float Ap[] = { -0.764f, -0.257f, -0.064f };
   float X[] = { 0.455f, -0.285f };
   int incX = -1;
   float Ap_expected[] = { -0.788367f, -0.218097f, -0.126108f };
   cblas_sspr(order, uplo, N, alpha, X, incX, Ap);
   {
     int i;
     for (i = 0; i < 3; i++) {
       gsl_test_rel(Ap[i], Ap_expected[i], flteps, "sspr(case 1426)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int N = 2;
   float alpha = -0.3f;
   float Ap[] = { -0.764f, -0.257f, -0.064f };
   float X[] = { 0.455f, -0.285f };
   int incX = -1;
   float Ap_expected[] = { -0.788367f, -0.218097f, -0.126108f };
   cblas_sspr(order, uplo, N, alpha, X, incX, Ap);
   {
     int i;
     for (i = 0; i < 3; i++) {
       gsl_test_rel(Ap[i], Ap_expected[i], flteps, "sspr(case 1427)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int N = 2;
   float alpha = -0.3f;
   float Ap[] = { -0.764f, -0.257f, -0.064f };
   float X[] = { 0.455f, -0.285f };
   int incX = -1;
   float Ap_expected[] = { -0.788367f, -0.218097f, -0.126108f };
   cblas_sspr(order, uplo, N, alpha, X, incX, Ap);
   {
     int i;
     for (i = 0; i < 3; i++) {
       gsl_test_rel(Ap[i], Ap_expected[i], flteps, "sspr(case 1428)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int N = 2;
   float alpha = -0.3f;
   float Ap[] = { -0.764f, -0.257f, -0.064f };
   float X[] = { 0.455f, -0.285f };
   int incX = -1;
   float Ap_expected[] = { -0.788367f, -0.218097f, -0.126108f };
   cblas_sspr(order, uplo, N, alpha, X, incX, Ap);
   {
     int i;
     for (i = 0; i < 3; i++) {
       gsl_test_rel(Ap[i], Ap_expected[i], flteps, "sspr(case 1429)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 121;
   int N = 2;
   double alpha = -1;
   double Ap[] = { 0.819, 0.175, -0.809 };
   double X[] = { -0.645, -0.222 };
   int incX = -1;
   double Ap_expected[] = { 0.769716, 0.03181, -1.225025 };
   cblas_dspr(order, uplo, N, alpha, X, incX, Ap);
   {
     int i;
     for (i = 0; i < 3; i++) {
       gsl_test_rel(Ap[i], Ap_expected[i], dbleps, "dspr(case 1430)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int N = 2;
   double alpha = -1;
   double Ap[] = { 0.819, 0.175, -0.809 };
   double X[] = { -0.645, -0.222 };
   int incX = -1;
   double Ap_expected[] = { 0.769716, 0.03181, -1.225025 };
   cblas_dspr(order, uplo, N, alpha, X, incX, Ap);
   {
     int i;
     for (i = 0; i < 3; i++) {
       gsl_test_rel(Ap[i], Ap_expected[i], dbleps, "dspr(case 1431)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int N = 2;
   double alpha = -1;
   double Ap[] = { 0.819, 0.175, -0.809 };
   double X[] = { -0.645, -0.222 };
   int incX = -1;
   double Ap_expected[] = { 0.769716, 0.03181, -1.225025 };
   cblas_dspr(order, uplo, N, alpha, X, incX, Ap);
   {
     int i;
     for (i = 0; i < 3; i++) {
       gsl_test_rel(Ap[i], Ap_expected[i], dbleps, "dspr(case 1432)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int N = 2;
   double alpha = -1;
   double Ap[] = { 0.819, 0.175, -0.809 };
   double X[] = { -0.645, -0.222 };
   int incX = -1;
   double Ap_expected[] = { 0.769716, 0.03181, -1.225025 };
   cblas_dspr(order, uplo, N, alpha, X, incX, Ap);
   {
     int i;
     for (i = 0; i < 3; i++) {
       gsl_test_rel(Ap[i], Ap_expected[i], dbleps, "dspr(case 1433)");
     }
   };
  };

}

max_stars_repo_name: utdsimmons/ohpc

max_issues_repo_name: utdsimmons/ohpc

max_forks_repo_name: utdsimmons/ohpc


text:
/*
FRACTAL - A program growing fractals to benchmark parallelization and drawing
libraries.

Copyright 2009-2021, Javier Burguete Tolosa.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

1. Redistributions of source code must retain the above copyright notice, this
  list of conditions and the following disclaimer.

2. Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY Javier Burguete Tolosa ``AS IS'' AND ANY EXPRESS
OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
EVENT SHALL Javier Burguete Tolosa OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

/**
 * \file main.c
 * \brief Source file to define the main function.
 * \author Javier Burguete Tolosa.
 * \copyright Copyright 2009-2021, Javier Burguete Tolosa.
 */
#define _GNU_SOURCE
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <unistd.h>
#include <libintl.h>
#include <gsl/gsl_rng.h>
#include <glib.h>
#include <png.h>
#include <ft2build.h>
#include FT_FREETYPE_H
#include <gtk/gtk.h>
#include <GL/glew.h>
#if HAVE_FREEGLUT
#include <GL/freeglut.h>
#elif HAVE_SDL
#include <SDL.h>
#elif HAVE_GLFW
#include <GLFW/glfw3.h>
#endif

#include "config.h"
#include "fractal.h"
#include "image.h"
#include "text.h"
#include "graphic.h"
#include "draw.h"
#include "simulator.h"

#if HAVE_SDL
SDL_Window *window;             ///< SDL window.
#elif HAVE_GLFW
GLFWwindow *window;             ///< GLFW window.
#endif

/**
 * Function to do one main loop iteration.
 */
void
main_iteration ()
{
  GMainContext *context;
  context = g_main_context_default ();
  while (g_main_context_pending (context))
    g_main_context_iteration (context, 0);
}

/**
 * Function to do the main loop.
 */
void
main_loop ()
{

#if HAVE_FREEGLUT

  // Passing the GTK+ signals to the FreeGLUT main loop
  glutIdleFunc (main_iteration);
  // Setting our draw resize function as the FreeGLUT reshape function
  glutReshapeFunc (draw_resize);
  // Setting our draw function as the FreeGLUT display function
  glutDisplayFunc (draw);
  // FreeGLUT main loop
  glutMainLoop ();

#elif HAVE_SDL

  SDL_Event event[1];
  while (1)
    {
      main_iteration ();
      while (SDL_PollEvent (event))
        {
          if (event->type == SDL_QUIT)
            return;
          if (event->type == SDL_WINDOWEVENT
              && event->window.event == SDL_WINDOWEVENT_SIZE_CHANGED)
            draw_resize (event->window.data1, event->window.data2);
        }
      draw ();
    }

#elif HAVE_GLFW

  while (!glfwWindowShouldClose (window))
    {
      main_iteration ();
      glfwPollEvents ();
      draw ();
    }

#else

  g_main_loop_run (dialog_simulator->loop);
  g_main_loop_unref (dialog_simulator->loop);

#endif
}

/**
 * Main function.
 *
 * \return 0 on success.
 */
int
main (int argn,                 ///< Arguments number.
      char **argc)              ///< Array of arguments.
{
  if (argn > 2)
    {
      printf ("Bad arguments number\n");
      return 1;
    }

// PARALELLIZING INIT
  nthreads = threads_number ();
// END

  // Initing locales
#if DEBUG
  printf ("Initing locales\n");
  fflush (stdout);
#endif
  bindtextdomain ("fractal", "./po");
  bind_textdomain_codeset ("fractal", "UTF-8");
  textdomain ("fractal");

  // Initing graphic window
#if HAVE_FREEGLUT

#if DEBUG
  printf ("Initing FreeGLUT window\n");
  fflush (stdout);
#endif
  glutInit (&argn, argc);
  glutInitDisplayMode (GLUT_DEPTH | GLUT_DOUBLE | GLUT_RGBA);
  glutInitWindowSize (window_width, window_height);
  glutCreateWindow ("fractal");

#elif HAVE_SDL

#if DEBUG
  printf ("Initing SDL window\n");
  fflush (stdout);
#endif
  SDL_Init (SDL_INIT_VIDEO);
  window = SDL_CreateWindow ("fractal",
                             SDL_WINDOWPOS_CENTERED, SDL_WINDOWPOS_CENTERED,
                             window_width, window_height,
                             SDL_WINDOW_RESIZABLE | SDL_WINDOW_OPENGL);
  if (!window)
    {
      printf ("ERROR! unable to create the window: %s\n", SDL_GetError ());
      return 1;
    }
  SDL_GL_SetAttribute (SDL_GL_CONTEXT_MAJOR_VERSION, 2);
  if (!SDL_GL_CreateContext (window))
    {
      printf ("ERROR! SDL_GL_CreateContext: %s\n", SDL_GetError ());
      return 1;
    }

#elif HAVE_GLFW

#if DEBUG
  printf ("Initing GLFW window\n");
  fflush (stdout);
#endif
  if (!glfwInit ())
    {
      printf ("ERROR! unable to init GLFW\n");
      return 1;
    }
  window
    = glfwCreateWindow (window_width, window_height, "fractal", NULL, NULL);
  if (!window)
    {
      printf ("ERROR! unable to open the window\n");
      glfwTerminate ();
      return 1;
    }
  glfwMakeContextCurrent (window);

#endif

  // Initing GTK+
#if DEBUG
  printf ("Initing GTK+\n");
  fflush (stdout);
#endif
#if !GTK4
  gtk_init (&argn, &argc);
#else
  gtk_init ();
#endif

  // Creating the main GTK+ window
#if DEBUG
  printf ("Creating simulator dialog\n");
  fflush (stdout);
#endif
  dialog_simulator_create ();

#if !HAVE_GTKGLAREA
  // Initing drawing data
#if DEBUG
  printf ("Initing drawing data\n");
  fflush (stdout);
#endif
  if (!graphic_init (graphic, "logo.png"))
    return 1;
#endif

  // Opening input file
  if (argn == 2 && !fractal_input (argc[1]))
    return 1;

  // Updating view
  if (argn == 2)
    fractal ();
  dialog_simulator_update ();

  // Main loop
#if DEBUG
  printf ("Main loop\n");
  fflush (stdout);
#endif
  main_loop ();

  // Freeing memory
#if HAVE_GLFW
  glfwDestroyWindow (window);
  glfwTerminate ();
#endif
  graphic_destroy (graphic);

  return 0;
}

max_stars_repo_name: jburguete/fractal

max_issues_repo_name: jburguete/fractal

max_forks_repo_name: jburguete/fractal


text:
//This takes univariate X, windows X with W, applies the 1D FFT to each frame,
//takes power, transforms with T mat to B cfs, and then compresses power with log.
//Then the cepstral coeffs (CCs) are obtained by the 1D DCT-II of each frame,
//followed by a lifter (weighting in the cepstral domain).

//This is equivalent to :                                                                                spectrogram -> get_ccs.
//which is equivalent to:                                                                                spectrogram -> dct -> lifter.
//which is equivalent to:                                             stft -> apply_spectrogram_T_mat -> pow_compress -> dct -> lifter.
//which is equivalent to:                     window_univar -> fft_squared -> apply_spectrogram_T_mat -> pow_compress -> dct -> lifter.
//which is equivalent to:             window_univar -> fft_hc -> hc_square -> apply_spectrogram_T_mat -> pow_compress -> dct -> lifter.
//which is equivalent to: frame_univar -> apply_win -> fft_hc -> hc_square -> apply_spectrogram_T_mat -> pow_compress -> dct -> lifter.

#include <stdio.h>
#include <float.h>
#include <math.h>
#include <cblas.h>
#include <fftw3.h>

#ifndef M_PIf
   #define M_PIf 3.141592653589793238462643383279502884f
#endif

#ifndef M_PI
   #define M_PI 3.141592653589793238462643383279502884
#endif

#ifdef __cplusplus
namespace ov {
extern "C" {
#endif

int mfccs_s (float *Y, const int iscolmajor, const int R, const int C, const float *X, const int N, const float *W, const int L, const float *H, const int B, const int dim, const int c0, const float stp, const int mn0, const int nfft, const float preg, const int ndct, const float Q, const int K);
int mfccs_d (double *Y, const int iscolmajor, const int R, const int C, const double *X, const int N, const double *W, const int L, const double *H, const int B, const int dim, const int c0, const double stp, const int mn0, const int nfft, const double preg, const int ndct, const double Q, const int K);

//int mfccs_s (float *Y, sigi Yi, const float *X, const sigi Xi, const float *W, const sigi Wi, const float *H, const sigi H, const int dim, const int c0, const float stp, const float preg, const int ndct, const float Q)
int mfccs_s (float *Y, const int iscolmajor, const int R, const int C, const float *X, const int N, const float *W, const int L, const float *H, const int B, const int dim, const int c0, const float stp, const int mn0, const int nfft, const float preg, const int ndct, const float Q, const int K)
{
    const float z = 0.0f, o = 1.0f;
    const float Q2 = 0.5f*Q, PQ = M_PIf/Q;
    //const int normalize_H = 0;
    const int F = nfft/2 + 1;   //num non-negative FFT freqs
    const int F2 = F-1+nfft%2;  //for square loop below
    const int Lpre = L/2;       //nsamps before center samp
    int ss = c0 - Lpre;         //start samp of current frame
    int r, c, f, b, k;
    float *X1, *Y1, *X1d, *Y1d; //, *freqs, *cfs, *H;
    fftwf_plan fft_plan, dct_plan;
    //struct timespec tic, toc;

    //Checks
    if (N<1) { fprintf(stderr,"error in mfccs_s: N (length X) must be positive\n"); return 1; }
    if (R<1) { fprintf(stderr,"error in mfccs_s: R (nrows Y) must be positive\n"); return 1; }
    if (C<1) { fprintf(stderr,"error in mfccs_s: C (ncols Y) must be positive\n"); return 1; }
    if (L<1) { fprintf(stderr,"error in mfccs_s: L (winlength) must be positive\n"); return 1; }
    if (c0<0) { fprintf(stderr,"error in mfccs_s: c0 (center samp of 1st frame) must be nonnegative\n"); return 1; }
    if (c0>N-1) { fprintf(stderr,"error in mfccs_s: c0 (center samp of 1st frame) must be < N (length X)\n"); return 1; }
    if (L>=N) { fprintf(stderr,"error in mfccs_s: L (winlength) must be < N (length X)\n"); return 1; }
    if (stp<=0.0f) { fprintf(stderr,"error in mfccs_s: stp (step size) must be positive\n"); return 1; }
    if (nfft<L) { fprintf(stderr,"error in mfccs_s: nfft must be >= L (winlength)\n"); return 1; }
    if (preg!=preg || preg<0.0f) { fprintf(stderr,"error in mfccs_s: preg must be nonnegative\n"); return 1; }
    if (Q<0.0f) { fprintf(stderr,"error in mfccs_s: Q must be positive (or 0 to skip lifter)\n"); return 1; }
    if (K<1) { fprintf(stderr,"error in mfccs_s: K must be positive\n"); return 1; }
    if (K>ndct) { fprintf(stderr,"error in mfccs_s: K must be <= ndct>\n"); return 1; }

    //Initialize fftwf for FFT
    X1 = fftwf_alloc_real((size_t)nfft);
    Y1 = fftwf_alloc_real((size_t)nfft);
    fft_plan = fftwf_plan_r2r_1d(nfft,X1,Y1,FFTW_R2HC,FFTW_ESTIMATE);
    if (!fft_plan) { fprintf(stderr,"error in mfccs_s: problem creating fftw plan for fft\n"); return 1; }
    cblas_scopy(nfft,&z,0,&X1[0],1); //zero-pad

    //Initialize fftwf for DCT
    X1d = fftwf_alloc_real((size_t)ndct);
    Y1d = fftwf_alloc_real((size_t)ndct);
    dct_plan = fftwf_plan_r2r_1d(ndct,X1d,Y1d,FFTW_REDFT10,FFTW_ESTIMATE);
    if (!dct_plan) { fprintf(stderr,"error in mfccs_s: problem creating fftw plan for dct\n"); return 1; }
    cblas_scopy(ndct,&z,0,&X1d[0],1); //zero-pad

    //Get frequencies and H
    //if (get_stft_freqs_s(&freqs[0],F,fs,nfft)) { fprintf(stderr,"error in mfccs_s: problem getting STFT freqs\n"); return 1; }
    //if (get_cfs_T_s(&cfs[0],B,fs,'mel')) { fprintf(stderr,"error in mfccs_s: problem getting Mel cfs\n"); return 1; }
    //if (get_spectrogram_T_mat_s(&H[0],iscolmajor,&freqs[0],F,&cfs[0],B,normalize_H))
    //{ fprintf(stderr,"error in mfccs_s: problem getting STFT freq -> Mel cfs transform matrix\n"); return 1; }

    //clock_gettime(CLOCK_REALTIME,&tic);

    //Initialize lifter
    //float *lift; int n;
    //if (!(lift=(float *)malloc((size_t)K*sizeof(float)))) { fprintf(stderr,"error in get_ccs_s: problem with malloc for lifter. "); perror("malloc"); return 1; }
    //for (k=0; k<K; k++) { lift[k] = fmaf(Q2,sinf(k*PQ),1.0f); }

    if (dim==0)
    {
        if (iscolmajor)
        {
            for (c=0; c<C; c++)
            {
                //Window
                ss = (int)(roundf(c*stp)) + c0 - Lpre;
                if (ss<0)
                {
                    cblas_scopy(-ss,&z,0,&X1[0],1);
                    cblas_ssbmv(CblasColMajor,CblasUpper,L+ss,0,1.0f,&W[-ss],1,&X[0],1,0.0f,&X1[-ss],1);
                }
                else if (ss+L<=N)
                {
                    cblas_ssbmv(CblasColMajor,CblasUpper,L,0,1.0f,&W[0],1,&X[ss],1,0.0f,&X1[0],1);
                }
                else if (ss<N)
                {
                    cblas_ssbmv(CblasColMajor,CblasUpper,N-ss,0,1.0f,&W[0],1,&X[ss],1,0.0f,&X1[0],1);
                    cblas_scopy(L-N+ss,&z,0,&X1[N-ss],1);
                }
                else
                {
                    cblas_scopy(L,&z,0,&X1[0],1);
                }

                //Subtract mean
                if (mn0) { cblas_saxpy(L,-cblas_sdot(L,&X1[0],1,&o,0)/L,&o,0,&X1[0],1); }

                //FFT
                fftwf_execute(fft_plan);

                //Power
                for (f=0; f<nfft; f++) { Y1[f] *= Y1[f]; }
                for (f=1; f<F2; f++) { Y1[f] += Y1[nfft-f]; }

                //Transform freq scale
                cblas_sgemv(CblasColMajor,CblasNoTrans,B,F,1.0f,&H[0],B,&Y1[0],1,0.0f,&X1d[0],1);

                //Power compress
                for (b=0; b<B; b++) { X1d[b] = logf(X1d[b]+preg); }

                //DCT
                fftwf_execute(dct_plan);
                cblas_scopy(K,&Y1d[0],1,&Y[c*K],1);
                //cblas_ssbmv(CblasColMajor,CblasUpper,K,0,1.0f,&lift[0],1,&Y1d[0],1,0.0f,&Y[c*K],1);
                //for (k=0; k<K; k++) { Y[k+c*K] = lift[k]*Y1d[k]; }
            }
            //Lifter
            if (Q>FLT_EPSILON) { for (k=0; k<K; k++) { cblas_sscal(C,fmaf(Q2,sinf(PQ*k),1.0f),&Y[k],K); } }
        }
        else
        {
            for (c=0; c<C; c++)
            {
                //Window
                ss = (int)(roundf(c*stp)) + c0 - Lpre;
                if (ss<0)
                {
                    cblas_scopy(-ss,&z,0,&X1[0],1);
                    cblas_ssbmv(CblasRowMajor,CblasUpper,L+ss,0,1.0f,&W[-ss],1,&X[0],1,0.0f,&X1[-ss],1);
                }
                else if (ss+L<=N)
                {
                    cblas_ssbmv(CblasRowMajor,CblasUpper,L,0,1.0f,&W[0],1,&X[ss],1,0.0f,&X1[0],1);
                }
                else if (ss<N)
                {
                    cblas_ssbmv(CblasRowMajor,CblasUpper,N-ss,0,1.0f,&W[0],1,&X[ss],1,0.0f,&X1[0],1);
                    cblas_scopy(L-N+ss,&z,0,&X1[N-ss],1);
                }
                else
                {
                    cblas_scopy(L,&z,0,&X1[0],1);
                }

                //Subtract mean
                if (mn0) { cblas_saxpy(L,-cblas_sdot(L,&X1[0],1,&o,0)/L,&o,0,&X1[0],1); }

                //FFT
                fftwf_execute(fft_plan);

                //Power
                for (f=0; f<nfft; f++) { Y1[f] *= Y1[f]; }
                for (f=1; f<F2; f++) { Y1[f] += Y1[nfft-f]; }

                //Transform freq scale
                cblas_sgemv(CblasRowMajor,CblasNoTrans,B,F,1.0f,&H[0],F,&Y1[0],1,0.0f,&X1d[0],1);

                //Power compress
                for (b=0; b<B; b++) { X1d[b] = logf(X1d[b]+preg); }

                //DCT
                fftwf_execute(dct_plan);
                cblas_scopy(K,&Y1d[0],1,&Y[c],C);
            }
            //Lifter
            if (Q>FLT_EPSILON) { for (k=0; k<K; k++) { cblas_sscal(C,fmaf(Q2,sinf(PQ*k),1.0f),&Y[k*C],1); } }
        }
    }
    else if (dim==1)
    {
        if (iscolmajor)
        {
            for (r=0; r<R; r++)
            {
                //Window
                ss = (int)(roundf(r*stp)) + c0 - Lpre;
                if (ss<0)
                {
                    cblas_scopy(-ss,&z,0,&X1[0],1);
                    cblas_ssbmv(CblasColMajor,CblasUpper,L+ss,0,1.0f,&W[-ss],1,&X[0],1,0.0f,&X1[-ss],1);
                }
                else if (ss+L<=N)
                {
                    cblas_ssbmv(CblasColMajor,CblasUpper,L,0,1.0f,&W[0],1,&X[ss],1,0.0f,&X1[0],1);
                }
                else if (ss<N)
                {
                    cblas_ssbmv(CblasColMajor,CblasUpper,N-ss,0,1.0f,&W[0],1,&X[ss],1,0.0f,&X1[0],1);
                    cblas_scopy(L-N+ss,&z,0,&X1[N-ss],1);
                }
                else
                {
                    cblas_scopy(L,&z,0,&X1[0],1);
                }

                //Subtract mean
                if (mn0) { cblas_saxpy(L,-cblas_sdot(L,&X1[0],1,&o,0)/L,&o,0,&X1[0],1); }

                //FFT
                fftwf_execute(fft_plan);

                //Power
                for (f=0; f<nfft; f++) { Y1[f] *= Y1[f]; }
                for (f=1; f<F2; f++) { Y1[f] += Y1[nfft-f]; }

                //Transform freq scale
                cblas_sgemv(CblasColMajor,CblasNoTrans,B,F,1.0f,&H[0],B,&Y1[0],1,0.0f,&X1d[0],1);

                //Power compress
                for (b=0; b<B; b++) { X1d[b] = logf(X1d[b]+preg); }

                //DCT
                fftwf_execute(dct_plan);
                cblas_scopy(K,&Y1d[0],1,&Y[r],R);
            }
            //Lifter
            if (Q>FLT_EPSILON) { for (k=0; k<K; k++) { cblas_sscal(R,fmaf(Q2,sinf(PQ*k),1.0f),&Y[k*R],1); } }
        }
        else
        {
            for (r=0; r<R; r++)
            {
                //Window
                ss = (int)(roundf(r*stp)) + c0 - Lpre;
                if (ss<0)
                {
                    cblas_scopy(-ss,&z,0,&X1[0],1);
                    cblas_ssbmv(CblasRowMajor,CblasUpper,L+ss,0,1.0f,&W[-ss],1,&X[0],1,0.0f,&X1[-ss],1);
                }
                else if (ss+L<=N)
                {
                    cblas_ssbmv(CblasRowMajor,CblasUpper,L,0,1.0f,&W[0],1,&X[ss],1,0.0f,&X1[0],1);
                }
                else if (ss<N)
                {
                    cblas_ssbmv(CblasRowMajor,CblasUpper,N-ss,0,1.0f,&W[0],1,&X[ss],1,0.0f,&X1[0],1);
                    cblas_scopy(L-N+ss,&z,0,&X1[N-ss],1);
                }
                else
                {
                    cblas_scopy(L,&z,0,&X1[0],1);
                }

                //Subtract mean
                if (mn0) { cblas_saxpy(L,-cblas_sdot(L,&X1[0],1,&o,0)/L,&o,0,&X1[0],1); }

                //FFT
                fftwf_execute(fft_plan);

                //Power
                for (f=0; f<nfft; f++) { Y1[f] *= Y1[f]; }
                for (f=1; f<F2; f++) { Y1[f] += Y1[nfft-f]; }

                //Transform freq scale
                cblas_sgemv(CblasRowMajor,CblasNoTrans,B,F,1.0f,&H[0],F,&Y1[0],1,0.0f,&X1d[0],1);

                //Power compress
                for (b=0; b<B; b++) { X1d[b] = logf(X1d[b]+preg); }

                //DCT
                fftwf_execute(dct_plan);
                cblas_scopy(K,&Y1d[0],1,&Y[r*K],1);
            }
            //Lifter
            if (Q>FLT_EPSILON) { for (k=0; k<K; k++) { cblas_sscal(R,fmaf(Q2,sinf(PQ*k),1.0f),&Y[k],K); } }
        }
    }
    else
    {
        fprintf(stderr,"error in mfccs_s: dim must be 0 or 1.\n"); return 1;
    }
    //clock_gettime(CLOCK_REALTIME,&toc); fprintf(stderr,"elapsed time = %f ms\n",(toc.tv_sec-tic.tv_sec)*1e3+(toc.tv_nsec-tic.tv_nsec)/1e6);

    //Exit
    fftwf_destroy_plan(fft_plan); fftwf_free(X1); fftwf_free(Y1);
    fftwf_destroy_plan(dct_plan); fftwf_free(X1d); fftwf_free(Y1d);
    return 0;
}

int mfccs_d (double *Y, const int iscolmajor, const int R, const int C, const double *X, const int N, const double *W, const int L, const double *H, const int B, const int dim, const int c0, const double stp, const int mn0, const int nfft, const double preg, const int ndct, const double Q, const int K)
{
    const double z = 0.0, o = 1.0;
    const double Q2 = 0.5*Q, PQ = M_PI/Q;
    const int F = nfft/2 + 1;   //num non-negative FFT freqs
    const int F2 = F-1+nfft%2;  //for square loop below
    const int Lpre = L/2;       //nsamps before center samp
    int ss = c0 - Lpre;         //start samp of current frame
    int r, c, f, b, k;
    double *X1, *Y1, *X1d, *Y1d;
    fftw_plan fft_plan, dct_plan;

    //Checks
    if (N<1) { fprintf(stderr,"error in mfccs_d: N (length X) must be positive\n"); return 1; }
    if (R<1) { fprintf(stderr,"error in mfccs_d: R (nrows Y) must be positive\n"); return 1; }
    if (C<1) { fprintf(stderr,"error in mfccs_d: C (ncols Y) must be positive\n"); return 1; }
    if (L<1) { fprintf(stderr,"error in mfccs_d: L (winlength) must be positive\n"); return 1; }
    if (c0<0) { fprintf(stderr,"error in mfccs_d: c0 (center samp of 1st frame) must be nonnegative\n"); return 1; }
    if (c0>N-1) { fprintf(stderr,"error in mfccs_d: c0 (center samp of 1st frame) must be < N (length X)\n"); return 1; }
    if (L>=N) { fprintf(stderr,"error in mfccs_d: L (winlength) must be < N (length X)\n"); return 1; }
    if (stp<=0.0) { fprintf(stderr,"error in mfccs_d: stp (step size) must be positive\n"); return 1; }
    if (nfft<L) { fprintf(stderr,"error in mfccs_d: nfft must be >= L (winlength)\n"); return 1; }
    if (preg!=preg || preg<0.0) { fprintf(stderr,"error in mfccs_d: preg must be nonnegative\n"); return 1; }
    if (Q<0.0) { fprintf(stderr,"error in mfccs_d: Q must be positive (or 0 to skip lifter)\n"); return 1; }
    if (K<1) { fprintf(stderr,"error in mfccs_d: K must be positive\n"); return 1; }
    if (K>ndct) { fprintf(stderr,"error in mfccs_d: K must be <= ndct>\n"); return 1; }

    //Initialize fftw for FFT
    X1 = fftw_alloc_real((size_t)nfft);
    Y1 = fftw_alloc_real((size_t)nfft);
    fft_plan = fftw_plan_r2r_1d(nfft,X1,Y1,FFTW_R2HC,FFTW_ESTIMATE);
    if (!fft_plan) { fprintf(stderr,"error in mfccs_d: problem creating fftw plan for fft\n"); return 1; }
    cblas_dcopy(nfft,&z,0,&X1[0],1); //zero-pad

    //Initialize fftw for DCT
    X1d = fftw_alloc_real((size_t)ndct);
    Y1d = fftw_alloc_real((size_t)ndct);
    dct_plan = fftw_plan_r2r_1d(ndct,X1d,Y1d,FFTW_REDFT10,FFTW_ESTIMATE);
    if (!dct_plan) { fprintf(stderr,"error in mfccs_d: problem creating fftw plan for dct\n"); return 1; }
    cblas_dcopy(ndct,&z,0,&X1d[0],1); //zero-pad

    if (dim==0)
    {
        if (iscolmajor)
        {
            for (c=0; c<C; c++)
            {
                //Window
                ss = (int)(round(c*stp)) + c0 - Lpre;
                if (ss<0)
                {
                    cblas_dcopy(-ss,&z,0,&X1[0],1);
                    cblas_dsbmv(CblasColMajor,CblasUpper,L+ss,0,1.0,&W[-ss],1,&X[0],1,0.0,&X1[-ss],1);
                }
                else if (ss+L<=N)
                {
                    cblas_dsbmv(CblasColMajor,CblasUpper,L,0,1.0,&W[0],1,&X[ss],1,0.0,&X1[0],1);
                }
                else if (ss<N)
                {
                    cblas_dsbmv(CblasColMajor,CblasUpper,N-ss,0,1.0,&W[0],1,&X[ss],1,0.0,&X1[0],1);
                    cblas_dcopy(L-N+ss,&z,0,&X1[N-ss],1);
                }
                else
                {
                    cblas_dcopy(L,&z,0,&X1[0],1);
                }

                //Subtract mean
                if (mn0) { cblas_daxpy(L,-cblas_ddot(L,&X1[0],1,&o,0)/L,&o,0,&X1[0],1); }

                //FFT
                fftw_execute(fft_plan);

                //Power
                for (f=0; f<nfft; f++) { Y1[f] *= Y1[f]; }
                for (f=1; f<F2; f++) { Y1[f] += Y1[nfft-f]; }

                //Transform freq scale
                cblas_dgemv(CblasColMajor,CblasNoTrans,B,F,1.0,&H[0],B,&Y1[0],1,0.0,&X1d[0],1);

                //Power compress
                for (b=0; b<B; b++) { X1d[b] = log(X1d[b]+preg); }

                //DCT
                fftw_execute(dct_plan);
                cblas_dcopy(K,&Y1d[0],1,&Y[c*K],1);
            }
            //Lifter
            if (Q>DBL_EPSILON) { for (k=0; k<K; k++) { cblas_dscal(C,fma(Q2,sin(PQ*k),1.0),&Y[k],K); } }
        }
        else
        {
            for (c=0; c<C; c++)
            {
                //Window
                ss = (int)(round(c*stp)) + c0 - Lpre;
                if (ss<0)
                {
                    cblas_dcopy(-ss,&z,0,&X1[0],1);
                    cblas_dsbmv(CblasRowMajor,CblasUpper,L+ss,0,1.0,&W[-ss],1,&X[0],1,0.0,&X1[-ss],1);
                }
                else if (ss+L<=N)
                {
                    cblas_dsbmv(CblasRowMajor,CblasUpper,L,0,1.0,&W[0],1,&X[ss],1,0.0,&X1[0],1);
                }
                else if (ss<N)
                {
                    cblas_dsbmv(CblasRowMajor,CblasUpper,N-ss,0,1.0,&W[0],1,&X[ss],1,0.0,&X1[0],1);
                    cblas_dcopy(L-N+ss,&z,0,&X1[N-ss],1);
                }
                else
                {
                    cblas_dcopy(L,&z,0,&X1[0],1);
                }

                //Subtract mean
                if (mn0) { cblas_daxpy(L,-cblas_ddot(L,&X1[0],1,&o,0)/L,&o,0,&X1[0],1); }

                //FFT
                fftw_execute(fft_plan);

                //Power
                for (f=0; f<nfft; f++) { Y1[f] *= Y1[f]; }
                for (f=1; f<F2; f++) { Y1[f] += Y1[nfft-f]; }

                //Transform freq scale
                cblas_dgemv(CblasRowMajor,CblasNoTrans,B,F,1.0,&H[0],F,&Y1[0],1,0.0,&X1d[0],1);

                //Power compress
                for (b=0; b<B; b++) { X1d[b] = log(X1d[b]+preg); }

                //DCT
                fftw_execute(dct_plan);
                cblas_dcopy(K,&Y1d[0],1,&Y[c],C);
            }
            //Lifter
            if (Q>DBL_EPSILON) { for (k=0; k<K; k++) { cblas_dscal(C,fma(Q2,sin(PQ*k),1.0),&Y[k*C],1); } }
        }
    }
    else if (dim==1)
    {
        if (iscolmajor)
        {
            for (r=0; r<R; r++)
            {
                //Window
                ss = (int)(round(r*stp)) + c0 - Lpre;
                if (ss<0)
                {
                    cblas_dcopy(-ss,&z,0,&X1[0],1);
                    cblas_dsbmv(CblasColMajor,CblasUpper,L+ss,0,1.0,&W[-ss],1,&X[0],1,0.0,&X1[-ss],1);
                }
                else if (ss+L<=N)
                {
                    cblas_dsbmv(CblasColMajor,CblasUpper,L,0,1.0,&W[0],1,&X[ss],1,0.0,&X1[0],1);
                }
                else if (ss<N)
                {
                    cblas_dsbmv(CblasColMajor,CblasUpper,N-ss,0,1.0,&W[0],1,&X[ss],1,0.0,&X1[0],1);
                    cblas_dcopy(L-N+ss,&z,0,&X1[N-ss],1);
                }
                else
                {
                    cblas_dcopy(L,&z,0,&X1[0],1);
                }

                //Subtract mean
                if (mn0) { cblas_daxpy(L,-cblas_ddot(L,&X1[0],1,&o,0)/L,&o,0,&X1[0],1); }

                //FFT
                fftw_execute(fft_plan);

                //Power
                for (f=0; f<nfft; f++) { Y1[f] *= Y1[f]; }
                for (f=1; f<F2; f++) { Y1[f] += Y1[nfft-f]; }

                //Transform freq scale
                cblas_dgemv(CblasColMajor,CblasNoTrans,B,F,1.0,&H[0],B,&Y1[0],1,0.0,&X1d[0],1);

                //Power compress
                for (b=0; b<B; b++) { X1d[b] = log(X1d[b]+preg); }

                //DCT
                fftw_execute(dct_plan);
                cblas_dcopy(K,&Y1d[0],1,&Y[r],R);
            }
            //Lifter
            if (Q>DBL_EPSILON) { for (k=0; k<K; k++) { cblas_dscal(R,fma(Q2,sin(PQ*k),1.0),&Y[k*R],1); } }
        }
        else
        {
            for (r=0; r<R; r++)
            {
                //Window
                ss = (int)(round(r*stp)) + c0 - Lpre;
                if (ss<0)
                {
                    cblas_dcopy(-ss,&z,0,&X1[0],1);
                    cblas_dsbmv(CblasRowMajor,CblasUpper,L+ss,0,1.0,&W[-ss],1,&X[0],1,0.0,&X1[-ss],1);
                }
                else if (ss+L<=N)
                {
                    cblas_dsbmv(CblasRowMajor,CblasUpper,L,0,1.0,&W[0],1,&X[ss],1,0.0,&X1[0],1);
                }
                else if (ss<N)
                {
                    cblas_dsbmv(CblasRowMajor,CblasUpper,N-ss,0,1.0,&W[0],1,&X[ss],1,0.0,&X1[0],1);
                    cblas_dcopy(L-N+ss,&z,0,&X1[N-ss],1);
                }
                else
                {
                    cblas_dcopy(L,&z,0,&X1[0],1);
                }

                //Subtract mean
                if (mn0) { cblas_daxpy(L,-cblas_ddot(L,&X1[0],1,&o,0)/L,&o,0,&X1[0],1); }

                //FFT
                fftw_execute(fft_plan);

                //Power
                for (f=0; f<nfft; f++) { Y1[f] *= Y1[f]; }
                for (f=1; f<F2; f++) { Y1[f] += Y1[nfft-f]; }

                //Transform freq scale
                cblas_dgemv(CblasRowMajor,CblasNoTrans,B,F,1.0,&H[0],F,&Y1[0],1,0.0,&X1d[0],1);

                //Power compress
                for (b=0; b<B; b++) { X1d[b] = log(X1d[b]+preg); }

                //DCT
                fftw_execute(dct_plan);
                cblas_dcopy(K,&Y1d[0],1,&Y[r*K],1);
            }
            //Lifter
            if (Q>DBL_EPSILON) { for (k=0; k<K; k++) { cblas_dscal(R,fma(Q2,sin(PQ*k),1.0),&Y[k],K); } }
        }
    }
    else
    {
        fprintf(stderr,"error in mfccs_d: dim must be 0 or 1.\n"); return 1;
    }

    //Exit
    fftw_destroy_plan(fft_plan); fftw_free(X1); fftw_free(Y1);
    fftw_destroy_plan(dct_plan); fftw_free(X1d); fftw_free(Y1d);
    return 0;
}

#ifdef __cplusplus
}
}
#endif


text:
/**
 *
 * @file core_dpotrf.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Hatem Ltaief
 * @author Mathieu Faverge
 * @author Jakub Kurzak
 * @date 2010-11-15
 * @generated d Tue Jan  7 11:44:45 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 * @ingroup CORE_double
 *
 *  CORE_dpotrf - Computes the Cholesky factorization of a symmetric positive definite
 *  (or Hermitian positive definite in the complex case) matrix A.
 *  The factorization has the form
 *
 *    \f[ A = \{_{L\times L^H, if uplo = PlasmaLower}^{U^H\times U, if uplo = PlasmaUpper} \f]
 *
 *  where U is an upper triangular matrix and L is a lower triangular matrix.
 *
 *******************************************************************************
 *
 * @param[in] uplo
 *          = PlasmaUpper: Upper triangle of A is stored;
 *          = PlasmaLower: Lower triangle of A is stored.
 *
 * @param[in] N
 *          The order of the matrix A. N >= 0.
 *

 * @param[in,out] A
 *          On entry, the symmetric positive definite (or Hermitian) matrix A.
 *          If uplo = PlasmaUpper, the leading N-by-N upper triangular part of A
 *          contains the upper triangular part of the matrix A, and the strictly
 *          lower triangular part of A is not referenced.
 *          If UPLO = 'L', the leading N-by-N lower triangular part of A
 *          contains the lower triangular part of the matrix A, and the strictly
 *          upper triangular part of A is not referenced.
 *          On exit, if return value = 0, the factor U or L from the Cholesky
 *          factorization A = U**T*U or A = L*L**T.
 *
 * @param[in] LDA
 *          The leading dimension of the array A. LDA >= max(1,N).
 *
 * @param[out] info
 *          - 0 on successful exit
 *          - <0 if -i, the i-th argument had an illegal value
 *          - >0 if i, the leading minor of order i of A is not positive
 *               definite, so the factorization could not be completed, and the
 *               solution has not been computed.
 *
 ******************************************************************************/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_dpotrf = PCORE_dpotrf
#define CORE_dpotrf PCORE_dpotrf
#endif
void CORE_dpotrf(PLASMA_enum uplo, int N, double *A, int LDA, int *info)
{
    *info = LAPACKE_dpotrf_work(
        LAPACK_COL_MAJOR,
        lapack_const(uplo),
        N, A, LDA );
}


text:
/**
 * @file batchf_zgemv.c
 *
 * Part of API test for Batched BLAS routines.
 *
 * @author  Samuel  D. Relton
 * @author  Pedro   V. Lara
 * @author  Mawussi Zounon
 * @date 2016-06-01
 *
 * @precisions normal z -> c d s
 *
 **/

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

#define COMPLEX

void batchf_zgemv(
	const enum BBLAS_TRANS trans,
	const int m, const int n,
	const BBLAS_Complex64_t alpha,
	const BBLAS_Complex64_t **arrayA, const int lda,
	const BBLAS_Complex64_t **arrayx, const int incx,
	const BBLAS_Complex64_t beta,
	BBLAS_Complex64_t **arrayy, const int incy,
	const int batch_count, int info)
{
	/* Local variables */
	int first_index = 0;
	int batch_iter = 0;
	char func_name[15] = "batchf_zgemv";

	/* Check input arguments */
	if (batch_count < 0)
	{
		xerbla_batch(func_name, BBLAS_ERR_BATCH_COUNT, -1);
	}

	if ((trans != BblasTrans) &&
		(trans != BblasNoTrans) &&
		(trans != BblasConjTrans))
	{
		xerbla_batch(func_name, BBLAS_ERR_TRANS, first_index);
		info = BBLAS_ERR_TRANS;
	}

	if (m < 0)
	{
		xerbla_batch(func_name, BBLAS_ERR_M, first_index);
		info = BBLAS_ERR_M;
	}

	if (n < 0)
	{
		xerbla_batch(func_name, BBLAS_ERR_N, first_index);
		info = BBLAS_ERR_N;
	}

	/* Column major */
	if ((lda < 1) && (lda < m))
	{
		xerbla_batch(func_name, BBLAS_ERR_LDA, first_index);
		info = BBLAS_ERR_LDA;
	}

	if (incx < 1)
	{
		xerbla_batch(func_name, BBLAS_ERR_INCX, first_index);
		info = BBLAS_ERR_INCX;
	}
	if (incy < 1)
	{
		xerbla_batch(func_name, BBLAS_ERR_INCY, first_index);
		info = BBLAS_ERR_INCY;
	}

	/* Call CBLAS */
	for (batch_iter = 0; batch_iter < batch_count; batch_iter++)
	{
		cblas_zgemv(
			BblasColMajor,
			trans,
			m, n,
			CBLAS_SADDR( alpha ),
			arrayA[batch_iter], lda,
			arrayx[batch_iter], incx,
			CBLAS_SADDR( beta ),
			arrayy[batch_iter], incy);
	} /* End fixed size for loop */
	/* Successful */
	info = BBLAS_SUCCESS;
}
#undef COMPLEX

max_stars_repo_name: mawussi/BBLAS-group