view Matrix.h @ 1:6422640a802f

first upload

#ifndef MatrixH
#define MatrixH

/*
  Matrix.cpp - matrix operations.

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

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

//--Similar balance----------------------------------------------------------
void BalanceSim(int n, double** A);

//--Choleski factorization---------------------------------------------------
int Choleski(int N, double** L, double** A);

//--matrix copy--------------------------------------------------------------
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);

//--matrix determinant calculation-------------------------------------------
double det(int N, double** A, int mode=0);
cdouble det(int N, cdouble** A, int mode=0);

//--power methods for solving eigenproblems----------------------------------
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);

//--Gaussian elimination for solving linear systems--------------------------
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 GECP: Gaussian elimination with maximal column pivoting for solving linear system Ax=b

  In: 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.
*/
template<class T, class Ta>int GECP(int N, T* x, Ta** A, T *b, Ta* logdet=0)
{
	if (logdet) *logdet=1E-302; int c, p, ip, *rp=new int[N]; for (int i=0; i<N; i++) rp[i]=i;
  Ta m;
  //Gaussian eliminating
  for (int i=0; i<N-1; i++)
  {
    p=i, ip=i+1;
    while (ip<N){if (fabs(A[rp[ip]][i])>fabs(A[rp[p]][i])) p=ip; ip++;}
    if (A[rp[p]][i]==0) {delete[] rp; return 1;}
    if (p!=i) {c=rp[i]; rp[i]=rp[p]; rp[p]=c;}
    for (int j=i+1; j<N; j++)
    {
      m=A[rp[j]][i]/A[rp[i]][i];
      A[rp[j]][i]=0;
      for (int k=i+1; k<N; k++) A[rp[j]][k]-=m*A[rp[i]][k];
      b[rp[j]]-=m*b[rp[i]];
    }
  }
  if (A[rp[N-1]][N-1]==0) {delete[] rp; return 1;}
  //backward substitution
  x[N-1]=b[rp[N-1]]/A[rp[N-1]][N-1];
  for (int i=N-2; i>=0; i--)
  {
    x[i]=b[rp[i]]; for (int j=i+1; j<N; j++) x[i]-=A[rp[i]][j]*x[j]; x[i]/=A[rp[i]][i];
  }
  if (logdet){*logdet=log(fabs(A[rp[0]][0])); for (int n=1; n<N; n++) *logdet+=log(fabs(A[rp[n]][n]));}
	delete[] rp;
	return 0;
}//GECP

//--inverse lower and upper triangular matrices------------------------------
double GILT(int N, double** A);
double GIUT(int N, double** A);

//--inverse matrix calculation with gaussian elimination---------------------
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);

//--Gaussian-Seidel method for solving linear systems------------------------
int GSI(int N, double* x0, double** A, double* b, double ep=1e-4, int maxiter=50);

//Reduction to upper Heissenberg matrix by elimination with pivoting
void Hessenb(int n, double** A);

//--Householder algorithm converting a matrix tridiagonal--------------------
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);

//--inner product------------------------------------------------------------
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 Inner: Inner product <xw, y>

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

  Returns inner product of xw and y.
*/
template<class Tx, class Tw>cdouble Inner(int N, Tx* x, Tw* w, cdouble* y)
{
  cdouble result=0;
  for (int i=0; i<N; i++) result+=(x[i]*w[i])**y[i];
  return result;
}//Inner
template<class Tx, class Tw>cdouble Inner(int N, Tx* x, Tw* w, double* y)
{
  cdouble result=0;
  for (int i=0; i<N; i++) result+=x[i]*w[i]*y[i];
  return result;
}//Inner

//--Jacobi iterative method for solving linear systems-----------------------
int JI(int N, double* x, double** A, double* b, double ep=1e-4, int maxiter=50);

//--LDL factorization of a symmetric matrix----------------------------------
int LDL(int N, double** L, double* d, double** A);

//--LQ factorization by Gram-Schmidt-----------------------------------------
void LQ_GS(int M, int N, double** A, double** L, double** Q);

//--1-stage Least-square solution of overdetermined linear system------------
/*
  template LSLinera: direct LS solution of A[M][N]x[N]=y[M], M>=N.

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

  Returns the log determinant of AtA. Contents of matrix A and vector y are unchanged on return.
*/
template <class T> T LSLinear(int M, int N, T* x, T** A, T* y)
{
	MList* mlist=new MList;
	T** AtA=MultiplyXtX(N, M, A, mlist);
	T* Aty=MultiplyXty(N, M, A, y, mlist);
	T result; GECP(N, x, AtA, Aty, &result);
	delete mlist;
	return result;
}//LSLinear

//--2-stage Least-square solution of overdetermined linear system------------
void LSLinear2(int M, int N, double* x, double** A, double* y);

//--LU factorization of a non-singular matrix--------------------------------
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);

//--LU factorization method for solving a linear system----------------------
void LU(int N, double* x, double** A, double* y, int* ind=0);

//--find maximum of vector---------------------------------------------------
int maxind(double* data, int from, int to);

//--matrix linear combination------------------------------------------------
/*
  template MultiAdd: matrix linear combination Z=X+aY

  In: matrices X[M][N], Y[M][N], scalar a
  Out: matrix Z[M][N]

  Returns pointer to Z. Z is created anew if Z=0 is specified on start.
*/
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)
{
  if (!Z){Allocate2(Tz, M, N, Z); if (List) List->Add(Z, 2);}
  for (int i=0; i<M; i++) for (int j=0; j<N; j++) Z[i][j]=X[i][j]+a*Y[i][j];
  return Z;
}//MultiAdd

/*
  template MultiAdd: vector linear combination z=x+ay

  In: vectors x[N], y[N], scalar a
  Out: vector z[N]

  Returns pointer to z. z is created anew if z=0 is specified on start.
*/
template<class Tz, class Tx, class Ty, class Ta> Tz* MultiAdd(int N, Tz* z, Tx* x, Ty* y, Ta a, MList* List=0)
{
  if (!z){z=new Tz[N]; if (List) List->Add(z, 1);}
  for (int i=0; i<N; i++) z[i]=x[i]+Ty(a)*y[i];
  return z;
}//MultiAdd
template<class Tx, class Ty, class Ta> Tx* MultiAdd(int N, Tx* x, Ty* y, Ta a, MList* List=0)
{
  return MultiAdd(N, (Tx*)0, x, y, a, List);
}//MultiAdd

//--matrix multiplication by constant----------------------------------------
/*
  template Multiply: matrix-constant multiplication Z=aX

  In: matrix X[M][N], scalar a
  Out: matrix Z[M][N]

  Returns pointer to Z. Z is created anew if Z=0 is specified on start.
*/
template<class Tz, class Tx, class Ta> Tz** Multiply(int M, int N, Tz** Z, Tx** X, Ta a, MList* List=0)
{
  if (!Z){Allocate2(Tz, M, N, Z); if (List) List->Add(Z, 2);}
  for (int i=0; i<M; i++) for (int j=0; j<N; j++) Z[i][j]=a*X[i][j];
  return Z;
}//Multiply

/*
  template Multiply: vector-constant multiplication z=ax

  In: matrix x[N], scalar a
  Out: matrix z[N]

  Returns pointer to z. z is created anew if z=0 is specified on start.
*/
template<class Tz, class Tx, class Ta> Tz* Multiply(int N, Tz* z, Tx* x, Ta a, MList* List=0)
{
  if (!z){z=new Tz[N]; if (List) List->Add(z, 1);}
  for (int i=0; i<N; i++) z[i]=x[i]*a;
  return z;
}//Multiply
template<class Tx, class Ta>Tx* Multiply(int N, Tx* x, Ta a, MList* List=0)
{
  return Multiply(N, (Tx*)0, x, a, List);
}//Multiply

//--matrix multiplication operations-----------------------------------------
/*
  macro Multiply_full: matrix-matrix multiplication z=xy and multiplication-accumulation z=z+xy.

  Each expansion of the macro defines three function templates; two are named $MULTIPLY and do matrix
  multiplication only; one is named $MULTIADD and accumulates the multiplicated result on top of a
  specified matrix. One of the two $MULTIPLY functions allows using a pre-allocated matrix to accept
  the matrix product, while the other directly allocates a new matrix and returns the pointer.

  Functions are named after their exact functions. For example, MultiplyXYc multiplies matrix X with
  the complex conjugate of matrix Y, where postfix "c" attched to Y stands for conjugate. Likewise,
  the postfix "t" stands for transpose, and "h" stnads for Hermitian (conjugate transpose).

  Three dimension arguments are needed by each function template. The first and last of the three are
  the number of rows and columns, respectively, of the product (output) matrix. The middle one is the
  "other", common dimension of both multiplier matrices.
*/
#define Multiply_full(MULTIPLY, MULTIADD, xx, yy) \
  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;}}
Multiply_full(MultiplyXY, MultiAddXY, x[n][m], y[m][k])     //z[N][K]=x[N][M]y[M][K], identical z and x or y not allowed
Multiply_full(MultiplyXYc, MultiAddXYc, x[n][m], *y[m][k])  //z[N][K]=x[N][M](y*)[M][K], identical z and x or y not allowed
Multiply_full(MultiplyXYt, MultiAddXYt, x[n][m], y[k][m])   //z[N][K]=x[N][M]Yt[M][K], identical z and x or y not allowed
Multiply_full(MultiplyXYh, MultiAddXYh, x[n][m], *y[k][m])  //z[N][K]=x[N][M](Yt*)[M][K], identical z and x or y not allowed

/*
  macro Multiply_square: square matrix multiplication z=xy and multiplication-accumulation z=z+xy.
  Identical z and x allowed for multiplication but not for multiplication-accumulation.

  Each expansion of the macro defines three function templates; two are named $MULTIPLY and do matrix
  multiplication only; one is named $MULTIADD and accumulates the multiplicated result on top of a
  specified matrix. One of the two $MULTIPLY functions allows using a pre-allocated matrix to accept
  the matrix product, while the other directly allocates a new matrix and returns the pointer.
*/
#define Multiply_square(MULTIPLY, MULTIADD) \
  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);}
Multiply_square(MultiplyXY, MultiAddXY)

/*
  macro Multiply_xx: matrix self multiplication z=xx and self-multiplication-accumulation z=z+xx.

  Each expansion of the macro defines three function templates; two are named $MULTIPLY and do matrix
  multiplication only; one is named $MULTIADD and accumulates the multiplicated result on top of a
  specified matrix. One of the two $MULTIPLY functions allows using a pre-allocated matrix to accept
  the matrix product, while the other directly allocates a new matrix and returns the pointer.
*/
#define Multiply_xx(MULTIPLY, MULTIADD, xx1, xx2, zzt) \
  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;}}
Multiply_xx(MultiplyXtX, MultiAddXtX, x[n][m], x[n][k], zz) //z[M][M]=xt[M][N]x[N][M], identical z and x NOT ALLOWED.
Multiply_xx(MultiplyXXt, MultiAddXXt, x[m][n], x[k][n], zz) //z[M][M]=X[M][N]xt[N][M], identical z and x NOT ALLOWED.
Multiply_xx(MultiplyXhX, MultiAddXhX, *x[n][m], x[n][k], *zz) //z[M][M]=(xt*)[M][N]x[N][M], identical z and x NOT ALLOWED.
Multiply_xx(MultiplyXXh, MultiAddXXh, x[m][n], *x[k][n], *zz) //z[M][M]=x[M][N](xt*)[N][M], identical z and x NOT ALLOWED.
Multiply_xx(MultiplyXcXt, MultiAddXcXt, *x[m][n], x[k][n], *zz) //z[M][M]=(x*)[M][N]xt[N][M], identical z and x NOT ALLOWED.

//matrix-vector multiplication routines
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);

//--matrix norms-------------------------------------------------------------
double Norm1(int N, double** A);

//--QR factorization---------------------------------------------------------
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);

//--QR factorization of a tridiagonal matrix---------------------------------
int QL(int N, double* d, double* sd, double** z);
int QR(int N, double **A, cdouble* ev);

//--QU factorization of a matrix---------------------------------------------
void QU(int N, double** Q, double** A);

//--Extract the real part----------------------------------------------------
double** Real(int M, int N, double** z, cdouble** x, MList* List);
double** Real(int M, int N, cdouble** x, MList* List);

//--Finding roots of real polynomials----------------------------------------
int Roots(int N, double* p, double* rr, double* ri);

//--Sor iterative method for solving linear systems--------------------------
int SorI(int N, double* x0, double** a, double* b, double w=1, double ep=1e-4, int maxiter=50);

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

//--matrix transpose operation-----------------------------------------------
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); //z[NM]=a[MN]'
double** transpose(int N, int M, double** a, MList* List=0);

//--rotation matrix converting between vectors-------------------------------
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); 

#endif