IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 1953


Ignore:
Timestamp:
Oct 4, 2004, 3:43:58 PM (22 years ago)
Author:
gusciora
Message:

Added F32 types for the VectorFit routine.

Location:
trunk/psLib/src
Files:
11 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/collections/psVector.h

    r1920 r1953  
    1 
    21/** @file  psVector.h
    32 *
    … …  
    1211 *  @author Ross Harman, MHPCC
    1312 *
    14  *  @version $Revision: 1.22 $ $Name: not supported by cvs2svn $
    15  *  @date $Date: 2004-09-28 23:26:48 $
     13 *  @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
     14 *  @date $Date: 2004-10-05 01:43:57 $
    1615 *
    1716 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii

  • trunk/psLib/src/dataManip/psConstants.h

    r1945 r1953  
    66 *  @author GLG, MHPCC
    77 *
    8  *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2004-10-03 23:35:47 $
     8 *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2004-10-05 01:43:58 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    1818#define RIGHT_SPLINE_DERIV 0.0
    1919
     20#define GEN_STATIC_RECYCLED_VECTOR(NAME, SIZE, TYPE) \
     21static psVector *NAME = NULL; \
     22NAME = psVectorRecycle(NAME, SIZE, TYPE); \
     23p_psMemSetPersistent(NAME, true); \
     24p_psMemSetPersistent(NAME->data.V, true); \
     25
     26
     27
     28
    2029/** Preprocessor macro to generate error on an incorrect type */
    2130#define PS_CHECK_VECTOR_TYPE(NAME, TYPE) \
    2231if (NAME->type.type != TYPE) { \
    23     psError(__func__,"Invalid operation: %s has incorrect type.", #NAME); \
     32    psError(__func__,"Unallowable operation: %s has incorrect type.", #NAME); \
    2433}
    2534
    … …  
    2736#define PS_CHECK_NULL_VECTOR(NAME) \
    2837if (NAME == NULL || NAME->data.V == NULL) { \
    29     psError(__func__,"Invalid operation: %s or its data is NULL.", #NAME); \
     38    psError(__func__,"Unallowable operation: %s or its data is NULL.", #NAME); \
    3039}
    3140
    … …  
    3443if (NAME == NULL || NAME->data.V == NULL) { \
    3544    if (ACTION == 0) { \
    36         psError(__func__,"Invalid operation: %s or its data is NULL.", #NAME); \
     45        psError(__func__,"Unallowable operation: %s or its data is NULL.", #NAME); \
    3746        return(NULL); \
    3847    } \
    3948    if (ACTION == 1) { \
    40         psAbort(__func__, "Invalid operation: %s or its data is NULL.", #NAME); \
     49        psAbort(__func__, "Unallowable operation: %s or its data is NULL.", #NAME); \
    4150    } \
    4251}
    … …  
    4554#define PS_CHECK_NULL_PTR(NAME) \
    4655if (NAME == NULL) { \
    47     psError(__func__,"Invalid operation: %s is NULL.", #NAME); \
     56    psError(__func__,"Unallowable operation: %s is NULL.", #NAME); \
    4857}
    4958
    … …  
    5160#define PS_CHECK_EMPTY_VECTOR(NAME) \
    5261if (NAME->n < 1) { \
    53     psError(__func__,"Invalid operation: %s has zero n value.", #NAME); \
     62    psError(__func__,"Unallowable operation: %s has zero n value.", #NAME); \
    5463}
    5564
    … …  
    5867if (NAME->n < 1) { \
    5968    if (ACTION == 0) { \
    60         psError(__func__,"Invalid operation: %s has zero n value.", #NAME); \
     69        psError(__func__,"Unallowable operation: %s has zero n value.", #NAME); \
    6170        return(NULL); \
    6271    } \
    6372    if (ACTION == 1) { \
    64         psAbort(__func__,"Invalid operation: %s has zero n value.", #NAME); \
     73        psAbort(__func__,"Unallowable operation: %s has zero n value.", #NAME); \
    6574    } \
    6675}
    … …  
    7584#define PS_CHECK_NULL_IMAGE(NAME)                                                           \
    7685if (NAME == NULL || NAME->data.V == NULL) {                                                         \
    77     psError(__func__,"Invalid operation: %s or its data is NULL.", #NAME);                          \
     86    psError(__func__,"Unallowable operation: %s or its data is NULL.", #NAME);                          \
    7887}
    7988
    … …  
    8190#define PS_CHECK_EMPTY_IMAGE(NAME)                                                           \
    8291if (NAME->numCols < 1 || NAME->numRows < 1) {                                                       \
    83     psError(__func__,"Invalid operation: %s has zero rows or columns (%dx%d).", #NAME,              \
     92    psError(__func__,"Unallowable operation: %s has zero rows or columns (%dx%d).", #NAME,              \
    8493            NAME->numCols, NAME->numRows);                                                          \
    8594}
    … …  
    8897#define PS_CHECK_NULL_1DPOLY(NAME)                                                          \
    8998if (NAME == NULL || NAME->coeff == NULL) {                                                         \
    90     psError(__func__,"Invalid operation: %s or its data is NULL.", #NAME);                          \
     99    psError(__func__,"Unallowable operation: %s or its data is NULL.", #NAME);                          \
    91100}
    92101
    … …  
    95104if (NAME == NULL || NAME->coeff == NULL) {                                                         \
    96105    if (ACTION == 0) { \
    97         psError(__func__,"Invalid operation: %s or its coeffs is NULL.", #NAME);                          \
     106        psError(__func__,"Unallowable operation: %s or its coeffs is NULL.", #NAME);                          \
    98107        return(NULL); \
    99108    } \
    100109    if (ACTION == 1) { \
    101         psAbort(__func__,"Invalid operation: %s or its coeffs is NULL.", #NAME);                          \
     110        psAbort(__func__,"Unallowable operation: %s or its coeffs is NULL.", #NAME);                          \
    102111    } \
    103112}

  • trunk/psLib/src/dataManip/psFunctions.c

    r1947 r1953  
    77 *  polynomials.  It also contains a Gaussian functions.
    88 *
    9  *  @version $Revision: 1.45 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2004-10-04 00:40:23 $
     9 *  @version $Revision: 1.46 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2004-10-05 01:43:58 $
    1111 *
    1212 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    19881988 
    19891989XXX: nobody asked for us to generate this routine.
     1990 
     1991XXX: add another psScalar argument for the result.
     1992 
     1993XXX: The VectorCopy routines seg fault when I declare range32 as static.
    19901994 *****************************************************************************/
    19911995psScalar *p_psVectorInterpolate(psVector *domain,
    … …  
    19941998                                psScalar *x)
    19951999{
     2000    psVector *range32 = NULL;
     2001    psVector *domain32 = NULL;
    19962002    psTrace(".psLib.dataManip.psFunctions.p_psVectorInterpolate", 4,
    19972003            "---- p_psVectorInterpolate() begin ----\n");
    … …  
    20002006            (domain->type.type != x->type.type)) {
    20012007        // XXX psError
    2002         printf("domain->type.type != range->type.type != x->type.type\n");
     2008        printf("p_psVectorInterpolate(): domain->type.type != range->type.type != x->type.type\n");
    20032009        printf("%d %d %d\n", domain->type.type, range->type.type, x->type.type);
     2010        exit(1);
    20042011    }
    20052012    if (domain->n != range->n) {
    20062013        // XXX psError
    20072014        printf("domain->n != range->n\n");
     2015        exit(1);
    20082016    }
    20092017    if (order > (domain->n - 1)) {
    20102018        // XXX psError: not enough data points for order-order interpolation.
    20112019        printf("not enough data points for order-order interpolation.\n");
     2020        exit(1);
    20122021    }
    20132022
    … …  
    20202029                             order,
    20212030                             x->data.F32), PS_TYPE_F32));
     2031    } else if (x->type.type == PS_TYPE_F64) {
     2032        // XXX: use recycled vectors here.
     2033        range32 = psVectorCopy(range32, range, PS_TYPE_F32);
     2034        domain32 = psVectorCopy(domain32, domain, PS_TYPE_F32);
     2035
     2036        psScalar *tmpScalar = psScalarAlloc((double)
     2037                                            p_ps1DInterpolateF32(domain32->data.F32,
     2038                                                                 range32->data.F32,
     2039                                                                 domain32->n,
     2040                                                                 order,
     2041                                                                 (float) x->data.F64), PS_TYPE_F64);
     2042        psFree(range32);
     2043        psFree(domain32);
     2044
     2045        psTrace(".psLib.dataManip.psFunctions.p_psVectorInterpolate", 4,
     2046                "---- p_psVectorInterpolate() end ----\n");
     2047        // XXX: Convert data type to F64?
     2048        return(tmpScalar);
     2049
    20222050    } else {
    20232051        // XXX psError: type not supported

  • trunk/psLib/src/dataManip/psMinimize.c

    r1945 r1953  
    99 *  @author George Gusciora, MHPCC
    1010 *
    11  *  @version $Revision: 1.51 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-10-03 23:35:47 $
     11 *  @version $Revision: 1.52 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-10-05 01:43:58 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    2424/* INCLUDE FILES                                                             */
    2525/*****************************************************************************/
    26 #include <stdlib.h>
    27 #include <stdio.h>
    28 #include <string.h>
    29 #include <stdarg.h>
    30 #include <float.h>
    31 #include <math.h>
    32 
    33 #include <gsl/gsl_multifit_nlin.h>
    34 #include <gsl/gsl_multimin.h>
    35 #include <gsl/gsl_rng.h>
    36 #include <gsl/gsl_randist.h>
    37 #include <gsl/gsl_vector.h>
    38 #include <gsl/gsl_blas.h>
    39 
    40 #include "psMemory.h"
    41 #include "psVector.h"
    42 #include "psImage.h"
    43 #include "psTrace.h"
    44 #include "psLogMsg.h"
    45 #include "psError.h"
    46 #include "psAbort.h"
    47 #include "psFunctions.h"
    4826#include "psMinimize.h"
    49 #include "psMatrix.h"
    50 #include "psConstants.h"
    51 
    5227/*****************************************************************************/
    5328/* DEFINE STATEMENTS                                                         */
    … …  
    655630}
    656631
     632#define GEN_STATIC_RECYCLED_VECTOR(NAME, SIZE, TYPE) \
     633static psVector *NAME = NULL; \
     634NAME = psVectorRecycle(NAME, SIZE, TYPE); \
     635p_psMemSetPersistent(NAME, true); \
     636p_psMemSetPersistent(NAME->data.V, true); \
     637
    657638
    658639/******************************************************************************
    … …  
    661642coefficients of that polynomial.
    662643 
    663 XXX: yErr is currently ignored.
    664  
    665 XXX: must add type F32 (currently F64 only).
    666 
    667 XXX: Use private name?
    668  *****************************************************************************/
     644XXX:
     645yErr is currently ignored.
     646 
     647XXX:
     648Use private name?
     649*****************************************************************************/
    669650psPolynomial1D *p_psVectorFitPolynomial1DCheby(psPolynomial1D* myPoly,
    670651        const psVector* restrict x,
    … …  
    675656    int k;
    676657    int n = x->n;
    677     psVector *f = psVectorAlloc(n, PS_TYPE_F64);
    678658    double fac;
    679659    double sum;
    680     // XXX: Use static memory here.
    681     psScalar *tmpScalar = psScalarAlloc(0.0, PS_TYPE_F32);
     660    GEN_STATIC_RECYCLED_VECTOR(f, n, PS_TYPE_F64);
    682661    psScalar *fScalar;
     662    psScalar tmpScalar;
     663    tmpScalar.type.type = PS_TYPE_F64;
    683664
    684665    // XXX: These assignments appear too simple to warrant code and
    … …  
    690671    double bpa = 0.5 * (max+min);  // 0
    691672
    692     // XXX: Eliminate this later by generating a F64 version of the
    693     // LaGrange interpolation routines.
    694     PS_VECTOR_F64_TO_F32(x, x32);
    695     PS_VECTOR_F64_TO_F32(y, y32);
    696 
    697673    // In this loop, we first calculate the values of X for which the
    698674    // Chebyshev polynomials are zero (see NR, section 5.4).  Then we
    … …  
    707683        double Y = cos(M_PI * (0.5 + ((float) i)) / ((float) n));
    708684        double X = (Y + bma + bpa) - 1.0;
    709         tmpScalar->data.F32 = (float) X;
     685        tmpScalar.data.F64 = X;
    710686
    711687        // We interpolate against are tabluated x,y vectors to determine the
    712688        // function value at X.
    713         fScalar = p_psVectorInterpolate((psVector *) x32, (psVector *) y32,
    714                                         3, tmpScalar);
    715 
    716         f->data.F64[i] = (double) fScalar->data.F32;
     689        fScalar = p_psVectorInterpolate((psVector *) x,
     690                                        (psVector *) y,
     691                                        3,
     692                                        &tmpScalar);
     693
     694        f->data.F64[i] = fScalar->data.F64;
    717695        psFree(fScalar);
    718696
    … …  
    724702    // We have the values for f() at the zero points, we now calculate the
    725703    // coefficients of the Chebyshev polynomial: NR 5.8.7.
     704
    726705    fac = 2.0/((float) n);
    727 
    728706    // XXX: is this loop bound correct?
    729707    for (j=0;j<myPoly->n;j++) {
    … …  
    737715    }
    738716
    739     // XXX: Must free memory.
    740     psFree(f);
    741     psFree(x32);
    742     psFree(y32);
    743     psFree(tmpScalar);
    744 
    745717    return(myPoly);
    746718}
    … …  
    750722polynomial of degree myPoly to the data points (x, y) and return the
    751723coefficients of that polynomial.
    752  
    753 XXX: must add type F32 (currently F64 only).
    754724 
    755725XXX: Use private name?
    … …  
    875845polynomial.
    876846 
    877 XXX: must add type F32 (currently F64 only).
     847XXX: type F32 is done via vector conversion only.
    878848 *****************************************************************************/
    879849psPolynomial1D* psVectorFitPolynomial1D(psPolynomial1D* myPoly,
    … …  
    882852                                        const psVector* restrict yErr)
    883853{
    884     bool mustFreeMyYErr = false;
    885     bool mustFreemyX = false;
    886854    int i;
    887855    psPolynomial1D *tmpPoly;
    888     psVector *myX = NULL;
    889     psVector *myYErr = NULL;
    890 
    891856    PS_CHECK_NULL_1DPOLY_ACTION(myPoly, 0);
    892857    PS_CHECK_NULL_VECTOR_ACTION(y, 0);
    893858    PS_CHECK_EMPTY_VECTOR_ACTION(y, 0);
     859    psVector *x64 = NULL;
     860    psVector *y64 = NULL;
     861    psVector *yErr64 = NULL;
     862    static psVector *x64Static = NULL;
     863    static psVector *y64Static = NULL;
     864    static psVector *yErr64Static = NULL;
     865
     866    if (y->type.type == PS_TYPE_F64) {
     867        y64 = (psVector *) y;
     868    } else if (y->type.type == PS_TYPE_F32) {
     869        y64Static = psVectorRecycle(y64Static, y->n, PS_TYPE_F64);
     870        p_psMemSetPersistent(y64Static, true);
     871        p_psMemSetPersistent(y64Static->data.V, true);
     872        y64 = y64Static;
     873    } else {
     874        // XXX: psError() bad type.
     875        psAbort(__func__, "Bad type for y64 (%d)", y->type.type);
     876    }
     877
    894878
    895879    // If yErr==NULL, set all errors equal.
    896880    if (yErr == NULL) {
    897         myYErr = psVectorAlloc(y->n, y->type.type);
    898         mustFreeMyYErr = true;
    899 
    900         if (y->type.type == PS_TYPE_F32) {
    901             for (i=0;i<myYErr->n;i++) {
    902                 myYErr->data.F32[i] = 1.0;
    903             }
    904         } else if (y->type.type == PS_TYPE_F64) {
    905             for (i=0;i<myYErr->n;i++) {
    906                 myYErr->data.F64[i] = 1.0;
    907             }
    908         }
     881        yErr64Static = psVectorRecycle(yErr64Static, y->n, PS_TYPE_F64);
     882        p_psMemSetPersistent(yErr64Static, true);
     883        p_psMemSetPersistent(yErr64Static->data.V, true);
     884
     885        for (i=0;i<yErr64Static->n;i++) {
     886            yErr64Static->data.F64[i] = 1.0;
     887        }
     888        yErr64 = yErr64Static;
    909889    } else {
    910         myYErr = (psVector *) yErr;
    911     }
    912 
    913     // If x==NULL, create an myX vector with x values set to (0:n), and if
     890        if (yErr->type.type == PS_TYPE_F64) {
     891            yErr64 = (psVector *) yErr;
     892        } else if (yErr->type.type == PS_TYPE_F32) {
     893            yErr64Static = psVectorRecycle(yErr64Static, yErr->n, PS_TYPE_F64);
     894            p_psMemSetPersistent(yErr64Static, true);
     895            p_psMemSetPersistent(yErr64Static->data.V, true);
     896        } else {
     897            // XXX: psError() bad type.
     898            psAbort(__func__, "Bad type for yErr64");
     899        }
     900    }
     901
     902    // If x==NULL, create an x64 vector with x values set to (0:n), and if
    914903    // this is a Chebyshev polynomial, we must scale to (-1:1).
    915 
    916904    // XXX: Verify that this is the correct action.
    917905    if (x == NULL) {
    918         myX = psVectorAlloc(y->n, y->type.type);
    919         mustFreemyX = true;
    920 
    921         if (y->type.type == PS_TYPE_F32) {
    922             if (myPoly->type == PS_POLYNOMIAL_ORD) {
    923                 for (i=0;i<myX->n;i++) {
    924                     myX->data.F32[i] = (float) i;
    925                 }
    926             } else if (myPoly->type == PS_POLYNOMIAL_CHEB) {
    927                 float min = 0.0;
    928                 float max = (float) (y->n - 1);
    929 
    930                 for (i=0;i<myX->n;i++) {
    931                     myX->data.F32[i] = (((float) i) - 0.5 * (min + max)) /
    932                                        (0.5 * (max - min));
    933                 }
    934             }
    935         } else if (y->type.type == PS_TYPE_F64) {
    936             if (myPoly->type == PS_POLYNOMIAL_ORD) {
    937                 for (i=0;i<myX->n;i++) {
    938                     myX->data.F64[i] = (float) i;
    939                 }
    940             } else if (myPoly->type == PS_POLYNOMIAL_CHEB) {
    941                 double min = 0.0;
    942                 double max = (double) (y->n - 1);
    943 
    944                 for (i=0;i<myX->n;i++) {
    945                     myX->data.F64[i] = (((float) i) - 0.5 * (min + max)) /
    946                                        (0.5 * (max - min));
    947                 }
     906        x64Static = psVectorRecycle(x64Static, y->n, PS_TYPE_F64);
     907        p_psMemSetPersistent(x64Static, true);
     908        p_psMemSetPersistent(x64Static->data.V, true);
     909        x64 = x64Static;
     910
     911        if (myPoly->type == PS_POLYNOMIAL_ORD) {
     912            for (i=0;i<x64->n;i++) {
     913                x64->data.F64[i] = (float) i;
     914            }
     915        } else if (myPoly->type == PS_POLYNOMIAL_CHEB) {
     916            double min = 0.0;
     917            double max = (double) (y->n - 1);
     918
     919            for (i=0;i<x64->n;i++) {
     920                x64->data.F64[i] = (((double) i) - 0.5 * (min + max)) /
     921                                   (0.5 * (max - min));
    948922            }
    949923        }
    950924    } else {
    951         myX = (psVector *) x;
    952     }
    953 
    954     PS_CHECK_VECTOR_SIZE_EQUAL(y, myX);
    955     PS_CHECK_VECTOR_SIZE_EQUAL(y, myYErr);
     925        if (x->type.type == PS_TYPE_F64) {
     926            x64 = (psVector *) x;
     927        } else if (x->type.type == PS_TYPE_F32) {
     928            x64Static = psVectorRecycle(x64Static, x->n, PS_TYPE_F64);
     929            p_psMemSetPersistent(x64Static, true);
     930            p_psMemSetPersistent(x64Static->data.V, true);
     931            x64 = x64Static;
     932        } else {
     933            // XXX: psError() bad type.
     934            psAbort(__func__, "Bad type for x64");
     935        }
     936    }
     937
     938    PS_CHECK_VECTOR_SIZE_EQUAL(y64, x64);
     939    PS_CHECK_VECTOR_SIZE_EQUAL(y64, yErr64);
    956940
    957941    // Call the appropriate vector fitting routine.
    958942    if (myPoly->type == PS_POLYNOMIAL_CHEB) {
    959         tmpPoly = p_psVectorFitPolynomial1DCheby(myPoly, myX, y, myYErr);
     943        tmpPoly = p_psVectorFitPolynomial1DCheby(myPoly, x64, y64, yErr64);
    960944    } else if (myPoly->type == PS_POLYNOMIAL_ORD) {
    961         tmpPoly = p_psVectorFitPolynomial1DOrd(myPoly, myX, y, myYErr);
     945        tmpPoly = p_psVectorFitPolynomial1DOrd(myPoly, x64, y64, yErr64);
    962946    } else {
    963947        // XXX: psErrorMsg()
    … …  
    966950
    967951    // Free any allocated memory.
    968     if (mustFreeMyYErr == true) {
    969         psFree(myYErr);
    970     }
    971     if (mustFreemyX == true) {
    972         psFree(myX);
    973     }
     952    /*
     953        if ((yErr == NULL) || (yErr->type.type != PS_TYPE_F64)) {
     954            psFree(yErr64);
     955        }
     956        if ((x == NULL) || (x->type.type != PS_TYPE_F64)) {
     957            psFree(x64);
     958        }
     959        if ((y == NULL) || (y->type.type != PS_TYPE_F64)) {
     960            psFree(y64);
     961        }
     962    */
    974963
    975964    return(myPoly);

  • trunk/psLib/src/dataManip/psMinimize.h

    r1846 r1953  
    88 *  @author George Gusciora, MHPCC
    99 *
    10  *  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
    11  *  @date $Date: 2004-09-22 01:30:21 $
     10 *  @version $Revision: 1.26 $ $Name: not supported by cvs2svn $
     11 *  @date $Date: 2004-10-05 01:43:58 $
    1212 *
    1313 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    1818#define PS_MINIMIZE_H
    1919
     20#include <stdio.h>
     21#include <float.h>
     22#include <math.h>
     23
     24#include "psMemory.h"
     25#include "psVector.h"
     26#include "psArray.h"
     27#include "psImage.h"
     28#include "psMatrix.h"
    2029#include "psFunctions.h"
    21 #include "psArray.h"
    22 #include "psVector.h"
    23 #include "psMatrix.h"
     30#include "psTrace.h"
     31#include "psLogMsg.h"
     32#include "psError.h"
     33#include "psAbort.h"
     34#include "psConstants.h"
     35
    2436
    2537/** \file psMinimize.h
    … …  
    4456                                    float tol);
    4557
    46 /** Derive a polynomial fit by chi^2 minimisation (analytically)
     58/** Derive a polynomial fit.
    4759 *
    4860 *  psVectorFitPolynomial1d returns the polynomial that best fits the

  • trunk/psLib/src/math/psConstants.h

    r1945 r1953  
    66 *  @author GLG, MHPCC
    77 *
    8  *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2004-10-03 23:35:47 $
     8 *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2004-10-05 01:43:58 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    1818#define RIGHT_SPLINE_DERIV 0.0
    1919
     20#define GEN_STATIC_RECYCLED_VECTOR(NAME, SIZE, TYPE) \
     21static psVector *NAME = NULL; \
     22NAME = psVectorRecycle(NAME, SIZE, TYPE); \
     23p_psMemSetPersistent(NAME, true); \
     24p_psMemSetPersistent(NAME->data.V, true); \
     25
     26
     27
     28
    2029/** Preprocessor macro to generate error on an incorrect type */
    2130#define PS_CHECK_VECTOR_TYPE(NAME, TYPE) \
    2231if (NAME->type.type != TYPE) { \
    23     psError(__func__,"Invalid operation: %s has incorrect type.", #NAME); \
     32    psError(__func__,"Unallowable operation: %s has incorrect type.", #NAME); \
    2433}
    2534
    … …  
    2736#define PS_CHECK_NULL_VECTOR(NAME) \
    2837if (NAME == NULL || NAME->data.V == NULL) { \
    29     psError(__func__,"Invalid operation: %s or its data is NULL.", #NAME); \
     38    psError(__func__,"Unallowable operation: %s or its data is NULL.", #NAME); \
    3039}
    3140
    … …  
    3443if (NAME == NULL || NAME->data.V == NULL) { \
    3544    if (ACTION == 0) { \
    36         psError(__func__,"Invalid operation: %s or its data is NULL.", #NAME); \
     45        psError(__func__,"Unallowable operation: %s or its data is NULL.", #NAME); \
    3746        return(NULL); \
    3847    } \
    3948    if (ACTION == 1) { \
    40         psAbort(__func__, "Invalid operation: %s or its data is NULL.", #NAME); \
     49        psAbort(__func__, "Unallowable operation: %s or its data is NULL.", #NAME); \
    4150    } \
    4251}
    … …  
    4554#define PS_CHECK_NULL_PTR(NAME) \
    4655if (NAME == NULL) { \
    47     psError(__func__,"Invalid operation: %s is NULL.", #NAME); \
     56    psError(__func__,"Unallowable operation: %s is NULL.", #NAME); \
    4857}
    4958
    … …  
    5160#define PS_CHECK_EMPTY_VECTOR(NAME) \
    5261if (NAME->n < 1) { \
    53     psError(__func__,"Invalid operation: %s has zero n value.", #NAME); \
     62    psError(__func__,"Unallowable operation: %s has zero n value.", #NAME); \
    5463}
    5564
    … …  
    5867if (NAME->n < 1) { \
    5968    if (ACTION == 0) { \
    60         psError(__func__,"Invalid operation: %s has zero n value.", #NAME); \
     69        psError(__func__,"Unallowable operation: %s has zero n value.", #NAME); \
    6170        return(NULL); \
    6271    } \
    6372    if (ACTION == 1) { \
    64         psAbort(__func__,"Invalid operation: %s has zero n value.", #NAME); \
     73        psAbort(__func__,"Unallowable operation: %s has zero n value.", #NAME); \
    6574    } \
    6675}
    … …  
    7584#define PS_CHECK_NULL_IMAGE(NAME)                                                           \
    7685if (NAME == NULL || NAME->data.V == NULL) {                                                         \
    77     psError(__func__,"Invalid operation: %s or its data is NULL.", #NAME);                          \
     86    psError(__func__,"Unallowable operation: %s or its data is NULL.", #NAME);                          \
    7887}
    7988
    … …  
    8190#define PS_CHECK_EMPTY_IMAGE(NAME)                                                           \
    8291if (NAME->numCols < 1 || NAME->numRows < 1) {                                                       \
    83     psError(__func__,"Invalid operation: %s has zero rows or columns (%dx%d).", #NAME,              \
     92    psError(__func__,"Unallowable operation: %s has zero rows or columns (%dx%d).", #NAME,              \
    8493            NAME->numCols, NAME->numRows);                                                          \
    8594}
    … …  
    8897#define PS_CHECK_NULL_1DPOLY(NAME)                                                          \
    8998if (NAME == NULL || NAME->coeff == NULL) {                                                         \
    90     psError(__func__,"Invalid operation: %s or its data is NULL.", #NAME);                          \
     99    psError(__func__,"Unallowable operation: %s or its data is NULL.", #NAME);                          \
    91100}
    92101
    … …  
    95104if (NAME == NULL || NAME->coeff == NULL) {                                                         \
    96105    if (ACTION == 0) { \
    97         psError(__func__,"Invalid operation: %s or its coeffs is NULL.", #NAME);                          \
     106        psError(__func__,"Unallowable operation: %s or its coeffs is NULL.", #NAME);                          \
    98107        return(NULL); \
    99108    } \
    100109    if (ACTION == 1) { \
    101         psAbort(__func__,"Invalid operation: %s or its coeffs is NULL.", #NAME);                          \
     110        psAbort(__func__,"Unallowable operation: %s or its coeffs is NULL.", #NAME);                          \
    102111    } \
    103112}

  • trunk/psLib/src/math/psMinimize.c

    r1945 r1953  
    99 *  @author George Gusciora, MHPCC
    1010 *
    11  *  @version $Revision: 1.51 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-10-03 23:35:47 $
     11 *  @version $Revision: 1.52 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-10-05 01:43:58 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    2424/* INCLUDE FILES                                                             */
    2525/*****************************************************************************/
    26 #include <stdlib.h>
    27 #include <stdio.h>
    28 #include <string.h>
    29 #include <stdarg.h>
    30 #include <float.h>
    31 #include <math.h>
    32 
    33 #include <gsl/gsl_multifit_nlin.h>
    34 #include <gsl/gsl_multimin.h>
    35 #include <gsl/gsl_rng.h>
    36 #include <gsl/gsl_randist.h>
    37 #include <gsl/gsl_vector.h>
    38 #include <gsl/gsl_blas.h>
    39 
    40 #include "psMemory.h"
    41 #include "psVector.h"
    42 #include "psImage.h"
    43 #include "psTrace.h"
    44 #include "psLogMsg.h"
    45 #include "psError.h"
    46 #include "psAbort.h"
    47 #include "psFunctions.h"
    4826#include "psMinimize.h"
    49 #include "psMatrix.h"
    50 #include "psConstants.h"
    51 
    5227/*****************************************************************************/
    5328/* DEFINE STATEMENTS                                                         */
    … …  
    655630}
    656631
     632#define GEN_STATIC_RECYCLED_VECTOR(NAME, SIZE, TYPE) \
     633static psVector *NAME = NULL; \
     634NAME = psVectorRecycle(NAME, SIZE, TYPE); \
     635p_psMemSetPersistent(NAME, true); \
     636p_psMemSetPersistent(NAME->data.V, true); \
     637
    657638
    658639/******************************************************************************
    … …  
    661642coefficients of that polynomial.
    662643 
    663 XXX: yErr is currently ignored.
    664  
    665 XXX: must add type F32 (currently F64 only).
    666 
    667 XXX: Use private name?
    668  *****************************************************************************/
     644XXX:
     645yErr is currently ignored.
     646 
     647XXX:
     648Use private name?
     649*****************************************************************************/
    669650psPolynomial1D *p_psVectorFitPolynomial1DCheby(psPolynomial1D* myPoly,
    670651        const psVector* restrict x,
    … …  
    675656    int k;
    676657    int n = x->n;
    677     psVector *f = psVectorAlloc(n, PS_TYPE_F64);
    678658    double fac;
    679659    double sum;
    680     // XXX: Use static memory here.
    681     psScalar *tmpScalar = psScalarAlloc(0.0, PS_TYPE_F32);
     660    GEN_STATIC_RECYCLED_VECTOR(f, n, PS_TYPE_F64);
    682661    psScalar *fScalar;
     662    psScalar tmpScalar;
     663    tmpScalar.type.type = PS_TYPE_F64;
    683664
    684665    // XXX: These assignments appear too simple to warrant code and
    … …  
    690671    double bpa = 0.5 * (max+min);  // 0
    691672
    692     // XXX: Eliminate this later by generating a F64 version of the
    693     // LaGrange interpolation routines.
    694     PS_VECTOR_F64_TO_F32(x, x32);
    695     PS_VECTOR_F64_TO_F32(y, y32);
    696 
    697673    // In this loop, we first calculate the values of X for which the
    698674    // Chebyshev polynomials are zero (see NR, section 5.4).  Then we
    … …  
    707683        double Y = cos(M_PI * (0.5 + ((float) i)) / ((float) n));
    708684        double X = (Y + bma + bpa) - 1.0;
    709         tmpScalar->data.F32 = (float) X;
     685        tmpScalar.data.F64 = X;
    710686
    711687        // We interpolate against are tabluated x,y vectors to determine the
    712688        // function value at X.
    713         fScalar = p_psVectorInterpolate((psVector *) x32, (psVector *) y32,
    714                                         3, tmpScalar);
    715 
    716         f->data.F64[i] = (double) fScalar->data.F32;
     689        fScalar = p_psVectorInterpolate((psVector *) x,
     690                                        (psVector *) y,
     691                                        3,
     692                                        &tmpScalar);
     693
     694        f->data.F64[i] = fScalar->data.F64;
    717695        psFree(fScalar);
    718696
    … …  
    724702    // We have the values for f() at the zero points, we now calculate the
    725703    // coefficients of the Chebyshev polynomial: NR 5.8.7.
     704
    726705    fac = 2.0/((float) n);
    727 
    728706    // XXX: is this loop bound correct?
    729707    for (j=0;j<myPoly->n;j++) {
    … …  
    737715    }
    738716
    739     // XXX: Must free memory.
    740     psFree(f);
    741     psFree(x32);
    742     psFree(y32);
    743     psFree(tmpScalar);
    744 
    745717    return(myPoly);
    746718}
    … …  
    750722polynomial of degree myPoly to the data points (x, y) and return the
    751723coefficients of that polynomial.
    752  
    753 XXX: must add type F32 (currently F64 only).
    754724 
    755725XXX: Use private name?
    … …  
    875845polynomial.
    876846 
    877 XXX: must add type F32 (currently F64 only).
     847XXX: type F32 is done via vector conversion only.
    878848 *****************************************************************************/
    879849psPolynomial1D* psVectorFitPolynomial1D(psPolynomial1D* myPoly,
    … …  
    882852                                        const psVector* restrict yErr)
    883853{
    884     bool mustFreeMyYErr = false;
    885     bool mustFreemyX = false;
    886854    int i;
    887855    psPolynomial1D *tmpPoly;
    888     psVector *myX = NULL;
    889     psVector *myYErr = NULL;
    890 
    891856    PS_CHECK_NULL_1DPOLY_ACTION(myPoly, 0);
    892857    PS_CHECK_NULL_VECTOR_ACTION(y, 0);
    893858    PS_CHECK_EMPTY_VECTOR_ACTION(y, 0);
     859    psVector *x64 = NULL;
     860    psVector *y64 = NULL;
     861    psVector *yErr64 = NULL;
     862    static psVector *x64Static = NULL;
     863    static psVector *y64Static = NULL;
     864    static psVector *yErr64Static = NULL;
     865
     866    if (y->type.type == PS_TYPE_F64) {
     867        y64 = (psVector *) y;
     868    } else if (y->type.type == PS_TYPE_F32) {
     869        y64Static = psVectorRecycle(y64Static, y->n, PS_TYPE_F64);
     870        p_psMemSetPersistent(y64Static, true);
     871        p_psMemSetPersistent(y64Static->data.V, true);
     872        y64 = y64Static;
     873    } else {
     874        // XXX: psError() bad type.
     875        psAbort(__func__, "Bad type for y64 (%d)", y->type.type);
     876    }
     877
    894878
    895879    // If yErr==NULL, set all errors equal.
    896880    if (yErr == NULL) {
    897         myYErr = psVectorAlloc(y->n, y->type.type);
    898         mustFreeMyYErr = true;
    899 
    900         if (y->type.type == PS_TYPE_F32) {
    901             for (i=0;i<myYErr->n;i++) {
    902                 myYErr->data.F32[i] = 1.0;
    903             }
    904         } else if (y->type.type == PS_TYPE_F64) {
    905             for (i=0;i<myYErr->n;i++) {
    906                 myYErr->data.F64[i] = 1.0;
    907             }
    908         }
     881        yErr64Static = psVectorRecycle(yErr64Static, y->n, PS_TYPE_F64);
     882        p_psMemSetPersistent(yErr64Static, true);
     883        p_psMemSetPersistent(yErr64Static->data.V, true);
     884
     885        for (i=0;i<yErr64Static->n;i++) {
     886            yErr64Static->data.F64[i] = 1.0;
     887        }
     888        yErr64 = yErr64Static;
    909889    } else {
    910         myYErr = (psVector *) yErr;
    911     }
    912 
    913     // If x==NULL, create an myX vector with x values set to (0:n), and if
     890        if (yErr->type.type == PS_TYPE_F64) {
     891            yErr64 = (psVector *) yErr;
     892        } else if (yErr->type.type == PS_TYPE_F32) {
     893            yErr64Static = psVectorRecycle(yErr64Static, yErr->n, PS_TYPE_F64);
     894            p_psMemSetPersistent(yErr64Static, true);
     895            p_psMemSetPersistent(yErr64Static->data.V, true);
     896        } else {
     897            // XXX: psError() bad type.
     898            psAbort(__func__, "Bad type for yErr64");
     899        }
     900    }
     901
     902    // If x==NULL, create an x64 vector with x values set to (0:n), and if
    914903    // this is a Chebyshev polynomial, we must scale to (-1:1).
    915 
    916904    // XXX: Verify that this is the correct action.
    917905    if (x == NULL) {
    918         myX = psVectorAlloc(y->n, y->type.type);
    919         mustFreemyX = true;
    920 
    921         if (y->type.type == PS_TYPE_F32) {
    922             if (myPoly->type == PS_POLYNOMIAL_ORD) {
    923                 for (i=0;i<myX->n;i++) {
    924                     myX->data.F32[i] = (float) i;
    925                 }
    926             } else if (myPoly->type == PS_POLYNOMIAL_CHEB) {
    927                 float min = 0.0;
    928                 float max = (float) (y->n - 1);
    929 
    930                 for (i=0;i<myX->n;i++) {
    931                     myX->data.F32[i] = (((float) i) - 0.5 * (min + max)) /
    932                                        (0.5 * (max - min));
    933                 }
    934             }
    935         } else if (y->type.type == PS_TYPE_F64) {
    936             if (myPoly->type == PS_POLYNOMIAL_ORD) {
    937                 for (i=0;i<myX->n;i++) {
    938                     myX->data.F64[i] = (float) i;
    939                 }
    940             } else if (myPoly->type == PS_POLYNOMIAL_CHEB) {
    941                 double min = 0.0;
    942                 double max = (double) (y->n - 1);
    943 
    944                 for (i=0;i<myX->n;i++) {
    945                     myX->data.F64[i] = (((float) i) - 0.5 * (min + max)) /
    946                                        (0.5 * (max - min));
    947                 }
     906        x64Static = psVectorRecycle(x64Static, y->n, PS_TYPE_F64);
     907        p_psMemSetPersistent(x64Static, true);
     908        p_psMemSetPersistent(x64Static->data.V, true);
     909        x64 = x64Static;
     910
     911        if (myPoly->type == PS_POLYNOMIAL_ORD) {
     912            for (i=0;i<x64->n;i++) {
     913                x64->data.F64[i] = (float) i;
     914            }
     915        } else if (myPoly->type == PS_POLYNOMIAL_CHEB) {
     916            double min = 0.0;
     917            double max = (double) (y->n - 1);
     918
     919            for (i=0;i<x64->n;i++) {
     920                x64->data.F64[i] = (((double) i) - 0.5 * (min + max)) /
     921                                   (0.5 * (max - min));
    948922            }
    949923        }
    950924    } else {
    951         myX = (psVector *) x;
    952     }
    953 
    954     PS_CHECK_VECTOR_SIZE_EQUAL(y, myX);
    955     PS_CHECK_VECTOR_SIZE_EQUAL(y, myYErr);
     925        if (x->type.type == PS_TYPE_F64) {
     926            x64 = (psVector *) x;
     927        } else if (x->type.type == PS_TYPE_F32) {
     928            x64Static = psVectorRecycle(x64Static, x->n, PS_TYPE_F64);
     929            p_psMemSetPersistent(x64Static, true);
     930            p_psMemSetPersistent(x64Static->data.V, true);
     931            x64 = x64Static;
     932        } else {
     933            // XXX: psError() bad type.
     934            psAbort(__func__, "Bad type for x64");
     935        }
     936    }
     937
     938    PS_CHECK_VECTOR_SIZE_EQUAL(y64, x64);
     939    PS_CHECK_VECTOR_SIZE_EQUAL(y64, yErr64);
    956940
    957941    // Call the appropriate vector fitting routine.
    958942    if (myPoly->type == PS_POLYNOMIAL_CHEB) {
    959         tmpPoly = p_psVectorFitPolynomial1DCheby(myPoly, myX, y, myYErr);
     943        tmpPoly = p_psVectorFitPolynomial1DCheby(myPoly, x64, y64, yErr64);
    960944    } else if (myPoly->type == PS_POLYNOMIAL_ORD) {
    961         tmpPoly = p_psVectorFitPolynomial1DOrd(myPoly, myX, y, myYErr);
     945        tmpPoly = p_psVectorFitPolynomial1DOrd(myPoly, x64, y64, yErr64);
    962946    } else {
    963947        // XXX: psErrorMsg()
    … …  
    966950
    967951    // Free any allocated memory.
    968     if (mustFreeMyYErr == true) {
    969         psFree(myYErr);
    970     }
    971     if (mustFreemyX == true) {
    972         psFree(myX);
    973     }
     952    /*
     953        if ((yErr == NULL) || (yErr->type.type != PS_TYPE_F64)) {
     954            psFree(yErr64);
     955        }
     956        if ((x == NULL) || (x->type.type != PS_TYPE_F64)) {
     957            psFree(x64);
     958        }
     959        if ((y == NULL) || (y->type.type != PS_TYPE_F64)) {
     960            psFree(y64);
     961        }
     962    */
    974963
    975964    return(myPoly);

  • trunk/psLib/src/math/psMinimize.h

    r1846 r1953  
    88 *  @author George Gusciora, MHPCC
    99 *
    10  *  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
    11  *  @date $Date: 2004-09-22 01:30:21 $
     10 *  @version $Revision: 1.26 $ $Name: not supported by cvs2svn $
     11 *  @date $Date: 2004-10-05 01:43:58 $
    1212 *
    1313 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    1818#define PS_MINIMIZE_H
    1919
     20#include <stdio.h>
     21#include <float.h>
     22#include <math.h>
     23
     24#include "psMemory.h"
     25#include "psVector.h"
     26#include "psArray.h"
     27#include "psImage.h"
     28#include "psMatrix.h"
    2029#include "psFunctions.h"
    21 #include "psArray.h"
    22 #include "psVector.h"
    23 #include "psMatrix.h"
     30#include "psTrace.h"
     31#include "psLogMsg.h"
     32#include "psError.h"
     33#include "psAbort.h"
     34#include "psConstants.h"
     35
    2436
    2537/** \file psMinimize.h
    … …  
    4456                                    float tol);
    4557
    46 /** Derive a polynomial fit by chi^2 minimisation (analytically)
     58/** Derive a polynomial fit.
    4759 *
    4860 *  psVectorFitPolynomial1d returns the polynomial that best fits the

  • trunk/psLib/src/math/psPolynomial.c

    r1947 r1953  
    77 *  polynomials.  It also contains a Gaussian functions.
    88 *
    9  *  @version $Revision: 1.45 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2004-10-04 00:40:23 $
     9 *  @version $Revision: 1.46 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2004-10-05 01:43:58 $
    1111 *
    1212 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    19881988 
    19891989XXX: nobody asked for us to generate this routine.
     1990 
     1991XXX: add another psScalar argument for the result.
     1992 
     1993XXX: The VectorCopy routines seg fault when I declare range32 as static.
    19901994 *****************************************************************************/
    19911995psScalar *p_psVectorInterpolate(psVector *domain,
    … …  
    19941998                                psScalar *x)
    19951999{
     2000    psVector *range32 = NULL;
     2001    psVector *domain32 = NULL;
    19962002    psTrace(".psLib.dataManip.psFunctions.p_psVectorInterpolate", 4,
    19972003            "---- p_psVectorInterpolate() begin ----\n");
    … …  
    20002006            (domain->type.type != x->type.type)) {
    20012007        // XXX psError
    2002         printf("domain->type.type != range->type.type != x->type.type\n");
     2008        printf("p_psVectorInterpolate(): domain->type.type != range->type.type != x->type.type\n");
    20032009        printf("%d %d %d\n", domain->type.type, range->type.type, x->type.type);
     2010        exit(1);
    20042011    }
    20052012    if (domain->n != range->n) {
    20062013        // XXX psError
    20072014        printf("domain->n != range->n\n");
     2015        exit(1);
    20082016    }
    20092017    if (order > (domain->n - 1)) {
    20102018        // XXX psError: not enough data points for order-order interpolation.
    20112019        printf("not enough data points for order-order interpolation.\n");
     2020        exit(1);
    20122021    }
    20132022
    … …  
    20202029                             order,
    20212030                             x->data.F32), PS_TYPE_F32));
     2031    } else if (x->type.type == PS_TYPE_F64) {
     2032        // XXX: use recycled vectors here.
     2033        range32 = psVectorCopy(range32, range, PS_TYPE_F32);
     2034        domain32 = psVectorCopy(domain32, domain, PS_TYPE_F32);
     2035
     2036        psScalar *tmpScalar = psScalarAlloc((double)
     2037                                            p_ps1DInterpolateF32(domain32->data.F32,
     2038                                                                 range32->data.F32,
     2039                                                                 domain32->n,
     2040                                                                 order,
     2041                                                                 (float) x->data.F64), PS_TYPE_F64);
     2042        psFree(range32);
     2043        psFree(domain32);
     2044
     2045        psTrace(".psLib.dataManip.psFunctions.p_psVectorInterpolate", 4,
     2046                "---- p_psVectorInterpolate() end ----\n");
     2047        // XXX: Convert data type to F64?
     2048        return(tmpScalar);
     2049
    20222050    } else {
    20232051        // XXX psError: type not supported

  • trunk/psLib/src/math/psSpline.c

    r1947 r1953  
    77 *  polynomials.  It also contains a Gaussian functions.
    88 *
    9  *  @version $Revision: 1.45 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2004-10-04 00:40:23 $
     9 *  @version $Revision: 1.46 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2004-10-05 01:43:58 $
    1111 *
    1212 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    19881988 
    19891989XXX: nobody asked for us to generate this routine.
     1990 
     1991XXX: add another psScalar argument for the result.
     1992 
     1993XXX: The VectorCopy routines seg fault when I declare range32 as static.
    19901994 *****************************************************************************/
    19911995psScalar *p_psVectorInterpolate(psVector *domain,
    … …  
    19941998                                psScalar *x)
    19951999{
     2000    psVector *range32 = NULL;
     2001    psVector *domain32 = NULL;
    19962002    psTrace(".psLib.dataManip.psFunctions.p_psVectorInterpolate", 4,
    19972003            "---- p_psVectorInterpolate() begin ----\n");
    … …  
    20002006            (domain->type.type != x->type.type)) {
    20012007        // XXX psError
    2002         printf("domain->type.type != range->type.type != x->type.type\n");
     2008        printf("p_psVectorInterpolate(): domain->type.type != range->type.type != x->type.type\n");
    20032009        printf("%d %d %d\n", domain->type.type, range->type.type, x->type.type);
     2010        exit(1);
    20042011    }
    20052012    if (domain->n != range->n) {
    20062013        // XXX psError
    20072014        printf("domain->n != range->n\n");
     2015        exit(1);
    20082016    }
    20092017    if (order > (domain->n - 1)) {
    20102018        // XXX psError: not enough data points for order-order interpolation.
    20112019        printf("not enough data points for order-order interpolation.\n");
     2020        exit(1);
    20122021    }
    20132022
    … …  
    20202029                             order,
    20212030                             x->data.F32), PS_TYPE_F32));
     2031    } else if (x->type.type == PS_TYPE_F64) {
     2032        // XXX: use recycled vectors here.
     2033        range32 = psVectorCopy(range32, range, PS_TYPE_F32);
     2034        domain32 = psVectorCopy(domain32, domain, PS_TYPE_F32);
     2035
     2036        psScalar *tmpScalar = psScalarAlloc((double)
     2037                                            p_ps1DInterpolateF32(domain32->data.F32,
     2038                                                                 range32->data.F32,
     2039                                                                 domain32->n,
     2040                                                                 order,
     2041                                                                 (float) x->data.F64), PS_TYPE_F64);
     2042        psFree(range32);
     2043        psFree(domain32);
     2044
     2045        psTrace(".psLib.dataManip.psFunctions.p_psVectorInterpolate", 4,
     2046                "---- p_psVectorInterpolate() end ----\n");
     2047        // XXX: Convert data type to F64?
     2048        return(tmpScalar);
     2049
    20222050    } else {
    20232051        // XXX psError: type not supported

  • trunk/psLib/src/mathtypes/psVector.h

    r1920 r1953  
    1 
    21/** @file  psVector.h
    32 *
    … …  
    1211 *  @author Ross Harman, MHPCC
    1312 *
    14  *  @version $Revision: 1.22 $ $Name: not supported by cvs2svn $
    15  *  @date $Date: 2004-09-28 23:26:48 $
     13 *  @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
     14 *  @date $Date: 2004-10-05 01:43:57 $
    1615 *
    1716 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
Note: See TracChangeset for help on using the changeset viewer.