Index: trunk/psLib/src/math/psMathUtils.c =================================================================== --- trunk/psLib/src/math/psMathUtils.c (revision 6305) +++ trunk/psLib/src/math/psMathUtils.c (revision 7132) @@ -1,8 +1,8 @@ /** @file psMathUtils.c * - * This file contains standard math rotines. + * This file contains standard math routines. * - * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-02-02 21:09:07 $ + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-05-18 01:20:51 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -48,31 +48,25 @@ /***************************************************************************** -vectorBinDisectTYPE(): This is a macro for a private function which takes as -input an array of data as well as a single value for that data. The input -vector values are assumed to be non-decreasing (v[i-1] <= v[i] for all i). -This routine does a binary disection of the vector and returns "i" such that -(v[i] <= x <= v[i+1). If x lies outside the range of v[], then this routine -prints a warning message and returns (-2 or -1). +This is a macro covering the various types for the below function. *****************************************************************************/ -#define FUNC_MACRO_VECTOR_BIN_DISECT(TYPE) \ -static psS32 vectorBinDisect##TYPE(ps##TYPE *bins, \ - psS32 numBins, \ - ps##TYPE x) \ -{ \ - psS32 min; \ - psS32 max; \ - psS32 mid; \ +#define VECTOR_BIN_DISECT_CASE(TYPE) \ +case PS_TYPE_##TYPE: { \ + ps##TYPE *bounds = bins->data.TYPE; \ + ps##TYPE value = x->data.TYPE; \ + long min; \ + long max; \ + long mid; \ psTrace(__func__, 4, "---- %s() begin ----\n", __func__); \ /* psTrace(__func__, 6, "Determining the bin for: %f\n", x); */\ - if (x < bins[0]) { \ + if (value < bounds[0]) { \ psLogMsg(__func__, PS_LOG_WARN, \ "vectorBinDisect%s(): ordinate %f is outside vector range (%f - %f).", \ - #TYPE, x, bins[0], bins[numBins-1]); \ + #TYPE, value, bounds[0], bounds[numBins-1]); \ return(-2); \ } \ - if (x > bins[numBins-1]) { \ + if (value > bounds[numBins-1]) { \ psLogMsg(__func__, PS_LOG_WARN, \ "vectorBinDisect%s(): ordinate %f is outside vector range (%f - %f).", \ - #TYPE, x, bins[0], bins[numBins-1]); \ + #TYPE, value, bounds[0], bounds[numBins-1]); \ return(-1); \ } \ @@ -82,11 +76,9 @@ mid = ((max+1)-min)/2; \ while (min != max) { \ - /*psTrace(__func__, 6, "(min, mid, max) is (%u, %u, %u): (x, bins[mid]) is (%f, %f)\n", min, mid, max, x, bins[mid]); */\ \ - if (x == bins[mid]) { \ - /*psTrace(__func__, 6, "%.2f should be in this range (%.2f - %.2f)\n", x, bins[mid], bins[mid+1]); */\ + if (value == bounds[mid]) { \ psTrace(__func__, 4, "---- %s(%d) end (1) ----\n", __func__, mid); \ return(mid); \ - } else if (x < bins[mid]) { \ + } else if (value < bounds[mid]) { \ max = mid-1; \ } else { \ @@ -95,26 +87,17 @@ mid = ((max+1)+min)/2; \ } \ - /*psTrace(__func__, 6, "%.2f should be in this range (%.2f - %.2f)\n", x, bins[min], bins[min+1]);*/ \ psTrace(__func__, 4, "---- %s(%d) end (2) ----\n", __func__, min); \ return(min); \ -} \ - -FUNC_MACRO_VECTOR_BIN_DISECT(S8) -FUNC_MACRO_VECTOR_BIN_DISECT(S16) -FUNC_MACRO_VECTOR_BIN_DISECT(S32) -FUNC_MACRO_VECTOR_BIN_DISECT(S64) -FUNC_MACRO_VECTOR_BIN_DISECT(U8) -FUNC_MACRO_VECTOR_BIN_DISECT(U16) -FUNC_MACRO_VECTOR_BIN_DISECT(U32) -FUNC_MACRO_VECTOR_BIN_DISECT(U64) -FUNC_MACRO_VECTOR_BIN_DISECT(F32) -FUNC_MACRO_VECTOR_BIN_DISECT(F64) +} /***************************************************************************** -p_psVectorBinDisect(): A wrapper to the above p_psVectorBinDisect(). +p_psVectorBinDisect(): This function takes as input an array of data as well as a single value for that data. +The input vector values are assumed to be non-decreasing (v[i-1] <= v[i] for all i). This routine does a +binary disection of the vector and returns "i" such that (v[i] <= x <= v[i+1). If x lies outside the range of +v[], then this routine prints a warning message and returns (-2 or -1). *****************************************************************************/ psS32 p_psVectorBinDisect( - psVector *bins, - psScalar *x) + const psVector *bins, + const psScalar *x) { PS_ASSERT_VECTOR_NON_NULL(bins, -4); @@ -122,31 +105,23 @@ PS_ASSERT_PTR_NON_NULL(x, -6); PS_ASSERT_PTR_TYPE_EQUAL(x, bins, -3); - char* strType; + long numBins = bins->n; // Number of bins switch (x->type.type) { - case PS_TYPE_U8: - return(vectorBinDisectU8(bins->data.U8, bins->n, x->data.U8)); - case PS_TYPE_U16: - return(vectorBinDisectU16(bins->data.U16, bins->n, x->data.U16)); - case PS_TYPE_U32: - return(vectorBinDisectU32(bins->data.U32, bins->n, x->data.U32)); - case PS_TYPE_U64: - return(vectorBinDisectU64(bins->data.U64, bins->n, x->data.U64)); - case PS_TYPE_S8: - return(vectorBinDisectS8(bins->data.S8, bins->n, x->data.S8)); - case PS_TYPE_S16: - return(vectorBinDisectS16(bins->data.S16, bins->n, x->data.S16)); - case PS_TYPE_S32: - return(vectorBinDisectS32(bins->data.S32, bins->n, x->data.S32)); - case PS_TYPE_S64: - return(vectorBinDisectS64(bins->data.S64, bins->n, x->data.S64)); - case PS_TYPE_F32: - return(vectorBinDisectF32(bins->data.F32, bins->n, x->data.F32)); - case PS_TYPE_F64: - return(vectorBinDisectF64(bins->data.F64, bins->n, x->data.F64)); - default: - PS_TYPE_NAME(strType, x->type.type); - psError(PS_ERR_BAD_PARAMETER_TYPE, "psLib type %s is not supported.", strType); - return -1; + VECTOR_BIN_DISECT_CASE(S8); + VECTOR_BIN_DISECT_CASE(S16); + VECTOR_BIN_DISECT_CASE(S32); + VECTOR_BIN_DISECT_CASE(S64); + VECTOR_BIN_DISECT_CASE(U8); + VECTOR_BIN_DISECT_CASE(U16); + VECTOR_BIN_DISECT_CASE(U32); + VECTOR_BIN_DISECT_CASE(U64); + VECTOR_BIN_DISECT_CASE(F32); + VECTOR_BIN_DISECT_CASE(F64); + default: { + char* strType; + PS_TYPE_NAME(strType, x->type.type); + psError(PS_ERR_BAD_PARAMETER_TYPE, "psLib type %s is not supported.", strType); + return -1; + } } return(-3); @@ -154,23 +129,9 @@ - -/***************************************************************************** -fullInterpolateTYPE(): This routine will take as input n-element psVectors -domain and range, and the x value, assumed to lie with the domain vector. It -produces as output the (order-1)-order LaGrange interpolated value of x. - -XXX: do we error check for non-distinct domain values? - -This routine will only be called inside this file. So, we do not have to -do PS_ASSERTS. - *****************************************************************************/ -#define FUNC_MACRO_FULL_INTERPOLATE(TYPE) \ -static psScalar *vectorInterpolate##TYPE( \ - psScalar *out, \ - psVector *domain, \ - psVector *range, \ - psS32 order, \ - psScalar *x) \ -{ \ +/************************************************************************************************************* +Helper macro for p_psVectorInterpolate: handles each of the types +*************************************************************************************************************/ +#define VECTOR_INTERPOLATE_CASE(TYPE) \ +case PS_TYPE_##TYPE: { \ psTrace(__func__, 4, "---- %s() begin %u-order.) (%d data points) ----\n", __func__, order, order+1); \ if (x->data.TYPE < domain->data.TYPE[0]) { \ @@ -213,16 +174,5 @@ psTrace(__func__, 4, "---- %s(....) end ----\n", __func__); \ return(out); \ -} \ - -FUNC_MACRO_FULL_INTERPOLATE(U8) -FUNC_MACRO_FULL_INTERPOLATE(U16) -FUNC_MACRO_FULL_INTERPOLATE(U32) -FUNC_MACRO_FULL_INTERPOLATE(U64) -FUNC_MACRO_FULL_INTERPOLATE(S8) -FUNC_MACRO_FULL_INTERPOLATE(S16) -FUNC_MACRO_FULL_INTERPOLATE(S32) -FUNC_MACRO_FULL_INTERPOLATE(S64) -FUNC_MACRO_FULL_INTERPOLATE(F32) -FUNC_MACRO_FULL_INTERPOLATE(F64) +} /***************************************************************************** @@ -236,8 +186,8 @@ psScalar *p_psVectorInterpolate( psScalar *out, - psVector *domain, - psVector *range, + const psVector *domain, + const psVector *range, psS32 order, - psScalar *x) + const psScalar *x) { psTrace(__func__, 4, "---- %s() begin ----\n", __func__); @@ -249,36 +199,27 @@ PS_ASSERT_PTR_TYPE_EQUAL(domain, range, NULL); PS_ASSERT_PTR_TYPE_EQUAL(domain, x, NULL); - if (out == NULL) { + if (!out) { out = psScalarAlloc(0, x->type.type); } else { PS_ASSERT_PTR_TYPE_EQUAL(domain, out, NULL); } - char* strType; switch (x->type.type) { - case PS_TYPE_U8: - return(vectorInterpolateU8(out, domain, range, order, x)); - case PS_TYPE_U16: - return(vectorInterpolateU16(out, domain, range, order, x)); - case PS_TYPE_U32: - return(vectorInterpolateU32(out, domain, range, order, x)); - case PS_TYPE_U64: - return(vectorInterpolateU64(out, domain, range, order, x)); - case PS_TYPE_S8: - return(vectorInterpolateS8(out, domain, range, order, x)); - case PS_TYPE_S16: - return(vectorInterpolateS16(out, domain, range, order, x)); - case PS_TYPE_S32: - return(vectorInterpolateS32(out, domain, range, order, x)); - case PS_TYPE_S64: - return(vectorInterpolateS64(out, domain, range, order, x)); - case PS_TYPE_F32: - return(vectorInterpolateF32(out, domain, range, order, x)); - case PS_TYPE_F64: - return(vectorInterpolateF64(out, domain, range, order, x)); - default: - PS_TYPE_NAME(strType, x->type.type); - psError(PS_ERR_BAD_PARAMETER_TYPE, "psLib type %s is not supported.", strType); - return NULL; + VECTOR_INTERPOLATE_CASE(U8); + VECTOR_INTERPOLATE_CASE(U16); + VECTOR_INTERPOLATE_CASE(U32); + VECTOR_INTERPOLATE_CASE(U64); + VECTOR_INTERPOLATE_CASE(S8); + VECTOR_INTERPOLATE_CASE(S16); + VECTOR_INTERPOLATE_CASE(S32); + VECTOR_INTERPOLATE_CASE(S64); + VECTOR_INTERPOLATE_CASE(F32); + VECTOR_INTERPOLATE_CASE(F64); + default: { + char* strType; + PS_TYPE_NAME(strType, x->type.type); + psError(PS_ERR_BAD_PARAMETER_TYPE, "psLib type %s is not supported.", strType); + return NULL; + } } @@ -287,34 +228,14 @@ /***************************************************************************** -These macros and functions define the following functions: - - p_psNormalizeVectorRange(myData, low, high) - -That assumes that the low/high arguments are PS_TYPE_F64; the vector myData -can be of any type. Arguments low/high will be converted to the appropriate -type and one of the type-specific functions below will be called: - - p_psNormalizeVectorRangeU8(myData, low, high) - p_psNormalizeVectorRangeU16(myData, low, high) - p_psNormalizeVectorRangeU32(myData, low, high) - p_psNormalizeVectorRangeU64(myData, low, high) - p_psNormalizeVectorRangeS8(myData, low, high) - p_psNormalizeVectorRangeS16(myData, low, high) - p_psNormalizeVectorRangeS32(myData, low, high) - p_psNormalizeVectorRangeS64(myData, low, high) - p_psNormalizeVectorRangeF32(myData, low, high) - p_psNormalizeVectorRangeF64(myData, low, high) +Helper macro for p_psNormalizeVectorRange: handles each of the types *****************************************************************************/ -#define PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(TYPE) \ -static void p_psNormalizeVectorRange##TYPE( \ - psVector* myData, \ - ps##TYPE outLow, \ - ps##TYPE outHigh) \ -{ \ - ps##TYPE min = (ps##TYPE) PS_MAX_##TYPE; \ - ps##TYPE max = (ps##TYPE) -PS_MAX_##TYPE; \ - psS32 i = 0; \ - \ - for (i = 0; i < myData->n; i++) { \ +#define NORMALIZE_VECTOR_RANGE_CASE(TYPE) \ +case PS_TYPE_##TYPE: { \ + ps##TYPE low = outLow; \ + ps##TYPE high = outHigh; \ + ps##TYPE min = (ps##TYPE)PS_MAX_##TYPE; \ + ps##TYPE max = (ps##TYPE)-PS_MAX_##TYPE; \ + \ + for (long i = 0; i < myData->n; i++) { \ if (myData->data.TYPE[i] < min) { \ min = myData->data.TYPE[i]; \ @@ -327,79 +248,53 @@ /* Ensure that max!=min before we divide by (max-min) */ \ if (max != min) { \ - for (i = 0; i < myData->n; i++) { \ - myData->data.TYPE[i] = (outLow + (myData->data.TYPE[i] - min) * \ - (outHigh - outLow) / (max - min)); \ + for (long i = 0; i < myData->n; i++) { \ + myData->data.TYPE[i] = (low + (myData->data.TYPE[i] - min) * \ + (high - low) / (max - min)); \ } \ } else { \ psLogMsg(__func__, PS_LOG_WARN, "WARNING: (max==min). Setting all elements to min.\n"); \ - for (i = 0; i < myData->n; i++) \ + for (long i = 0; i < myData->n; i++) \ { \ \ - myData->data.TYPE[i] = outLow; \ + myData->data.TYPE[i] = low; \ \ } \ } \ + break; \ } \ -PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U8) -PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U16) -PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U32) -PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U64) -PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S8) -PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S16) -PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S32) -PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S64) -PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(F32) -PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(F64) - -psS32 p_psNormalizeVectorRange(psVector* myData, - psF64 outLow, - psF64 outHigh) +/************************************************************************************************************* +p_psNormalizeVectorRange(myData, low, high): this function normalises the vector (myData) to the range +low: +high. +*************************************************************************************************************/ +bool p_psNormalizeVectorRange(psVector* myData, + psF64 outLow, + psF64 outHigh) { - PS_ASSERT_VECTOR_NON_NULL(myData, -1); - char* strType; - + PS_ASSERT_VECTOR_NON_NULL(myData, false); psTrace(__func__, 4, "---- %s() begin ----\n", __func__); + switch (myData->type.type) { - case PS_TYPE_U8: - p_psNormalizeVectorRangeU8(myData, (psU8) outLow, (psU8) outHigh); - break; - case PS_TYPE_U16: - p_psNormalizeVectorRangeU16(myData, (psU16) outLow, (psU16) outHigh); - break; - case PS_TYPE_U32: - p_psNormalizeVectorRangeU32(myData, (psU32) outLow, (psU32) outHigh); - break; - case PS_TYPE_U64: - p_psNormalizeVectorRangeU64(myData, (psU64) outLow, (psU64) outHigh); - break; - case PS_TYPE_S8: - p_psNormalizeVectorRangeS8(myData, (psS8) outLow, (psS8) outHigh); - break; - case PS_TYPE_S16: - p_psNormalizeVectorRangeS16(myData, (psS16) outLow, (psS16) outHigh); - break; - case PS_TYPE_S32: - p_psNormalizeVectorRangeS32(myData, (psS32) outLow, (psS32) outHigh); - break; - case PS_TYPE_S64: - p_psNormalizeVectorRangeS64(myData, (psS64) outLow, (psS64) outHigh); - break; - case PS_TYPE_F32: - p_psNormalizeVectorRangeF32(myData, (psF32) outLow, (psF32) outHigh); - break; - case PS_TYPE_F64: - p_psNormalizeVectorRangeF64(myData, (psF64) outLow, (psF64) outHigh); - break; - case PS_TYPE_C32: - case PS_TYPE_C64: - default: - PS_TYPE_NAME(strType, myData->type.type); - psError(PS_ERR_BAD_PARAMETER_TYPE, "psLib type %s is not supported.", strType); - break; + NORMALIZE_VECTOR_RANGE_CASE(U8); + NORMALIZE_VECTOR_RANGE_CASE(U16); + NORMALIZE_VECTOR_RANGE_CASE(U32); + NORMALIZE_VECTOR_RANGE_CASE(U64); + NORMALIZE_VECTOR_RANGE_CASE(S8); + NORMALIZE_VECTOR_RANGE_CASE(S16); + NORMALIZE_VECTOR_RANGE_CASE(S32); + NORMALIZE_VECTOR_RANGE_CASE(S64); + NORMALIZE_VECTOR_RANGE_CASE(F32); + NORMALIZE_VECTOR_RANGE_CASE(F64); + default: { + char* strType; + PS_TYPE_NAME(strType, myData->type.type); + psError(PS_ERR_BAD_PARAMETER_TYPE, "psLib type %s is not supported.", strType); + break; + } } psTrace(__func__, 4, "---- %s() end ----\n", __func__); - return(0); -} - - + return true; +} + +