Index: trunk/psLib/src/dataManip/psMinimize.c =================================================================== --- trunk/psLib/src/dataManip/psMinimize.c (revision 1945) +++ trunk/psLib/src/dataManip/psMinimize.c (revision 1953) @@ -9,6 +9,6 @@ * @author George Gusciora, MHPCC * - * @version $Revision: 1.51 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-10-03 23:35:47 $ + * @version $Revision: 1.52 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-10-05 01:43:58 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -24,30 +24,5 @@ /* INCLUDE FILES */ /*****************************************************************************/ -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include - -#include "psMemory.h" -#include "psVector.h" -#include "psImage.h" -#include "psTrace.h" -#include "psLogMsg.h" -#include "psError.h" -#include "psAbort.h" -#include "psFunctions.h" #include "psMinimize.h" -#include "psMatrix.h" -#include "psConstants.h" - /*****************************************************************************/ /* DEFINE STATEMENTS */ @@ -655,4 +630,10 @@ } +#define GEN_STATIC_RECYCLED_VECTOR(NAME, SIZE, TYPE) \ +static psVector *NAME = NULL; \ +NAME = psVectorRecycle(NAME, SIZE, TYPE); \ +p_psMemSetPersistent(NAME, true); \ +p_psMemSetPersistent(NAME->data.V, true); \ + /****************************************************************************** @@ -661,10 +642,10 @@ coefficients of that polynomial. -XXX: yErr is currently ignored. - -XXX: must add type F32 (currently F64 only). - -XXX: Use private name? - *****************************************************************************/ +XXX: +yErr is currently ignored. + +XXX: +Use private name? +*****************************************************************************/ psPolynomial1D *p_psVectorFitPolynomial1DCheby(psPolynomial1D* myPoly, const psVector* restrict x, @@ -675,10 +656,10 @@ int k; int n = x->n; - psVector *f = psVectorAlloc(n, PS_TYPE_F64); double fac; double sum; - // XXX: Use static memory here. - psScalar *tmpScalar = psScalarAlloc(0.0, PS_TYPE_F32); + GEN_STATIC_RECYCLED_VECTOR(f, n, PS_TYPE_F64); psScalar *fScalar; + psScalar tmpScalar; + tmpScalar.type.type = PS_TYPE_F64; // XXX: These assignments appear too simple to warrant code and @@ -690,9 +671,4 @@ double bpa = 0.5 * (max+min); // 0 - // XXX: Eliminate this later by generating a F64 version of the - // LaGrange interpolation routines. - PS_VECTOR_F64_TO_F32(x, x32); - PS_VECTOR_F64_TO_F32(y, y32); - // In this loop, we first calculate the values of X for which the // Chebyshev polynomials are zero (see NR, section 5.4). Then we @@ -707,12 +683,14 @@ double Y = cos(M_PI * (0.5 + ((float) i)) / ((float) n)); double X = (Y + bma + bpa) - 1.0; - tmpScalar->data.F32 = (float) X; + tmpScalar.data.F64 = X; // We interpolate against are tabluated x,y vectors to determine the // function value at X. - fScalar = p_psVectorInterpolate((psVector *) x32, (psVector *) y32, - 3, tmpScalar); - - f->data.F64[i] = (double) fScalar->data.F32; + fScalar = p_psVectorInterpolate((psVector *) x, + (psVector *) y, + 3, + &tmpScalar); + + f->data.F64[i] = fScalar->data.F64; psFree(fScalar); @@ -724,6 +702,6 @@ // We have the values for f() at the zero points, we now calculate the // coefficients of the Chebyshev polynomial: NR 5.8.7. + fac = 2.0/((float) n); - // XXX: is this loop bound correct? for (j=0;jn;j++) { @@ -737,10 +715,4 @@ } - // XXX: Must free memory. - psFree(f); - psFree(x32); - psFree(y32); - psFree(tmpScalar); - return(myPoly); } @@ -750,6 +722,4 @@ polynomial of degree myPoly to the data points (x, y) and return the coefficients of that polynomial. - -XXX: must add type F32 (currently F64 only). XXX: Use private name? @@ -875,5 +845,5 @@ polynomial. -XXX: must add type F32 (currently F64 only). +XXX: type F32 is done via vector conversion only. *****************************************************************************/ psPolynomial1D* psVectorFitPolynomial1D(psPolynomial1D* myPoly, @@ -882,82 +852,96 @@ const psVector* restrict yErr) { - bool mustFreeMyYErr = false; - bool mustFreemyX = false; int i; psPolynomial1D *tmpPoly; - psVector *myX = NULL; - psVector *myYErr = NULL; - PS_CHECK_NULL_1DPOLY_ACTION(myPoly, 0); PS_CHECK_NULL_VECTOR_ACTION(y, 0); PS_CHECK_EMPTY_VECTOR_ACTION(y, 0); + psVector *x64 = NULL; + psVector *y64 = NULL; + psVector *yErr64 = NULL; + static psVector *x64Static = NULL; + static psVector *y64Static = NULL; + static psVector *yErr64Static = NULL; + + if (y->type.type == PS_TYPE_F64) { + y64 = (psVector *) y; + } else if (y->type.type == PS_TYPE_F32) { + y64Static = psVectorRecycle(y64Static, y->n, PS_TYPE_F64); + p_psMemSetPersistent(y64Static, true); + p_psMemSetPersistent(y64Static->data.V, true); + y64 = y64Static; + } else { + // XXX: psError() bad type. + psAbort(__func__, "Bad type for y64 (%d)", y->type.type); + } + // If yErr==NULL, set all errors equal. if (yErr == NULL) { - myYErr = psVectorAlloc(y->n, y->type.type); - mustFreeMyYErr = true; - - if (y->type.type == PS_TYPE_F32) { - for (i=0;in;i++) { - myYErr->data.F32[i] = 1.0; - } - } else if (y->type.type == PS_TYPE_F64) { - for (i=0;in;i++) { - myYErr->data.F64[i] = 1.0; - } - } + yErr64Static = psVectorRecycle(yErr64Static, y->n, PS_TYPE_F64); + p_psMemSetPersistent(yErr64Static, true); + p_psMemSetPersistent(yErr64Static->data.V, true); + + for (i=0;in;i++) { + yErr64Static->data.F64[i] = 1.0; + } + yErr64 = yErr64Static; } else { - myYErr = (psVector *) yErr; - } - - // If x==NULL, create an myX vector with x values set to (0:n), and if + if (yErr->type.type == PS_TYPE_F64) { + yErr64 = (psVector *) yErr; + } else if (yErr->type.type == PS_TYPE_F32) { + yErr64Static = psVectorRecycle(yErr64Static, yErr->n, PS_TYPE_F64); + p_psMemSetPersistent(yErr64Static, true); + p_psMemSetPersistent(yErr64Static->data.V, true); + } else { + // XXX: psError() bad type. + psAbort(__func__, "Bad type for yErr64"); + } + } + + // If x==NULL, create an x64 vector with x values set to (0:n), and if // this is a Chebyshev polynomial, we must scale to (-1:1). - // XXX: Verify that this is the correct action. if (x == NULL) { - myX = psVectorAlloc(y->n, y->type.type); - mustFreemyX = true; - - if (y->type.type == PS_TYPE_F32) { - if (myPoly->type == PS_POLYNOMIAL_ORD) { - for (i=0;in;i++) { - myX->data.F32[i] = (float) i; - } - } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { - float min = 0.0; - float max = (float) (y->n - 1); - - for (i=0;in;i++) { - myX->data.F32[i] = (((float) i) - 0.5 * (min + max)) / - (0.5 * (max - min)); - } - } - } else if (y->type.type == PS_TYPE_F64) { - if (myPoly->type == PS_POLYNOMIAL_ORD) { - for (i=0;in;i++) { - myX->data.F64[i] = (float) i; - } - } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { - double min = 0.0; - double max = (double) (y->n - 1); - - for (i=0;in;i++) { - myX->data.F64[i] = (((float) i) - 0.5 * (min + max)) / - (0.5 * (max - min)); - } + x64Static = psVectorRecycle(x64Static, y->n, PS_TYPE_F64); + p_psMemSetPersistent(x64Static, true); + p_psMemSetPersistent(x64Static->data.V, true); + x64 = x64Static; + + if (myPoly->type == PS_POLYNOMIAL_ORD) { + for (i=0;in;i++) { + x64->data.F64[i] = (float) i; + } + } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { + double min = 0.0; + double max = (double) (y->n - 1); + + for (i=0;in;i++) { + x64->data.F64[i] = (((double) i) - 0.5 * (min + max)) / + (0.5 * (max - min)); } } } else { - myX = (psVector *) x; - } - - PS_CHECK_VECTOR_SIZE_EQUAL(y, myX); - PS_CHECK_VECTOR_SIZE_EQUAL(y, myYErr); + if (x->type.type == PS_TYPE_F64) { + x64 = (psVector *) x; + } else if (x->type.type == PS_TYPE_F32) { + x64Static = psVectorRecycle(x64Static, x->n, PS_TYPE_F64); + p_psMemSetPersistent(x64Static, true); + p_psMemSetPersistent(x64Static->data.V, true); + x64 = x64Static; + } else { + // XXX: psError() bad type. + psAbort(__func__, "Bad type for x64"); + } + } + + PS_CHECK_VECTOR_SIZE_EQUAL(y64, x64); + PS_CHECK_VECTOR_SIZE_EQUAL(y64, yErr64); // Call the appropriate vector fitting routine. if (myPoly->type == PS_POLYNOMIAL_CHEB) { - tmpPoly = p_psVectorFitPolynomial1DCheby(myPoly, myX, y, myYErr); + tmpPoly = p_psVectorFitPolynomial1DCheby(myPoly, x64, y64, yErr64); } else if (myPoly->type == PS_POLYNOMIAL_ORD) { - tmpPoly = p_psVectorFitPolynomial1DOrd(myPoly, myX, y, myYErr); + tmpPoly = p_psVectorFitPolynomial1DOrd(myPoly, x64, y64, yErr64); } else { // XXX: psErrorMsg() @@ -966,10 +950,15 @@ // Free any allocated memory. - if (mustFreeMyYErr == true) { - psFree(myYErr); - } - if (mustFreemyX == true) { - psFree(myX); - } + /* + if ((yErr == NULL) || (yErr->type.type != PS_TYPE_F64)) { + psFree(yErr64); + } + if ((x == NULL) || (x->type.type != PS_TYPE_F64)) { + psFree(x64); + } + if ((y == NULL) || (y->type.type != PS_TYPE_F64)) { + psFree(y64); + } + */ return(myPoly);