libalberta-dev 3.0.1-1build1 / usr / include / alberta / alberta_util.h

/*******************************************************************************
 *
 * ALBERTA_UTIL:  tools for messages, memory allocation, parameters, etc.
 *
 * file: alberta_util.h
 *
 * description: public header file of the ALBERTA_UTIL package
 *
 *--------------------------------------------------------------------------
 *
 *  authors:   Alfred Schmidt
 *             Zentrum fuer Technomathematik
 *             Fachbereich 3 Mathematik/Informatik
 *             Universitaet Bremen
 *             Bibliothekstr. 2
 *             D-28359 Bremen, Germany
 *
 *             Kunibert G. Siebert
 *             Institut fuer Mathematik
 *             Universitaet Augsburg
 *             Universitaetsstr. 14
 *             D-86159 Augsburg, Germany
 *
 *
 *  www.alberta-fem.de
 *
 *  (c) by A. Schmidt and K.G. Siebert (1996-2003),
 *         C.-J. Heine (2002-2009)
 *
 ******************************************************************************/

#ifndef _ALBERTA_UTIL_H_
#define _ALBERTA_UTIL_H_

#include <sys/types.h>
#include <stdlib.h>
#include <stddef.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <float.h>
#include <stdarg.h>
#ifndef HAVE_OBSTACK_H
# define HAVE_OBSTACK_H 1
#endif
#if HAVE_OBSTACK_H
# include <obstack.h>
#else
/* Use our private copy */
# include <alberta/obstack.h>
#endif

#line 46 "../../alberta_util/src/alberta_util.h.in.in"

#ifdef __cplusplus
extern "C" {
#elif 0
}
#endif

#ifdef HAVE_STDBOOL_H
# include <stdbool.h>
# line 56 "../../alberta_util/src/alberta_util.h.in.in"
#else
# ifndef HAVE__BOOL
#  ifdef __cplusplus
typedef bool _Bool;
#  else
#   define _Bool signed char
#  endif
# endif
# ifndef __cplusplus
#  define bool _Bool
#  define false 0
#  define true 1
# else /* __cplusplus */
#  define bool bool
#  define false false
#  define true true
# endif /* __cplusplus */
# define __bool_true_false_are_defined 1
#endif

#if defined(__GNUC__)
# define __ATTRIBUTE_UNUSED__ __attribute__((unused))
# define __LIKELY__(arg)   __builtin_expect(arg, true)
# define __UNLIKELY__(arg) __builtin_expect(arg, false)
#else
# define __LIKELY__(arg)       (arg)
# define __UNLIKELY__(arg)     (arg)
# define __UNUSED_ATTRIBUTE__  /* nothing */
#endif

#if !ALBERTA_DEBUG\
 && defined(__GNUC__) && __GNUC__ >= 4 && __GNUC_MINOR__ >= 1\
 && !defined(__clang__) /* clang is a bloody liar ;) */

# define __FLATTEN_ATTRIBUTE__       __attribute__((flatten))
# define __FORCE_INLINE_ATTRIBUTE__  __attribute__((always_inline))
#else
# define __FLATTEN_ATTRIBUTE__      /* nothing */
# define __FORCE_INLINE_ATTRIBUTE__ /* nothing */
#endif


#ifndef CPP_CONCAT
/* Concatenation of preprocessor tokens; the recursion into the
 * underscored versions servers as evaluation loop.
 */
# define ___CPP_CONCAT(A, B) A##B
# define __CPP_CONCAT(A, B)  ___CPP_CONCAT(A, B)
# define _CPP_CONCAT(A, B)    __CPP_CONCAT(A, B)
# define CPP_CONCAT(A, B)    _CPP_CONCAT(A, B)
# define CPP_CONCAT3(A, B, C) CPP_CONCAT(A, CPP_CONCAT(B, C))
# define CPP_CONCAT4(A, B, C, D) CPP_CONCAT(A, CPP_CONCAT(B, CPP_CONCAT(C, D)))
#endif

#ifndef CPP_STRINGIZE
/* Turning of preprocessor tokens into string constants. */
# define ___CPP_STRINGIZE(A) #A
# define __CPP_STRINGIZE(A)  ___CPP_STRINGIZE(A)
# define _CPP_STRINGIZE(A)   __CPP_STRINGIZE(A)
# define CPP_STRINGIZE(A)    _CPP_STRINGIZE(A)
#endif

#ifndef nil
#define nil           NULL
#endif

#ifndef MAX
#define MAX(a, b)     ((a) > (b) ? (a) : (b))
#define MIN(a, b)     ((a) < (b) ? (a) : (b))
#endif
#ifndef ABS
#define ABS(a)        ((a) >= 0 ? (a) : -(a))
#endif

#ifndef SQR
#define SQR(a)        ((a)*(a))
#endif

#ifndef M_E
#define M_E		2.7182818284590452354
#endif
#ifndef M_LOG2E
#define M_LOG2E		1.4426950408889634074
#endif
#ifndef M_LOG10E
#define M_LOG10E	0.43429448190325182765
#endif
#ifndef M_LN2
#define M_LN2		0.69314718055994530942
#endif
#ifndef M_LN10
#define M_LN10		2.30258509299404568402
#endif
#ifndef M_PI
#define M_PI		3.14159265358979323846
#endif
#ifndef M_PI_2
#define M_PI_2		1.57079632679489661923
#endif
#ifndef M_PI_4
#define M_PI_4		0.78539816339744830962
#endif
#ifndef M_1_PI
#define M_1_PI		0.31830988618379067154
#endif
#ifndef M_2_PI
#define M_2_PI		0.63661977236758134308
#endif
#ifndef M_2_SQRTPI
#define M_2_SQRTPI	1.12837916709551257390
#endif
#ifndef M_SQRT2
#define M_SQRT2		1.41421356237309504880
#endif
#ifndef M_SQRT1_2
#define M_SQRT1_2	0.70710678118654752440
#endif

/*--------------------------------------------------------------------------*/
/*  Definition of basic data types      				    */
/*--------------------------------------------------------------------------*/

#ifndef REAL_DEFINED
#define REAL_DEFINED

#define USE_DOUBLE_AS_REAL  1

#if USE_DOUBLE_AS_REAL

typedef double              REAL;

#define REAL_MANT_DIG     DBL_MANT_DIG
#define REAL_EPSILON      DBL_EPSILON
#define REAL_DIG          DBL_DIG
#define REAL_MIN_EXP      DBL_MIN_EXP
#define REAL_MIN          DBL_MIN
#define REAL_MIN_10_EXP   DBL_MIN_10_EXP
#define REAL_MAX_EXP      DBL_MAX_EXP
#define REAL_MAX          DBL_MAX
#define REAL_MAX_10_EXP   DBL_MAX_10_EXP

#else  /*  USE_DOUBLE_AS_REAL  */

typedef float             REAL;

#define REAL_MANT_DIG     FLT_MANT_DIG
#define REAL_EPSILON      FLT_EPSILON
#define REAL_DIG          FLT_DIG
#define REAL_MIN_EXP      FLT_MIN_EXP
#define REAL_MIN          FLT_MIN
#define REAL_MIN_10_EXP   FLT_MIN_10_EXP
#define REAL_MAX_EXP      FLT_MAX_EXP
#define REAL_MAX          FLT_MAX
#define REAL_MAX_10_EXP   FLT_MAX_10_EXP

#endif  /*  USE_DOUBLE_AS_REAL  */

#ifndef LARGE
# define LARGE         REAL_MAX
#endif

#endif  /*  REAL_DEFINED  */

#ifndef U_CHAR_DEFINED
# define U_CHAR_DEFINED

typedef unsigned char       U_CHAR;
typedef signed char         S_CHAR;
#endif

#ifndef FLAGS_DEFINED
# define FLAGS_DEFINED

typedef unsigned long       FLAGS;
#endif

#ifndef BITS_64_DEFINED
# define BITS_64_DEFINED

typedef FLAGS BITS_64[64/(8*sizeof(FLAGS))];
#endif

#ifndef BITS_128_DEFINED
# define BITS_128_DEFINED

typedef FLAGS BITS_128[128/(8*sizeof(FLAGS))];
#endif

#ifndef BITS_256_DEFINED
# define BITS_256_DEFINED

typedef FLAGS BITS_256[256/(8*sizeof(FLAGS))];
#endif

#ifndef SCRATCH_MEM_DEFINED
# define SCATCH_MEM_DEFINED 1

typedef struct obstack SCRATCH_MEM[1];
typedef struct obstack *SCRATCH_MEM_PTR; /* A reference to an existing pool */

# line 249 "../../alberta_util/src/alberta_util.h.in.in"
# define obstack_chunk_alloc ALBERTA_OBSTACK_CHUNK_ALLOC
# define obstack_chunk_free  ALBERTA_OBSTACK_CHUNK_FREE

# define SCRATCH_MEM_INIT(handle) obstack_init((handle))
# define SCRATCH_MEM_ALLOC(handle, n_elem, type)			\
  (type *)obstack_alloc((handle), (n_elem)*sizeof(type))
# define SCRATCH_MEM_CALLOC(handle, n_elem, type)			\
  (type *)memset(SCRATCH_MEM_ALLOC(handle, n_elem, type),		\
		 0, (n_elem)*sizeof(type))
# define SCRATCH_MEM_ZAP(handle)			\
  {						\
    struct obstack _AI_obst = *(handle);	\
    obstack_free(&_AI_obst, NULL);		\
  }
# define SCRATCH_MEM_CPY(to, from) *(to) = *(from)
#endif /* SCRATCH_MEM_DEFINED */

#ifndef WORKSPACE_DEFINED
# define WORKSPACE_DEFINED
typedef struct workspace    WORKSPACE;

struct workspace
{
  size_t  size;
  void    *work;
};
#endif

#ifndef DBL_LIST_NODE_DEFINED
# define DBL_LIST_NODE_DEFINED
/* A simple, abstracted doubly linked list, see alberta_util_inlines.h
 * for inline functions and support macros.
 */
typedef struct dbl_list_node {
  struct dbl_list_node *next, *prev;
} DBL_LIST_NODE;
#endif

/*--------------------------------------------------------------------------*/
/*  functions and macros for messages	                         	    */
/*--------------------------------------------------------------------------*/

#ifdef __GNUC__
# define ALBERTA_UNUSED(proto) proto __attribute__((unused))
# define ALBERTA_DEFUNUSED(proto) ALBERTA_UNUSED(proto); proto
# define ALBERTA_NORETURN(proto) proto __attribute__((noreturn))
# define ALBERTA_DEFNORETURN(proto) ALBERTA_NORETURN(proto); proto
#else
# define ALBERTA_UNUSED(proto) proto
# define ALBERTA_DEFUNUSED(proto) proto
# define ALBERTA_NORETURN(proto) proto
# define ALBERTA_DEFNORETURN(proto) proto
#endif

/*--------------------------------------------------------------------------*/
/* Some of the printing functions are declared to return type int only for  */
/* the purpose of defining the macros below without "if" statements!        */
/*--------------------------------------------------------------------------*/

int  print_msg(const char *format, ...);
int  print_error_msg(const char *format, ...);
ALBERTA_NORETURN(void print_error_msg_exit(const char *format, ...));
int  print_error_funcname(const char *name, const char *file, int line);
void print_warn_msg(const char *format, ...);
void print_warn_funcname(const char *name, const char *file, int line);
int  print_funcname(const char *name);
void alberta_print_int_vec(const char *s, const int *vec, int no);
void alberta_print_real_vec(const char *s, const REAL *vec, int no);

void change_msg_out(FILE *fp);
void open_msg_file(const char *filename, const char *type);
void change_error_out(FILE *fp);
void open_error_file(const char *filename, const char *type);
const char *generate_filename(const char *, const char *, int);

#define FUNCNAME(nn) static const char *funcName __ATTRIBUTE_UNUSED__ = nn

#define ERROR(...)							\
  ( print_error_funcname(funcName ? funcName : __func__, __FILE__, __LINE__), \
    print_error_msg(__VA_ARGS__) )
#define ERROR_EXIT(...)							\
  ( print_error_funcname(funcName ? funcName : __func__, __FILE__, __LINE__), \
    print_error_msg_exit(__VA_ARGS__) )
#define WARNING(...)					\
  ( print_warn_funcname(funcName ? funcName : __func__, __FILE__, __LINE__), \
    print_warn_msg(__VA_ARGS__) )

#define TEST(test, ...)						\
  do {								\
    if (!(test)) {						\
      print_error_funcname(funcName				\
			   ? funcName : __func__, __FILE__, __LINE__);	\
      print_error_msg(__VA_ARGS__);				\
    }								\
  } while (0)
#define TEST_EXIT(test, ...)					\
  do {								\
    if (!(test)) {						\
      print_error_funcname(funcName				\
			   ? funcName : __func__, __FILE__, __LINE__);	\
      print_error_msg_exit(__VA_ARGS__);			\
    }								\
  } while (0)

/*--------------------------------------------------------------------------*/
/* Some testing macros useful while debugging.                              */
/*--------------------------------------------------------------------------*/

#if ALBERTA_DEBUG==1
#define DEBUG_TEST(test, ...)						\
  do {									\
    if (!(test)) {							\
      print_error_funcname(funcName					\
			   ? funcName : __func__, __FILE__, __LINE__);	\
      print_error_msg(__VA_ARGS__);					\
    }									\
  } while (0)
#define DEBUG_TEST_EXIT(test, ...)					\
  do {									\
    if (!(test)) {							\
      print_error_funcname(funcName					\
			   ? funcName : __func__, __FILE__, __LINE__);	\
      print_error_msg_exit(__VA_ARGS__);				\
    }									\
  } while (0)
#else
#define DEBUG_TEST(test, ...)      do { (void)funcName; } while (0)
#define DEBUG_TEST_EXIT(test, ...) do { (void)funcName; } while (0)
#endif


#define MSG print_funcname(funcName ? funcName : __func__), print_msg

#define INFO(info, noinfo, ...)						\
  do {									\
    if (msg_info&&(MIN(msg_info, (info))>=(noinfo))) {			\
      print_funcname(funcName ? funcName : __func__); print_msg(__VA_ARGS__); \
    }									\
  } while (0)

#define PRINT_INFO(info, noinfo, ...)			\
  do {							\
    if (msg_info&&(MIN(msg_info, (info))>=(noinfo))) {	\
      print_msg(__VA_ARGS__);				\
    }							\
  } while (0)

#define PRINT_INT_VEC						\
  print_funcname(funcName ? funcName : __func__), alberta_print_int_vec
#define PRINT_REAL_VEC							\
  print_funcname(funcName ? funcName : __func__), alberta_print_real_vec

extern void alberta_wait(const char *, int);
#define WAIT					\
  alberta_wait(funcName ? funcName : __func__, msg_wait)
#define WAIT_REALLY				\
  alberta_wait(funcName ? funcName : __func__, 1)

/*--------------------------------------------------------------------------*/
/*  Definition of extern variables      				    */
extern bool          msg_wait;                       /* declared in msg.c   */
extern const char    *funcName;                      /* declared in msg.c   */
extern int           msg_info;                       /* declared in msg.c   */

/*--------------------------------------------------------------------------*/
/*  Macros and proto-types for memory allocation			    */
/*--------------------------------------------------------------------------*/

void **alberta_matrix(int nr, int nc, size_t elsize,
		      const char *, const char *, int);
void free_alberta_matrix(void **ptr, int nr, int nc, size_t elsize);
void clear_alberta_matrix(void **ptr, int nr, int nc, size_t elsize);
void ***alberta_3array(int nr, int nc, int nd, size_t size,
		       const char *fct, const char *file, int line);
void free_alberta_3array(void ***array, int nr, int nc, int nd, size_t size);
void clear_alberta_3array(void ***array, int nr, int nc, int nd, size_t size);

void ****alberta_4array(int n0, int n1, int n2, int n3, size_t size,
			const char *fct, const char *file, int line);
void free_alberta_4array(void ****array,
			 int n0, int n1, int n2, int n3, size_t size);
void clear_alberta_4array(void ****array,
			  int n0, int n1, int n2, int n3, size_t size);

WORKSPACE *get_workspace(size_t, const char *, const char *, int);
#define GET_WORKSPACE(s)						\
  get_workspace(s, funcName ? funcName : __func__, __FILE__, __LINE__)

WORKSPACE  *realloc_workspace(WORKSPACE *ws, size_t new_size,
			      const char *fct, const char *file, int line);
static inline
WORKSPACE *__REALLOC_WORKSPACE(WORKSPACE *ws, size_t sz,
			       const char *fct, const char *file, int line)
{
  if (ws == NULL || ws->size < sz) {
    return realloc_workspace(
      ws, sz, fct, __FILE__, __LINE__);
  } else {
    return ws;
  }
}
#define REALLOC_WORKSPACE(ws, sz)					\
  __REALLOC_WORKSPACE((ws), (sz),					\
		      funcName ? funcName : __func__, __FILE__, __LINE__)

void clear_workspace(WORKSPACE *ws);
#define CLEAR_WORKSPACE(ws) clear_workspace(ws)
void free_workspace(WORKSPACE *ws);
#define FREE_WORKSPACE(ws) free_workspace(ws)

/* If the size of a block is not known then alberta_free() and
 * alberta_realloc() accept the "special" size ~0 to indicate
 * this. The only purpose of passing the size to realloc() and free()
 * is to keep track of the memory consumption such that
 * print_mem_use() prints the correct statistics. In principle we
 * should get rid of this nonsense; on any decent system memory
 * allocation statistics can be obtained by more reliable
 * methods. -- cH.
 */
#define MEM_UNKNOWN_SIZE ((size_t)~0)

#ifndef ALBERTA_EFENCE
# define ALBERTA_EFENCE       0
#endif
#ifndef ALBERTA_DUMA
# define ALBERTA_DUMA       0
#endif
#ifndef ALBERTA_ALLOC_RECORD
# define ALBERTA_ALLOC_RECORD 0
#endif
#ifndef ALBERTA_ALLOC_MALLOC
# define ALBERTA_ALLOC_MALLOC 0
#endif

#if !ALBERTA_DEBUG
# undef  ALBERTA_EFENCE
# define ALBERTA_EFENCE 0
# undef  ALBERTA_DUMA
# define ALBERTA_DUMA 0
# undef  ALBERTA_ALLOC_RECORD
# define ALBERTA_ALLOC_RECORD 0
#endif

#if ALBERTA_EFENCE || ALBERTA_DUMA || ALBERTA_ALLOC_RECORD
# undef ALBERTA_ALLOC_MALLOC
# define ALBERTA_ALLOC_MALLOC 1
#endif

#if ALBERTA_ALLOC_MALLOC
# define alberta_alloc(size, fct, file, line)	\
  ((void)(fct), (void)(file), (void)(line), malloc(size))
# define alberta_realloc(ptr, old_size, new_size, fct, file, line)	\
  ((void)(fct), (void)(file), (void)(line), realloc(ptr, new_size))
# define alberta_calloc(size, elsize, fct, file, line) \
  ((void)(fct), (void)(file), (void)(line), calloc(size, elsize))
# define alberta_free(ptr, size)                 \
  ((void)(size), free(ptr))
# define alberta_obstack_chunk_alloc                   malloc
# define alberta_obstack_chunk_free                    free
void print_mem_use(void);
#else
void *alberta_alloc(size_t size, const char *fct, const char *file, int line);
void *alberta_realloc(void *ptr, size_t old_size, size_t new_size,
		      const char *fct, const char *file, int line);
void *alberta_calloc(size_t size, size_t elsize, const char *fct,
		     const char *file, int line);
void alberta_free(void *ptr, size_t size);
void *alberta_obstack_chunk_alloc(size_t size);
void alberta_obstack_chunk_free(void *ptr);
void print_mem_use(void);
#endif

#if ALBERTA_ALLOC_RECORD
/* proto-types */
extern void *
malloc_record__(size_t size, const char *file, const char *func, int line);
extern void free_record__(void *ptr);
extern char *strdup_record__(const char *string,
			     const char *file, const char *func, int line);
extern void *
realloc_record__(void *ptr, size_t size,
		 const char *file, const char *func, int line);
extern void *
calloc_record__(size_t nmemb, size_t size,
		const char *file, const char *func, int line);
extern void print_alloc_records(int tail);

# undef alberta_alloc
# define alberta_alloc(size, fct, file, line)	\
  malloc_record__(size, file, fct, line)
# undef alberta_realloc
# define alberta_realloc(ptr, old_size, new_size, fct, file, line)	\
  realloc_record__(ptr, new_size, file, fct, line)
# undef alberta_calloc
# define alberta_calloc(size, elsize, fct, file, line)	\
  calloc_record__(size, elsize, file, fct, line)
# undef alberta_free
# define alberta_free(ptr, size)		\
  free_record__(ptr)
# undef print_mem_use
# define print_mem_use() print_alloc_records(-1)

/* The following two must be functions and obey the calling
 * conventions of the obstack implementation.
 */
# undef alberta_obstack_chunk_alloc
static inline void *alberta_obstack_chunk_alloc(size_t size)
{
  return malloc_record__(size, __func__, __FILE__, __LINE__);
}
# undef alberta_obstack_chunk_free
static inline void alberta_obstack_chunk_free(void *ptr)
{
  free_record__(ptr);
}

# undef free
# define free(ptr) free_record__((ptr))
# undef malloc
# define malloc(sz) malloc_record__((sz), __FILE__, __func__, __LINE__)
# undef calloc
# define calloc(n, sz) calloc_record__((n), (sz), __FILE__, __func__, __LINE__)
# undef realloc
# define realloc(op, sz)					\
  realloc_record__((op), (sz), __FILE__, __func__, __LINE__)
# undef strdup
# define strdup(s) strdup_record__(s, __FILE__, __func__, __LINE__)

#endif /* end of ALBERTA_ALLOC_RECORD stuff */

#if ALBERTA_EFENCE

# include <efence.h>
# line 575 "../../alberta_util/src/alberta_util.h.in.in"

/* There are just so many of them and efence does not define it */
# undef strdup
# define strdup(s)				\
  ({char *_strdup_s = malloc(strlen(s)+1);	\
    strcpy(_strdup_s, s);			\
    _strdup_s;					\
  })
#endif /* end of EFENCE stuff */

#if ALBERTA_DUMA

# include <duma.h>
# line 594 "../../alberta_util/src/alberta_util.h.in.in"

#endif /* end of DUMA stuff */

/* standard allocation */
#define MEM_ALLOC(n, type)						\
  ((type*)alberta_alloc((n)*sizeof(type),				\
			funcName ? funcName : __func__, __FILE__, __LINE__))
#define MEM_REALLOC(ptr, old_n, new_n, type)				\
  ((type*)alberta_realloc((void *)ptr, (old_n)*sizeof(type),		\
			  (new_n)*sizeof(type),				\
			  funcName ? funcName : __func__, __FILE__, __LINE__))
#define MEM_CALLOC(n, type)						\
  ((type*)alberta_calloc((n), sizeof(type),				\
			 funcName ? funcName : __func__, __FILE__, __LINE__))
#define MEM_FREE(ptr, n, type) alberta_free((void *)ptr, (n)*sizeof(type))
#define ALBERTA_OBSTACK_CHUNK_ALLOC alberta_obstack_chunk_alloc
#define ALBERTA_OBSTACK_CHUNK_FREE  alberta_obstack_chunk_free

/* Genereate a self destructing hanble. */
static inline SCRATCH_MEM_PTR SCRATCH_MEM_GET(void)
{
  SCRATCH_MEM_PTR scrptr;
  SCRATCH_MEM     scratch;

  SCRATCH_MEM_INIT(scratch);

  scrptr = (SCRATCH_MEM_PTR)SCRATCH_MEM_ALLOC(scratch, 1, SCRATCH_MEM);
  SCRATCH_MEM_CPY(scrptr, scratch);
  return scrptr;
}

#define VECMAT_ALLOC(nr, nc, stride, type)			\
  ((type**)alberta_matrix((nr), (nc), (stride)*sizeof(type),	\
			  funcName ? funcName : __func__, __FILE__, __LINE__))
#define VECMAT_FREE(ptr, nr, nc, stride, type)				\
  free_alberta_matrix((void **)ptr, (nr), (nc), (stride)*sizeof(type))
#define VECMAT_CLEAR(ptr, nr, nc, stride, type)				\
  clear_alberta_matrix((void **)ptr, (nr), (nc), (stride)*sizeof(type))
#define MAT_ALLOC(nr, nc, type)     VECMAT_ALLOC(nr, nc, 1, type)
#define MAT_FREE(ptr, nr, nc, type) VECMAT_FREE(ptr, nr, nc, 1, type)
#define MAT_CLEAR(ptr, nr, nc, type) VECMAT_CLEAR(ptr, nc, nc, 1, type)

#define ARRAY3_ALLOC(nr, nc, nd, type)					\
  ((type ***)alberta_3array((nr), (nc), (nd), sizeof(type),		\
			    funcName ? funcName : __func__, __FILE__, __LINE__))
#define ARRAY3_FREE(ptr, nr, nc, nd, type)				\
  free_alberta_3array((void ***)(ptr), (nr), (nc), (nd), sizeof(type))
#define ARRAY3_CLEAR(ptr, nr, nc, nd, type)	\
  clear_alberta_3array((void ***)(ptr), (nr), (nc), (nd), sizeof(type))

#define ARRAY4_ALLOC(n0, n1, n2, n3, type)				\
  ((type ****)alberta_4array((n0), (n1), (n2), (n3), sizeof(type),	\
			    funcName ? funcName : __func__, __FILE__, __LINE__))
#define ARRAY4_FREE(ptr, n0, n1, n2, n3, type)				\
  free_alberta_4array((void ****)(ptr), (n0), (n1), (n2), (n3), sizeof(type))
#define ARRAY4_CLEAR(ptr, n0, n1, n2, n3, type)				\
  clear_alberta_4array((void ****)(ptr), (n0), (n1), (n2), (n3), sizeof(type))

#define NAME(struct)							\
  ((struct) ? ((struct)->name ? (struct)->name : #struct "->name unknown") : \
   #struct " pointer to NULL")

#define GET_STRUCT(struct, from)					\
  TEST_EXIT(((from)  && (struct = (from)->struct)) || ((struct) = NULL), \
	    (from)							\
	    ? "GET_STRUCT: %s->%s == NULL\n"				\
	    : "GET_STRUCT: %s == NULL%s\n",				\
	    #from, (from) ? #struct : "")

/*--------------------------------------------------------------------------*/
/* functions for handling parameters                                        */
/*--------------------------------------------------------------------------*/

void add_parameter(int p, const char *key, const char *par);
void Add_parameter(int p, const char *key, const char *par,
		   const char *fname, const char *file, int line);
#define ADD_PARAMETER(p, key, param)				\
  Add_parameter(p, key, param,					\
		funcName ? funcName : __func__, __FILE__, __LINE__)

void init_parameters(int p, const char *filename);
void init_parameters_cpp(int p, const char *fn, const char *flags);

int get_parameter(int flag, const char *key, const char *format, ...);
int init_param_func_name(const char *, const char *, int call_line);

#define GET_PARAMETER							\
  (!init_param_func_name(funcName ? funcName : __func__, __FILE__, __LINE__) \
   ? 0 : get_parameter)
void save_parameters(const char *file, int info);
void print_parameters(void);

/*--------------------------------------------------------------------------*/
/*  generating filenames                                                    */
/*--------------------------------------------------------------------------*/

const char *generate_filename(const char *path, const char *fn, int ntime);

/* A poor man's SIGFPE exception catcher. If unmask == 1, then SIGFPE
 * is unmasked and caught, if possible. If unmask == 0 floating point
 * exception are not caught (this is the default for most systems).
 */
extern void sigfpe_init(bool unmask);


/*--------------------------------------------------------------------------*/
/*  library of Orthogonal Error Methods                                     */
/*                                                                          */
/*  most routines are the C-version of OFM-lib by Willy Doerfler            */
/*                                                                          */
/*  author: Kunibert G. Siebert                                             */
/*          Institut fuer Mathematik                                        */
/*          Universitaet Augsburg                                           */
/*          Universitaetsstr. 14                                            */
/*          D-86159 Augsburg, Germany                                       */
/*                                                                          */
/*          http://scicomp.math.uni-augsburg.de/Siebert/                    */
 /*                                                                          */
 /* (c) by K.G. Siebert (2000-2003)                                          */
 /*--------------------------------------------------------------------------*/

 /*--------------------------------------------------------------------------*/
 /*---  data structures for orthogonal error methods  -----------------------*/
 /*--------------------------------------------------------------------------*/

/*******************************************************************************
 *  solver identification for interface to OEM-lib
 ******************************************************************************/
typedef enum {
  NoSolver = 0,
  BiCGStab = 1,
  CG       = 2,
  GMRes    = 3,
  ODir     = 4,
  ORes     = 5,
  TfQMR    = 6,
  GMRes_k  = 7,
  SymmLQ   = 8
} OEM_SOLVER;

 /* Apply a matrix to vector, either directly or by solving the equation
  * M^{-1} u = rhs
  */
typedef int (*OEM_MV_FCT)(void *data, int dim, const REAL *rhs, REAL *u);

typedef struct oem_data OEM_DATA;
struct oem_data
{
  OEM_MV_FCT mat_vec;
  void       *mat_vec_data;
  OEM_MV_FCT mat_vec_T;
  void       *mat_vec_T_data;
  void       (*left_precon)(void *, int, REAL *);
  void       *left_precon_data;
  void       (*right_precon)(void *, int, REAL *);
  void       *right_precon_data;

  REAL       (*scp)(void *, int, const REAL *, const REAL *);
  void       *scp_data;

  WORKSPACE  *ws;

  REAL       tolerance;
  int        restart;
  int        max_iter;
  int        info;

  int        terminate_reason;

  REAL       initial_residual;
  REAL       residual;
};

/*--------------------------------------------------------------------------*/
/*--  implemented solvers  -------------------------------------------------*/
/*--------------------------------------------------------------------------*/

int oem_bicgstab(OEM_DATA *oem, int dim, const REAL *rhs, REAL *u);
int oem_cg(OEM_DATA *oem, int dim, const REAL *rhs, REAL *u);
int oem_gmres(OEM_DATA *oem, int dim, const REAL *rhs, REAL *u);
int oem_gmres_k(OEM_DATA *oem, int dim, const REAL *rhs, REAL *u);
int oem_odir(OEM_DATA *oem, int dim, const REAL *rhs, REAL *u);
int oem_ores(OEM_DATA *oem, int dim, const REAL *rhs, REAL *u);
int oem_tfqmr(OEM_DATA *oem, int dim, const REAL *rhs, REAL *u);
int oem_symmlq(OEM_DATA *oem, int dim, const REAL *rhs, REAL *u);

/*--------------------------------------------------------------------------*/
/*  function call:                                                          */
/*    oem_?(oem_data, dim, b, x);                                           */
/*                                                                          */
/*    input:  oem_data:  above described data structure                     */
/*            dim:       dimension of the linear system                     */
/*            b:         vector storing the right hand side of the sytem    */
/*            x:         vector stroing the initial guess                   */
/*                                                                          */
/*    output: x:         solution of the linear system                      */
/*            oem_data->residual:  residual of the equation                 */
/*                                                                          */
/*    return value:      number of iterations                               */
/*--------------------------------------------------------------------------*/

/*--------------------------------------------------------------------------*/
/*  Solver for linear saddle point problem                                  */
/*--------------------------------------------------------------------------*/

/* y += factor A x */
typedef void (*OEM_GEMV_FCT)(void *data,
			     REAL factor,
			     int dim_x, const REAL *x, int dim_y, REAL *y);

typedef struct oem_sp_data OEM_SP_DATA;
struct oem_sp_data
{
  OEM_MV_FCT   solve_Auf;
  void         *solve_Auf_data;  /*--  solve A u = f  in X  ------------------*/
  OEM_GEMV_FCT B;
  void         *B_data;          /*--  compute B p ---------------------------*/
  OEM_GEMV_FCT Bt;
  void         *Bt_data;         /*--  compute B^* u -------------------------*/
  OEM_MV_FCT   project;
  void         *project_data;    /*--  project B^* u to Y  -------------------*/
  int          (*precon)(void *ud,
			 int dimY, const REAL *g_Btu, const REAL *r, REAL *Cr);
  void         *precon_data;     /*--  compute C r in Y  ---------------------*/

  WORKSPACE    *ws;              /*--  workspace  ----------------------------*/

  REAL         tolerance;
  int          restart;
  int          max_iter;
  int          info;

  REAL         initial_residual;
  REAL         residual;
};

/*--------------------------------------------------------------------------*/
/*--  implemented solvers  -------------------------------------------------*/
/*--------------------------------------------------------------------------*/

int oem_spcg(OEM_SP_DATA *data, int dimX, const REAL *f, REAL *u, int dimY,
	     const REAL *g, REAL *p);

/*--------------------------------------------------------------------------*/
/*  library of solvers for nonlinear systems                                */
/*                                                                          */
/*  some routines are the C-version of NGL-lib by Willy Doerfler            */
/*     http://www.mathematik.uni-kl.de/~doerfler                            */
 /*                                                                          */
 /*  author: Kunibert G. Siebert                                             */
 /*          Institut f"ur Angewandte Mathematik                             */
 /*          Albert-Ludwigs-Universit"at Freiburg                            */
 /*          Hermann-Herder-Str. 10                                          */
 /*          79104 Freiburg                                                  */
 /*          Germany                                                         */
 /*                                                                          */
 /* email:  kunibert@mathematik.uni-freiburg.de                              */
 /*                                                                          */
 /* (c) by K.G. Siebert (2000)                                               */
 /*                                                                          */
 /*--------------------------------------------------------------------------*/

typedef struct nls_data NLS_DATA;
struct nls_data
{
  void       (*update)(void *update_data,
		       int dim, const REAL *x_k, bool ud_DF, REAL *Fx);
  void       *update_data;
  int        (*solve)(void *solve_data, int dim, const REAL *rhs, REAL *d);
  void       *solve_data;
  REAL       (*norm)(void *norm_data, int dim, const REAL *residual);
  void       *norm_data;

  WORKSPACE  *ws;

  REAL       tolerance;
  int        restart;
  int        max_iter;
  int        info;

  REAL       initial_residual;
  REAL       residual;
};

/*--------------------------------------------------------------------------*/
/*  functions:                                                              */
/*--------------------------------------------------------------------------*/
/*    update:  calculate F(x) and DF(x)                                     */
/*    call:                                                                 */
/*             update(update_data, dim, x, ud_DF, Fx)                       */
/*    arguments:                                                            */
/*             update_data: pointer to user data                            */
/*             dim:         dimension of the nonlinear system               */
/*             x:           vector storing the actual iterate               */
/*             ud_DF:       ud_DF==1: update DF(x)                          */
/*             Fx:          Fx!=NULL:  update F(x) and store it in Fx       */
/*--------------------------------------------------------------------------*/
/*    solve:   solve the linearized system DF(x) d = b for the correction   */
/*    call:                                                                 */
/*             solve(solve_data, dim, b, d)                                 */
/*             solve_data:  pointer to user data                            */
/*             dim:         dimension of the nonlinear system               */
/*             b:           right hand side of the linearized system        */
/*             d:           initial guess and solution of the system        */
/*                                                                          */
/*             return value is the number of iterations used by an          */
/*             iterative solver                                             */
/*--------------------------------------------------------------------------*/
/*    norm:    compute the norm of the residual                             */
/*             if norm==NULL, the Eucledian norm of the coefficient vector  */
/*             in R^dim is used.                                            */
/*    call:                                                                 */
/*             norm(norm_data, dim, d)                                      */
/*             norm_data:   pointer to user data                            */
/*             dim:         dimension of the nonlinear system               */
/*             d:           coefficient vector of the residual              */
/*                                                                          */
/*             return value is the norm of the residual                     */
/*--------------------------------------------------------------------------*/
/*  ws: workspace                                                           */
/*--------------------------------------------------------------------------*/

int nls_newton(NLS_DATA *, int, REAL *);
int nls_newton_ds(NLS_DATA *, int, REAL *);
int nls_newton_fs(NLS_DATA *, int, REAL *);
int nls_newton_br(NLS_DATA *, REAL, int, REAL *);

int nls_banach(NLS_DATA *, int, REAL *);

/* More or less simplistic vector iteration to compute the spectral
 * condition number of a positive (semi-) definite matrix.
 */
extern REAL vector_iteration(OEM_MV_FCT Av, void *Av_data,
			     REAL *x0, const REAL **kernel, int ker_dim,
			     int dim, REAL tol, int max_iter, int info);
extern REAL matrix_condition(OEM_MV_FCT AIv, OEM_DATA *AIv_data,
			     const REAL **kernel, int ker_dim, int dim,
			     REAL tol, int max_iter, int info);

/*--------------------------------------------------------------------------*/
/*-- Simple Gauss elimination with total pivot search for  -----------------*/
/*-- non-degenerate quadratic matrices.                    -----------------*/
/*--------------------------------------------------------------------------*/
extern void square_gauss(REAL *M, REAL *b, REAL *x, int n, int nbcol);

/*--------------------------------------------------------------------------*/
/*--  inline wrappers to FORTRAN BLAS functions  ---------------------------*/
/*--------------------------------------------------------------------------*/

#include "alberta_util_inlines.h"
#line 940 "../../alberta_util/src/alberta_util.h.in.in"

#ifdef __cplusplus
}
#endif

#endif  /* !_ALBERTA_UTIL_H_ */

/*
 * Local Variables: ***
 * mode: c ***
 * End: ***
 */