xue@11: /*
xue@11:     Harmonic sinusoidal modelling and tools
xue@11: 
xue@11:     C++ code package for harmonic sinusoidal modelling and relevant signal processing.
xue@11:     Centre for Digital Music, Queen Mary, University of London.
xue@11:     This file copyright 2011 Wen Xue.
xue@11: 
xue@11:     This program is free software; you can redistribute it and/or
xue@11:     modify it under the terms of the GNU General Public License as
xue@11:     published by the Free Software Foundation; either version 2 of the
xue@11:     License, or (at your option) any later version.
xue@11: */
xue@1: #ifndef MatrixH
xue@1: #define MatrixH
xue@1: 
Chris@5: /**
Chris@5:   \file matrix.h - matrix operations.
xue@1: 
xue@1:   Matrices are accessed by double pointers as MATRIX[Y][X], where Y is the row index.
xue@1: */
xue@1: 
xue@1: #include "xcomplex.h"
xue@1: #include "arrayalloc.h"
xue@1: 
xue@1: //--Similar balance----------------------------------------------------------
xue@1: void BalanceSim(int n, double** A);
xue@1: 
xue@1: //--Choleski factorization---------------------------------------------------
xue@1: int Choleski(int N, double** L, double** A);
xue@1: 
xue@1: //--matrix copy--------------------------------------------------------------
xue@1: double** Copy(int M, int N, double** Z, double** A, MList* List=0);
xue@1:   double** Copy(int M, int N, double** A, MList* List=0);
xue@1:   double** Copy(int N, double** Z, double ** A, MList* List=0);
xue@1:   double** Copy(int N, double ** A, MList* List=0);
xue@1: cdouble** Copy(int M, int N, cdouble** Z, cdouble** A, MList* List=0);
xue@1:   cdouble** Copy(int M, int N, cdouble** A, MList* List=0);
xue@1:   cdouble** Copy(int N, cdouble** Z, cdouble** A, MList* List=0);
xue@1:   cdouble** Copy(int N, cdouble** A, MList* List=0);
xue@1: double* Copy(int N, double* z, double* a, MList* List=0);
xue@1:   double* Copy(int N, double* a, MList* List=0);
xue@1: cdouble* Copy(int N, cdouble* z, cdouble* a, MList* List=0);
xue@1:   cdouble* Copy(int N, cdouble* a, MList* List=0);
xue@1: 
xue@1: //--matrix determinant calculation-------------------------------------------
xue@1: double det(int N, double** A, int mode=0);
xue@1: cdouble det(int N, cdouble** A, int mode=0);
xue@1: 
xue@1: //--power methods for solving eigenproblems----------------------------------
xue@1: int EigPower(int N, double& l, double* x, double** A, double ep=1e-6, int maxiter=50);
xue@1: int EigPowerA(int N, double& l, double* x, double** A, double ep=1e-6, int maxiter=50);
xue@1: int EigPowerI(int N, double& l, double* x, double** A, double ep=1e-6, int maxiter=50);
xue@1: int EigPowerS(int N, double& l, double* x, double** A, double ep=1e-6, int maxiter=50);
xue@1: int EigPowerWielandt(int N, double& m, double* u, double l, double* v, double** A, double ep=1e-06, int maxiter=50);
xue@1: int EigenValues(int N, double** A, cdouble* ev);
xue@1: int EigSym(int N, double** A, double* d, double** Q);
xue@1: 
xue@1: //--Gaussian elimination for solving linear systems--------------------------
xue@1: int GEB(int N, double* x, double** A, double* b);
xue@1: int GESCP(int N, double* x, double** A, double*b);
xue@1: void GExL(int N, double* x, double** L, double* a);
xue@1: void GExLAdd(int N, double* x, double** L, double* a);
xue@1: void GExL1(int N, double* x, double** L, double a);
xue@1: void GExL1Add(int N, double* x, double** L, double a);
xue@1: 
xue@1: /*
xue@1:   template GECP: Gaussian elimination with maximal column pivoting for solving linear system Ax=b
xue@1: 
xue@1:   In: matrix A[N][N], vector b[N]
xue@1:   Out: vector x[N]
xue@1: 
xue@1:   Returns 0 if successful. Contents of matrix A and vector b are destroyed on return.
xue@1: */
xue@1: template<class T, class Ta>int GECP(int N, T* x, Ta** A, T *b, Ta* logdet=0)
xue@1: {
xue@1: 	if (logdet) *logdet=1E-302; int c, p, ip, *rp=new int[N]; for (int i=0; i<N; i++) rp[i]=i;
xue@1:   Ta m;
xue@1:   //Gaussian eliminating
xue@1:   for (int i=0; i<N-1; i++)
xue@1:   {
xue@1:     p=i, ip=i+1;
xue@1:     while (ip<N){if (fabs(A[rp[ip]][i])>fabs(A[rp[p]][i])) p=ip; ip++;}
xue@1:     if (A[rp[p]][i]==0) {delete[] rp; return 1;}
xue@1:     if (p!=i) {c=rp[i]; rp[i]=rp[p]; rp[p]=c;}
xue@1:     for (int j=i+1; j<N; j++)
xue@1:     {
xue@1:       m=A[rp[j]][i]/A[rp[i]][i];
xue@1:       A[rp[j]][i]=0;
xue@1:       for (int k=i+1; k<N; k++) A[rp[j]][k]-=m*A[rp[i]][k];
xue@1:       b[rp[j]]-=m*b[rp[i]];
xue@1:     }
xue@1:   }
xue@1:   if (A[rp[N-1]][N-1]==0) {delete[] rp; return 1;}
xue@1:   //backward substitution
xue@1:   x[N-1]=b[rp[N-1]]/A[rp[N-1]][N-1];
xue@1:   for (int i=N-2; i>=0; i--)
xue@1:   {
xue@1:     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];
xue@1:   }
xue@1:   if (logdet){*logdet=log(fabs(A[rp[0]][0])); for (int n=1; n<N; n++) *logdet+=log(fabs(A[rp[n]][n]));}
xue@1: 	delete[] rp;
xue@1: 	return 0;
xue@1: }//GECP
xue@1: 
xue@1: //--inverse lower and upper triangular matrices------------------------------
xue@1: double GILT(int N, double** A);
xue@1: double GIUT(int N, double** A);
xue@1: 
xue@1: //--inverse matrix calculation with gaussian elimination---------------------
xue@1: double GICP(int N, double** X, double** A);
xue@1: double GICP(int N, double** A);
xue@1: cdouble GICP(int N, cdouble** A);
xue@1: double GISCP(int N, double** X, double** A);
xue@1: double GISCP(int N, double** A);
xue@1: 
xue@1: //--Gaussian-Seidel method for solving linear systems------------------------
xue@1: int GSI(int N, double* x0, double** A, double* b, double ep=1e-4, int maxiter=50);
xue@1: 
xue@1: //Reduction to upper Heissenberg matrix by elimination with pivoting
xue@1: void Hessenb(int n, double** A);
xue@1: 
xue@1: //--Householder algorithm converting a matrix tridiagonal--------------------
xue@1: void HouseHolder(int N, double** A);
xue@1: void HouseHolder(int N, double** T, double** Q, double** A);
xue@1: void HouseHolder(int n, double **A, double* d, double* sd);
xue@1: 
xue@1: //--inner product------------------------------------------------------------
xue@1: double Inner(int N, double* x, double* y);
xue@1: cdouble Inner(int N, double* x, cdouble* y);
xue@1: cdouble Inner(int N, cdouble* x, cdouble* y);
xue@1: cdouble Inner(int N, cfloat* x, cdouble* y);
xue@1: cfloat Inner(int N, cfloat* x, cfloat* y);
xue@1: double Inner(int M, int N, double** X, double** Y);
xue@1: 
xue@1: /*
xue@1:   template Inner: Inner product <xw, y>
xue@1: 
xue@1:   In: vectors x[N], w[N] and y[N]
xue@1: 
xue@1:   Returns inner product of xw and y.
xue@1: */
xue@1: template<class Tx, class Tw>cdouble Inner(int N, Tx* x, Tw* w, cdouble* y)
xue@1: {
xue@1:   cdouble result=0;
xue@1:   for (int i=0; i<N; i++) result+=(x[i]*w[i])**y[i];
xue@1:   return result;
xue@1: }//Inner
xue@1: template<class Tx, class Tw>cdouble Inner(int N, Tx* x, Tw* w, double* y)
xue@1: {
xue@1:   cdouble result=0;
xue@1:   for (int i=0; i<N; i++) result+=x[i]*w[i]*y[i];
xue@1:   return result;
xue@1: }//Inner
xue@1: 
xue@1: //--Jacobi iterative method for solving linear systems-----------------------
xue@1: int JI(int N, double* x, double** A, double* b, double ep=1e-4, int maxiter=50);
xue@1: 
xue@1: //--LDL factorization of a symmetric matrix----------------------------------
xue@1: int LDL(int N, double** L, double* d, double** A);
xue@1: 
xue@1: //--LQ factorization by Gram-Schmidt-----------------------------------------
xue@1: void LQ_GS(int M, int N, double** A, double** L, double** Q);
xue@1: 
xue@1: //--1-stage Least-square solution of overdetermined linear system------------
xue@1: /*
xue@1:   template LSLinera: direct LS solution of A[M][N]x[N]=y[M], M>=N.
xue@1: 
xue@1:   In: matrix A[M][N], vector y[M]
xue@1:   Out: vector x[N]
xue@1: 
xue@1:   Returns the log determinant of AtA. Contents of matrix A and vector y are unchanged on return.
xue@1: */
xue@1: template <class T> T LSLinear(int M, int N, T* x, T** A, T* y)
xue@1: {
xue@1: 	MList* mlist=new MList;
xue@1: 	T** AtA=MultiplyXtX(N, M, A, mlist);
xue@1: 	T* Aty=MultiplyXty(N, M, A, y, mlist);
xue@1: 	T result; GECP(N, x, AtA, Aty, &result);
xue@1: 	delete mlist;
xue@1: 	return result;
xue@1: }//LSLinear
xue@1: 
xue@1: //--2-stage Least-square solution of overdetermined linear system------------
xue@1: void LSLinear2(int M, int N, double* x, double** A, double* y);
xue@1: 
xue@1: //--LU factorization of a non-singular matrix--------------------------------
xue@1: int LU_Direct(int mode, int N, double* diag, double** a);
xue@1: int LU(int N, double** l, double** u, double** a);  void LU_DiagL(int N, double** l, double* diagl, double** u, double** a);
xue@1: int LU_PD(int N, double** b);
xue@1: int LU_PD(int N, double** b, double* c);
xue@1: double LUCP(double **a, int N, int *ind, int& parity, int mode=1);
xue@1: 
xue@1: //--LU factorization method for solving a linear system----------------------
xue@1: void LU(int N, double* x, double** A, double* y, int* ind=0);
xue@1: 
xue@1: //--find maximum of vector---------------------------------------------------
xue@1: int maxind(double* data, int from, int to);
xue@1: 
xue@1: //--matrix linear combination------------------------------------------------
xue@1: /*
xue@1:   template MultiAdd: matrix linear combination Z=X+aY
xue@1: 
xue@1:   In: matrices X[M][N], Y[M][N], scalar a
xue@1:   Out: matrix Z[M][N]
xue@1: 
xue@1:   Returns pointer to Z. Z is created anew if Z=0 is specified on start.
xue@1: */
xue@1: 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)
xue@1: {
xue@1:   if (!Z){Allocate2(Tz, M, N, Z); if (List) List->Add(Z, 2);}
xue@1:   for (int i=0; i<M; i++) for (int j=0; j<N; j++) Z[i][j]=X[i][j]+a*Y[i][j];
xue@1:   return Z;
xue@1: }//MultiAdd
xue@1: 
xue@1: /*
xue@1:   template MultiAdd: vector linear combination z=x+ay
xue@1: 
xue@1:   In: vectors x[N], y[N], scalar a
xue@1:   Out: vector z[N]
xue@1: 
xue@1:   Returns pointer to z. z is created anew if z=0 is specified on start.
xue@1: */
xue@1: template<class Tz, class Tx, class Ty, class Ta> Tz* MultiAdd(int N, Tz* z, Tx* x, Ty* y, Ta a, MList* List=0)
xue@1: {
xue@1:   if (!z){z=new Tz[N]; if (List) List->Add(z, 1);}
xue@1:   for (int i=0; i<N; i++) z[i]=x[i]+Ty(a)*y[i];
xue@1:   return z;
xue@1: }//MultiAdd
xue@1: template<class Tx, class Ty, class Ta> Tx* MultiAdd(int N, Tx* x, Ty* y, Ta a, MList* List=0)
xue@1: {
xue@1:   return MultiAdd(N, (Tx*)0, x, y, a, List);
xue@1: }//MultiAdd
xue@1: 
xue@1: //--matrix multiplication by constant----------------------------------------
xue@1: /*
xue@1:   template Multiply: matrix-constant multiplication Z=aX
xue@1: 
xue@1:   In: matrix X[M][N], scalar a
xue@1:   Out: matrix Z[M][N]
xue@1: 
xue@1:   Returns pointer to Z. Z is created anew if Z=0 is specified on start.
xue@1: */
xue@1: template<class Tz, class Tx, class Ta> Tz** Multiply(int M, int N, Tz** Z, Tx** X, Ta a, MList* List=0)
xue@1: {
xue@1:   if (!Z){Allocate2(Tz, M, N, Z); if (List) List->Add(Z, 2);}
xue@1:   for (int i=0; i<M; i++) for (int j=0; j<N; j++) Z[i][j]=a*X[i][j];
xue@1:   return Z;
xue@1: }//Multiply
xue@1: 
xue@1: /*
xue@1:   template Multiply: vector-constant multiplication z=ax
xue@1: 
xue@1:   In: matrix x[N], scalar a
xue@1:   Out: matrix z[N]
xue@1: 
xue@1:   Returns pointer to z. z is created anew if z=0 is specified on start.
xue@1: */
xue@1: template<class Tz, class Tx, class Ta> Tz* Multiply(int N, Tz* z, Tx* x, Ta a, MList* List=0)
xue@1: {
xue@1:   if (!z){z=new Tz[N]; if (List) List->Add(z, 1);}
xue@1:   for (int i=0; i<N; i++) z[i]=x[i]*a;
xue@1:   return z;
xue@1: }//Multiply
xue@1: template<class Tx, class Ta>Tx* Multiply(int N, Tx* x, Ta a, MList* List=0)
xue@1: {
xue@1:   return Multiply(N, (Tx*)0, x, a, List);
xue@1: }//Multiply
xue@1: 
xue@1: //--matrix multiplication operations-----------------------------------------
xue@1: /*
xue@1:   macro Multiply_full: matrix-matrix multiplication z=xy and multiplication-accumulation z=z+xy.
xue@1: 
xue@1:   Each expansion of the macro defines three function templates; two are named $MULTIPLY and do matrix
xue@1:   multiplication only; one is named $MULTIADD and accumulates the multiplicated result on top of a
xue@1:   specified matrix. One of the two $MULTIPLY functions allows using a pre-allocated matrix to accept
xue@1:   the matrix product, while the other directly allocates a new matrix and returns the pointer.
xue@1: 
xue@1:   Functions are named after their exact functions. For example, MultiplyXYc multiplies matrix X with
xue@1:   the complex conjugate of matrix Y, where postfix "c" attched to Y stands for conjugate. Likewise,
xue@1:   the postfix "t" stands for transpose, and "h" stnads for Hermitian (conjugate transpose).
xue@1: 
xue@1:   Three dimension arguments are needed by each function template. The first and last of the three are
xue@1:   the number of rows and columns, respectively, of the product (output) matrix. The middle one is the
xue@1:   "other", common dimension of both multiplier matrices.
xue@1: */
xue@1: #define Multiply_full(MULTIPLY, MULTIADD, xx, yy) \
xue@1:   template<class Tx, class Ty>Tx** MULTIPLY(int N, int M, int K, Tx** z, Tx** x, Ty** y, MList* List=0){ \
xue@1: 		if (!z) {Allocate2(Tx, N, K, z); if (List) List->Add(z, 2);} \
xue@1:     for (int n=0; n<N; n++) for (int k=0; k<K; k++){ \
xue@1:       Tx zz=0; for (int m=0; m<M; m++) zz+=xx*yy; z[n][k]=zz;} return z;} \
xue@1:   template<class Tx, class Ty>Tx** MULTIPLY(int N, int M, int K, Tx** x, Ty** y, MList* List=0){ \
xue@1:     Tx** Allocate2(Tx, N, K, z); if (List) List->Add(z, 2); \
xue@1:     for (int n=0; n<N; n++) for (int k=0; k<K; k++){ \
xue@1:       Tx zz=0; for (int m=0; m<M; m++) zz+=xx*yy; z[n][k]=zz;} return z;} \
xue@1:   template<class Tx, class Ty>void MULTIADD(int N, int M, int K, Tx** z, Tx** x, Ty** y){ \
xue@1:     for (int n=0; n<N; n++) for (int k=0; k<K; k++){ \
xue@1:       Tx zz=0; for (int m=0; m<M; m++) zz+=xx*yy; z[n][k]+=zz;}}
xue@1: 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
xue@1: 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
xue@1: 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
xue@1: 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
xue@1: 
xue@1: /*
xue@1:   macro Multiply_square: square matrix multiplication z=xy and multiplication-accumulation z=z+xy.
xue@1:   Identical z and x allowed for multiplication but not for multiplication-accumulation.
xue@1: 
xue@1:   Each expansion of the macro defines three function templates; two are named $MULTIPLY and do matrix
xue@1:   multiplication only; one is named $MULTIADD and accumulates the multiplicated result on top of a
xue@1:   specified matrix. One of the two $MULTIPLY functions allows using a pre-allocated matrix to accept
xue@1:   the matrix product, while the other directly allocates a new matrix and returns the pointer.
xue@1: */
xue@1: #define Multiply_square(MULTIPLY, MULTIADD) \
xue@1:   template<class T>T** MULTIPLY(int N, T** z, T** x, T** y, MList* List=0){ \
xue@1:     if (!z){Allocate2(T, N, N, z); if (List) List->Add(z, 2);} \
xue@1:     if (z!=x) MULTIPLY(N, N, N, z, x, y); else{ \
xue@1:       T* tmp=new T[N]; int sizeN=sizeof(T)*N; for (int i=0; i<N; i++){ \
xue@1:         MULTIPLY(1, N, N, &tmp, &x[i], y); memcpy(z[i], tmp, sizeN);} delete[] tmp;} \
xue@1:     return z;} \
xue@1:   template<class T>T** MULTIPLY(int N, T** x, T** y, MList* List=0){return MULTIPLY(N, N, N, x, y, List);}\
xue@1:   template<class T>void MULTIADD(int N, T** z, T** x, T** y){MULTIADD(N, N, N, z, x, y);}
xue@1: Multiply_square(MultiplyXY, MultiAddXY)
xue@1: 
xue@1: /*
xue@1:   macro Multiply_xx: matrix self multiplication z=xx and self-multiplication-accumulation z=z+xx.
xue@1: 
xue@1:   Each expansion of the macro defines three function templates; two are named $MULTIPLY and do matrix
xue@1:   multiplication only; one is named $MULTIADD and accumulates the multiplicated result on top of a
xue@1:   specified matrix. One of the two $MULTIPLY functions allows using a pre-allocated matrix to accept
xue@1:   the matrix product, while the other directly allocates a new matrix and returns the pointer.
xue@1: */
xue@1: #define Multiply_xx(MULTIPLY, MULTIADD, xx1, xx2, zzt) \
xue@1:   template<class T>T** MULTIPLY(int M, int N, T** z, T** x, MList* List=0){ \
xue@1:     if (!z){Allocate2(T, M, M, z); if (List) List->Add(z, 2);} \
xue@1:     for (int m=0; m<M; m++) for (int k=0; k<=m; k++){ \
xue@1:       T zz=0; for (int n=0; n<N; n++) zz+=xx1*xx2; z[m][k]=zz; if (m!=k) z[k][m]=zzt;} \
xue@1:     return z;} \
xue@1:   template<class T>T** MULTIPLY(int M, int N, T ** x, MList* List=0){ \
xue@1:     T** Allocate2(T, M, M, z); if (List) List->Add(z, 2); \
xue@1:     for (int m=0; m<M; m++) for (int k=0; k<=m; k++){ \
xue@1:       T zz=0; for (int n=0; n<N; n++) zz+=xx1*xx2; z[m][k]=zz; if (m!=k) z[k][m]=zzt;} \
xue@1:     return z;} \
xue@1:   template<class T>void MULTIADD(int M, int N, T** z, T** x){ \
xue@1:     for (int m=0; m<M; m++) for (int k=0; k<=m; k++){ \
xue@1:       T zz=0; for (int n=0; n<N; n++) zz+=xx1*xx2; z[m][k]+=zz; if (m!=k) z[k][m]+=zzt;}}
xue@1: 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.
xue@1: 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.
xue@1: 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.
xue@1: 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.
xue@1: 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.
xue@1: 
xue@1: //matrix-vector multiplication routines
xue@1: double* MultiplyXy(int M, int N, double* z, double** x, double* y, MList* List=0);
xue@1: double* MultiplyXy(int M, int N, double** x, double* y, MList* List=0);
xue@1: cdouble* MultiplyXy(int M, int N, cdouble* z, cdouble** x, cdouble* y, MList* List=0);
xue@1: cdouble* MultiplyXy(int M, int N, cdouble** x, cdouble* y, MList* List=0);
xue@1: cdouble* MultiplyXy(int M, int N, cdouble* z, double** x, cdouble* y, MList* List=0);
xue@1: cdouble* MultiplyXy(int M, int N, double** x, cdouble* y, MList* List=0);
xue@1: double* MultiplyxY(int M, int N, double* zt, double* xt, double** y, MList* List=0);
xue@1: double* MultiplyxY(int M, int N, double* xt, double** y, MList* List=0);
xue@1: cdouble* MultiplyxY(int M, int N, cdouble* zt, cdouble* xt, cdouble** y, MList* List=0);
xue@1: cdouble* MultiplyxY(int M, int N, cdouble* xt, cdouble** y, MList* List=0);
xue@1: double* MultiplyXty(int M, int N, double* z, double** x, double* y, MList* List=0);
xue@1: double* MultiplyXty(int M, int N, double** x, double* y, MList* List=0);
xue@1: cdouble* MultiplyXty(int M, int N, cdouble* z, cdouble** x, cdouble* y, MList* List=0);
xue@1: cdouble* MultiplyXty(int M, int N, cdouble** x, cdouble* y, MList* List=0);
xue@1: cdouble* MultiplyXhy(int M, int N, cdouble* z, cdouble** x, cdouble* y, MList* List=0);
xue@1: cdouble* MultiplyXhy(int M, int N, cdouble** x, cdouble* y, MList* List=0);
xue@1: cdouble* MultiplyXcy(int M, int N, cdouble* z, cdouble** x, cfloat* y, MList* List=0);
xue@1: cdouble* MultiplyXcy(int M, int N, cdouble** x, cfloat* y, MList* List=0);
xue@1: cdouble* MultiplyXcy(int M, int N, cdouble* z, cdouble** x, cdouble* y, MList* List=0);
xue@1: cdouble* MultiplyXcy(int M, int N, cdouble** x, cdouble* y, MList* List=0);
xue@1: double* MultiplyxYt(int M, int N, double* zt, double* xt, double** y, MList* MList=0);
xue@1: double* MultiplyxYt(int M, int N, double* xt, double** y, MList* MList=0);
xue@1: 
xue@1: //--matrix norms-------------------------------------------------------------
xue@1: double Norm1(int N, double** A);
xue@1: 
xue@1: //--QR factorization---------------------------------------------------------
xue@1: void QR_GS(int M, int N, double** A, double** Q, double** R);
xue@1: void QR_householder(int M, int N, double** A, double** Q, double** R);
xue@1: 
xue@1: //--QR factorization of a tridiagonal matrix---------------------------------
xue@1: int QL(int N, double* d, double* sd, double** z);
xue@1: int QR(int N, double **A, cdouble* ev);
xue@1: 
xue@1: //--QU factorization of a matrix---------------------------------------------
xue@1: void QU(int N, double** Q, double** A);
xue@1: 
xue@1: //--Extract the real part----------------------------------------------------
xue@1: double** Real(int M, int N, double** z, cdouble** x, MList* List);
xue@1: double** Real(int M, int N, cdouble** x, MList* List);
xue@1: 
xue@1: //--Finding roots of real polynomials----------------------------------------
xue@1: int Roots(int N, double* p, double* rr, double* ri);
xue@1: 
xue@1: //--Sor iterative method for solving linear systems--------------------------
xue@1: int SorI(int N, double* x0, double** a, double* b, double w=1, double ep=1e-4, int maxiter=50);
xue@1: 
xue@1: //--Submatrix----------------------------------------------------------------
xue@1: void SetSubMatrix(double** z, double** x, int Y1, int Y, int X1, int X);
xue@1: void SetSubMatrix(cdouble** z, cdouble** x, int Y1, int Y, int X1, int X);
xue@1: cdouble** SubMatrix(cdouble** z, cdouble** x, int Y1, int Y, int X1, int X, MList* List=0);
xue@1: cdouble** SubMatrix(cdouble** x, int Y1, int Y, int X1, int X, MList* List=0);
xue@1: cdouble* SubVector(cdouble* z, cdouble* x, int X1, int X, MList* List=0);
xue@1: cdouble* SubVector(cdouble* x, int X1, int X, MList* List=0);
xue@1: 
xue@1: //--matrix transpose operation-----------------------------------------------
xue@1: void transpose(int N, double** a);
xue@1: void transpose(int N, cdouble** a);
xue@1: double** transpose(int N, int M, double** ta, double** a, MList* List=0); //z[NM]=a[MN]'
xue@1: double** transpose(int N, int M, double** a, MList* List=0);
xue@1: 
xue@1: //--rotation matrix converting between vectors-------------------------------
xue@1: double** Unitary(int N, double** P, double* x, double* y, MList* List=0);
xue@1: double** Unitary(int N, double* x, double* y, MList* List=0);
xue@1: cdouble** Unitary(int N, cdouble** P, cdouble* x, cdouble* y, MList* List=0);
xue@1: cdouble** Unitary(int N, cdouble* x, cdouble* y, MList* List=0); 
xue@1: 
xue@1: #endif