Changeset 7132 for trunk/psLib/src/math/psMathUtils.c
- Timestamp:
- May 17, 2006, 3:20:51 PM (20 years ago)
- File:
-
- 1 edited
-
trunk/psLib/src/math/psMathUtils.c (modified) (11 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psLib/src/math/psMathUtils.c
r6305 r7132 1 1 /** @file psMathUtils.c 2 2 * 3 * This file contains standard math ro tines.3 * This file contains standard math routines. 4 4 * 5 * @version $Revision: 1. 1$ $Name: not supported by cvs2svn $6 * @date $Date: 2006-0 2-02 21:09:07$5 * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ 6 * @date $Date: 2006-05-18 01:20:51 $ 7 7 * 8 8 * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii … … 48 48 49 49 /***************************************************************************** 50 vectorBinDisectTYPE(): This is a macro for a private function which takes as 51 input an array of data as well as a single value for that data. The input 52 vector values are assumed to be non-decreasing (v[i-1] <= v[i] for all i). 53 This routine does a binary disection of the vector and returns "i" such that 54 (v[i] <= x <= v[i+1). If x lies outside the range of v[], then this routine 55 prints a warning message and returns (-2 or -1). 50 This is a macro covering the various types for the below function. 56 51 *****************************************************************************/ 57 #define FUNC_MACRO_VECTOR_BIN_DISECT(TYPE) \ 58 static psS32 vectorBinDisect##TYPE(ps##TYPE *bins, \ 59 psS32 numBins, \ 60 ps##TYPE x) \ 61 { \ 62 psS32 min; \ 63 psS32 max; \ 64 psS32 mid; \ 52 #define VECTOR_BIN_DISECT_CASE(TYPE) \ 53 case PS_TYPE_##TYPE: { \ 54 ps##TYPE *bounds = bins->data.TYPE; \ 55 ps##TYPE value = x->data.TYPE; \ 56 long min; \ 57 long max; \ 58 long mid; \ 65 59 psTrace(__func__, 4, "---- %s() begin ----\n", __func__); \ 66 60 /* psTrace(__func__, 6, "Determining the bin for: %f\n", x); */\ 67 if ( x < bins[0]) { \61 if (value < bounds[0]) { \ 68 62 psLogMsg(__func__, PS_LOG_WARN, \ 69 63 "vectorBinDisect%s(): ordinate %f is outside vector range (%f - %f).", \ 70 #TYPE, x, bins[0], bins[numBins-1]); \64 #TYPE, value, bounds[0], bounds[numBins-1]); \ 71 65 return(-2); \ 72 66 } \ 73 if ( x > bins[numBins-1]) { \67 if (value > bounds[numBins-1]) { \ 74 68 psLogMsg(__func__, PS_LOG_WARN, \ 75 69 "vectorBinDisect%s(): ordinate %f is outside vector range (%f - %f).", \ 76 #TYPE, x, bins[0], bins[numBins-1]); \70 #TYPE, value, bounds[0], bounds[numBins-1]); \ 77 71 return(-1); \ 78 72 } \ … … 82 76 mid = ((max+1)-min)/2; \ 83 77 while (min != max) { \ 84 /*psTrace(__func__, 6, "(min, mid, max) is (%u, %u, %u): (x, bins[mid]) is (%f, %f)\n", min, mid, max, x, bins[mid]); */\85 78 \ 86 if (x == bins[mid]) { \ 87 /*psTrace(__func__, 6, "%.2f should be in this range (%.2f - %.2f)\n", x, bins[mid], bins[mid+1]); */\ 79 if (value == bounds[mid]) { \ 88 80 psTrace(__func__, 4, "---- %s(%d) end (1) ----\n", __func__, mid); \ 89 81 return(mid); \ 90 } else if ( x < bins[mid]) { \82 } else if (value < bounds[mid]) { \ 91 83 max = mid-1; \ 92 84 } else { \ … … 95 87 mid = ((max+1)+min)/2; \ 96 88 } \ 97 /*psTrace(__func__, 6, "%.2f should be in this range (%.2f - %.2f)\n", x, bins[min], bins[min+1]);*/ \98 89 psTrace(__func__, 4, "---- %s(%d) end (2) ----\n", __func__, min); \ 99 90 return(min); \ 100 } \ 101 102 FUNC_MACRO_VECTOR_BIN_DISECT(S8) 103 FUNC_MACRO_VECTOR_BIN_DISECT(S16) 104 FUNC_MACRO_VECTOR_BIN_DISECT(S32) 105 FUNC_MACRO_VECTOR_BIN_DISECT(S64) 106 FUNC_MACRO_VECTOR_BIN_DISECT(U8) 107 FUNC_MACRO_VECTOR_BIN_DISECT(U16) 108 FUNC_MACRO_VECTOR_BIN_DISECT(U32) 109 FUNC_MACRO_VECTOR_BIN_DISECT(U64) 110 FUNC_MACRO_VECTOR_BIN_DISECT(F32) 111 FUNC_MACRO_VECTOR_BIN_DISECT(F64) 91 } 112 92 113 93 /***************************************************************************** 114 p_psVectorBinDisect(): A wrapper to the above p_psVectorBinDisect(). 94 p_psVectorBinDisect(): This function takes as input an array of data as well as a single value for that data. 95 The input vector values are assumed to be non-decreasing (v[i-1] <= v[i] for all i). This routine does a 96 binary disection of the vector and returns "i" such that (v[i] <= x <= v[i+1). If x lies outside the range of 97 v[], then this routine prints a warning message and returns (-2 or -1). 115 98 *****************************************************************************/ 116 99 psS32 p_psVectorBinDisect( 117 psVector *bins,118 psScalar *x)100 const psVector *bins, 101 const psScalar *x) 119 102 { 120 103 PS_ASSERT_VECTOR_NON_NULL(bins, -4); … … 122 105 PS_ASSERT_PTR_NON_NULL(x, -6); 123 106 PS_ASSERT_PTR_TYPE_EQUAL(x, bins, -3); 124 char* strType;107 long numBins = bins->n; // Number of bins 125 108 126 109 switch (x->type.type) { 127 case PS_TYPE_U8: 128 return(vectorBinDisectU8(bins->data.U8, bins->n, x->data.U8)); 129 case PS_TYPE_U16: 130 return(vectorBinDisectU16(bins->data.U16, bins->n, x->data.U16)); 131 case PS_TYPE_U32: 132 return(vectorBinDisectU32(bins->data.U32, bins->n, x->data.U32)); 133 case PS_TYPE_U64: 134 return(vectorBinDisectU64(bins->data.U64, bins->n, x->data.U64)); 135 case PS_TYPE_S8: 136 return(vectorBinDisectS8(bins->data.S8, bins->n, x->data.S8)); 137 case PS_TYPE_S16: 138 return(vectorBinDisectS16(bins->data.S16, bins->n, x->data.S16)); 139 case PS_TYPE_S32: 140 return(vectorBinDisectS32(bins->data.S32, bins->n, x->data.S32)); 141 case PS_TYPE_S64: 142 return(vectorBinDisectS64(bins->data.S64, bins->n, x->data.S64)); 143 case PS_TYPE_F32: 144 return(vectorBinDisectF32(bins->data.F32, bins->n, x->data.F32)); 145 case PS_TYPE_F64: 146 return(vectorBinDisectF64(bins->data.F64, bins->n, x->data.F64)); 147 default: 148 PS_TYPE_NAME(strType, x->type.type); 149 psError(PS_ERR_BAD_PARAMETER_TYPE, "psLib type %s is not supported.", strType); 150 return -1; 110 VECTOR_BIN_DISECT_CASE(S8); 111 VECTOR_BIN_DISECT_CASE(S16); 112 VECTOR_BIN_DISECT_CASE(S32); 113 VECTOR_BIN_DISECT_CASE(S64); 114 VECTOR_BIN_DISECT_CASE(U8); 115 VECTOR_BIN_DISECT_CASE(U16); 116 VECTOR_BIN_DISECT_CASE(U32); 117 VECTOR_BIN_DISECT_CASE(U64); 118 VECTOR_BIN_DISECT_CASE(F32); 119 VECTOR_BIN_DISECT_CASE(F64); 120 default: { 121 char* strType; 122 PS_TYPE_NAME(strType, x->type.type); 123 psError(PS_ERR_BAD_PARAMETER_TYPE, "psLib type %s is not supported.", strType); 124 return -1; 125 } 151 126 } 152 127 return(-3); … … 154 129 155 130 156 157 /***************************************************************************** 158 fullInterpolateTYPE(): This routine will take as input n-element psVectors 159 domain and range, and the x value, assumed to lie with the domain vector. It 160 produces as output the (order-1)-order LaGrange interpolated value of x. 161 162 XXX: do we error check for non-distinct domain values? 163 164 This routine will only be called inside this file. So, we do not have to 165 do PS_ASSERTS. 166 *****************************************************************************/ 167 #define FUNC_MACRO_FULL_INTERPOLATE(TYPE) \ 168 static psScalar *vectorInterpolate##TYPE( \ 169 psScalar *out, \ 170 psVector *domain, \ 171 psVector *range, \ 172 psS32 order, \ 173 psScalar *x) \ 174 { \ 131 /************************************************************************************************************* 132 Helper macro for p_psVectorInterpolate: handles each of the types 133 *************************************************************************************************************/ 134 #define VECTOR_INTERPOLATE_CASE(TYPE) \ 135 case PS_TYPE_##TYPE: { \ 175 136 psTrace(__func__, 4, "---- %s() begin %u-order.) (%d data points) ----\n", __func__, order, order+1); \ 176 137 if (x->data.TYPE < domain->data.TYPE[0]) { \ … … 213 174 psTrace(__func__, 4, "---- %s(....) end ----\n", __func__); \ 214 175 return(out); \ 215 } \ 216 217 FUNC_MACRO_FULL_INTERPOLATE(U8) 218 FUNC_MACRO_FULL_INTERPOLATE(U16) 219 FUNC_MACRO_FULL_INTERPOLATE(U32) 220 FUNC_MACRO_FULL_INTERPOLATE(U64) 221 FUNC_MACRO_FULL_INTERPOLATE(S8) 222 FUNC_MACRO_FULL_INTERPOLATE(S16) 223 FUNC_MACRO_FULL_INTERPOLATE(S32) 224 FUNC_MACRO_FULL_INTERPOLATE(S64) 225 FUNC_MACRO_FULL_INTERPOLATE(F32) 226 FUNC_MACRO_FULL_INTERPOLATE(F64) 176 } 227 177 228 178 /***************************************************************************** … … 236 186 psScalar *p_psVectorInterpolate( 237 187 psScalar *out, 238 psVector *domain,239 psVector *range,188 const psVector *domain, 189 const psVector *range, 240 190 psS32 order, 241 psScalar *x)191 const psScalar *x) 242 192 { 243 193 psTrace(__func__, 4, "---- %s() begin ----\n", __func__); … … 249 199 PS_ASSERT_PTR_TYPE_EQUAL(domain, range, NULL); 250 200 PS_ASSERT_PTR_TYPE_EQUAL(domain, x, NULL); 251 if ( out == NULL) {201 if (!out) { 252 202 out = psScalarAlloc(0, x->type.type); 253 203 } else { 254 204 PS_ASSERT_PTR_TYPE_EQUAL(domain, out, NULL); 255 205 } 256 char* strType;257 206 258 207 switch (x->type.type) { 259 case PS_TYPE_U8: 260 return(vectorInterpolateU8(out, domain, range, order, x)); 261 case PS_TYPE_U16: 262 return(vectorInterpolateU16(out, domain, range, order, x)); 263 case PS_TYPE_U32: 264 return(vectorInterpolateU32(out, domain, range, order, x)); 265 case PS_TYPE_U64: 266 return(vectorInterpolateU64(out, domain, range, order, x)); 267 case PS_TYPE_S8: 268 return(vectorInterpolateS8(out, domain, range, order, x)); 269 case PS_TYPE_S16: 270 return(vectorInterpolateS16(out, domain, range, order, x)); 271 case PS_TYPE_S32: 272 return(vectorInterpolateS32(out, domain, range, order, x)); 273 case PS_TYPE_S64: 274 return(vectorInterpolateS64(out, domain, range, order, x)); 275 case PS_TYPE_F32: 276 return(vectorInterpolateF32(out, domain, range, order, x)); 277 case PS_TYPE_F64: 278 return(vectorInterpolateF64(out, domain, range, order, x)); 279 default: 280 PS_TYPE_NAME(strType, x->type.type); 281 psError(PS_ERR_BAD_PARAMETER_TYPE, "psLib type %s is not supported.", strType); 282 return NULL; 208 VECTOR_INTERPOLATE_CASE(U8); 209 VECTOR_INTERPOLATE_CASE(U16); 210 VECTOR_INTERPOLATE_CASE(U32); 211 VECTOR_INTERPOLATE_CASE(U64); 212 VECTOR_INTERPOLATE_CASE(S8); 213 VECTOR_INTERPOLATE_CASE(S16); 214 VECTOR_INTERPOLATE_CASE(S32); 215 VECTOR_INTERPOLATE_CASE(S64); 216 VECTOR_INTERPOLATE_CASE(F32); 217 VECTOR_INTERPOLATE_CASE(F64); 218 default: { 219 char* strType; 220 PS_TYPE_NAME(strType, x->type.type); 221 psError(PS_ERR_BAD_PARAMETER_TYPE, "psLib type %s is not supported.", strType); 222 return NULL; 223 } 283 224 } 284 225 … … 287 228 288 229 /***************************************************************************** 289 These macros and functions define the following functions: 290 291 p_psNormalizeVectorRange(myData, low, high) 292 293 That assumes that the low/high arguments are PS_TYPE_F64; the vector myData 294 can be of any type. Arguments low/high will be converted to the appropriate 295 type and one of the type-specific functions below will be called: 296 297 p_psNormalizeVectorRangeU8(myData, low, high) 298 p_psNormalizeVectorRangeU16(myData, low, high) 299 p_psNormalizeVectorRangeU32(myData, low, high) 300 p_psNormalizeVectorRangeU64(myData, low, high) 301 p_psNormalizeVectorRangeS8(myData, low, high) 302 p_psNormalizeVectorRangeS16(myData, low, high) 303 p_psNormalizeVectorRangeS32(myData, low, high) 304 p_psNormalizeVectorRangeS64(myData, low, high) 305 p_psNormalizeVectorRangeF32(myData, low, high) 306 p_psNormalizeVectorRangeF64(myData, low, high) 230 Helper macro for p_psNormalizeVectorRange: handles each of the types 307 231 *****************************************************************************/ 308 #define PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(TYPE) \ 309 static void p_psNormalizeVectorRange##TYPE( \ 310 psVector* myData, \ 311 ps##TYPE outLow, \ 312 ps##TYPE outHigh) \ 313 { \ 314 ps##TYPE min = (ps##TYPE) PS_MAX_##TYPE; \ 315 ps##TYPE max = (ps##TYPE) -PS_MAX_##TYPE; \ 316 psS32 i = 0; \ 317 \ 318 for (i = 0; i < myData->n; i++) { \ 232 #define NORMALIZE_VECTOR_RANGE_CASE(TYPE) \ 233 case PS_TYPE_##TYPE: { \ 234 ps##TYPE low = outLow; \ 235 ps##TYPE high = outHigh; \ 236 ps##TYPE min = (ps##TYPE)PS_MAX_##TYPE; \ 237 ps##TYPE max = (ps##TYPE)-PS_MAX_##TYPE; \ 238 \ 239 for (long i = 0; i < myData->n; i++) { \ 319 240 if (myData->data.TYPE[i] < min) { \ 320 241 min = myData->data.TYPE[i]; \ … … 327 248 /* Ensure that max!=min before we divide by (max-min) */ \ 328 249 if (max != min) { \ 329 for ( i = 0; i < myData->n; i++) { \330 myData->data.TYPE[i] = ( outLow + (myData->data.TYPE[i] - min) * \331 ( outHigh - outLow) / (max - min)); \250 for (long i = 0; i < myData->n; i++) { \ 251 myData->data.TYPE[i] = (low + (myData->data.TYPE[i] - min) * \ 252 (high - low) / (max - min)); \ 332 253 } \ 333 254 } else { \ 334 255 psLogMsg(__func__, PS_LOG_WARN, "WARNING: (max==min). Setting all elements to min.\n"); \ 335 for ( i = 0; i < myData->n; i++) \256 for (long i = 0; i < myData->n; i++) \ 336 257 { \ 337 258 \ 338 myData->data.TYPE[i] = outLow; \259 myData->data.TYPE[i] = low; \ 339 260 \ 340 261 } \ 341 262 } \ 263 break; \ 342 264 } \ 343 265 344 PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U8) 345 PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U16) 346 PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U32) 347 PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U64) 348 PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S8) 349 PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S16) 350 PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S32) 351 PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S64) 352 PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(F32) 353 PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(F64) 354 355 psS32 p_psNormalizeVectorRange(psVector* myData, 356 psF64 outLow, 357 psF64 outHigh) 266 /************************************************************************************************************* 267 p_psNormalizeVectorRange(myData, low, high): this function normalises the vector (myData) to the range 268 low: 269 high. 270 *************************************************************************************************************/ 271 bool p_psNormalizeVectorRange(psVector* myData, 272 psF64 outLow, 273 psF64 outHigh) 358 274 { 359 PS_ASSERT_VECTOR_NON_NULL(myData, -1); 360 char* strType; 361 275 PS_ASSERT_VECTOR_NON_NULL(myData, false); 362 276 psTrace(__func__, 4, "---- %s() begin ----\n", __func__); 277 363 278 switch (myData->type.type) { 364 case PS_TYPE_U8: 365 p_psNormalizeVectorRangeU8(myData, (psU8) outLow, (psU8) outHigh); 366 break; 367 case PS_TYPE_U16: 368 p_psNormalizeVectorRangeU16(myData, (psU16) outLow, (psU16) outHigh); 369 break; 370 case PS_TYPE_U32: 371 p_psNormalizeVectorRangeU32(myData, (psU32) outLow, (psU32) outHigh); 372 break; 373 case PS_TYPE_U64: 374 p_psNormalizeVectorRangeU64(myData, (psU64) outLow, (psU64) outHigh); 375 break; 376 case PS_TYPE_S8: 377 p_psNormalizeVectorRangeS8(myData, (psS8) outLow, (psS8) outHigh); 378 break; 379 case PS_TYPE_S16: 380 p_psNormalizeVectorRangeS16(myData, (psS16) outLow, (psS16) outHigh); 381 break; 382 case PS_TYPE_S32: 383 p_psNormalizeVectorRangeS32(myData, (psS32) outLow, (psS32) outHigh); 384 break; 385 case PS_TYPE_S64: 386 p_psNormalizeVectorRangeS64(myData, (psS64) outLow, (psS64) outHigh); 387 break; 388 case PS_TYPE_F32: 389 p_psNormalizeVectorRangeF32(myData, (psF32) outLow, (psF32) outHigh); 390 break; 391 case PS_TYPE_F64: 392 p_psNormalizeVectorRangeF64(myData, (psF64) outLow, (psF64) outHigh); 393 break; 394 case PS_TYPE_C32: 395 case PS_TYPE_C64: 396 default: 397 PS_TYPE_NAME(strType, myData->type.type); 398 psError(PS_ERR_BAD_PARAMETER_TYPE, "psLib type %s is not supported.", strType); 399 break; 279 NORMALIZE_VECTOR_RANGE_CASE(U8); 280 NORMALIZE_VECTOR_RANGE_CASE(U16); 281 NORMALIZE_VECTOR_RANGE_CASE(U32); 282 NORMALIZE_VECTOR_RANGE_CASE(U64); 283 NORMALIZE_VECTOR_RANGE_CASE(S8); 284 NORMALIZE_VECTOR_RANGE_CASE(S16); 285 NORMALIZE_VECTOR_RANGE_CASE(S32); 286 NORMALIZE_VECTOR_RANGE_CASE(S64); 287 NORMALIZE_VECTOR_RANGE_CASE(F32); 288 NORMALIZE_VECTOR_RANGE_CASE(F64); 289 default: { 290 char* strType; 291 PS_TYPE_NAME(strType, myData->type.type); 292 psError(PS_ERR_BAD_PARAMETER_TYPE, "psLib type %s is not supported.", strType); 293 break; 294 } 400 295 } 401 296 psTrace(__func__, 4, "---- %s() end ----\n", __func__); 402 return (0);403 } 404 405 297 return true; 298 } 299 300
Note:
See TracChangeset
for help on using the changeset viewer.
