IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 1008


Ignore:
Timestamp:
Jun 11, 2004, 3:33:16 PM (22 years ago)
Author:
desonia
Message:

bug fixes found while testing psFFT routines.

Location:
trunk/psLib/src
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/Makefile.Globals

    r806 r1008  
    55##  Assumptions:    Variable "prefix" already defined
    66##
    7 ##  $Revision: 1.9 $  $Name: not supported by cvs2svn $
    8 ##  $Date: 2004-05-28 20:26:49 $
     7##  $Revision: 1.10 $  $Name: not supported by cvs2svn $
     8##  $Date: 2004-06-12 01:33:16 $
    99##
    1010##  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    2626ifndef includedir
    2727    export includedir = $(prefix)/include
     28endif
     29
     30ifndef testbindir
     31    export testbindir = $(prefix)/bin/test
    2832endif
    2933

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

    r945 r1008  
    88 *  @author Ross Harman, MHPCC
    99 *
    10  *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
    11  *  @date $Date: 2004-06-09 18:20:01 $
     10 *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
     11 *  @date $Date: 2004-06-12 01:33:16 $
    1212 *
    1313 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    145145    }
    146146
    147     in->data.V = psRealloc(in->data.V,nalloc*PSELEMTYPE_SIZEOF(elemType));
     147    in->data.V = psRealloc(in->data.V,nalloc*PSELEMTYPE_SIZEOF(type));
    148148
    149149    in->n = 0;
    150     in->type.type = elemType;
     150    in->type.type = type;
    151151    in->nalloc = nalloc;
    152152

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

    r831 r1008  
    55 *  @author Robert DeSonia, MHPCC
    66 *
    7  *  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2004-06-02 23:29:21 $
     7 *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2004-06-12 01:33:16 $
    99 *
    1010 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    448448    numElements = in->n;
    449449
    450     out = psVectorRecycle(out, numElements, PS_TYPE_C32);
     450    out = psVectorRecycle(out, PS_TYPE_C32, numElements);
     451    out->n = numElements;
    451452
    452453    if (type == PS_TYPE_F32) {
    … …  
    504505        psLogMsg(__func__,PS_LOG_WARN,"Real portion of a non-Complex type called called for. "
    505506                 "Just a vector copy was performed.");
    506         out = psVectorRecycle(out,numElements,type);
     507        out = psVectorRecycle(out,type,numElements);
     508        out->n = numElements;
    507509        memcpy(out->data.V,in->data.V,numElements*PSELEMTYPE_SIZEOF(type));
    508510        return out;
    … …  
    513515        psC32* inVec = in->data.C32;
    514516
    515         out = psVectorRecycle(out,numElements,PS_TYPE_F32);
     517        out = psVectorRecycle(out,PS_TYPE_F32,numElements);
     518        out->n = numElements;
    516519        outVec = out->data.F32;
    517520
    … …  
    547550        psLogMsg(__func__,PS_LOG_WARN,"Imaginary portion of a non-Complex type called for. "
    548551                 "A zeroed vector was returned.");
    549         out = psVectorRecycle(out,numElements,type);
     552        out = psVectorRecycle(out,type,numElements);
     553        out->n = numElements;
    550554        memset(out->data.V,0,PSELEMTYPE_SIZEOF(type)*numElements);
    551555        return out;
    … …  
    556560        psC32* inVec = in->data.C32;
    557561
    558         out = psVectorRecycle(out,numElements, PS_TYPE_F32);
     562        out = psVectorRecycle(out,PS_TYPE_F32,numElements);
     563        out->n = numElements;
    559564        outVec = out->data.F32;
    560565
    … …  
    579584
    580585    if (real == NULL || imag == NULL) {
     586        psVectorFree(out);
    581587        return NULL;
    582588    }
    … …  
    591597    if (imag->type.type != type) {
    592598        psError(__func__,"The inputs to psVectorComplex must be the same type.");
     599        psVectorFree(out);
    593600        return NULL;
    594601    }
    … …  
    596603    if (PS_IS_PSELEMTYPE_COMPLEX(type)) {
    597604        psError(__func__,"The inputs to psVectorComplex can not be complex.");
     605        psVectorFree(out);
    598606        return NULL;
    599607    }
    … …  
    604612        psF32* imagVec = imag->data.F32;
    605613
    606         out = psVectorRecycle(out,numElements, PS_TYPE_C32);
     614        out = psVectorRecycle(out,PS_TYPE_C32,numElements);
     615        out->n = numElements;
    607616        outVec = out->data.C32;
    608617
    … …  
    639648                 "Vector copy was performed instead.");
    640649
    641         out = psVectorRecycle(out,numElements,type);
     650        out = psVectorRecycle(out,type,numElements);
     651        out->n = numElements;
    642652        memcpy(out->data.V,in->data.V,PSELEMTYPE_SIZEOF(type)*numElements);
    643653        return out;
    … …  
    648658        psC32* inVec = in->data.C32;
    649659
    650         out = psVectorRecycle(out,numElements,PS_TYPE_C32);
     660        out = psVectorRecycle(out,PS_TYPE_C32, numElements);
     661        out->n = numElements;
    651662        outVec = out->data.C32;
    652663
    … …  
    667678{
    668679    psElemType type;
    669     unsigned int numElements;
    670     int numElementsSquared;
     680    unsigned int outNumElements;
     681    unsigned int inNumElements;
     682    unsigned int inHalfNumElements;
     683    unsigned int inNumElementsSquared;
    671684
    672685    if (in == NULL) {
    673         return NULL;
    674     }
    675 
    676     type = in->type.type;
    677     numElements = in->n;
    678     numElementsSquared = numElements*numElements;
     686        psVectorFree(out);
     687        return NULL;
     688    }
     689
     690    type = in->type.type;
     691    inNumElements = in->n;
     692    inNumElementsSquared = inNumElements*inNumElements;
     693    inHalfNumElements = inNumElements/2;
     694    outNumElements = inHalfNumElements+1;
    679695
    680696    /* if not a complex number, this is not implemented */
    681697    if (! PS_IS_PSELEMTYPE_COMPLEX(type)) {
    682698        psError(__func__,"Power Spectrum for non-complex inputs is not implemented.");
     699        psVectorFree(out);
    683700        return NULL;
    684701    }
    … …  
    687704        psF32* outVec;
    688705        psC32* inVec = in->data.C32;
    689         psF32 real;
    690         psF32 imag;
    691 
    692 
    693         out = psVectorRecycle(out,numElements,PS_TYPE_F32);
     706        psF32 inAbs1;
     707        psF32 inAbs2;
     708
     709        out = psVectorRecycle(out,PS_TYPE_F32,outNumElements);
     710        out->n = outNumElements;
    694711        outVec = out->data.F32;
    695712
    696         for (unsigned int i=0;i<numElements;i++) {
    697             real = crealf(inVec[i]);
    698             imag = cimagf(inVec[i]);
    699             outVec[i] = (real*real+imag*imag)/numElementsSquared;
    700         }
     713        // from ADD: P_0 = |C_0|^2/N^2
     714        inAbs1 = cabsf(inVec[0]);
     715        outVec[0] = inAbs1*inAbs1/inNumElementsSquared;
     716
     717        // from ADD: P_j = (|C_j|^2+|C_N-j|^2)/N^2, where j = 1,2,...,(N/2-1)
     718        for (unsigned int i=1;i<inHalfNumElements;i++) {
     719            inAbs1 = cabsf(inVec[i]);
     720            inAbs2 = cabsf(inVec[inNumElements-i]);
     721            outVec[i] = (inAbs1*inAbs1+inAbs2*inAbs2)/inNumElementsSquared;
     722        }
     723
     724        // from ADD: P_N/2 = |C_N/2|^2/N^2
     725        inAbs1 = cabsf(inVec[inHalfNumElements]);
     726        outVec[inHalfNumElements] = inAbs1*inAbs1/inNumElementsSquared;
    701727    } else {
    702728        psError(__func__,"Can not power spectrum for given vector type (%d).",

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

    r831 r1008  
    55 *  @author Robert DeSonia, MHPCC
    66 *
    7  *  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2004-06-02 23:29:21 $
     7 *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2004-06-12 01:33:16 $
    99 *
    1010 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    448448    numElements = in->n;
    449449
    450     out = psVectorRecycle(out, numElements, PS_TYPE_C32);
     450    out = psVectorRecycle(out, PS_TYPE_C32, numElements);
     451    out->n = numElements;
    451452
    452453    if (type == PS_TYPE_F32) {
    … …  
    504505        psLogMsg(__func__,PS_LOG_WARN,"Real portion of a non-Complex type called called for. "
    505506                 "Just a vector copy was performed.");
    506         out = psVectorRecycle(out,numElements,type);
     507        out = psVectorRecycle(out,type,numElements);
     508        out->n = numElements;
    507509        memcpy(out->data.V,in->data.V,numElements*PSELEMTYPE_SIZEOF(type));
    508510        return out;
    … …  
    513515        psC32* inVec = in->data.C32;
    514516
    515         out = psVectorRecycle(out,numElements,PS_TYPE_F32);
     517        out = psVectorRecycle(out,PS_TYPE_F32,numElements);
     518        out->n = numElements;
    516519        outVec = out->data.F32;
    517520
    … …  
    547550        psLogMsg(__func__,PS_LOG_WARN,"Imaginary portion of a non-Complex type called for. "
    548551                 "A zeroed vector was returned.");
    549         out = psVectorRecycle(out,numElements,type);
     552        out = psVectorRecycle(out,type,numElements);
     553        out->n = numElements;
    550554        memset(out->data.V,0,PSELEMTYPE_SIZEOF(type)*numElements);
    551555        return out;
    … …  
    556560        psC32* inVec = in->data.C32;
    557561
    558         out = psVectorRecycle(out,numElements, PS_TYPE_F32);
     562        out = psVectorRecycle(out,PS_TYPE_F32,numElements);
     563        out->n = numElements;
    559564        outVec = out->data.F32;
    560565
    … …  
    579584
    580585    if (real == NULL || imag == NULL) {
     586        psVectorFree(out);
    581587        return NULL;
    582588    }
    … …  
    591597    if (imag->type.type != type) {
    592598        psError(__func__,"The inputs to psVectorComplex must be the same type.");
     599        psVectorFree(out);
    593600        return NULL;
    594601    }
    … …  
    596603    if (PS_IS_PSELEMTYPE_COMPLEX(type)) {
    597604        psError(__func__,"The inputs to psVectorComplex can not be complex.");
     605        psVectorFree(out);
    598606        return NULL;
    599607    }
    … …  
    604612        psF32* imagVec = imag->data.F32;
    605613
    606         out = psVectorRecycle(out,numElements, PS_TYPE_C32);
     614        out = psVectorRecycle(out,PS_TYPE_C32,numElements);
     615        out->n = numElements;
    607616        outVec = out->data.C32;
    608617
    … …  
    639648                 "Vector copy was performed instead.");
    640649
    641         out = psVectorRecycle(out,numElements,type);
     650        out = psVectorRecycle(out,type,numElements);
     651        out->n = numElements;
    642652        memcpy(out->data.V,in->data.V,PSELEMTYPE_SIZEOF(type)*numElements);
    643653        return out;
    … …  
    648658        psC32* inVec = in->data.C32;
    649659
    650         out = psVectorRecycle(out,numElements,PS_TYPE_C32);
     660        out = psVectorRecycle(out,PS_TYPE_C32, numElements);
     661        out->n = numElements;
    651662        outVec = out->data.C32;
    652663
    … …  
    667678{
    668679    psElemType type;
    669     unsigned int numElements;
    670     int numElementsSquared;
     680    unsigned int outNumElements;
     681    unsigned int inNumElements;
     682    unsigned int inHalfNumElements;
     683    unsigned int inNumElementsSquared;
    671684
    672685    if (in == NULL) {
    673         return NULL;
    674     }
    675 
    676     type = in->type.type;
    677     numElements = in->n;
    678     numElementsSquared = numElements*numElements;
     686        psVectorFree(out);
     687        return NULL;
     688    }
     689
     690    type = in->type.type;
     691    inNumElements = in->n;
     692    inNumElementsSquared = inNumElements*inNumElements;
     693    inHalfNumElements = inNumElements/2;
     694    outNumElements = inHalfNumElements+1;
    679695
    680696    /* if not a complex number, this is not implemented */
    681697    if (! PS_IS_PSELEMTYPE_COMPLEX(type)) {
    682698        psError(__func__,"Power Spectrum for non-complex inputs is not implemented.");
     699        psVectorFree(out);
    683700        return NULL;
    684701    }
    … …  
    687704        psF32* outVec;
    688705        psC32* inVec = in->data.C32;
    689         psF32 real;
    690         psF32 imag;
    691 
    692 
    693         out = psVectorRecycle(out,numElements,PS_TYPE_F32);
     706        psF32 inAbs1;
     707        psF32 inAbs2;
     708
     709        out = psVectorRecycle(out,PS_TYPE_F32,outNumElements);
     710        out->n = outNumElements;
    694711        outVec = out->data.F32;
    695712
    696         for (unsigned int i=0;i<numElements;i++) {
    697             real = crealf(inVec[i]);
    698             imag = cimagf(inVec[i]);
    699             outVec[i] = (real*real+imag*imag)/numElementsSquared;
    700         }
     713        // from ADD: P_0 = |C_0|^2/N^2
     714        inAbs1 = cabsf(inVec[0]);
     715        outVec[0] = inAbs1*inAbs1/inNumElementsSquared;
     716
     717        // from ADD: P_j = (|C_j|^2+|C_N-j|^2)/N^2, where j = 1,2,...,(N/2-1)
     718        for (unsigned int i=1;i<inHalfNumElements;i++) {
     719            inAbs1 = cabsf(inVec[i]);
     720            inAbs2 = cabsf(inVec[inNumElements-i]);
     721            outVec[i] = (inAbs1*inAbs1+inAbs2*inAbs2)/inNumElementsSquared;
     722        }
     723
     724        // from ADD: P_N/2 = |C_N/2|^2/N^2
     725        inAbs1 = cabsf(inVec[inHalfNumElements]);
     726        outVec[inHalfNumElements] = inAbs1*inAbs1/inNumElementsSquared;
    701727    } else {
    702728        psError(__func__,"Can not power spectrum for given vector type (%d).",

  • trunk/psLib/src/fft/psVectorFFT.c

    r831 r1008  
    55 *  @author Robert DeSonia, MHPCC
    66 *
    7  *  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2004-06-02 23:29:21 $
     7 *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2004-06-12 01:33:16 $
    99 *
    1010 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    448448    numElements = in->n;
    449449
    450     out = psVectorRecycle(out, numElements, PS_TYPE_C32);
     450    out = psVectorRecycle(out, PS_TYPE_C32, numElements);
     451    out->n = numElements;
    451452
    452453    if (type == PS_TYPE_F32) {
    … …  
    504505        psLogMsg(__func__,PS_LOG_WARN,"Real portion of a non-Complex type called called for. "
    505506                 "Just a vector copy was performed.");
    506         out = psVectorRecycle(out,numElements,type);
     507        out = psVectorRecycle(out,type,numElements);
     508        out->n = numElements;
    507509        memcpy(out->data.V,in->data.V,numElements*PSELEMTYPE_SIZEOF(type));
    508510        return out;
    … …  
    513515        psC32* inVec = in->data.C32;
    514516
    515         out = psVectorRecycle(out,numElements,PS_TYPE_F32);
     517        out = psVectorRecycle(out,PS_TYPE_F32,numElements);
     518        out->n = numElements;
    516519        outVec = out->data.F32;
    517520
    … …  
    547550        psLogMsg(__func__,PS_LOG_WARN,"Imaginary portion of a non-Complex type called for. "
    548551                 "A zeroed vector was returned.");
    549         out = psVectorRecycle(out,numElements,type);
     552        out = psVectorRecycle(out,type,numElements);
     553        out->n = numElements;
    550554        memset(out->data.V,0,PSELEMTYPE_SIZEOF(type)*numElements);
    551555        return out;
    … …  
    556560        psC32* inVec = in->data.C32;
    557561
    558         out = psVectorRecycle(out,numElements, PS_TYPE_F32);
     562        out = psVectorRecycle(out,PS_TYPE_F32,numElements);
     563        out->n = numElements;
    559564        outVec = out->data.F32;
    560565
    … …  
    579584
    580585    if (real == NULL || imag == NULL) {
     586        psVectorFree(out);
    581587        return NULL;
    582588    }
    … …  
    591597    if (imag->type.type != type) {
    592598        psError(__func__,"The inputs to psVectorComplex must be the same type.");
     599        psVectorFree(out);
    593600        return NULL;
    594601    }
    … …  
    596603    if (PS_IS_PSELEMTYPE_COMPLEX(type)) {
    597604        psError(__func__,"The inputs to psVectorComplex can not be complex.");
     605        psVectorFree(out);
    598606        return NULL;
    599607    }
    … …  
    604612        psF32* imagVec = imag->data.F32;
    605613
    606         out = psVectorRecycle(out,numElements, PS_TYPE_C32);
     614        out = psVectorRecycle(out,PS_TYPE_C32,numElements);
     615        out->n = numElements;
    607616        outVec = out->data.C32;
    608617
    … …  
    639648                 "Vector copy was performed instead.");
    640649
    641         out = psVectorRecycle(out,numElements,type);
     650        out = psVectorRecycle(out,type,numElements);
     651        out->n = numElements;
    642652        memcpy(out->data.V,in->data.V,PSELEMTYPE_SIZEOF(type)*numElements);
    643653        return out;
    … …  
    648658        psC32* inVec = in->data.C32;
    649659
    650         out = psVectorRecycle(out,numElements,PS_TYPE_C32);
     660        out = psVectorRecycle(out,PS_TYPE_C32, numElements);
     661        out->n = numElements;
    651662        outVec = out->data.C32;
    652663
    … …  
    667678{
    668679    psElemType type;
    669     unsigned int numElements;
    670     int numElementsSquared;
     680    unsigned int outNumElements;
     681    unsigned int inNumElements;
     682    unsigned int inHalfNumElements;
     683    unsigned int inNumElementsSquared;
    671684
    672685    if (in == NULL) {
    673         return NULL;
    674     }
    675 
    676     type = in->type.type;
    677     numElements = in->n;
    678     numElementsSquared = numElements*numElements;
     686        psVectorFree(out);
     687        return NULL;
     688    }
     689
     690    type = in->type.type;
     691    inNumElements = in->n;
     692    inNumElementsSquared = inNumElements*inNumElements;
     693    inHalfNumElements = inNumElements/2;
     694    outNumElements = inHalfNumElements+1;
    679695
    680696    /* if not a complex number, this is not implemented */
    681697    if (! PS_IS_PSELEMTYPE_COMPLEX(type)) {
    682698        psError(__func__,"Power Spectrum for non-complex inputs is not implemented.");
     699        psVectorFree(out);
    683700        return NULL;
    684701    }
    … …  
    687704        psF32* outVec;
    688705        psC32* inVec = in->data.C32;
    689         psF32 real;
    690         psF32 imag;
    691 
    692 
    693         out = psVectorRecycle(out,numElements,PS_TYPE_F32);
     706        psF32 inAbs1;
     707        psF32 inAbs2;
     708
     709        out = psVectorRecycle(out,PS_TYPE_F32,outNumElements);
     710        out->n = outNumElements;
    694711        outVec = out->data.F32;
    695712
    696         for (unsigned int i=0;i<numElements;i++) {
    697             real = crealf(inVec[i]);
    698             imag = cimagf(inVec[i]);
    699             outVec[i] = (real*real+imag*imag)/numElementsSquared;
    700         }
     713        // from ADD: P_0 = |C_0|^2/N^2
     714        inAbs1 = cabsf(inVec[0]);
     715        outVec[0] = inAbs1*inAbs1/inNumElementsSquared;
     716
     717        // from ADD: P_j = (|C_j|^2+|C_N-j|^2)/N^2, where j = 1,2,...,(N/2-1)
     718        for (unsigned int i=1;i<inHalfNumElements;i++) {
     719            inAbs1 = cabsf(inVec[i]);
     720            inAbs2 = cabsf(inVec[inNumElements-i]);
     721            outVec[i] = (inAbs1*inAbs1+inAbs2*inAbs2)/inNumElementsSquared;
     722        }
     723
     724        // from ADD: P_N/2 = |C_N/2|^2/N^2
     725        inAbs1 = cabsf(inVec[inHalfNumElements]);
     726        outVec[inHalfNumElements] = inAbs1*inAbs1/inNumElementsSquared;
    701727    } else {
    702728        psError(__func__,"Can not power spectrum for given vector type (%d).",

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

    r945 r1008  
    88 *  @author Ross Harman, MHPCC
    99 *
    10  *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
    11  *  @date $Date: 2004-06-09 18:20:01 $
     10 *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
     11 *  @date $Date: 2004-06-12 01:33:16 $
    1212 *
    1313 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    145145    }
    146146
    147     in->data.V = psRealloc(in->data.V,nalloc*PSELEMTYPE_SIZEOF(elemType));
     147    in->data.V = psRealloc(in->data.V,nalloc*PSELEMTYPE_SIZEOF(type));
    148148
    149149    in->n = 0;
    150     in->type.type = elemType;
     150    in->type.type = type;
    151151    in->nalloc = nalloc;
    152152
Note: See TracChangeset for help on using the changeset viewer.