ncigt-fil
Functions
Linear System Solvers

standard numerical solutions to linear system problems More...

Functions

double * alloc1double (size_t n1)
 
double ** alloc2double (size_t n1, size_t n2)
 
void cgsolv (COMPLEX *A, int n, COMPLEX *b, COMPLEX *x, int k_max, double thresh)
 
void compute_pinv (double **u, double w[], double **v, short n, short m, double **S_inv)
 
int compute_svd (double **a, int m, int n, double w[], double **v)
 
void eig_jacobi (double **a, double d[], double **v, int n)
 
void eig_sort (double d[], double **v, int n)
 
void free1double (double *p)
 
void free2double (double **p)
 
int lsv (COMPLEX *A, int m, int n, COMPLEX *x0, double *sigma, COMPLEX *v, float thresh)
 
int rank_one_est (COMPLEX *C, float *mask, int Nx, int Ny, double *d, COMPLEX *u, COMPLEX *v, int tid)
 
void svd_backsubstitute (double **u, double w[], double **v, int n, int m, double b[], double x[])
 
void svd_sort (double **u, double *w, double **v, int m, int n)
 

Detailed Description

standard numerical solutions to linear system problems

A linear of system of equations can be ordered as

\[ Ax = b \]

where A is a system matrix of size m-by-n, b is a vector of measured data, and x is the solution vector—which after solution describes a linear combination of values from A to form b.

Function Documentation

double* alloc1double ( size_t  n1)

1D vector memory allocation for SVD routines

double** alloc2double ( size_t  n1,
size_t  n2 
)

2D array allocation, e.g. A[][], for SVD routines.

void cgsolv ( COMPLEX *  A,
int  n,
COMPLEX *  b,
COMPLEX *  x,
int  k_max,
double  thresh 
)

Generic CG solver

Solves the problem Ax = b using the method of conjugate gradients. A must be conjugate symmetric, which implies that x and b are vectors of equal length, n.

Design based on the practical CG algorithm in ``Matrix Computations'' by Golub and Van Loan.

One can specify either a fixed number of iterations, or a residual error stopping criterion. To use the default values, (k_max = 20, thresh = 1e-9), declare '0' as the input values for k_max and thresh.

Parameters
Asystem matrix, must be complex-conjugate symmetric
nlength of data vector
bobjective vector
xsolution
k_maxmaximum number of iterations
threshresidual error threshold
void compute_pinv ( double **  u,
double  w[],
double **  v,
short  n,
short  m,
double **  S_inv 
)

Computes the inverse of the sensitivity matrix, using the output of compute_svd(). It is given by $ v $ times the singular values $ w $ and the traspose of $ u $. This function is simmilar to svd_backsubstitute(), but does not multiply by b. It returns the inverse matrix instead of a solution given by the inverse times b.

Parameters
uleft singular vectors
warray of singular values
vright singular vectors, transposed
nnumber of rows in v
mnumber of rows in u
S_invPseudo-inverse of A
int compute_svd ( double **  a,
int  m,
int  n,
double  w[],
double **  v 
)

Computes the singular value decomposition (SVD) of a matrix, $ A = U \Sigma V^H $, where $ U $ and $ V $ are matrices with orthogonal columns, and $ \Sigma $ is a diagonal matrix holding the singular values, $ \sigma_i $, associated with the singular vectors.

The most common use for an SVD is to compute a psuedo-inverse of the matrix $ A $, in order to solve the linear system $ Ax = b $. Once the SVD has been computed, the solution can be formed via

\[ x = V \Sigma^{-1} U^T b \]

This can be perfomed computationally using the function svd_backsubstitute().

In the case of small singular values, e.g. $ 0 < \sigma_i < 1 $, computing the inverse of $ \Sigma $ will induce excessively large contributions from the associated singular vectors. Thus, it is common to specify a threshold—discarding singular values and associated vectors below the threshold to limit such effects. To peform this in C programs, one needs to set all singular values below the threshold to zero.

Note, this function accepts only real-valued input. To find the SVD of a complex-valued matrix, one can split the real and imaginary components as

\[ A x = (\mathrm{re}\{A\} + j \mathrm{im}\{A\}) x = b \]

\[ \left[ \begin{array}{cc} \mathrm{re}\{A\} & -\mathrm{im}\{A\} \\ \mathrm{im}\{A\} & \mathrm{re}\{A\} \end{array} \right] \left[ \begin{array}{c} \mathrm{re}\{x\} \\ \mathrm{im}\{x\} \end{array} \right] = \left[ \begin{array}{c} \mathrm{re}\{b\} \\ \mathrm{im}\{b\} \end{array} \right] \]

Usage:

  double **A;        
  double **U;        
  double *W;         
  double **V;        
  double **A_inv;
  A = alloc2double(m,n);
  V = alloc2double(n,n);
  W = alloc1double(n);
  A_inv = alloc2double(n,m);
  compute_svd( A, m, n, W, V );
  U = A;
  compute_pinv( U, W, V, n, m, A_inv);
  free2double(A);
  free2double(V);
  free1double(W);
  free2double(A_inv);

Credits: from CWP/SU. Similar to code in Netlib's EISPACK and CLAPACK. See also discussions in NR in C

Parameters
ainput matrix to decompose, will be overwritten by the left signular vectors
mnumber of rows in a
nnumber of columns in a
wan array of output singular values
va matrix output of the right signular vectors, transposed
void eig_jacobi ( double **  a,
double  d[],
double **  v,
int  n 
)

Find eigenvalues and corresponding eigenvectors via the Jacobi algorithm for symmetric matrices

Credits: from CWP/SU.

Parameters
aa real-valued symmetric input matrix to decompose
dan output array of eigen values
van output array of eigen vectors
nthe size (rows,cols) of the input matrix
void eig_sort ( double  d[],
double **  v,
int  n 
)

sort eigenvalues and corresponding eigenvectors in descending order

Credits: derived from CWP/SU.

Parameters
darray of eigen-values
varray of eigen-vectors
nsize of eigen system
void free1double ( double *  p)

free memory allocated by alloc1double().

void free2double ( double **  p)

free memory allocated by alloc2double().

int lsv ( COMPLEX *  A,
int  m,
int  n,
COMPLEX *  x0,
double *  sigma,
COMPLEX *  v,
float  thresh 
)

Estimate of the dominant Left Singular Vector

The function employs an iterative algorithm that solves

\[ \max_x \| A^H x \| \]

to estimate the left singular vector associated with the largest singular value. ( from this vector, the singular value and associated right singular vector can also be determined. )

If A is tall-and-thin, i.e. the number of rows, m, is much larger than the number of columns, n, then this algorithm is significantly faster than computing the SVD directly.

The function is most useful when only a rank-one estimate of the matrix is desired.

Parameters
Athe input matrix.
mnumber of rows in input matrix.
nnumber of columns in input matrix. Note, it is assumed that m>>n
x0the output vector. size: m x 1 (needs to be pre-allocated)
thresh[opt] convergent threshold (default=1e-4)
int rank_one_est ( COMPLEX *  C,
float *  mask,
int  Nx,
int  Ny,
double *  d,
COMPLEX *  u,
COMPLEX *  v,
int  tid 
)

uses Lanczos algorithm to determine singular vectors associated with the largest signular value.

Parameters
Cinput data
Nxnumber of rows in input data
Nynumber of columns in input data
dthe largest signular value. (Set to NULL if undesired)
uthe right singular vector associated with the largest singular value. (Set to NULL if undesired)
vthe left signular vector associated with the largest singular value. (Set to NULL if undesired)
tidthread id for the Lanczos algorithm code. Set to zero is not using threads.
void svd_backsubstitute ( double **  u,
double  w[],
double **  v,
int  n,
int  m,
double  b[],
double  x[] 
)

Perform back substitution of the singular value decomposittion, operation on the objective, $ b $, to compute the solution, $ x $, of the linear system of equations.

Back-substition can be performed much more quickly that explicitly computing the pseudo-inverse. However, if the inverse will be used to operate on multiple data sets (as in UNFOLD-SENSE), it is preferable to compute the inverse once and then apply it to all time frames in the series.

Credits: from CWP/SU. Similar to code in Netlib's EISPACK and CLAPACK. See also discussions in Numerical Recipes in C.

Parameters
uLeft singular vectors
wArray of singular values
vright signular vectors
nnumber of rows in u
mnumber of columns in v
bobjective
xoutput
void svd_sort ( double **  u,
double *  w,
double **  v,
int  m,
int  n 
)

Sort the singular values and corresponding singular vectors in descending order

Assumes the input of the singular value decomposition of a matrix a[][] of the form: a[][] = u[][] w[] vt[][]

Credits: Nils Maercklin, GeoForschungsZentrum (GFZ) Potsdam, Germany, 2001.

Parameters
uleft singular vectors
warray of singular values
vright singular vectors, transposed