Index: /trunk/psLib/src/dataManip/psMatrixVectorArithmetic.c =================================================================== --- /trunk/psLib/src/dataManip/psMatrixVectorArithmetic.c (revision 1084) +++ /trunk/psLib/src/dataManip/psMatrixVectorArithmetic.c (revision 1085) @@ -29,6 +29,6 @@ * @author Ross Harman, MHPCC * - * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-06-24 18:53:42 $ + * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-06-25 00:31:12 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -76,10 +76,66 @@ /*****************************************************************************/ +// Local version of min function (cast to double for complex numbers) #define MIN(A,B)(((double)(A)<(double)(B))?(A):(B)) +// Local version of max function (cast to double for complex numbers) #define MAX(A,B)(((double)(A)>(double)(B))?(A):(B)) - -// VECTOR_XXXX operations +// Conversion for degrees to radians +#define D2R M_PI/360.0 + +// Conversion for radians to degrees +#define R2D 360.0/M_PI + +// Binary SCALAR_XXXX operations +#define SCALAR_SCALAR(OUT,IN1,OP,IN2,TYPE) \ +{ \ + ps##TYPE *o = NULL; \ + ps##TYPE *i1 = NULL; \ + ps##TYPE *i2 = NULL; \ + o = &((psScalar*)OUT)->data.TYPE; \ + i1 = &((psScalar*)IN1)->data.TYPE; \ + i2 = &((psScalar*)IN2)->data.TYPE; \ + *o = OP; \ +} + +#define SCALAR_VECTOR(OUT,IN1,OP,IN2,TYPE) \ +{ \ + int i = 0; \ + int npt = 0; \ + ps##TYPE *o = NULL; \ + ps##TYPE *i1 = NULL; \ + ps##TYPE *i2 = NULL; \ + npt = ((psVector*)IN1)->n; \ + o = ((psVector*)OUT)->data.TYPE; \ + i1 = &((psScalar*)IN1)->data.TYPE; \ + i2 = ((psVector*)IN2)->data.TYPE; \ + for (i=0; i < npt; i++, o++, i2++) { \ + *o = OP; \ + } \ +} + +#define SCALAR_IMAGE(OUT,IN1,OP,IN2,TYPE) \ +{ \ + int i = 0; \ + int j = 0; \ + int numRows = 0; \ + int numCols = 0; \ + ps##TYPE *o = NULL; \ + ps##TYPE *i1 = NULL; \ + ps##TYPE *i2 = NULL; \ + numRows = ((psImage*)IN1)->numRows; \ + numCols = ((psImage*)IN1)->numCols; \ + for(j = 0; j < numCols; j++) { \ + o = ((psImage*)OUT)->data.TYPE[j]; \ + i1 = &((psScalar*)IN1)->data.TYPE; \ + i2 = ((psImage*)IN2)->data.TYPE[j]; \ + for(i = 0; i < numRows; i++, o++, i2++) { \ + *o = OP; \ + } \ + } \ +} + +// Binary VECTOR_XXXX operations #define VECTOR_SCALAR(OUT,IN1,OP,IN2,TYPE) \ { \ @@ -109,6 +165,6 @@ n2 = ((psVector*)IN2)->n; \ if(n1 != n2) { \ - psError(__func__, ": Inconsistent element count: %d vs %d", n1, n2); \ - return OUT; \ + psError(__func__, ": Inconsistent element count: %d vs %d", n1, n2); \ + return OUT; \ } \ o = ((psVector*)OUT)->data.TYPE; \ @@ -138,6 +194,6 @@ if(dim1 == PS_DIMEN_VECTOR) { /* Regular vectors */ \ if(n1!=numRows2) { \ - psError(__func__, ": Inconsistent element count: %d vs %d", n1, numRows2); \ - return OUT; \ + psError(__func__, ": Inconsistent element count: %d vs %d", n1, numRows2); \ + return OUT; \ } \ \ @@ -156,16 +212,16 @@ } \ \ - for(j = 0; j < numRows2; j++) { \ - o = ((psImage*)OUT)->data.TYPE[j]; \ - i1 = ((psVector*)IN1)->data.TYPE; \ - i2 = ((psImage*)IN2)->data.TYPE[j]; \ - for(i = 0; i < numCols2; i++, o++, i1++, i2++) { \ - *o = OP; \ - } \ - } \ - } \ -} - -// IMAGE_XXXX operations + for(j = 0; j < numRows2; j++) { \ + o = ((psImage*)OUT)->data.TYPE[j]; \ + i1 = ((psVector*)IN1)->data.TYPE; \ + i2 = ((psImage*)IN2)->data.TYPE[j]; \ + for(i = 0; i < numCols2; i++, o++, i1++, i2++) { \ + *o = OP; \ + } \ + } \ + } \ +} + +// Binary IMAGE_XXXX operations #define IMAGE_SCALAR(OUT,IN1,OP,IN2,TYPE) \ { \ @@ -207,6 +263,6 @@ if(dim2 == PS_DIMEN_VECTOR) { /* Regular vectors */ \ if(n2!=numRows1) { \ - psError(__func__, ": Inconsistent element count: %d vs %d", n2, numRows1); \ - return OUT; \ + psError(__func__, ": Inconsistent element count: %d vs %d", n2, numRows1); \ + return OUT; \ } \ \ @@ -225,12 +281,12 @@ } \ \ - for(j = 0; j < numRows1; j++) { \ - o = ((psImage*)OUT)->data.TYPE[j]; \ - i1 = ((psVector*)IN2)->data.TYPE; \ - i2 = ((psImage*)IN1)->data.TYPE[j]; \ - for(i = 0; i < numCols1; i++, o++, i2++, i1++) { \ - *o = OP; \ - } \ - } \ + for(j = 0; j < numRows1; j++) { \ + o = ((psImage*)OUT)->data.TYPE[j]; \ + i1 = ((psVector*)IN2)->data.TYPE; \ + i2 = ((psImage*)IN1)->data.TYPE[j]; \ + for(i = 0; i < numCols1; i++, o++, i2++, i1++) { \ + *o = OP; \ + } \ + } \ } \ } @@ -252,7 +308,7 @@ numCols2 = ((psImage*)IN2)->numCols; \ if(numRows1!=numRows2 || numCols1!=numCols2) { \ - psError(__func__, ": Inconsistent element count: numRows: %d vs %d numCols: %d vs %d", numRows1, \ - numRows2, numCols1, numCols2); \ - return OUT; \ + psError(__func__, ": Inconsistent element count: numRows: %d vs %d numCols: %d vs %d", numRows1, \ + numRows2, numCols1, numCols2); \ + return OUT; \ } \ for(j = 0; j < numRows1; j++) { \ @@ -266,53 +322,5 @@ } -// SCALAR_XXXX operations -#define SCALAR_SCALAR(OUT,IN1,OP,IN2,TYPE) \ -{ \ - ps##TYPE *o = NULL; \ - ps##TYPE *i1 = NULL; \ - ps##TYPE *i2 = NULL; \ - o = &((psScalar*)OUT)->data.TYPE; \ - i1 = &((psScalar*)IN1)->data.TYPE; \ - i2 = &((psScalar*)IN2)->data.TYPE; \ - *o = OP; \ -} - -#define SCALAR_VECTOR(OUT,IN1,OP,IN2,TYPE) \ -{ \ - int i = 0; \ - int npt = 0; \ - ps##TYPE *o = NULL; \ - ps##TYPE *i1 = NULL; \ - ps##TYPE *i2 = NULL; \ - npt = ((psVector*)IN1)->n; \ - o = ((psVector*)OUT)->data.TYPE; \ - i1 = &((psScalar*)IN1)->data.TYPE; \ - i2 = ((psVector*)IN2)->data.TYPE; \ - for (i=0; i < npt; i++, o++, i2++) { \ - *o = OP; \ - } \ -} - -#define SCALAR_IMAGE(OUT,IN1,OP,IN2,TYPE) \ -{ \ - int i = 0; \ - int j = 0; \ - int numRows = 0; \ - int numCols = 0; \ - ps##TYPE *o = NULL; \ - ps##TYPE *i1 = NULL; \ - ps##TYPE *i2 = NULL; \ - numRows = ((psImage*)IN1)->numRows; \ - numCols = ((psImage*)IN1)->numCols; \ - for(j = 0; j < numCols; j++) { \ - o = ((psImage*)OUT)->data.TYPE[j]; \ - i1 = &((psScalar*)IN1)->data.TYPE; \ - i2 = ((psImage*)IN2)->data.TYPE[j]; \ - for(i = 0; i < numRows; i++, o++, i2++) { \ - *o = OP; \ - } \ - } \ -} - +// Preprocessor macro function to create arithmetic function based on input type #define LEVEL1(DIM1,DIM2,OUT,IN1,OP,IN2) \ switch (IN1->type) { \ @@ -357,4 +365,5 @@ } +// Preprocessor macro function to create arithmetic function operation name #define LEVEL2(DIM1,DIM2,OUT,IN1,OP,IN2) \ if(!strncmp(OP, "+", 1)) { \ @@ -373,7 +382,7 @@ } \ } else if(!strncmp(OP, "min", 3)) { \ - LEVEL1(DIM1,DIM2,OUT,IN1,MIN(*i1,*i2),IN2); \ + LEVEL1(DIM1,DIM2,OUT,IN1,MIN(*i1,*i2),IN2); \ } else if(!strncmp(OP, "max", 3)) { \ - LEVEL1(DIM1,DIM2,OUT,IN1,MAX(*i1,*i2),IN2); \ + LEVEL1(DIM1,DIM2,OUT,IN1,MAX(*i1,*i2),IN2); \ } else { \ psError(__func__, ": Line %d - Invalid operator: %s", __LINE__, #OP); \ @@ -407,4 +416,9 @@ if(psType2 == NULL) { psError(__func__, ": Line %d - Null in2 argument", __LINE__); + return out; + } + + if(op == NULL) { + psError(__func__, ": Line %d - Null op argument", __LINE__); return out; } @@ -465,2 +479,271 @@ return out; } + +// Unary SCALAR operations +#define SCALAR(OUT,IN,OP,TYPE) \ +{ \ + ps##TYPE *o = NULL; \ + ps##TYPE *i1 = NULL; \ + o = &((psScalar*)OUT)->data.TYPE; \ + i1 = &((psScalar*)IN)->data.TYPE; \ + *o = OP; \ +} + +// Unary IMAGE operations +#define VECTOR(OUT,IN,OP,TYPE) \ +{ \ + int i = 0; \ + int nIn = 0; \ + int nOut = 0; \ + ps##TYPE *o = NULL; \ + ps##TYPE *i1 = NULL; \ + nIn = ((psVector*)IN)->n; \ + nOut = ((psVector*)OUT)->n; \ + if(nIn != nOut) { \ + psError(__func__, ": Inconsistent element count: %d vs %d", nIn, nOut); \ + return OUT; \ + } \ + o = ((psVector*)OUT)->data.TYPE; \ + i1 = ((psVector*)IN)->data.TYPE; \ + for(i = 0; i < nIn; i++, o++, i1++) { \ + *o = OP; \ + } \ +} + +// Unary IMAGE operations +#define IMAGE(OUT,IN,OP,TYPE) \ +{ \ + int i = 0; \ + int j = 0; \ + int numRowsIn = 0; \ + int numColsIn = 0; \ + int numRowsOut = 0; \ + int numColsOut = 0; \ + ps##TYPE *o = NULL; \ + ps##TYPE *i1 = NULL; \ + numRowsIn = ((psImage*)IN)->numRows; \ + numColsIn = ((psImage*)IN)->numCols; \ + numRowsOut = ((psImage*)OUT)->numRows; \ + numColsOut = ((psImage*)OUT)->numCols; \ + if(numRowsIn!=numRowsOut || numColsIn!=numColsOut) { \ + psError(__func__, ": Inconsistent element count: numRows: %d vs %d numCols: %d vs %d", numRowsIn, \ + numRowsOut, numColsIn, numColsOut); \ + return OUT; \ + } \ + for(j = 0; j < numRowsIn; j++) { \ + o = ((psImage*)OUT)->data.TYPE[j]; \ + i1 = ((psImage*)IN)->data.TYPE[j]; \ + for(i = 0; i < numColsIn; i++, o++, i1++) { \ + *o = OP; \ + } \ + } \ +} + + +// Preprocessor macro function to create arithmetic function based on input type +#define UNARY_TYPE(DIM,OUT,IN,OP) \ +switch (IN->type) { \ +case PS_TYPE_S8: \ + DIM(OUT,IN,OP,S8); \ + break; \ +case PS_TYPE_U8: \ + DIM(OUT,IN,OP,U8); \ + break; \ +case PS_TYPE_S16: \ + DIM(OUT,IN,OP,S16); \ + break; \ +case PS_TYPE_U16: \ + DIM(OUT,IN,OP,U16); \ + break; \ +case PS_TYPE_S32: \ + DIM(OUT,IN,OP,S32); \ + break; \ +case PS_TYPE_U32: \ + DIM(OUT,IN,OP,U32); \ + break; \ +case PS_TYPE_S64: \ + DIM(OUT,IN,OP,S64); \ + break; \ +case PS_TYPE_U64: \ + DIM(OUT,IN,OP,U64); \ + break; \ +case PS_TYPE_F32: \ + DIM(OUT,IN,OP,F32); \ + break; \ +case PS_TYPE_F64: \ + DIM(OUT,IN,OP,F64); \ + break; \ +case PS_TYPE_C32: \ + DIM(OUT,IN,OP,C32); \ + break; \ +case PS_TYPE_C64: \ + DIM(OUT,IN,OP,C64); \ + break; \ +default: \ + psError(__func__, ": Line %d - Invalid PS_TYPE: %d", __LINE__, IN->type); \ +} + + +// Preprocessor macro function to create arithmetic function operation name +#define UNARY_OP(DIM,OUT,IN,OP) \ +if(!strncmp(OP, "abs", 3)) { \ + if(PS_IS_PSELEMTYPE_COMPLEX(IN->type)) { \ + UNARY_TYPE(DIM,OUT,IN,cabs(*i1)); \ + } else { \ + UNARY_TYPE(DIM,OUT,IN,fabs(*i1)); \ + } \ +} else if(!strncmp(OP, "exp", 3)) { \ + if(PS_IS_PSELEMTYPE_COMPLEX(IN->type)) { \ + UNARY_TYPE(DIM,OUT,IN,exp(*i1)); \ + } else { \ + UNARY_TYPE(DIM,OUT,IN,cexp(*i1)); \ + } \ +} else if(!strncmp(OP, "ln", 2)) { \ + if(PS_IS_PSELEMTYPE_COMPLEX(IN->type)) { \ + UNARY_TYPE(DIM,OUT,IN,log(*i1)); \ + } else { \ + UNARY_TYPE(DIM,OUT,IN,clog(*i1)); \ + } \ +} else if(!strncmp(OP, "ten", 3)) { \ + if(PS_IS_PSELEMTYPE_COMPLEX(IN->type)) { \ + UNARY_TYPE(DIM,OUT,IN,exp10(*i1)); \ + } else { \ + UNARY_TYPE(DIM,OUT,IN,cpow(*i1,10.0)); \ + } \ +} else if(!strncmp(OP, "log", 3)) { \ + if(PS_IS_PSELEMTYPE_COMPLEX(IN->type)) { \ + UNARY_TYPE(DIM,OUT,IN,log10(*i1)); \ + } else { \ + UNARY_TYPE(DIM,OUT,IN,clog10(*i1)); \ + } \ +} else if(!strncmp(OP, "sin", 3)) { \ + if(PS_IS_PSELEMTYPE_COMPLEX(IN->type)) { \ + UNARY_TYPE(DIM,OUT,IN,sin(*i1)); \ + } else { \ + UNARY_TYPE(DIM,OUT,IN,csin(*i1)); \ + } \ +} else if(!strncmp(OP, "dsin", 4)) { \ + if(PS_IS_PSELEMTYPE_COMPLEX(IN->type)) { \ + UNARY_TYPE(DIM,OUT,IN,sin(*i1*D2R)); \ + } else { \ + UNARY_TYPE(DIM,OUT,IN,csin(*i1*D2R)); \ + } \ +} else if(!strncmp(OP, "cos", 3)) { \ + if(PS_IS_PSELEMTYPE_COMPLEX(IN->type)) { \ + UNARY_TYPE(DIM,OUT,IN,cos(*i1)); \ + } else { \ + UNARY_TYPE(DIM,OUT,IN,ccos(*i1)); \ + } \ +} else if(!strncmp(OP, "dcos", 4)) { \ + if(PS_IS_PSELEMTYPE_COMPLEX(IN->type)) { \ + UNARY_TYPE(DIM,OUT,IN,cos(*i1*D2R)); \ + } else { \ + UNARY_TYPE(DIM,OUT,IN,ccos(*i1*D2R)); \ + } \ +} else if(!strncmp(OP, "tan", 3)) { \ + if(PS_IS_PSELEMTYPE_COMPLEX(IN->type)) { \ + UNARY_TYPE(DIM,OUT,IN,tan(*i1)); \ + } else { \ + UNARY_TYPE(DIM,OUT,IN,ctan(*i1)); \ + } \ +} else if(!strncmp(OP, "dtan", 4)) { \ + if(PS_IS_PSELEMTYPE_COMPLEX(IN->type)) { \ + UNARY_TYPE(DIM,OUT,IN,tan(*i1*D2R)); \ + } else { \ + UNARY_TYPE(DIM,OUT,IN,ctan(*i1*D2R)); \ + } \ +} else if(!strncmp(OP, "asin", 4)) { \ + if(PS_IS_PSELEMTYPE_COMPLEX(IN->type)) { \ + UNARY_TYPE(DIM,OUT,IN,asin(*i1)); \ + } else { \ + UNARY_TYPE(DIM,OUT,IN,casin(*i1)); \ + } \ +} else if(!strncmp(OP, "dasin", 5)) { \ + if(PS_IS_PSELEMTYPE_COMPLEX(IN->type)) { \ + UNARY_TYPE(DIM,OUT,IN,R2D*asin(*i1)); \ + } else { \ + UNARY_TYPE(DIM,OUT,IN,R2D*casin(*i1)); \ + } \ +} else if(!strncmp(OP, "acos", 4)) { \ + if(PS_IS_PSELEMTYPE_COMPLEX(IN->type)) { \ + UNARY_TYPE(DIM,OUT,IN,acos(*i1)); \ + } else { \ + UNARY_TYPE(DIM,OUT,IN,cacos(*i1)); \ + } \ +} else if(!strncmp(OP, "dacos", 5)) { \ + if(PS_IS_PSELEMTYPE_COMPLEX(IN->type)) { \ + UNARY_TYPE(DIM,OUT,IN,R2D*acos(*i1)); \ + } else { \ + UNARY_TYPE(DIM,OUT,IN,R2D*cacos(*i1)); \ + } \ +} else if(!strncmp(OP, "atan", 4)) { \ + if(PS_IS_PSELEMTYPE_COMPLEX(IN->type)) { \ + UNARY_TYPE(DIM,OUT,IN,atan(*i1)); \ + } else { \ + UNARY_TYPE(DIM,OUT,IN,catan(*i1)); \ + } \ +} else if(!strncmp(OP, "datan", 5)) { \ + if(PS_IS_PSELEMTYPE_COMPLEX(IN->type)) { \ + UNARY_TYPE(DIM,OUT,IN,R2D*atan(*i1)); \ + } else { \ + UNARY_TYPE(DIM,OUT,IN,R2D*catan(*i1)); \ + } \ +} else { \ + psError(__func__, ": Line %d - Invalid operator: %s", __LINE__, #OP); \ +} + +void *psUnaryOp(void *out, void *in, char *op) +{ + psDimen dimIn = 0; + psDimen dimOut = 0; + psElemType elTypeIn = 0; + psElemType elTypeOut = 0; + psType *psTypeIn = NULL; + psType *psTypeOut = NULL; + + psTypeOut = (psType*)out; + if(psTypeOut == NULL) { + psError(__func__, ": Line %d - Null out argument", __LINE__); + return out; + } + + psTypeIn = (psType*)in; + if(psTypeIn == NULL) { + psError(__func__, ": Line %d - Null in argument", __LINE__); + return out; + } + + if(op == NULL) { + psError(__func__, ": Line %d - Null op argument", __LINE__); + return out; + } + + dimIn = psTypeIn->dimen; + dimOut = psTypeOut->dimen; + elTypeIn = psTypeIn->type; + elTypeOut = psTypeOut->type; + + if(elTypeIn!=elTypeOut) { + psError(__func__, ": Line %d - Element types for arguments inconsistent: (%d, %d)", __LINE__, + elTypeIn, elTypeOut); + return out; + } + + if(dimIn==PS_DIMEN_OTHER || dimOut==PS_DIMEN_OTHER) { + psError(__func__, ": Line %d - PS_DIMEN_OTHER not allowed for arguments: (%d, %d)", __LINE__, + dimIn, dimOut); + return out; + } + + if(dimIn == PS_DIMEN_SCALAR) { + UNARY_OP(SCALAR,out,psTypeIn,op); // scalar + } else if(dimIn==PS_DIMEN_VECTOR || dimIn==PS_DIMEN_TRANSV) { + UNARY_OP(VECTOR,out,psTypeIn,op); // vector + } else if(dimIn == PS_DIMEN_IMAGE) { + UNARY_OP(IMAGE,out,psTypeIn,op); // image + } else { + psError(__func__, ": Line %d - Invalid dimensionality for in arg: %d", __LINE__, dimIn); + } + + return out; +}