/****************************************************************************/ /* 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 */