diff --git a/lasp/c/CMakeLists.txt b/lasp/c/CMakeLists.txt index 58cf81d..84910ff 100644 --- a/lasp/c/CMakeLists.txt +++ b/lasp/c/CMakeLists.txt @@ -5,7 +5,7 @@ endif(!LASP_DEBUG) add_library(lasp_lib lasp_fft.c - lasp_math.c + lasp_mat.c lasp_math_raw.c lasp_alg.c lasp_assert.c diff --git a/lasp/c/lasp_alg.c b/lasp/c/lasp_alg.c index a51e300..6947b3e 100644 --- a/lasp/c/lasp_alg.c +++ b/lasp/c/lasp_alg.c @@ -8,10 +8,12 @@ #include "lasp_alg.h" void cmv_dot(const cmat* A,const vc* restrict x,vc* restrict b){ + dbgassert(A && x && b,NULLPTRDEREF); + assert_vx(b); + assert_vx(x); + dbgassert(A->n_cols == x->n_rows,SIZEINEQUAL); + dbgassert(A->n_rows == b->n_rows,SIZEINEQUAL); - dbgassert(A->n_rows == b->size,SIZEINEQUAL); - dbgassert(A->n_cols == x->size,SIZEINEQUAL); - #if LASP_USE_BLAS == 1 dbgassert(false,"Untested function. Is not functional for strides"); /* typedef enum CBLAS_ORDER {CblasRowMajor=101, CblasColMajor=102} CBLAS_ORDER; */ diff --git a/lasp/c/lasp_alg.h b/lasp/c/lasp_alg.h index 6509f6e..1443576 100644 --- a/lasp/c/lasp_alg.h +++ b/lasp/c/lasp_alg.h @@ -8,7 +8,7 @@ #pragma once #ifndef LASP_ALG_H #define LASP_ALG_H -#include "lasp_math.h" +#include "lasp_mat.h" /** * Compute the dot product of two vectors of floats @@ -19,8 +19,10 @@ */ static inline d vd_dot(const vd * a,const vd* b) { dbgassert(a && b,NULLPTRDEREF); - dbgassert(a->size == b->size,SIZEINEQUAL); - return d_dot(getvdval(a,0),getvdval(b,0),a->size); + assert_vx(a); + assert_vx(b); + dbgassert(a->n_rows == b->n_rows,SIZEINEQUAL); + return d_dot(getvdval(a,0),getvdval(b,0),a->n_rows); } /** @@ -93,11 +95,10 @@ static inline void cmat_add_cmat(cmat* x,cmat* y,c fac) { */ static inline void vd_elem_prod(vd* result,const vd* a,const vd* b) { dbgassert(result && a && b,NULLPTRDEREF); - dbgassert(result->size==a->size,SIZEINEQUAL); - dbgassert(b->size==a->size,SIZEINEQUAL); + assert_equalsize(a,b); d_elem_prod_d(getvdval(result,0), getvdval(a,0), - getvdval(b,0),a->size); + getvdval(b,0),a->n_rows); } /** * Compute the element-wise (Hadamard) product of a and b, and store @@ -110,12 +111,13 @@ static inline void vd_elem_prod(vd* result,const vd* a,const vd* b) { static inline void vc_hadamard(vc* result,const vc* a,const vc* b) { fsTRACE(15); dbgassert(result && a && b,NULLPTRDEREF); - dbgassert(result->size==a->size,SIZEINEQUAL); - dbgassert(b->size==a->size,SIZEINEQUAL); + assert_equalsize(a,b); + assert_equalsize(a,result); + c_hadamard(getvcval(result,0), getvcval(a,0), getvcval(b,0), - a->size); + a->n_rows); check_overflow_vx(*result); check_overflow_vx(*a); diff --git a/lasp/c/lasp_aps.c b/lasp/c/lasp_aps.c index 537710c..75b83a5 100644 --- a/lasp/c/lasp_aps.c +++ b/lasp/c/lasp_aps.c @@ -64,7 +64,7 @@ void AvPowerSpectra_free(AvPowerSpectra* aps) { dFifo_free(aps->fifo); dmat_free(&aps->buffer); - us nweight = aps->weighting.size; + us nweight = aps->weighting.n_rows; if(nweight > 0) { for(us blockno = 0; blockno < nweight; blockno++) { cmat_free(&aps->ps_storage[blockno]); @@ -127,7 +127,7 @@ AvPowerSpectra* AvPowerSpectra_alloc(const us nfft, aps->oldest_block = 0; if(weighting) { - us nweight = weighting->size; + us nweight = weighting->n_rows; iVARTRACE(15,nweight); /* Allocate vectors and matrices */ aps->ps_storage = a_malloc(nweight*sizeof(cmat)); @@ -137,12 +137,12 @@ AvPowerSpectra* AvPowerSpectra_alloc(const us nfft, } /* Allocate space and copy weighting coefficients */ - aps->weighting = vd_alloc(weighting->size); + aps->weighting = vd_alloc(weighting->n_rows); vd_copy(&aps->weighting,weighting); } else { TRACE(15,"no weighting"); - aps->weighting.size = 0; + aps->weighting.n_rows = 0; } aps->ps_result = cmat_alloc(nfft/2+1,nchannels*nchannels); @@ -186,7 +186,7 @@ static void AvPowerSpectra_addBlock(AvPowerSpectra* aps, ps_single); vd weighting = aps->weighting; - us nweight = weighting.size; + us nweight = weighting.n_rows; if(nweight == 0) { diff --git a/lasp/c/lasp_aps.h b/lasp/c/lasp_aps.h index e7215ba..973037a 100644 --- a/lasp/c/lasp_aps.h +++ b/lasp/c/lasp_aps.h @@ -9,7 +9,7 @@ #ifndef LASP_APS_H #define LASP_APS_H #include "lasp_types.h" -#include "lasp_math.h" +#include "lasp_mat.h" #include "lasp_window.h" typedef enum { diff --git a/lasp/c/lasp_decimation.c b/lasp/c/lasp_decimation.c index af78ad6..c444b3b 100644 --- a/lasp/c/lasp_decimation.c +++ b/lasp/c/lasp_decimation.c @@ -62,7 +62,7 @@ Decimator* Decimator_create(us nchannels,DEC_FAC df) { dec->nchannels = nchannels; dec->dec_fac = filter->dec_fac; - dmat h = dmat_foreign(filter->ntaps,1,filter->h); + dmat h = dmat_foreign_data(filter->ntaps,1,filter->h,false); for(us channelno=0;channelnofbs[channelno] = FilterBank_create(&h,DEC_FFT_LEN); diff --git a/lasp/c/lasp_decimation.h b/lasp/c/lasp_decimation.h index 207a4dc..4dcfc01 100644 --- a/lasp/c/lasp_decimation.h +++ b/lasp/c/lasp_decimation.h @@ -10,7 +10,7 @@ #ifndef LASP_DECIMATION_H #define LASP_DECIMATION_H #include "lasp_types.h" -#include "lasp_math.h" +#include "lasp_mat.h" typedef struct Decimator_s Decimator; diff --git a/lasp/c/lasp_dfifo.h b/lasp/c/lasp_dfifo.h index 587d2aa..98f198f 100644 --- a/lasp/c/lasp_dfifo.h +++ b/lasp/c/lasp_dfifo.h @@ -9,7 +9,7 @@ #ifndef LASP_DFIFO_H #define LASPDFIFO_H #include "lasp_types.h" -#include "lasp_math.h" +#include "lasp_mat.h" typedef struct dFifo_s dFifo; diff --git a/lasp/c/lasp_fft.c b/lasp/c/lasp_fft.c index 162c0c9..f92374a 100644 --- a/lasp/c/lasp_fft.c +++ b/lasp/c/lasp_fft.c @@ -54,11 +54,11 @@ void Fft_ifft_single(const Fft* fft,const vc* freqdata,vd* result) { fsTRACE(15); dbgassert(fft && freqdata && result,NULLPTRDEREF); const us nfft = fft->nfft; - dbgassert(result->size == nfft, + dbgassert(result->n_rows == nfft, "Invalid size for time data rows." " Should be equal to nfft"); - dbgassert(freqdata->size == (nfft/2+1),"Invalid number of rows in" + dbgassert(freqdata->n_rows == (nfft/2+1),"Invalid number of rows in" " result array"); @@ -113,11 +113,13 @@ void Fft_fft_single(const Fft* fft,const vd* timedata,vc* result) { dbgassert(fft && timedata && result,NULLPTRDEREF); const us nfft = fft->nfft; - dbgassert(timedata->size == nfft, + assert_vx(timedata); + assert_vx(result); + dbgassert(timedata->n_rows == nfft, "Invalid size for time data rows." " Should be equal to nfft"); - dbgassert(result->size == (nfft/2+1),"Invalid number of rows in" + dbgassert(result->n_rows == (nfft/2+1),"Invalid number of rows in" " result array"); d* result_ptr = (d*) getvcval(result,0); diff --git a/lasp/c/lasp_fft.h b/lasp/c/lasp_fft.h index cf4506d..5da9e1d 100644 --- a/lasp/c/lasp_fft.h +++ b/lasp/c/lasp_fft.h @@ -9,7 +9,7 @@ #ifndef LASP_FFT_H #define LASP_FFT_H #include "lasp_types.h" -#include "lasp_math.h" +#include "lasp_mat.h" /** * Perform forward FFT's on real time data. diff --git a/lasp/c/lasp_filterbank.c b/lasp/c/lasp_filterbank.c index 10a1a89..dcfff9f 100644 --- a/lasp/c/lasp_filterbank.c +++ b/lasp/c/lasp_filterbank.c @@ -100,9 +100,7 @@ dmat FilterBank_filter(FilterBank* fb, const us nfilters = fb->filters.n_cols; /* Push samples to the input fifo */ - dmat samples = dmat_foreign_vd((vd*) x); - dFifo_push(fb->input_fifo,&samples); - dmat_free(&samples); + dFifo_push(fb->input_fifo,x); dmat input_block = dmat_alloc(nfft,1); diff --git a/lasp/c/lasp_filterbank.h b/lasp/c/lasp_filterbank.h index 17f40bc..be53344 100644 --- a/lasp/c/lasp_filterbank.h +++ b/lasp/c/lasp_filterbank.h @@ -13,7 +13,7 @@ #ifndef LASP_FILTERBANK_H #define LASP_FILTERBANK_H #include "lasp_types.h" -#include "lasp_math.h" +#include "lasp_mat.h" typedef struct FilterBank_s FilterBank; /** diff --git a/lasp/c/lasp_math.c b/lasp/c/lasp_mat.c similarity index 61% rename from lasp/c/lasp_math.c rename to lasp/c/lasp_mat.c index ff6e415..6c39a15 100644 --- a/lasp/c/lasp_math.c +++ b/lasp/c/lasp_mat.c @@ -1,4 +1,4 @@ -// lasp_math.c +// lasp_mat.c // // Author: J.A. de Jong -ASCEE // @@ -6,10 +6,8 @@ // ////////////////////////////////////////////////////////////////////// #define TRACERPLUS (-10) -#include "lasp_math.h" - +#include "lasp_mat.h" #include "lasp_assert.h" -#include "lasp_math.h" #include "lasp_tracer.h" #include @@ -50,35 +48,6 @@ void print_cmat(const cmat* m) { } feTRACE(50); } -void print_vc(const vc* m) { - fsTRACE(50); - size_t row; - - for(row=0;rowsize;row++){ - c val = *getvcval(m,row); - - d rval = creal(val); - d ival = cimag(val); - indent_trace(); - printf("%c%2.2e%c%2.2ei ",rval< 0 ?'-': ' ', d_abs(rval),ival<0 ? '-' : '+',d_abs(ival) ) ; - printf("\n"); - - } - feTRACE(50); -} -void print_vd(const vd* m) { - fsTRACE(50); - size_t row; - iVARTRACE(20,m->size); - for(row=0;rowsize;row++){ - - d rval = *getvdval(m,row); - indent_trace(); - printf("%c%2.2e ",rval< 0 ? '\r': ' ',rval); - printf("\n"); - } - feTRACE(50); -} #endif diff --git a/lasp/c/lasp_math.h b/lasp/c/lasp_mat.h similarity index 70% rename from lasp/c/lasp_math.h rename to lasp/c/lasp_mat.h index bf83baf..f4e31c9 100644 --- a/lasp/c/lasp_math.h +++ b/lasp/c/lasp_mat.h @@ -1,4 +1,4 @@ -// lasp_math.h +// lasp_mat.h // // Author: J.A. de Jong - ASCEE // @@ -6,29 +6,13 @@ // copying of matrices and vectors. ////////////////////////////////////////////////////////////////////// #pragma once -#ifndef LASP_MATH_H -#define LASP_MATH_H +#ifndef LASP_MAT_H +#define LASP_MAT_H #include "lasp_math_raw.h" #include "lasp_alloc.h" #include "lasp_tracer.h" #include "lasp_assert.h" -/// Vector of floating point numbers -typedef struct { - us size; - bool _foreign_data; - d* _data; /**< Pointer set if data storage is - intern. If this is set to zero, the - vector is a sub-vector. */ -} vd; -/// Vector of complex floating point numbers -typedef struct { - us size; - bool _foreign_data; - c* _data; /**< Pointer set if data storage is - intern. If this is set to zero, the - vector is a sub-vector. */ -} vc; /// Dense matrix of floating point values typedef struct { us n_rows; @@ -46,10 +30,20 @@ typedef struct { c* _data; } cmat; -#define setvecval(vec,index,val) \ - dbgassert((((us) index) <= (vec)->size),OUTOFBOUNDSVEC); \ - (vec)->_data[index] = val; +typedef dmat vd; +typedef cmat vc; +#define assert_equalsize(a,b) \ + dbgassert((a)->n_rows == (b)->n_rows,SIZEINEQUAL); \ + dbgassert((a)->n_cols == (b)->n_cols,SIZEINEQUAL); + +#define is_vx(vx) ((vx)->n_cols == 1) +#define assert_vx(vx) dbgassert(is_vx(vx),"Not a vector!") + +#define setvecval(vec,index,val) \ + assert_vx(vec); \ + dbgassert((((us) index) < (vec)->n_rows),OUTOFBOUNDSVEC); \ + (vec)->_data[index] = val; #define setmatval(mat,row,col,val) \ dbgassert((((us) row) <= mat->n_rows),OUTOFBOUNDSMATR); \ @@ -62,22 +56,6 @@ typedef struct { * @param mat The vector * @param row The row */ -static inline d* getvdval(const vd* vec,us row){ - dbgassert(row < vec->size,OUTOFBOUNDSVEC); - return &vec->_data[row]; -} - -/** - * Return pointer to a value from a complex vector - * - * @param mat The vector - * @param row The row - */ -static inline c* getvcval(const vc* vec,us row){ - dbgassert(row < vec->size,OUTOFBOUNDSVEC); - return &vec->_data[row]; -} - /** * Return a value from a matrix of floating points * @@ -91,7 +69,6 @@ static inline d* getdmatval(const dmat* mat,us row,us col){ dbgassert(col < mat->n_cols,OUTOFBOUNDSMATC); return &mat->_data[col*mat->stride+row]; } - /** * Return a value from a matrix of complex floating points * @@ -106,17 +83,29 @@ static inline c* getcmatval(const cmat* mat,const us row,const us col){ return &mat->_data[col*mat->stride+row]; } +static inline d* getvdval(const vd* vec,us row){ + dbgassert(vec,NULLPTRDEREF); + assert_vx(vec); + return getdmatval(vec,row,0); +} + +/** + * Return pointer to a value from a complex vector + * + * @param mat The vector + * @param row The row + */ +static inline c* getvcval(const vc* vec,us row){ + dbgassert(vec,NULLPTRDEREF); + assert_vx(vec); + return getcmatval(vec,row,0); +} + + + #ifdef LASP_DEBUG #define OVERFLOW_MAGIC_NUMBER (-10e-45) -#define check_overflow_vx(vx) \ - TRACE(15,"Checking overflow " #vx); \ - if(!(vx)._foreign_data) { \ - dbgassert((vx)._data[(vx).size] == OVERFLOW_MAGIC_NUMBER, \ - "Buffer overflow detected on" #vx ); \ - } \ - - #define check_overflow_xmat(xmat) \ TRACE(15,"Checking overflow " #xmat); \ if(!(xmat)._foreign_data) { \ @@ -125,31 +114,13 @@ static inline c* getcmatval(const cmat* mat,const us row,const us col){ "Buffer overflow detected on" #xmat ); \ } \ +#define check_overflow_vx check_overflow_xmat + #else #define check_overflow_vx(vx) #define check_overflow_xmat(xmat) #endif -/** - * Sets all values in a vector to the value - * - * @param vec the vector to set - * @param value - */ -static inline void vd_set(vd* vec, d value){ - d_set(vec->_data,value,vec->size); -} - -/** - * Sets all values in a vector to the value - * - * @param vec the vector to set - * @param value - */ -static inline void vc_set(vc* vec,const c value){ - c_set(vec->_data,value,vec->size); -} - /** * Sets all values in a matrix to the value * @@ -164,7 +135,7 @@ static inline void dmat_set(dmat* mat,const d value){ } } } - +#define vd_set dmat_set /** * Sets all values in a matrix to the value @@ -180,49 +151,8 @@ static inline void cmat_set(cmat* mat,const c value){ } } } +#define vc_set cmat_set -/** - * Allocate data for a float vector. - * - * @param size Size of the vector - * - * @return vd with allocated data - */ -static inline vd vd_alloc(us size) { - vd result = { size, NULL,NULL}; - #ifdef LASP_DEBUG - result._data = (d*) a_malloc((size+1)*sizeof(d)); - result._data[size] = OVERFLOW_MAGIC_NUMBER; - #else - result._data = (d*) a_malloc(size*sizeof(d)); - #endif // LASP_DEBUG - result._foreign_data = false; - #ifdef LASP_DEBUG - vd_set(&result,NAN); - #endif // LASP_DEBUG - return result; -} -/** - * Allocate data for a complex vector. - * - * @param size Size of the vector - * - * @return vc with allocated data - */ -static inline vc vc_alloc(us size) { - vc result = { size, NULL, NULL}; - #ifdef LASP_DEBUG - result._data = (c*) a_malloc((size+1)*sizeof(c)); - result._data[size] = OVERFLOW_MAGIC_NUMBER; - #else - result._data = (c*) a_malloc(size*sizeof(c)); - #endif // LASP_DEBUG - result._foreign_data = false; - #ifdef LASP_DEBUG - vc_set(&result,NAN+I*NAN); - #endif // LASP_DEBUG - return result; -} /** * Allocate data for a matrix of floating points * @@ -250,7 +180,6 @@ static inline dmat dmat_alloc(us n_rows, return result; } - /** * Allocate a matrix of complex floating points * @@ -276,6 +205,30 @@ static inline cmat cmat_alloc(const us n_rows, #endif // LASP_DEBUG return result; } + +/** + * Allocate data for a float vector. + * + * @param size Size of the vector + * + * @return vd with allocated data + */ +static inline vd vd_alloc(us size) { + return dmat_alloc(size,1); +} + +/** + * Allocate data for a complex vector. + * + * @param size Size of the vector + * + * @return vc with allocated data + */ +static inline vc vc_alloc(us size) { + return cmat_alloc(size,1); +} + + /** * Creates a dmat from foreign data. Does not copy the data, but only * initializes the row pointers. Assumes column-major ordering for the @@ -288,32 +241,62 @@ static inline cmat cmat_alloc(const us n_rows, * * @return */ -static inline dmat dmat_foreign(const us n_rows, - const us n_cols, - d* data) { - - dbgassert(data,NULLPTRDEREF); - dmat result = {n_rows,n_cols,true,n_rows,data}; - return result; -} -static inline dmat dmat_foreign_vd(vd* vector) { - dbgassert(vector,NULLPTRDEREF); - dmat result = {vector->size,1,true,vector->size,vector->_data}; +static inline dmat dmat_foreign(dmat* other) { + dbgassert(other,NULLPTRDEREF); + dmat result = {other->n_rows, + other->n_cols, + true, + other->stride, + other->_data}; return result; } /** - * Create vd from foreign data + * Create a dmat from foreign data. Assumes the stride of the data is + * n_rows. * - * @param size Size of the vector - * @param data Pointer to data + * @param n_rows Number of rows + * @param n_cols Number of columns + * @param data Pointer to data storage * - * @return + * @return dmat */ -static inline vd vd_foreign(const us size,d* data) { +static inline dmat dmat_foreign_data(us n_rows, + us n_cols, + d* data, + bool own_data) { + dbgassert(data,NULLPTRDEREF); - vd result = {size,true,data}; + dmat result = {n_rows, + n_cols, + !own_data, + n_rows, + data}; return result; } +/** + * Create a cmat from foreign data. Assumes the stride of the data is + * n_rows. + * + * @param n_rows Number of rows + * @param n_cols Number of columns + * @param data Pointer to data storage + * + * @return dmat + */ +static inline cmat cmat_foreign_data(us n_rows, + us n_cols, + c* data, + bool own_data) { + + dbgassert(data,NULLPTRDEREF); + cmat result = {n_rows, + n_cols, + !own_data, + n_rows, + data}; + return result; +} + /** * Creates a cmat from foreign data. Does not copy the data, but only * initializes the row pointers. Assumes column-major ordering for the @@ -326,34 +309,18 @@ static inline vd vd_foreign(const us size,d* data) { * * @return */ -static inline cmat cmat_foreign(const us n_rows, - const us n_cols, - c* data) { - dbgassert(data,NULLPTRDEREF); - cmat result = {n_rows,n_cols,true,n_rows,data}; +static inline cmat cmat_foreign(cmat* other) { + dbgassert(other,NULLPTRDEREF); + cmat result = {other->n_rows, + other->n_cols, + true, + other->stride, + other->_data}; return result; } -/** - * Free's data of a vector. Is safe to run on sub-vecs as well, to - * make API consistent. (Only free's data if data pointer is set) - * - * @param f Vector to free - */ -static inline void vd_free(vd* f) { - dbgassert(f,NULLPTRDEREF); - if(!(f->_foreign_data)) a_free(f->_data); -} -/** - * Free's data of a vector. Is safe to run on sub-vecs as well, to - * make API consistent. (Only free's data if data pointer is set) - * - * @param f Vector to free - */ -static inline void vc_free(vc* f) { - dbgassert(f,NULLPTRDEREF); - if(!(f->_foreign_data)) a_free(f->_data); -} + + /** * Free's data of dmat. Safe to run on sub-matrices as well. * @@ -363,6 +330,8 @@ static inline void dmat_free(dmat* m) { dbgassert(m,NULLPTRDEREF); if(!(m->_foreign_data)) a_free(m->_data); } +#define vd_free dmat_free + /** * Free's data of dmat. Safe to run on sub-matrices as well. * @@ -372,7 +341,7 @@ static inline void cmat_free(cmat* m) { dbgassert(m,NULLPTRDEREF); if(!(m->_foreign_data)) a_free(m->_data); } - +#define vc_free cmat_free /** * Copy some rows from one matrix to another @@ -457,40 +426,6 @@ static inline cmat cmat_submat(cmat* parent, return result; } -/** - * Copy contents of one vector to another - * - * @param to : Vector to write to - * @param from : Vector to read from - */ -static inline void vd_copy(vd* to,const vd* from) { - dbgassert(to && from,NULLPTRDEREF); - dbgassert(to->size==from->size,SIZEINEQUAL); - d_copy(to->_data,from->_data,to->size,1,1); -} -static inline void vd_copy_rows(vd* to, - const vd* from, - const us startrow_to, - const us startrow_from, - const us nrows) { - dbgassert(to && from,NULLPTRDEREF); - dbgassert(startrow_from+nrows <= from->size,OUTOFBOUNDSMATR); - dbgassert(startrow_to+nrows <= to->size,OUTOFBOUNDSMATR); - d_copy(&to->_data[startrow_to], - &from->_data[startrow_from], - nrows,1,1); -} -/** - * Copy contents of one vector to another - * - * @param to : Vector to write to - * @param from : Vector to read from - */ -static inline void vc_copy(vc* to,const vc* from) { - dbgassert(to && from,NULLPTRDEREF); - dbgassert(to->size==from->size,SIZEINEQUAL); - c_copy(to->_data,from->_data,to->size); -} /** * Copy contents of one matrix to another. Sizes should be equal * @@ -524,6 +459,31 @@ static inline void cmat_copy(cmat* to,const cmat* from) { } } +/** + * Copy contents of one vector to another + * + * @param to : Vector to write to + * @param from : Vector to read from + */ +static inline void vd_copy(vd* to,const vd* from) { + dbgassert(to && from,NULLPTRDEREF); + assert_vx(to); + assert_vx(from); + dmat_copy(to,from); +} +/** + * Copy contents of one vector to another + * + * @param to : Vector to write to + * @param from : Vector to read from + */ +static inline void vc_copy(vc* to,const vc* from) { + dbgassert(to && from,NULLPTRDEREF); + assert_vx(to); + assert_vx(from); + cmat_copy(to,from); +} + /** * Get a reference to a column of a matrix as a vector @@ -534,7 +494,7 @@ static inline void cmat_copy(cmat* to,const cmat* from) { * @return vector with reference to column */ static inline vd dmat_column(dmat* x,us col) { - vd res = {x->n_rows,true,getdmatval(x,0,col)}; + vd res = {x->n_rows,1,true,x->stride,getdmatval(x,0,col)}; return res; } @@ -547,7 +507,7 @@ static inline vd dmat_column(dmat* x,us col) { * @return vector with reference to column */ static inline vc cmat_column(cmat* x,us col) { - vc res = {x->n_rows,true,getcmatval(x,0,col)}; + vc res = {x->n_rows,1,true,x->stride,getcmatval(x,0,col)}; return res; } @@ -557,11 +517,14 @@ static inline vc cmat_column(cmat* x,us col) { * @param a * @param b */ -static inline void vc_conj(vc* a,const vc* b) { +static inline void cmat_conj(cmat* a,const cmat* b) { fsTRACE(15); dbgassert(a && b,NULLPTRDEREF); - dbgassert(a->size == b->size,SIZEINEQUAL); - carray_conj(a->_data,b->_data,a->size); + dbgassert(a->n_cols == b->n_cols,SIZEINEQUAL); + dbgassert(a->n_rows == b->n_rows,SIZEINEQUAL); + for(us col=0;coln_cols;col++) { + carray_conj(getcmatval(a,0,col),getcmatval(b,0,col),a->n_rows); + } feTRACE(15); } @@ -570,7 +533,7 @@ static inline void vc_conj(vc* a,const vc* b) { * * @param x */ -static inline void cmat_conj(cmat* x) { +static inline void cmat_conj_inplace(cmat* x) { dbgassert(x,NULLPTRDEREF); for(us col=0;coln_cols;col++) { c_conj_inplace(getcmatval(x,0,col),x->n_rows); @@ -580,14 +543,15 @@ static inline void cmat_conj(cmat* x) { #ifdef LASP_DEBUG void print_cmat(const cmat* m); -void print_vc(const vc* m); -void print_vd(const vd* m); void print_dmat(const dmat* m); +#define print_vc print_cmat +#define print_vd print_dmat + #else #define print_cmat(m) #define print_vc(m) #define print_dmat(m) #endif -#endif // LASP_MATH_H +#endif // LASP_MAT_H ////////////////////////////////////////////////////////////////////// diff --git a/lasp/c/lasp_ps.c b/lasp/c/lasp_ps.c index ead4cda..f78ee68 100644 --- a/lasp/c/lasp_ps.c +++ b/lasp/c/lasp_ps.c @@ -148,7 +148,7 @@ void PowerSpectra_compute(const PowerSpectra* ps, vc j_vec = cmat_column(&fft_work,j); /* Compute the conjugate of spectra j */ - vc_conj(&j_vec_conj,&j_vec); + cmat_conj(&j_vec_conj,&j_vec); check_overflow_xmat(fft_work); diff --git a/lasp/c/lasp_python.h b/lasp/c/lasp_python.h index 8eb64e0..286a286 100644 --- a/lasp/c/lasp_python.h +++ b/lasp/c/lasp_python.h @@ -8,7 +8,6 @@ #pragma once #ifndef LASP_PYTHON_H #define LASP_PYTHON_H -#define TRACERPLUS (-5) #include #ifdef LASP_DOUBLE_PRECISION #define LASP_NUMPY_FLOAT_TYPE NPY_FLOAT64 diff --git a/lasp/c/lasp_signals.h b/lasp/c/lasp_signals.h index 7f5ea78..71205f0 100644 --- a/lasp/c/lasp_signals.h +++ b/lasp/c/lasp_signals.h @@ -8,7 +8,7 @@ #pragma once #ifndef LASP_SIGNALS_H #define LASP_SIGNALS_H -#include "lasp_math.h" +#include "lasp_mat.h" /** * Compute the signal power, that is \f$ \frac{1}{N} \sum_{i=0}^{N-1} @@ -19,10 +19,10 @@ */ static inline d signal_power(vd* signal) { d res = 0; - for(us i=0;isize;i++) { + for(us i=0;in_rows;i++) { res+= d_pow(*getvdval(signal,i),2); } - res /= signal->size; + res /= signal->n_rows; return res; } diff --git a/lasp/c/lasp_window.c b/lasp/c/lasp_window.c index e3769f4..00cd797 100644 --- a/lasp/c/lasp_window.c +++ b/lasp/c/lasp_window.c @@ -59,7 +59,10 @@ static d bartlett(us n,us N) { } int window_create(const WindowType wintype,vd* result,d* win_pow) { fsTRACE(15); - us nfft = result->size; + dbgassert(result && win_pow,NULLPTRDEREF); + assert_vx(result); + + us nfft = result->n_rows; d (*win_fun)(us,us); switch (wintype) { case Hann: { diff --git a/lasp/c/lasp_window.h b/lasp/c/lasp_window.h index 3a507fa..4ef81ad 100644 --- a/lasp/c/lasp_window.h +++ b/lasp/c/lasp_window.h @@ -8,7 +8,7 @@ #pragma once #ifndef LASP_WINDOW_H #define LASP_WINDOW_H -#include "lasp_math.h" +#include "lasp_mat.h" typedef enum { Hann = 0, diff --git a/lasp/config.pxi.in b/lasp/config.pxi.in index 1af75a4..70d1e49 100644 --- a/lasp/config.pxi.in +++ b/lasp/config.pxi.in @@ -29,7 +29,7 @@ cdef extern from "lasp_tracer.h": void fsTRACE(int) void feTRACE(int) void clearScreen() -cdef extern from "lasp_math.h": +cdef extern from "lasp_mat.h": ctypedef struct dmat: us n_cols us n_rows @@ -37,16 +37,17 @@ cdef extern from "lasp_math.h": bint _foreign_data ctypedef struct cmat: pass - ctypedef struct vd: - pass - ctypedef struct vc: - pass - dmat dmat_foreign(us n_rows, - us n_cols, - d* data) - cmat cmat_foreign(us n_rows, - us n_cols, - c* data) + ctypedef cmat vc + ctypedef dmat vd + + dmat dmat_foreign_data(us n_rows, + us n_cols, + d* data, + bint own_data) + cmat cmat_foreign_data(us n_rows, + us n_cols, + c* data, + bint own_data) cmat cmat_alloc(us n_rows,us n_cols) dmat dmat_alloc(us n_rows,us n_cols) diff --git a/lasp/wrappers.pyx b/lasp/wrappers.pyx index b84c8b6..ad26785 100644 --- a/lasp/wrappers.pyx +++ b/lasp/wrappers.pyx @@ -50,13 +50,15 @@ cdef class Fft: # result[:,:] = np.nan+1j*np.nan cdef c[::1,:] result_view = result - cdef cmat r = cmat_foreign(result.shape[0], - result.shape[1], - &result_view[0,0]) + cdef cmat r = cmat_foreign_data(result.shape[0], + result.shape[1], + &result_view[0,0], + False) - cdef dmat t = dmat_foreign(timedata.shape[0], - timedata.shape[1], - &timedata[0,0]) + cdef dmat t = dmat_foreign_data(timedata.shape[0], + timedata.shape[1], + &timedata[0,0], + False) Fft_fft(self._fft,&t,&r) @@ -74,18 +76,20 @@ cdef class Fft: # result[:,:] = np.nan+1j*np.nan - cdef cmat f = cmat_foreign(freqdata.shape[0], - freqdata.shape[1], - &freqdata[0,0]) + cdef cmat f = cmat_foreign_data(freqdata.shape[0], + freqdata.shape[1], + &freqdata[0,0], + False) timedata = np.empty((nfft,nchannels), dtype=NUMPY_FLOAT_TYPE, order='F') cdef d[::1,:] timedata_view = timedata - cdef dmat t = dmat_foreign(timedata.shape[0], - timedata.shape[1], - &timedata_view[0,0]) + cdef dmat t = dmat_foreign_data(timedata.shape[0], + timedata.shape[1], + &timedata_view[0,0], + False) Fft_ifft(self._fft,&f,&t) @@ -142,9 +146,10 @@ cdef class PowerSpectra: cmat result_mat - td = dmat_foreign(nfft, - nchannels, - &timedata[0,0]) + td = dmat_foreign_data(nfft, + nchannels, + &timedata[0,0], + False) # The array here is created in such a way that the strides # increase with increasing dimension. This is required for @@ -159,9 +164,10 @@ cdef class PowerSpectra: cdef c[::1,:,:] result_view = result - result_mat = cmat_foreign(nfft//2+1, - nchannels*nchannels, - &result_view[0,0,0]) + result_mat = cmat_foreign_data(nfft//2+1, + nchannels*nchannels, + &result_view[0,0,0], + False) @@ -208,8 +214,8 @@ cdef class AvPowerSpectra: cdef vd weighting_vd cdef vd* weighting_ptr = NULL if(weighting.size != 0): - weighting_vd = vd_foreign(weighting.size, - &weighting[0]) + weighting_vd = dmat_foreign_data(weighting.n_rows,1, + &weighting[0],False) weighting_ptr = &weighting_vd self.aps = AvPowerSpectra_alloc(nfft, @@ -242,9 +248,10 @@ cdef class AvPowerSpectra: if nchannels != self.nchannels: raise RuntimeError('Invalid number of channels') - td = dmat_foreign(nsamples, - nchannels, - &timedata[0,0]) + td = dmat_foreign_data(nsamples, + nchannels, + &timedata[0,0], + False) result_ptr = AvPowerSpectra_addTimeData(self.aps, &td) @@ -260,9 +267,10 @@ cdef class AvPowerSpectra: dtype = NUMPY_COMPLEX_TYPE, order='F') cdef c[::1,:,:] result_view = result - cdef cmat res = cmat_foreign(self.nfft//2+1, - nchannels*nchannels, - &result_view[0,0,0]) + cdef cmat res = cmat_foreign_data(self.nfft//2+1, + nchannels*nchannels, + &result_view[0,0,0], + True) # Copy result cmat_copy(&res,result_ptr) @@ -282,7 +290,11 @@ cdef class FilterBank: cdef: c_FilterBank* fb def __cinit__(self,d[::1,:] h, us nfft): - cdef dmat hmat = dmat_foreign(h.shape[0],h.shape[1],&h[0,0]) + cdef dmat hmat = dmat_foreign_data(h.shape[0], + h.shape[1], + &h[0,0], + True) + self.fb = FilterBank_create(&hmat,nfft) dmat_free(&hmat) if not self.fb: @@ -293,7 +305,9 @@ cdef class FilterBank: FilterBank_free(self.fb) def filter_(self,d[:] input_): - cdef vd input_vd = vd_foreign(input_.size,&input_[0]) + cdef dmat input_vd = dmat_foreign_data(input_.shape[0],1, + &input_[0],False) + cdef dmat output = FilterBank_filter(self.fb,&input_vd) # Steal the pointer from output @@ -326,9 +340,10 @@ cdef class Decimator: def decimate(self,d[::1,:] samples): assert samples.shape[1] == self.nchannels,'Invalid number of channels' - cdef dmat d_samples = dmat_foreign(samples.shape[0], - samples.shape[1], - &samples[0,0]) + cdef dmat d_samples = dmat_foreign_data(samples.shape[0], + samples.shape[1], + &samples[0,0], + False) cdef dmat res = Decimator_decimate(self.dec,&d_samples) result = dmat_to_ndarray(&res,True) @@ -337,4 +352,4 @@ cdef class Decimator: def __dealloc__(self): if self.dec != NULL: - Decimator_free(self.dec) \ No newline at end of file + Decimator_free(self.dec) diff --git a/test/test_bf.c b/test/test_bf.c index a97cb8f..5375747 100644 --- a/test/test_bf.c +++ b/test/test_bf.c @@ -5,7 +5,7 @@ // Description: // ////////////////////////////////////////////////////////////////////// -#include "lasp_math.h" +#include "lasp_mat.h" int main() { diff --git a/test/test_math.c b/test/test_math.c index a04eceb..404eeac 100644 --- a/test/test_math.c +++ b/test/test_math.c @@ -7,7 +7,7 @@ ////////////////////////////////////////////////////////////////////// #include "lasp_fft.h" -#include "lasp_math.h" +#include "lasp_mat.h" #include "lasp_tracer.h" #include "lasp_alg.h"