libalberta2-dev 2.0.1-5 / usr / include / alberta / alberta_inlines.h

#ifndef _ALBERT_INLINES_H_
#define _ALBERT_INLINES_H_

/*--------------------------------------------------------------------------*/
/* ALBERT:   an Adaptive multi Level finite element toolbox using           */
/*           Bisectioning refinement and Error control by Residual          */
/*           Techniques                                                     */
/*                                                                          */
/* file: albert_inlines.h                                                   */
/*                                                                          */
/*                                                                          */
/* description: Blas-like inline functions for REAL_Ds and REAL_DDs         */
/*                                                                          */
/*--------------------------------------------------------------------------*/
/*                                                                          */
/*  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                                    */
/*                                                                          */
/*             Claus-Justus Heine                                           */
/*             Abteilung fuer Angewandte Mathematik                         */
/*             Albert-Ludwigs-Universitaet Freiburg                         */
/*             Hermann-Herder-Str. 10                                       */
/*             D-79104 Freiburg im Breisgau, Germany                        */
/*                                                                          */
/*  http://www.mathematik.uni-freiburg.de/IAM/ALBERT                        */
/*                                                                          */
/*  (c) by A. Schmidt, K.G. Siebert, C.-J. Heine (1996-2003)                */
/*                                                                          */
/*--------------------------------------------------------------------------*/

#include "alberta.h" /* essentially a no-op when included from alberta.h */

#ifndef DIM_OF_WORLD
# error Need to know the dimension of the World :)
#endif

/* multiple invocations of macro-arguments can be harmful, if the macro
 * argument is, e.g., a function-call.
 *
 * NOTE: as DIM_OF_WORLD is a constant, the C-compiler should unroll all 
 * loops when compiling with optimizations, so there should be no need
 * for hand-unrolling, except in some simple 1D cases.
 *
 * Also, all modern compilers do function inlining, so the
 * function-call  over-head is _not_ a problem.
 *
 * Note: the function may be nested, they return the address of the
 * _modified_ operand. So AXPY(a, AX(b, x), y) is valid.
 */

/* The following functions are defined here:
 *
 *   AX(a, x)              --  x *= a (alias SCAL_DOW is also defined)
 *   AXEY(a, x, y)         --  y  = a x
 *   AXPBY(a, x, b, y, z)  --  z  = a x + by
 *   AXPBYP(a, x, b, y, z) --  z += a x + by
 *   AXPY(a, x, y)         --  y += a x
 *   COPY(src, dst)        --  dst := src
 *   DIST(x, y)            --  sqrt(DST2(x, y))
 *   DST2(x, y)            --  SCP(x-y, x-y)
 *   NRM2(x)               --  SCP(x, x)
 *   NORM(x)               --  sqrt(NRM2(x))
 *   MTV(m, v, b)          --  b += m^t v
 *   MV(m, v, b)           --  b += m v   (m is a matrix)
 *   MDIV(m, v, b)         --  scale v by the inverse of the diagonal -> b
 *   SCP(x, y)             --  <x, y>
 *   SET(val, x)           --  x[i] = val, i=1, ..., DOW
 *   WEDGE(x, y, n)        --  n  = x /\ y             in 3D
 *   WEDGE(x, y)           --  x0 * y1 - x1 * y0       in 2D
 *
 * The actual function named is generated by adding a _DOW() suffix.
 *
 * Prefix        Version
 * none          REAL_D
 * M             REAL_DD
 * DM            diagonal matrix, diagonal stored in REAL_D vector
 * SM            symmetric matrix, data type REAL_DDS (albert.h)
 *
 * Further:
 * Macros EXPAND and FORMAT (with named pre- and suffixes) for easier 
 * print-out of REAL_D and REAL_DD, use like this:
 *
 * printf("text"MFORMAT_DOW"more text\n", MEXPAND_DOW(m));
 *
 */

#ifndef C_CONST_IS_BRAIN_DAMAGED
# define C_CONST_IS_BRAIN_DAMAGED 1
#endif
#if C_CONST_IS_BRAIN_DAMAGED
# define CCIBD_CONST /* */
#else
# define CCIBD_CONST const /* actually, this will never happen ... :( */
#endif

#define SCAL_DOW(a, x) AX_DOW(a, x)
static inline REAL *AX_DOW(REAL a, REAL_D x)
{
  int i;

  for (i = 0; i < DIM_OF_WORLD; i++) {
    x[i] *= a;
  }
  return x;
}

#define MSCAL_DOW(a, m) MAX_DOW(a, m)
static inline REAL_D *MAX_DOW(REAL a, REAL_DD m)
{
  int i;

  for (i = 0; i < DIM_OF_WORLD; i++) {
    AX_DOW(a, m[i]);
  }
  return m;
}

#define DMSCAL_DOW(a, m) DMAX_DOW(a, m)
#define DMAX_DOW(a, m) AX_DOW(a, m)

#define SMSCAL_DOW(a, m) SMAX_DOW(a, m)
static inline REAL_DDS *SMAX_DOW(REAL a, REAL_DDS *m)
{
  int i;

  AX_DOW(a, m->row0);
  for (i = 0; i < DIM_OF_WORLD - 1; i++) {
    m->row1[i] *= a;
  }
  for (i = 0; i < DIM_OF_WORLD - 2; i++) {
    m->row2[i] *= a;
  }
  return m;
}

static inline REAL *AXEY_DOW(REAL a, const REAL_D x, REAL_D y)
{
  int i;

  for (i = 0; i < DIM_OF_WORLD; i++) {
    y[i] = a * x[i];
  }
  return y;
}

static inline REAL_D *MAXEY_DOW(REAL a, CCIBD_CONST REAL_DD x, REAL_DD y)
{
  int i;

  for (i = 0; i < DIM_OF_WORLD; i++) {
    AXEY_DOW(a, x[i], y[i]);
  }
  return y;
}

#define DMAXEY_DOW(a, x, y) AXEY_DOW(a, x, y)

static inline REAL_DDS *SMAXEY_DOW(REAL a, const REAL_DDS *x, REAL_DDS *y)
{
  int i;

  AXEY_DOW(a, x->row0, y->row0);
  for (i = 0; i < DIM_OF_WORLD - 1; i++) {
    y->row1[i] = a*x->row1[i];
  }
  for (i = 0; i < DIM_OF_WORLD - 2; i++) {
    y->row2[i] = a*x->row2[i];
  }
  return y;
}

static inline REAL *AXPY_DOW(REAL a, const REAL_D x, REAL_D y)
{
  int i;

  for (i = 0; i < DIM_OF_WORLD; i++) {
    y[i] += a * x[i];
  }
  return y;
}

static inline REAL_D *MAXPY_DOW(REAL a, CCIBD_CONST REAL_DD x, REAL_DD y)
{
  int i;

  for (i = 0; i < DIM_OF_WORLD; i++) {
    AXPY_DOW(a, x[i], y[i]);
  }
  return y;
}

/* same as above, but add the transposed matrix to y */
static inline REAL_D *MAXTPY_DOW(REAL a, CCIBD_CONST REAL_DD x, REAL_DD y)
{
  int i, j;

  for (i = 0; i < DIM_OF_WORLD; i++) {
    for (j = 0; j < DIM_OF_WORLD; j++) {
      y[i][j] += a*x[j][i];
    }
  }
  return y;
}

#define DMAXPY_DOW(a, x, y) AXPY_DOW(a, x, y)
#define DMAXTPY_DOW(a, x, y) AXPY_DOW(a, x, y) /* transpose of diagonal matrix :) */

static inline REAL_DDS *SMAXPY_DOW(REAL a, const REAL_DDS *x, REAL_DDS *y)
{
  int i;

  AXPY_DOW(a, x->row0, y->row0);
  for (i = 0; i < DIM_OF_WORLD - 1; i++) {
    y->row1[i] += a*x->row1[i];
  }
  for (i = 0; i < DIM_OF_WORLD - 2; i++) {
    y->row2[i] += a*x->row2[i];
  }
  return y;
}

#define SMAXTPY_DOW(a, x, y) SMAXPY_DOW(a, x, y) /* transpose of symmetric matrix :) */

static inline REAL *AXPBY_DOW(REAL a, const REAL_D x, REAL b, const REAL_D y,
			      REAL_D z)
{
  int i;

  for (i = 0; i < DIM_OF_WORLD; i++) {
    z[i] = b*y[i] + a * x[i];
  }
  return z;
}

static inline REAL_D *MAXPBY_DOW(REAL a, CCIBD_CONST REAL_DD x,
				 REAL b, CCIBD_CONST REAL_DD y,
				 REAL_DD z)
{
  int i;

  for (i = 0; i < DIM_OF_WORLD; i++) {
    AXPBY_DOW(a, x[i], b, y[i], z[i]);
  }
  return z;
}

#define DMAXPBY_DOW(a, x, b, y, z) AXPBY_DOW(a, x, b, y, z)

static inline REAL_DDS *SMAXPBY_DOW(REAL a, const REAL_DDS *x,
				    REAL b, const REAL_DDS *y,
				    REAL_DDS *z)
{
  int i;

  AXPBY_DOW(a, x->row0, b, y->row0, z->row0);
  for (i = 0; i < DIM_OF_WORLD - 1; i++) {
    z->row1[i] = a*x->row1[i] + b*y->row1[i];
  }
  for (i = 0; i < DIM_OF_WORLD - 2; i++) {
    z->row2[i] = a*x->row2[i] + b*y->row2[i];
  }
  return z;
}

static inline REAL *AXPBYP_DOW(REAL a, const REAL_D x, REAL b, const REAL_D y,
			       REAL_D z)
{
  int i;

  for (i = 0; i < DIM_OF_WORLD; i++) {
    z[i] += b*y[i] + a * x[i];
  }
  return z;
}

static inline REAL_D *MAXPBYP_DOW(REAL a, CCIBD_CONST REAL_DD x,
				  REAL b, CCIBD_CONST REAL_DD y,
				  REAL_DD z)
{
  int i;

  for (i = 0; i < DIM_OF_WORLD; i++) {
    AXPBYP_DOW(a, x[i], b, y[i], z[i]);
  }
  return z;
}

#define DMAXPBYP_DOW(a, x, b, y, z) AXPBYP_DOW(a, x, b, y, z)

static inline REAL_DDS *SMAXPBYP_DOW(REAL a, const REAL_DDS *x,
				     REAL b, const REAL_DDS *y,
				     REAL_DDS *z)
{
  int i;

  AXPBYP_DOW(a, x->row0, b, y->row0, z->row0);
  for (i = 0; i < DIM_OF_WORLD - 1; i++) {
    z->row1[i] += a*x->row1[i] + b*y->row1[i];
  }
  for (i = 0; i < DIM_OF_WORLD - 2; i++) {
    z->row2[i] += a*x->row2[i] + b*y->row2[i];
  }
  return z;
}

static inline REAL *COPY_DOW(const REAL_D x, REAL_D y)
{
  memcpy(y, x, sizeof(REAL_D));
  return y;
}

static inline REAL_D *MCOPY_DOW(CCIBD_CONST REAL_DD x, REAL_DD y)
{
  memcpy(y, x, sizeof(REAL_DD));
  return y;
}

#define DMCOPY_DOW(src, dst) COPY_DOW(src, dst)

static inline REAL_DDS *SMCOPY_DOW(CCIBD_CONST REAL_DDS *x, REAL_DDS *y)
{
  memcpy(y, x, sizeof(REAL_DDS));
  return y;
}


static inline REAL DST2_DOW(const REAL_D x, const REAL_D y)
{
# if DIM_OF_WORLD == 1
  return SQR(ABS(x[0] - y[0]));
# else
  int i;
  REAL accu;

  accu = SQR(x[0] - y[0]);
  for (i = 1; i < DIM_OF_WORLD; i++) {
    accu += SQR(x[i] - y[i]);
  }
  return accu;
# endif
}

static inline REAL MDST2_DOW(CCIBD_CONST REAL_DD a, CCIBD_CONST REAL_DD b)
{
  int  i;
  REAL res = 0.0;

  for (i = 0; i < DIM_OF_WORLD; i++) {
    res += DST2_DOW(a[i], b[i]);
  }
  return res;
}

#define DMDST2_DOW(x, y) DST2_DOW(x, y)

static inline REAL SMDST2_DOW(const REAL_DDS *a, const REAL_DDS *b)
{
  int  i;
  REAL res = 0.0;

  res = SQR(a->row0[0] - b->row0[0]);
  for (i = 1; i < DIM_OF_WORLD; i++) {
    res += 2.0*SQR(a->row0[i] - b->row0[i]);
  }
#if DIM_OF_WORLD > 1
  res += SQR(a->row1[0] - b->row1[0]);
  for (i = 1; i < DIM_OF_WORLD-1; i++) {
    res += 2.0*SQR(a->row1[i] - b->row1[i]);
  }
#endif
#if DIM_OF_WORLD > 2
  res += SQR(a->row2[0] - b->row2[0]);
  for (i = 1; i < DIM_OF_WORLD-2; i++) {
    res += 2.0*SQR(a->row2[i] - b->row2[i]);
  }
#endif
  return res;
}

static inline REAL NRM2_DOW(const REAL_D x)
{
  int i;
  REAL accu;

  accu = SQR(x[0]);
  for (i = 1; i < DIM_OF_WORLD; i++) {
    accu += SQR(x[i]);
  }
  return accu;
}

static inline REAL MNRM2_DOW(CCIBD_CONST REAL_DD m)
{
  int  i;
  REAL res = 0.0;

  for (i = 0; i < DIM_OF_WORLD; i++) {
    res += NRM2_DOW(m[i]);
  }
  return res;
}

#define DMNRM2_DOW(m) NRM2_DOW(x)

static inline REAL SMNRM2_DOW(const REAL_DDS *a)
{
  int  i;
  REAL res = 0.0;

  res = SQR(a->row0[0]);
  for (i = 1; i < DIM_OF_WORLD; i++) {
    res += 2.0*SQR(a->row0[i]);
  }
#if DIM_OF_WORLD > 1
  res += SQR(a->row1[0]);
  for (i = 1; i < DIM_OF_WORLD-1; i++) {
    res += 2.0*SQR(a->row1[i]);
  }
#endif
#if DIM_OF_WORLD > 2
  res += SQR(a->row2[0]);
  for (i = 1; i < DIM_OF_WORLD-2; i++) {
    res += 2.0*SQR(a->row2[i]);
  }
#endif
  return res;
}

static inline REAL SCP_DOW(const REAL_D x, const REAL_D y)
{
  REAL res;
  int  i;

  res = x[0] * y[0];
  for (i = 1; i < DIM_OF_WORLD; i++) {
    res += x[i]*y[i];
    
  }
  return res;
}

/* symmetric case, no need to load all values, should result in faster
 * code (less registers needed)
 */
static inline REAL *SMV_DOW(const REAL_DDS *m, const REAL_D v, REAL_D b)
{
  b[0] += SCP_DOW(m->row0, v);
#if DIM_OF_WORLD == 2
  b[1] +=  m->row0[1]*v[0] + m->row1[0]*v[1];
#endif
#if DIM_OF_WORLD == 3
  b[1] +=  m->row0[1]*v[0] + m->row1[0]*v[1] + m->row1[1]*v[2];
  b[2] +=  m->row0[2]*v[0] + m->row1[1]*v[1] + m->row2[0]*v[2];
#endif
  return b;
}

static inline REAL *MTV_DOW(CCIBD_CONST REAL_DD m, const REAL_D v, REAL_D b)
{
  int i, j;

  for (i = 0; i < DIM_OF_WORLD; i++) {
    for (j = 0; j < DIM_OF_WORLD; j++) {
      b[i] += m[j][i] * v[j];
    }
  }
  return b;
}

#define SMTV_DOW(m, v, b) SMV_DOW(m, v, b)
#define DMTV_DOW(m, v, b) DMV_DOW(m, v, b)

static inline REAL *MDIV_DOW(CCIBD_CONST REAL_DD m, const REAL_D v, REAL_D b)
{
  int i;

  for (i = 0; i < DIM_OF_WORLD; i++) {
    b[i] = v[i] / m[i][i];
  }
  return b;
}

static inline REAL *SMDIV_DOW(const REAL_DDS *m, const REAL_D v, REAL_D b)
{
  b[0] = v[0] / m->row0[0];
#if DIM_OF_WORLD == 2
  b[1] = v[1] / m->row1[0];
#endif
#if DIM_OF_WORLD == 3
  b[1] = v[2] / m->row2[0];
#endif
  return b;
}

static inline REAL *DMDIV_DOW(const REAL_D m, const REAL_D y, REAL_D r)
{
  int i;
  for (i = 0; i < DIM_OF_WORLD; i++) {
    r[i] = y[i] / m[i];
  }
  return r;
}

static inline REAL *DMV_DOW(const REAL_D x, const REAL_D y, REAL_D r)
{
  int i;
  for (i = 0; i < DIM_OF_WORLD; i++) {
    r[i] += x[i]*y[i];
  }
  return r;
}

static inline REAL *MV_DOW(CCIBD_CONST REAL_DD m, const REAL_D v, REAL_D b)
{
  int i;

  for (i = 0; i < DIM_OF_WORLD; i++) {
    b[i] += SCP_DOW(m[i], v);
  }
  return b;
}

static inline REAL *MGEMV_DOW(REAL a, CCIBD_CONST REAL_DD m,
			      const REAL_D v, REAL beta,
			      REAL_D b)
{
  int i;

  for (i = 0; i < DIM_OF_WORLD; i++) {
    b[i] = beta*b[i] + a * SCP_DOW(m[i], v);
  }
  return b;
}

static inline REAL *SMGEMV_DOW(REAL a, const REAL_DDS *m,
			       const REAL_D v, REAL beta,
			       REAL_D b)
{
  b[0] = beta*b[0] + a*SCP_DOW(m->row0, v);
#if DIM_OF_WORLD == 2
  b[1] = beta*b[1] + a*(m->row0[1]*v[0] + m->row1[0]*v[1]);
#endif
#if DIM_OF_WORLD == 3
  b[1] = beta*b[1] + a*(m->row0[1]*v[0] + m->row1[0]*v[1] + m->row1[1]*v[2]);
  b[2] = beta*b[2] + a*(m->row0[2]*v[0] + m->row1[1]*v[1] + m->row2[0]*v[2]);
#endif
  return b;
}

static inline REAL *DMGEMV_DOW(REAL a, const REAL_D x, const REAL_D y,
			       REAL beta, REAL_D r)
{
  int i;
  for (i = 0; i < DIM_OF_WORLD; i++) {
    r[i] = beta*r[i] + a*x[i]*y[i];
  }
  return r;
}

static inline REAL MSCP_DOW(CCIBD_CONST REAL_DD x, CCIBD_CONST REAL_DD y)
{
  REAL res;
  int  i;

  res = SCP_DOW(x[0], y[0]);
  for (i = 1; i < DIM_OF_WORLD; i++) {
    res += SCP_DOW(x[i], y[i]);    
  }
  return res;
}

#define DMSCP_DOW(x, y) SCP_DOW(x, y)

static inline REAL SMSCP_DOW(const REAL_DDS *x, const REAL_DDS *y)
{
#if DIM_OF_WORLD > 1
  int i;
#endif
  REAL res;

  res = SCP_DOW(x->row0, y->row1);
#if DIM_OF_WORLD > 1
  res += x->row1[0]*y->row1[0];
  for (i = 1; i < DIM_OF_WORLD - 1; i++) {
    res += 2.0*x->row1[i]*y->row1[i];
  }
#endif
#if DIM_OF_WORLD > 2
  res += x->row2[0]*y->row2[0];
  for (i = 1; i < DIM_OF_WORLD - 2; i++) {
    res += 2.0*x->row2[i]*y->row2[i];
  }
#endif
  return res;
}

static inline REAL *SET_DOW(REAL val, REAL_D x)
{
  int i;
  for (i = 0; i < DIM_OF_WORLD; i++) {
    x[i] = val;
  }
  return x;
}

static inline REAL_D *MSET_DOW(REAL val, REAL_DD m)
{
  int i, j;

  for (i = 0; i < DIM_OF_WORLD; i++) {
    m[i][i] = val;
    for (j = i+1; j < DIM_OF_WORLD; j++) {
      m[j][i] = m[i][j] = 0.0;
    }
  }
  return m;
}

#define DMSET_DOW(val, m) SET_DOW(val, m)

static inline REAL_DDS *SMSET_DOW(REAL val, REAL_DDS *m)
{
  int i;

  m->row0[0] = val;
  for (i = 1; i < DIM_OF_WORLD; i++) {
    m->row0[i] = 0.0;
  }
#if DIM_OF_WORLD > 1
  m->row1[0] = val;
  for (i = 1; i < DIM_OF_WORLD - 1; i++) {
    m->row1[i] = val;
  }
#endif
#if DIM_OF_WORLD > 1
  m->row2[0] = val;
  for (i = 1; i < DIM_OF_WORLD - 2; i++) {
    m->row2[i] = val;
  }
#endif

  return m;
}

#if DIM_OF_WORLD == 2
static inline REAL WEDGE_DOW(const REAL_D a, const REAL_D b)
{
  return a[0]*b[1] - a[1]*b[0];
}
#endif

#if DIM_OF_WORLD == 3
static inline REAL *WEDGE_DOW(const REAL_D a, const REAL_D b, REAL_D r)
{
  r[0] = a[1]*b[2] - a[2]*b[1];
  r[1] = a[2]*b[0] - a[0]*b[2];
  r[2] = a[0]*b[1] - a[1]*b[0];
  return r;
}
#endif

#define MAT_SWITCH_TYPE(type, body_f, body_s, body_d)		\
switch (type) {							\
 case dowbm_full: body_f; break;				\
 case dowbm_symm: body_s; break;				\
 case dowbm_diag: body_d; break;				\
 default: ERROR_EXIT("Unknown DOWBM_TYPE (%d)\n", type);	\
}

/* BODY(F, CAST, PRE, SUF) is supposed to be a "multiplex" macro where
 * BLAS routines are accessed via F##AXPY(..., CAST PRE var##SUF, ...)
 */
#define MAT_EMIT_BODY_SWITCH(type)			\
	MAT_SWITCH_TYPE(type,				\
			MAT_BODY(M, (REAL_D *),, full),	\
			MAT_BODY(SM, ,& , symm),	\
			MAT_BODY(DM, , , diag))

/* defines where only DOW == 1 plays a special role */
# if DIM_OF_WORLD == 1
#  define DIST_DOW(x,y)    ABS((x)[0]-(y)[0])
#  define NORM_DOW(x)      ABS((x)[0])
#  define MNRM_DOW(m)      ABS((m)[0][0])
#  define DMNRM_DOW(m)     NRM_DOW(m)
#  define SMNRM_DOW(m)     ABS((m)->row0[0])
#  define MDIST_DOW(a,b)   ABS((a)[0][0] - (b)[0][0])
#  define DMDIST_DOW(a,b)  DIST_DOW(a, b)
#  define SDMDIST_DOW(a,b) ABS((a)->row0[0] - (b)->row0[0])
# else
#  define NORM_DOW(x)     sqrt(NRM2_DOW(x))
#  define DIST_DOW(x,y)   sqrt(DST2_DOW(x, y))
#  define MNORM_DOW(m)    sqrt(MNRM2_DOW(m))
#  define DMNORM_DOW(m)   sqrt(DMNRM2_DOW(m))
#  define SMNORM_DOW(m)   sqrt(DMNRM2_DOW(m))
#  define MDIST_DOW(a,b)  sqrt(MDST2_DOW(a, b))
#  define DMDIST_DOW(a,b) sqrt(DMDST2_DOW(a, b))
#  define SMDIST_DOW(a,b) sqrt(SMDST2_DOW(a, b))
# endif

/* defines different for all DOWs */
# if   DIM_OF_WORLD == 1
#  define REAL_DDS_ENT(a,i,j) (a).row0[0]
#  define EXPAND_DOW(x)  (x)[0]
#  define FORMAT_DOW     "%10.5le"
#  define MEXPAND_DOW(m) (m)[0][0]
#  define MFORMAT_DOW    FORMAT_DOW
#  define SMEXPAND_DOW(m) (m)->row0[0]
#  define SMFORMAT_DOW    FORMAT_DOW
#  define DMEXPAND_DOW(m) EXPAND_DOW(m)
#  define DMFORMAT_DOW    FORMAT_DOW
# elif DIM_OF_WORLD == 2
#  define REAL_DDS_ENT(a,i,j) ((i) == 0 ? (a).row0[j] : (a).row1[j-1])
#  define EXPAND_DOW(x)  (x)[0], (x)[1]
#  define FORMAT_DOW     "[%10.5le, %10.5le]"
#  define MEXPAND_DOW(m) (m)[0][0], (m)[0][1], (m)[1][0], (m)[1][1]
#  define MFORMAT_DOW    "["FORMAT_DOW", "FORMAT_DOW"]"
#  define DMEXPAND_DOW(m) EXPAND_DOW(m)
#  define DMFORMAT_DOW    FORMAT_DOW
#  define SMEXPAND_DOW(m) (m)->row0[0], (m)->row0[1], (m)->row0[1], (m)->row1[0]
#  define SMFORMAT_DOW    MFORMAT_DOW
# elif DIM_OF_WORLD == 3
#  define REAL_DDS_ENT(a,i,j) \
  ((i) == 0 ? (a).row0[j] : ((i) == 1 ? (a).row1[j-1] : (a).row2[0]))
#  define EXPAND_DOW(x)  (x)[0], (x)[1], (x)[2]
#  define FORMAT_DOW     "[%10.5le, %10.5le, %10.5le]"
#  define MEXPAND_DOW(m) \
  (m)[0][0], (m)[0][1], (m)[0][2], \
  (m)[1][0], (m)[1][1], (m)[1][2], \
  (m)[2][0], (m)[2][1], (m)[2][2]
#  define MFORMAT_DOW    "["FORMAT_DOW", "FORMAT_DOW", "FORMAT_DOW"]"
#  define DMEXPAND_DOW(m) EXPAND_DOW(m)
#  define DMFORMAT_DOW    FORMAT_DOW
#  define SMEXPAND_DOW(m)			\
  (m)->row0[0], (m)->row0[1], (m)->row0[2],	\
  (m)->row0[1], (m)->row1[0], (m)->row1[1],	\
  (m)->row0[2], (m)->row1[1], (m)->row2[0]
#  define SMFORMAT_DOW    MFORMAT_DOW
# endif

#endif /* _ALBERT_INLINES_H_ */