source: proiecte/NBody/Tree codes/vectmath.h @ 170

/****************************************************************************/
/* VECTMATH.H: include file for vector/matrix operations.                   */
/* Copyright (c) 1999 by Joshua E. Barnes, Tokyo, JAPAN.                    */
/****************************************************************************/

#ifndef _vectmath_h
#define _vectmath_h

#include "vectdefs.h"

/*
 * Vector operations.
 */

#define CLRV(v)                 /* CLeaR Vector */                      \
{                                                                       \
    int _i;                                                             \
    for (_i = 0; _i < NDIM; _i++)                                       \
        (v)[_i] = 0.0;                                                  \
}

#define UNITV(v,j)              /* UNIT Vector */                       \
{                                                                       \
    int _i;                                                             \
    for (_i = 0; _i < NDIM; _i++)                                       \
        (v)[_i] = (_i == (j) ? 1.0 : 0.0);                              \
}

#define SETV(v,u)               /* SET Vector */                        \
{                                                                       \
    int _i;                                                             \
    for (_i = 0; _i < NDIM; _i++)                                       \
        (v)[_i] = (u)[_i];                                              \
}

#if defined(THREEDIM)

#define ADDV(v,u,w)             /* ADD Vector */                        \
{                                                                       \
    (v)[0] = (u)[0] + (w)[0];                                           \
    (v)[1] = (u)[1] + (w)[1];                                           \
    (v)[2] = (u)[2] + (w)[2];                                           \
}

#define SUBV(v,u,w)             /* SUBtract Vector */                   \
{                                                                       \
    (v)[0] = (u)[0] - (w)[0];                                           \
    (v)[1] = (u)[1] - (w)[1];                                           \
    (v)[2] = (u)[2] - (w)[2];                                           \
}

#define MULVS(v,u,s)            /* MULtiply Vector by Scalar */         \
{                                                                       \
    (v)[0] = (u)[0] * s;                                                \
    (v)[1] = (u)[1] * s;                                                \
    (v)[2] = (u)[2] * s;                                                \
}

#else

#define ADDV(v,u,w)             /* ADD Vector */                        \
{                                                                       \
    int _i;                                                             \
    for (_i = 0; _i < NDIM; _i++)                                       \
        (v)[_i] = (u)[_i] + (w)[_i];                                    \
}

#define SUBV(v,u,w)             /* SUBtract Vector */                   \
{                                                                       \
    int _i;                                                             \
    for (_i = 0; _i < NDIM; _i++)                                       \
        (v)[_i] = (u)[_i] - (w)[_i];                                    \
}

#define MULVS(v,u,s)            /* MULtiply Vector by Scalar */         \
{                                                                       \
    int _i;                                                             \
    for (_i = 0; _i < NDIM; _i++)                                       \
        (v)[_i] = (u)[_i] * (s);                                        \
}

#endif

#define DIVVS(v,u,s)            /* DIVide Vector by Scalar */           \
{                                                                       \
    int _i;                                                             \
    for (_i = 0; _i < NDIM; _i++)                                       \
        (v)[_i] = (u)[_i] / (s);                                        \
}

#if defined(THREEDIM)

#define DOTVP(s,v,u)            /* DOT Vector Product */                \
{                                                                       \
    (s) = (v)[0]*(u)[0] + (v)[1]*(u)[1] + (v)[2]*(u)[2];                \
}

#else

#define DOTVP(s,v,u)            /* DOT Vector Product */                \
{                                                                       \
    int _i;                                                             \
    (s) = 0.0;                                                          \
    for (_i = 0; _i < NDIM; _i++)                                       \
        (s) += (v)[_i] * (u)[_i];                                       \
}

#endif

#define ABSV(s,v)               /* ABSolute value of a Vector */        \
{                                                                       \
    real _tmp;                                                          \
    int _i;                                                             \
    _tmp = 0.0;                                                         \
    for (_i = 0; _i < NDIM; _i++)                                       \
        _tmp += (v)[_i] * (v)[_i];                                      \
    (s) = rsqrt(_tmp);                                                  \
}

#define DISTV(s,u,v)            /* DISTance between Vectors */          \
{                                                                       \
    real _tmp;                                                          \
    int _i;                                                             \
    _tmp = 0.0;                                                         \
    for (_i = 0; _i < NDIM; _i++)                                       \
        _tmp += ((u)[_i]-(v)[_i]) * ((u)[_i]-(v)[_i]);                  \
    (s) = rsqrt(_tmp);                                                  \
}

#if defined(TWODIM)

#define CROSSVP(s,v,u)          /* CROSS Vector Product */              \
{                                                                       \
    (s) = (v)[0]*(u)[1] - (v)[1]*(u)[0];                                \
}

#endif

#if defined(THREEDIM)

#define CROSSVP(v,u,w)          /* CROSS Vector Product */              \
{                                                                       \
    (v)[0] = (u)[1]*(w)[2] - (u)[2]*(w)[1];                             \
    (v)[1] = (u)[2]*(w)[0] - (u)[0]*(w)[2];                             \
    (v)[2] = (u)[0]*(w)[1] - (u)[1]*(w)[0];                             \
}

#endif

/*
 * Matrix operations.
 */

#define CLRM(p)                 /* CLeaR Matrix */                      \
{                                                                       \
    int _i, _j;                                                         \
    for (_i = 0; _i < NDIM; _i++)                                       \
        for (_j = 0; _j < NDIM; _j++)                                   \
            (p)[_i][_j] = 0.0;                                          \
}

#define SETMI(p)                /* SET Matrix to Identity */            \
{                                                                       \
    int _i, _j;                                                         \
    for (_i = 0; _i < NDIM; _i++)                                       \
        for (_j = 0; _j < NDIM; _j++)                                   \
            (p)[_i][_j] = (_i == _j ? 1.0 : 0.0);                       \
}

#define SETM(p,q)               /* SET Matrix */                        \
{                                                                       \
    int _i, _j;                                                         \
    for (_i = 0; _i < NDIM; _i++)                                       \
        for (_j = 0; _j < NDIM; _j++)                                   \
            (p)[_i][_j] = (q)[_i][_j];                                  \
}

#define TRANM(p,q)              /* TRANspose Matrix */                  \
{                                                                       \
    int _i, _j;                                                         \
    for (_i = 0; _i < NDIM; _i++)                                       \
        for (_j = 0; _j < NDIM; _j++)                                   \
            (p)[_i][_j] = (q)[_j][_i];                                  \
}

#define ADDM(p,q,r)             /* ADD Matrix */                        \
{                                                                       \
    int _i, _j;                                                         \
    for (_i = 0; _i < NDIM; _i++)                                       \
        for (_j = 0; _j < NDIM; _j++)                                   \
            (p)[_i][_j] = (q)[_i][_j] + (r)[_i][_j];                    \
}

#define SUBM(p,q,r)             /* SUBtract Matrix */                   \
{                                                                       \
    int _i, _j;                                                         \
    for (_i = 0; _i < NDIM; _i++)                                       \
        for (_j = 0; _j < NDIM; _j++)                                   \
            (p)[_i][_j] = (q)[_i][_j] - (r)[_i][_j];                    \
}

#define MULM(p,q,r)             /* Multiply Matrix */                   \
{                                                                       \
    int _i, _j, _k;                                                     \
    for (_i = 0; _i < NDIM; _i++)                                       \
        for (_j = 0; _j < NDIM; _j++) {                                 \
            (p)[_i][_j] = 0.0;                                          \
            for (_k = 0; _k < NDIM; _k++)                               \
                (p)[_i][_j] += (q)[_i][_k] * (r)[_k][_j];               \
        }                                                               \
}

#define MULMS(p,q,s)            /* MULtiply Matrix by Scalar */         \
{                                                                       \
    int _i, _j;                                                         \
    for (_i = 0; _i < NDIM; _i++)                                       \
        for (_j = 0; _j < NDIM; _j++)                                   \
            (p)[_i][_j] = (q)[_i][_j] * (s);                            \
}

#define DIVMS(p,q,s)            /* DIVide Matrix by Scalar */           \
{                                                                       \
    int _i, _j;                                                         \
    for (_i = 0; _i < NDIM; _i++)                                       \
        for (_j = 0; _j < NDIM; _j++)                                   \
            (p)[_i][_j] = (q)[_i][_j] / (s);                            \
}

#define MULMV(v,p,u)            /* MULtiply Matrix by Vector */         \
{                                                                       \
    int _i, _j;                                                         \
    for (_i = 0; _i < NDIM; _i++) {                                     \
        (v)[_i] = 0.0;                                                  \
        for (_j = 0; _j < NDIM; _j++)                                   \
            (v)[_i] += (p)[_i][_j] * (u)[_j];                           \
    }                                                                   \
}

#define OUTVP(p,v,u)            /* OUTer Vector Product */              \
{                                                                       \
    int _i, _j;                                                         \
    for (_i = 0; _i < NDIM; _i++)                                       \
        for (_j = 0; _j < NDIM; _j++)                                   \
            (p)[_i][_j] = (v)[_i] * (u)[_j];                            \
}

#define TRACEM(s,p)             /* TRACE of Matrix */                   \
{                                                                       \
    int _i;                                                             \
    (s) = 0.0;                                                          \
    for (_i = 0.0; _i < NDIM; _i++)                                     \
        (s) += (p)[_i][_i];                                             \
}

/*
 * Enhancements for tree codes.
 */

#if defined(THREEDIM)

#define DOTPSUBV(s,v,u,w)       /* SUB Vectors, form DOT Prod */        \
{                                                                       \
    (v)[0] = (u)[0] - (w)[0];    (s)  = (v)[0] * (v)[0];                \
    (v)[1] = (u)[1] - (w)[1];    (s) += (v)[1] * (v)[1];                \
    (v)[2] = (u)[2] - (w)[2];    (s) += (v)[2] * (v)[2];                \
}

#define DOTPMULMV(s,v,p,u)      /* MUL Mat by Vect, form DOT Prod */    \
{                                                                       \
    DOTVP(v[0], p[0], u);    (s)  = (v)[0] * (u)[0];                    \
    DOTVP(v[1], p[1], u);    (s) += (v)[1] * (u)[1];                    \
    DOTVP(v[2], p[2], u);    (s) += (v)[2] * (u)[2];                    \
}

#define ADDMULVS(v,u,s)         /* MUL Vect by Scalar, ADD to vect */   \
{                                                                       \
    (v)[0] += (u)[0] * (s);                                             \
    (v)[1] += (u)[1] * (s);                                             \
    (v)[2] += (u)[2] * (s);                                             \
}

#define ADDMULVS2(v,u,s,w,r)    /* 2 times MUL V by S, ADD to vect */   \
{                                                                       \
    (v)[0] += (u)[0] * (s) + (w)[0] * (r);                              \
    (v)[1] += (u)[1] * (s) + (w)[1] * (r);                              \
    (v)[2] += (u)[2] * (s) + (w)[2] * (r);                              \
}

#endif

/*
 * Misc. impure operations.
 */

#define SETVS(v,s)              /* SET Vector to Scalar */              \
{                                                                       \
    int _i;                                                             \
    for (_i = 0; _i < NDIM; _i++)                                       \
        (v)[_i] = (s);                                                  \
}

#define ADDVS(v,u,s)            /* ADD Vector and Scalar */             \
{                                                                       \
    int _i;                                                             \
    for (_i = 0; _i < NDIM; _i++)                                       \
        (v)[_i] = (u)[_i] + (s);                                        \
}

#define SETMS(p,s)              /* SET Matrix to Scalar */              \
{                                                                       \
    int _i, _j;                                                         \
    for (_i = 0; _i < NDIM; _i++)                                       \
        for (_j = 0; _j < NDIM; _j++)                                   \
            (p)[_i][_j] = (s);                                          \
}

#endif  /* ! _vectmath_h */