/*******************************************************************************
 * ALBERTA:  an Adaptive multi Level finite element toolbox using
 *           Bisectioning refinement and Error control by Residual
 *           Techniques for scientific Applications
 *
 * file:     el_vec.h
 *
 * description: support routines for operation on element vector and matrices
 *
 *******************************************************************************
 *
 *  this file's author(s): Claus-Justus Heine
 *                         Abteilung fuer Angewandte Mathematik
 *                         Albert-Ludwigs-Universitaet Freiburg
 *                         Hermann-Herder-Str. 10
 *                         D-79104 Freiburg im Breisgau, Germany
 *
 *  http://www.alberta-fem.de
 *
 *  (c) by C.-J. Heine (2009)
 *
 ******************************************************************************/

#ifndef _ALBERTA_EL_VEC_H_
#define _ALBERTA_EL_VEC_H_

#include "alberta.h"

/* Initializer for an uninitialized element-vector. */
#define DEFUN_INIT_EL_VEC(VECNAME)			\
  static inline EL_##VECNAME##_VEC *			\
  CPP_CONCAT3(init_el_, VECNAME##_name, _vec)(		\
    EL_##VECNAME##_VEC *vec, int size, int size_max)	\
  {							\
    vec->n_components     = size;			\
    vec->n_components_max = size_max;			\
    vec->reserved         = 1;				\
    DBL_LIST_INIT(&vec->chain);				\
    return vec;						\
  }							\
  struct _AI_semicolon_dummy

DEFUN_INIT_EL_VEC(INT);
DEFUN_INIT_EL_VEC(DOF);
DEFUN_INIT_EL_VEC(REAL);
DEFUN_INIT_EL_VEC(REAL_DD);
/*DEFUN_INIT_EL_VEC(REAL_D);*/
DEFUN_INIT_EL_VEC(UCHAR);
DEFUN_INIT_EL_VEC(SCHAR);
DEFUN_INIT_EL_VEC(PTR);

static inline EL_REAL_D_VEC *
init_el_real_d_vec(EL_REAL_D_VEC *vec, int size, int size_max)
{
  vec->n_components     = size;
  vec->n_components_max = size_max;
  vec->reserved         = DIM_OF_WORLD;
  DBL_LIST_INIT(&vec->chain);

  return vec;
}

#define EL_VEC_SIZEOF(VECNAME, _size_max)	\
  (((_size_max) - 1)				\
   *						\
   sizeof(EL_##VECNAME##_VEC_TYPE)		\
   +						\
   sizeof(EL_##VECNAME##_VEC))

/* An automatic element vector, beware that this is just an
 * uninitialized memory region if _init == false.
 */
#define DEF_EL_VEC_VAR(VECNAME, name, _size, _size_max, _init)	\
  char _AI_##name##_space[EL_VEC_SIZEOF(VECNAME, (_size_max))];	\
  EL_##VECNAME##_VEC *name =					\
    (_init)							\
    ? CPP_CONCAT3(init_el_, VECNAME##_name, _vec)(		\
      (EL_##VECNAME##_VEC *)_AI_##name##_space,			\
      (_size), (_size_max))					\
    : (EL_##VECNAME##_VEC *)_AI_##name##_space

/* Automatic element vector of fixed maximal size. The macro
 * "size_max" needs to be constant.
 */
#define DEF_EL_VEC_CONST(VECNAME, name, _size, _size_max)		\
  struct {								\
      EL_##VECNAME##_VEC el_vec;					\
      EL_##VECNAME##_VEC_TYPE data[MAX(1, (_size_max)-1)];		\
  } _AI_##name = {							\
    {									\
      (_size),     /* n_components */					\
      (_size_max), /* n_components_max */				\
      DBL_LIST_INITIALIZER(_AI_##name.el_vec.chain), /* chain */	\
      ((sizeof(EL_##VECNAME##_VEC_TYPE) + sizeof(REAL) - 1)		\
       / sizeof(REAL)), /* reserved (stride) */				\
    },									\
  };									\
  EL_##VECNAME##_VEC *name = &_AI_##name.el_vec

/* Allocate and initialize a single, un-chained element vector. */
#define ALLOC_EL_VEC(VECNAME, _size, _size_max)		\
  CPP_CONCAT3(init_el_, VECNAME##_name, _vec)(		\
    (EL_##VECNAME##_VEC *)MEM_ALLOC(			\
      EL_VEC_SIZEOF(VECNAME, (_size_max)), char),	\
    (_size), (_size_max))

/* <<< get_dof_indices() */

static inline EL_DOF_VEC *get_dof_indices(EL_DOF_VEC *dofs,
					  const FE_SPACE *fe_space,
					  const EL *el)
{
  const BAS_FCTS *bas_fcts = fe_space->bas_fcts;
  const FE_SPACE *fe_chain;

  if (dofs) {
    CHAIN_DO(fe_space, const FE_SPACE) {
      GET_DOF_INDICES(fe_space->bas_fcts, el, fe_space->admin, dofs->vec);
      dofs->n_components = fe_space->bas_fcts->n_bas_fcts;
      dofs = CHAIN_NEXT(dofs, EL_DOF_VEC);
    } CHAIN_WHILE(fe_space, const FE_SPACE);
  } else {
    EL_DOF_VEC *chain_dofs;

    dofs = (EL_DOF_VEC *)GET_DOF_INDICES(bas_fcts, el, fe_space->admin, NULL);
    dofs->n_components = bas_fcts->n_bas_fcts;
    DBL_LIST_INIT(&dofs->chain);
    CHAIN_FOREACH(fe_chain, fe_space, const FE_SPACE) {
      chain_dofs = (EL_DOF_VEC *)
	GET_DOF_INDICES(fe_chain->bas_fcts, el, fe_chain->admin, NULL);
      chain_dofs->n_components = fe_chain->bas_fcts->n_bas_fcts;
      CHAIN_ADD_TAIL(dofs, chain_dofs);
    }
  }
  return dofs;
}

/* >>> */

/* <<< get_bound() */

static inline EL_BNDRY_VEC *get_bound(EL_BNDRY_VEC *bndry,
				      const BAS_FCTS *bas_fcts,
				      const EL_INFO *el_info)

{
  const BAS_FCTS *bfcts_chain;

  if (bndry) {
    CHAIN_DO(bas_fcts, const BAS_FCTS) {
      GET_BOUND(bas_fcts, el_info, bndry->vec);
      bndry->n_components = bas_fcts->n_bas_fcts;
      bndry = CHAIN_NEXT(bndry, EL_BNDRY_VEC);
    } CHAIN_WHILE(bas_fcts, const BAS_FCTS);
  } else {
    EL_BNDRY_VEC *chain_bndry;

    bndry = (EL_BNDRY_VEC *)GET_BOUND(bas_fcts, el_info, NULL);
    bndry->n_components = bas_fcts->n_bas_fcts;
    DBL_LIST_INIT(&bndry->chain);
    CHAIN_FOREACH(bfcts_chain, bas_fcts, const BAS_FCTS) {
      chain_bndry = (EL_BNDRY_VEC *)GET_BOUND(bfcts_chain, el_info, NULL);
      chain_bndry->n_components = bfcts_chain->n_bas_fcts;
      CHAIN_ADD_TAIL(bndry, chain_bndry);
    }
  }
  return bndry;
}

/* >>> */

/* <<< interpol() */

static inline
void el_interpol(EL_REAL_VEC *coeff,
		 const EL_INFO *el_info, int wall,
		 const EL_INT_VEC *indices,
		 LOC_FCT_AT_QP f, void *ud,
		 const BAS_FCTS *bas_fcts)
{
  const BAS_FCTS *bfcts_chain;
  const int *ind = NULL;
  int n_ind = -1;

  if (indices) {
    ind   = indices->vec;
    n_ind = indices->n_components;
  }
  INTERPOL(bas_fcts, coeff, el_info, wall, n_ind, ind, f, ud);
  CHAIN_FOREACH(bfcts_chain, bas_fcts, const BAS_FCTS) {
    if (indices) {
      indices = CHAIN_NEXT(indices, EL_INT_VEC);
      ind     = indices->vec;
      n_ind   = indices->n_components;
    }
    coeff   = CHAIN_NEXT(coeff, EL_REAL_VEC);
    INTERPOL(bfcts_chain, coeff, el_info, wall, n_ind, ind, f, ud);
  }
}

/* >>> */

/* <<< interpol_dow() */

static inline
void el_interpol_dow(EL_REAL_VEC_D *coeff,
		     const EL_INFO *el_info, int wall,
		     const EL_INT_VEC *indices,
		     LOC_FCT_D_AT_QP f, void *f_data,
		     const BAS_FCTS *bas_fcts)
{
  const BAS_FCTS *bfcts_chain;
  const int *ind = NULL;
  int n_ind = -1;

  if (indices) {
    ind   = indices->vec;
    n_ind = indices->n_components;
  }
  INTERPOL_DOW(bas_fcts, coeff, el_info, wall, n_ind, ind, f, f_data);
  CHAIN_FOREACH(bfcts_chain, bas_fcts, const BAS_FCTS) {
    if (indices) {
      indices = CHAIN_NEXT(indices, EL_INT_VEC);
      ind     = indices->vec;
      n_ind   = indices->n_components;
    }
    coeff   = CHAIN_NEXT(coeff, EL_REAL_VEC_D);
    INTERPOL_DOW(bfcts_chain, coeff, el_info, wall, n_ind, ind, f, f_data);
  }
}

/* >>> */

/* <<< dirichlet_map() */

static inline void __dirichlet_map(S_CHAR *bound,
				   const BNDRY_FLAGS *bndry_bits,
				   int n_components,
				   const BNDRY_FLAGS mask)
{
  int i;
  BNDRY_FLAGS all;

  if (mask == NULL) {
    BNDRY_FLAGS_ALL(all);
    mask = all;
  }

  for (i = 0; i < n_components; i++) {
    if (BNDRY_FLAGS_IS_INTERIOR(bndry_bits[i])) {
      bound[i] = INTERIOR;
      continue;
    }
    if (BNDRY_FLAGS_IS_PARTOF(bndry_bits[i], mask)) {
      bound[i] = DIRICHLET;
    } else {
      bound[i] = INTERIOR;
    }
  }
}

static inline void dirichlet_map_single(EL_SCHAR_VEC *bound,
					const EL_BNDRY_VEC *bndry_bits,
					const BNDRY_FLAGS mask)
{
  bound->n_components = bndry_bits->n_components;

  __dirichlet_map(bound->vec, bndry_bits->vec, bound->n_components, mask);
}

static inline void dirichlet_map(EL_SCHAR_VEC *bound,
				 const EL_BNDRY_VEC *bndry_bits,
				 const BNDRY_FLAGS mask)
{
  const EL_BNDRY_VEC *bndry_chain;

  dirichlet_map_single(bound, bndry_bits, mask);
  CHAIN_FOREACH(bndry_chain, bndry_bits, EL_BNDRY_VEC) {
    bound = CHAIN_NEXT(bound, EL_SCHAR_VEC);
    dirichlet_map_single(bound, bndry_chain, mask);
  }
}

/* >>> */

/* <<< filling el-vecs with data */

#define DEFUN_FILL_EL_VEC(VECNAME, vecname)				\
  static inline const EL_##VECNAME##_VEC *				\
  fill_el_##vecname##_vec(EL_##VECNAME##_VEC *el_vec,			\
			  EL *el,					\
			  const DOF_##VECNAME##_VEC *dof_vec)		\
  {									\
    const FE_SPACE *fe_space = dof_vec->fe_space;			\
    const FE_SPACE *fe_chain;						\
									\
    if (el_vec) {							\
      CHAIN_DO(fe_space, const FE_SPACE) {				\
	fe_space->bas_fcts->get_##vecname##_vec(el_vec->vec, el, dof_vec); \
	el_vec->n_components = fe_space->bas_fcts->n_bas_fcts;		\
	el_vec  = CHAIN_NEXT(el_vec, EL_##VECNAME##_VEC);		\
	dof_vec = CHAIN_NEXT(dof_vec, DOF_##VECNAME##_VEC);		\
      } CHAIN_WHILE(fe_space, const FE_SPACE);				\
    } else {								\
      EL_##VECNAME##_VEC *chain_vec;					\
									\
      el_vec = (EL_##VECNAME##_VEC *)					\
	fe_space->bas_fcts->get_##vecname##_vec(NULL, el, dof_vec);	\
      el_vec->n_components = fe_space->bas_fcts->n_bas_fcts;		\
      DBL_LIST_INIT(&el_vec->chain);					\
      CHAIN_FOREACH(fe_chain, fe_space, const FE_SPACE) {		\
	dof_vec = CHAIN_NEXT(dof_vec, DOF_##VECNAME##_VEC);		\
	chain_vec = (EL_##VECNAME##_VEC *)				\
	  fe_chain->bas_fcts->get_##vecname##_vec(NULL, el, dof_vec);	\
	chain_vec->n_components = fe_chain->bas_fcts->n_bas_fcts;	\
	CHAIN_ADD_TAIL(el_vec, chain_vec);				\
      }									\
    }									\
    return el_vec;							\
  }									\
  struct _AI_semicolon_dummy

DEFUN_FILL_EL_VEC(INT, int);
DEFUN_FILL_EL_VEC(REAL, real);
DEFUN_FILL_EL_VEC(REAL_D, real_d);
DEFUN_FILL_EL_VEC(REAL_DD, real_dd);
DEFUN_FILL_EL_VEC(UCHAR, uchar);
DEFUN_FILL_EL_VEC(SCHAR, schar);

static inline
const EL_REAL_VEC_D *
fill_el_real_vec_d(EL_REAL_VEC_D *el_vec, EL *el, const DOF_REAL_VEC_D *dof_vec)
{
  const FE_SPACE *fe_space = dof_vec->fe_space;
  const FE_SPACE *fe_chain;

  if (el_vec) {
    CHAIN_DO(fe_space, const FE_SPACE) {
      fe_space->bas_fcts->get_real_vec_d(el_vec->vec, el, dof_vec);
      el_vec->n_components = fe_space->bas_fcts->n_bas_fcts;
      el_vec  = CHAIN_NEXT(el_vec, EL_REAL_VEC_D);
      dof_vec = CHAIN_NEXT(dof_vec, DOF_REAL_VEC_D);
    } CHAIN_WHILE(fe_space, const FE_SPACE);
  } else {
    EL_REAL_VEC_D *el_vec_chain;

    el_vec = (EL_REAL_VEC_D *)
      fe_space->bas_fcts->get_real_vec_d(NULL, el, dof_vec);
    el_vec->n_components = fe_space->bas_fcts->n_bas_fcts;
    DBL_LIST_INIT(&el_vec->chain);
    CHAIN_FOREACH(fe_chain, fe_space, const FE_SPACE) {
      dof_vec = CHAIN_NEXT(dof_vec, DOF_REAL_VEC_D);
      el_vec_chain =
	(EL_REAL_VEC_D *)fe_chain->bas_fcts->get_real_vec_d(NULL, el, dof_vec);
      el_vec_chain->n_components = fe_chain->bas_fcts->n_bas_fcts;
      CHAIN_ADD_TAIL(el_vec, el_vec_chain);
    }
  }

  return el_vec;
}

/* >>> */

/* <<< f = c * f + (a A + b B) u, result is an element-vector */

/* We have two kinds of functions: those starting with two underscores
 * are for internal use, they assume fixed data-types and ignore any
 * chain attached to the respective objects.
 *
 * The naming scheme rather treats the beasts as local operators; the
 * different kinds are distinguished by RangeDim/DomainDim
 * combinations.
 *
 */

/* <<< SCAL x SCAL (RANGE x DOMAIN) */

static inline void __el_bi_mat_vec(REAL a, const EL_MATRIX *A,
				   REAL b, const EL_MATRIX *B,
				   const EL_REAL_VEC *u_h,
				   REAL c, EL_REAL_VEC *f_h)
{
  int i, j;
  REAL sum;

  if (A && B) {
    for (i = 0; i < A->n_row; i++) {
      sum = 0.0;
      for (j = 0; j < A->n_col; j++) {
	sum += (a * A->data.real[i][j] + b * B->data.real[i][j]) * u_h->vec[j];
      }
      f_h->vec[i] = c * f_h->vec[i] + sum;
    }
  } else {
    if (B) {
      A = B;
      a = b;
      B = NULL;
      b = 0.0;
    }
    for (i = 0; i < A->n_row; i++) {
      sum = 0.0;
      for (j = 0; j < A->n_col; j++) {
	sum += a * A->data.real[i][j] * u_h->vec[j];
      }
      f_h->vec[i] = c * f_h->vec[i] + sum;
    }
  }
}

__FLATTEN_ATTRIBUTE__
static inline
EL_REAL_VEC *el_bi_mat_vec(REAL a,
			   const EL_MATRIX *A,
			   REAL b,
			   const EL_MATRIX *B,
			   const EL_REAL_VEC *u_h,
			   REAL c,
			   EL_REAL_VEC *f_h)
{
  const EL_MATRIX *A_chain;

  if (A == NULL) {
    A = B;
    a = b;
    B = NULL;
    b = 0.0;
  }
  COL_CHAIN_DO(A, const EL_MATRIX) {
    __el_bi_mat_vec(a, A, b, B, u_h, c, f_h);
    ROW_CHAIN_FOREACH(A_chain, A, const EL_MATRIX) {
      B = B ? ROW_CHAIN_NEXT(B, const EL_MATRIX) : NULL;
      u_h = CHAIN_NEXT(u_h, const EL_REAL_VEC);
      __el_bi_mat_vec(a, A_chain, b, B, u_h, 1.0, f_h);
    }
    /* roll-over to start of this row */
    B = B ? ROW_CHAIN_NEXT(B, const EL_MATRIX) : NULL;
    u_h = CHAIN_NEXT(u_h, const EL_REAL_VEC);
    /* next row */
    f_h = CHAIN_NEXT(f_h, EL_REAL_VEC);
    B = B ? COL_CHAIN_NEXT(B, const EL_MATRIX) : NULL;
  } COL_CHAIN_WHILE(A, const EL_MATRIX);

  return f_h;
}

/* >>> */

/* <<< SCAL x DOW (RANGE x DOMAIN) */

/* <<< REAL x REAL_D, fixed stride */

static inline void __el_bi_mat_vec_rrd(REAL a, const EL_MATRIX *A,
				       REAL b, const EL_MATRIX *B,
				       const EL_REAL_D_VEC *u_h,
				       REAL c, EL_REAL_VEC *f_h)
{
  int i, j;
  REAL sum;

  /* A and B must be either REAL or REAL_D matrices */
  if (A && B) {
    if (A->type == MATENT_REAL) {
      if (B->type == MATENT_REAL) {
	for (j = 0; j < A->n_col; j++) {
	  sum = SUM_DOW(u_h->vec[j]);
	  for (i = 0; i < A->n_row; i++) {
	    f_h->vec[i] =
	      c * f_h->vec[i]
	      +
	      (a * A->data.real[i][j] + b * B->data.real[i][j]) * sum;
	  }
	}
	return;
      } else {
	for (j = 0; j < A->n_col; j++) {
	  REAL u_sum = SUM_DOW(u_h->vec[j]);
	  for (i = 0; i < A->n_row; i++) {
	    f_h->vec[i] =
	      c * f_h->vec[i]
	      +
	      a * A->data.real[i][j] * u_sum
	      +
	      b * SCP_DOW(B->data.real_d[i][j], u_h->vec[j]);
	  }
	}
	return;
      }
    } else {
      if (B->type == MATENT_REAL) {
	__el_bi_mat_vec_rrd(b, B, a, A, u_h, c, f_h);
      } else {
	for (j = 0; j < A->n_col; j++) {
	  for (i = 0; i < A->n_row; i++) {
	    f_h->vec[i] =
	      c * f_h->vec[i]
	      +
	      a * SCP_DOW(A->data.real_d[i][j], u_h->vec[j])
	      +
	      b * SCP_DOW(B->data.real_d[i][j], u_h->vec[j]);
	  }
	}
      }
      return;
    }
  } else if (A) {
    if (A->type == MATENT_REAL) {
      for (j = 0; j < A->n_col; j++) {
	sum = SUM_DOW(u_h->vec[j]);
	for (i = 0; i < A->n_row; i++) {
	  f_h->vec[i] = c * f_h->vec[i] + a * A->data.real[i][j] * sum;
	}
      }
    } else {
      for (i = 0; i < A->n_row; i++) {
	for (j = 0; j < A->n_col; j++) {
	  f_h->vec[i] =
	    c * f_h->vec[i] + a * SCP_DOW(A->data.real_d[i][j], u_h->vec[j]);
	}
      }
    }
    return;
  }
}

__FLATTEN_ATTRIBUTE__
static inline
EL_REAL_VEC *el_bi_mat_vec_rrd(REAL a,
			       const EL_MATRIX *A,
			       REAL b,
			       const EL_MATRIX *B,
			       const EL_REAL_D_VEC *u_h,
			       REAL c,
			       EL_REAL_VEC *f_h)
{
  const EL_MATRIX *A_chain;

  if (A == NULL) {
    A = B;
    a = b;
    B = NULL;
    b = 0.0;
  }
  COL_CHAIN_DO(A, const EL_MATRIX) {
    __el_bi_mat_vec_rrd(a, A, b, B, u_h, c, f_h);
    ROW_CHAIN_FOREACH(A_chain, A, const EL_MATRIX) {
      B = B ? ROW_CHAIN_NEXT(B, const EL_MATRIX) : NULL;
      u_h = CHAIN_NEXT(u_h, const EL_REAL_D_VEC);
      __el_bi_mat_vec_rrd(a, A_chain, b, B, u_h, 1.0, f_h);
    }
    /* roll over to start of this row */
    B = B ? ROW_CHAIN_NEXT(B, const EL_MATRIX) : NULL;
    u_h = CHAIN_NEXT(u_h, const EL_REAL_D_VEC);
    /* next row */
    f_h = CHAIN_NEXT(f_h, EL_REAL_VEC);
    B = B ? COL_CHAIN_NEXT(B, const EL_MATRIX) : NULL;
  } COL_CHAIN_WHILE(A, const EL_MATRIX);

  return f_h;
}

/* >>> */

/* <<< REAL x REAL_D, vaiable stride */

static inline
void __el_bi_mat_vec_scl_dow(REAL a,
			     const EL_MATRIX *A,
			     REAL b,
			     const EL_MATRIX *B,
			     const EL_REAL_VEC_D *u_h,
			     REAL c,
			     EL_REAL_VEC *f_h)
{
  if (u_h->stride != 1) {
    __el_bi_mat_vec_rrd(a, A, b, B, (const EL_REAL_D_VEC *)u_h, c, f_h);
  } else {
    __el_bi_mat_vec(a, A, b, B, (const EL_REAL_VEC *)u_h, c, f_h);
  }
}

__FLATTEN_ATTRIBUTE__
static inline
EL_REAL_VEC *el_bi_mat_vec_scl_dow(REAL a,
				   const EL_MATRIX *A,
				   REAL b,
				   const EL_MATRIX *B,
				   const EL_REAL_VEC_D *u_h,
				   REAL c,
				   EL_REAL_VEC *f_h)
{
  const EL_MATRIX *A_chain;
  if (A == NULL) {
    A = B;
    a = b;
    B = NULL;
    b = 0.0;
  }
  COL_CHAIN_DO(A, const EL_MATRIX) {
    __el_bi_mat_vec_scl_dow(a, A, b, B, u_h, c, f_h);
    ROW_CHAIN_FOREACH(A_chain, A, const EL_MATRIX) {
      B = B ? ROW_CHAIN_NEXT(B, const EL_MATRIX) : NULL;
      u_h = CHAIN_NEXT(u_h, const EL_REAL_VEC_D);
      __el_bi_mat_vec_scl_dow(a, A_chain, b, B, u_h, 1.0, f_h);
    }
    /* roll-over to start of row */
    B = B ? ROW_CHAIN_NEXT(B, const EL_MATRIX) : NULL;
    u_h = CHAIN_NEXT(u_h, const EL_REAL_VEC_D);
    /* next row */
    B = B ? COL_CHAIN_NEXT(B, const EL_MATRIX) : NULL;
    f_h = CHAIN_NEXT(f_h, EL_REAL_VEC);
  } COL_CHAIN_WHILE(A, const EL_MATRIX);

  return f_h;
}

/* >>> */

/* >>> */

/* <<< DOW x SCAL (RANGE x DOMAIN) */

/* <<< READ_D x REAL, fixed stride */

static inline void
__el_bi_mat_vec_rdr(REAL a, const EL_MATRIX *A,
		    REAL b, const EL_MATRIX *B,
		    const EL_REAL_VEC *u_h,
		    REAL c, EL_REAL_D_VEC *f_h)
{
  int i, j;
  REAL val;

  if (A && B) {
    if (A->type == MATENT_REAL) {
      if (B->type == MATENT_REAL) {
	for (i = 0; i < A->n_row; i++) {
	  val = 0.0;
	  for (j = 0; j < A->n_col; j++) {
	    val +=
	      (a * A->data.real[i][j] + b * B->data.real[i][j]) * u_h->vec[j];
	  }
	  for (j = 0; j < DIM_OF_WORLD; j++) {
	    f_h->vec[i][j] = c * f_h->vec[i][j] + val;
	  }
	}
	return;
      } else {
	for (i = 0; i < A->n_row; i++) {
	  val = 0.0;
	  for (j = 0; j < A->n_col; j++) {
	    val += a * A->data.real[i][j] * u_h->vec[j];
	  }
	  for (j = 0; j < DIM_OF_WORLD; j++) {
	    f_h->vec[i][j] = c * f_h->vec[i][j] + val;
	  }
	  for (j = 0; j < A->n_col; j++) {
	    AXPY_DOW(b * u_h->vec[j], B->data.real_d[i][j], f_h->vec[i]);
	  }
	}
	return;
      }
    } else {
      if (B->type == MATENT_REAL) {
	__el_bi_mat_vec_rdr(b, B, a, A, u_h, c, f_h);
      } else {
	for (i = 0; i < A->n_row; i++) {
	  for (j = 0; j < DIM_OF_WORLD; j++) {
	    f_h->vec[i][j] = c * f_h->vec[i][j];
	  }
	  for (j = 0; j < A->n_col; j++) {
	    AXPY_DOW(a * u_h->vec[j], A->data.real_d[i][j], f_h->vec[i]);
	    AXPY_DOW(b * u_h->vec[j], B->data.real_d[i][j], f_h->vec[i]);
	  }
	}
      }
      return;
    }
  } else if (A) {
    if (A->type == MATENT_REAL) {
      for (i = 0; i < A->n_row; i++) {
	val = 0.0;
	for (j = 0; j < A->n_col; j++) {
	  val += a * A->data.real[i][j] * u_h->vec[j];
	}
	for (j = 0; j < DIM_OF_WORLD; j++) {
	  f_h->vec[i][j] = c * f_h->vec[i][j] + val;
	}
      }
    } else {
      for (i = 0; i < A->n_row; i++) {
	for (j = 0; j < A->n_col; j++) {
	  AXPBY_DOW(a * u_h->vec[j], A->data.real_d[i][j], c, f_h->vec[i],
		    f_h->vec[i]);
	}
      }
    }
    return;
  }
}

__FLATTEN_ATTRIBUTE__
static inline
EL_REAL_D_VEC *el_bi_mat_vec_rdr(REAL a,
				 const EL_MATRIX *A,
				 REAL b,
				 const EL_MATRIX *B,
				 const EL_REAL_VEC *u_h,
				 REAL c,
				 EL_REAL_D_VEC *f_h)
{
  const EL_MATRIX *A_chain;

  if (A == NULL) {
    A = B;
    a = b;
    B = NULL;
    b = 0.0;
  }
  COL_CHAIN_DO(A, const EL_MATRIX) {
    __el_bi_mat_vec_rdr(a, A, b, B, u_h, c, f_h);
    ROW_CHAIN_FOREACH(A_chain, A, const EL_MATRIX) {
      B = B ? ROW_CHAIN_NEXT(B, const EL_MATRIX) : NULL;
      u_h = CHAIN_NEXT(u_h, const EL_REAL_VEC);
      __el_bi_mat_vec_rdr(a, A, b, B, u_h, 1.0, f_h);
    }
    /* roll-over to start of row */
    B = B ? ROW_CHAIN_NEXT(B, const EL_MATRIX) : NULL;
    u_h = CHAIN_NEXT(u_h, const EL_REAL_VEC);
    /* next row */
    B = B ? COL_CHAIN_NEXT(B, const EL_MATRIX) : NULL;
    f_h = CHAIN_NEXT(f_h, EL_REAL_D_VEC);
  } COL_CHAIN_WHILE(A, const EL_MATRIX);

  return f_h;
}

/* >>> */

/* <<< REAL_D x REAL, variable stride */

static inline
void __el_bi_mat_vec_dow_scl(REAL a,
			     const EL_MATRIX *A,
			     REAL b,
			     const EL_MATRIX *B,
			     const EL_REAL_VEC *u_h,
			     REAL c,
			     EL_REAL_VEC_D *f_h)
{
  if (f_h->stride != 1) {
    __el_bi_mat_vec_rdr(a, A, b, B, u_h, c, (EL_REAL_D_VEC *)f_h);
  } else {
    __el_bi_mat_vec(a, A, b, B, u_h, c, (EL_REAL_VEC *)f_h);
  }
}

__FLATTEN_ATTRIBUTE__
static inline
EL_REAL_VEC_D *el_bi_mat_vec_dow_scl(REAL a,
				     const EL_MATRIX *A,
				     REAL b,
				     const EL_MATRIX *B,
				     const EL_REAL_VEC *u_h,
				     REAL c,
				     EL_REAL_VEC_D *f_h)
{
  const EL_MATRIX *A_chain;

  if (A == NULL) {
    A = B;
    a = b;
    B = NULL;
    b = 0.0;
  }
  COL_CHAIN_DO(A, const EL_MATRIX) {
    __el_bi_mat_vec_dow_scl(a, A, b, B, u_h, c, f_h);
    ROW_CHAIN_FOREACH(A_chain, A, const EL_MATRIX) {
      B = B ? ROW_CHAIN_NEXT(B, const EL_MATRIX) : NULL;
      u_h = CHAIN_NEXT(u_h, const EL_REAL_VEC);
      __el_bi_mat_vec_dow_scl(a, A, b, B, u_h, 1.0, f_h);
    }
    /* roll-over to start of row */
    B = B ? ROW_CHAIN_NEXT(B, const EL_MATRIX) : NULL;
    u_h = CHAIN_NEXT(u_h, const EL_REAL_VEC);
    /* next row */
    B = B ? COL_CHAIN_NEXT(B, const EL_MATRIX) : NULL;
    f_h = CHAIN_NEXT(f_h, EL_REAL_VEC_D);
  } COL_CHAIN_WHILE(A, const EL_MATRIX);

  return f_h;
}

/* >>> */

/* >>> */

/* <<< DOW x DOW (RANGE x DOMAIN) */

/* <<< fixed stride */

/* Multiply an element matrix by an element vector and store the
 * result inside an element vector.
 */
static inline
void __el_bi_mat_vec_d(REAL a, const EL_MATRIX *A,
		       REAL b, const EL_MATRIX *B,
		       const EL_REAL_D_VEC *u_h,
		       REAL c, EL_REAL_D_VEC *f_h)
{
  int i, j;
  REAL *sum;

  if (A && B) {
#undef MAT_BI_BODY
#define MAT_BI_BODY(PFX1, CC1, C1, S1, TYPE1,		\
		    PFX2, CC2, C2, S2, TYPE2)		\
    for (i = 0; i < A->n_row; i++) {			\
      sum = f_h->vec[i];				\
      for (j = 0; j < A->n_col; j++) {			\
	PFX1##PFX2##BIMV_DOW(a, CC1 A->data.S1[i][j],	\
			     b, CC2 B->data.S2[i][j],	\
			     u_h->vec[j],		\
			     c, sum);			\
      }							\
    }
    MAT_EMIT_BI_BODY_SWITCH(A->type, B->type);
  } else {
    if (A == NULL) {
      A = B;
      a = b;
      B = NULL;
      b = 0.0;
    }
#undef MAT_BODY
#define MAT_BODY(PFX, CC, C, S, TYPE)					\
    for (i = 0; i < A->n_row; i++) {					\
      sum = f_h->vec[i];						\
      for (j = 0; j < A->n_col; j++) {					\
	PFX##GEMV_DOW(a, CC A->data.S[i][j], u_h->vec[j], c, sum);	\
      }									\
    }
    MAT_EMIT_BODY_SWITCH(A->type);
  }
}

__FLATTEN_ATTRIBUTE__
static inline
EL_REAL_D_VEC *el_bi_mat_vec_d(REAL a,
			       const EL_MATRIX *A,
			       REAL b,
			       const EL_MATRIX *B,
			       const EL_REAL_D_VEC *u_h,
			       REAL c,
			       EL_REAL_D_VEC *f_h)
{
  const EL_MATRIX *A_chain;

  if (A == NULL) {
    A = B;
    a = b;
    B = NULL;
    b = 0.0;
  }
  COL_CHAIN_DO(A, const EL_MATRIX) {
    __el_bi_mat_vec_d(a, A, b, B, u_h, c, f_h);
    ROW_CHAIN_FOREACH(A_chain, A, const EL_MATRIX) {
      /* next column */
      B = B ? ROW_CHAIN_NEXT(B, const EL_MATRIX) : NULL;
      u_h = CHAIN_NEXT(u_h, const EL_REAL_D_VEC);
      __el_bi_mat_vec_d(a, A, b, B, u_h, 1.0, f_h);
    }
    /* roll-over to start of row */
    B = B ? ROW_CHAIN_NEXT(B, const EL_MATRIX) : NULL;
    u_h = CHAIN_NEXT(u_h, const EL_REAL_D_VEC);
    /* next row */
    B = B ? COL_CHAIN_NEXT(B, const EL_MATRIX) : NULL;
    f_h = CHAIN_NEXT(f_h, EL_REAL_D_VEC);
  } COL_CHAIN_WHILE(A, const EL_MATRIX);

  return f_h;
}

/* >>> */

/* <<< REAL_D x REAL_D, variable stride */

static inline
void __el_bi_mat_vec_dow(REAL a,
			 const EL_MATRIX *A,
			 REAL b,
			 const EL_MATRIX *B,
			 const EL_REAL_VEC_D *u_h,
			 REAL c,
			 EL_REAL_VEC_D *f_h)
{
  if (f_h->stride == 1) {
    __el_bi_mat_vec_scl_dow(a, A, b, B, u_h, c, (EL_REAL_VEC *)f_h);
  } else if (u_h->stride == 1) {
    __el_bi_mat_vec_dow_scl(a, A, b, B, (const EL_REAL_VEC *)u_h, c, f_h);
  } else {
    __el_bi_mat_vec_d(a, A, b, B, (const EL_REAL_D_VEC *)u_h,
		      c, (EL_REAL_D_VEC *)f_h);
  }
}

__FLATTEN_ATTRIBUTE__
static inline
EL_REAL_VEC_D *el_bi_mat_vec_dow(REAL a,
				 const EL_MATRIX *A,
				 REAL b,
				 const EL_MATRIX *B,
				 const EL_REAL_VEC_D *u_h,
				 REAL c,
				 EL_REAL_VEC_D *f_h)
{
  const EL_MATRIX *A_chain;

  if (A == NULL) {
    A = B;
    a = b;
    B = NULL;
    b = 0.0;
  }
  COL_CHAIN_DO(A, const EL_MATRIX) {
    __el_bi_mat_vec_dow(a, A, b, B, u_h, c, f_h);
    ROW_CHAIN_FOREACH(A_chain, A, const EL_MATRIX) {
      /* next column */
      B = B ? ROW_CHAIN_NEXT(B, const EL_MATRIX) : NULL;
      u_h = CHAIN_NEXT(u_h, const EL_REAL_VEC_D);
      __el_bi_mat_vec_dow(a, A, b, B, u_h, 1.0, f_h);
    }
    /* roll-over to start of row */
    B = B ? ROW_CHAIN_NEXT(B, const EL_MATRIX) : NULL;
    u_h = CHAIN_NEXT(u_h, const EL_REAL_VEC_D);
    /* next row */
    f_h = CHAIN_NEXT(f_h, EL_REAL_VEC_D);
    B = B ? COL_CHAIN_NEXT(B, const EL_MATRIX) : NULL;
  } COL_CHAIN_WHILE(A, const EL_MATRIX);

  return f_h;
}

/* >>> */

/* >>> */

/* <<< f = c * f + a A u, f += a A u versions */

/* <<< SCL x SCL */

static inline EL_REAL_VEC *
el_gen_mat_vec(REAL a, const EL_MATRIX *A, const EL_REAL_VEC *u_h,
	       REAL b, EL_REAL_VEC *f_h)
{
  return el_bi_mat_vec(a, A, 0.0, NULL, u_h, b, f_h);
}

static inline EL_REAL_VEC *
el_mat_vec(REAL a, const EL_MATRIX *A, const EL_REAL_VEC *u_h, EL_REAL_VEC *f_h)
{
  return el_gen_mat_vec(a, A, u_h, 1.0, f_h);
}

/* >>> */

/* <<< SCL x DOW */

static inline EL_REAL_VEC *
el_gen_mat_vec_rrd(REAL a, const EL_MATRIX *A, const EL_REAL_D_VEC *u_h,
		   REAL b, EL_REAL_VEC *f_h)
{
  return el_bi_mat_vec_rrd(a, A, 0.0, NULL, u_h, b, f_h);
}

static inline EL_REAL_VEC *
el_mat_vec_rrd(REAL a, const EL_MATRIX *A, const EL_REAL_D_VEC *u_h,
	       EL_REAL_VEC *f_h)
{
  return el_gen_mat_vec_rrd(a, A, u_h, 1.0, f_h);
}

static inline EL_REAL_VEC *
el_gen_mat_vec_scl_dow(REAL a, const EL_MATRIX *A, const EL_REAL_VEC_D *u_h,
		       REAL b, EL_REAL_VEC *f_h)
{
  return el_bi_mat_vec_scl_dow(a, A, 0.0, NULL, u_h, b, f_h);
}

static inline EL_REAL_VEC *
el_mat_vec_scl_dow(REAL a, const EL_MATRIX *A, const EL_REAL_VEC_D *u_h,
		   EL_REAL_VEC *f_h)
{
  return el_gen_mat_vec_scl_dow(a, A, u_h, 1.0, f_h);
}

/* >>> */

/* <<< DOW x SCL */

static inline EL_REAL_D_VEC *
el_gen_mat_vec_rdr(REAL a, const EL_MATRIX *A, const EL_REAL_VEC *u_h,
		   REAL b, EL_REAL_D_VEC *f_h)
{
  return el_bi_mat_vec_rdr(a, A, 0.0, NULL, u_h, b, f_h);
}

static inline EL_REAL_D_VEC *
el_mat_vec_rdr(REAL a, const EL_MATRIX *A, const EL_REAL_VEC *u_h,
	       EL_REAL_D_VEC *f_h)
{
  return el_gen_mat_vec_rdr(a, A, u_h, 1.0, f_h);
}

static inline EL_REAL_VEC_D *
el_gen_mat_vec_dow_scl(REAL a, const EL_MATRIX *A, const EL_REAL_VEC *u_h,
		       REAL b, EL_REAL_VEC_D *f_h)
{
  return el_bi_mat_vec_dow_scl(a, A, 0.0, NULL, u_h, b, f_h);
}

static inline EL_REAL_VEC_D *
el_mat_vec_dow_scl(REAL a, const EL_MATRIX *A, const EL_REAL_VEC *u_h,
		   EL_REAL_VEC_D *f_h)
{
  return el_gen_mat_vec_dow_scl(a, A, u_h, 1.0, f_h);
}

/* >>> */

/* <<< DOW x DOW */

static inline EL_REAL_D_VEC *
el_gen_mat_vec_d(REAL a, const EL_MATRIX *A, const EL_REAL_D_VEC *u_h,
		 REAL b, EL_REAL_D_VEC *f_h)
{
  return el_bi_mat_vec_d(a, A, 0.0, NULL, u_h, b, f_h);
}

static inline EL_REAL_D_VEC *
el_mat_vec_d(REAL a, const EL_MATRIX *A, const EL_REAL_D_VEC *u_h,
	     EL_REAL_D_VEC *f_h)
{
  return el_gen_mat_vec_d(a, A, u_h, 1.0, f_h);
}

static inline EL_REAL_VEC_D *
el_gen_mat_vec_dow(REAL a, const EL_MATRIX *A, const EL_REAL_VEC_D *u_h,
		   REAL b, EL_REAL_VEC_D *f_h)
{
  return el_bi_mat_vec_dow(a, A, 0.0, NULL, u_h, b, f_h);
}

static inline EL_REAL_VEC_D *
el_mat_vec_dow(REAL a, const EL_MATRIX *A, const EL_REAL_VEC_D *u_h,
	       EL_REAL_VEC_D *f_h)
{
  return el_gen_mat_vec_dow(a, A, u_h, 1.0, f_h);
}

/* >>> */

/* >>> */

/* >>> */

/* <<< f = c * f + (a A + b B) u, result is a global DOF-vector */

/* <<< SCAL-SCAL */

static inline
void __bi_mat_el_vec(REAL a, const EL_MATRIX *A,
		     REAL b, const EL_MATRIX *B,
		     const EL_REAL_VEC *u_h,
		     REAL c, DOF_REAL_VEC *f_h,
		     const EL_DOF_VEC *dof,
		     const EL_SCHAR_VEC *bnd)
{
  int i, j;
  REAL sum;

  if (A && B) {
    for (i = 0; i < A->n_row; i++) {
      if (bnd == NULL || bnd->vec[i] < DIRICHLET) {
	sum = 0.0;
	for (j = 0; j < A->n_col; j++) {
	  sum +=
	    (a * A->data.real[i][j] + b * B->data.real[i][j]) * u_h->vec[j];
	}
	f_h->vec[dof->vec[i]] = c * f_h->vec[dof->vec[i]] + sum;
      }
    }
  } else {
    if (B) {
      A = B;
      a = b;
      B = NULL;
      b = 0.0;
    }
    for (i = 0; i < A->n_row; i++) {
      if (bnd == NULL || bnd->vec[i] < DIRICHLET) {
	sum = 0.0;
	for (j = 0; j < A->n_col; j++) {
	  sum += a * A->data.real[i][j] * u_h->vec[j];
	}
	f_h->vec[dof->vec[i]] = c * f_h->vec[dof->vec[i]] + sum;
      }
    }
  }
}

__FLATTEN_ATTRIBUTE__
static inline
void bi_mat_el_vec(REAL a,
		   const EL_MATRIX *A,
		   REAL b,
		   const EL_MATRIX *B,
		   const EL_REAL_VEC *u_h,
		   REAL c,
		   DOF_REAL_VEC *f_h,
		   const EL_DOF_VEC *dof,
		   const EL_SCHAR_VEC *bnd)
{
  if (A == NULL) {
    A = B;
    a = b;
    B = NULL;
    b = 0.0;
  }
  COL_CHAIN_DO(A, const EL_MATRIX) {
    ROW_CHAIN_DO(A, const EL_MATRIX) {

      __bi_mat_el_vec(a, A, b, B, u_h, c, f_h, dof, bnd);

      B = B ? ROW_CHAIN_NEXT(B, const EL_MATRIX) : NULL;
      u_h = CHAIN_NEXT(u_h, const EL_REAL_VEC);
    } ROW_CHAIN_WHILE(A, const EL_MATRIX);
    f_h   = CHAIN_NEXT(f_h, DOF_REAL_VEC);
    dof   = CHAIN_NEXT(dof, EL_DOF_VEC);
    bnd   = bnd ? CHAIN_NEXT(bnd, EL_SCHAR_VEC) : NULL;
    B     = B ? COL_CHAIN_NEXT(B, const EL_MATRIX) : NULL;
  } COL_CHAIN_WHILE(A, const EL_MATRIX);
}

/* >>> */

/* <<< SCAL-DOW version */

/* <<< fixed stride */

static inline
void __bi_mat_el_vec_rrd(REAL a, const EL_MATRIX *A,
			 REAL b, const EL_MATRIX *B,
			 const EL_REAL_D_VEC *u_h,
			 REAL c, DOF_REAL_VEC *f_h,
			 const EL_DOF_VEC *dof,
			 const EL_SCHAR_VEC *bnd)
{
  int i, j;
  REAL sum;

  /* A and B must be either REAL or REAL_D matrices */
  if (A && B) {
    if (A->type == MATENT_REAL) {
      if (B->type == MATENT_REAL) {
	for (j = 0; j < A->n_col; j++) {
	  sum = SUM_DOW(u_h->vec[j]);
	  for (i = 0; i < A->n_row; i++) {
	    if (bnd == NULL || bnd->vec[i] < DIRICHLET) {
	      f_h->vec[dof->vec[i]] =
		c * f_h->vec[dof->vec[i]]
		+
		(a * A->data.real[i][j] + b * B->data.real[i][j]) * sum;
	    }
	  }
	}
	return;
      } else {
	for (j = 0; j < A->n_col; j++) {
	  REAL u_sum = SUM_DOW(u_h->vec[j]);
	  for (i = 0; i < A->n_row; i++) {
	    if (bnd == NULL || bnd->vec[i] < DIRICHLET) {
	      f_h->vec[dof->vec[i]] =
		c * f_h->vec[dof->vec[i]]
		+
		a * A->data.real[i][j] * u_sum
		+
		b * SCP_DOW(B->data.real_d[i][j], u_h->vec[j]);
	    }
	  }
	}
	return;
      }
    } else {
      if (B->type == MATENT_REAL) {
	__bi_mat_el_vec_rrd(b, B, a, A, u_h, c, f_h, dof, bnd);
      } else {
	for (j = 0; j < A->n_col; j++) {
	  for (i = 0; i < A->n_row; i++) {
	    if (bnd == NULL || bnd->vec[i] < DIRICHLET) {
	      f_h->vec[dof->vec[i]] =
		c * f_h->vec[dof->vec[i]]
		+
		a * SCP_DOW(A->data.real_d[i][j], u_h->vec[j])
		+
		b * SCP_DOW(B->data.real_d[i][j], u_h->vec[j]);
	    }
	  }
	}
      }
      return;
    }
  } else if (A) {
    if (A->type == MATENT_REAL) {
      for (j = 0; j < A->n_col; j++) {
	sum = SUM_DOW(u_h->vec[j]);
	for (i = 0; i < A->n_row; i++) {
	  if (bnd == NULL || bnd->vec[i] < DIRICHLET) {
	    f_h->vec[dof->vec[i]] =
	      c * f_h->vec[dof->vec[i]] + + a * A->data.real[i][j] * sum;
	  }
	}
      }
    } else {
      for (i = 0; i < A->n_row; i++) {
	for (j = 0; j < A->n_col; j++) {
	  if (bnd == NULL || bnd->vec[i] < DIRICHLET) {
	    f_h->vec[dof->vec[i]] =
	      c * f_h->vec[dof->vec[i]]
	      +
	      a * SCP_DOW(A->data.real_d[i][j], u_h->vec[j]);
	  }
	}
      }
    }
    return;
  }
}

__FLATTEN_ATTRIBUTE__
static inline
void bi_mat_el_vec_rrd(REAL a,
		       const EL_MATRIX *A,
		       REAL b,
		       const EL_MATRIX *B,
		       const EL_REAL_D_VEC *u_h,
		       REAL c,
		       DOF_REAL_VEC *f_h,
		       const EL_DOF_VEC *dof,
		       const EL_SCHAR_VEC *bnd)
{
  if (A == NULL) {
    A = B;
    a = b;
    B = NULL;
    b = 0.0;
  }
  COL_CHAIN_DO(A, const EL_MATRIX) {
    ROW_CHAIN_DO(A, const EL_MATRIX) {

      __bi_mat_el_vec_rrd(a, A, b, B, u_h, c, f_h, dof, bnd);

      B = B ? ROW_CHAIN_NEXT(B, const EL_MATRIX) : NULL;
      u_h = CHAIN_NEXT(u_h, const EL_REAL_D_VEC);
    } ROW_CHAIN_WHILE(A, const EL_MATRIX);
    f_h   = CHAIN_NEXT(f_h, DOF_REAL_VEC);
    dof   = CHAIN_NEXT(dof, EL_DOF_VEC);
    bnd   = bnd ? CHAIN_NEXT(bnd, EL_SCHAR_VEC) : NULL;
    B     = B ? COL_CHAIN_NEXT(B, const EL_MATRIX) : NULL;
  } COL_CHAIN_WHILE(A, const EL_MATRIX);
}

/* >>> */

/* <<< variable stride version */

static inline
void __bi_mat_el_vec_scl_dow(REAL a, const EL_MATRIX *A,
			     REAL b, const EL_MATRIX *B,
			     const EL_REAL_VEC_D *u_h,
			     REAL c, DOF_REAL_VEC *f_h,
			     const EL_DOF_VEC *dof,
			     const EL_SCHAR_VEC *bnd)
{
  if (u_h->stride != 1) {
    __bi_mat_el_vec_rrd(a, A, b, B, (const EL_REAL_D_VEC *)u_h,
			c, f_h, dof, bnd);
  } else {
    __bi_mat_el_vec(a, A, b, B, (const EL_REAL_VEC *)u_h,
		    c, f_h, dof, bnd);
  }
}

__FLATTEN_ATTRIBUTE__
static inline
void bi_mat_el_vec_scl_dow(REAL a,
			   const EL_MATRIX *A,
			   REAL b,
			   const EL_MATRIX *B,
			   const EL_REAL_VEC_D *u_h,
			   REAL c,
			   DOF_REAL_VEC *f_h,
			   const EL_DOF_VEC *dof,
			   const EL_SCHAR_VEC *bnd)
{
  if (A == NULL) {
    A = B;
    a = b;
    B = NULL;
    b = 0.0;
  }
  COL_CHAIN_DO(A, const EL_MATRIX) {
    ROW_CHAIN_DO(A, const EL_MATRIX) {

      __bi_mat_el_vec_scl_dow(a, A, b, B, u_h, c, f_h, dof, bnd);

      B = B ? ROW_CHAIN_NEXT(B, const EL_MATRIX) : NULL;
      u_h = CHAIN_NEXT(u_h, const EL_REAL_VEC_D);
    } ROW_CHAIN_WHILE(A, const EL_MATRIX);
    f_h   = CHAIN_NEXT(f_h, DOF_REAL_VEC);
    dof   = CHAIN_NEXT(dof, EL_DOF_VEC);
    bnd   = bnd ? CHAIN_NEXT(bnd, EL_SCHAR_VEC) : NULL;
    B     = B ? COL_CHAIN_NEXT(B, const EL_MATRIX) : NULL;
  } COL_CHAIN_WHILE(A, const EL_MATRIX);
}

/* >>> */

/* >>> */

/* <<< DOW-SCAL version */

/* <<< fixed stride */

static inline
void __bi_mat_el_vec_rdr(REAL a, const EL_MATRIX *A,
			 REAL b, const EL_MATRIX *B,
			 const EL_REAL_VEC *u_h,
			 REAL c, DOF_REAL_D_VEC *f_h,
			 const EL_DOF_VEC *dof,
			 const EL_SCHAR_VEC *bnd)
{
  int i, j;
  REAL val;

  if (A && B) {
    if (A->type == MATENT_REAL) {
      if (B->type == MATENT_REAL) {
	for (i = 0; i < A->n_row; i++) {
	  if (bnd == NULL || bnd->vec[i] < DIRICHLET) {
	    val = 0.0;
	    for (j = 0; j < A->n_col; j++) {
	      val +=
		(a * A->data.real[i][j] + b * B->data.real[i][j]) * u_h->vec[j];
	    }
	    for (j = 0; j < DIM_OF_WORLD; j++) {
	      f_h->vec[dof->vec[i]][j] = c * f_h->vec[dof->vec[i]][j] + val;
	    }
	  }
	}
	return;
      } else {
	for (i = 0; i < A->n_row; i++) {
	  if (bnd == NULL || bnd->vec[i] < DIRICHLET) {
	    val = 0.0;
	    for (j = 0; j < A->n_col; j++) {
	      val += a * A->data.real[i][j] * u_h->vec[j];
	    }
	    for (j = 0; j < DIM_OF_WORLD; j++) {
	      f_h->vec[dof->vec[i]][j] = c * f_h->vec[dof->vec[i]][j] + val;
	    }
	    for (j = 0; j < A->n_col; j++) {
	      AXPY_DOW(b * u_h->vec[j], B->data.real_d[i][j],
		       f_h->vec[dof->vec[i]]);
	    }
	  }
	}
	return;
      }
    } else {
      if (B->type == MATENT_REAL) {
	__bi_mat_el_vec_rdr(b, B, a, A, u_h, c, f_h, dof, bnd);
      } else {
	for (i = 0; i < A->n_row; i++) {
	  if (bnd == NULL || bnd->vec[i] < DIRICHLET) {
	    for (j = 0; j < DIM_OF_WORLD; j++) {
	      f_h->vec[dof->vec[i]][j] = c * f_h->vec[dof->vec[i]][j];
	    }
	    for (j = 0; j < A->n_col; j++) {
	      AXPY_DOW(a * u_h->vec[j], A->data.real_d[i][j],
		       f_h->vec[dof->vec[i]]);
	      AXPY_DOW(b * u_h->vec[j], B->data.real_d[i][j],
		       f_h->vec[dof->vec[i]]);
	    }
	  }
	}
      }
    }
  } else if (A) {
    if (A->type == MATENT_REAL) {
      for (i = 0; i < A->n_row; i++) {
	if (bnd == NULL || bnd->vec[i] < DIRICHLET) {
	  val = 0.0;
	  for (j = 0; j < A->n_col; j++) {
	    val += a * A->data.real[i][j] * u_h->vec[j];
	  }
	  for (j = 0; j < DIM_OF_WORLD; j++) {
	    f_h->vec[dof->vec[i]][j] = c * f_h->vec[dof->vec[i]][j] + val;
	  }
	}
      }
    } else {
      for (i = 0; i < A->n_row; i++) {
	if (bnd == NULL || bnd->vec[i] < DIRICHLET) {
	  for (j = 0; j < A->n_col; j++) {
	    AXPBY_DOW(a * u_h->vec[j], A->data.real_d[i][j],
		      c, f_h->vec[dof->vec[i]],
		      f_h->vec[dof->vec[i]]);
	  }
	}
      }
    }
    return;
  }
}

static inline void bi_mat_el_vec_rdr(REAL a,
				     const EL_MATRIX *A,
				     REAL b,
				     const EL_MATRIX *B,
				     const EL_REAL_VEC *u_h,
				     REAL c,
				     DOF_REAL_D_VEC *f_h,
				     const EL_DOF_VEC *dof,
				     const EL_SCHAR_VEC *bnd)
{
  if (A == NULL) {
    A = B;
    a = b;
    B = NULL;
    b = 0.0;
  }
  COL_CHAIN_DO(A, const EL_MATRIX) {
    ROW_CHAIN_DO(A, const EL_MATRIX) {

      __bi_mat_el_vec_rdr(a, A, b, B, u_h, c, f_h, dof, bnd);

      B = B ? ROW_CHAIN_NEXT(B, const EL_MATRIX) : NULL;
      u_h = CHAIN_NEXT(u_h, const EL_REAL_VEC);
    } ROW_CHAIN_WHILE(A, const EL_MATRIX);
    f_h   = CHAIN_NEXT(f_h, DOF_REAL_D_VEC);
    dof   = CHAIN_NEXT(dof, EL_DOF_VEC);
    bnd   = bnd ? CHAIN_NEXT(bnd, EL_SCHAR_VEC) : NULL;
    B     = B ? COL_CHAIN_NEXT(B, const EL_MATRIX) : NULL;
  } COL_CHAIN_WHILE(A, const EL_MATRIX);
}

/* >>> */

/* <<< variable stride version */

static inline
void __bi_mat_el_vec_dow_scl(REAL a, const EL_MATRIX *A,
			     REAL b, const EL_MATRIX *B,
			     const EL_REAL_VEC *u_h,
			     REAL c, DOF_REAL_VEC_D *f_h,
			     const EL_DOF_VEC *dof,
			     const EL_SCHAR_VEC *bnd)
{
  if (f_h->stride != 1) {
    __bi_mat_el_vec_rdr(a, A, b, B, u_h, c, (DOF_REAL_D_VEC *)f_h, dof, bnd);
  } else {
    __bi_mat_el_vec(a, A, b, B, u_h, c, (DOF_REAL_VEC *)f_h, dof, bnd);
  }
}

__FLATTEN_ATTRIBUTE__
static inline
void bi_mat_el_vec_dow_scl(REAL a,
			   const EL_MATRIX *A,
			   REAL b,
			   const EL_MATRIX *B,
			   const EL_REAL_VEC *u_h,
			   REAL c,
			   DOF_REAL_VEC_D *f_h,
			   const EL_DOF_VEC *dof,
			   const EL_SCHAR_VEC *bnd)
{
  if (A == NULL) {
    A = B;
    a = b;
    B = NULL;
    b = 0.0;
  }
  COL_CHAIN_DO(A, const EL_MATRIX) {
    ROW_CHAIN_DO(A, const EL_MATRIX) {

      __bi_mat_el_vec_dow_scl(a, A, b, B, u_h, c, f_h, dof, bnd);

      B = B ? ROW_CHAIN_NEXT(B, const EL_MATRIX) : NULL;
      u_h = CHAIN_NEXT(u_h, const EL_REAL_VEC);
    } ROW_CHAIN_WHILE(A, const EL_MATRIX);
    f_h   = CHAIN_NEXT(f_h, DOF_REAL_VEC_D);
    dof   = CHAIN_NEXT(dof, EL_DOF_VEC);
    bnd   = bnd ? CHAIN_NEXT(bnd, EL_SCHAR_VEC) : NULL;
    B     = B ? COL_CHAIN_NEXT(B, const EL_MATRIX) : NULL;
  } COL_CHAIN_WHILE(A, const EL_MATRIX);
}

/* >>> */

/* >>> */

/* <<< DOW-DOW version */

/* <<< fixed stride */

/* Multiply an element matrix by an element vector and store the
 * result inside an element vector.
 */
static inline
void __bi_mat_el_vec_d(REAL a, const EL_MATRIX *A,
		       REAL b, const EL_MATRIX *B,
		       const EL_REAL_D_VEC *u_h,
		       REAL c, DOF_REAL_D_VEC *f_h,
		       const EL_DOF_VEC *dof,
		       const EL_SCHAR_VEC *bnd)
{
  int i, j;
  REAL *sum;

  if (A && B) {
#undef MAT_BI_BODY
#define MAT_BI_BODY(PFX1, CC1, C1, S1, TYPE1,		\
		    PFX2, CC2, C2, S2, TYPE2)		\
    for (i = 0; i < A->n_row; i++) {			\
      if (bnd == NULL || bnd->vec[i] < DIRICHLET) {	\
	sum = f_h->vec[dof->vec[i]];			\
	for (j = 0; j < A->n_col; j++) {		\
	  PFX1##PFX2##BIMV_DOW(a, CC1 A->data.S1[i][j],	\
			       b, CC2 B->data.S2[i][j],	\
			       u_h->vec[j],		\
			       c, sum);			\
	}						\
      }							\
    }
    MAT_EMIT_BI_BODY_SWITCH(A->type, B->type);
  } else {
    if (A == NULL) {
      A = B;
      a = b;
      B = NULL;
      b = 0.0;
    }
#undef MAT_BODY
#define MAT_BODY(PFX, CC, C, S, TYPE)					\
    for (i = 0; i < A->n_row; i++) {					\
      if (bnd == NULL || bnd->vec[i] < DIRICHLET) {			\
	sum = f_h->vec[dof->vec[i]];					\
	for (j = 0; j < A->n_col; j++) {				\
	  PFX##GEMV_DOW(a, CC A->data.S[i][j], u_h->vec[j], c, sum);	\
	}								\
      }									\
    }
    MAT_EMIT_BODY_SWITCH(A->type);
  }
}

__FLATTEN_ATTRIBUTE__
static inline
void bi_mat_el_vec_d(REAL a,
		     const EL_MATRIX *A,
		     REAL b,
		     const EL_MATRIX *B,
		     const EL_REAL_D_VEC *u_h,
		     REAL c,
		     DOF_REAL_D_VEC *f_h,
		     const EL_DOF_VEC *dof,
		     const EL_SCHAR_VEC *bnd)
{
  if (A == NULL) {
    A = B;
    a = b;
    B = NULL;
    b = 0.0;
  }
  COL_CHAIN_DO(A, const EL_MATRIX) {
    ROW_CHAIN_DO(A, const EL_MATRIX) {

      __bi_mat_el_vec_d(a, A, b, B, u_h, c, f_h, dof, bnd);

      B = B ? ROW_CHAIN_NEXT(B, const EL_MATRIX) : NULL;
      u_h = CHAIN_NEXT(u_h, const EL_REAL_D_VEC);
    } ROW_CHAIN_WHILE(A, const EL_MATRIX);
    f_h = CHAIN_NEXT(f_h, DOF_REAL_D_VEC);
    dof = CHAIN_NEXT(dof, EL_DOF_VEC);
    bnd = bnd ? CHAIN_NEXT(bnd, EL_SCHAR_VEC) : NULL;
    B   = B ? COL_CHAIN_NEXT(B, const EL_MATRIX) : NULL;
  } COL_CHAIN_WHILE(A, const EL_MATRIX);
}

/* >>> */

/* <<< variable stride version */

static inline
void __bi_mat_el_vec_dow(REAL a, const EL_MATRIX *A,
			 REAL b, const EL_MATRIX *B,
			 const EL_REAL_VEC_D *u_h,
			 REAL c, DOF_REAL_VEC_D *f_h,
			 const EL_DOF_VEC *dof,
			 const EL_SCHAR_VEC *bnd)
{
  if (f_h->stride == 1) {
    __bi_mat_el_vec_scl_dow(a, A, b, B, u_h,
			    c, (DOF_REAL_VEC *)f_h, dof, bnd);
  } else if (u_h->stride == 1) {
    __bi_mat_el_vec_dow_scl(a, A, b, B, (const EL_REAL_VEC *)u_h,
			    c, f_h, dof, bnd);
  } else {
    __bi_mat_el_vec_d(a, A, b, B, (const EL_REAL_D_VEC *)u_h,
		      c, (DOF_REAL_D_VEC *)f_h, dof, bnd);
  }
}

__FLATTEN_ATTRIBUTE__
static inline
void bi_mat_el_vec_dow(REAL a,
		       const EL_MATRIX *A,
		       REAL b,
		       const EL_MATRIX *B,
		       const EL_REAL_VEC_D *u_h,
		       REAL c,
		       DOF_REAL_VEC_D *f_h,
		       const EL_DOF_VEC *dof,
		       const EL_SCHAR_VEC *bnd)
{
  if (A == NULL) {
    A = B;
    a = b;
    B = NULL;
    b = 0.0;
  }
  COL_CHAIN_DO(A, const EL_MATRIX) {
    ROW_CHAIN_DO(A, const EL_MATRIX) {

      __bi_mat_el_vec_dow(a, A, b, B, u_h, c, f_h, dof, bnd);

      B = B ? ROW_CHAIN_NEXT(B, const EL_MATRIX) : NULL;
      u_h = CHAIN_NEXT(u_h, const EL_REAL_VEC_D);
    } ROW_CHAIN_WHILE(A, const EL_MATRIX);
    f_h = CHAIN_NEXT(f_h, DOF_REAL_VEC_D);
    dof = CHAIN_NEXT(dof, EL_DOF_VEC);
    bnd = bnd ? CHAIN_NEXT(bnd, EL_SCHAR_VEC) : NULL;
    B   = B ? COL_CHAIN_NEXT(B, const EL_MATRIX) : NULL;
  } COL_CHAIN_WHILE(A, const EL_MATRIX);
}

/* >>> */

/* >>> */

/* <<< f = c * f + a A u, f += a A u versions */

/* <<< SCL x SCL */

static inline void
gen_mat_el_vec(REAL a, const EL_MATRIX *A, const EL_REAL_VEC *u_h,
	       REAL b, DOF_REAL_VEC *f_h,
	       const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bnd)
{
  bi_mat_el_vec(a, A, 0.0, NULL, u_h, b, f_h, dof, bnd);
}

static inline void
mat_el_vec(REAL a, const EL_MATRIX *A, const EL_REAL_VEC *u_h,
	   DOF_REAL_VEC *f_h, const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bnd)
{
  gen_mat_el_vec(a, A, u_h, 1.0, f_h, dof, bnd);
}

/* >>> */

/* <<< DOW x DOW */

static inline void
gen_mat_el_vec_d(REAL a, const EL_MATRIX *A, const EL_REAL_D_VEC *u_h,
		 REAL b, DOF_REAL_D_VEC *f_h,
		 const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bnd)
{
  bi_mat_el_vec_d(a, A, 0.0, NULL, u_h, b, f_h, dof, bnd);
}

static inline void
mat_el_vec_d(REAL a, const EL_MATRIX *A, const EL_REAL_D_VEC *u_h,
	     DOF_REAL_D_VEC *f_h,
	     const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bnd)
{
  gen_mat_el_vec_d(a, A, u_h, 1.0, f_h, dof, bnd);
}

static inline void
gen_mat_el_vec_dow(REAL a, const EL_MATRIX *A, const EL_REAL_VEC_D *u_h,
		   REAL b, DOF_REAL_VEC_D *f_h,
		   const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bnd)
{
  bi_mat_el_vec_dow(a, A, 0.0, NULL, u_h, b, f_h, dof, bnd);
}

static inline void
mat_el_vec_dow(REAL a, const EL_MATRIX *A, const EL_REAL_VEC_D *u_h,
	       DOF_REAL_VEC_D *f_h,
	       const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bnd)
{
  gen_mat_el_vec_dow(a, A, u_h, 1.0, f_h, dof, bnd);
}

/* >>> */

/* <<< SCL x DOW */

static inline void
gen_mat_el_vec_rrd(REAL a, const EL_MATRIX *A, const EL_REAL_D_VEC *u_h,
		   REAL b, DOF_REAL_VEC *f_h,
		   const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bnd)
{
  bi_mat_el_vec_rrd(a, A, 0.0, NULL, u_h, b, f_h, dof, bnd);
}

static inline void
mat_el_vec_rrd(REAL a, const EL_MATRIX *A, const EL_REAL_D_VEC *u_h,
	       DOF_REAL_VEC *f_h,
	       const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bnd)
{
  gen_mat_el_vec_rrd(a, A, u_h, 1.0, f_h, dof, bnd);
}

static inline void
gen_mat_el_vec_scl_dow(REAL a, const EL_MATRIX *A, const EL_REAL_VEC_D *u_h,
		       REAL b, DOF_REAL_VEC *f_h,
		       const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bnd)
{
  bi_mat_el_vec_scl_dow(a, A, 0.0, NULL, u_h, b, f_h, dof, bnd);
}

static inline void
mat_el_vec_scl_dow(REAL a, const EL_MATRIX *A, const EL_REAL_VEC_D *u_h,
		   DOF_REAL_VEC *f_h,
		   const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bnd)
{
  gen_mat_el_vec_scl_dow(a, A, u_h, 1.0, f_h, dof, bnd);
}

/* >>> */

/* <<< DOW x SCL */

static inline void
gen_mat_el_vec_rdr(REAL a, const EL_MATRIX *A, const EL_REAL_VEC *u_h,
		   REAL b, DOF_REAL_D_VEC *f_h,
		   const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bnd)
{
  bi_mat_el_vec_rdr(a, A, 0.0, NULL, u_h, b, f_h, dof, bnd);
}

static inline void
mat_el_vec_rdr(REAL a, const EL_MATRIX *A, const EL_REAL_VEC *u_h,
	       DOF_REAL_D_VEC *f_h,
	       const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bnd)
{
  gen_mat_el_vec_rdr(a, A, u_h, 1.0, f_h, dof, bnd);
}

static inline void
gen_mat_el_vec_dow_scl(REAL a, const EL_MATRIX *A, const EL_REAL_VEC *u_h,
		       REAL b, DOF_REAL_VEC_D *f_h,
		       const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bnd)
{
  bi_mat_el_vec_dow_scl(a, A, 0.0, NULL, u_h, b, f_h, dof, bnd);
}

static inline void
mat_el_vec_dow_scl(REAL a, const EL_MATRIX *A, const EL_REAL_VEC *u_h,
		   DOF_REAL_VEC_D *f_h,
		   const EL_DOF_VEC *dof, const EL_SCHAR_VEC *bnd)
{
  gen_mat_el_vec_dow_scl(a, A, u_h, 1.0, f_h, dof, bnd);
}

/* >>> */

/* >>> */

/* >>> */

/* <<< some blas functions (incomplete) */

/* <<< axpy, axey, set for EL_MATRIX */

static inline
EL_MATRIX *__el_mat_set(REAL a, EL_MATRIX *res)
{
  int i, j;

#undef MAT_BODY
#define MAT_BODY(PFX, CC, C, S, TYPE)		\
  for (i = 0; i < res->n_row; i++) {		\
    for (j = 0; j < res->n_col; j++) {		\
      PFX##SET_DOW(a, res->data.S[i][j]);	\
    }						\
  }
  MAT_EMIT_BODY_SWITCH(res->type);

  return res;
}

static inline EL_MATRIX *el_mat_set(REAL a, EL_MATRIX *res)
{
  COL_CHAIN_DO(res, EL_MATRIX) {
    ROW_CHAIN_DO(res, EL_MATRIX) {
      __el_mat_set(a, res);
    } ROW_CHAIN_WHILE(res, EL_MATRIX);
  } COL_CHAIN_WHILE(res, EL_MATRIX);

  return res;
}

static inline
EL_MATRIX *__el_mat_axey(REAL a, const EL_MATRIX *A, EL_MATRIX *res)
{
  int i, j;

#undef MAT_TRI_BODY
#define MAT_TRI_BODY(PFX1, CC1, C1, S1, TYPE1,				\
		     PFX2, CC2, C2, S2, TYPE2)				\
  for (i = 0; i < A->n_row; i++) {					\
    for (j = 0; j < A->n_col; j++) {					\
      PFX1##PFX2##AXEY_DOW(a, CC2 A->data.S2[i][j], res->data.S1[i][j]); \
    }									\
  }
  MAT_EMIT_TRI_BODY_SWITCH(A->type, res->type);

  return res;
}

static inline EL_MATRIX *
el_mat_axey(REAL a, const EL_MATRIX *A, EL_MATRIX *res)
{
  COL_CHAIN_DO(res, EL_MATRIX) {
    ROW_CHAIN_DO(res, EL_MATRIX) {
      __el_mat_axey(a, A, res);
      A = ROW_CHAIN_NEXT(A, const EL_MATRIX);
    } ROW_CHAIN_WHILE(res, EL_MATRIX);
    A = COL_CHAIN_NEXT(A, const EL_MATRIX);
  } COL_CHAIN_WHILE(res, EL_MATRIX);

  return res;
}

static inline
EL_MATRIX *__el_mat_axpy(REAL a, const EL_MATRIX *A, EL_MATRIX *res)
{
  int i, j;

#undef MAT_TRI_BODY
#define MAT_TRI_BODY(PFX1, CC1, C1, S1, TYPE1,				\
		     PFX2, CC2, C2, S2, TYPE2)				\
  for (i = 0; i < A->n_row; i++) {					\
    for (j = 0; j < A->n_col; j++) {					\
      PFX1##PFX2##AXPY_DOW(a, CC2 A->data.S2[i][j], res->data.S1[i][j]); \
    }									\
  }
  MAT_EMIT_TRI_BODY_SWITCH(res->type, A->type);

  return res;
}

static inline EL_MATRIX *
el_mat_axpy(REAL a, const EL_MATRIX *A, EL_MATRIX *res)
{
  COL_CHAIN_DO(res, EL_MATRIX) {
    ROW_CHAIN_DO(res, EL_MATRIX) {
      __el_mat_axpy(a, A, res);
      A = ROW_CHAIN_NEXT(A, const EL_MATRIX);
    } ROW_CHAIN_WHILE(res, EL_MATRIX);
    A = COL_CHAIN_NEXT(A, const EL_MATRIX);
  } COL_CHAIN_WHILE(res, EL_MATRIX);

  return res;
}

static inline
EL_MATRIX *M_el_mat_axpby(REAL a, const EL_MATRIX *X,
			  REAL b, const EL_MATRIX *Y, EL_MATRIX *res)
{

  int i, j;

  switch (X->type) {
  case MATENT_REAL_DD:
    switch (Y->type) {
    case MATENT_REAL_DD:
      for (i = 0; i < res->n_row; i++) {
	for (j = 0; j < res->n_col; j++) {
	  MMMAXPBY_DOW(a, (const REAL_D *) X->data.real_dd[i][j],
		       b, (const REAL_D *) Y->data.real_dd[i][j],
		       res->data.real_dd[i][j]);
	}
      }
      break;
    case MATENT_REAL_D:
      for (i = 0; i < res->n_row; i++) {
	for (j = 0; j < res->n_col; j++) {
	  MMDMAXPBY_DOW(a, (const REAL_D *) X->data.real_dd[i][j],
			b, Y->data.real_d[i][j], res->data.real_dd[i][j]);
	}
      }
      break;
    case MATENT_REAL:
      for (i = 0; i < res->n_row; i++) {
	for (j = 0; j < res->n_col; j++) {
	  MMSCMAXPBY_DOW(a, (const REAL_D *) X->data.real_dd[i][j],
			 b, Y->data.real[i][j], res->data.real_dd[i][j]);
	}
      }
      break;
    default:
      ERROR_EXIT("!!! unknown MATENT_TYPE: %d", Y->type);
      break;
    }
    break;
  case MATENT_REAL_D:
    if (Y->type == MATENT_REAL_D) {
      for (i = 0; i < res->n_row; i++) {
	for (j = 0; j < res->n_col; j++) {
	  MDMDMAXPBY_DOW(a, X->data.real_d[i][j],
			 b, Y->data.real_d[i][j], res->data.real_dd[i][j]);
	}
      }
    } else if (Y->type == MATENT_REAL) {
      for (i = 0; i < res->n_row; i++) {
	for (j = 0; j < res->n_col; j++) {
	  MDMSCMAXPBY_DOW(a, X->data.real_d[i][j],
			  b, Y->data.real[i][j], res->data.real_dd[i][j]);
	}
      }
    }
    break;
  case MATENT_REAL:
    if (res->type == MATENT_REAL) {
      for (i = 0; i < res->n_row; i++) {
	for (j = 0; j < res->n_col; j++) {
	  MSCMSCMAXPBY_DOW(a, X->data.real[i][j],
			   b, Y->data.real[i][j], res->data.real_dd[i][j]);
	}
      }
    }
    break;
  default:
    ERROR_EXIT("!!! unknown MATENT_TYPE: %d", X->type);
    break;
  }

  return res;
}

static inline
EL_MATRIX *DM_el_mat_axpby(REAL a, const EL_MATRIX *X,
			   REAL b, const EL_MATRIX *Y, EL_MATRIX *res)
{

  int i, j;

  switch (X->type) {

  case MATENT_REAL_D:
    if (Y->type == MATENT_REAL_D) {
      for (i = 0; i < res->n_row; i++) {
	for (j = 0; j < res->n_col; j++) {
	  DMDMDMAXPBY_DOW(a, X->data.real_d[i][j],
			  b, Y->data.real_d[i][j], res->data.real_d[i][j]);
	}
      }
    } else if (Y->type == MATENT_REAL) {
      for (i = 0; i < res->n_row; i++) {
	for (j = 0; j < res->n_col; j++) {
	  DMDMSCMAXPBY_DOW(a, X->data.real_d[i][j],
			   b, Y->data.real[i][j], res->data.real_d[i][j]);
	}
      }
    }
    break;
  case MATENT_REAL:
    if (Y->type == MATENT_REAL) {
      for (i = 0; i < res->n_row; i++) {
	for (j = 0; j < res->n_col; j++) {
	  DMSCMSCMAXPBY_DOW(a, X->data.real[i][j],
			    b, Y->data.real[i][j], res->data.real_d[i][j]);
	}
      }
    }
    break;
  default:
    ERROR_EXIT("!!! unknown MATENT_TYPE: %d", X->type);
    break;
  }

  return res;
}

static inline
EL_MATRIX *SCM_el_mat_axpby(REAL a, const EL_MATRIX *X,
			    REAL b, const EL_MATRIX *Y, EL_MATRIX *res)
{
  int i, j;

  for (i = 0; i < res->n_row; i++) {
    for (j = 0; j < res->n_col; j++) {
      SCMSCMSCMAXPBY_DOW(a, X->data.real[i][j],
			 b, Y->data.real[i][j], res->data.real[i][j]);
    }
  }

  return res;
}

static inline
EL_MATRIX *__el_mat_axpby(REAL a, const EL_MATRIX *X,
			  REAL b, const EL_MATRIX *Y, EL_MATRIX *res)
{
  FUNCNAME("el_mat_axpby");


  REAL tmp;
  const EL_MATRIX *tmp_mat;


  TEST_EXIT(res,
	    "please get_el_matrix(result,...) bevore using el_mat_axpby!!!\n");
  TEST_EXIT(res->type >= X->type,
	    "MATENT_TYPE of EL_MATRIX *result have to be equal "
	    "or greater then MATENT_TYPE of EL_MATRIX *X!!!\n");
  TEST_EXIT(res->type >= Y->type,
	    "MATENT_TYPE of EL_MATRIX *result have to be equal "
	    "or greater then MATENT_TYPE of EL_MATRIX *Y!!!\n");


  if (Y->type > X->type) {
    tmp_mat = X;
    tmp = a;
    X = Y;
    a = b;
    Y = tmp_mat;
    b = tmp;
  }

  switch (res->type) {
  case MATENT_REAL_DD:
    res = M_el_mat_axpby(a, X, b, Y, res);
    break;
  case MATENT_REAL_D:
    res = DM_el_mat_axpby(a, X, b, Y, res);
    break;
  case MATENT_REAL:
    res = SCM_el_mat_axpby(a, X, b, Y, res);
    break;
  default:
    ERROR_EXIT("!!! unknown MATENT_TYPE: %d", res->type);
    break;
  }

  return (res);

}

static inline EL_MATRIX *
el_mat_axpby(REAL a, const EL_MATRIX *X,
	     REAL b, const EL_MATRIX *Y, EL_MATRIX *res)
{
  COL_CHAIN_DO(res, EL_MATRIX) {
    ROW_CHAIN_DO(res, EL_MATRIX) {
      __el_mat_axpby(a, X, b, Y, res);
      X = ROW_CHAIN_NEXT(X, const EL_MATRIX);
      Y = ROW_CHAIN_NEXT(Y, const EL_MATRIX);
    } ROW_CHAIN_WHILE(res, EL_MATRIX);
    X = COL_CHAIN_NEXT(X, const EL_MATRIX);
    Y = COL_CHAIN_NEXT(Y, const EL_MATRIX);
  } COL_CHAIN_WHILE(res, EL_MATRIX);

  return res;
}

/* >>> */

/* <<< axpy() for element vectors  */

static inline
void __el_real_vec_axpy(REAL a, const EL_REAL_VEC *x, EL_REAL_VEC *y)
{
  FUNCNAME("__el_real_vec_axpy");
  int i;
  DEBUG_TEST_EXIT(x->n_components == y->n_components,
		  "element vector dimensions do not match, x: %d, y: %d.\n",
		  x->n_components, y->n_components);
  for (i = 0; i < x->n_components; i++) {
    y->vec[i] += a * x->vec[i];
  }
}

__FLATTEN_ATTRIBUTE__
static inline
void el_real_vec_axpy(REAL a, const EL_REAL_VEC *x, EL_REAL_VEC *y)
{
  CHAIN_DO(x, const EL_REAL_VEC) {
    __el_real_vec_axpy(a, x, y);
    y = CHAIN_NEXT(y, EL_REAL_VEC);
  } CHAIN_WHILE(x, const EL_REAL_VEC);
}

static inline
void __el_real_d_vec_axpy(REAL a, const EL_REAL_D_VEC *x, EL_REAL_D_VEC *y)
{
  FUNCNAME("__el_real_d_vec_axpy");
  int i;
  DEBUG_TEST_EXIT(x->n_components == y->n_components,
		  "element vector dimensions do not match, x: %d, y: %d.\n",
		  x->n_components, y->n_components);
  for (i = 0; i < x->n_components; i++) {
    AXPY_DOW(a, x->vec[i], y->vec[i]);
  }
}

__FLATTEN_ATTRIBUTE__
static inline
void el_real_d_vec_axpy(REAL a, const EL_REAL_D_VEC *x, EL_REAL_D_VEC *y)
{
  CHAIN_DO(x, const EL_REAL_D_VEC) {
    __el_real_d_vec_axpy(a, x, y);
    y = CHAIN_NEXT(y, EL_REAL_D_VEC);
  } CHAIN_WHILE(x, const EL_REAL_D_VEC);
}

__FLATTEN_ATTRIBUTE__
static inline
void el_real_vec_d_axpy(REAL a, const EL_REAL_VEC_D *x, EL_REAL_VEC_D *y)
{
  CHAIN_DO(x, const EL_REAL_VEC_D) {
    if (x->stride == 1) {
      __el_real_vec_axpy(a, (const EL_REAL_VEC *)x, (EL_REAL_VEC *)y);
    } else {
      __el_real_d_vec_axpy(a, (const EL_REAL_D_VEC *)x, (EL_REAL_D_VEC *)y);
    }
    y = CHAIN_NEXT(y, EL_REAL_VEC_D);
  } CHAIN_WHILE(x, const EL_REAL_VEC_D);
}

/* >>> */

/* <<< axey() for element vectors  */

static inline
void __el_real_vec_axey(REAL a, const EL_REAL_VEC *x, EL_REAL_VEC *y)
{
  FUNCNAME("__el_real_vec_axey");
  int i;
  DEBUG_TEST_EXIT(x->n_components == y->n_components,
		  "element vector dimensions do not match, x: %d, y: %d.\n",
		  x->n_components, y->n_components);
  for (i = 0; i < x->n_components; i++) {
    y->vec[i] = a * x->vec[i];
  }
}

__FLATTEN_ATTRIBUTE__
static inline
void el_real_vec_axey(REAL a, const EL_REAL_VEC *x, EL_REAL_VEC *y)
{
  CHAIN_DO(x, const EL_REAL_VEC) {
    __el_real_vec_axey(a, x, y);
    y = CHAIN_NEXT(y, EL_REAL_VEC);
  } CHAIN_WHILE(x, const EL_REAL_VEC);
}

static inline
void __el_real_d_vec_axey(REAL a, const EL_REAL_D_VEC *x, EL_REAL_D_VEC *y)
{
  FUNCNAME("__el_real_d_vec_axey");
  int i;
  DEBUG_TEST_EXIT(x->n_components == y->n_components,
		  "element vector dimensions do not match, x: %d, y: %d.\n",
		  x->n_components, y->n_components);
  for (i = 0; i < x->n_components; i++) {
    AXEY_DOW(a, x->vec[i], y->vec[i]);
  }
}

__FLATTEN_ATTRIBUTE__
static inline
void el_real_d_vec_axey(REAL a, const EL_REAL_D_VEC *x, EL_REAL_D_VEC *y)
{
  CHAIN_DO(x, const EL_REAL_D_VEC) {
    __el_real_d_vec_axey(a, x, y);
    y = CHAIN_NEXT(y, EL_REAL_D_VEC);
  } CHAIN_WHILE(x, const EL_REAL_D_VEC);
}

__FLATTEN_ATTRIBUTE__
static inline
void el_real_vec_d_axey(REAL a, const EL_REAL_VEC_D *x, EL_REAL_VEC_D *y)
{
  CHAIN_DO(x, const EL_REAL_VEC_D) {
    if (x->stride == 1) {
      __el_real_vec_axey(a, (const EL_REAL_VEC *)x, (EL_REAL_VEC *)y);
    } else {
      __el_real_d_vec_axey(a, (const EL_REAL_D_VEC *)x, (EL_REAL_D_VEC *)y);
    }
    y = CHAIN_NEXT(y, EL_REAL_VEC_D);
  } CHAIN_WHILE(x, const EL_REAL_VEC_D);
}

/* >>> */

/* <<< set() for element vectors  */

static inline
void __el_real_vec_set(REAL a, EL_REAL_VEC *y)
{
  int i;
  for (i = 0; i < y->n_components; i++) {
    y->vec[i] = a;
  }
}

__FLATTEN_ATTRIBUTE__
static inline
void el_real_vec_set(REAL a, EL_REAL_VEC *y)
{
  CHAIN_DO(y, EL_REAL_VEC) {
    __el_real_vec_set(a, y);
  } CHAIN_WHILE(y, EL_REAL_VEC);
}

static inline
void __el_real_d_vec_set(REAL a, EL_REAL_D_VEC *y)
{
  int i;
  for (i = 0; i < y->n_components; i++) {
    SET_DOW(a, y->vec[i]);
  }
}

__FLATTEN_ATTRIBUTE__
static inline
void el_real_d_vec_set(REAL a, EL_REAL_D_VEC *y)
{
  CHAIN_DO(y, EL_REAL_D_VEC) {
    __el_real_d_vec_set(a, y);
  } CHAIN_WHILE(y, EL_REAL_D_VEC);
}

__FLATTEN_ATTRIBUTE__
static inline
void el_real_vec_d_set(REAL a, EL_REAL_VEC_D *y)
{
  CHAIN_DO(y, EL_REAL_VEC_D) {
    if (y->stride == 1) {
      __el_real_vec_set(a, (EL_REAL_VEC *)y);
    } else {
      __el_real_d_vec_set(a, (EL_REAL_D_VEC *)y);
    }
  } CHAIN_WHILE(y, EL_REAL_VEC_D);
}

/* >>> */

/* >>> */

#endif /* _ALBERTA_EL_VEC_H_ */