Changed internals of xmat to use strides instead of column pointers. Some other bugfixes

This commit is contained in:
Anne de Jong 2018-02-12 20:44:07 +01:00 committed by Anne de Jong
parent b31866974e
commit 54e9eeca37
9 changed files with 167 additions and 218 deletions

View File

@ -16,9 +16,7 @@ void cmv_dot(const cmat* A,const vc* restrict x,vc* restrict b){
dbgassert(A->n_cols == x->size,SIZEINEQUAL); dbgassert(A->n_cols == x->size,SIZEINEQUAL);
#if ASCEE_USE_BLAS == 1 #if ASCEE_USE_BLAS == 1
dbgassert(false,"Untested function");
dbgassert(A->_data,"Matrix-vector product only works for allocated data matrices");
/* typedef enum CBLAS_ORDER {CblasRowMajor=101, CblasColMajor=102} CBLAS_ORDER; */ /* typedef enum CBLAS_ORDER {CblasRowMajor=101, CblasColMajor=102} CBLAS_ORDER; */
/* typedef enum CBLAS_TRANSPOSE {CblasNoTrans=111, CblasTrans=112, CblasConjTrans=113, CblasConjNoTrans=114} CBLAS_TRANSPOSE; */ /* typedef enum CBLAS_TRANSPOSE {CblasNoTrans=111, CblasTrans=112, CblasConjTrans=113, CblasConjNoTrans=114} CBLAS_TRANSPOSE; */
/* /*
@ -44,10 +42,10 @@ void cmv_dot(const cmat* A,const vc* restrict x,vc* restrict b){
(d*) &alpha, /* alpha */ (d*) &alpha, /* alpha */
(d*) A->_data, /* A */ (d*) A->_data, /* A */
A->n_rows, /* lda */ A->n_rows, /* lda */
(d*) x->ptr, /* */ (d*) x->_data, /* */
1, 1,
(d*) &beta, /* beta */ (d*) &beta, /* beta */
(d*) b->ptr, (d*) b->_data,
1); 1);

View File

@ -18,8 +18,9 @@
* @return dot product as float * @return dot product as float
*/ */
static inline d vd_dot(const vd * a,const vd* b) { static inline d vd_dot(const vd * a,const vd* b) {
dbgassert(a && b,NULLPTRDEREF);
dbgassert(a->size == b->size,SIZEINEQUAL); dbgassert(a->size == b->size,SIZEINEQUAL);
return d_dot(a->ptr,b->ptr,a->size); return d_dot(getvdval(a,0),getvdval(b,0),a->size);
} }
/** /**
@ -30,14 +31,10 @@ static inline d vd_dot(const vd * a,const vd* b) {
*/ */
static inline void dmat_scale(dmat* y,const c fac){ static inline void dmat_scale(dmat* y,const c fac){
dbgassert(y,NULLPTRDEREF); dbgassert(y,NULLPTRDEREF);
if(likely(y->_data)) { for(us col=0;col<y->n_cols;col++) {
d_scale(y->_data,fac,y->n_cols*y->n_rows); d_scale(getdmatval(y,0,col),fac,y->n_rows);
}
else {
for(us col=0;col<y->n_cols;col++) {
d_scale(y->col_ptrs[col],fac,y->n_rows);
}
} }
} }
/** /**
* y = fac * y * y = fac * y
@ -48,13 +45,8 @@ static inline void dmat_scale(dmat* y,const c fac){
static inline void cmat_scale(cmat* y,const c fac){ static inline void cmat_scale(cmat* y,const c fac){
dbgassert(y,NULLPTRDEREF); dbgassert(y,NULLPTRDEREF);
dbgassert(y,NULLPTRDEREF); dbgassert(y,NULLPTRDEREF);
if(likely(y->_data)) { for(us col=0;col<y->n_cols;col++) {
c_scale(y->_data,fac,y->n_cols*y->n_rows); c_scale(getcmatval(y,0,col),fac,y->n_rows);
}
else {
for(us col=0;col<y->n_cols;col++) {
c_scale(y->col_ptrs[col],fac,y->n_rows);
}
} }
} }
@ -69,13 +61,9 @@ static inline void dmat_add_dmat(dmat* x,dmat* y,d fac) {
dbgassert(x && y,NULLPTRDEREF); dbgassert(x && y,NULLPTRDEREF);
dbgassert(x->n_cols == y->n_cols,SIZEINEQUAL); dbgassert(x->n_cols == y->n_cols,SIZEINEQUAL);
dbgassert(x->n_rows == y->n_rows,SIZEINEQUAL); dbgassert(x->n_rows == y->n_rows,SIZEINEQUAL);
if(likely(x->_data && y->_data)) { for(us col=0;col<y->n_cols;col++) {
d_add_to(x->_data,y->_data,fac,x->n_cols*x->n_rows); d_add_to(getdmatval(x,0,col),
} getdmatval(y,0,col),fac,x->n_rows);
else {
for(us col=0;col<y->n_cols;col++) {
d_add_to(x->col_ptrs[col],y->col_ptrs[col],fac,x->n_rows);
}
} }
} }
/** /**
@ -86,17 +74,12 @@ static inline void dmat_add_dmat(dmat* x,dmat* y,d fac) {
* @param[in] fac * @param[in] fac
*/ */
static inline void cmat_add_cmat(cmat* x,cmat* y,c fac) { static inline void cmat_add_cmat(cmat* x,cmat* y,c fac) {
// dbgassert(x && y,NULLPTRDEREF); dbgassert(x && y,NULLPTRDEREF);
// dbgassert(x->n_cols == y->n_cols,SIZEINEQUAL); dbgassert(x->n_cols == y->n_cols,SIZEINEQUAL);
// dbgassert(x->n_rows == y->n_rows,SIZEINEQUAL); dbgassert(x->n_rows == y->n_rows,SIZEINEQUAL);
// if(likely(x->data && y->data)) {
// TRACE(15,"Scale whole");
// c_add_to(x->data,y->data,fac,x->n_cols*x->n_rows);
// }
// else {
for(us col=0;col<y->n_cols;col++) { for(us col=0;col<y->n_cols;col++) {
TRACE(15,"Scale columns"); c_add_to(getcmatval(x,0,col),
c_add_to(x->col_ptrs[col],y->col_ptrs[col],fac,x->n_rows); getcmatval(y,0,col),fac,x->n_rows);
} }
} }
@ -112,7 +95,9 @@ static inline void vd_elem_prod(vd* result,const vd* a,const vd* b) {
dbgassert(result && a && b,NULLPTRDEREF); dbgassert(result && a && b,NULLPTRDEREF);
dbgassert(result->size==a->size,SIZEINEQUAL); dbgassert(result->size==a->size,SIZEINEQUAL);
dbgassert(b->size==a->size,SIZEINEQUAL); dbgassert(b->size==a->size,SIZEINEQUAL);
d_elem_prod_d(result->ptr,a->ptr,b->ptr,a->size); d_elem_prod_d(getvdval(result,0),
getvdval(a,0),
getvdval(b,0),a->size);
} }
/** /**
* Compute the element-wise (Hadamard) product of a and b, and store * Compute the element-wise (Hadamard) product of a and b, and store
@ -127,7 +112,11 @@ static inline void vc_hadamard(vc* result,const vc* a,const vc* b) {
dbgassert(result && a && b,NULLPTRDEREF); dbgassert(result && a && b,NULLPTRDEREF);
dbgassert(result->size==a->size,SIZEINEQUAL); dbgassert(result->size==a->size,SIZEINEQUAL);
dbgassert(b->size==a->size,SIZEINEQUAL); dbgassert(b->size==a->size,SIZEINEQUAL);
c_hadamard(result->ptr,a->ptr,b->ptr,a->size); c_hadamard(getvcval(result,0),
getvcval(a,0),
getvcval(b,0),
a->size);
check_overflow_vx(*result); check_overflow_vx(*result);
check_overflow_vx(*a); check_overflow_vx(*a);
check_overflow_vx(*b); check_overflow_vx(*b);
@ -153,9 +142,9 @@ static inline void cmat_hadamard(cmat* result,
dbgassert(b->n_cols==a->n_cols,SIZEINEQUAL); dbgassert(b->n_cols==a->n_cols,SIZEINEQUAL);
for(us col=0;col<result->n_cols;col++) { for(us col=0;col<result->n_cols;col++) {
c_hadamard(result->col_ptrs[col], c_hadamard(getcmatval(result,0,col),
a->col_ptrs[col], getcmatval(a,0,col),
b->col_ptrs[col], getcmatval(b,0,col),
a->n_rows); a->n_rows);
} }

View File

@ -9,12 +9,18 @@
#ifndef ASCEE_ALLOC_H #ifndef ASCEE_ALLOC_H
#define ASCEE_ALLOC_H #define ASCEE_ALLOC_H
#include <malloc.h> #include <malloc.h>
#include "ascee_tracer.h"
/** /**
* Reserved words for memory allocation. Can be changed to something * Reserved words for memory allocation. Can be changed to something
* else when required. For example for debugging purposes. * else when required. For example for debugging purposes.
*/ */
#define a_malloc malloc static inline void* a_malloc(size_t nbytes) {
void* ptr = malloc(nbytes);
if(!ptr) {
FATAL("Memory allocation failed. Exiting");
}
return ptr;
}
#define a_free free #define a_free free
#define a_realloc realloc #define a_realloc realloc

View File

@ -72,7 +72,7 @@ void print_vd(const vd* m) {
iVARTRACE(20,m->size); iVARTRACE(20,m->size);
for(row=0;row<m->size;row++){ for(row=0;row<m->size;row++){
d rval = m->ptr[row]; d rval = *getvdval(m,row);
indent_trace(); indent_trace();
printf("%c%2.2e ",rval< 0 ? '\r': ' ',rval); printf("%c%2.2e ",rval< 0 ? '\r': ' ',rval);
printf("\n"); printf("\n");

View File

@ -16,8 +16,7 @@
/// Vector of floating point numbers /// Vector of floating point numbers
typedef struct { typedef struct {
us size; us size;
d* ptr; /**< This pointer points to the data bool _foreign_data;
of this vector */
d* _data; /**< Pointer set if data storage is d* _data; /**< Pointer set if data storage is
intern. If this is set to zero, the intern. If this is set to zero, the
vector is a sub-vector. */ vector is a sub-vector. */
@ -25,8 +24,7 @@ typedef struct {
/// Vector of complex floating point numbers /// Vector of complex floating point numbers
typedef struct { typedef struct {
us size; us size;
c* ptr; /**< This pointer points to the data bool _foreign_data;
of this vector */
c* _data; /**< Pointer set if data storage is c* _data; /**< Pointer set if data storage is
intern. If this is set to zero, the intern. If this is set to zero, the
vector is a sub-vector. */ vector is a sub-vector. */
@ -35,33 +33,101 @@ typedef struct {
typedef struct { typedef struct {
us n_rows; us n_rows;
us n_cols; us n_cols;
d** col_ptrs; bool _foreign_data;
us stride;
d* _data; d* _data;
} dmat; } dmat;
/// Dense matrix of complex floating point values /// Dense matrix of complex floating point values
typedef struct { typedef struct {
us n_rows; us n_rows;
us n_cols; us n_cols;
c** col_ptrs; bool _foreign_data;
us stride;
c* _data; c* _data;
} cmat; } cmat;
#define setvecval(vec,index,val) \
dbgassert((((us) index) <= (vec)->size),OUTOFBOUNDSVEC); \
(vec)->_data[index] = val;
#define setmatval(mat,row,col,val) \
dbgassert((((us) row) <= mat->n_rows),OUTOFBOUNDSMATR); \
dbgassert((((us) col) <= mat->n_cols),,OUTOFBOUNDSMATC); \
(mat)->data[(col)*(mat)->stride+(row)] = val;
/**
* Return pointer to a value from a vector
*
* @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
*
* @param mat The matrix
* @param row The row
* @param col The column
*/
static inline d* getdmatval(const dmat* mat,us row,us col){
dbgassert(mat,NULLPTRDEREF);
dbgassert(row < mat->n_rows,OUTOFBOUNDSMATR);
dbgassert(col < mat->n_cols,OUTOFBOUNDSMATC);
return &mat->_data[col*mat->stride+row];
}
/**
* Return a value from a matrix of complex floating points
*
* @param mat The matrix
* @param row The row
* @param col The column
*/
static inline c* getcmatval(const cmat* mat,const us row,const us col){
dbgassert(mat,NULLPTRDEREF);
dbgassert(row < mat->n_rows,OUTOFBOUNDSMATR);
dbgassert(col < mat->n_cols,OUTOFBOUNDSMATC);
return &mat->_data[col*mat->stride+row];
}
#ifdef ASCEE_DEBUG #ifdef ASCEE_DEBUG
#define OVERFLOW_MAGIC_NUMBER (-10e-45) #define OVERFLOW_MAGIC_NUMBER (-10e-45)
#define check_overflow_vx(vx) \ #define check_overflow_vx(vx) \
TRACE(15,"Checking overflow " #vx); \ TRACE(15,"Checking overflow " #vx); \
if((vx)._data) { \ if(!(vx)._foreign_data) { \
dbgassert((vx)._data[(vx).size] == OVERFLOW_MAGIC_NUMBER, \ dbgassert((vx)._data[(vx).size] == OVERFLOW_MAGIC_NUMBER, \
"Buffer overflow detected on" #vx ); \ "Buffer overflow detected on" #vx ); \
} } \
else { \
WARN("Cannot check overflow on foreign buffer"); \
}
#define check_overflow_xmat(xmat) \ #define check_overflow_xmat(xmat) \
TRACE(15,"Checking overflow " #xmat); \ TRACE(15,"Checking overflow " #xmat); \
if((xmat)._data) { \ if(!(xmat)._foreign_data) { \
for(us _overflow=0;_overflow<(xmat).n_cols;_overflow++) \ dbgassert((xmat)._data[((xmat).n_cols-1)*(xmat).stride+(xmat).n_rows] \
dbgassert((xmat)._data[(xmat).n_rows*(xmat).n_cols+ \ == OVERFLOW_MAGIC_NUMBER, \
_overflow] == OVERFLOW_MAGIC_NUMBER, \ "Buffer overflow detected on" #xmat ); \
"Buffer overflow detected on" #xmat ); \ } \
else { \
WARN("Cannot check overflow on foreign buffer"); \
} }
#else #else
@ -76,7 +142,7 @@ typedef struct {
* @param value * @param value
*/ */
static inline void vd_set(vd* vec, d value){ static inline void vd_set(vd* vec, d value){
d_set(vec->ptr,value,vec->size); d_set(vec->_data,value,vec->size);
} }
/** /**
@ -86,7 +152,7 @@ static inline void vd_set(vd* vec, d value){
* @param value * @param value
*/ */
static inline void vc_set(vc* vec,const c value){ static inline void vc_set(vc* vec,const c value){
c_set(vec->ptr,value,vec->size); c_set(vec->_data,value,vec->size);
} }
/** /**
@ -98,7 +164,7 @@ static inline void vc_set(vc* vec,const c value){
static inline void dmat_set(dmat* mat,const d value){ static inline void dmat_set(dmat* mat,const d value){
dbgassert(mat,NULLPTRDEREF); dbgassert(mat,NULLPTRDEREF);
for(us col=0;col<mat->n_cols;col++) { for(us col=0;col<mat->n_cols;col++) {
d_set(mat->col_ptrs[col],value,mat->n_rows); d_set(getdmatval(mat,0,col),value,mat->n_rows);
} }
} }
@ -112,7 +178,7 @@ static inline void dmat_set(dmat* mat,const d value){
static inline void cmat_set(cmat* mat,const c value){ static inline void cmat_set(cmat* mat,const c value){
dbgassert(mat,NULLPTRDEREF); dbgassert(mat,NULLPTRDEREF);
for(us col=0;col<mat->n_cols;col++) { for(us col=0;col<mat->n_cols;col++) {
c_set(mat->col_ptrs[col],value,mat->n_rows); c_set(getcmatval(mat,0,col),value,mat->n_rows);
} }
} }
@ -126,14 +192,12 @@ static inline void cmat_set(cmat* mat,const c value){
static inline vd vd_alloc(us size) { static inline vd vd_alloc(us size) {
vd result = { size, NULL,NULL}; vd result = { size, NULL,NULL};
#ifdef ASCEE_DEBUG #ifdef ASCEE_DEBUG
result._data = (d*) a_malloc((size+1)*sizeof(d)); result._data = (d*) a_malloc((size+1)*sizeof(d));
result._data[size] = OVERFLOW_MAGIC_NUMBER; result._data[size] = OVERFLOW_MAGIC_NUMBER;
#else #else
result._data = (d*) a_malloc(size*sizeof(d)); result._data = (d*) a_malloc(size*sizeof(d));
#endif // ASCEE_DEBUG #endif // ASCEE_DEBUG
result._foreign_data = false;
result.ptr = result._data;
dbgassert(result._data,ALLOCFAILED);
#ifdef ASCEE_DEBUG #ifdef ASCEE_DEBUG
vd_set(&result,NAN); vd_set(&result,NAN);
#endif // ASCEE_DEBUG #endif // ASCEE_DEBUG
@ -154,8 +218,7 @@ static inline vc vc_alloc(us size) {
#else #else
result._data = (c*) a_malloc(size*sizeof(c)); result._data = (c*) a_malloc(size*sizeof(c));
#endif // ASCEE_DEBUG #endif // ASCEE_DEBUG
dbgassert(result._data,ALLOCFAILED); result._foreign_data = false;
result.ptr = result._data;
#ifdef ASCEE_DEBUG #ifdef ASCEE_DEBUG
vc_set(&result,NAN+I*NAN); vc_set(&result,NAN+I*NAN);
#endif // ASCEE_DEBUG #endif // ASCEE_DEBUG
@ -172,26 +235,15 @@ static inline vc vc_alloc(us size) {
*/ */
static inline dmat dmat_alloc(us n_rows, static inline dmat dmat_alloc(us n_rows,
us n_cols) { us n_cols) {
dmat result = { n_rows, n_cols, NULL, NULL}; dmat result = { n_rows, n_cols, false, n_rows, NULL};
/**
* Here storage is allocated for both the data, as well as the
* column pointers. In debug mode, extra memory is allocated to
* check for possible buffer overflows.
*/
#ifdef ASCEE_DEBUG #ifdef ASCEE_DEBUG
result._data = (d*) a_malloc((n_rows*n_cols+n_cols)*sizeof(d)); result._data = (d*) a_malloc((n_rows*n_cols+1)*sizeof(d));
for(us i=0;i<n_cols;i++) result._data[n_rows*n_cols] = OVERFLOW_MAGIC_NUMBER;
result._data[n_rows*n_cols+i] = OVERFLOW_MAGIC_NUMBER;
#else #else
result._data = (d*) a_malloc((n_rows*n_cols+1)*sizeof(d)); result._data = (d*) a_malloc((n_rows*n_cols)*sizeof(d));
#endif // ASCEE_DEBUG #endif // ASCEE_DEBUG
dbgassert(result._data,ALLOCFAILED);
result.col_ptrs = a_malloc(n_cols*sizeof(d*));
for(us col=0;col<n_cols;col++) {
result.col_ptrs[col] = &result._data[n_rows*col];
}
#ifdef ASCEE_DEBUG #ifdef ASCEE_DEBUG
dmat_set(&result,NAN); dmat_set(&result,NAN);
#endif // ASCEE_DEBUG #endif // ASCEE_DEBUG
@ -211,23 +263,14 @@ static inline dmat dmat_alloc(us n_rows,
*/ */
static inline cmat cmat_alloc(const us n_rows, static inline cmat cmat_alloc(const us n_rows,
const us n_cols) { const us n_cols) {
cmat result = { n_rows, n_cols, NULL, NULL}; cmat result = { n_rows, n_cols, false, n_rows, NULL};
/**
* Here storage is allocated for both the data, as well as the
* column pointers. In debug mode, extra memory is allocated to
* check for possible buffer overflows.
*/
#ifdef ASCEE_DEBUG #ifdef ASCEE_DEBUG
result._data = (c*) a_malloc((n_rows*n_cols+n_cols)*sizeof(c)); result._data = (c*) a_malloc((n_rows*n_cols+1)*sizeof(c));
for(us i=0;i<n_cols;i++) result._data[n_rows*n_cols] = OVERFLOW_MAGIC_NUMBER;
result._data[n_rows*n_cols+i] = OVERFLOW_MAGIC_NUMBER;
#else #else
result._data = (c*) a_malloc((n_rows*n_cols+1)*sizeof(c)); result._data = (c*) a_malloc((n_rows*n_cols)*sizeof(c));
#endif // ASCEE_DEBUG #endif // ASCEE_DEBUG
result.col_ptrs = a_malloc(n_cols*sizeof(c*));
for(us col=0;col<n_cols;col++) {
result.col_ptrs[col] = &result._data[n_rows*col];
}
#ifdef ASCEE_DEBUG #ifdef ASCEE_DEBUG
cmat_set(&result,NAN+I*NAN); cmat_set(&result,NAN+I*NAN);
@ -238,7 +281,7 @@ static inline cmat cmat_alloc(const us n_rows,
* Creates a dmat from foreign data. Does not copy the data, but only * Creates a dmat from foreign data. Does not copy the data, but only
* initializes the row pointers. Assumes column-major ordering for the * initializes the row pointers. Assumes column-major ordering for the
* data. Please do not keep this one alive after the data has been * data. Please do not keep this one alive after the data has been
* destroyed. * destroyed. Assumes the column stride equals to n_rows.
* *
* @param n_rows Number of rows * @param n_rows Number of rows
* @param n_cols Number of columns * @param n_cols Number of columns
@ -251,21 +294,14 @@ static inline dmat dmat_foreign(const us n_rows,
d* data) { d* data) {
dbgassert(data,NULLPTRDEREF); dbgassert(data,NULLPTRDEREF);
dmat result = {n_rows,n_cols,NULL,NULL}; dmat result = {n_rows,n_cols,true,n_rows,data};
d** colptrs = malloc(sizeof(d*)*n_cols);
dbgassert(colptrs,ALLOCFAILED);
result.col_ptrs = colptrs;
for(us i=0;i<n_cols;i++) {
colptrs[i] = &data[i*n_rows];
}
return result; return result;
} }
/** /**
* Creates a cmat from foreign data. Does not copy the data, but only * Creates a cmat from foreign data. Does not copy the data, but only
* initializes the row pointers. Assumes column-major ordering for the * initializes the row pointers. Assumes column-major ordering for the
* data. Please do not keep this one alive after the data has been * data. Please do not keep this one alive after the data has been
* destroyed. * destroyed. Assumes the column stride equals to n_rows.
* *
* @param n_rows * @param n_rows
* @param n_cols * @param n_cols
@ -277,13 +313,7 @@ static inline cmat cmat_foreign(const us n_rows,
const us n_cols, const us n_cols,
c* data) { c* data) {
dbgassert(data,NULLPTRDEREF); dbgassert(data,NULLPTRDEREF);
cmat result = {n_rows,n_cols,NULL,NULL}; cmat result = {n_rows,n_cols,true,n_rows,data};
c** colptrs = malloc(sizeof(c*)*n_cols);
dbgassert(colptrs,ALLOCFAILED);
result.col_ptrs = colptrs;
for(us i=0;i<n_cols;i++) {
colptrs[i] = &data[i*n_rows];
}
return result; return result;
} }
@ -295,7 +325,7 @@ static inline cmat cmat_foreign(const us n_rows,
*/ */
static inline void vd_free(vd* f) { static inline void vd_free(vd* f) {
dbgassert(f,NULLPTRDEREF); dbgassert(f,NULLPTRDEREF);
if(likely(f->_data)) a_free(f->_data); if(!(f->_foreign_data)) a_free(f->_data);
} }
/** /**
* Free's data of a vector. Is safe to run on sub-vecs as well, to * Free's data of a vector. Is safe to run on sub-vecs as well, to
@ -305,7 +335,7 @@ static inline void vd_free(vd* f) {
*/ */
static inline void vc_free(vc* f) { static inline void vc_free(vc* f) {
dbgassert(f,NULLPTRDEREF); dbgassert(f,NULLPTRDEREF);
if(likely(f->_data)) a_free(f->_data); if(!(f->_foreign_data)) a_free(f->_data);
} }
/** /**
* Free's data of dmat. Safe to run on sub-matrices as well. * Free's data of dmat. Safe to run on sub-matrices as well.
@ -313,9 +343,8 @@ static inline void vc_free(vc* f) {
* @param m Matrix to free * @param m Matrix to free
*/ */
static inline void dmat_free(dmat* m) { static inline void dmat_free(dmat* m) {
if(m->_data) a_free(m->_data); dbgassert(m,NULLPTRDEREF);
dbgassert(m->col_ptrs,NULLPTRDEREF); if(!(m->_foreign_data)) a_free(m->_data);
a_free(m->col_ptrs);
} }
/** /**
* Free's data of dmat. Safe to run on sub-matrices as well. * Free's data of dmat. Safe to run on sub-matrices as well.
@ -323,70 +352,10 @@ static inline void dmat_free(dmat* m) {
* @param m Matrix to free * @param m Matrix to free
*/ */
static inline void cmat_free(cmat* m) { static inline void cmat_free(cmat* m) {
if(m->_data) a_free(m->_data); dbgassert(m,NULLPTRDEREF);
dbgassert(m->col_ptrs,NULLPTRDEREF); if(!(m->_foreign_data)) a_free(m->_data);
a_free(m->col_ptrs);
} }
#define setvecval(vec,index,val) \
dbgassert((((us) index) <= (vec)->size),OUTOFBOUNDSVEC); \
(vec)->ptr[index] = val;
#define setmatval(mat,row,col,val) \
dbgassert((((us) row) <= mat->n_rows),OUTOFBOUNDSMATR); \
dbgassert((((us) col) <= mat->n_cols),,OUTOFBOUNDSMATC); \
(mat)->data[(col)*(mat)->n_rows+(row)] = val;
/**
* Return pointer to a value from a vector
*
* @param mat The vector
* @param row The row
*/
static inline d* getvdval(const vd* vec,us row){
dbgassert(row < vec->size,OUTOFBOUNDSVEC);
return &vec->ptr[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->ptr[row];
}
/**
* Return a value from a matrix of floating points
*
* @param mat The matrix
* @param row The row
* @param col The column
*/
static inline d* getdmatval(const dmat* mat,us row,us col){
dbgassert(mat,NULLPTRDEREF);
dbgassert(row < mat->n_rows,OUTOFBOUNDSMATR);
dbgassert(col < mat->n_cols,OUTOFBOUNDSMATC);
return &mat->col_ptrs[col][row];
}
/**
* Return a value from a matrix of complex floating points
*
* @param mat The matrix
* @param row The row
* @param col The column
*/
static inline c* getcmatval(const cmat* mat,const us row,const us col){
dbgassert(mat,NULLPTRDEREF);
dbgassert(row < mat->n_rows,OUTOFBOUNDSMATR);
dbgassert(col < mat->n_cols,OUTOFBOUNDSMATC);
return &mat->col_ptrs[col][row];
}
/** /**
* Copy some rows from one matrix to another * Copy some rows from one matrix to another
@ -432,16 +401,9 @@ static inline dmat dmat_submat(const dmat* parent,
dbgassert(n_rows+startrow <= parent->n_rows,OUTOFBOUNDSMATR); dbgassert(n_rows+startrow <= parent->n_rows,OUTOFBOUNDSMATR);
dbgassert(n_cols+startcol <= parent->n_cols,OUTOFBOUNDSMATC); dbgassert(n_cols+startcol <= parent->n_cols,OUTOFBOUNDSMATC);
d** col_ptrs = malloc(sizeof(d*)*n_cols); dmat result = { n_rows,n_cols,true,n_rows-startrow,
dbgassert(col_ptrs,ALLOCFAILED); getdmatval(parent,startrow,startcol)};
for(us col=0;col<n_cols;col++) {
col_ptrs[col] = getdmatval(parent,
startrow,
startcol+col);
}
dmat result = { n_rows,n_cols,NULL,NULL};
result.col_ptrs = col_ptrs;
return result; return result;
} }
/** /**
@ -465,15 +427,9 @@ static inline cmat cmat_submat(cmat* parent,
dbgassert(n_rows+startrow <= parent->n_rows,OUTOFBOUNDSMATR); dbgassert(n_rows+startrow <= parent->n_rows,OUTOFBOUNDSMATR);
dbgassert(n_cols+startcol <= parent->n_cols,OUTOFBOUNDSMATC); dbgassert(n_cols+startcol <= parent->n_cols,OUTOFBOUNDSMATC);
c** col_ptrs = malloc(sizeof(c*)*n_cols); cmat result = { n_rows,n_cols,true,n_rows-startrow,
dbgassert(col_ptrs,ALLOCFAILED); getcmatval(parent,startrow,startcol)};
for(us col=0;col<n_cols;col++) {
col_ptrs[col] = getcmatval(parent,
startrow,
startcol+col);
}
cmat result = { n_rows,n_cols,col_ptrs,NULL};
return result; return result;
} }
@ -486,7 +442,7 @@ static inline cmat cmat_submat(cmat* parent,
static inline void vd_copy(vd* to,const vd* from) { static inline void vd_copy(vd* to,const vd* from) {
dbgassert(to && from,NULLPTRDEREF); dbgassert(to && from,NULLPTRDEREF);
dbgassert(to->size==from->size,SIZEINEQUAL); dbgassert(to->size==from->size,SIZEINEQUAL);
d_copy(to->ptr,from->ptr,to->size); d_copy(to->_data,from->_data,to->size);
} }
static inline void vd_copy_rows(vd* to, static inline void vd_copy_rows(vd* to,
const vd* from, const vd* from,
@ -496,8 +452,8 @@ static inline void vd_copy_rows(vd* to,
dbgassert(to && from,NULLPTRDEREF); dbgassert(to && from,NULLPTRDEREF);
dbgassert(startrow_from+nrows <= from->size,OUTOFBOUNDSMATR); dbgassert(startrow_from+nrows <= from->size,OUTOFBOUNDSMATR);
dbgassert(startrow_to+nrows <= to->size,OUTOFBOUNDSMATR); dbgassert(startrow_to+nrows <= to->size,OUTOFBOUNDSMATR);
d_copy(&to->ptr[startrow_to], d_copy(&to->_data[startrow_to],
&from->ptr[startrow_from], &from->_data[startrow_from],
nrows); nrows);
} }
/** /**
@ -509,7 +465,7 @@ static inline void vd_copy_rows(vd* to,
static inline void vc_copy(vc* to,const vc* from) { static inline void vc_copy(vc* to,const vc* from) {
dbgassert(to && from,NULLPTRDEREF); dbgassert(to && from,NULLPTRDEREF);
dbgassert(to->size==from->size,SIZEINEQUAL); dbgassert(to->size==from->size,SIZEINEQUAL);
c_copy(to->ptr,from->ptr,to->size); c_copy(to->_data,from->_data,to->size);
} }
/** /**
* Copy contents of one matrix to another. Sizes should be equal * Copy contents of one matrix to another. Sizes should be equal
@ -522,7 +478,9 @@ static inline void dmat_copy(dmat* to,const dmat* from) {
dbgassert(to->n_rows==from->n_rows,SIZEINEQUAL); dbgassert(to->n_rows==from->n_rows,SIZEINEQUAL);
dbgassert(to->n_cols==from->n_cols,SIZEINEQUAL); dbgassert(to->n_cols==from->n_cols,SIZEINEQUAL);
for(us col=0;col<to->n_cols;col++) { for(us col=0;col<to->n_cols;col++) {
d_copy(to->col_ptrs[col],from->col_ptrs[col],to->n_rows); d_copy(getdmatval(to,0,col),
getdmatval(from,0,col),
to->n_rows);
} }
} }
/** /**
@ -536,7 +494,9 @@ static inline void cmat_copy(cmat* to,const cmat* from) {
dbgassert(to->n_rows==from->n_rows,SIZEINEQUAL); dbgassert(to->n_rows==from->n_rows,SIZEINEQUAL);
dbgassert(to->n_cols==from->n_cols,SIZEINEQUAL); dbgassert(to->n_cols==from->n_cols,SIZEINEQUAL);
for(us col=0;col<to->n_cols;col++) { for(us col=0;col<to->n_cols;col++) {
c_copy(to->col_ptrs[col],from->col_ptrs[col],to->n_rows); c_copy(getcmatval(to,0,col),
getcmatval(from,0,col),
to->n_rows);
} }
} }
@ -550,10 +510,7 @@ static inline void cmat_copy(cmat* to,const cmat* from) {
* @return vector with reference to column * @return vector with reference to column
*/ */
static inline vd dmat_column(dmat* x,us col) { static inline vd dmat_column(dmat* x,us col) {
vd res; vd res = {x->n_rows,true,getdmatval(x,0,col)};
res.size = x->n_rows;
res.ptr = getdmatval(x,0,col);
res._data = NULL;
return res; return res;
} }
@ -566,10 +523,7 @@ static inline vd dmat_column(dmat* x,us col) {
* @return vector with reference to column * @return vector with reference to column
*/ */
static inline vc cmat_column(cmat* x,us col) { static inline vc cmat_column(cmat* x,us col) {
vc res; vc res = {x->n_rows,true,getcmatval(x,0,col)};
res.size = x->n_rows;
res.ptr = getcmatval(x,0,col);
res._data = NULL;
return res; return res;
} }
@ -583,7 +537,7 @@ static inline void vc_conj(vc* a,const vc* b) {
fsTRACE(15); fsTRACE(15);
dbgassert(a && b,NULLPTRDEREF); dbgassert(a && b,NULLPTRDEREF);
dbgassert(a->size == b->size,SIZEINEQUAL); dbgassert(a->size == b->size,SIZEINEQUAL);
carray_conj(a->ptr,b->ptr,a->size); carray_conj(a->_data,b->_data,a->size);
feTRACE(15); feTRACE(15);
} }
@ -595,7 +549,7 @@ static inline void vc_conj(vc* a,const vc* b) {
static inline void cmat_conj(cmat* x) { static inline void cmat_conj(cmat* x) {
dbgassert(x,NULLPTRDEREF); dbgassert(x,NULLPTRDEREF);
for(us col=0;col<x->n_cols;col++) { for(us col=0;col<x->n_cols;col++) {
c_conj_inplace(x->col_ptrs[col],x->n_rows); c_conj_inplace(getcmatval(x,0,col),x->n_rows);
} }
} }

View File

@ -10,6 +10,7 @@
#define ASCEE_TRACER_H #define ASCEE_TRACER_H
#include <string.h> #include <string.h>
#include <stdio.h> #include <stdio.h>
#include <stdlib.h>
static inline void clearScreen() { static inline void clearScreen() {
printf("\033c\n"); printf("\033c\n");

View File

@ -33,7 +33,7 @@ Fft* Fft_alloc(const us nfft) {
fft->fft_work = vd_alloc(2*nfft+15); fft->fft_work = vd_alloc(2*nfft+15);
fft->fft_result = vd_alloc(nfft); fft->fft_result = vd_alloc(nfft);
npy_rffti(nfft,fft->fft_work.ptr); npy_rffti(nfft,getvdval(&fft->fft_work,0));
check_overflow_vx(fft->fft_work); check_overflow_vx(fft->fft_work);
/* print_vd(&fft->fft_work); */ /* print_vd(&fft->fft_work); */
@ -70,7 +70,8 @@ void Fft_fft_single(const Fft* fft,const vd* timedata,vc* result) {
vd_copy(&fft_result,timedata); vd_copy(&fft_result,timedata);
/* Perform fft */ /* Perform fft */
npy_rfftf(nfft,fft_result.ptr,fft->fft_work.ptr); npy_rfftf(nfft,getvdval(&fft_result,0),
getvdval(&fft->fft_work,0));
/* Fftpack stores the data a bit strange, the resulting array /* Fftpack stores the data a bit strange, the resulting array
* has the DC value at 0,the first cosine at 1, the first sine * has the DC value at 0,the first cosine at 1, the first sine

View File

@ -35,7 +35,7 @@ int main() {
vc res2 = vc_alloc(3); vc res2 = vc_alloc(3);
c_hadamard(res2.ptr,vc1.ptr,vc2.ptr,3); c_hadamard(res2._data,vc1._data,vc2._data,3);
print_vc(&res2); print_vc(&res2);

View File

@ -31,7 +31,7 @@ int main() {
vc c1 = vc_alloc(5); vc c1 = vc_alloc(5);
/* vc_set(&c1,10); */ /* vc_set(&c1,10); */
/* c_add_to(c1.ptr,a.ptr,1,3); */ /* c_add_to(c1.ptr,a.ptr,1,3); */
c_hadamard(c1.ptr,a,b.ptr,5); c_hadamard(c1._data,a,b._data,5);