Merged master

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

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; */

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);
+    assert_vx(a);
-    return d_dot(getvdval(a,0),getvdval(b,0),a->size);
+    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);
+    assert_equalsize(a,b);
-    dbgassert(b->size==a->size,SIZEINEQUAL);
     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);
+    assert_equalsize(a,b);
-    dbgassert(b->size==a->size,SIZEINEQUAL);
+    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);

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) {
 

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 {

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;channelno<nchannels;channelno++) {
         dec->fbs[channelno] = FilterBank_create(&h,DEC_FFT_LEN);

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;
 

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;
 

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);

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.

lasp/c/lasp_filterbank.c

@@ -68,12 +68,6 @@ FilterBank* FilterBank_create(const dmat* h,
 
     dmat_free(&temp);
 
-    /* Push initial output of zeros to the output fifo */
-    dmat initial_output = dmat_alloc(fb->P_m_1,nfilters);
-    dmat_set(&initial_output,0);
-    dFifo_push(fb->output_fifo,&initial_output);
-    dmat_free(&initial_output);
-    
     feTRACE(15);
     return fb;
 }
@@ -100,9 +94,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,x);
-    dFifo_push(fb->input_fifo,&samples);
-    dmat_free(&samples);
 
     dmat input_block = dmat_alloc(nfft,1);
 

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;
 
 /** 

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 <math.h>
@@ -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;row<m->size;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;row<m->size;row++){
-	
-        d rval = *getvdval(m,row);
-        indent_trace();
-        printf("%c%2.2e ",rval< 0 ? '\r': ' ',rval);
-        printf("\n");
-    }
-    feTRACE(50);
-}
 
 #endif
 

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
+#ifndef LASP_MAT_H
-#define LASP_MATH_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)                              \
+typedef dmat vd;
-    dbgassert((((us) index) <= (vec)->size),OUTOFBOUNDSVEC);  \
+typedef cmat vc;
-    (vec)->_data[index] = val;
 
+#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,
+static inline dmat dmat_foreign(dmat* other) {
-                                const us n_cols,
+    dbgassert(other,NULLPTRDEREF);
-                                d* data) {
+    dmat result = {other->n_rows,
-
+                   other->n_cols,
-    dbgassert(data,NULLPTRDEREF);
+                   true,
-    dmat result = {n_rows,n_cols,true,n_rows,data};
+                   other->stride,
-    return result;
+                   other->_data};
-}
-static inline dmat dmat_foreign_vd(vd* vector) {
-    dbgassert(vector,NULLPTRDEREF);
-    dmat result = {vector->size,1,true,vector->size,vector->_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 n_rows Number of rows
- * @param data Pointer to data
+ * @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,
+static inline cmat cmat_foreign(cmat* other) {
-                                const us n_cols,
+    dbgassert(other,NULLPTRDEREF);
-                                c* data) {
+    cmat result = {other->n_rows,
-    dbgassert(data,NULLPTRDEREF);
+                   other->n_cols,
-    cmat result = {n_rows,n_cols,true,n_rows,data};
+                   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);
+    dbgassert(a->n_cols == b->n_cols,SIZEINEQUAL);
-    carray_conj(a->_data,b->_data,a->size);
+    dbgassert(a->n_rows == b->n_rows,SIZEINEQUAL);
+    for(us col=0;col<a->n_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;col<x->n_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
 //////////////////////////////////////////////////////////////////////

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);
 

lasp/c/lasp_python.h

@@ -8,7 +8,6 @@
 #pragma once
 #ifndef LASP_PYTHON_H
 #define LASP_PYTHON_H
-#define TRACERPLUS (-5)
 #include <numpy/ndarrayobject.h>
 #ifdef LASP_DOUBLE_PRECISION
 #define LASP_NUMPY_FLOAT_TYPE NPY_FLOAT64

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;i<signal->size;i++) {
+    for(us i=0;i<signal->n_rows;i++) {
         res+= d_pow(*getvdval(signal,i),2);
     }
-    res /= signal->size;
+    res /= signal->n_rows;
     return res;
 }
 

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: {

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,

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:
+    ctypedef cmat vc
-        pass
+    ctypedef dmat vd
-    ctypedef struct vc:
+
-        pass
+    dmat dmat_foreign_data(us n_rows,
-    dmat dmat_foreign(us n_rows,
+                           us n_cols,
-                      us n_cols,
+                           d* data,
-                      d* data)
+                           bint own_data)
-    cmat cmat_foreign(us n_rows,
+    cmat cmat_foreign_data(us n_rows,
-                      us n_cols,
+                           us n_cols,
-                      c* data)
+                           c* data,
+                           bint own_data)
 
     cmat cmat_alloc(us n_rows,us n_cols)
     dmat dmat_alloc(us n_rows,us n_cols)

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],
+        cdef cmat r = cmat_foreign_data(result.shape[0],
-                                   result.shape[1],
+                                        result.shape[1],
-                                   &result_view[0,0])
+                                        &result_view[0,0],
+                                        False)
 
-        cdef dmat t = dmat_foreign(timedata.shape[0],
+        cdef dmat t = dmat_foreign_data(timedata.shape[0],
-                                   timedata.shape[1],
+                                        timedata.shape[1],
-                                   &timedata[0,0])
+                                        &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],
+        cdef cmat f = cmat_foreign_data(freqdata.shape[0],
-                                   freqdata.shape[1],
+                                        freqdata.shape[1],
-                                   &freqdata[0,0])
+                                        &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],
+        cdef dmat t = dmat_foreign_data(timedata.shape[0],
-                                   timedata.shape[1],
+                                        timedata.shape[1],
-                                   &timedata_view[0,0])
+                                        &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,
+        td = dmat_foreign_data(nfft,
-                          nchannels,
+                               nchannels,
-                          &timedata[0,0])
+                               &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,
+        result_mat = cmat_foreign_data(nfft//2+1,
-                                  nchannels*nchannels,
+                                       nchannels*nchannels,
-                                  &result_view[0,0,0])
+                                       &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_vd = dmat_foreign_data(weighting.n_rows,1,
-                                      &weighting[0])
+                                             &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,
+        td = dmat_foreign_data(nsamples,
-                          nchannels,
+                               nchannels,
-                          &timedata[0,0])
+                               &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,
+        cdef cmat res = cmat_foreign_data(self.nfft//2+1,
-                                     nchannels*nchannels,
+                                          nchannels*nchannels,
-                                     &result_view[0,0,0])
+                                          &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],
+        cdef dmat d_samples = dmat_foreign_data(samples.shape[0],
-                                           samples.shape[1],
+                                                samples.shape[1],
-                                           &samples[0,0])
+                                                &samples[0,0],
+                                                False)
         
         cdef dmat res = Decimator_decimate(self.dec,&d_samples)
         result = dmat_to_ndarray(&res,True)
@@ -339,7 +354,6 @@ cdef class Decimator:
         if self.dec != NULL:
             Decimator_free(self.dec)
 
-
 cdef extern from "lasp_sp_lowpass.h":
     ctypedef struct c_SPLowpass "SPLowpass"
     c_SPLowpass* SPLowpass_create(d fs,d tau)

test/test_bf.c

@@ -5,7 +5,7 @@
 // Description:
 //
 //////////////////////////////////////////////////////////////////////
-#include "lasp_math.h"
+#include "lasp_mat.h"
 
 int main() {
 

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"