#include "xcomplex.h"
#include "arrayalloc.h"

Include dependency graph for matrix.h:

This graph shows which files directly or indirectly include this file:

Macros
#define	Multiply_full(MULTIPLY, MULTIADD, xx, yy)

#define	Multiply_square(MULTIPLY, MULTIADD)

#define	Multiply_xx(MULTIPLY, MULTIADD, xx1, xx2, zzt)

Functions
void	BalanceSim (int n, double **A)

int	Choleski (int N, double L, double A)

double **	Copy (int M, int N, double Z, double A, MList *List=0)

double **	Copy (int M, int N, double *A, MList List=0)

double **	Copy (int N, double Z, double A, MList *List=0)

double **	Copy (int N, double *A, MList List=0)

cdouble **	Copy (int M, int N, cdouble Z, cdouble A, MList *List=0)

cdouble **	Copy (int M, int N, cdouble *A, MList List=0)

cdouble **	Copy (int N, cdouble Z, cdouble A, MList *List=0)

cdouble **	Copy (int N, cdouble *A, MList List=0)

double *	Copy (int N, double z, double a, MList *List=0)

double *	Copy (int N, double a, MList List=0)

cdouble *	Copy (int N, cdouble z, cdouble a, MList *List=0)

cdouble *	Copy (int N, cdouble a, MList List=0)

double	det (int N, double **A, int mode=0)

cdouble	det (int N, cdouble **A, int mode=0)

int	EigPower (int N, double &l, double x, double *A, double ep=1e-6, int maxiter=50)

int	EigPowerA (int N, double &l, double x, double *A, double ep=1e-6, int maxiter=50)

int	EigPowerI (int N, double &l, double x, double *A, double ep=1e-6, int maxiter=50)

int	EigPowerS (int N, double &l, double x, double *A, double ep=1e-6, int maxiter=50)

int	EigPowerWielandt (int N, double &m, double u, double l, double v, double **A, double ep=1e-06, int maxiter=50)

int	EigenValues (int N, double *A, cdouble ev)

int	EigSym (int N, double *A, double d, double **Q)

int	GEB (int N, double x, double A, double b)

int	GESCP (int N, double x, double A, double b)

void	GExL (int N, double x, double L, double a)

void	GExLAdd (int N, double x, double L, double a)

void	GExL1 (int N, double x, double *L, double a)

void	GExL1Add (int N, double x, double *L, double a)

template<class T , class Ta >
int	GECP (int N, T x, Ta A, T b, Ta *logdet=0)

double	GILT (int N, double **A)

double	GIUT (int N, double **A)

double	GICP (int N, double X, double A)

double	GICP (int N, double **A)

cdouble	GICP (int N, cdouble **A)

double	GISCP (int N, double X, double A)

double	GISCP (int N, double **A)

int	GSI (int N, double x0, double A, double b, double ep=1e-4, int maxiter=50)

void	Hessenb (int n, double **A)

void	HouseHolder (int N, double **A)

void	HouseHolder (int N, double T, double Q, double **A)

void	HouseHolder (int n, double *A, double d, double *sd)

double	Inner (int N, double x, double y)

cdouble	Inner (int N, double x, cdouble y)

cdouble	Inner (int N, cdouble x, cdouble y)

cdouble	Inner (int N, cfloat x, cdouble y)

cfloat	Inner (int N, cfloat x, cfloat y)

double	Inner (int M, int N, double X, double Y)

template<class Tx , class Tw >
cdouble	Inner (int N, Tx x, Tw w, cdouble *y)

template<class Tx , class Tw >
cdouble	Inner (int N, Tx x, Tw w, double *y)

int	JI (int N, double x, double A, double b, double ep=1e-4, int maxiter=50)

int	LDL (int N, double *L, double d, double **A)

void	LQ_GS (int M, int N, double A, double L, double **Q)

template<class T >
T	LSLinear (int M, int N, T x, T A, T y)

void	LSLinear2 (int M, int N, double x, double A, double y)

int	LU_Direct (int mode, int N, double diag, double *a)

int	LU (int N, double l, double u, double **a)

void	LU_DiagL (int N, double *l, double diagl, double u, double a)

int	LU_PD (int N, double **b)

int	LU_PD (int N, double *b, double c)

double	LUCP (double *a, int N, int ind, int &parity, int mode=1)

void	LU (int N, double x, double A, double y, int *ind=0)

int	maxind (double *data, int from, int to)

template<class Tz , class Tx , class Ty , class Ta >
Tz **	MultiAdd (int M, int N, Tz Z, Tx X, Ty *Y, Ta a, MList List=0)

template<class Tz , class Tx , class Ty , class Ta >
Tz *	MultiAdd (int N, Tz z, Tx x, Ty y, Ta a, MList List=0)

template<class Tx , class Ty , class Ta >
Tx *	MultiAdd (int N, Tx x, Ty y, Ta a, MList *List=0)

template<class Tz , class Tx , class Ta >
Tz **	Multiply (int M, int N, Tz Z, Tx X, Ta a, MList *List=0)

template<class Tz , class Tx , class Ta >
Tz *	Multiply (int N, Tz z, Tx x, Ta a, MList *List=0)

template<class Tx , class Ta >
Tx *	Multiply (int N, Tx x, Ta a, MList List=0)

	Multiply_full (MultiplyXY, MultiAddXY, x[n][m], y[m][k]) Multiply_full(MultiplyXYc

*y[m][k]	Multiply_full (MultiplyXYt, MultiAddXYt, x[n][m], y[k][m]) Multiply_full(MultiplyXYh

	Multiply_square (MultiplyXY, MultiAddXY) Multiply_xx(MultiplyXtX

zz	Multiply_xx (MultiplyXXt, MultiAddXXt, x[m][n], x[k][n], zz) Multiply_xx(MultiplyXhX

zz *zz	Multiply_xx (MultiplyXXh, MultiAddXXh, x[m][n],x[k][n],zz) Multiply_xx(MultiplyXcXt

zz zz zz double *	MultiplyXy (int M, int N, double z, double x, double y, MList *List=0)

double *	MultiplyXy (int M, int N, double *x, double y, MList *List=0)

cdouble *	MultiplyXy (int M, int N, cdouble z, cdouble x, cdouble y, MList *List=0)

cdouble *	MultiplyXy (int M, int N, cdouble *x, cdouble y, MList *List=0)

cdouble *	MultiplyXy (int M, int N, cdouble z, double x, cdouble y, MList *List=0)

cdouble *	MultiplyXy (int M, int N, double *x, cdouble y, MList *List=0)

double *	MultiplyxY (int M, int N, double zt, double xt, double *y, MList List=0)

double *	MultiplyxY (int M, int N, double xt, double y, MList List=0)

cdouble *	MultiplyxY (int M, int N, cdouble zt, cdouble xt, cdouble *y, MList List=0)

cdouble *	MultiplyxY (int M, int N, cdouble xt, cdouble y, MList List=0)

double *	MultiplyXty (int M, int N, double z, double x, double y, MList *List=0)

double *	MultiplyXty (int M, int N, double *x, double y, MList *List=0)

cdouble *	MultiplyXty (int M, int N, cdouble z, cdouble x, cdouble y, MList *List=0)

cdouble *	MultiplyXty (int M, int N, cdouble *x, cdouble y, MList *List=0)

cdouble *	MultiplyXhy (int M, int N, cdouble z, cdouble x, cdouble y, MList *List=0)

cdouble *	MultiplyXhy (int M, int N, cdouble *x, cdouble y, MList *List=0)

cdouble *	MultiplyXcy (int M, int N, cdouble z, cdouble x, cfloat y, MList *List=0)

cdouble *	MultiplyXcy (int M, int N, cdouble *x, cfloat y, MList *List=0)

cdouble *	MultiplyXcy (int M, int N, cdouble z, cdouble x, cdouble y, MList *List=0)

cdouble *	MultiplyXcy (int M, int N, cdouble *x, cdouble y, MList *List=0)

double *	MultiplyxYt (int M, int N, double zt, double xt, double *y, MList MList=0)

double *	MultiplyxYt (int M, int N, double xt, double y, MList MList=0)

double	Norm1 (int N, double **A)

void	QR_GS (int M, int N, double A, double Q, double **R)

void	QR_householder (int M, int N, double A, double Q, double **R)

int	QL (int N, double d, double sd, double **z)

int	QR (int N, double *A, cdouble ev)

void	QU (int N, double Q, double A)

double **	Real (int M, int N, double z, cdouble x, MList *List)

double **	Real (int M, int N, cdouble *x, MList List)

int	Roots (int N, double p, double rr, double *ri)

int	SorI (int N, double x0, double a, double b, double w=1, double ep=1e-4, int maxiter=50)

void	SetSubMatrix (double z, double x, int Y1, int Y, int X1, int X)

void	SetSubMatrix (cdouble z, cdouble x, int Y1, int Y, int X1, int X)

cdouble **	SubMatrix (cdouble z, cdouble x, int Y1, int Y, int X1, int X, MList *List=0)

cdouble **	SubMatrix (cdouble *x, int Y1, int Y, int X1, int X, MList List=0)

cdouble *	SubVector (cdouble z, cdouble x, int X1, int X, MList *List=0)

cdouble *	SubVector (cdouble x, int X1, int X, MList List=0)

void	transpose (int N, double **a)

void	transpose (int N, cdouble **a)

double **	transpose (int N, int M, double ta, double a, MList *List=0)

double **	transpose (int N, int M, double *a, MList List=0)

double **	Unitary (int N, double *P, double x, double y, MList List=0)

double **	Unitary (int N, double x, double y, MList *List=0)

cdouble **	Unitary (int N, cdouble *P, cdouble x, cdouble y, MList List=0)

cdouble **	Unitary (int N, cdouble x, cdouble y, MList *List=0)

Variables
	MultiAddXYc

	x [n][m]

*y[m][k]	MultiAddXYh

	MultiAddXtX

zz	MultiAddXhX

zz *zz	MultiAddXcXt

Detailed Description

matrix operations.

Matrices are accessed by double pointers as MATRIX[Y][X], where Y is the row index.

Macro Definition Documentation

#define Multiply_full	(	MULTIPLY,
		MULTIADD,
		xx,
		yy
	)

Value:

template<class Tx, class Ty>Tx** MULTIPLY(int N, int M, int K, Tx** z, Tx** x, Ty** y, MList* List=0){ \
                if (!z) {Allocate2(Tx, N, K, z); if (List) List->Add(z, 2);} \
    for (int n=0; n<N; n++) for (int k=0; k<K; k++){ \
      Tx zz=0; for (int m=0; m<M; m++) zz+=xx*yy; z[n][k]=zz;} return z;} \
  template<class Tx, class Ty>Tx** MULTIPLY(int N, int M, int K, Tx** x, Ty** y, MList* List=0){ \
    Tx** Allocate2(Tx, N, K, z); if (List) List->Add(z, 2); \
    for (int n=0; n<N; n++) for (int k=0; k<K; k++){ \
      Tx zz=0; for (int m=0; m<M; m++) zz+=xx*yy; z[n][k]=zz;} return z;} \
  template<class Tx, class Ty>void MULTIADD(int N, int M, int K, Tx** z, Tx** x, Ty** y){ \
    for (int n=0; n<N; n++) for (int k=0; k<K; k++){ \
      Tx zz=0; for (int m=0; m<M; m++) zz+=xx*yy; z[n][k]+=zz;}}

#define Multiply_square	(	MULTIPLY,
		MULTIADD
	)

Value:

template<class T>T** MULTIPLY(int N, T** z, T** x, T** y, MList* List=0){ \
    if (!z){Allocate2(T, N, N, z); if (List) List->Add(z, 2);} \
    if (z!=x) MULTIPLY(N, N, N, z, x, y); else{ \
      T* tmp=new T[N]; int sizeN=sizeof(T)*N; for (int i=0; i<N; i++){ \
        MULTIPLY(1, N, N, &tmp, &x[i], y); memcpy(z[i], tmp, sizeN);} delete[] tmp;} \
    return z;} \
  template<class T>T** MULTIPLY(int N, T** x, T** y, MList* List=0){return MULTIPLY(N, N, N, x, y, List);}\
  template<class T>void MULTIADD(int N, T** z, T** x, T** y){MULTIADD(N, N, N, z, x, y);}

#define Multiply_xx	(	MULTIPLY,
		MULTIADD,
		xx1,
		xx2,
		zzt
	)

Value:

template<class T>T** MULTIPLY(int M, int N, T** z, T** x, MList* List=0){ \
    if (!z){Allocate2(T, M, M, z); if (List) List->Add(z, 2);} \
    for (int m=0; m<M; m++) for (int k=0; k<=m; k++){ \
      T zz=0; for (int n=0; n<N; n++) zz+=xx1*xx2; z[m][k]=zz; if (m!=k) z[k][m]=zzt;} \
    return z;} \
  template<class T>T** MULTIPLY(int M, int N, T ** x, MList* List=0){ \
    T** Allocate2(T, M, M, z); if (List) List->Add(z, 2); \
    for (int m=0; m<M; m++) for (int k=0; k<=m; k++){ \
      T zz=0; for (int n=0; n<N; n++) zz+=xx1*xx2; z[m][k]=zz; if (m!=k) z[k][m]=zzt;} \
    return z;} \
  template<class T>void MULTIADD(int M, int N, T** z, T** x){ \
    for (int m=0; m<M; m++) for (int k=0; k<=m; k++){ \
      T zz=0; for (int n=0; n<N; n++) zz+=xx1*xx2; z[m][k]+=zz; if (m!=k) z[k][m]+=zzt;}}

Function Documentation

void BalanceSim	(	int	n,
		double **	A
	)

function BalanceSim: applies a similarity transformation to matrix a so that a is "balanced". This is used by various eigenvalue evaluation routines.

In: matrix A[n][n] Out: balanced matrix a

No return value.

int Choleski	(	int	N,
		double **	L,
		double **	A
	)

function Choleski: Choleski factorization A=LL', where L is lower triangular. The symmetric matrix A[N][N] is positive definite iff A can be factored as LL', where L is lower triangular with nonzero diagonl entries.

In: matrix A[N][N] Out: mstrix L[N][N].

Returns 0 if successful. On return content of matrix a is not changed.

double** Copy	(	int	M,
		int	N,
		double **	Z,
		double **	A,
		MList *	List
	)

function Copy: duplicate the matrix A as matrix Z.

In: matrix A[M][N] Out: matrix Z[M][N]

Returns pointer to Z. Z is created anew if Z=0 is supplied on start.

double* Copy	(	int	N,
		double *	z,
		double *	a,
		MList *	List
	)

function Copy: duplicating vector a as vector z

In: vector a[N] Out: vector z[N]

Returns pointer to z. z is created anew is z=0 is specified on start.

double det	(	int	N,
		double **	A,
		int	mode
	)

function det: computes determinant by Gaussian elimination method with column pivoting

In: matrix A[N][N]

Returns det(A). On return content of matrix A is unchanged if mode=0.

int EigenValues	(	int	N,
		double **	A,
		cdouble *	ev
	)

function EigenValues: solves for eigenvalues of general system

In: matrix A[N][N] Out: eigenvalues ev[N]

Returns 0 if successful. Content of matrix A is destroyed on return.

int EigPower	(	int	N,
		double &	l,
		double *	x,
		double **	A,
		double	ep,
		int	maxiter
	)

function EigPower: power method for solving dominant eigenvalue and eigenvector

In: matrix A[N][N], initial arbitrary vector x[N]. Out: eigenvalue l, eigenvector x[N].

Returns 0 is successful. Content of matrix A is unchangd on return. Initial x[N] must not be zero.

int EigPowerA	(	int	N,
		double &	l,
		double *	x,
		double **	A,
		double	ep,
		int	maxiter
	)

function EigPowerA: EigPower with Aitken acceleration

In: matrix A[N][N], initial arbitrary vector x[N]. Out: eigenvalue l, eigenvector x[N].

Returns 0 is successful. Content of matrix A is unchangd on return. Initial x[N] must not be zero.

int EigPowerI	(	int	N,
		double &	l,
		double *	x,
		double **	A,
		double	ep,
		int	maxiter
	)

function EigPowerI: Inverse power method for solving the eigenvalue given an approximate non-zero eigenvector.

In: matrix A[N][N], approximate eigenvector x[N]. Out: eigenvalue l, eigenvector x[N].

Returns 0 is successful. Content of matrix A is unchangd on return. Initial x[N] must not be zero.

int EigPowerS	(	int	N,
		double &	l,
		double *	x,
		double **	A,
		double	ep,
		int	maxiter
	)

function EigPowerS: symmetric power method for solving the dominant eigenvalue with its eigenvector

In: matrix A[N][N], initial arbitrary vector x[N]. Out: eigenvalue l, eigenvector x[N].

Returns 0 is successful. Content of matrix A is unchangd on return. Initial x[N] must not be zero.

int EigPowerWielandt	(	int	N,
		double &	m,
		double *	u,
		double	l,
		double *	v,
		double **	A,
		double	ep,
		int	maxiter
	)

function EigPowerWielandt: Wielandt's deflation algorithm for solving a second dominant eigenvalue and eigenvector (m,u) given the dominant eigenvalue and eigenvector (l,v).

In: matrix A[N][N], first eigenvalue l with eigenvector v[N] Out: second eigenvalue m with eigenvector u

Returns 0 if successful. Content of matrix A is unchangd on return. Initial u[N] must not be zero.

int EigSym	(	int	N,
		double **	A,
		double *	d,
		double **	Q
	)

function EigSym: Solves real symmetric eigensystem A

In: matrix A[N][N] Out: eigenvalues d[N], transform matrix Q[N][N], so that diag(d)=Q'AQ, A=Q diag(d) Q', AQ=Q diag(d)

Returns 0 if successful. Content of matrix A is unchanged on return.

int GEB	(	int	N,
		double *	x,
		double **	A,
		double *	b
	)

function GEB: Gaussian elimination with backward substitution for solving linear system Ax=b.

In: coefficient matrix A[N][N], vector b[N] Out: vector x[N]

Returns 0 if successful. Contents of matrix A and vector b are destroyed on return.

int GESCP	(	int	N,
		double *	x,
		double **	A,
		double *	b
	)

function GESCP: Gaussian elimination with scaled column pivoting for solving linear system Ax=b

In: matrix A[N][N], vector b[N] Out: vector x[N]

Returns 0 is successful. Contents of matrix A and vector b are destroyed on return.

void GExL	(	int	N,
		double *	x,
		double **	L,
		double *	a
	)

function GExL: solves linear system xL=a, L being lower-triangular. This is used in LU factorization for solving linear systems.

In: lower-triangular matrix L[N][N], vector a[N] Out: vector x[N]

No return value. Contents of matrix L and vector a are unchanged at return.

void GExL1	(	int	N,
		double *	x,
		double **	L,
		double	a
	)

function GExL1: solves linear system xL=(0, 0, ..., 0, a)', L being lower-triangular.

In: lower-triangular matrix L[N][N], a Out: vector x[N]

No return value. Contents of matrix L and vector a are unchanged at return.

void GExL1Add	(	int	N,
		double *	x,
		double **	L,
		double	a
	)

function GExL1Add: solves linear system *L=(0, 0, ..., 0, a)', L being lower-triangular, and add the solution * to x[].

In: lower-triangular matrix L[N][N], vector a Out: updated vector x[N]

No return value. Contents of matrix L and vector a are unchanged at return.

void GExLAdd	(	int	N,
		double *	x,
		double **	L,
		double *	a
	)

function GExLAdd: solves linear system *L=a, L being lower-triangular, and add the solution * to x[].

In: lower-triangular matrix L[N][N], vector a[N] Out: updated vector x[N]

No return value. Contents of matrix L and vector a are unchanged at return.

double GICP	(	int	N,
		double **	A
	)

function GICP: matrix inverse using Gaussian elimination with column pivoting: inv(A)->A.

In: matrix A[N][N] Out: matrix A[N][N]

Returns the determinant of the inverse matrix, 0 on failure.

double GILT	(	int	N,
		double **	A
	)

function GILT: inv(lower trangular of A)->lower trangular of A

In: matrix A[N][N] Out: matrix A[N][N]

Returns the determinant of the lower trangular of A

double GISCP	(	int	N,
		double **	X,
		double **	A
	)

function GISCP: wrapper function that does not overwrite input matrix A: inv(A)->X.

In: matrix A[N][N] Out: matrix X[N][N]

Returns the determinant of the inverse matrix, 0 on failure.

double GISCP	(	int	N,
		double **	A
	)

function GISCP: matrix inverse using Gaussian elimination w. scaled column pivoting: inv(A)->A.

In: matrix A[N][N] Out: matrix A[N][N]

Returns the determinant of the inverse matrix, 0 on failure.

double GIUT	(	int	N,
		double **	A
	)

function GIUT: inv(upper trangular of A)->upper trangular of A

In: matrix A[N][N] Out: matrix A[N][N]

Returns the determinant of the upper trangular of A

int GSI	(	int	N,
		double *	x0,
		double **	A,
		double *	b,
		double	ep,
		int	maxiter
	)

function GSI: Gaussian-Seidel iterative algorithm for solving linear system Ax=b. Breaks down if any Aii=0, like the Jocobi method JI(...).

Gaussian-Seidel iteration is x(k)=(D-L)^(-1)(Ux(k-1)+b), where D is diagonal, L is lower triangular, U is upper triangular and A=L+D+U.

In: matrix A[N][N], vector b[N], initial vector x0[N] Out: vector x0[N]

Returns 0 is successful. Contents of matrix A and vector b remain unchanged on return.

void Hessenb	(	int	N,
		double **	A
	)

function Hessenb: reducing a square matrix A to upper Hessenberg form

In: matrix A[N][N] Out: matrix A[N][N], in upper Hessenberg form

No return value.

void HouseHolder	(	int	N,
		double **	A
	)

function HouseHolder: house holder method converting a symmetric matrix into a tridiagonal symmetric matrix, or a non-symmetric matrix into an upper-Hessenberg matrix, using similarity transformation.

In: matrix A[N][N] Out: matrix A[N][N] after transformation

No return value.

void HouseHolder	(	int	N,
		double **	T,
		double **	Q,
		double **	A
	)

function HouseHolder: house holder transformation T=Q'AQ or A=QTQ', where T is tridiagonal and Q is unitary i.e. QQ'=I.

In: matrix A[N][N] Out: matrix tridiagonal matrix T[N][N] and unitary matrix Q[N][N]

No return value. Identical A and T allowed. Content of matrix A is unchanged if A!=T.

void HouseHolder	(	int	N,
		double **	A,
		double *	d,
		double *	sd
	)

function HouseHolder: nr version of householder method for transforming symmetric matrix A to QTQ', where T is tridiagonal and Q is orthonormal.

In: matrix A[N][N] Out: A[N][N]: now containing Q d[N]: containing diagonal elements of T sd[N]: containing subdiagonal elements of T as sd[1:N-1].

No return value.

double Inner	(	int	N,
		double *	x,
		double *	y
	)

function Inner: inner product z=y'x

In: vectors x[N], y[N]

Returns inner product of x and y.

double Inner	(	int	M,
		int	N,
		double **	X,
		double **	Y
	)

function Inner: inner product z=tr(Y'X)

In: matrices X[M][N], Y[M][N]

Returns inner product of X and Y.

int JI	(	int	N,
		double *	x0,
		double **	A,
		double *	b,
		double	ep,
		int	maxiter
	)

function JI: Jacobi interative algorithm for solving linear system Ax=b Breaks down if A[i][i]=0 for any i. Reorder A so that this does not happen.

Jacobi iteration is x(k)=D^(-1)((L+U)x(k-1)+b), D is diagonal, L is lower triangular, U is upper triangular and A=L+D+U.

In: matrix A[N][N], vector b[N], initial vector x0[N] Out: vector x0[N]

Returns 0 if successful. Contents of matrix A and vector b are unchanged on return.

int LDL	(	int	N,
		double **	L,
		double *	d,
		double **	A
	)

function LDL: LDL' decomposition A=LDL', where L is lower triangular and D is diagonal identical l and a allowed.

The symmetric matrix A is positive definite iff A can be factorized as LDL', where L is lower triangular with ones on its diagonal and D is diagonal with positive diagonal entries.

If a symmetric matrix A can be reduced by Gaussian elimination without row interchanges, then it can be factored into LDL', where L is lower triangular with ones on its diagonal and D is diagonal with non-zero diagonal entries.

In: matrix A[N][N] Out: lower triangular matrix L[N][N], vector d[N] containing diagonal elements of D

Returns 0 if successful. Content of matrix A is unchanged on return.

void LQ_GS	(	int	M,
		int	N,
		double **	A,
		double **	L,
		double **	Q
	)

function LQ_GS: LQ decomposition using Gram-Schmidt method

In: matrix A[M][N], M<=N Out: matrices L[M][M], Q[M][N]

No return value.

void LSLinear2	(	int	M,
		int	N,
		double *	x,
		double **	A,
		double *	y
	)

function LSLinear2: 2-dtage LS solution of A[M][N]x[N][1]=y[M][1], M>=N. Use of this function requires the submatrix A[N][N] be invertible.

In: matrix A[M][N], vector y[M], M>=N. Out: vector x[N].

No return value. Contents of matrix A and vector y are unchanged on return.

int LU	(	int	N,
		double **	L,
		double **	U,
		double **	A
	)

function LU: LU decomposition A=LU, where L is lower triangular with diagonal entries 1 and U is upper triangular.

LU is possible if A can be reduced by Gaussian elimination without row interchanges.

In: matrix A[N][N] Out: matrices L[N][N] and U[N][N], subject to input values of L and U: if L euqals NULL, L is not returned if U equals NULL or A, U is returned in A, s.t. A is modified if L equals A, L is returned in A, s.t. A is modified if L equals U, L and U are returned in the same matrix when L and U are returned in the same matrix, diagonal of L (all 1) is not returned

Returns 0 if successful.

void LU	(	int	N,
		double *	x,
		double **	A,
		double *	y,
		int *	ind
	)

function LU: Solving linear system Ax=y by LU factorization

In: matrix A[N][N], vector y[N] Out: x[N]

No return value. On return A contains its LU factorization (with pivoting, diag mode 1), y remains unchanged.

int LU_Direct	(	int	mode,
		int	N,
		double *	diag,
		double **	A
	)

function LU_Direct: LU factorization A=LU.

In: matrix A[N][N], vector diag[N] specifying main diagonal of L or U, according to mode (0=LDiag, 1=UDiag). Out: matrix A[N][N] now containing L and U.

Returns 0 if successful.

int LU_PD	(	int	N,
		double **	b
	)

function LU_PD: LU factorization for pentadiagonal A=LU

In: pentadiagonal matrix A[N][N] stored in a compact format, i.e. A[i][j]->b[i-j, j] the main diagonal is b[0][0]~b[0][N-1] the 1st upper subdiagonal is b[-1][1]~b[-1][N-1] the 2nd upper subdiagonal is b[-2][2]~b[-2][N-1] the 1st lower subdiagonal is b[1][0]~b[1][N-2] the 2nd lower subdiagonal is b[2][0]~b[2][N-3]

Out: L[N][N] and U[N][N], main diagonal of L being all 1 (probably), stored in a compact format in b[-2:2][N].

Returns 0 if successful.

int LU_PD	(	int	N,
		double **	b,
		double *	c
	)

function LU_PD: solve pentadiagonal system Ax=c

In: pentadiagonal matrix A[N][N] stored in a compact format in b[-2:2][N], vector c[N] Out: vector c now containing x.

Returns 0 if successful. On return b is in the LU form.

double LUCP	(	double **	A,
		int	N,
		int *	ind,
		int &	parity,
		int	mode
	)

function LUCP: LU decomposition A=LU with column pivoting

In: matrix A[N][N] Out: matrix A[N][N] now holding L and U by L_U[i][j]=A[ind[i]][j], where L_U hosts L and U according to mode: mode=0: L diag=abs(U diag), U diag as return mode=1: L diag=1, U diag as return mode=2: U diag=1, L diag as return

Returns the determinant of A.

int maxind	(	double *	data,
		int	from,
		int	to
	)

function maxind: returns the index of the maximal value of data[from:(to-1)].

In: vector data containing at least $to entries. Out: the index to the maximal entry of data[from:(to-1)]

Returns the index to the maximal value.

int QL	(	int	N,
		double *	d,
		double *	sd,
		double **	z
	)

function QL: QL method for solving tridiagonal symmetric matrix eigenvalue problem.

In: A[N][N]: tridiagonal symmetric matrix stored in d[N] and sd[] arranged so that d[0:n-1] contains the diagonal elements of A, sd[0]=0, sd[1:n-1] contains the subdiagonal elements of A. z[N][N]: pre-transform matrix z[N][N] compatible with HouseHolder() routine. Out: d[N]: the eigenvalues of A z[N][N] the eigenvectors of A.

Returns 0 if successful. sd[] should have storage for at least N+1 entries.

int QR	(	int	N,
		double **	A,
		cdouble *	ev
	)

function QR: nr version of QR method for solving upper Hessenberg system A. This is compatible with Hessenb method.

In: matrix A[N][N] Out: vector ev[N] of eigenvalues

Returns 0 on success. Content of matrix A is destroyed on return.

void QR_GS	(	int	M,
		int	N,
		double **	A,
		double **	Q,
		double **	R
	)

function QR_GS: QR decomposition A=QR using Gram-Schmidt method

In: matrix A[M][N], M>=N Out: Q[M][N], R[N][N]

No return value.

void QR_householder	(	int	M,
		int	N,
		double **	A,
		double **	Q,
		double **	R
	)

function QR_householder: QR decomposition using householder transform

In: A[M][N], M>=N Out: Q[M][M], R[M][N]

No return value.

void QU	(	int	N,
		double **	Q,
		double **	A
	)

function QU: Unitary decomposition A=QU, where Q is unitary and U is upper triangular

In: matrix A[N][N] Out: matrices Q[N][N], A[n][n] now containing U

No return value.

double** Real	(	int	M,
		int	N,
		double **	z,
		cdouble **	x,
		MList *	List
	)

function Real: extracts the real part of matrix X

In: matrix x[M][N]; Out: matrix z[M][N]

Returns pointer to z. z is created anew if z=0 is specified on start.

void SetSubMatrix	(	double **	z,
		double **	x,
		int	Y1,
		int	Y,
		int	X1,
		int	X
	)

function SetSubMatrix: copy matrix x[Y][X] into matrix z at (Y1, X1).

In: matrix x[Y][X], matrix z with dimensions no less than [Y+Y1][X+X1] Out: matrix z, updated.

No return value.

int SorI	(	int	N,
		double *	x0,
		double **	a,
		double *	b,
		double	w,
		double	ep,
		int	maxiter
	)

function SorI: Sor iteration algorithm for solving linear system Ax=b.

Sor method is an extension of the Gaussian-Siedel method, with the latter equivalent to the former with w set to 1. The Sor iteration is given by x(k)=(D-wL)^(-1)(((1-w)D+wU)x(k-1)+wb), where 0<w<2, D is diagonal, L is lower triangular, U is upper triangular and A=L+D+U. Sor method converges if A is positive definite.

In: matrix A[N][N], vector b[N], initial vector x0[N] Out: vector x0[N]

Returns 0 if successful. Contents of matrix A and vector b are unchanged on return.

cdouble** SubMatrix	(	cdouble **	z,
		cdouble **	x,
		int	Y1,
		int	Y,
		int	X1,
		int	X,
		MList *	List
	)

function SubMatrix: extract a submatrix of x at (Y1, X1) to z[Y][X].

In: matrix x of dimensions no less than [Y+Y1][X+X1] Out: matrix z[Y][X].

Returns pointer to z. z is created anew if z=0 is specifid on start.

cdouble* SubVector	(	cdouble *	z,
		cdouble *	x,
		int	X1,
		int	X,
		MList *	List
	)

function SubVector: extract a subvector of x at X1 to z[X].

In: vector x no shorter than X+X1. Out: vector z[X].

Returns pointer to z. z is created anew if z=0 is specifid on start.

void transpose	(	int	N,
		double **	a
	)

function transpose: matrix transpose: A'->A

In: matrix a[N][N] Out: matrix a[N][N] after transpose

No return value.

double** transpose	(	int	N,
		int	M,
		double **	ta,
		double **	a,
		MList *	List
	)

function transpose: matrix transpose: A'->Z

In: matrix a[M][N] Out: matrix z[N][M]

Returns pointer to z. z is created anew if z=0 is specifid on start.

double** Unitary	(	int	N,
		double **	P,
		double *	x,
		double *	y,
		MList *	List
	)

function Unitary: given x & y s.t. |x|=|y|, find unitary matrix P s.t. Px=y. P is given in closed form as I-(x-y)(x-y)'/(x-y)'x

In: vectors x[N] and y[N] Out: matrix P[N][N]

Returns pointer to P. P is created anew if P=0 is specified on start.

HSM and tools

Macros

Functions

Variables

Detailed Description

Macro Definition Documentation

Function Documentation