text:
#ifndef __PROJ_GAUSS_MIXTURES_H__
#define __PROJ_GAUSS_MIXTURES_H__

/* include */
#include <stdbool.h>
#include <float.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_linalg.h>

/* variable declarations and definitions */

double halflogtwopi; /* constant used in calculation */

int nthreads;
int tid;

struct gaussian{
  double alpha;
  gsl_vector *mm;
  gsl_matrix *VV;
};

struct datapoint{
  gsl_vector *ww;
  gsl_matrix *SS;
  gsl_matrix *RR;
  double logweight;
};

struct gaussian * newgaussians, * startnewgaussians;
gsl_matrix * qij;
gsl_permutation * p;
gsl_vector * wminusRm, * TinvwminusRm;
gsl_matrix * Tij,* Tij_inv, * VRT, * VRTTinv, * Rtrans, * I;
struct modelbs{
  gsl_vector *bbij;
  gsl_matrix *BBij;
};

struct modelbs * bs;

FILE *logfile,*convlogfile;

gsl_rng * randgen; /* global random number generator */

/* function declarations */

static inline bool bovy_isfin(double x){ /* true if x is finite */
  return (x > DBL_MAX || x < -DBL_MAX) ? false : true;}
void minmax(gsl_matrix * q, int row, bool isrow, double * min, double * max);
double logsum(gsl_matrix * q, int row, bool isrow);
double normalize_row(gsl_matrix * q, int row,bool isrow,bool noweight, double weight);
void bovy_randvec(gsl_vector * eps, int d, double length); /* returns random vector */
double bovy_det(gsl_matrix * A);/* determinant of matrix A */
void calc_splitnmerge(struct datapoint * data,int N,struct gaussian * gaussians, int K, gsl_matrix * qij, int * snmhierarchy);
void splitnmergegauss(struct gaussian * gaussians,int K, gsl_matrix * qij, int j, int k, int l);
void proj_EM_step(struct datapoint * data, int N, struct gaussian * gaussians, int K,bool * fixamp, bool * fixmean, bool * fixcovar, double * avgloglikedata, bool likeonly, double w,bool noproj, bool diagerrs, bool noweight);
void proj_EM(struct datapoint * data, int N, struct gaussian * gaussians, int K,bool * fixamp, bool * fixmean, bool * fixcovar, double * avgloglikedata, double tol,long long int maxiter, bool likeonly, double w,bool keeplog, FILE *logfile,FILE *tmplogfile, bool noproj, bool diagerrs, bool noweight);
void proj_gauss_mixtures(struct datapoint * data, int N, struct gaussian * gaussians, int K,bool * fixamp, bool * fixmean, bool * fixcovar, double * avgloglikedata, double tol,long long int maxiter, bool likeonly, double w, int splitnmerge, bool keeplog, FILE *logfile,FILE *convlogfile, bool noproj, bool diagerrs, bool noweight);
void calc_qstarij(double * qstarij, gsl_matrix * qij, int partial_indx[3]);

#endif /* proj_gauss_mixtures.h */

max_stars_repo_name: surbut/mashr

max_issues_repo_name: surbut/mashr

max_forks_repo_name: surbut/mashr


text:
/**
  * Copyright 2016 José Manuel Abuín Mosquera <josemanuel.abuin@usc.es>
  *
  * This file is part of Matrix Market Suite.
  *
  * Matrix Market Suite is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
  * the Free Software Foundation, either version 3 of the License, or
  * (at your option) any later version.
  *
  * Matrix Market Suite is distributed in the hope that it will be useful,
  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  * GNU General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with Matrix Market Suite. If not, see <http://www.gnu.org/licenses/>.
  */

#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>

#include <cblas.h>

#include "DMxV.h"
#include "basic.h"

void usageDMxV(){

	fprintf(stderr, "\n");
	fprintf(stderr, "Usage: MM-Suite DMxV [options] <input-matrix> <input-vector>\n");
	fprintf(stderr, "\nInput/output options:\n\n");
	fprintf(stderr, "       -o STR        Output file name. Default: stdout\n");
	fprintf(stderr, "       -r            Input format is row per line. Default: False\n");
	fprintf(stderr, "\nParameters options:\n\n");
	fprintf(stderr, "       -a DOUBLE     Alpha. Default: 1.0\n");
	fprintf(stderr, "       -b DOUBLE     Beta. Default: 0.0\n");
	fprintf(stderr, "\nPerformance options:\n\n");
	fprintf(stderr, "       -t INT        Number of threads to use in OpenBLAS. Default: 1\n");
	fprintf(stderr, "\n");

}

int DMxV(int argc, char *argv[]) {

	int 			ret_code = 1;
	int 			option;

	unsigned long 		*II;
	unsigned long 		*J;
	double 			*values;

	unsigned long 		M;
	unsigned long 		N;
	unsigned long long 	nz;

	double 			*vectorValues;
	unsigned long 		M_Vector;
	unsigned long 		N_Vector;
	unsigned long long 	nz_vector;

	char			*outputFileName = NULL;

	char			*inputMatrixFile = NULL;
	char			*inputVectorFile = NULL;
	char			*outputVectorFile = NULL;

	int			inputFormatRow = 0;
	int			basicOps = 0;
	int			numThreads = 1;

	double			alpha = 1.0;
	double			beta = 0.0;

	while ((option = getopt(argc, argv,"ero:b:a:t:")) >= 0) {
		switch (option) {
			case 'o' :
				//free(outputFileName);

				outputFileName = (char *) malloc(sizeof(char)*strlen(optarg)+1);
				strcpy(outputFileName,optarg);

				break;

			case 'r':
				inputFormatRow = 1;
				break;

			case 'e':
				basicOps = 1;
				break;

			case 'b':
				beta = atof(optarg);
				break;

			case 'a':
				alpha = atof(optarg);
				break;

			case 't':
				numThreads = atoi(optarg);
				break;

			default: break;
		}

	}

	if ((optind + 3 != argc) && (optind + 2 != argc)) {
		usageDMxV();
		return 0;
	}

	openblas_set_num_threads(numThreads);

	if(optind + 3 == argc) { //We have an output vector

		outputVectorFile = (char *)malloc(sizeof(char)*strlen(argv[optind+2])+1);
		strcpy(outputVectorFile,argv[optind+2]);
	}

	if(outputFileName == NULL) {
		outputFileName = (char *) malloc(sizeof(char)*7);
		sprintf(outputFileName,"stdout");
	}

	inputMatrixFile = (char *)malloc(sizeof(char)*strlen(argv[optind])+1);

	if(inputMatrixFile == NULL) {
		fprintf(stderr, "[%s] Error reserving memory for input matrix file name\n",__func__);
		return 0;
	}

	inputVectorFile = (char *)malloc(sizeof(char)*strlen(argv[optind+1])+1);

	if(inputVectorFile == NULL) {
		fprintf(stderr, "[%s] Error reserving memory for input vector file name\n",__func__);
		return 0;
	}

	strcpy(inputMatrixFile,argv[optind]);
	strcpy(inputVectorFile,argv[optind+1]);

	//Read matrix
	if(inputFormatRow){

		if(!readDenseCoordinateMatrixRowLine(inputMatrixFile,&II,&J,&values,&M,&N,&nz)){
			usageDMxV();
			fprintf(stderr, "[%s] Can not read Matrix\n",__func__);
			return 0;
		}

	}
	else {
		if(!readDenseCoordinateMatrix(inputMatrixFile,&II,&J,&values,&M,&N,&nz)){
			usageDMxV();
			fprintf(stderr, "[%s] Can not read Matrix\n",__func__);
			return 0;
		}
	}

	//Read vector
	if(!readDenseVector(inputVectorFile, &vectorValues,&M_Vector,&N_Vector,&nz_vector)){
		usageDMxV();
		fprintf(stderr, "[%s] Can not read Vector\n",__func__);
		return 0;
	}

	/*
	void cblas_dgemv(const enum CBLAS_ORDER order,
                 const enum CBLAS_TRANSPOSE TransA, 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);
                 */

        double *result=(double *) malloc(nz_vector * sizeof(double));

        //Read output vector if any
	if(outputVectorFile != NULL) {
		if(!readDenseVector(outputVectorFile, &result,&M_Vector,&N_Vector,&nz_vector)){
			usageDMxV();
			fprintf(stderr, "[%s] Can not read Vector %s\n",__func__, outputVectorFile);
			return 0;
		}
	}

	//cblas_dgemv(CblasColMajor,CblasNoTrans,M,N,1.0,values,N,vectorValues,1,0.0,result,1);
	double t_real = realtime();
	if(basicOps){
		mms_dgemv(M, N, alpha, values, vectorValues,beta, result);
	}
	else{
		cblas_dgemv(CblasRowMajor,CblasNoTrans,M,N,alpha,values,N,vectorValues,1,beta,result,1);
	}

	fprintf(stderr, "\n[%s] Time spent in cblas_dgemv: %.6f sec\n", __func__, realtime() - t_real);

	writeDenseVector(outputFileName, result,M_Vector,N_Vector,nz_vector);

	return ret_code;
}

ext: c


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

#include <stdio.h>
#include <fpu_control.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_ieee_utils.h>

  /* Handle libc5, where _FPU_SETCW is not available, suggested by
     OKUJI Yoshinori <okuji@gnu.org> and Evgeny Stambulchik
     <fnevgeny@plasma-gate.weizmann.ac.il> */

#ifndef _FPU_SETCW
#include <i386/fpu_control.h>
#define _FPU_SETCW(cw) __setfpucw(cw)
#endif

int
gsl_ieee_set_mode (int precision, int rounding, int exception_mask)
{
  unsigned short mode = 0 ;

  switch (precision)
    {
    case GSL_IEEE_SINGLE_PRECISION:
      mode |= _FPU_SINGLE ;
      break ;
    case GSL_IEEE_DOUBLE_PRECISION:
      mode |= _FPU_DOUBLE ;
      break ;
    case GSL_IEEE_EXTENDED_PRECISION:
      mode |= _FPU_EXTENDED ;
      break ;
    default:
      mode |= _FPU_EXTENDED ;
    }

  switch (rounding)
    {
    case GSL_IEEE_ROUND_TO_NEAREST:
      mode |= _FPU_RC_NEAREST ;
      break ;
    case GSL_IEEE_ROUND_DOWN:
      mode |= _FPU_RC_DOWN ;
      break ;
    case GSL_IEEE_ROUND_UP:
      mode |= _FPU_RC_UP ;
      break ;
    case GSL_IEEE_ROUND_TO_ZERO:
      mode |= _FPU_RC_ZERO ;
      break ;
    default:
      mode |= _FPU_RC_NEAREST ;
    }

  if (exception_mask & GSL_IEEE_MASK_INVALID)
    mode |= _FPU_MASK_IM ;

  if (exception_mask & GSL_IEEE_MASK_DENORMALIZED)
    mode |= _FPU_MASK_DM ;

  if (exception_mask & GSL_IEEE_MASK_DIVISION_BY_ZERO)
    mode |= _FPU_MASK_ZM ;

  if (exception_mask & GSL_IEEE_MASK_OVERFLOW)
    mode |= _FPU_MASK_OM ;

  if (exception_mask & GSL_IEEE_MASK_UNDERFLOW)
    mode |= _FPU_MASK_UM ;

  if (exception_mask & GSL_IEEE_TRAP_INEXACT)
    {
      mode &= ~ _FPU_MASK_PM ;
    }
  else
    {
      mode |= _FPU_MASK_PM ;
    }

  _FPU_SETCW(mode) ;

  return GSL_SUCCESS ;
}


text:
/* AUTORIGHTS
Copyright (C) 2007 Princeton University

This file is part of Ferret Toolkit.

Ferret Toolkit is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software Foundation,
Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_rng.h>

#include <cass.h>
#include <hash_table.h>
#include <arena.h>

#define	MAX_INT	2147483647

#define DEFAULT_PRIME	1017881
//#define DEFAULT_PRIME	 507217
//#define DEFAULT_PRIME	 432251
//#define DEFAULT_PRIME	 300301
//#define DEFAULT_PRIME	 128311
//#define DEFAULT_PRIME	  10007

MemArena	*hasharena;

typedef struct nodeT {
	int				val;
	struct nodeT	*next;
} nodeT;

typedef struct bucketT {
	int		cnt, idx;
	nodeT	*nodes;
} bucketT;

struct hashTable_t {
	int		size;
	int		cnt_item; 	// number of items
	int		cnt_bucket;	// number of non-empty buckets
	int		key_size;	// size of hash key (number of integers)
	hashT	rndVec;
	bucketT	**buckets;
};

void
hashTable_init(void)
{
	hasharena = mkmemarena(NULL, NULL, NULL, 0);
}

void
hashT_print(hashT key, int key_size, const char *tag, FILE *fp)
{
	int		i;

	fprintf(fp, "[%d", key[0]);
	for (i=1; i<key_size; i++)
		fprintf(fp, ",%d", key[i]);
	fprintf(fp, "]%s", tag);
}

nodeT *
nodeT_new(int val)
{
	nodeT	*ret = (nodeT *)memarenamalloc(hasharena, sizeof(nodeT));

	if (ret) {
		ret->val = val;
		ret->next = NULL;
	}

	return ret;
}

bucketT *
bucketT_new(int idx)
{
	bucketT	*ret = (bucketT *)memarenamalloc(hasharena, sizeof(bucketT));

	if (ret) {
		ret->idx = idx;
		ret->cnt = 0;
		ret->nodes = NULL;
	}

	return ret;
}

void
bucketT_insert(bucketT *bucket, int val)
{
	nodeT	*node = nodeT_new(val);
	if (bucket->cnt > 0)
		node->next = bucket->nodes;

	bucket->nodes = node;
	bucket->cnt++;
}

void
bucketT_free(bucketT *bucket)
{
#ifdef NOTDEF  // arena takes care of freeing the bucket
	nodeT	*p, *q;

	if (bucket) {
		p = bucket->nodes;
		while (p != NULL) {
			q = p->next;
			free(p);
			p = q;
		}

		free(bucket);
	}
#endif
}

int
bucketT_dump(bucketT *bucket, CASS_FILE *out)
{
	int		ret;
	nodeT	*p = bucket->nodes;

	ret = cass_write_uint32(&(bucket->cnt), 1, out);
	assert(ret == 1);

	ret = 0;
	while (p != NULL) {
		ret += cass_write_uint32(&(p->val), 1, out);
		p = p->next;
	}

	assert(ret == bucket->cnt);
	return 0;
}

bucketT *
bucketT_load(int idx, CASS_FILE *in)
{
	int			i, ret, cnt, val;
	bucketT		*bucket = bucketT_new(idx);

	ret = cass_read_uint32(&cnt, 1, in);
	assert(ret == 1);

	ret = 0;
	for (i=0; i<cnt; i++) {
		ret += cass_read_uint32(&val, 1, in);
		bucketT_insert(bucket, val);
	}

	assert(ret == cnt);

	return bucket;
}

hashTable_t *
hashTable_new(int key_size)
{
	int				i;
	hashTable_t		*ret = (hashTable_t *)malloc(sizeof(hashTable_t));

	const gsl_rng_type	*T;
	gsl_rng				*r;

	ret->key_size = key_size;
	ret->rndVec = (int *)calloc(key_size, sizeof(int));

	gsl_rng_env_setup();
	T = gsl_rng_default;
	r = gsl_rng_alloc(T);

	for (i=0; i<key_size; i++)
		ret->rndVec[i] = (int)(gsl_rng_uniform(r) * MAX_INT);

	gsl_rng_free(r);

	ret->size = DEFAULT_PRIME;
	ret->buckets = (bucketT **)calloc(ret->size, sizeof(bucketT *));
	for (i=0; i<ret->size; i++)
		ret->buckets[i] = NULL;

	ret->cnt_item = 0;
	ret->cnt_bucket = 0;

	return ret;
}

int
hashTable_free(hashTable_t *hash_table)
{
	int		i;

	if (hash_table != NULL) {
		for (i=0; i<hash_table->size; i++) {
			if (hash_table->buckets[i] != NULL)
				bucketT_free(hash_table->buckets[i]);
		}

		free(hash_table->buckets);
		free(hash_table->rndVec);
		free(hash_table);
	}

	return 0;
}

int
hashTable_dump(hashTable_t *hash_table, CASS_FILE *out)
{
	int			i, ret;

	ret = cass_write_uint32(&(hash_table->key_size), 1, out);
	ret += cass_write_uint32(&(hash_table->cnt_item), 1, out);
	ret += cass_write_uint32(&(hash_table->cnt_bucket), 1, out);
	assert(ret == 3);

	ret = cass_write_uint32(hash_table->rndVec,
		hash_table->key_size, out);
	assert(ret == hash_table->key_size);

	for (i=0; i<hash_table->size; i++)
		if (hash_table->buckets[i] != NULL) {
			ret = cass_write_uint32(&i, 1, out);
			assert(ret == 1);

			ret = bucketT_dump(hash_table->buckets[i], out);
			if (ret != 0)
				return ret;
		}

	return 0;
}

hashTable_t *
hashTable_load(CASS_FILE *in)
{
	int				i, x, ret;
	int				 key_size, cnt_item, cnt_bucket;
	hashTable_t		*hash_table = type_calloc(hashTable_t, 1);

	ret = cass_read_uint32(&key_size, 1, in);
	ret += cass_read_uint32(&cnt_item, 1, in);
	ret += cass_read_uint32(&cnt_bucket, 1, in);
	assert(ret == 3);

	hash_table->size = DEFAULT_PRIME;
	hash_table->key_size = key_size;
	hash_table->cnt_item = cnt_item;
	hash_table->cnt_bucket = cnt_bucket;

	hash_table->rndVec = (int *)calloc(key_size, sizeof(int));
	ret = cass_read_uint32(hash_table->rndVec, key_size, in);
	assert(ret == key_size);

	hash_table->buckets = (bucketT **)calloc(hash_table->size,
			sizeof(bucketT *));
	for (i=0; i<cnt_bucket; i++) {
		ret = cass_read_uint32(&x, 1, in);
		assert(ret == 1);

		hash_table->buckets[x] = bucketT_load(x, in);
		if (hash_table->buckets[x] == NULL)
			return NULL;
	}

	return hash_table;
}

int
hashTable_hash(hashTable_t *table, hashT key)
{
	int		i, k, sum = 0;

	for (i=0; i<table->key_size; i++) {
		//printf("%d\t", key[i]);
		sum += table->rndVec[i] * key[i];
	}

	k = abs(sum % (table->size));
	//printf("%d\n", k);

	return k;
}

void
hashTable_insert(hashTable_t *table, hashT key, int val)
{
	int			idx;
	bucketT		*bucket;

	idx = hashTable_hash(table, key);
	bucket = table->buckets[idx];
	if (bucket == NULL) {
		bucket = bucketT_new(idx);
		table->buckets[idx] = bucket;
		table->cnt_bucket++;
	}
	bucketT_insert(bucket, val);
	table->cnt_item++;
}

bucketT *
hashTable_get_bucket(hashTable_t *table, hashT key)
{
	int		idx = hashTable_hash(table, key);

	return table->buckets[idx];
}

int
hashTable_search(hashTable_t *table, bitmap_t *map, hashT key)
{
	int			i;
	nodeT		*p;
	int			idx = hashTable_hash(table, key);
	bucketT		*bucket = table->buckets[idx];

	if (bucket) {
		p = bucket->nodes;
		for (i=0; i<bucket->cnt; i++) {
			bitmap_insert(map, p->val);
			p = p->next;
		}
	}

	return 0;
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
// --------------------------------------------------------------------------
//                   OpenMS -- Open-Source Mass Spectrometry
// --------------------------------------------------------------------------
// Copyright The OpenMS Team -- Eberhard Karls University Tuebingen,
// ETH Zurich, and Freie Universitaet Berlin 2002-2013.
//
// This software is released under a three-clause BSD license:
//  * Redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//  * 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.
//  * Neither the name of any author or any participating institution
//    may be used to endorse or promote products derived from this software
//    without specific prior written permission.
// For a full list of authors, refer to the file AUTHORS.
// --------------------------------------------------------------------------
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "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 ANY OF THE AUTHORS OR THE CONTRIBUTING
// INSTITUTIONS 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.
//
// --------------------------------------------------------------------------
// $Maintainer: $
// $Authors: Hendrik Weisser $
// --------------------------------------------------------------------------

#ifndef OPENMS_ANALYSIS_MAPMATCHING_TRANSFORMATIONMODEL_H
#define OPENMS_ANALYSIS_MAPMATCHING_TRANSFORMATIONMODEL_H

#include <OpenMS/DATASTRUCTURES/Param.h>

#include <gsl/gsl_bspline.h>
#include <gsl/gsl_interp.h>

namespace OpenMS
{
  /**
       @brief Base class for transformation models

       Implements the identity (no transformation). Parameters and data are ignored.

       @ingroup MapAlignment
  */
  class OPENMS_DLLAPI TransformationModel
  {
public:
    /// Coordinate pair
    typedef std::pair<DoubleReal, DoubleReal> DataPoint;
    /// Vector of coordinate pairs
    typedef std::vector<DataPoint> DataPoints;

    /// Constructor
    TransformationModel() {}

    /// Alternative constructor (derived classes should implement this one!)
    TransformationModel(const TransformationModel::DataPoints &,
                        const Param &) :
      params_() {}

    /// Destructor
    virtual ~TransformationModel() {}

    /// Evaluates the model at the given value
    virtual DoubleReal evaluate(const DoubleReal value) const
    {
      return value;
    }

    /// Gets the (actual) parameters
    void getParameters(Param & params) const
    {
      params = params_;
    }

    /// Gets the default parameters
    static void getDefaultParameters(Param & params)
    {
      params.clear();
    }

protected:
    /// Parameters
    Param params_;
  };

  /**
       @brief Linear model for transformations

       The model can be inferred from data or specified using explicit parameters. If data is given, a least squares fit is used to find the model parameters (slope and intercept). Depending on parameter @p symmetric_regression, a normal regression (@e y on @e x) or symmetric regression (@f$ y - x @f$ on @f$ y + x @f$) is performed.

       Without data, the model can be specified by giving the parameters @p slope and @p intercept explicitly.

       @ingroup MapAlignment
  */
  class OPENMS_DLLAPI TransformationModelLinear :
    public TransformationModel
  {
public:
    /**
        @brief Constructor

        @exception IllegalArgument is thrown if neither data points nor explicit parameters (slope/intercept) are given.
    */
    TransformationModelLinear(const DataPoints & data, const Param & params);

    /// Destructor
    ~TransformationModelLinear();

    /// Evaluates the model at the given value
    virtual DoubleReal evaluate(const DoubleReal value) const;

    using TransformationModel::getParameters;

    /// Gets the "real" parameters
    void getParameters(DoubleReal & slope, DoubleReal & intercept) const;

    /// Gets the default parameters
    static void getDefaultParameters(Param & params);

    /**
         @brief Computes the inverse

         @exception DivisionByZero is thrown if the slope is zero.
    */
    void invert();

protected:
    /// Parameters of the linear model
    DoubleReal slope_, intercept_;
    /// Was the model estimated from data?
    bool data_given_;
    /// Use symmetric regression?
    bool symmetric_;
  };

  /**
       @brief Interpolation model for transformations

       Between the data points, the interpolation uses the neighboring points. Outside the range spanned by the points, we extrapolate using a line through the first and the last point.

       Different types of interpolation (controlled by the parameter @p interpolation_type) are supported: "linear", "polynomial", "cspline", and "akima". Note that the number of required data points may differ between types.

       @ingroup MapAlignment
  */
  class OPENMS_DLLAPI TransformationModelInterpolated :
    public TransformationModel
  {
public:
    /**
         @brief Constructor

         @exception IllegalArgument is thrown if there are not enough data points or if an unknown interpolation type is given.
    */
    TransformationModelInterpolated(const DataPoints & data,
                                    const Param & params);

    /// Destructor
    ~TransformationModelInterpolated();

    /// Evaluates the model at the given value
    DoubleReal evaluate(const DoubleReal value) const;

    /// Gets the default parameters
    static void getDefaultParameters(Param & params);

protected:
    /// Data coordinates
    std::vector<double> x_, y_;
    /// Number of data points
    size_t size_;
    /// Look-up accelerator
    gsl_interp_accel * acc_;
    /// Interpolation function
    gsl_interp * interp_;
    /// Linear model for extrapolation
    TransformationModelLinear * lm_;
  };

  /**
       @brief B-spline model for transformations

       In the range of the data points, the transformation is evaluated from a cubic smoothing spline fit to the points. The number of breakpoints is given as a parameter (@p num_breakpoints). Outside of this range, linear extrapolation through the last point with the slope of the spline at that point is used.

       Positioning of the breakpoints is controlled by the parameter @p break_positions. Valid choices are "uniform" (equidistant spacing on the data range) and "quantiles" (equal numbers of data points in every interval).

       @ingroup MapAlignment
  */
  class OPENMS_DLLAPI TransformationModelBSpline :
    public TransformationModel
  {
public:
    /**
         @brief Constructor

         @exception IllegalArgument is thrown if not enough data points are given or if the required parameter @p num_breakpoints is missing.
    */
    TransformationModelBSpline(const DataPoints & data, const Param & params);

    /// Destructor
    ~TransformationModelBSpline();

    /// Evaluates the model at the given value
    DoubleReal evaluate(const DoubleReal value) const;

    /// Gets the default parameters
    static void getDefaultParameters(Param & params);

protected:
    /**
        @brief Finds quantile values

        @param x Data vector to find quantiles in
        @param quantiles Quantiles to find (values between 0 and 1)
        @param results Resulting quantiles (vector must be already allocated to the correct size!)
    */
    void getQuantiles_(const gsl_vector * x, const std::vector<double> &
                       quantiles, gsl_vector * results);

    /// Computes the B-spline fit
    void computeFit_();

    /// Computes the linear extrapolation
    void computeLinear_(const double pos, double & slope, double & offset,
                        double & sd_err);

    /// Vectors for B-spline computation
    gsl_vector * x_, * y_, * w_, * bsplines_, * coeffs_;
    /// Covariance matrix
    gsl_matrix * cov_;
    /// B-spline workspace
    gsl_bspline_workspace * workspace_;
    /// Number of data points and coefficients
    size_t size_, ncoeffs_;
    // First/last breakpoint
    double xmin_, xmax_;
    /// Parameters for linear extrapolation
    double slope_min_, slope_max_, offset_min_, offset_max_;
    /// Fitting errors of linear extrapolation
    double sd_err_left_, sd_err_right_;
  };

} // end of namespace OpenMS

#endif // OPENMS_ANALYSIS_MAPMATCHING_TRANSFORMATIONMODEL_H

ext: h

max_stars_repo_name: open-ms/all-svn-branches

max_issues_repo_name: kreinert/OpenMS

max_forks_repo_name: kreinert/OpenMS


text:
#pragma once

#define NOMINMAX
#define NODRAWTEXT
#define NOGDI
#define NOBITMAP
#define NOMCX
#define NOSERVICE
#define NOHELP

#include <iostream>
#include <vector>
#include <string>
#include <string_view>
#include <sstream>
#include <variant>
#include <codecvt>
#include <locale>
#include <filesystem>
#include <fstream>
#include <chrono>
#include <functional>
#include <numeric>

#ifndef _WIN32
#include <wsl/winadapter.h>
#include "directml_guids.h"
#else
#include <Windows.h>
#endif
#include <wil/result.h>
#include <wrl/client.h>
#include <gsl/gsl>
#include <rapidjson/document.h>
#include <rapidjson/istreamwrapper.h>
#include <fmt/format.h>

#ifdef _GAMING_XBOX_SCARLETT // Xbox Series X/S - Use GDK
#include <d3d12_xs.h>
#include <d3dx12_xs.h>
#include <d3d12shader_xs.h>
#include <dxcapi_xs.h>
#include <pix3.h>
using IAdapter = IDXGIAdapter;
#elif _GAMING_XBOX_XBOXONE // XboxOne - Use GDK
#include <d3d12_x.h>
#include <d3dx12_x.h>
#include <d3d12shader_x.h>
#include <dxcapi_x.h>
#include <pix3.h>
using IAdapter = IDXGIAdapter;
#else // Desktop/PC - Use DirectX-Headers
#include <directx/d3d12.h>
#include <directx/d3dx12.h>
#include <directx/dxcore.h>
#ifndef DXCOMPILER_NONE
#include <directx/d3d12shader.h>
#endif
#include <dxcapi.h>
#include <WinPixEventRuntime/pix3.h>

#define IGraphicsUnknown IUnknown
#define IID_GRAPHICS_PPV_ARGS IID_PPV_ARGS
using IAdapter = IDXCoreAdapter;
#endif

#include <DirectML.h>
#include "DirectMLX.h"

#include "Logging.h"

max_stars_repo_name: miaobin/DirectML

max_issues_repo_name: miaobin/DirectML

max_forks_repo_name: miaobin/DirectML


text:
/**
 *
 * @file core_ssyrfb.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
 * @date 2010-11-15
 * @generated s Tue Jan  7 11:44:48 2014
 *
 **/
#include <lapacke.h>
#include "common.h"
#undef COMPLEX
#define REAL

/***************************************************************************//**
 *
 * @ingroup CORE_float
 *
 *  CORE_ssyrfb overwrites the symmetric complex N-by-N tile C with
 *
 *    Q**T*C*Q
 *
 *  where Q is a complex unitary matrix defined as the product of k
 *  elementary reflectors
 *
 *    Q = H(1) H(2) . . . H(k)
 *
 *  as returned by CORE_sgeqrt. Only PlasmaLower supported!
 *
 *******************************************************************************
 *
 * @param[in] uplo
 *         @arg PlasmaLower : the upper part of the symmetric matrix C
 *                            is not referenced.
 *         @arg PlasmaUpper : the lower part of the symmetric matrix C
 *                            is not referenced (not supported).
 *
 * @param[in] n
 *         The number of rows/columns of the tile C.  N >= 0.
 *
 * @param[in] k
 *         The number of elementary reflectors whose product defines
 *         the matrix Q. K >= 0.
 *
 * @param[in] ib
 *         The inner-blocking size.  IB >= 0.
 *
 * @param[in] nb
 *         The blocking size.  NB >= 0.
 *
 * @param[in] A
 *         The i-th column must contain the vector which defines the
 *         elementary reflector H(i), for i = 1,2,...,k, as returned by
 *         CORE_sgeqrt in the first k columns of its array argument A.
 *
 * @param[in] lda
 *         The leading dimension of the array A. LDA >= max(1,N).
 *
 * @param[in] T
 *         The IB-by-K triangular factor T of the block reflector.
 *         T is upper triangular by block (economic storage);
 *         The rest of the array is not referenced.
 *
 * @param[in] ldt
 *         The leading dimension of the array T. LDT >= IB.
 *
 * @param[in,out] C
 *         On entry, the symmetric N-by-N tile C.
 *         On exit, C is overwritten by Q**T*C*Q.
 *
 * @param[in] ldc
 *         The leading dimension of the array C. LDC >= max(1,M).
 *
 * @param[in,out] WORK
 *         On exit, if INFO = 0, WORK(1) returns the optimal LDWORK.
 *
 * @param[in] ldwork
 *         The dimension of the array WORK. LDWORK >= max(1,N);
 *
 *******************************************************************************
 *
 * @return
 *          \retval PLASMA_SUCCESS successful exit
 *          \retval <0 if -i, the i-th argument had an illegal value
 *
 ******************************************************************************/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_ssyrfb = PCORE_ssyrfb
#define CORE_ssyrfb PCORE_ssyrfb
#define CORE_sormlq PCORE_sormlq
#define CORE_sormqr PCORE_sormqr
int  CORE_sormlq(PLASMA_enum side, PLASMA_enum trans,
                 int M, int N, int IB, int K,
                 const float *V, int LDV,
                 const float *T, int LDT,
                 float *C, int LDC,
                 float *WORK, int LDWORK);
int  CORE_sormqr(PLASMA_enum side, PLASMA_enum trans,
                 int M, int N, int K, int IB,
                 const float *V, int LDV,
                 const float *T, int LDT,
                 float *C, int LDC,
                 float *WORK, int LDWORK);
#endif
int CORE_ssyrfb( PLASMA_enum uplo, int n,
                 int k, int ib, int nb,
                 const float *A, int lda,
                 const float *T, int ldt,
                 float *C, int ldc,
                 float *WORK, int ldwork )
{
    float tmp;
    int i, j;

    /* Check input arguments */
    if ((uplo != PlasmaUpper) && (uplo != PlasmaLower)) {
        coreblas_error(1, "Illegal value of uplo");
        return -1;
    }
    if (n < 0) {
        coreblas_error(2, "Illegal value of n");
        return -2;
    }
    if (k < 0) {
        coreblas_error(3, "Illegal value of k");
        return -3;
    }
    if (ib < 0) {
        coreblas_error(4, "Illegal value of ib");
        return -4;
    }
    if (nb < 0) {
        coreblas_error(5, "Illegal value of nb");
        return -5;
    }
    if ( (lda < max(1,n)) && (n > 0) ) {
        coreblas_error(7, "Illegal value of lda");
        return -7;
    }
    if ( (ldt < max(1,ib)) && (ib > 0) ) {
        coreblas_error(9, "Illegal value of ldt");
        return -9;
    }
    if ( (ldc < max(1,n)) && (n > 0) ) {
        coreblas_error(11, "Illegal value of ldc");
        return -11;
    }

    if (uplo == PlasmaLower) {
        /* Rebuild the symmetric block: WORK <- C */
        for (j = 0; j < n; j++) {
            *(WORK + j + j * ldwork) =  *(C + ldc*j + j);
            for (i = j+1; i < n; i++){
                tmp = *(C + i + j*ldc);
                *(WORK + i + j * ldwork) = tmp;
                *(WORK + j + i * ldwork) = ( tmp );
            }
        }

        /* Left */
        CORE_sormqr(PlasmaLeft, PlasmaTrans, n, n, k, ib,
                    A, lda, T, ldt, WORK, ldwork, WORK+nb*ldwork, ldwork);
        /* Right */
        CORE_sormqr(PlasmaRight, PlasmaNoTrans, n, n, k, ib,
                    A, lda, T, ldt, WORK, ldwork, WORK+nb*ldwork, ldwork);

        /*
         * Copy back the final result to the lower part of C
         */
        LAPACKE_slacpy_work( LAPACK_COL_MAJOR, lapack_const(PlasmaLower), n, n, WORK, ldwork, C, ldc );
    }
    else {
        /* Rebuild the symmetric block: WORK <- C */
        for (j = 0; j < n; j++) {
            for (i = 0; i < j; i++){
                tmp = *(C + i + j*ldc);
                *(WORK + i + j * ldwork) = tmp;
                *(WORK + j + i * ldwork) = ( tmp );
            }
            *(WORK + j + j * ldwork) =  *(C + ldc*j + j);
        }

        /* Right */
        CORE_sormlq(PlasmaRight, PlasmaTrans, n, n, k, ib,
                    A, lda, T, ldt, WORK, ldwork, WORK+nb*ldwork, ldwork);
        /* Left */
        CORE_sormlq(PlasmaLeft, PlasmaNoTrans, n, n, k, ib,
                    A, lda, T, ldt, WORK, ldwork, WORK+nb*ldwork, ldwork);

        /*
         * Copy back the final result to the upper part of C
         */
        LAPACKE_slacpy_work( LAPACK_COL_MAJOR, lapack_const(PlasmaUpper), n, n, WORK, ldwork, C, ldc );
    }
    return 0;
}


text:
#include <string.h>
#include <gbpLib.h>
#include <gbpMisc.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>

void invert_square_matrix(double *matrix_in, int size, double *inverse_out) {
    gsl_matrix *     matrix;
    gsl_matrix *     LU_decomp;
    gsl_matrix *     inverse;
    gsl_permutation *permutation;
    int              signum;
    gsl_matrix_view  m;
    double *         matrix_copy;
    int              i_x, i_y, i_M;

    // Allocate some RAM and set some pointers
    matrix_copy = (double *)SID_malloc(sizeof(double) * size * size);
    memcpy(matrix_copy, matrix_in, sizeof(double) * size * size);
    m           = gsl_matrix_view_array(matrix_copy, size, size);
    matrix      = &m.matrix;
    LU_decomp   = matrix;
    inverse     = gsl_matrix_alloc(size, size);
    permutation = gsl_permutation_alloc(size);

    // Perform LU decomposition.  Result ends-up in matrix,
    //   which is pointed to by LU_decomp
    gsl_linalg_LU_decomp(matrix, permutation, &signum);
    gsl_linalg_LU_invert(LU_decomp, permutation, inverse);

    // Copy result into the out-going array
    for(i_y = 0, i_M = 0; i_y < size; i_y++) {
        for(i_x = 0; i_x < size; i_x++, i_M++) {
            inverse_out[i_M] = gsl_matrix_get(inverse, i_x, i_y);
        }
    }

    // Test that the matrix*inverse=identity
    /*
    for(i_x=0;i_x<size;i_x++){
       for(i_y=0;i_y<size;i_y++,i_M++){
          double test;
          int j_M;
          test=0.;
          for(j_M=0;j_M<size;j_M++) test+=matrix_in[i_y*size+j_M]*inverse_out[j_M*size+i_x];
          fprintf(stderr,"%le ",test);
       }
       fprintf(stderr,"\n");
    }
    */

    // Free some RAM
    gsl_matrix_free(inverse);
    gsl_permutation_free(permutation);
    SID_free(SID_FARG matrix_copy);
}

max_stars_repo_name: gbpoole/gbpCode

max_issues_repo_name: gbpoole/gbpCode

max_forks_repo_name: gbpoole/gbpCode


text:
/*
 * Implement Heap sort -- direct and indirect sorting
 * Based on descriptions in Sedgewick "Algorithms in C"
 *
 * Copyright (C) 1999  Thomas Walter
 *
 * 18 February 2000: Modified for GSL by Brian Gough
 *
 * This is free software; you can redistribute it and/or modify it
 * under the terms of the GNU General Public License as published by the
 * Free Software Foundation; either version 3, or (at your option) any
 * later version.
 *
 * This source is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
 * for more details.
 */

#include <config.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_sort_vector.h>

#define BASE_LONG_DOUBLE
#include "templates_on.h"
#include "sortvecind_source.c"
#include "templates_off.h"
#undef  BASE_LONG_DOUBLE

#define BASE_DOUBLE
#include "templates_on.h"
#include "sortvecind_source.c"
#include "templates_off.h"
#undef  BASE_DOUBLE

#define BASE_FLOAT
#include "templates_on.h"
#include "sortvecind_source.c"
#include "templates_off.h"
#undef  BASE_FLOAT

#define BASE_ULONG
#include "templates_on.h"
#include "sortvecind_source.c"
#include "templates_off.h"
#undef  BASE_ULONG

#define BASE_LONG
#include "templates_on.h"
#include "sortvecind_source.c"
#include "templates_off.h"
#undef  BASE_LONG

#define BASE_UINT
#include "templates_on.h"
#include "sortvecind_source.c"
#include "templates_off.h"
#undef  BASE_UINT

#define BASE_INT
#include "templates_on.h"
#include "sortvecind_source.c"
#include "templates_off.h"
#undef  BASE_INT

#define BASE_USHORT
#include "templates_on.h"
#include "sortvecind_source.c"
#include "templates_off.h"
#undef  BASE_USHORT

#define BASE_SHORT
#include "templates_on.h"
#include "sortvecind_source.c"
#include "templates_off.h"
#undef  BASE_SHORT

#define BASE_UCHAR
#include "templates_on.h"
#include "sortvecind_source.c"
#include "templates_off.h"
#undef  BASE_UCHAR

#define BASE_CHAR
#include "templates_on.h"
#include "sortvecind_source.c"
#include "templates_off.h"
#undef  BASE_CHAR

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#if !defined(PETIGA_H)
#define PETIGA_H

/*
#include <petscconf.h>
#undef  static inline
#define static inline static __inline
*/

#include <petsc.h>

#if PETSC_VERSION_LT(3,7,0)
#error "PetIGA requires PETSc 3.7 or higher"
#endif

#include <petsc/private/petscimpl.h>

#if PETSC_VERSION_LT(3,8,0)
#  ifdef PETSC_FUNCTION_NAME
#    undef  __FUNCT__
#    define __FUNCT__ PETSC_FUNCTION_NAME
#  endif
#endif

#if PETSC_VERSION_LT(3,17,0)
#undef  SETERRQ
#define SETERRQ(comm,ierr,...) return PetscError(comm,__LINE__,PETSC_FUNCTION_NAME,__FILE__,ierr,PETSC_ERROR_INITIAL,__VA_ARGS__)
#endif

typedef ISLocalToGlobalMapping LGMap;
#define LGMap LGMap

/* ---------------------------------------------------------------- */

typedef struct _p_IGA         *IGA;
typedef struct _n_IGAAxis     *IGAAxis;
typedef struct _n_IGARule     *IGARule;
typedef struct _n_IGABasis    *IGABasis;
typedef struct _n_IGAForm     *IGAForm;
typedef struct _n_IGAElement  *IGAElement;
typedef struct _n_IGAPoint    *IGAPoint;
typedef struct _n_IGAProbe    *IGAProbe;

/* ---------------------------------------------------------------- */

struct _n_IGAAxis {
  PetscInt refct;
  /**/
  PetscInt   p; /* polynomial order    */
  PetscInt   m; /* last knot index     */
  PetscReal *U; /* knot vector         */
  /**/
  PetscBool  periodic; /* periodicity  */
  PetscInt   nnp,nel;  /* bases, spans */
  PetscInt   *span;    /* span indices */
};

PETSC_EXTERN PetscErrorCode IGAAxisCreate(IGAAxis *axis);
PETSC_EXTERN PetscErrorCode IGAAxisDestroy(IGAAxis *axis);
PETSC_EXTERN PetscErrorCode IGAAxisReset(IGAAxis axis);
PETSC_EXTERN PetscErrorCode IGAAxisReference(IGAAxis axis);
PETSC_EXTERN PetscErrorCode IGAAxisCopy(IGAAxis base,IGAAxis axis);
PETSC_EXTERN PetscErrorCode IGAAxisDuplicate(IGAAxis base,IGAAxis *axis);
PETSC_EXTERN PetscErrorCode IGAAxisSetPeriodic(IGAAxis axis,PetscBool periodic);
PETSC_EXTERN PetscErrorCode IGAAxisGetPeriodic(IGAAxis axis,PetscBool *periodic);
PETSC_EXTERN PetscErrorCode IGAAxisSetDegree(IGAAxis axis,PetscInt p);
PETSC_EXTERN PetscErrorCode IGAAxisGetDegree(IGAAxis axis,PetscInt *p);
PETSC_EXTERN PetscErrorCode IGAAxisSetKnots(IGAAxis axis,PetscInt m,const PetscReal U[]);
PETSC_EXTERN PetscErrorCode IGAAxisGetKnots(IGAAxis axis,PetscInt *m,PetscReal *U[]);
PETSC_EXTERN PetscErrorCode IGAAxisGetLimits(IGAAxis axis,PetscReal *Ui,PetscReal *Uf);
PETSC_EXTERN PetscErrorCode IGAAxisGetSizes(IGAAxis axis,PetscInt *nel,PetscInt *nnp);
PETSC_EXTERN PetscErrorCode IGAAxisGetSpans(IGAAxis axis,PetscInt *nel,PetscInt *spans[]);
PETSC_EXTERN PetscErrorCode IGAAxisInit(IGAAxis axis,PetscInt p,PetscInt m,const PetscReal U[]);
PETSC_EXTERN PetscErrorCode IGAAxisInitBreaks(IGAAxis axis,PetscInt nu,const PetscReal u[],PetscInt C);
PETSC_EXTERN PetscErrorCode IGAAxisInitUniform(IGAAxis axis,PetscInt N,PetscReal Ui,PetscReal Uf,PetscInt C);
PETSC_EXTERN PetscErrorCode IGAAxisSetUp(IGAAxis axis);

typedef enum {
  IGA_RULE_LEGENDRE=0, /* Gauss-Legendre         */
  IGA_RULE_LOBATTO,    /* Gauss-Lobatto          */
  IGA_RULE_REDUCED,    /* Reduced Gauss-Legendre */
  IGA_RULE_USER        /* User-defined           */
} IGARuleType;

PETSC_EXTERN const char *const IGARuleTypes[];

struct _n_IGARule {
  PetscInt refct;
  /**/
  IGARuleType type;   /* rule type */
  /**/
  PetscInt  nqp;      /* number of quadrature points */
  PetscReal *point;   /* [nqp] quadrature points  */
  PetscReal *weight;  /* [nqp] quadrature weights */
};

PETSC_EXTERN PetscErrorCode IGARuleCreate(IGARule *rule);
PETSC_EXTERN PetscErrorCode IGARuleDestroy(IGARule *rule);
PETSC_EXTERN PetscErrorCode IGARuleReset(IGARule rule);
PETSC_EXTERN PetscErrorCode IGARuleReference(IGARule rule);
PETSC_EXTERN PetscErrorCode IGARuleCopy(IGARule base,IGARule rule);
PETSC_EXTERN PetscErrorCode IGARuleDuplicate(IGARule base,IGARule *rule);
PETSC_EXTERN PetscErrorCode IGARuleSetType(IGARule rule,IGARuleType type);
PETSC_EXTERN PetscErrorCode IGARuleSetSize(IGARule rule,PetscInt nqp);
PETSC_EXTERN PetscErrorCode IGARuleSetUp(IGARule rule);
PETSC_EXTERN PetscErrorCode IGARuleSetRule(IGARule rule,PetscInt q,const PetscReal x[],const PetscReal w[]);
PETSC_EXTERN PetscErrorCode IGARuleGetRule(IGARule rule,PetscInt *q,PetscReal *x[],PetscReal *w[]);

typedef enum {
  IGA_BASIS_BSPLINE=0, /* B-Spline basis functions */
  IGA_BASIS_BERNSTEIN, /* Bernstein polynomials (C^0 B-Spline basis) */
  IGA_BASIS_LAGRANGE,  /* Lagrange polynomials on Newton-Cotes points */
  IGA_BASIS_SPECTRAL   /* Lagrange polynomials on Gauss-Lobatto points */
} IGABasisType;

PETSC_EXTERN const char *const IGABasisTypes[];

struct _n_IGABasis {
  PetscInt refct;
  /**/
  IGABasisType type;  /* basis type */
  /**/
  PetscInt  nel;      /* number of elements */
  PetscInt  nqp;      /* number of quadrature points */
  PetscInt  nen;      /* number of local basis functions */

  PetscInt  *offset;  /* [nel] basis offset */
  PetscReal *detJac;  /* [nel] element length */
  PetscReal *weight;  /* [nel][nqp] quadrature weight */
  PetscReal *point;   /* [nel][nqp] quadrature point */
  PetscReal *value;   /* [nel][nqp][nen][5] basis derivatives */

  PetscReal  bnd_detJac;
  PetscReal  bnd_weight;
  PetscReal  bnd_point[2];
  PetscReal *bnd_value[2]; /* [nen][5] */
};

PETSC_EXTERN PetscErrorCode IGABasisCreate(IGABasis *basis);
PETSC_EXTERN PetscErrorCode IGABasisDestroy(IGABasis *basis);
PETSC_EXTERN PetscErrorCode IGABasisReset(IGABasis basis);
PETSC_EXTERN PetscErrorCode IGABasisReference(IGABasis basis);
PETSC_EXTERN PetscErrorCode IGABasisSetType(IGABasis basis,IGABasisType type);
PETSC_EXTERN PetscErrorCode IGABasisInitQuadrature (IGABasis basis,IGAAxis axis,IGARule rule);
PETSC_EXTERN PetscErrorCode IGABasisInitCollocation(IGABasis basis,IGAAxis axis);

/* ---------------------------------------------------------------- */

typedef PetscErrorCode (*IGAFormExact)(IGAPoint point,PetscInt k,PetscScalar V[],void *ctx);
typedef PetscErrorCode (*IGAFormScalar)(IGAPoint point,const PetscScalar U[],PetscInt n,PetscScalar S[],void *ctx);
typedef PetscErrorCode (*IGAFormVector)(IGAPoint point,PetscScalar F[],void *ctx);
typedef PetscErrorCode (*IGAFormMatrix)(IGAPoint point,PetscScalar K[],void *ctx);
typedef PetscErrorCode (*IGAFormSystem)(IGAPoint point,PetscScalar K[],PetscScalar F[],void *ctx);
typedef PetscErrorCode (*IGAFormFunction)(IGAPoint point,
                                          const PetscScalar U[],
                                          PetscScalar F[],void *ctx);
typedef PetscErrorCode (*IGAFormJacobian)(IGAPoint point,
                                          const PetscScalar U[],
                                          PetscScalar J[],void *ctx);
typedef PetscErrorCode (*IGAFormIFunction)(IGAPoint point,
                                           PetscReal a,const PetscScalar V[],
                                           PetscReal t,const PetscScalar U[],
                                           PetscScalar F[],void *ctx);
typedef PetscErrorCode (*IGAFormIJacobian)(IGAPoint point,
                                           PetscReal a,const PetscScalar V[],
                                           PetscReal t,const PetscScalar U[],
                                           PetscScalar J[],void *ctx);
typedef PetscErrorCode (*IGAFormI2Function)(IGAPoint point,
                                            PetscReal a,const PetscScalar A[],
                                            PetscReal v,const PetscScalar V[],
                                            PetscReal t,const PetscScalar U[],
                                            PetscScalar F[],void *ctx);
typedef PetscErrorCode (*IGAFormI2Jacobian)(IGAPoint point,
                                            PetscReal a,const PetscScalar A[],
                                            PetscReal v,const PetscScalar V[],
                                            PetscReal t,const PetscScalar U[],
                                            PetscScalar J[],void *ctx);
typedef PetscErrorCode (*IGAFormIEFunction)(IGAPoint point,
                                            PetscReal a,const PetscScalar V[],
                                            PetscReal t,const PetscScalar U[],
                                            PetscReal t0,const PetscScalar U0[],
                                            PetscScalar F[],void *ctx);
typedef PetscErrorCode (*IGAFormIEJacobian)(IGAPoint point,
                                            PetscReal a,const PetscScalar V[],
                                            PetscReal t,const PetscScalar U[],
                                            PetscReal t0,const PetscScalar U0[],
                                            PetscScalar J[],void *ctx);
typedef PetscErrorCode (*IGAFormRHSFunction)(IGAPoint point,
                                             PetscReal t,const PetscScalar U[],
                                             PetscScalar F[],void *ctx);
typedef PetscErrorCode (*IGAFormRHSJacobian)(IGAPoint point,
                                             PetscReal t,const PetscScalar U[],
                                             PetscScalar J[],void *ctx);

PETSC_EXTERN PetscErrorCode IGAFormJacobianFD(IGAPoint p,
                                              const PetscScalar U[],
                                              PetscScalar J[],void *ctx);
PETSC_EXTERN PetscErrorCode IGAFormIJacobianFD(IGAPoint p,
                                               PetscReal s,const PetscScalar V[],
                                               PetscReal t,const PetscScalar U[],
                                               PetscScalar J[],void *ctx);
PETSC_EXTERN PetscErrorCode IGAFormIEJacobianFD(IGAPoint p,
                                                PetscReal s, const PetscScalar V[],
                                                PetscReal t, const PetscScalar U[],
                                                PetscReal t0,const PetscScalar U0[],
                                                PetscScalar J[],void *ctx);
PETSC_EXTERN PetscErrorCode IGAFormI2JacobianFD(IGAPoint p,
                                                PetscReal a,const PetscScalar A[],
                                                PetscReal v,const PetscScalar V[],
                                                PetscReal t,const PetscScalar U[],
                                                PetscScalar J[],void *ctx);
PETSC_EXTERN PetscErrorCode IGAFormRHSJacobianFD(IGAPoint p,
                                                 PetscReal t,const PetscScalar U[],
                                                 PetscScalar F[],void *ctx);

typedef struct _IGAFormBC *IGAFormBC;
struct _IGAFormBC {
  PetscInt    count;
  PetscInt    field[64];
  PetscScalar value[64];
};

typedef struct _IGAFormOps *IGAFormOps;
struct _IGAFormOps {
  /**/
  IGAFormVector     Vector;
  void              *VecCtx;
  IGAFormMatrix     Matrix;
  void              *MatCtx;
  IGAFormSystem     System;
  void              *SysCtx;
  /**/
  IGAFormFunction   Function;
  void              *FunCtx;
  IGAFormJacobian   Jacobian;
  void              *JacCtx;
  /**/
  IGAFormIFunction  IFunction;
  IGAFormI2Function I2Function;
  void              *IFunCtx;
  IGAFormIJacobian  IJacobian;
  IGAFormI2Jacobian I2Jacobian;
  void              *IJacCtx;
  /**/
  IGAFormIEFunction IEFunction;
  void              *IEFunCtx;
  IGAFormIEJacobian IEJacobian;
  void              *IEJacCtx;
  /**/
  IGAFormRHSFunction RHSFunction;
  void               *RHSFunCtx;
  IGAFormRHSJacobian RHSJacobian;
  void               *RHSJacCtx;
};

struct _n_IGAForm {
  PetscInt refct;
  /**/
  IGAFormOps ops;
  PetscInt   dof;
  IGAFormBC  value[3][2];
  IGAFormBC  load [3][2];
  PetscBool  visit[3][2];
};

PETSC_EXTERN PetscErrorCode IGAGetForm(IGA iga,IGAForm *form);
PETSC_EXTERN PetscErrorCode IGASetForm(IGA iga,IGAForm form);

PETSC_EXTERN PetscErrorCode IGAFormCreate(IGAForm *form);
PETSC_EXTERN PetscErrorCode IGAFormDestroy(IGAForm *form);
PETSC_EXTERN PetscErrorCode IGAFormReset(IGAForm form);
PETSC_EXTERN PetscErrorCode IGAFormReference(IGAForm form);

PETSC_EXTERN PetscErrorCode IGAFormSetBoundaryValue(IGAForm form,PetscInt axis,PetscInt side,PetscInt field,PetscScalar value);
PETSC_EXTERN PetscErrorCode IGAFormSetBoundaryLoad (IGAForm form,PetscInt axis,PetscInt side,PetscInt field,PetscScalar value);
PETSC_EXTERN PetscErrorCode IGAFormSetBoundaryForm (IGAForm form,PetscInt axis,PetscInt side,PetscBool flag);
PETSC_EXTERN PetscErrorCode IGAFormClearBoundary   (IGAForm form,PetscInt axis,PetscInt side);

PETSC_EXTERN PetscErrorCode IGAFormSetVector     (IGAForm form,IGAFormVector      Vector,     void *ctx);
PETSC_EXTERN PetscErrorCode IGAFormSetMatrix     (IGAForm form,IGAFormMatrix      Matrix,     void *ctx);
PETSC_EXTERN PetscErrorCode IGAFormSetSystem     (IGAForm form,IGAFormSystem      System,     void *ctx);
PETSC_EXTERN PetscErrorCode IGAFormSetFunction   (IGAForm form,IGAFormFunction    Function,   void *ctx);
PETSC_EXTERN PetscErrorCode IGAFormSetJacobian   (IGAForm form,IGAFormJacobian    Jacobian,   void *ctx);
PETSC_EXTERN PetscErrorCode IGAFormSetIFunction  (IGAForm form,IGAFormIFunction   IFunction,  void *ctx);
PETSC_EXTERN PetscErrorCode IGAFormSetIJacobian  (IGAForm form,IGAFormIJacobian   IJacobian,  void *ctx);
PETSC_EXTERN PetscErrorCode IGAFormSetI2Function (IGAForm form,IGAFormI2Function  IFunction,  void *ctx);
PETSC_EXTERN PetscErrorCode IGAFormSetI2Jacobian (IGAForm form,IGAFormI2Jacobian  IJacobian,  void *ctx);
PETSC_EXTERN PetscErrorCode IGAFormSetIEFunction (IGAForm form,IGAFormIEFunction  IEFunction, void *ctx);
PETSC_EXTERN PetscErrorCode IGAFormSetIEJacobian (IGAForm form,IGAFormIEJacobian  IEJacobian, void *ctx);
PETSC_EXTERN PetscErrorCode IGAFormSetRHSFunction(IGAForm form,IGAFormRHSFunction RHSFunction,void *ctx);
PETSC_EXTERN PetscErrorCode IGAFormSetRHSJacobian(IGAForm form,IGAFormRHSJacobian RHSJacobian,void *ctx);

PETSC_EXTERN PetscErrorCode IGASetFixTable(IGA iga,Vec table);
PETSC_EXTERN PetscErrorCode IGASetBoundaryValue(IGA iga,PetscInt axis,PetscInt side,PetscInt field,PetscScalar value);
PETSC_EXTERN PetscErrorCode IGASetBoundaryLoad (IGA iga,PetscInt axis,PetscInt side,PetscInt field,PetscScalar value);
PETSC_EXTERN PetscErrorCode IGASetBoundaryForm (IGA iga,PetscInt axis,PetscInt side,PetscBool flag);

PETSC_EXTERN PetscErrorCode IGASetFormVector     (IGA iga,IGAFormVector      Vector,     void *ctx);
PETSC_EXTERN PetscErrorCode IGASetFormMatrix     (IGA iga,IGAFormMatrix      Matrix,     void *ctx);
PETSC_EXTERN PetscErrorCode IGASetFormSystem     (IGA iga,IGAFormSystem      System,     void *ctx);
PETSC_EXTERN PetscErrorCode IGASetFormFunction   (IGA iga,IGAFormFunction    Function,   void *ctx);
PETSC_EXTERN PetscErrorCode IGASetFormJacobian   (IGA iga,IGAFormJacobian    Jacobian,   void *ctx);
PETSC_EXTERN PetscErrorCode IGASetFormIFunction  (IGA iga,IGAFormIFunction   IFunction,  void *ctx);
PETSC_EXTERN PetscErrorCode IGASetFormIJacobian  (IGA iga,IGAFormIJacobian   IJacobian,  void *ctx);
PETSC_EXTERN PetscErrorCode IGASetFormI2Function (IGA iga,IGAFormI2Function  IFunction,  void *ctx);
PETSC_EXTERN PetscErrorCode IGASetFormI2Jacobian (IGA iga,IGAFormI2Jacobian  IJacobian,  void *ctx);
PETSC_EXTERN PetscErrorCode IGASetFormIEFunction (IGA iga,IGAFormIEFunction  IEFunction, void *ctx);
PETSC_EXTERN PetscErrorCode IGASetFormIEJacobian (IGA iga,IGAFormIEJacobian  IEJacobian, void *ctx);
PETSC_EXTERN PetscErrorCode IGASetFormRHSFunction(IGA iga,IGAFormRHSFunction RHSFunction,void *ctx);
PETSC_EXTERN PetscErrorCode IGASetFormRHSJacobian(IGA iga,IGAFormRHSJacobian RHSJacobian,void *ctx);

/* ---------------------------------------------------------------- */

typedef struct _IGAOps *IGAOps;
struct _IGAOps {
  PetscErrorCode (*create)(IGA);
  PetscErrorCode (*destroy)(IGA);
  PetscErrorCode (*view)(IGA);
  PetscErrorCode (*setup)(IGA);
  PetscErrorCode (*setfromoptions)(IGA);
};

struct _p_IGA {
  PETSCHEADER(struct _IGAOps);
  PetscBool  setup;
  PetscInt   setupstage;
  PetscBool  collocation;

  VecType    vectype;
  MatType    mattype;
  char       **fieldname;

  PetscInt  dim;   /* parametric dimension of the function space*/
  PetscInt  dof;   /* number of degrees of freedom per node */
  PetscInt  order; /* maximum derivative order */

  IGAAxis     axis[3];
  IGARule     rule[3];
  IGABasis    basis[3];
  IGAForm     form;

  IGAElement  iterator;

  PetscBool   rational;
  PetscInt    geometry;
  PetscInt    property;
  PetscReal   *rationalW;
  PetscReal   *geometryX;
  PetscScalar *propertyA;

  PetscBool   fixtable;
  PetscScalar *fixtableU;

  PetscInt  proc_sizes[3];
  PetscInt  proc_ranks[3];

  PetscInt  elem_sizes[3];
  PetscInt  elem_start[3];
  PetscInt  elem_width[3];

  PetscInt  geom_sizes[3];
  PetscInt  geom_lstart[3];
  PetscInt  geom_lwidth[3];
  PetscInt  geom_gstart[3];
  PetscInt  geom_gwidth[3];

  PetscInt  node_shift[3];
  PetscInt  node_sizes[3];
  PetscInt  node_lstart[3];
  PetscInt  node_lwidth[3];
  PetscInt  node_gstart[3];
  PetscInt  node_gwidth[3];

  AO          ao;
  LGMap       lgmap;
  PetscLayout map;
  VecScatter  g2l,l2g,l2l;
  PetscInt    nwork;
  Vec         vwork[16];
  Vec         natural;
  VecScatter  n2g,g2n;

  DM elem_dm;
  DM geom_dm;
  DM node_dm;
  DM draw_dm;
};

PETSC_EXTERN PetscClassId IGA_CLASSID;
#define IGA_FILE_CLASSID 1211299

PETSC_EXTERN PetscErrorCode IGAInitializePackage(void);
PETSC_EXTERN PetscErrorCode IGAFinalizePackage(void);
PETSC_EXTERN PetscErrorCode IGARegisterAll(void);

PETSC_EXTERN PetscErrorCode IGACreate(MPI_Comm comm,IGA *iga);
PETSC_EXTERN PetscErrorCode IGADestroy(IGA *iga);
PETSC_EXTERN PetscErrorCode IGAReset(IGA iga);
PETSC_EXTERN PetscErrorCode IGASetUp(IGA iga);
PETSC_EXTERN PetscErrorCode IGAView(IGA iga,PetscViewer viewer);

PETSC_EXTERN PetscErrorCode IGAGetOptionsPrefix(IGA iga,const char *prefix[]);
PETSC_EXTERN PetscErrorCode IGASetOptionsPrefix(IGA iga,const char prefix[]);
PETSC_EXTERN PetscErrorCode IGAPrependOptionsPrefix(IGA iga,const char prefix[]);
PETSC_EXTERN PetscErrorCode IGAAppendOptionsPrefix(IGA iga,const char prefix[]);
PETSC_EXTERN PetscErrorCode IGASetFromOptions(IGA iga);
PETSC_EXTERN PetscErrorCode IGAViewFromOptions(IGA iga,PetscObject obj,const char name[]);
PETSC_EXTERN PetscErrorCode IGAOptionsAlias(const char name[],const char defval[],const char alias[]);

PETSC_EXTERN PetscErrorCode IGALoad(IGA iga,PetscViewer viewer);
PETSC_EXTERN PetscErrorCode IGASave(IGA iga,PetscViewer viewer);
PETSC_EXTERN PetscErrorCode IGARead(IGA iga,const char filename[]);
PETSC_EXTERN PetscErrorCode IGAWrite(IGA iga,const char filename[]);
PETSC_EXTERN PetscErrorCode IGAPrint(IGA iga,PetscViewer viewer);
PETSC_EXTERN PetscErrorCode IGADraw(IGA iga,PetscViewer viewer);

PETSC_EXTERN PetscErrorCode IGASetGeometryDim(IGA iga,PetscInt dim);
PETSC_EXTERN PetscErrorCode IGAGetGeometryDim(IGA iga,PetscInt *dim);
PETSC_EXTERN PetscErrorCode IGALoadGeometry(IGA iga,PetscViewer viewer);
PETSC_EXTERN PetscErrorCode IGASaveGeometry(IGA iga,PetscViewer viewer);

PETSC_EXTERN PetscErrorCode IGASetPropertyDim(IGA iga,PetscInt dim);
PETSC_EXTERN PetscErrorCode IGAGetPropertyDim(IGA iga,PetscInt *dim);
PETSC_EXTERN PetscErrorCode IGALoadProperty(IGA iga,PetscViewer viewer);
PETSC_EXTERN PetscErrorCode IGASaveProperty(IGA iga,PetscViewer viewer);

PETSC_EXTERN PetscErrorCode IGALoadVec(IGA iga,Vec vec,PetscViewer viewer);
PETSC_EXTERN PetscErrorCode IGASaveVec(IGA iga,Vec vec,PetscViewer viewer);
PETSC_EXTERN PetscErrorCode IGADrawVec(IGA iga,Vec vec,PetscViewer viewer);
PETSC_EXTERN PetscErrorCode IGAReadVec(IGA iga,Vec vec,const char filename[]);
PETSC_EXTERN PetscErrorCode IGAWriteVec(IGA iga,Vec vec,const char filename[]);
PETSC_EXTERN PetscErrorCode IGADrawVecVTK(IGA iga,Vec vec,const char filename[]);

PETSC_EXTERN PetscErrorCode IGASetDim(IGA iga,PetscInt dim);
PETSC_EXTERN PetscErrorCode IGAGetDim(IGA iga,PetscInt *dim);
PETSC_EXTERN PetscErrorCode IGASetDof(IGA iga,PetscInt dof);
PETSC_EXTERN PetscErrorCode IGAGetDof(IGA iga,PetscInt *dof);
PETSC_EXTERN PetscErrorCode IGASetName(IGA iga,const char name[]);
PETSC_EXTERN PetscErrorCode IGAGetName(IGA iga,const char *name[]);
PETSC_EXTERN PetscErrorCode IGASetFieldName(IGA iga,PetscInt field,const char name[]);
PETSC_EXTERN PetscErrorCode IGAGetFieldName(IGA iga,PetscInt field,const char *name[]);
PETSC_EXTERN PetscErrorCode IGASetOrder(IGA iga,PetscInt order);
PETSC_EXTERN PetscErrorCode IGAGetOrder(IGA iga,PetscInt *order);
PETSC_EXTERN PetscErrorCode IGASetProcessors(IGA iga,PetscInt i,PetscInt processors);
PETSC_EXTERN PetscErrorCode IGASetBasisType(IGA iga,PetscInt i,IGABasisType type);
PETSC_EXTERN PetscErrorCode IGASetRuleType(IGA iga,PetscInt i,IGARuleType type);
PETSC_EXTERN PetscErrorCode IGASetRuleSize(IGA iga,PetscInt i,PetscInt nqp);
PETSC_EXTERN PetscErrorCode IGASetQuadrature(IGA iga,PetscInt i,PetscInt q);
PETSC_EXTERN PetscErrorCode IGASetUseCollocation(IGA iga,PetscBool collocation);

PETSC_EXTERN PetscErrorCode IGAGetComm(IGA iga,MPI_Comm *comm);
PETSC_EXTERN PetscErrorCode IGAGetAxis(IGA iga,PetscInt i,IGAAxis *axis);
PETSC_EXTERN PetscErrorCode IGAGetRule(IGA iga,PetscInt i,IGARule *rule);
PETSC_EXTERN PetscErrorCode IGAGetBasis(IGA iga,PetscInt i,IGABasis *basis);

PETSC_EXTERN PetscErrorCode IGACreateDMDA(IGA iga,PetscInt bs,
                                          const PetscInt gsizes[],
                                          const PetscInt lsizes[],
                                          const PetscBool periodic[],
                                          PetscBool stencil_box,
                                          PetscInt stencil_width,
                                          DM *dm);
PETSC_EXTERN PetscErrorCode IGACreateElemDM(IGA iga,PetscInt bs,DM *dm);
PETSC_EXTERN PetscErrorCode IGACreateGeomDM(IGA iga,PetscInt bs,DM *dm);
PETSC_EXTERN PetscErrorCode IGACreateNodeDM(IGA iga,PetscInt bs,DM *dm);
PETSC_EXTERN PetscErrorCode IGAGetElemDM(IGA iga,DM *dm);
PETSC_EXTERN PetscErrorCode IGAGetGeomDM(IGA iga,DM *dm);
PETSC_EXTERN PetscErrorCode IGAGetNodeDM(IGA iga,DM *dm);

PETSC_EXTERN PetscErrorCode IGASetVecType(IGA iga,const VecType vectype);
PETSC_EXTERN PetscErrorCode IGASetMatType(IGA iga,const MatType mattype);

PETSC_EXTERN PetscErrorCode IGACreateVec(IGA iga,Vec *vec);
PETSC_EXTERN PetscErrorCode IGACreateMat(IGA iga,Mat *mat);

PETSC_EXTERN PetscErrorCode IGACreateCoordinates(IGA iga,Vec *coords);
PETSC_EXTERN PetscErrorCode IGACreateRigidBody(IGA iga,MatNullSpace *nsp);

PETSC_EXTERN PetscErrorCode IGACreateLocalVec(IGA iga, Vec *lvec);
PETSC_EXTERN PetscErrorCode IGAGetLocalVec(IGA iga,Vec *lvec);
PETSC_EXTERN PetscErrorCode IGARestoreLocalVec(IGA iga,Vec *lvec);
PETSC_EXTERN PetscErrorCode IGAGlobalToLocalBegin(IGA iga,Vec gvec,Vec lvec,InsertMode addv);
PETSC_EXTERN PetscErrorCode IGAGlobalToLocalEnd  (IGA iga,Vec gvec,Vec lvec,InsertMode addv);
PETSC_EXTERN PetscErrorCode IGAGlobalToLocal     (IGA iga,Vec gvec,Vec lvec,InsertMode addv);
PETSC_EXTERN PetscErrorCode IGALocalToGlobalBegin(IGA iga,Vec lvec,Vec gvec,InsertMode addv);
PETSC_EXTERN PetscErrorCode IGALocalToGlobalEnd  (IGA iga,Vec lvec,Vec gvec,InsertMode addv);
PETSC_EXTERN PetscErrorCode IGALocalToGlobal     (IGA iga,Vec lvec,Vec gvec,InsertMode addv);
PETSC_EXTERN PetscErrorCode IGALocalToLocalBegin (IGA iga,Vec gvec,Vec lvec,InsertMode addv);
PETSC_EXTERN PetscErrorCode IGALocalToLocalEnd   (IGA iga,Vec gvec,Vec lvec,InsertMode addv);
PETSC_EXTERN PetscErrorCode IGALocalToLocal      (IGA iga,Vec gvec,Vec lvec,InsertMode addv);

PETSC_EXTERN PetscErrorCode IGAGetNaturalVec(IGA iga,Vec *nvec);
PETSC_EXTERN PetscErrorCode IGANaturalToGlobal(IGA iga,Vec nvec,Vec gvec);
PETSC_EXTERN PetscErrorCode IGAGlobalToNatural(IGA iga,Vec gvec,Vec nvec);

PETSC_EXTERN PetscErrorCode IGAGetLocalVecArray(IGA iga,Vec gvec,Vec *lvec,const PetscScalar *array[]);
PETSC_EXTERN PetscErrorCode IGARestoreLocalVecArray(IGA iga,Vec gvec,Vec *lvec,const PetscScalar *array[]);

PETSC_EXTERN PetscErrorCode IGAClone(IGA iga,PetscInt dof,IGA *newiga);

#define DMIGA "iga"
PETSC_EXTERN PetscErrorCode DMIGACreate(IGA iga,DM *dm);
PETSC_EXTERN PetscErrorCode DMIGASetIGA(DM dm,IGA iga);
PETSC_EXTERN PetscErrorCode DMIGAGetIGA(DM dm,IGA *iga);

/* ---------------------------------------------------------------- */

struct _n_IGAElement {
  PetscInt refct;
  /**/
  PetscInt  start[3];
  PetscInt  width[3];
  PetscInt  sizes[3];
  PetscInt  ID[3];
  /**/
  PetscBool atboundary;
  PetscInt  boundary_id;
  /**/
  PetscInt count;
  PetscInt index;
  /**/
  PetscInt neq;
  PetscInt nen;
  PetscInt dof;
  PetscInt dim;
  PetscInt nsd;
  PetscInt npd;

  PetscInt    *mapping;   /*[nen]      */
  PetscInt    *rowmap;    /*[neq]      */
  PetscInt    *colmap;    /*[nen]      */
  PetscBool   rational;
  PetscReal   *rationalW; /*[nen]      */
  PetscBool   geometry;
  PetscReal   *geometryX; /*[nen][nsd] */
  PetscBool   property;
  PetscScalar *propertyA; /*[nen][npd] */

  PetscInt  nqp;
  PetscInt  sqp[3];

  PetscReal *weight;   /*   [nqp] */
  PetscReal *detJac;   /*   [nqp] */

  PetscReal *basis[5]; /*0: [nqp][nen] */
                       /*1: [nqp][nen][dim] */
                       /*2: [nqp][nen][dim][dim] */
                       /*3: [nqp][nen][dim][dim][dim] */
                       /*4: [nqp][nen][dim][dim][dim][dim] */

  PetscReal *shape[5]; /*0: [nqp][nen] */
                       /*1: [nqp][nen][nsd] */
                       /*2: [nqp][nen][nsd][nsd] */
                       /*3: [nqp][nen][nsd][nsd][nsd] */
                       /*4: [nqp][nen][nsd][nsd][nsd][nsd] */

  PetscReal *mapU[5];  /*0: [nqp][dim] */
                       /*1: [nqp][dim][nsd] */
                       /*2: [nqp][dim][nsd][nsd] */
                       /*3: [nqp][dim][nsd][nsd][nsd] */
                       /*4: [nqp][dim][nsd][nsd][nsd][nsd] */

  PetscReal *mapX[5];  /*0: [nqp][nsd] */
                       /*1: [nqp][nsd][dim] */
                       /*2: [nqp][nsd][dim][dim] */
                       /*3: [nqp][nsd][dim][dim][dim] */
                       /*4: [nqp][nsd][dim][dim][dim][dim] */

  PetscReal *normal;   /*   [nqp][nsd] */
  PetscReal *detX;     /*   [nqp] */
  PetscReal *detS;     /*   [nqp] */

  IGA      parent;
  IGAPoint iterator;

  PetscBool   collocation;

  PetscInt     nfix;
  PetscInt    *ifix;
  PetscScalar *vfix;
  PetscScalar *ufix;

  PetscInt     nflux;
  PetscInt    *iflux;
  PetscScalar *vflux;

  PetscInt    nval;
  PetscScalar *wval[8];
  PetscInt    nvec;
  PetscScalar *wvec[8];
  PetscInt    nmat;
  PetscScalar *wmat[4];

};

PETSC_EXTERN PetscErrorCode IGAElementCreate(IGAElement *element);
PETSC_EXTERN PetscErrorCode IGAElementDestroy(IGAElement *element);
PETSC_EXTERN PetscErrorCode IGAElementReference(IGAElement element);
PETSC_EXTERN PetscErrorCode IGAElementReset(IGAElement element);
PETSC_EXTERN PetscErrorCode IGAElementInit(IGAElement element,IGA iga);

PETSC_EXTERN PetscErrorCode IGAGetElement(IGA iga,IGAElement *element);
PETSC_EXTERN PetscErrorCode IGABeginElement(IGA iga,IGAElement *element);
PETSC_EXTERN PetscBool      IGANextElement(IGA iga,IGAElement element);
PETSC_EXTERN PetscErrorCode IGAEndElement(IGA iga,IGAElement *element);
PETSC_EXTERN PetscBool      IGAElementNextForm(IGAElement element,PetscBool visit[][2]);
PETSC_EXTERN PetscErrorCode IGAElementGetPoint(IGAElement element,IGAPoint *point);
PETSC_EXTERN PetscErrorCode IGAElementBeginPoint(IGAElement element,IGAPoint *point);
PETSC_EXTERN PetscBool      IGAElementNextPoint(IGAElement element,IGAPoint point);
PETSC_EXTERN PetscErrorCode IGAElementEndPoint(IGAElement element,IGAPoint *point);

PETSC_EXTERN PetscErrorCode IGAElementBuildClosure(IGAElement element);
PETSC_EXTERN PetscErrorCode IGAElementBuildTabulation(IGAElement element);

PETSC_EXTERN PetscErrorCode IGAElementGetParent(IGAElement element,IGA *parent);
PETSC_EXTERN PetscErrorCode IGAElementGetIndex(IGAElement element,PetscInt *index);
PETSC_EXTERN PetscErrorCode IGAElementGetCount(IGAElement element,PetscInt *count);
PETSC_EXTERN PetscErrorCode IGAElementGetSizes(IGAElement element,PetscInt *neq,PetscInt *nen,PetscInt *dof);
PETSC_EXTERN PetscErrorCode IGAElementGetClosure(IGAElement element,PetscInt *nen,const PetscInt *mapping[]);
PETSC_EXTERN PetscErrorCode IGAElementGetIndices(IGAElement element,
                                                 PetscInt *neq,const PetscInt *rowmap[],
                                                 PetscInt *nen,const PetscInt *colmap[]);

PETSC_EXTERN PetscErrorCode IGAElementGetWorkVec(IGAElement element,PetscScalar *V[]);
PETSC_EXTERN PetscErrorCode IGAElementGetWorkMat(IGAElement element,PetscScalar *M[]);
PETSC_EXTERN PetscErrorCode IGAElementGetValues(IGAElement element,const PetscScalar arrayU[],PetscScalar *U[]);

PETSC_EXTERN PetscErrorCode IGAElementBuildFix(IGAElement element);
PETSC_EXTERN PetscErrorCode IGAElementDelValues(IGAElement element,PetscScalar V[]);
PETSC_EXTERN PetscErrorCode IGAElementFixValues(IGAElement element,PetscScalar U[]);
PETSC_EXTERN PetscErrorCode IGAElementFixSystem(IGAElement element,PetscScalar K[],PetscScalar F[]);
PETSC_EXTERN PetscErrorCode IGAElementFixFunction(IGAElement element,PetscScalar F[]);
PETSC_EXTERN PetscErrorCode IGAElementFixJacobian(IGAElement element,PetscScalar J[]);

PETSC_EXTERN PetscErrorCode IGAElementAssembleVec(IGAElement element,const PetscScalar F[],Vec vec);
PETSC_EXTERN PetscErrorCode IGAElementAssembleMat(IGAElement element,const PetscScalar K[],Mat mat);

/* ---------------------------------------------------------------- */

struct _n_IGAPoint {
  PetscInt refct;
  /**/
  PetscBool atboundary;
  PetscInt  boundary_id;
  /**/
  PetscInt count;
  PetscInt index;
  /**/
  PetscInt neq;
  PetscInt nen;
  PetscInt dof;
  PetscInt dim;
  PetscInt nsd;
  PetscInt npd;

  PetscReal   *rational;/*  [nen] */
  PetscReal   *geometry;/*  [nen][nsd] */
  PetscScalar *property;/*  [nen][npd] */

  PetscReal *weight;   /*   [1] */
  PetscReal *detJac;   /*   [1] */

  PetscReal *point;    /*   [dim] */
  PetscReal *normal;   /*   [nsd] */

  PetscReal *basis[5]; /*0: [nen] */
                       /*1: [nen][dim] */
                       /*2: [nen][dim][dim] */
                       /*3: [nen][dim][dim][dim] */
                       /*4: [nen][dim][dim][dim][dim] */

  PetscReal *shape[5]; /*0: [nen]  */
                       /*1: [nen][nsd] */
                       /*2: [nen][nsd][nsd] */
                       /*3: [nen][nsd][nsd][nsd] */
                       /*3: [nen][nsd][nsd][nsd][nsd] */

  PetscReal *mapU[5];  /*0: [dim] */
                       /*1: [dim][nsd] */
                       /*2: [dim][nsd][nsd] */
                       /*3: [dim][nsd][nsd][nsd] */
                       /*4: [dim][nsd][nsd][nsd][nsd] */

  PetscReal *mapX[5];  /*0: [nsd] */
                       /*1: [nsd][dim] */
                       /*2: [nsd][dim][dim] */
                       /*3: [nsd][dim][dim][dim] */
                       /*4: [nsd][dim][dim][dim][dim] */

  PetscReal *detX;     /*   [1] */
  PetscReal *detS;     /*   [1] */

  IGAElement parent;

  PetscInt    nvec;
  PetscScalar *wvec[8];
  PetscInt    nmat;
  PetscScalar *wmat[4];
};

PETSC_EXTERN PetscErrorCode IGAPointCreate(IGAPoint *point);
PETSC_EXTERN PetscErrorCode IGAPointDestroy(IGAPoint *point);
PETSC_EXTERN PetscErrorCode IGAPointReference(IGAPoint point);
PETSC_EXTERN PetscErrorCode IGAPointReset(IGAPoint point);
PETSC_EXTERN PetscErrorCode IGAPointInit(IGAPoint point,IGAElement element);

PETSC_EXTERN PetscErrorCode IGAPointGetParent(IGAPoint point,IGAElement *element);
PETSC_EXTERN PetscErrorCode IGAPointGetIndex(IGAPoint point,PetscInt *index);
PETSC_EXTERN PetscErrorCode IGAPointGetCount(IGAPoint point,PetscInt *count);
PETSC_EXTERN PetscErrorCode IGAPointAtBoundary(IGAPoint point,PetscInt *axis,PetscInt *side);
PETSC_EXTERN PetscErrorCode IGAPointGetSizes(IGAPoint point,PetscInt *neq,PetscInt *nen,PetscInt *dof);
PETSC_EXTERN PetscErrorCode IGAPointGetDims(IGAPoint point,PetscInt *dim,PetscInt *nsd,PetscInt *npd);
PETSC_EXTERN PetscErrorCode IGAPointGetQuadrature(IGAPoint point,PetscReal *weigth,PetscReal *detJac);
PETSC_EXTERN PetscErrorCode IGAPointGetBasisFuns(IGAPoint point,PetscInt der,const PetscReal *basisfuns[]);
PETSC_EXTERN PetscErrorCode IGAPointGetShapeFuns(IGAPoint point,PetscInt der,const PetscReal *shapefuns[]);

PETSC_EXTERN PetscErrorCode IGAPointFormPoint(IGAPoint p,PetscReal x[]);
PETSC_EXTERN PetscErrorCode IGAPointFormGeomMap(IGAPoint p,PetscReal x[]);
PETSC_EXTERN PetscErrorCode IGAPointFormGradGeomMap(IGAPoint p,PetscReal F[]);
PETSC_EXTERN PetscErrorCode IGAPointFormInvGradGeomMap(IGAPoint p,PetscReal G[]);

PETSC_EXTERN PetscErrorCode IGAPointEvaluate (IGAPoint p,PetscInt ider,const PetscScalar U[],PetscScalar u[]);
PETSC_EXTERN PetscErrorCode IGAPointFormValue(IGAPoint p,const PetscScalar U[],PetscScalar u[]);
PETSC_EXTERN PetscErrorCode IGAPointFormGrad (IGAPoint p,const PetscScalar U[],PetscScalar u[]);
PETSC_EXTERN PetscErrorCode IGAPointFormHess (IGAPoint p,const PetscScalar U[],PetscScalar u[]);
PETSC_EXTERN PetscErrorCode IGAPointFormDel2 (IGAPoint p,const PetscScalar U[],PetscScalar u[]);
PETSC_EXTERN PetscErrorCode IGAPointFormDer3 (IGAPoint p,const PetscScalar U[],PetscScalar u[]);
PETSC_EXTERN PetscErrorCode IGAPointFormDer4 (IGAPoint p,const PetscScalar U[],PetscScalar u[]);

PETSC_EXTERN PetscErrorCode IGAPointGetWorkVec(IGAPoint point,PetscScalar *V[]);
PETSC_EXTERN PetscErrorCode IGAPointGetWorkMat(IGAPoint point,PetscScalar *M[]);
PETSC_EXTERN PetscErrorCode IGAPointAddArray(IGAPoint point,PetscInt n,const PetscScalar a[],PetscScalar A[]);
PETSC_EXTERN PetscErrorCode IGAPointAddVec(IGAPoint point,const PetscScalar f[],PetscScalar F[]);
PETSC_EXTERN PetscErrorCode IGAPointAddMat(IGAPoint point,const PetscScalar k[],PetscScalar K[]);

/* ---------------------------------------------------------------- */

struct _n_IGAProbe {
  PetscInt refct;
  /**/
  IGA iga;
  Vec gvec;
  Vec lvec;
  /**/
  PetscInt  order;
  PetscBool collective;
  PetscBool offprocess;
  /**/
  PetscInt  dim;
  PetscInt  nsd;
  PetscInt  dof;
  PetscInt  p[3];
  PetscReal *U[3];
  PetscInt  n[3];
  PetscInt  s[3];
  PetscInt  w[3];
  const PetscReal   *arrayW;
  const PetscReal   *arrayX;
  const PetscScalar *arrayA;
  /**/
  PetscInt    nen;
  PetscInt    *map;
  PetscReal   *W;
  PetscReal   *X;
  PetscScalar *A;
  /**/
  PetscReal point[3];
  PetscInt  ID[3];
  PetscReal *BD[3];
  PetscReal *basis[5]; /*0: [nen] */
                       /*1: [nen][dim] */
                       /*2: [nen][dim][dim] */
                       /*3: [nen][dim][dim][dim] */
                       /*4: [nen][dim][dim][dim][dim] */
  PetscReal *shape[5]; /*0: [nen] */
                       /*1: [nen][nsd] */
                       /*2: [nen][nsd][nsd] */
                       /*3: [nen][nsd][nsd][nsd] */
                       /*4: [nen][nsd][nsd][nsd][nsd] */
  PetscReal *mapU[5];  /*0: [dim] */
                       /*1: [dim][nsd] */
                       /*2: [dim][nsd][nsd] */
                       /*3: [dim][nsd][nsd][nsd] */
                       /*4: [dim][nsd][nsd][nsd][nsd] */
  PetscReal *mapX[5];  /*0: [nsd] */
                       /*1: [nsd][dim] */
                       /*2: [nsd][dim][dim] */
                       /*3: [nsd][dim][dim][dim] */
                       /*4: [nsd][dim][dim][dim][dim] */
  PetscReal *detX;
};

PETSC_EXTERN PetscErrorCode IGAProbeCreate(IGA iga,Vec A,IGAProbe *prb);
PETSC_EXTERN PetscErrorCode IGAProbeDestroy(IGAProbe *prb);
PETSC_EXTERN PetscErrorCode IGAProbeReference(IGAProbe prb);

PETSC_EXTERN PetscErrorCode IGAProbeSetOrder(IGAProbe prb,PetscInt order);
PETSC_EXTERN PetscErrorCode IGAProbeSetCollective(IGAProbe prb,PetscBool collective);

PETSC_EXTERN PetscErrorCode IGAProbeSetVec(IGAProbe prb,Vec A);
PETSC_EXTERN PetscErrorCode IGAProbeSetPoint(IGAProbe prb,const PetscReal u[]);

PETSC_EXTERN PetscErrorCode IGAProbeGeomMap  (IGAProbe prb,PetscReal x[]);
PETSC_EXTERN PetscErrorCode IGAProbeEvaluate (IGAProbe prb,PetscInt der,PetscScalar A[]);
PETSC_EXTERN PetscErrorCode IGAProbeFormValue(IGAProbe prb,PetscScalar A[]);
PETSC_EXTERN PetscErrorCode IGAProbeFormGrad (IGAProbe prb,PetscScalar A[]);
PETSC_EXTERN PetscErrorCode IGAProbeFormHess (IGAProbe prb,PetscScalar A[]);
PETSC_EXTERN PetscErrorCode IGAProbeFormDer3 (IGAProbe prb,PetscScalar A[]);
PETSC_EXTERN PetscErrorCode IGAProbeFormDer4 (IGAProbe prb,PetscScalar A[]);

/* ---------------------------------------------------------------- */

PETSC_EXTERN PetscLogEvent IGA_FormScalar;
PETSC_EXTERN PetscLogEvent IGA_FormVector;
PETSC_EXTERN PetscLogEvent IGA_FormMatrix;
PETSC_EXTERN PetscLogEvent IGA_FormSystem;
PETSC_EXTERN PetscLogEvent IGA_FormFunction;
PETSC_EXTERN PetscLogEvent IGA_FormJacobian;
PETSC_EXTERN PetscLogEvent IGA_FormIFunction;
PETSC_EXTERN PetscLogEvent IGA_FormIJacobian;

PETSC_EXTERN PetscErrorCode IGAComputeErrorNorm(IGA iga,PetscInt k,
                                                Vec vecU,IGAFormExact Exact,
                                                PetscReal enorm[],void *ctx);

PETSC_EXTERN PetscErrorCode IGAComputeScalar(IGA iga,Vec U,
                                             PetscInt n,PetscScalar S[],
                                             IGAFormScalar Scalar,void *ctx);

#define PCIGAEBE "igaebe"
#define PCIGABBB "igabbb"

PETSC_EXTERN PetscErrorCode IGAComputeVector(IGA iga,Vec B);
PETSC_EXTERN PetscErrorCode IGAComputeMatrix(IGA iga,Mat A);
PETSC_EXTERN PetscErrorCode IGAComputeSystem(IGA iga,Mat A,Vec B);

PETSC_EXTERN PetscErrorCode IGAComputeFunction(IGA iga,Vec U,Vec F);
PETSC_EXTERN PetscErrorCode IGAComputeJacobian(IGA iga,Vec U,Mat J);

PETSC_EXTERN PetscErrorCode IGAComputeIFunction(IGA iga,
                                                PetscReal a,Vec V,
                                                PetscReal t,Vec U,
                                                Vec F);
PETSC_EXTERN PetscErrorCode IGAComputeIJacobian(IGA iga,
                                                PetscReal a,Vec V,
                                                PetscReal t,Vec U,
                                                Mat J);
PETSC_EXTERN PetscErrorCode IGAComputeIEFunction(IGA iga,
                                                 PetscReal a, Vec V,
                                                 PetscReal t, Vec U,
                                                 PetscReal t0,Vec U0,
                                                 Vec F);
PETSC_EXTERN PetscErrorCode IGAComputeIEJacobian(IGA iga,
                                                 PetscReal a, Vec V,
                                                 PetscReal t, Vec U,
                                                 PetscReal t0,Vec U0,
                                                 Mat J);
PETSC_EXTERN PetscErrorCode IGAComputeRHSFunction(IGA iga,
                                                  PetscReal t,Vec U,
                                                  Vec F);
PETSC_EXTERN PetscErrorCode IGAComputeRHSJacobian(IGA iga,
                                                  PetscReal t,Vec U,
                                                  Mat J);
PETSC_EXTERN PetscErrorCode IGAComputeI2Function(IGA iga,
                                                 PetscReal a,Vec vecA,
                                                 PetscReal v,Vec vecV,
                                                 PetscReal t,Vec vecU,
                                                 Vec vecF);
PETSC_EXTERN PetscErrorCode IGAComputeI2Jacobian(IGA iga,
                                                 PetscReal a,Vec vecA,
                                                 PetscReal v,Vec vecV,
                                                 PetscReal t,Vec vecU,
                                                 Mat matJ);

PETSC_EXTERN PetscErrorCode KSPSetIGA(KSP ksp,IGA iga);
PETSC_EXTERN PetscErrorCode SNESSetIGA(SNES snes,IGA iga);
PETSC_EXTERN PetscErrorCode TSSetIGA(TS ts,IGA iga);

PETSC_EXTERN PetscErrorCode IGACreateKSP(IGA iga,KSP *ksp);
PETSC_EXTERN PetscErrorCode IGACreateSNES(IGA iga,SNES *snes);
PETSC_EXTERN PetscErrorCode IGACreateTS(IGA iga,TS *ts);
PETSC_EXTERN PetscErrorCode IGACreateTS2(IGA iga, TS *ts);

/* ---------------------------------------------------------------- */

PETSC_EXTERN PetscErrorCode IGAKSPFormRHS(KSP,Vec,void*);
PETSC_EXTERN PetscErrorCode IGAKSPFormOperators(KSP,Mat,Mat,void*);
PETSC_EXTERN PetscErrorCode IGASNESFormFunction(SNES,Vec,Vec,void*);
PETSC_EXTERN PetscErrorCode IGASNESFormJacobian(SNES,Vec,Mat,Mat,void*);
PETSC_EXTERN PetscErrorCode IGATSFormIFunction(TS,PetscReal,Vec,Vec,Vec,void*);
PETSC_EXTERN PetscErrorCode IGATSFormIJacobian(TS,PetscReal,Vec,Vec,PetscReal,Mat,Mat,void*);
PETSC_EXTERN PetscErrorCode IGATSFormI2Function(TS,PetscReal,Vec,Vec,Vec,Vec,void*);
PETSC_EXTERN PetscErrorCode IGATSFormI2Jacobian(TS,PetscReal,Vec,Vec,Vec,PetscReal,PetscReal,Mat,Mat,void*);

PETSC_EXTERN PetscErrorCode IGAPreparePCMG(IGA iga,PC pc);
PETSC_EXTERN PetscErrorCode IGAPreparePCBDDC(IGA iga,PC pc);
PETSC_EXTERN PetscErrorCode IGAPreparePCH2OPUS(IGA iga,PC pc);
PETSC_EXTERN PetscErrorCode IGASetOptionsHandlerPC(PC pc);
PETSC_EXTERN PetscErrorCode IGASetOptionsHandlerKSP(KSP ksp);
PETSC_EXTERN PetscErrorCode IGASetOptionsHandlerSNES(SNES snes);
PETSC_EXTERN PetscErrorCode IGASetOptionsHandlerTS(TS ts);

/* ---------------------------------------------------------------- */

#if defined(PETSC_USE_DEBUG)
#define IGACheckSetUp(iga,arg) do {                                     \
    if (PetscUnlikely(!(iga)->setup))                                   \
      SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,                 \
              "Must call IGASetUp() on argument %D \"%s\" before %s()", \
              (arg),#iga,PETSC_FUNCTION_NAME);                          \
    } while (0)
#define IGACheckSetUpStage(iga,arg,stg) do {                            \
    if (PetscUnlikely((iga)->setupstage<(stg))) IGACheckSetUp(iga,arg); \
    } while (0)
#else
#define IGACheckSetUp(iga,arg)          do {} while (0)
#define IGACheckSetUpStage(iga,arg,stg) do {} while (0)
#endif
#define IGACheckSetUpStage1(iga,arg) IGACheckSetUpStage(iga,arg,1)
#define IGACheckSetUpStage2(iga,arg) IGACheckSetUpStage(iga,arg,2)

#if defined(PETSC_USE_DEBUG)
#define IGACheckFormOp(iga,arg,FormOp) do {            \
    if (!iga->form->ops->FormOp)                       \
      SETERRQ(((PetscObject)iga)->comm,PETSC_ERR_USER, \
              "Must call IGASetForm%s() "              \
              "on argument %D \"%s\" before %s()",     \
              #FormOp,(arg),#iga,PETSC_FUNCTION_NAME); \
    } while (0)
#else
#define IGACheckFormOp(iga,arg,FormOp) do {} while (0)
#endif

/* ---------------------------------------------------------------- */

PETSC_EXTERN PetscErrorCode IGAOptionsAlias(const char name[],const char defval[],const char alias[]);
PETSC_EXTERN PetscErrorCode IGAOptionsDefault(const char prefix[],const char name[],const char value[]);
PETSC_EXTERN PetscErrorCode IGAOptionsReject(const char prefix[],const char name[]);

PETSC_EXTERN PetscEnum      IGAGetOptEnum(const char prefix[],const char name[],const char *const elist[],PetscEnum defval);
PETSC_EXTERN const char*    IGAGetOptString(const char prefix[],const char name[],const char defval[]);
PETSC_EXTERN PetscBool      IGAGetOptBool(const char prefix[],const char name[],PetscBool defval);
PETSC_EXTERN PetscInt       IGAGetOptInt(const char prefix[],const char name[],PetscInt defval);
PETSC_EXTERN PetscReal      IGAGetOptReal(const char prefix[],const char name[],PetscReal defval);
PETSC_EXTERN PetscScalar    IGAGetOptScalar(const char prefix[],const char name[],PetscScalar defval);

/* ---------------------------------------------------------------- */

#if PETSC_VERSION_LT(3,8,0)
PETSC_EXTERN PetscErrorCode TSSetMaxSteps(TS,PetscInt);
PETSC_EXTERN PetscErrorCode TSGetMaxSteps(TS,PetscInt*);
PETSC_EXTERN PetscErrorCode TSSetMaxTime(TS,PetscReal);
PETSC_EXTERN PetscErrorCode TSGetMaxTime(TS,PetscReal*);
#endif
#define TSALPHA1 TSALPHA
PETSC_EXTERN PetscErrorCode TSAlphaUseAdapt(TS,PetscBool);
#else

/* ---------------------------------------------------------------- */

#endif/*PETIGA_H*/

max_stars_repo_name: dalcinl/PetIGA

max_issues_repo_name: dalcinl/PetIGA

max_forks_repo_name: dalcinl/PetIGA


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

#ifndef __GSL_VECTOR_INT_H__
#define __GSL_VECTOR_INT_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 <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_inline.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_block_int.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

typedef struct
{
  size_t size;
  size_t stride;
  int *data;
  gsl_block_int *block;
  int owner;
}
gsl_vector_int;

typedef struct
{
  gsl_vector_int vector;
} _gsl_vector_int_view;

typedef _gsl_vector_int_view gsl_vector_int_view;

typedef struct
{
  gsl_vector_int vector;
} _gsl_vector_int_const_view;

typedef const _gsl_vector_int_const_view gsl_vector_int_const_view;

/* Allocation */

GSL_FUN gsl_vector_int *gsl_vector_int_alloc (const size_t n);
GSL_FUN gsl_vector_int *gsl_vector_int_calloc (const size_t n);

GSL_FUN gsl_vector_int *gsl_vector_int_alloc_from_block (gsl_block_int * b,
                                                     const size_t offset,
                                                     const size_t n,
                                                     const size_t stride);

GSL_FUN gsl_vector_int *gsl_vector_int_alloc_from_vector (gsl_vector_int * v,
                                                      const size_t offset,
                                                      const size_t n,
                                                      const size_t stride);

GSL_FUN void gsl_vector_int_free (gsl_vector_int * v);

/* Views */

GSL_FUN _gsl_vector_int_view
gsl_vector_int_view_array (int *v, size_t n);

GSL_FUN _gsl_vector_int_view
gsl_vector_int_view_array_with_stride (int *base,
                                         size_t stride,
                                         size_t n);

GSL_FUN _gsl_vector_int_const_view
gsl_vector_int_const_view_array (const int *v, size_t n);

GSL_FUN _gsl_vector_int_const_view
gsl_vector_int_const_view_array_with_stride (const int *base,
                                               size_t stride,
                                               size_t n);

GSL_FUN _gsl_vector_int_view
gsl_vector_int_subvector (gsl_vector_int *v,
                            size_t i,
                            size_t n);

GSL_FUN _gsl_vector_int_view
gsl_vector_int_subvector_with_stride (gsl_vector_int *v,
                                        size_t i,
                                        size_t stride,
                                        size_t n);

GSL_FUN _gsl_vector_int_const_view
gsl_vector_int_const_subvector (const gsl_vector_int *v,
                                  size_t i,
                                  size_t n);

GSL_FUN _gsl_vector_int_const_view
gsl_vector_int_const_subvector_with_stride (const gsl_vector_int *v,
                                              size_t i,
                                              size_t stride,
                                              size_t n);

/* Operations */

GSL_FUN void gsl_vector_int_set_zero (gsl_vector_int * v);
GSL_FUN void gsl_vector_int_set_all (gsl_vector_int * v, int x);
GSL_FUN int gsl_vector_int_set_basis (gsl_vector_int * v, size_t i);

GSL_FUN int gsl_vector_int_fread (FILE * stream, gsl_vector_int * v);
GSL_FUN int gsl_vector_int_fwrite (FILE * stream, const gsl_vector_int * v);
GSL_FUN int gsl_vector_int_fscanf (FILE * stream, gsl_vector_int * v);
GSL_FUN int gsl_vector_int_fprintf (FILE * stream, const gsl_vector_int * v,
                              const char *format);

GSL_FUN int gsl_vector_int_memcpy (gsl_vector_int * dest, const gsl_vector_int * src);

GSL_FUN int gsl_vector_int_reverse (gsl_vector_int * v);

GSL_FUN int gsl_vector_int_swap (gsl_vector_int * v, gsl_vector_int * w);
GSL_FUN int gsl_vector_int_swap_elements (gsl_vector_int * v, const size_t i, const size_t j);

GSL_FUN int gsl_vector_int_max (const gsl_vector_int * v);
GSL_FUN int gsl_vector_int_min (const gsl_vector_int * v);
GSL_FUN void gsl_vector_int_minmax (const gsl_vector_int * v, int * min_out, int * max_out);

GSL_FUN size_t gsl_vector_int_max_index (const gsl_vector_int * v);
GSL_FUN size_t gsl_vector_int_min_index (const gsl_vector_int * v);
GSL_FUN void gsl_vector_int_minmax_index (const gsl_vector_int * v, size_t * imin, size_t * imax);

GSL_FUN int gsl_vector_int_add (gsl_vector_int * a, const gsl_vector_int * b);
GSL_FUN int gsl_vector_int_sub (gsl_vector_int * a, const gsl_vector_int * b);
GSL_FUN int gsl_vector_int_mul (gsl_vector_int * a, const gsl_vector_int * b);
GSL_FUN int gsl_vector_int_div (gsl_vector_int * a, const gsl_vector_int * b);
GSL_FUN int gsl_vector_int_scale (gsl_vector_int * a, const int x);
GSL_FUN int gsl_vector_int_add_constant (gsl_vector_int * a, const double x);
GSL_FUN int gsl_vector_int_axpby (const int alpha, const gsl_vector_int * x, const int beta, gsl_vector_int * y);
GSL_FUN int gsl_vector_int_sum (const gsl_vector_int * a);

GSL_FUN int gsl_vector_int_equal (const gsl_vector_int * u,
                            const gsl_vector_int * v);

GSL_FUN int gsl_vector_int_isnull (const gsl_vector_int * v);
GSL_FUN int gsl_vector_int_ispos (const gsl_vector_int * v);
GSL_FUN int gsl_vector_int_isneg (const gsl_vector_int * v);
GSL_FUN int gsl_vector_int_isnonneg (const gsl_vector_int * v);

GSL_FUN INLINE_DECL int gsl_vector_int_get (const gsl_vector_int * v, const size_t i);
GSL_FUN INLINE_DECL void gsl_vector_int_set (gsl_vector_int * v, const size_t i, int x);
GSL_FUN INLINE_DECL int * gsl_vector_int_ptr (gsl_vector_int * v, const size_t i);
GSL_FUN INLINE_DECL const int * gsl_vector_int_const_ptr (const gsl_vector_int * v, const size_t i);

#ifdef HAVE_INLINE

INLINE_FUN
int
gsl_vector_int_get (const gsl_vector_int * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_VAL ("index out of range", GSL_EINVAL, 0);
    }
#endif
  return v->data[i * v->stride];
}

INLINE_FUN
void
gsl_vector_int_set (gsl_vector_int * v, const size_t i, int x)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_VOID ("index out of range", GSL_EINVAL);
    }
#endif
  v->data[i * v->stride] = x;
}

INLINE_FUN
int *
gsl_vector_int_ptr (gsl_vector_int * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (int *) (v->data + i * v->stride);
}

INLINE_FUN
const int *
gsl_vector_int_const_ptr (const gsl_vector_int * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (const int *) (v->data + i * v->stride);
}
#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_VECTOR_INT_H__ */

max_stars_repo_name: mgreter/astrometrylib

max_issues_repo_name: zzpwahaha/VimbaCamJILA

max_forks_repo_name: zzpwahaha/VimbaCamJILA


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include <complex.h>
#include <math.h>

#include <fftw3.h>

#include <gsl/gsl_sf_legendre.h>

#include "fsht.h"
#include "util.h"

int fshtinit (shdata *dat, int deg, int ntheta, int nphi, int measfft) {
	int ierr;
	complex double *fftbuf;

	/* This is actually one more than the maximum degree. */
	dat->deg = deg;

	/* The number of theta samples must at least equal the SH degree. */
	if (ntheta < deg) dat->ntheta = deg;
	else dat->ntheta = ntheta;

	/* At least this many phi values are required for fast FFT evaluation. */
	if (2 * deg - 1 > nphi) dat->nphi = 2 * deg - 1;
	else dat->nphi = nphi;

	/* Allocate the theta points and weights. */
	dat->theta = calloc (2 * dat->ntheta, sizeof(double));
	dat->weights = dat->theta + dat->ntheta;

	/* Find the Legendre-Gauss quadrature points. */
	gauleg (dat->ntheta, dat->theta, dat->weights);

	/* Temporarily allocate an FFT data buffer for planning. */
	fftbuf = fftw_malloc (dat->ntheta * dat->nphi * sizeof(complex double));

	/* Plan the forward and inverse transforms. */
	dat->fplan = fftw_plan_many_dft (1, &(dat->nphi), dat->ntheta, fftbuf,
			&(dat->nphi), 1, dat->nphi, fftbuf, &(dat->nphi), 1,
			dat->nphi, FFTW_FORWARD, measfft ? FFTW_MEASURE : FFTW_ESTIMATE);

	dat->bplan = fftw_plan_many_dft (1, &(dat->nphi), dat->ntheta, fftbuf,
			&(dat->nphi), 1, dat->nphi, fftbuf, &(dat->nphi), 1,
			dat->nphi, FFTW_BACKWARD, measfft ? FFTW_MEASURE : FFTW_ESTIMATE);

	/* The FFT buffer is no longer necessary, and will be reallocated
	 * on-the-fly when it is needed later. */
	fftw_free (fftbuf);

	return 0;
}

int fshtfree (shdata *dat) {
	if (dat->theta) free (dat->theta);

	fftw_destroy_plan (dat->fplan);
	fftw_destroy_plan (dat->bplan);

	dat->ntheta = dat->nphi = dat->deg = 0;

	return 0;
}

/* Scale the components of the spherical harmonic coefficients to relate the
 * far-field signature to the spherical harmonic expansion. The default (for
 * non-negative sgn) properly scales SH coefficients AFTER a forward
 * transform (angular to spherical). If sgn is negative, it scales the SH
 * coefficients BEFORE an inverse transform (spherical to angular). */
int shscale (complex double *samp, shdata *dat, int sgn) {
	complex double cscale[4] = { I, -1.0, -I, 1.0 };
	int i, j, off, idx;

	/* For negative sgn, flip the signs of the imaginary multipliers. */
	if (sgn < 0) {
		cscale[0] = -I;
		cscale[2] = I;
	}

	for (i = 0; i < dat->deg; ++i) {
		idx = i % 4;
		off = i * dat->nphi;

		/* Scale the zero order for all degrees. */
		samp[off] *= cscale[idx];

		/* Scale the nonzero orders for all degrees. */
		for (j = 1; j <= i; ++j) {
			samp[off + j] *= cscale[idx];
			samp[off + dat->nphi - j] *= cscale[idx];
		}
	}

	return 0;
}

/* Forward spherical harmonic transform: take samples of the function in theta
 * and phi into SH coefficients. */
int ffsht (complex double *samp, shdata *dat, int maxdeg) {
	int i, j, k, aoff, npk, dm1, deg = maxdeg;
	long lgn, lgi;
	double pc, scale, *lgvals;
	complex double *beta, *fftbuf;

	/* Copy the samples to the FFT buffer and perform the FFT. */
	fftbuf = fftw_malloc (dat->ntheta * dat->nphi * sizeof(complex double));
	memcpy (fftbuf, samp, dat->ntheta * dat->nphi * sizeof(complex double));

	/* Perform an in-place FFT of the buffer. */
	fftw_execute_dft (dat->fplan, fftbuf, fftbuf);

	/* Zero out the input samples, to prepare for storage of coefficients. */
	memset (samp, 0, dat->ntheta * dat->nphi * sizeof(complex double));

	beta = fftbuf;

	/* If no maximum degree is specified, use the default. */
	if (maxdeg < 1) deg = dat->deg;

	dm1 = deg - 1;

	/* Create storage for all Legendre polynomials. */
	lgn = gsl_sf_legendre_array_n (dm1);
	lgvals = malloc (lgn * sizeof(double));

	for (i = 0; i < dat->ntheta; ++i) {
		/* The scale factor that will be required. */
		scale = 2 * M_PI * dat->weights[i] / dat->nphi;

		/* Build the Legendre polynomials that we need.
		 * Don't use Condon-Shortley phase factor. */
		gsl_sf_legendre_array_e (GSL_SF_LEGENDRE_SPHARM, dm1, dat->theta[i], 1., lgvals);

		/* Handle m = 0 for all degrees. */
		for (j = 0; j < deg; ++j) {
			lgi = gsl_sf_legendre_array_index(j, 0);
			samp[j * dat->nphi] += scale * beta[0] * lgvals[lgi];
		}

		/* Handle nonzero orders for all relevant degrees. */
		for (k = 1; k < deg; ++k) {
			npk = dat->nphi - k;
			for (j = k; j < deg; ++j) {
				aoff = j * dat->nphi;
				lgi = gsl_sf_legendre_array_index(j, k);
				pc = scale * lgvals[lgi];
				/* The positive-order coefficient. */
				samp[aoff + k] += pc * beta[k];
				/* The negative-order coefficient. */
				samp[aoff + npk] += pc * beta[npk];
			}
		}

		beta += dat->nphi;
	}

	/* The FFT buffer is no longer required. */
	fftw_free (fftbuf);

	return deg;
}

/* Inverse spherical harmonic transform: take SH coefficients to sample of the
 * function in theta and phi. */
int ifsht (complex double *samp, shdata *dat, int maxdeg) {
	int i, j, k, aoff, npk, dm1, n, deg = maxdeg;
	long lgn, lgi;
	double *lgvals;
	complex double *beta, *fftbuf;

	/* The non-polar samples. */
	n = dat->ntheta * dat->nphi;

	/* Zero out the FFT buffer to prepare for evaluation of function. */
	fftbuf = fftw_malloc (n * sizeof(complex double));
	memset (fftbuf, 0, n * sizeof(complex double));

	beta = fftbuf;

	/* Use the default degree if no maximum is specified. */
	if (maxdeg < 1) deg = dat->deg;

	dm1 = deg - 1;

	/* Create storage for all Legendre polynomials. */
	lgn = gsl_sf_legendre_array_n(dm1);
	lgvals = malloc (lgn * sizeof(double));

	for (i = 0; i < dat->ntheta; ++i) {
		/* Build the Legendre polynomials that we need. */
		/* Don't use the Condon-Shortley phase factor. */
		gsl_sf_legendre_array_e (GSL_SF_LEGENDRE_SPHARM, dm1, dat->theta[i], 1., lgvals);

		for (j = 0; j < deg; ++j) {
			lgi = gsl_sf_legendre_array_index(j, 0);
			beta[0] += samp[j * dat->nphi] * lgvals[lgi];
		}

		for (k = 1; k < deg; ++k) {
			npk = dat->nphi - k;
			for (j = k; j < deg; ++j) {
				aoff = j * dat->nphi;
				lgi = gsl_sf_legendre_array_index(j, k);
				/* The positive-order coefficient. */
				beta[k] += samp[aoff + k] * lgvals[lgi];
				/* The negative-order coefficient. */
				beta[npk] += samp[aoff + npk] * lgvals[lgi];
			}
		}

		/* Move the next theta value in the FFT array. */
		beta += dat->nphi;
	}

	/* Perform the inverse FFT and copy the values to the storage area. */
	fftw_execute_dft (dat->bplan, fftbuf, fftbuf);
	memcpy (samp, fftbuf, n * sizeof(complex double));

	/* Eliminate the FFT buffer. */
	fftw_free (fftbuf);
	free (lgvals);

	return deg;
}

max_stars_repo_name: ahesford/fastsphere

max_issues_repo_name: ahesford/fastsphere

max_forks_repo_name: ahesford/fastsphere