Changeset 3026

Timestamp:

Jan 17, 2005, 12:17:29 PM (22 years ago)

Author:

desonia

Message:

changed psMatrixLUD to create the perm parameter according to the machine
size of 'size_t', which is used by GSL.

Location:

trunk/psLib

Files:

• src/dataManip/psFunctions.c (modified) (2 diffs)
• src/dataManip/psMatrix.c (modified) (5 diffs)
• src/dataManip/psMatrix.h (modified) (2 diffs)
• src/dataManip/psMinimize.c (modified) (9 diffs)
• src/math/psMatrix.c (modified) (5 diffs)
• src/math/psMatrix.h (modified) (2 diffs)
• src/math/psMinimize.c (modified) (9 diffs)
• src/math/psPolynomial.c (modified) (2 diffs)
• src/math/psSpline.c (modified) (2 diffs)
• test/dataIO (modified) (1 prop)
• test/dataIO/.cvsignore (modified) (1 diff)
• test/dataManip/tst_psMatrix03.c (modified) (4 diffs)
• test/fileUtils (modified) (1 prop)
• test/fileUtils/.cvsignore (modified) (1 diff)

trunk/psLib/src/dataManip/psFunctions.c

@@ -7 +7 @@
  *  polynomials.  It also contains a Gaussian functions.
  *
- *  @version $Revision: 1.81 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-01-14 23:40:24 $
+ *  @version $Revision: 1.82 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-01-17 22:17:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -1906 +1906 @@
         (tmp->spline)[i] = psPolynomial1DAlloc(order+1, PS_POLYNOMIAL_ORD);
     }
+
+    // This should be set by the psVectorFitSpline1D()
+    tmp->p_psDeriv2 = NULL;
 
     tmp->domains = (psF32 *) psAlloc((bounds->n) * sizeof(psF32));

trunk/psLib/src/dataManip/psMatrix.c

@@ -20 +20 @@
  *  @author Ross Harman, MHPCC
  *
- *  @version $Revision: 1.21 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-12-11 00:16:38 $
+ *  @version $Revision: 1.22 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-01-17 22:17:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -175 +175 @@
 /*****************************************************************************/
 
-psImage* psMatrixLUD(psImage* outImage, psVector* outPerm, psImage* inImage)
+psImage* psMatrixLUD(psImage* outImage, psVector** outPerm, psImage* inImage)
 {
     psS32 signum = 0;
@@ -190 +190 @@
     PS_CHECK_POINTERS(inImage, outImage, outImage);
     PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, outImage);
-    PS_VECTOR_CHECK_NULL_GENERAL(outPerm, psMatrixLUD_EXIT);
-    PS_CHECK_DIMEN_AND_TYPE(outPerm, PS_DIMEN_VECTOR, outImage);
+    PS_PTR_CHECK_NULL_GENERAL(outPerm, psMatrixLUD_EXIT);
     psImageRecycle(outImage, inImage->numCols, inImage->numRows, inImage->type.type);
-    psVectorRecycle(outPerm, inImage->numRows, inImage->type.type);
     PS_CHECK_SQUARE(inImage, outImage);
     PS_CHECK_SQUARE(outImage, outImage);
@@ -203 +201 @@
     // Initialize GSL data
     perm.size = numCols;
-    outPerm->n = numCols;
-    perm.data = outPerm->data.V;
+    if (sizeof(size_t) == 4) {
+        *outPerm = psVectorRecycle(*outPerm, numCols, PS_TYPE_S32);
+    } else if (sizeof(size_t) == 8) {
+        *outPerm = psVectorRecycle(*outPerm, numCols, PS_TYPE_S64);
+    } else {
+        psError(PS_ERR_UNKNOWN, true,
+                "Failed to allocate the permutation vector; "
+                "could not determine the cooresponding data type.");
+        psMatrixLUD_EXIT;
+    }
+
+    (*outPerm)->n = numCols;
+    perm.data = (*outPerm)->data.V;
     lu = gsl_matrix_alloc(numRows, numCols);
 
@@ -244 +253 @@
     PS_CHECK_DIMEN_AND_TYPE(inVector, PS_DIMEN_VECTOR, outVector);
     PS_VECTOR_CHECK_NULL(inPerm, outVector);
-    PS_CHECK_DIMEN_AND_TYPE(inPerm, PS_DIMEN_VECTOR, outVector);
     psVectorRecycle(outVector, inImage->numRows, inImage->type.type);
 

trunk/psLib/src/dataManip/psMatrix.h

@@ -22 +22 @@
  *  @author Ross Harman, MHPCC
  *
- *  @version $Revision: 1.13 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-12-09 20:51:22 $
+ *  @version $Revision: 1.14 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-01-17 22:17:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -45 +45 @@
 psImage* psMatrixLUD(
     psImage* outImage,                 ///< Image to return, or NULL.
-    psVector* outPerm,                 ///< Output permutation vector used by psMatrixLUSolve.
+    psVector** outPerm,                ///< Output permutation vector used by psMatrixLUSolve.
     psImage* inImage                   ///< Image to decompose.
 );

trunk/psLib/src/dataManip/psMinimize.c

@@ -9 +9 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.99 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-01-14 23:27:56 $
+ *  @version $Revision: 1.100 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-01-17 22:17:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -294 +294 @@
         This can not be implemented until SDR states what order spline should be
         created.
-     
+
         Should we error if mySPline is not NULL?
     */
@@ -616 +616 @@
     psS32 p;
     psVector *beta = psVectorAlloc(numParams, PS_TYPE_F64);
-    psVector *perm = psVectorAlloc(numParams, PS_TYPE_F64);
+    psVector *perm = NULL;
 
     psVector *paramDeltasF64 = psVectorAlloc(numParams, PS_TYPE_F64);
@@ -738 +738 @@
         // XXX: How do we know if these functions were successful?
         //
-        aOut = psMatrixLUD(aOut, perm, A);
+        aOut = psMatrixLUD(aOut, &perm, A);
         paramDeltasF64 = psMatrixLUSolve(paramDeltasF64, aOut, beta, perm);
 
@@ -936 +936 @@
     coeffs = psVectorAlloc(polyOrder, PS_TYPE_F64);
     X = psVectorAlloc(x->n, PS_TYPE_F64);
-    outPerm = psVectorAlloc(polyOrder, PS_TYPE_F64);
     xSums = psVectorAlloc(1 + 2 * polyOrder, PS_TYPE_F64);
 
@@ -943 +942 @@
         B->data.F64[i] = 0.0;
         coeffs->data.F64[i] = 0.0;
-        outPerm->data.F64[i] = 0.0;
         for (j = 0; j < polyOrder; j++) {
             A->data.F64[i][j] = 0.0;
@@ -987 +985 @@
 
     // XXX: How do we know if these routines were successful?
-    ALUD = psMatrixLUD(ALUD, outPerm, A);
+    ALUD = psMatrixLUD(ALUD, &outPerm, A);
     coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
 
@@ -1139 +1137 @@
 Algorithm:
  
-XXX completely ad hoc:  
+XXX completely ad hoc:
 start with the user-supplied starting parameter and
 call that b.  Calculate a/c as a fractional amount smaller/larger than b.
 Repeat this process until a local minimum is found.
  
-XXX:  
-new algorithm:  
+XXX:
+new algorithm:
 start at x=0, expand in one direction until the function
 decreases.  Then you have two points in the bracket.  Keep going until it
@@ -1151 +1149 @@
 direction.
  
-XXX: 
+XXX:
 This is F32 only.
  
-XXX: 
+XXX:
 output bracket vector should be an input as well.
 *****************************************************************************/

trunk/psLib/src/math/psMatrix.c

@@ -20 +20 @@
  *  @author Ross Harman, MHPCC
  *
- *  @version $Revision: 1.21 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-12-11 00:16:38 $
+ *  @version $Revision: 1.22 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-01-17 22:17:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -175 +175 @@
 /*****************************************************************************/
 
-psImage* psMatrixLUD(psImage* outImage, psVector* outPerm, psImage* inImage)
+psImage* psMatrixLUD(psImage* outImage, psVector** outPerm, psImage* inImage)
 {
     psS32 signum = 0;
@@ -190 +190 @@
     PS_CHECK_POINTERS(inImage, outImage, outImage);
     PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, outImage);
-    PS_VECTOR_CHECK_NULL_GENERAL(outPerm, psMatrixLUD_EXIT);
-    PS_CHECK_DIMEN_AND_TYPE(outPerm, PS_DIMEN_VECTOR, outImage);
+    PS_PTR_CHECK_NULL_GENERAL(outPerm, psMatrixLUD_EXIT);
     psImageRecycle(outImage, inImage->numCols, inImage->numRows, inImage->type.type);
-    psVectorRecycle(outPerm, inImage->numRows, inImage->type.type);
     PS_CHECK_SQUARE(inImage, outImage);
     PS_CHECK_SQUARE(outImage, outImage);
@@ -203 +201 @@
     // Initialize GSL data
     perm.size = numCols;
-    outPerm->n = numCols;
-    perm.data = outPerm->data.V;
+    if (sizeof(size_t) == 4) {
+        *outPerm = psVectorRecycle(*outPerm, numCols, PS_TYPE_S32);
+    } else if (sizeof(size_t) == 8) {
+        *outPerm = psVectorRecycle(*outPerm, numCols, PS_TYPE_S64);
+    } else {
+        psError(PS_ERR_UNKNOWN, true,
+                "Failed to allocate the permutation vector; "
+                "could not determine the cooresponding data type.");
+        psMatrixLUD_EXIT;
+    }
+
+    (*outPerm)->n = numCols;
+    perm.data = (*outPerm)->data.V;
     lu = gsl_matrix_alloc(numRows, numCols);
 
@@ -244 +253 @@
     PS_CHECK_DIMEN_AND_TYPE(inVector, PS_DIMEN_VECTOR, outVector);
     PS_VECTOR_CHECK_NULL(inPerm, outVector);
-    PS_CHECK_DIMEN_AND_TYPE(inPerm, PS_DIMEN_VECTOR, outVector);
     psVectorRecycle(outVector, inImage->numRows, inImage->type.type);
 

trunk/psLib/src/math/psMatrix.h

@@ -22 +22 @@
  *  @author Ross Harman, MHPCC
  *
- *  @version $Revision: 1.13 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-12-09 20:51:22 $
+ *  @version $Revision: 1.14 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-01-17 22:17:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -45 +45 @@
 psImage* psMatrixLUD(
     psImage* outImage,                 ///< Image to return, or NULL.
-    psVector* outPerm,                 ///< Output permutation vector used by psMatrixLUSolve.
+    psVector** outPerm,                ///< Output permutation vector used by psMatrixLUSolve.
     psImage* inImage                   ///< Image to decompose.
 );

trunk/psLib/src/math/psMinimize.c

@@ -9 +9 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.99 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-01-14 23:27:56 $
+ *  @version $Revision: 1.100 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-01-17 22:17:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -294 +294 @@
         This can not be implemented until SDR states what order spline should be
         created.
-     
+
         Should we error if mySPline is not NULL?
     */
@@ -616 +616 @@
     psS32 p;
     psVector *beta = psVectorAlloc(numParams, PS_TYPE_F64);
-    psVector *perm = psVectorAlloc(numParams, PS_TYPE_F64);
+    psVector *perm = NULL;
 
     psVector *paramDeltasF64 = psVectorAlloc(numParams, PS_TYPE_F64);
@@ -738 +738 @@
         // XXX: How do we know if these functions were successful?
         //
-        aOut = psMatrixLUD(aOut, perm, A);
+        aOut = psMatrixLUD(aOut, &perm, A);
         paramDeltasF64 = psMatrixLUSolve(paramDeltasF64, aOut, beta, perm);
 
@@ -936 +936 @@
     coeffs = psVectorAlloc(polyOrder, PS_TYPE_F64);
     X = psVectorAlloc(x->n, PS_TYPE_F64);
-    outPerm = psVectorAlloc(polyOrder, PS_TYPE_F64);
     xSums = psVectorAlloc(1 + 2 * polyOrder, PS_TYPE_F64);
 
@@ -943 +942 @@
         B->data.F64[i] = 0.0;
         coeffs->data.F64[i] = 0.0;
-        outPerm->data.F64[i] = 0.0;
         for (j = 0; j < polyOrder; j++) {
             A->data.F64[i][j] = 0.0;
@@ -987 +985 @@
 
     // XXX: How do we know if these routines were successful?
-    ALUD = psMatrixLUD(ALUD, outPerm, A);
+    ALUD = psMatrixLUD(ALUD, &outPerm, A);
     coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
 
@@ -1139 +1137 @@
 Algorithm:
  
-XXX completely ad hoc:  
+XXX completely ad hoc:
 start with the user-supplied starting parameter and
 call that b.  Calculate a/c as a fractional amount smaller/larger than b.
 Repeat this process until a local minimum is found.
  
-XXX:  
-new algorithm:  
+XXX:
+new algorithm:
 start at x=0, expand in one direction until the function
 decreases.  Then you have two points in the bracket.  Keep going until it
@@ -1151 +1149 @@
 direction.
  
-XXX: 
+XXX:
 This is F32 only.
  
-XXX: 
+XXX:
 output bracket vector should be an input as well.
 *****************************************************************************/

trunk/psLib/src/math/psPolynomial.c

@@ -7 +7 @@
  *  polynomials.  It also contains a Gaussian functions.
  *
- *  @version $Revision: 1.81 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-01-14 23:40:24 $
+ *  @version $Revision: 1.82 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-01-17 22:17:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -1906 +1906 @@
         (tmp->spline)[i] = psPolynomial1DAlloc(order+1, PS_POLYNOMIAL_ORD);
     }
+
+    // This should be set by the psVectorFitSpline1D()
+    tmp->p_psDeriv2 = NULL;
 
     tmp->domains = (psF32 *) psAlloc((bounds->n) * sizeof(psF32));

trunk/psLib/src/math/psSpline.c

@@ -7 +7 @@
  *  polynomials.  It also contains a Gaussian functions.
  *
- *  @version $Revision: 1.81 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-01-14 23:40:24 $
+ *  @version $Revision: 1.82 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-01-17 22:17:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -1906 +1906 @@
         (tmp->spline)[i] = psPolynomial1DAlloc(order+1, PS_POLYNOMIAL_ORD);
     }
+
+    // This should be set by the psVectorFitSpline1D()
+    tmp->p_psDeriv2 = NULL;
 
     tmp->domains = (psF32 *) psAlloc((bounds->n) * sizeof(psF32));

trunk/psLib/test/dataIO

@@ -3 +3 @@
 tst_psFits
 multi.fits
+table.fits

@@ -3 +3 @@
 tst_psFits
 multi.fits
+table.fits

trunk/psLib/test/dataIO/.cvsignore

@@ -3 +3 @@
 tst_psFits
 multi.fits
+table.fits

trunk/psLib/test/dataManip/tst_psMatrix03.c

@@ -14 +14 @@
  *  @author  Ross Harman, MHPCC
  *
- *  @version $Revision: 1.12 $  $Name: not supported by cvs2svn $
- *  @date  $Date: 2004-12-10 20:00:06 $
+ *  @version $Revision: 1.13 $  $Name: not supported by cvs2svn $
+ *  @date  $Date: 2005-01-17 22:17:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -108 +108 @@
     inVector->n = 3;
     luImage32 = (psImage*)psImageAlloc(3, 3, PS_TYPE_F32);
-    perm32 = (psVector*)psVectorAlloc(3, PS_TYPE_F32);
     outVector32 = (psVector*)psVectorAlloc(3, PS_TYPE_F32);
     inVector32 = (psVector*)psVectorAlloc(3, PS_TYPE_F32);
@@ -131 +130 @@
     printPositiveTestHeader(stdout, "psMatrix", "Calculate LU matrix");
     tempImage = luImage;
-    luImage = psMatrixLUD(luImage, perm, inImage);
+    luImage = psMatrixLUD(luImage, &perm, inImage);
     CHECK_MATRIX(luImage);
     if(luImage->type.dimen != PS_DIMEN_IMAGE) {
@@ -140 +139 @@
 
     tempImage32 = luImage32;
-    luImage32 = psMatrixLUD(luImage32, perm32, inImage32);
+    luImage32 = psMatrixLUD(luImage32, &perm32, inImage32);
     CHECK_MATRIX(luImage32);
     if(luImage32->type.dimen != PS_DIMEN_IMAGE) {

trunk/psLib/test/fileUtils

@@ -3 +3 @@
 tst_psFits
 multi.fits
+table.fits

@@ -3 +3 @@
 tst_psFits
 multi.fits
+table.fits

trunk/psLib/test/fileUtils/.cvsignore

@@ -3 +3 @@
 tst_psFits
 multi.fits
+table.fits

@@ -3 +3 @@
 tst_psFits
 multi.fits
+table.fits

@@ -3 +3 @@
 tst_psFits
 multi.fits
+table.fits