Changeset 1963 for trunk/psLib/src/math/psMinimize.c

Timestamp:

Oct 5, 2004, 12:47:21 PM (22 years ago)

Author:

gusciora

Message:

Added macros in psConstants.h to check for NULL and mis-typed vectors.

File:

• trunk/psLib/src/math/psMinimize.c (modified) (15 diffs)

trunk/psLib/src/math/psMinimize.c

@@ -9 +9 @@
  *  @author George Gusciora, MHPCC
  *
- *  @version $Revision: 1.52 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-10-05 01:43:58 $
+ *  @version $Revision: 1.53 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-10-05 22:47:21 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -56 +56 @@
  
 XXX: change name
+XXX: Use a static vector.
  *****************************************************************************/
 void p_psBuildSums1D(double x,
@@ -63 +64 @@
     int i = 0;
     double xSum = 0.0;
+
     if (sums == NULL) {
         sums = psVectorAlloc(polyOrder, PS_TYPE_F64);
+    }
+
+    if (polyOrder > sums->n) {
+        sums = psVectorRealloc(sums, polyOrder);
     }
 
@@ -77 +83 @@
 VectorNormalizeGen(): this routine returns a psVector with "x" elements.  The
 values of the vector will be scaled uniformly between -1.0 and 1.0.
+ 
+XXX: use a static vector.
 *****************************************************************************/
 psVector* VectorNormalizeGen(int x)
@@ -100 +108 @@
  
 XXX: This algorithm is derived from the Numerical Recipes.
+XXX: use recyvled vectors for internal data.
  *****************************************************************************/
 float *CalculateSecondDerivs(const psVector* restrict x,        ///< Ordinates (or NULL to just use the indices)
@@ -119 +128 @@
     float qn;
 
-    if (x == NULL) {
-        X = (float *) psAlloc(n * sizeof(float));
-        for(i=0;i<n;i++) {
-            X[i] = (float) i;
-        }
-        mustFreeX = true;
-    }
-
     // XXX: The second derivatives at the endpoints, undefined in the SDR,
     // are set in psConstants.h: LEFT_SPLINE_DERIV, RIGHT_SPLINE_DERIV.
@@ -167 +168 @@
     return(derivs2);
 }
-
 
 /******************************************************************************
@@ -212 +212 @@
 }
 
-
 /*****************************************************************************/
 /* FUNCTION IMPLEMENTATION - PUBLIC                                          */
@@ -258 +257 @@
     int i;
     float slope;
-
-    if (mySpline == NULL) {
-        //XXX psErrorMsg()
-    }
-    if (y == NULL) {
-        //XXX psErrorMsg()
-    }
+    psVector *x32 = NULL;
+    psVector *y32 = NULL;
+    psVector *yErr32 = NULL;
+    static psVector *x32Static = NULL;
+    static psVector *y32Static = NULL;
+    static psVector *yErr32Static = NULL;
+
+    PS_CHECK_NULL_PTR_RETURN_NULL(mySpline);
+    PS_CHECK_NULL_PTR_RETURN_NULL(y);
+    VECTOR_CHECK_TYPE_F32_OR_F64_RETURN_NULL(y);
+    CONVERT_VECTOR_F64_TO_F32(y, y32, y32Static);
+
+    // If yErr==NULL, set all errors equal.
+    if (yErr == NULL) {
+        GEN_PERSIST_YERR_VEC_F32(yErr32Static, y->n);
+        yErr32 = yErr32Static;
+    } else {
+        VECTOR_CHECK_TYPE_F32_OR_F64_RETURN_NULL(yErr);
+        CONVERT_VECTOR_F64_TO_F32(yErr, yErr32, yErr32Static);
+    }
+
+    // If x==NULL, create an x32 vector with x values set to (0:n).
+    if (x == NULL) {
+        GEN_PERSIST_X_INDEX_VEC_F32(x32Static, y->n);
+        x32 = x32Static;
+    } else {
+        VECTOR_CHECK_TYPE_F32_OR_F64_RETURN_NULL(x);
+        CONVERT_VECTOR_F64_TO_F32(x, x32, x32Static);
+    }
+    PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_NULL(x, y);
+    PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_NULL(yErr, y);
 
     if (y->n != (1 + mySpline->n)) {
-        psLogMsg(__func__, PS_LOG_WARN,
-                 "data size / spline size mismatch (%d %d)\n",
-                 y->n, mySpline->n);
-        // XXX: psErrorMsg()
+        psError(__func__, "data size / spline size mismatch (%d %d)\n",
+                y->n, mySpline->n);
+        return(NULL);
     }
 
@@ -295 +317 @@
     // Check if these are cubic splines (n==4).  If not, psError.
     if (4 != (mySpline->spline[0])->n) {
-        psLogMsg(__func__, PS_LOG_WARN,
-                 "Don't know how to generate %d-order splines.",
-                 (mySpline->spline[0])->n-1);
-        // XXX: psErrorMsg()
+        psError(__func__, "Don't know how to generate %d-order splines.",
+                (mySpline->spline[0])->n-1);
         return(NULL);
     }
@@ -315 +335 @@
         H = x->data.F32[i+1] - x->data.F32[i];
         psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
-                "X data (%f - %f) (%f)\n", x->data.F32[i], x->data.F32[i+1], H);
+                "x data (%f - %f) (%f)\n", x->data.F32[i], x->data.F32[i+1], H);
         //
         // ******** Calculate 0-order term ********
@@ -630 +650 @@
 }
 
-#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); \
-
-
 /******************************************************************************
 p_psVectorFitPolynomial1DCheb():  This routine will fit a Chebyshev
@@ -854 +867 @@
     int i;
     psPolynomial1D *tmpPoly;
-    PS_CHECK_NULL_1DPOLY_ACTION(myPoly, 0);
-    PS_CHECK_NULL_VECTOR_ACTION(y, 0);
-    PS_CHECK_EMPTY_VECTOR_ACTION(y, 0);
+    PS_CHECK_NULL_1DPOLY_RETURN_NULL(myPoly);
+    PS_CHECK_NULL_VECTOR_RETURN_NULL(y);
+    PS_CHECK_EMPTY_VECTOR_RETURN_NULL(y);
     psVector *x64 = NULL;
     psVector *y64 = NULL;
@@ -864 +877 @@
     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);
-    }
-
+    VECTOR_CHECK_TYPE_F32_OR_F64_RETURN_NULL(y);
+    CONVERT_VECTOR_F32_TO_F64(y, y64, y64Static);
 
     // If yErr==NULL, set all errors equal.
     if (yErr == NULL) {
-        yErr64Static = psVectorRecycle(yErr64Static, y->n, PS_TYPE_F64);
-        p_psMemSetPersistent(yErr64Static, true);
-        p_psMemSetPersistent(yErr64Static->data.V, true);
-
-        for (i=0;i<yErr64Static->n;i++) {
-            yErr64Static->data.F64[i] = 1.0;
-        }
+        GEN_PERSIST_YERR_VEC_F64(yErr64Static, y->n);
         yErr64 = yErr64Static;
     } else {
-        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.
+        VECTOR_CHECK_TYPE_F32_OR_F64_RETURN_NULL(yErr);
+        CONVERT_VECTOR_F32_TO_F64(yErr, yErr64, yErr64Static);
+    }
+
+    // If x==NULL, create an x64 vector with x values set to (0:n).
     if (x == NULL) {
-        x64Static = psVectorRecycle(x64Static, y->n, PS_TYPE_F64);
-        p_psMemSetPersistent(x64Static, true);
-        p_psMemSetPersistent(x64Static->data.V, true);
+        GEN_PERSIST_X_INDEX_VEC_F64(x64Static, y->n);
+        if (myPoly->type == PS_POLYNOMIAL_CHEB) {
+            p_psNormalizeVector(x64Static);
+        }
         x64 = x64Static;
-
-        if (myPoly->type == PS_POLYNOMIAL_ORD) {
-            for (i=0;i<x64->n;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;i<x64->n;i++) {
-                x64->data.F64[i] = (((double) i) - 0.5 * (min + max)) /
-                                   (0.5 * (max - min));
-            }
-        }
     } else {
-        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);
+        VECTOR_CHECK_TYPE_F32_OR_F64_RETURN_NULL(x);
+        CONVERT_VECTOR_F32_TO_F64(x, x64, x64Static);
+    }
+    PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_NULL(x64, y64);
+    PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_NULL(yErr64, y64);
 
     // Call the appropriate vector fitting routine.
     if (myPoly->type == PS_POLYNOMIAL_CHEB) {
-        tmpPoly = p_psVectorFitPolynomial1DCheby(myPoly, x64, y64, yErr64);
+        p_psVectorFitPolynomial1DCheby(myPoly, x64, y64, yErr64);
     } else if (myPoly->type == PS_POLYNOMIAL_ORD) {
         tmpPoly = p_psVectorFitPolynomial1DOrd(myPoly, x64, y64, yErr64);
@@ -948 +912 @@
         return(NULL);
     }
-
-    // Free any allocated memory.
-    /*
-        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);