smath.h File Reference

Macros
#define	AOS_LIB_TYPE
#define	AOS_SMAT(A)   s##A

Typedefs
typedef float(*	sminsearch_fun) (const float *x, void *info)

Functions
smat *	snew (long nx, long ny)
smat *	snew_file (long nx, long ny, const char *keywords, const char *format,...)
smat *	snew_do (long nx, long ny, float *p, mem_t *mem)
smat *	smat_cast (const void *A)
int	sinit (smat **A, long nx, long ny)
void	sfree_content (smat *A)
void	sfree_do (smat *A)
smat *	sref (const smat *in)
smat *	sref_reshape (const smat *in, long nx, long ny)
smat *	srefcols (const smat *in, long icol, long ncol)
void	scols (smat *out, const smat *in, long icol, long ncol)
int	sresize (smat *A, long nx, long ny)
smat *	ssub (const smat *in, long sx, long nx, long sy, long ny)
smat *	scat (const smat *in1, const smat *in2, int dim)
smat *	sdup (const smat *in)
void	szero (smat *A)
int	szerocol (smat *A, int icol)
uint32_t	shash (const smat *A, uint32_t key)
void	scp (smat **out0, const smat *in)
smat *	strans (const smat *A)
void	sset (smat *A, const float val)
void	sshow (const smat *A, const char *format,...)
void	svecperm (smat *out, const smat *in, const long *perm)
void	svecpermi (smat *out, const smat *in, const long *perm)
int	sflip (smat *A, int axis)
float	svecsum (const float *__restrict p, long np)
float	ssum (const smat *A)
void	snormalize_sum (smat *A, float sum)
void	snormalize_sumabs (smat *A, float sum)
float	smean (const smat *A)
float	strace (const smat *A)
void	ssort (smat *A, int ascend)
void	svecmaxmin (const float *__restrict p, long N, float *max, float *min)
float	svecdot (const float *__restrict p1, const float *__restrict p2, const float *__restrict p3, long n)
void	svecnormalize_sum (float *__restrict p, long nloc, float norm)
void	svecnormalize_sumabs (float *__restrict p, long nloc, float norm)
float	svecmaxabs (const float *__restrict p, long n)
void	smaxmin (const smat *A, float *max, float *min)
float	smax (const smat *A)
float	smaxabs (const smat *A)
float	smin (const smat *A)
float	ssumabs (const smat *in)
float	ssumsq (const smat *in)
float	ssumdiffsq (const smat *A, const smat *B)
void	sfftshift (smat *A)
int	scpcorner2center (smat *A, const smat *B)
int	sshift (smat **B0, const smat *A, int sx, int sy)
scell *	scell_cast (const cell *A)
scell *	scellnew2 (const scell *A)
scell *	scellnew3 (long nx, long ny, long *nnx, long *nny)
scell *	scellnew_same (long nx, long ny, long mx, long my)
scell *	scellnew_file (long nx, long ny, long *nnx, long *nny, const char *keywords, const char *format,...)
scell *	scellnewsame_file (long nx, long ny, long mx, long my, const char *keywords, const char *format,...)
void	scellzero (scell *dc)
void	scellset (scell *dc, float alpha)
scell *	scelltrans (const scell *A)
scell *	scellref (const scell *in)
scell *	scelldup (const scell *in)
void	scellcp (scell **out0, const scell *in)
scell *	scellreduce (const scell *A, int dim)
scell *	scellcat (const scell *A, const scell *B, int dim)
void	scellcat2 (scell **A, const scell *B, int dim)
scell *	scellcat_each (const scell *A, const scell *B, int dim)
scell *	s2cellref (const smat *A, long *dims, long ndim)
void	s2cell (scell **B, const smat *A, const scell *ref)
smat *	scell_col (scell *input, long icol)
float	scellsum (const scell *A)
smat *	scellsum_each (const scell *A)
uint32_t	scellhash (const scell *A, uint32_t key)
int	sisnan (const smat *A)
float	snorm (const smat *in)
float	sstd (const smat *in)
void	srandu (smat *A, const float mean, rand_t *rstat)
void	srandn (smat *A, const float sigma, rand_t *rstat)
void	sscale (smat *A, float w)
float	sdot (const smat *A, const smat *B)
float	swdot (const float *a, const smat *w, const float *b)
void	scwm (smat *A, const smat *B)
void	scwm2 (smat *A, const smat *B1, float wt1, const smat *B2, float wt2)
void	scwm3 (smat *__restrict A, const smat *__restrict W, const smat *__restrict B1, const float wt1, const smat *__restrict B2, const float wt2)
void	scwmcol (smat *__restrict A, const smat *__restrict B)
void	scwmrow (smat *__restrict A, const smat *__restrict B)
void	scwdiv (smat *B, const smat *A, float value)
void	smulvec (float *__restrict y, const smat *__restrict A, const float *__restrict x, const float alpha)
void	smm (smat **C0, const float beta, const smat *A, const smat *B, const char trans[2], const float alpha)
void	sinvspd_inplace (smat *A)
smat *	sinvspd (const smat *A)
void	sinv_inplace (smat *A)
smat *	sinv (const smat *A)
smat *	schol (const smat *A)
smat *	smcc (const smat *A, const smat *wt)
smat *	simcc (const smat *A, const smat *wt)
smat *	stmcc (const smat *A, const smat *wt)
smat *	spinv2 (const smat *A, const_anyarray W, float thres)
smat *	spinv (const smat *A, const_anyarray W)
float	sdiff (const smat *A, const smat *B)
int	scircle (smat *A, float cx, float cy, float dx, float dy, float r, float val)
int	srectangle (smat *A, float cx, float cy, float rx, float ry, float theta, float val)
void	srotvec (smat *A, const float theta)
void	srotvect (smat *A, const float theta)
void	srotvecnn (smat **B0, const smat *A, float theta)
void	smulvec3 (float *y, const smat *A, const float *x)
void	scorr (smat **corr, const smat *A, const smat *B)
void	spara3 (float *grad, const smat *corr)
void	scog (float *grad, const smat *i0, float offsetx, float offsety, float thres, float bkgrnd, float flux)
void	sshift2center (smat *A, float offsetx, float offsety)
int	sclip (smat *A, float min, float max)
void	sgramschmidt (smat *Mod, float *amp)
void	smuldiag (smat *A, const smat *s)
void	scwpow (smat *A, float power)
void	scwexp (smat *A, float alpha)
void	scwpow_thres (smat *A, float power, float thres)
void	ssvd (smat **U, smat **Sdiag, smat **VT, const smat *A)
void	ssvd_cache (smat **U, smat **Sdiag, smat **VT, const smat *A)
void	ssvd_pow (smat *A, float power, float thres)
void	sexpm (smat **out, float alpha, const smat *A, float beta)
void	spolyval (smat *A, smat *p)
void	saddI (smat *A, float val)
void	sadd (smat **B0, float bc, const smat *A, const float ac)
void	sadd_relax (smat **B0, float bc, const smat *A, const float ac)
void	saddcol (smat *B, long icol, float bc, const smat *A, const float ac)
void	sadds (smat *A, const float ac)
smat *	slogspace (float emin, float emax, long n)
smat *	slinspace (float min, float dx, long n)
smat *	sinterp1 (const smat *xin, const smat *yin, const smat *xnew, float y0)
smat *	sinterp1_2 (const smat *xyin, const smat *xnew)
smat *	sinterp1linear (const smat *xin, const smat *yin, const smat *xnew, float y0)
smat *	sinterp1log (const smat *xin, const smat *yin, const smat *xnew, float y0)
void	sblend (smat *__restrict A, smat *__restrict B, int overlap)
void	shistfill (smat **out, const smat *A, float center, float spacing, int n)
smat *	sspline_prep (smat *x, smat *y)
smat *	sspline_eval (smat *coeff, smat *x, smat *xnew)
smat *	sspline (smat *x, smat *y, smat *xnew)
void	scwlog10 (smat *A)
void	scwlog (smat *A)
void	sembed (smat *__restrict A, const smat *__restrict B, const float theta)
void	sembedd (smat *__restrict A, const smat *__restrict B, const float theta)
float	sfwhm (smat *A)
void	sgauss_fit (float *mr, float *ma, float *mb, float *theta, smat *A, float thres)
float	sgauss_width (smat *A, float thres)
float	sfwhm_gauss (smat *A)
smat *	senc (const smat *A, const smat *dvec, int type, int nthread)
int	sminsearch (float *x, int nmod, float ftol, int nmax, sminsearch_fun fun, void *info)
void	sbessik (float x, float xnu, float *ri, float *rk, float *rip, float *rkp)
float	strapz (const smat *x, const smat *y)
float	scellnorm (const scell *in)
void	scellscale (anyarray A, float w)
void	scelldropempty (scell **A0, int dim)
float	scelldot (const scell *A, const scell *B)
void	scellcwm (scell *B, const scell *A)
scell *	scellinvspd (scell *A)
scell *	scellinv (scell *A)
scell *	scellinvspd_each (scell *A)
scell *	scellpinv (const scell *A, const_anyarray W, float thres)
void	scellsvd_pow (scell *A, float power, float thres)
void	scellcwpow (scell *A, float power)
void	scelldropzero (scell *B, float thres)
float	scelldiff (const scell *A, const scell *B)
int	scellclip (scell *Ac, float min, float max)
scell *	scellsub (const scell *in, long sx, long nx, long sy, long ny)
scell *	sbspline_prep (smat *x, smat *y, smat *z)
smat *	sbspline_eval (scell *coeff, smat *x, smat *y, smat *xnew, smat *ynew)
void	smaprot (anyarray A_, real theta)
void	swritedata (file_t *fp, const smat *A, long ncol)
smat *	sreaddata (file_t *fp, header_t *header)
smat *	snew_mmap (long nx, long ny, const char *keywords, const char *format,...)
scell *	scellnew_mmap (long nx, long ny, long *nnx, long *nny, const char *keywords, const char *format,...)
scell *	scellnewsame_mmap (long nx, long ny, long mx, long my, const char *keywords, const char *format,...)
smat *	sread_mmap (const char *format,...)
scell *	scellread_mmap (const char *format,...)
ssp *	sspnew (long nx, long ny, long nzmax)
ssp *	sspref (ssp *A)
ssp *	sspdup (const ssp *A)
ssp *	sspnew2 (const ssp *A)
ssp *	sspnewrandu (long nx, long ny, const float mean, float fill, rand_t *rstat)
void	sspsetnzmax (ssp *sp, long nzmax)
ssp *	ssp_cast (const cell *A)
void	sspfree_do (ssp *sp)
void	sspdisp (const ssp *sp)
int	sspcheck (const ssp *sp)
void	sspscale (ssp *A, const float beta)
void	sspscalex (ssp *A, const smat *xs)
void	sspscaley (ssp *A, const smat *ys)
sspcell *	sspcellref (const sspcell *in)
sspcell *	sspcell_cast (const cell *A)
void	sspcellscale (sspcell *A, const float beta)
ssp *	sspnewdiag (long N, float *vec, float alpha)
smat *	sspdiag (const ssp *A)
void	sspmuldiag (ssp *__restrict A, const float *w, float alpha)
void	sspmv (smat *y, const ssp *A, const smat *__restrict x, char trans, float alpha)
void	sspmulcreal (float *__restrict y, const ssp *A, const fcomplex *__restrict x, float alpha)
void	sspmm (smat **yout, const ssp *A, const smat *x, const char trans[2], const float alpha)
void	smulsp (smat **yout, const smat *x, const ssp *A, const char trans[2], const float alpha)
float	sspwdot (const smat *y, const ssp *A, const smat *x)
float	sspcellwdot (const scell *y, const sspcell *A, const scell *x)
void	scellmm (panyarray C0, const_anyarray A, const_anyarray B, const char trans[2], const float alpha)
scell *	scellmm2 (const scell *A, const scell *B, const char trans[2])
ssp *	s2sp (smat *A, float thres)
void	sspfull (smat **out0, const ssp *A, const char trans, const float alpha)
void	sspcellfull (scell **out0, const sspcell *A, const char trans, const float alpha)
ssp *	sspadd2 (const ssp *A, float a, const ssp *B, float b)
void	sspadd (ssp **A0, float alpha, const ssp *B, float beta)
void	scelladd (panyarray pA, float ac, const_anyarray B, float bc)
void	scelladdsp (scell **pA, float ac, const sspcell *B, float bc)
void	sspcelladd (sspcell **pA, float ac, const sspcell *B, float bc)
void	sspaddI (ssp *A0, float alpha)
void	scelladdI (anyarray A, float alpha)
ssp *	ssptrans (const ssp *A)
void	sspconj (ssp *)
ssp *	sspmulsp (const ssp *A, const ssp *B, const char trans[2])
void	sspmulsp2 (ssp **C0, const ssp *A, const ssp *B, const char trans[2], const float scale)
sspcell *	sspcelltrans (const sspcell *spc)
ssp *	sspcat (const ssp *A, const ssp *B, int type)
ssp *	sspcell2sp (const sspcell *A)
smat *	sspsum (const ssp *A, int col)
smat *	sspsumabs (const ssp *A, int col)
void	sspclean (ssp *A)
void	sspdroptol (ssp *A, float thres)
void	sspcelldroptol (sspcell *A, float thres)
lmat *	sspdropemptycol (ssp *A)
void	sspsort (ssp *A)
void	sspcellsort (sspcell *A)
void	sspsym (ssp **A)
void	sspcellsym (sspcell **A)
ssp *	sspconvolvop (smat *A)
ssp *	sspperm (ssp *A, int reverse, long *pcol, long *prow)
ssp *	sspinvbdiag (const ssp *A, long bs)
scell *	sspblockextract (const ssp *A, long bs)
smat *	scell2m (const_anyarray A)
void	sspwritedata (file_t *fp, const ssp *sp)
ssp *	sspreaddata (file_t *fp, header_t *header)

Function Documentation

snew()

smat* snew	(	long	nx,
		long	ny
	)

Create a new SX(mat). Initialized to zero.

snew_file()

smat* snew_file	(	long	nx,
		long	ny,
		const char *	keywords,
		const char *	format,
			...
	)

Calls SX(new) with a filename to be saved to.

smat_cast()

smat* smat_cast

(

const void *

A

)

check and cast an object to matrix

sinit()

int sinit	(	smat **	A,
		long	nx,
		long	ny
	)

check the size of matrix if exist. Otherwise create it. content is not zeroed if already exists.

• empty matrix is initialized.
• existing matrix is not resized.

sfree_content()

void sfree_content

(

smat *

A

)

free content of a matrix object.

sfree_do()

void sfree_do

(

smat *

A

)

free a matrix object.

sref()

smat* sref

(

const smat *

in

)

creat a reference to an existing matrix. header is duplicated if exists.

sref_reshape()

smat* sref_reshape	(	const smat *	in,
		long	nx,
		long	ny
	)

create an new reference another with different shape.

• total number of elements must be preserved
• set ny to 1 or 0 to return a vector.

srefcols()

smat* srefcols	(	const smat *	in,
		long	icol,
		long	ncol
	)

creat a new matrix referencing columns in existing matrix. reference counted.

scols()

void scols	(	smat *	out,
		const smat *	in,
		long	icol,
		long	ncol
	)

Override a stack matrix struct with pointers to columns of another matrix. Does not add reference count the original data.

sresize()

int sresize	(	smat *	A,
		long	nx,
		long	ny
	)

Resize a matrix by adding or removing columns or rows. Data is kept whenever possible.

• 0 in dimension preserves original value.

ssub()

smat* ssub	(	const smat *	in,
		long	sx,
		long	nx,
		long	sy,
		long	ny
	)

Create a new sub matrix of nx*ny starting from(sx,sy) and copy the data.

scat()

smat* scat	(	const smat *	in1,
		const smat *	in2,
		int	dim
	)

Concatenate two matrixes into 1 along "dim" (1 for x, 2 for y)

sdup()

smat* sdup

(

const smat *

in

)

duplicate a matrix

szero()

void szero

(

smat *

A

)

Set elements to zero.

szerocol()

int szerocol	(	smat *	A,
		int	icol
	)

Set column elements to zero.

shash()

uint32_t shash	(	const smat *	A,
		uint32_t	key
	)

Compute the hashlittle

scp()

void scp	(	smat **	out0,
		const smat *	in
	)

copy the values from one matrix to another.

strans()

smat* strans

(

const smat *

A

)

transpose a matrix

sshow()

void sshow	(	const smat *	A,
		const char *	format,
			...
	)

display a matrix.

svecperm()

void svecperm	(	smat *	out,
		const smat *	in,
		const long *	perm
	)

Permute the vector so that out(:)=in(p);

svecpermi()

void svecpermi	(	smat *	out,
		const smat *	in,
		const long *	perm
	)

Reverse permute the vector so that out(p)=in(:);

sflip()

int sflip	(	smat *	A,
		int	axis
	)

Flip the matrix along the set axis. 0 to flip both, 1 along x, 2 along y.

ssum()

float ssum

(

const smat *

A

)

Calculate sum of all the elements in A.

smean()

float smean

(

const smat *

A

)

Calculate average of all the elements in A.

strace()

float strace

(

const smat *

A

)

compute the trace (sum of diagonal elements)

ssort()

void ssort	(	smat *	A,
		int	ascend
	)

Sort all columns of A, in ascending order if ascend is non zero, otherwise in descending order.

smax()

float smax

(

const smat *

A

)

find the maximum value

smaxabs()

float smaxabs

(

const smat *

A

)

find the maximum of abs of all elements

smin()

float smin

(

const smat *

A

)

find the minimum value

ssumabs()

float ssumabs

(

const smat *

A

)

compute the sum of abs of all elements

ssumsq()

float ssumsq

(

const smat *

A

)

compute the sum of A.*A

ssumdiffsq()

float ssumdiffsq	(	const smat *	A,
		const smat *	B
	)

compute the sum of (A-B)^2

sfftshift()

void sfftshift

(

smat *

A

)

shift frequency components by n/2

scpcorner2center()

int scpcorner2center	(	smat *	A,
		const smat *	B
	)

reorder B and embed/crop into center of A .

4 * * 3
* * * *
* * * *
2 * * 1

to

1 2
3 4

sshift()

int sshift	(	smat **	B0,
		const smat *	A,
		int	sx,
		int	sy
	)

cyclic shift A by nx and ny to B.

4   3     1   2

2   1 to  3   4

scell_cast()

scell* scell_cast

(

const cell *

A

)

cast a cell object to actual cell after checking.

scellnew2()

scell* scellnew2

(

const scell *

A

)

create an new cell similar to A in shape. When a cell is empty, it is created with an emtry matrix and cannot be overriden.

scellnew3()

scell* scellnew3	(	long	nx,
		long	ny,
		long *	nnx,
		long *	nny
	)

Create an new cell with matrix specified. Each block is stored continuously in memory.

scellnew_same()

scell* scellnew_same	(	long	nx,
		long	ny,
		long	mx,
		long	my
	)

Create an new cell with matrix specified. Each block is stored continuously in memory.

scellnew_file()

scell* scellnew_file	(	long	nx,
		long	ny,
		long *	nnx,
		long *	nny,
		const char *	keywords,
		const char *	format,
			...
	)

Calls cellnew3 with a filename to be saved to.

scellnewsame_file()

scell* scellnewsame_file	(	long	nx,
		long	ny,
		long	mx,
		long	my,
		const char *	keywords,
		const char *	format,
			...
	)

Calls cellnew_same with a filename to be saved to.

scellzero()

void scellzero

(

scell *

dc

)

Setting all elements of a cell to zero.

scelltrans()

scell* scelltrans

(

const scell *

A

)

transpose a cell object

scellref()

scell* scellref

(

const scell *

in

)

creat a cell reference an existing cell by referencing the elements.

scelldup()

scell* scelldup

(

const scell *

in

)

duplicate a SX(cell) object. Not implemented for lmat or imat as SX(cellcp) is defined in spmm.c

scellcp()

void scellcp	(	scell **	A_,
		const scell *	B_
	)

Copy B to A.

Takes parameters of matrix, sparse matrix, or cell array of them.

scellreduce()

scell* scellreduce	(	const scell *	A,
		int	dim
	)

reduce nx*ny cell matrix to 1*ny if dim=1 and nx*1 if dim=2 by merging the cells.

scellcat()

scell* scellcat	(	const scell *	A,
		const scell *	B,
		int	dim
	)

concatenate two cell matrices along dimenstion 'dim'.

scellcat2()

void scellcat2	(	scell **	A,
		const scell *	B,
		int	dim
	)

concatenate two cell matrices along dimenstion 'dim'.

scellcat_each()

scell* scellcat_each	(	const scell *	A,
		const scell *	B,
		int	dim
	)

concatenate coresponding elements of each SX(cell). They must have the same shape.

s2cellref()

scell* s2cellref	(	const smat *	A,
		long *	dims,
		long	ndim
	)

convert a vector to cell using dimensions specified in dims. Reference the vector

s2cell()

void s2cell	(	scell **	B,
		const smat *	A,
		const scell *	ref
	)

make A a cell array using shape information from ref if *B is NULL.

Notice that empty cells may be replaced by zeros if it is within valid column or row.

scell_col()

smat* scell_col	(	scell *	input,
		long	icol
	)

input is nsa*ncol cell. each cell has npix=nx*ny elements. Extract icol of cell as npix*nsa array.

scellsum()

float scellsum

(

const scell *

A

)

create sum of all the elements in A.

scellsum_each()

smat* scellsum_each

(

const scell *

A

)

return sum of each cell in a SX(mat)

scellhash()

uint32_t scellhash	(	const scell *	A,
		uint32_t	key
	)

Compute the hashlittle

sisnan()

int sisnan

(

const smat *

A

)

Check for NaN in elements

snorm()

float snorm

(

const smat *

A

)

compute the norm(2) of A

sstd()

float sstd

(

const smat *

A

)

compute the standard deviation of A

sdot()

float sdot	(	const smat *	A,
		const smat *	B
	)

compute the dot product of A and B. (inner product)

scwmcol()

void scwmcol	(	smat *__restrict	A,
		const smat *__restrict	B
	)

Compute component wise multiply A=A.*(B1*wt1+B2*wt2) component-wise multiply of three matrices. A=A.*W.*(B1*wt1+B2*wt2) Component-wise multiply each column of A with B A(:,i)=A(:,i).*B;

scwmrow()

void scwmrow	(	smat *__restrict	A,
		const smat *__restrict	B
	)

component-wise multiply of columns of A with combination of B1 and B2: A(:,i)=A(:,i)*(B1*wt1+B2*wt2); component wise multiply of 2d complex matrix A,W and 1d vector B. A(:,i)=A(:,i).*W(:,i).*B; Component wise multiply each row of A with B. A(i,:)=A(i,:)*B

sinvspd_inplace()

void sinvspd_inplace

(

smat *

A

)

inplace invert a small square SPD matrix using lapack dposv_, usually (A'*w*A). by solving Ax=I; copy x to A. dposv_ modifies A also. be careful

sinvspd()

smat* sinvspd

(

const smat *

A

)

out of place version of dinvspd_inplace

sinv_inplace()

void sinv_inplace

(

smat *

A

)

inplace invert a general square matrix using lapack dgesv_

sinv()

smat* sinv

(

const smat *

A

)

out of place version of dinv

schol()

smat* schol

(

const smat *

A

)

Compute the cholesky decomposition of a symmetric semi-definit dense matrix.

smcc()

smat* smcc	(	const smat *	A,
		const smat *	wt
	)

compute (A'*W*A); where diag(W)=wt

simcc()

smat* simcc	(	const smat *	A,
		const smat *	wt
	)

compute inv(dmcc(A, wt))

stmcc()

smat* stmcc	(	const smat *	A,
		const smat *	wt
	)

compute (A*W*A'); where diag(W)=wt

spinv()

smat* spinv	(	const smat *	A,
		const_anyarray	W
	)

A convenient wrapper.

sdiff()

float sdiff	(	const smat *	A,
		const smat *	B
	)

compute the relative difference betwee two vectors. sqrt(||A-B||/||A||) using norm2. for debugging purpose.

scorr()

void scorr	(	smat **	pout,
		const smat *	A,
		const smat *	B
	)

Compute the correlation matrix.

sblend()

void sblend	(	smat *__restrict	A,
		smat *__restrict	B,
		int	overlap
	)

blend B into center of A with width of overlap. The center (size is B->nx-overlap, B->ny-overlap) of A is replaced by center of B . The overlapping area is blended

sspline_prep()

smat* sspline_prep	(	smat *	x,
		smat *	y
	)

1D Cubic spline interpolation preparation. if x has only 1 column: x is the coordinate. y is the function value. if x has two columns: first column is the coordinate, y is null.

It is upto the user to make sure that the coordinate is increasingly ordered and evenly spaced .

If the values of a function \(f(x)\) and its derivative are know at x=0, and x=1 (normalized coordinate), then the function can be interpolated on the interval [0,1] using a third degree polynomial. This is called cubic interpolation. The formula of this polynomial can be easily derived.

A third degree polynomial and its derivative:

\[ f(x)=ax^3+bx^2+cx+d \]

\[ f(x)=3ax^3+2bx+c \]

The coefficients can be derived from the value and derivatives:

\begin{eqnarray*} a&=&2f(0)-2f(1)+f^\prime (0)+f^\prime(1)\\ b&=&-3f(0)+3f(1)-2f^\prime(0)-f^\prime(0)\\ c&=&f^\prime(0)\\ d&=&f(0)\\ \end{eqnarray*}

the derivatives can be computed as

\begin{eqnarray*} f^\prime(0)&=&\frac{f(1)-f(-1)}{2}\\ f^\prime(1)&=&\frac{f(2)-f(0)}{2}\\ \end{eqnarray*}

so we have the formula

\begin{eqnarray*} a&=&-0.5 f(-1) + 1.5 f(0) - 1.5 f(1) + 0.5 f(2)\\ b&=& f(-1) - 2.5 f(0) + 2 f(1) - 0.5 f(2)\\ c&=&-0.5 f(-1) + 0.5 f(1) \\ d&=& f(0) \\ \end{eqnarray*}

for the boundary pints, replace

\[f^\prime(0)=(f(1)-f(-1))/2\]

by

\[f^\prime(0)=(f(1)-f(0))\]

Otehr type of boundaries are handled in the same way.

see www.paulinternet.nl/?page=bicubicx

sspline_eval()

smat* sspline_eval	(	smat *	coeff,
		smat *	x,
		smat *	xnew
	)

Evluate the cubic spline represented by nx5 matrix coeff, at location array xnew.

sspline()

smat* sspline	(	smat *	x,
		smat *	y,
		smat *	xnew
	)

Do 1D cubic spline all at once by calling SX(spline_prep) and SX(spline_evald)

scwlog10()

void scwlog10

(

smat *

A

)

Do a component wise log10 on each element of A.

scwlog()

void scwlog

(

smat *

A

)

Do a component wise log10 on each element of A.

sfwhm()

float sfwhm

(

smat *

A

)

Calculate number of pixels having values larger than or equal to half of maximum and convert to diameter.

For more accurate calculation of Gaussian beams, use the sqrt(2*log(2))*SX(gauss_width).

sfwhm_gauss()

float sfwhm_gauss

(

smat *

A

)

Use gauss fit to compute fwhm

senc()

smat* senc	(	const smat *	psf,
		const smat *	dvec,
		int	type,
		int	nthread
	)

Compute the enclosed energy or azimuthal average of a.

Parameters

psf	The input array
dvec	The diameter for enclosed energy, or radius for azimuthal average
type	The type. -1: azimuthal average, 0: within a square, 1: within a circle, 2: within a slit
nthread	Number of CPU threads to use.

scellnorm()

float scellnorm

(

const scell *

A

)

compute norm2.

scelldropempty()

void scelldropempty	(	scell **	A0,
		int	dim
	)

drop empty rows or columns. size of *A0 is changed.

scelldot()

float scelldot	(	const scell *	A,
		const scell *	B
	)

Compute the dot produce of two dcell.

scellcwm()

void scellcwm	(	scell *	B,
		const scell *	A
	)

Component wise multiply of two dcell B=A.*B*alpha

scellinvspd()

scell* scellinvspd

(

scell *

A

)

Inplace Invert a SPD matrix. It is treated as a block matrix

scellinv()

scell* scellinv

(

scell *

A

)

Inplace Invert a matrix. It is treated as a block matrix.

scellinvspd_each()

scell* scellinvspd_each

(

scell *

A

)

invert each component of the dcell. Each cell is treated as an individual matrix.

scelldiff()

float scelldiff	(	const scell *	A,
		const scell *	B
	)

compute ||A-B||/||A|| use mean.

scellsub()

scell* scellsub	(	const scell *	in,
		long	sx,
		long	nx,
		long	sy,
		long	ny
	)

Create a new sub cell matrix of nx*ny starting from(sx,sy)

sbspline_prep()

scell* sbspline_prep	(	smat *	x,
		smat *	y,
		smat *	z
	)

2D cubic spline interpolation preparation. x is the x coordinate vector of the 2-d grid. y is the y coordinate vector of the 2-d grid. z is defined on the 2-d grid. It is upto the user to make sure that the coordinate is increasingly ordered and evenly spaced .

The boundaries are handled in the same way is SX(spline). i.e. replace

\[f^\prime(0)=(f(1)-f(-1))/2\]

by

\[f^\prime(0)=(f(1)-f(0))\]

Otehr type of boundaries are handled in the same way.

sbspline_eval()

smat* sbspline_eval	(	scell *	coeff,
		smat *	x,
		smat *	y,
		smat *	xnew,
		smat *	ynew
	)

Evaluate 2D cubic spline at location defined 2-d arrays by xnew, ynew

smaprot()

void smaprot	(	anyarray	A_,
		real	theta
	)

Rotate a or multiple 2d array ccw by angle

swritedata()

void swritedata	(	file_t *	fp,
		const smat *	A,
		long	ncol
	)

Routines for input/output to bin/fits file. Function to write dense matrix data into a file pointer. Generally used by library developer

sreaddata()

smat* sreaddata	(	file_t *	fp,
		header_t *	header
	)

Function to read dense matrix into memory from file pointer. Generally used by library developer.

snew_mmap()

smat* snew_mmap	(	long	nx,
		long	ny,
		const char *	keywords,
		const char *	format,
			...
	)

Create a new SX(mat) matrix object, mmapped from file. The file is truncated if already exists in rw mode.

scellnew_mmap()

scell* scellnew_mmap	(	long	nx,
		long	ny,
		long *	nnx,
		long *	nny,
		const char *	keywords,
		const char *	format,
			...
	)

Create a new dense matrix cell object, mmapped from file. The file is truncated if already exists.

scellnewsame_mmap()

scell* scellnewsame_mmap	(	long	nx,
		long	ny,
		long	mx,
		long	my,
		const char *	keywords,
		const char *	format,
			...
	)

Create a new dense matrix cell object, with identical blocks, mmapped from file. be aware that the data is not 8-byte aligned. The file is truncated if already exists.

sread_mmap()

smat* sread_mmap	(	const char *	format,
			...
	)

Map the file to memory in read only, shared mode.

scellread_mmap()

scell* scellread_mmap	(	const char *	format,
			...
	)

Map the file to memory in read only, shared mode.

sspnew()

ssp* sspnew	(	long	nx,
		long	ny,
		long	nzmax
	)

Create a nx*ny sparse matrix with memory for nmax max elements allocated.

sspref()

ssp* sspref

(

ssp *

A

)

reference a sparse object.

sspdup()

ssp* sspdup

(

const ssp *

A

)

copy a sparse matrix to another.

sspnew2()

ssp* sspnew2

(

const ssp *

A

)

Create a new sparse matrix of the same size as A.

sspsetnzmax()

void sspsetnzmax	(	ssp *	sp,
		long	nzmax
	)

resize a sparse matrix

ssp_cast()

ssp* ssp_cast

(

const cell *

A

)

Cast a pointer to sparse array.

sspfree_do()

void sspfree_do

(

ssp *

sp

)

free a sparse matrix

sspdisp()

void sspdisp

(

const ssp *

sp

)

Display a sparse array

sspcheck()

int sspcheck

(

const ssp *

sp

)

Check a sparse array for wrong orders. Return 1 if is lower triangle, 2 if upper triangle, and 3 if diagonal.

sspscalex()

void sspscalex	(	ssp *	A,
		const smat *	xs
	)

inplace scale sparse matrix elements.

sspscaley()

void sspscaley	(	ssp *	A,
		const smat *	ys
	)

inplace scale sparse matrix elements.

sspcellref()

sspcell* sspcellref

(

const sspcell *

A

)

Reference a spcell array.

sspcell_cast()

sspcell* sspcell_cast

(

const cell *

A

)

cast a cell object to SX(spcell) after checking

sspdiag()

smat* sspdiag

(

const ssp *

A

)

Extract diagonal element of A and return

sspwdot()

float sspwdot	(	const smat *	y,
		const ssp *	A,
		const smat *	x
	)

Multiply two vectors with weighting by sparse matrix. return y'*(A*x)

sspcellwdot()

float sspcellwdot	(	const scell *	y,
		const sspcell *	A,
		const scell *	x
	)

Multiply two cell arrays with weighting by sparse matrix

scellmm2()

scell* scellmm2	(	const scell *	A,
		const scell *	B,
		const char	trans[2]
	)

a different interface for multiplying cells.

ssptrans()

ssp* ssptrans

(

const ssp *

A

)

Transpose a sparse array

sspconj()

void sspconj

(

ssp *

A

)

Take conjugation elementwise.

sspmulsp()

ssp* sspmulsp	(	const ssp *	A,
		const ssp *	B,
		const char	trans[2]
	)

Multiply two sparse arrays and return the result op(A)*op(B)

sspcelltrans()

sspcell* sspcelltrans

(

const sspcell *

spc

)

Transpose a sparse cell

sspcat()

ssp* sspcat	(	const ssp *	A,
		const ssp *	B,
		int	dim
	)

Concatenate two sparse array along dim dimension

sspcell2sp()

ssp* sspcell2sp

(

const sspcell *

A

)

Concatenate a dspcell to sparse array

sspsum()

smat* sspsum	(	const ssp *	A,
		int	dim
	)

Sum elements of sparse array along dimension dim

sspsumabs()

smat* sspsumabs	(	const ssp *	A,
		int	col
	)

Sum abs of elements of sparse array along dimension dim

sspclean()

void sspclean

(

ssp *

A

)

Clean up a sparse array by dropping zeros

sspdropemptycol()

lmat* sspdropemptycol

(

ssp *

A

)

Drop empty columns

sspsort()

void sspsort

(

ssp *

A

)

Make sure the elements are sorted correctly. Does not change the location of data. can be done without harm.

sspcellsort()

void sspcellsort

(

sspcell *

A

)

Make sure the elements are sorted correctly.

sspsym()

void sspsym

(

ssp **

A

)

symmetricize a sparse matrix and drop values below a threshold.

sspcellsym()

void sspcellsym

(

sspcell **

A

)

symmetricize a sparse cell and drop values below a threshold.

sspconvolvop()

ssp* sspconvolvop

(

smat *

A

)

Create a sparse convolution operator C with C(i,j)=A(i-j); A must be very sparse with only a view non-zero value otherwise C will be too full.

sspperm()

ssp* sspperm	(	ssp *	A,
		int	reverse,
		long *	pcol,
		long *	prow
	)

Permute rows and columns of sparse matrix A;

sspinvbdiag()

ssp* sspinvbdiag	(	const ssp *	A,
		long	bs
	)

Invert a SPD sparse matrix that is block diagonal with block sizes of bs.

sspblockextract()

scell* sspblockextract	(	const ssp *	A,
		long	bs
	)

Extrat the diagonal blocks of size bs into cell arrays.

scell2m()

smat* scell2m

(

const_anyarray

A_

)

Convert dense or sparse matrix cell to matrix.

sspwritedata()

void sspwritedata	(	file_t *	fp,
		const ssp *	sp
	)

Function to write sparse matrix data into file pointed using a file pointer. Generally used by library developer. We do not convert data during saving, but rather do the conversion during reading.

sspreaddata()

ssp* sspreaddata	(	file_t *	fp,
		header_t *	header
	)

Function to read sparse matrix data from file pointer into memory. Used by library developer.