Index: trunk/psLib/src/dataManip/psMatrixVectorArithmetic.c =================================================================== --- trunk/psLib/src/dataManip/psMatrixVectorArithmetic.c (revision 1065) +++ trunk/psLib/src/dataManip/psMatrixVectorArithmetic.c (revision 1082) @@ -29,6 +29,6 @@ * @author Ross Harman, MHPCC * - * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-06-21 21:04:24 $ + * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-06-24 18:53:42 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -76,7 +76,7 @@ /*****************************************************************************/ -#define MIN(A,B) (((A) < (B)) ? (A) : (B)) - -#define MAX(A,B) (((A) > (B)) ? (A) : (B)) +#define MIN(A,B)(((double)(A)<(double)(B))?(A):(B)) + +#define MAX(A,B)(((double)(A)>(double)(B))?(A):(B)) @@ -85,13 +85,13 @@ { \ int i = 0; \ - int npt = 0; \ - ps##TYPE *o = NULL; \ - ps##TYPE *i1 = NULL; \ - ps##TYPE *i2 = NULL; \ - npt = ((psVector*)IN1)->n; \ + int n1 = 0; \ + ps##TYPE *o = NULL; \ + ps##TYPE *i1 = NULL; \ + ps##TYPE *i2 = NULL; \ + n1 = ((psVector*)IN1)->n; \ o = ((psVector*)OUT)->data.TYPE; \ i1 = ((psVector*)IN1)->data.TYPE; \ i2 = &((psScalar*)IN2)->data.TYPE; \ - for (i=0; i < npt; i++, o++, i1++) { \ + for (i=0; i < n1; i++, o++, i1++) { \ *o = OP; \ } \ @@ -101,17 +101,69 @@ { \ int i = 0; \ - int npt = 0; \ - ps##TYPE *o = NULL; \ - ps##TYPE *i1 = NULL; \ - ps##TYPE *i2 = NULL; \ - npt = ((psVector*)IN1)->n; \ + int n1 = 0; \ + int n2 = 0; \ + ps##TYPE *o = NULL; \ + ps##TYPE *i1 = NULL; \ + ps##TYPE *i2 = NULL; \ + n1 = ((psVector*)IN1)->n; \ + n2 = ((psVector*)IN2)->n; \ + if(n1 != n2) { \ + psError(__func__, ": Inconsistent element count: %d vs %d", n1, n2); \ + return OUT; \ + } \ o = ((psVector*)OUT)->data.TYPE; \ i1 = ((psVector*)IN1)->data.TYPE; \ i2 = ((psVector*)IN2)->data.TYPE; \ - for (i=0; i < npt; i++, o++, i1++, i2++) { \ + for (i=0; i < n1; i++, o++, i1++, i2++) { \ *o = OP; \ } \ } +#define VECTOR_IMAGE(OUT,IN1,OP,IN2,TYPE) \ +{ \ + int i = 0; \ + int j = 0; \ + int n1 = 0; \ + int numRows2 = 0; \ + int numCols2 = 0; \ + psDimen dim1 = 0; \ + ps##TYPE *o = NULL; \ + ps##TYPE *i1 = NULL; \ + ps##TYPE *i2 = NULL; \ + n1 = ((psVector*)IN1)->n; \ + dim1 = ((psVector*)IN1)->type.dimen; \ + numRows2 = ((psImage*)IN2)->numRows; \ + numCols2 = ((psImage*)IN2)->numCols; \ + \ + if(dim1 == PS_DIMEN_VECTOR) { /* Regular vectors */ \ + if(n1!=numRows2) { \ + psError(__func__, ": Inconsistent element count: %d vs %d", n1, numRows2); \ + return OUT; \ + } \ + \ + i1 = ((psVector*)IN1)->data.TYPE; \ + for(j = 0; j < numRows2; j++, i1++) { \ + o = ((psImage*)OUT)->data.TYPE[j]; \ + i2 = ((psImage*)IN2)->data.TYPE[j]; \ + for(i = 0; i < numCols2; i++, o++, i2++) { \ + *o = OP; \ + } \ + } \ + } else { /* Transposed vectors */ \ + if(n1!=numCols2) { \ + psError(__func__, ": Inconsistent element count: %d vs %d", n1, numCols2); \ + return OUT; \ + } \ + \ + 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 @@ -127,9 +179,9 @@ numRows = ((psImage*)IN1)->numRows; \ numCols = ((psImage*)IN1)->numCols; \ - for(j = 0; j < numCols; j++) { \ + for(j = 0; j < numRows; j++) { \ o = ((psImage*)OUT)->data.TYPE[j]; \ i1 = ((psImage*)IN1)->data.TYPE[j]; \ i2 = &((psScalar*)IN2)->data.TYPE; \ - for(i = 0; i < numRows; i++, o++, i1++) { \ + for(i = 0; i < numCols; i++, o++, i1++) { \ *o = OP; \ } \ @@ -137,20 +189,76 @@ } +#define IMAGE_VECTOR(OUT,IN1,OP,IN2,TYPE) \ +{ \ + int i = 0; \ + int j = 0; \ + int n2 = 0; \ + int numRows1 = 0; \ + int numCols1 = 0; \ + psDimen dim2 = 0; \ + ps##TYPE *o = NULL; \ + ps##TYPE *i1 = NULL; \ + ps##TYPE *i2 = NULL; \ + n2 = ((psVector*)IN2)->n; \ + dim2 = ((psVector*)IN2)->type.dimen; \ + numRows1 = ((psImage*)IN1)->numRows; \ + numCols1 = ((psImage*)IN1)->numCols; \ + \ + if(dim2 == PS_DIMEN_VECTOR) { /* Regular vectors */ \ + if(n2!=numRows1) { \ + psError(__func__, ": Inconsistent element count: %d vs %d", n2, numRows1); \ + return OUT; \ + } \ + \ + i2 = ((psVector*)IN2)->data.TYPE; \ + for(j = 0; j < numRows1; j++, i1++) { \ + o = ((psImage*)OUT)->data.TYPE[j]; \ + i1 = ((psImage*)IN1)->data.TYPE[j]; \ + for(i = 0; i < numCols1; i++, o++, i1++) { \ + *o = OP; \ + } \ + } \ + } else { /* Transposed vectors */ \ + if(n2!=numCols1) { \ + psError(__func__, ": Inconsistent element count: %d vs %d", n2, numCols1); \ + return OUT; \ + } \ + \ + 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; \ + } \ + } \ + } \ +} + #define IMAGE_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++) { \ + int numRows1 = 0; \ + int numCols1 = 0; \ + int numRows2 = 0; \ + int numCols2 = 0; \ + ps##TYPE *o = NULL; \ + ps##TYPE *i1 = NULL; \ + ps##TYPE *i2 = NULL; \ + numRows1 = ((psImage*)IN1)->numRows; \ + numCols1 = ((psImage*)IN1)->numCols; \ + numRows2 = ((psImage*)IN2)->numRows; \ + 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; \ + } \ + for(j = 0; j < numRows1; j++) { \ o = ((psImage*)OUT)->data.TYPE[j]; \ i1 = ((psImage*)IN1)->data.TYPE[j]; \ i2 = ((psImage*)IN2)->data.TYPE[j]; \ - for(i = 0; i < numRows; i++, o++, i1++, i2++) { \ + for(i = 0; i < numCols1; i++, o++, i1++, i2++) { \ *o = OP; \ } \ @@ -161,5 +269,4 @@ #define SCALAR_SCALAR(OUT,IN1,OP,IN2,TYPE) \ { \ - int i = 0; \ ps##TYPE *o = NULL; \ ps##TYPE *i1 = NULL; \ @@ -187,5 +294,4 @@ } -// SCALAR_XXXX operations #define SCALAR_IMAGE(OUT,IN1,OP,IN2,TYPE) \ { \ @@ -242,10 +348,8 @@ break; \ case PS_TYPE_C32: \ - psError(__func__, ": Line %d - C32 support not implemented", __LINE__); \ - /*DIM1##_##DIM2(OUT,IN1,OP,IN2,C32); */ \ + DIM1##_##DIM2(OUT,IN1,OP,IN2,C32); \ break; \ case PS_TYPE_C64: \ - psError(__func__, ": Line %d - C64 support not implemented", __LINE__); \ - /*DIM1##_##DIM2(OUT,IN1,OP,IN2,C64); */ \ + DIM1##_##DIM2(OUT,IN1,OP,IN2,C64); \ break; \ default: \ @@ -257,5 +361,5 @@ LEVEL1(DIM1,DIM2,OUT,IN1,*i1 + *i2,IN2); \ } else if(!strncmp(OP, "-", 1)) { \ - LEVEL1(DIM1, DIM2,OUT,IN1,*i1 - *i2,IN2); \ + LEVEL1(DIM1,DIM2,OUT,IN1,*i1 - *i2,IN2); \ } else if(!strncmp(OP, "*", 1)) { \ LEVEL1(DIM1,DIM2,OUT,IN1,*i1 * *i2,IN2); \ @@ -263,9 +367,13 @@ LEVEL1(DIM1,DIM2,OUT,IN1,*i1 / *i2,IN2); \ } else if(!strncmp(OP, "^", 1)) { \ - LEVEL1(DIM1,DIM2,OUT,IN1,pow(*i1,*i2),IN2); \ + if(PS_IS_PSELEMTYPE_COMPLEX(IN1->type)) { \ + LEVEL1(DIM1,DIM2,OUT,IN1,cpow(*i1,*i2),IN2); \ + } else { \ + LEVEL1(DIM1,DIM2,OUT,IN1,pow(*i1,*i2),IN2); \ + } \ } 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); \ @@ -277,4 +385,7 @@ psDimen dim2 = 0; psDimen dimOut = 0; + psElemType elType1 = 0; + psElemType elType2 = 0; + psElemType elTypeOut = 0; psType *psType1 = NULL; psType *psType2 = NULL; @@ -302,35 +413,47 @@ dim2 = psType2->dimen; dimOut = psTypeOut->dimen; - - if(dimOut == PS_DIMEN_OTHER) { - psError(__func__, ": Line %d - Invalid dimensionality for out arg: %d", __LINE__, dimOut); + elType1 = psType1->type; + elType2 = psType2->type; + elTypeOut = psTypeOut->type; + + if(elType1!=elType2 && elType1!=elTypeOut) { + psError(__func__, ": Line %d - Element types for arguments inconsistent: (%d, %d, %d)", __LINE__, + elType1, elType2, elTypeOut); return out; } - if(dim1 == PS_DIMEN_VECTOR) { - if(dim2 == PS_DIMEN_VECTOR) { - LEVEL2(VECTOR,VECTOR,out,psType1,op,psType2); // vector op vector + if(dim1==PS_DIMEN_OTHER || dim2==PS_DIMEN_OTHER || dimOut==PS_DIMEN_OTHER) { + psError(__func__, ": Line %d - PS_DIMEN_OTHER not allowed for arguments: (%d, %d, %d)", __LINE__, + dim1, dim2, dimOut); + return out; + } + + if(dim1 == PS_DIMEN_SCALAR) { + if(dim2 == PS_DIMEN_SCALAR) { + LEVEL2(SCALAR,SCALAR,out,psType1,op,psType2); // scalar op scalar + } else if(dim2==PS_DIMEN_VECTOR || dim2==PS_DIMEN_TRANSV) { + LEVEL2(SCALAR,VECTOR,out,psType1,op,psType2); // scalar op vector } else if(dim2 == PS_DIMEN_IMAGE) { - //LEVEL2(); // vector op image - } else if(dim2 == PS_DIMEN_SCALAR) { - LEVEL2(VECTOR,VECTOR,out,psType1,op,psType2); // vector op scalar + LEVEL2(SCALAR,IMAGE,out,psType1,op,psType2); // scalar op image + } else { + psError(__func__, ": Line %d - Invalid dimensionality for in2 arg: %d", __LINE__, dim2); + } + } else if(dim1==PS_DIMEN_VECTOR || dim1==PS_DIMEN_TRANSV) { + if(dim2 == PS_DIMEN_SCALAR) { + LEVEL2(VECTOR,SCALAR,out,psType1,op,psType2); // vector op scalar + } else if(dim2==PS_DIMEN_VECTOR || dim2==PS_DIMEN_TRANSV) { + LEVEL2(VECTOR,VECTOR,out,psType1,op,psType2); // vector op vector + } else if(dim2 == PS_DIMEN_IMAGE) { + LEVEL2(VECTOR,IMAGE,out,psType1,op,psType2); // vector op image } else { psError(__func__, ": Line %d - Invalid dimensionality for in2 arg: %d", __LINE__, dim2); } } else if(dim1 == PS_DIMEN_IMAGE) { - if(dim2 == PS_DIMEN_VECTOR) { - //LEVEL2(); // image op vector + if(dim2 == PS_DIMEN_SCALAR) { + LEVEL2(IMAGE,SCALAR,out,psType1,op,psType2); // image op scalar + } else if(dim2==PS_DIMEN_VECTOR || dim2==PS_DIMEN_TRANSV) { + LEVEL2(IMAGE,VECTOR,out,psType1,op,psType2); // image op vector } else if(dim2 == PS_DIMEN_IMAGE) { - LEVEL2(IMAGE,IMAGE,out,psType1,op,psType2); // image op image - } else if(dim2 == PS_DIMEN_SCALAR) { - LEVEL2(IMAGE,SCALAR,out,psType1,op,psType2); // image op scalar - } else { - psError(__func__, ": Line %d - Invalid dimensionality for in2 arg: %d", __LINE__, dim2); - } - } else if(dim1 == PS_DIMEN_SCALAR) { - if(dim2 == PS_DIMEN_VECTOR) { - LEVEL2(SCALAR,VECTOR,out,psType1,op,psType2); // scalar op vector - } else if(dim2 == PS_DIMEN_IMAGE) { - LEVEL2(SCALAR,IMAGE,out,psType1,op,psType2); // scalar op image + LEVEL2(IMAGE,IMAGE,out,psType1,op,psType2); // image op image } else { psError(__func__, ": Line %d - Invalid dimensionality for in2 arg: %d", __LINE__, dim2); @@ -342,6 +465,2 @@ return out; } - - - -