comparison 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