comparison matrix.cpp @ 5:5f3c32dc6e17

* Adjust comment syntax to permit Doxygen to generate HTML documentation; add Doxyfile

 //---------------------------------------------------------------------------
 #include <math.h>
 #include <memory.h>
 #include "matrix.h"
-//---------------------------------------------------------------------------
-/*
+
+/** \file matrix.h */
+
+//---------------------------------------------------------------------------
+
+/**
   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

     }
   }
 }//BalanceSim
 
 //---------------------------------------------------------------------------
-/*
+/**
   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]

 }//Choleski
 
 //---------------------------------------------------------------------------
 //matrix duplication routines
 
-/*
+/**
   function Copy: duplicate the matrix A as matrix Z.
 
   In: matrix A[M][N]
   Out: matrix Z[M][N]
 

 cdouble** Copy(int N, cdouble** A, MList* List){return Copy(N, N, 0, A, List);}
 
 //---------------------------------------------------------------------------
 //vector duplication routines
 
-/*
+/**
   function Copy: duplicating vector a as vector z
 
   In: vector a[N]
   Out: vector z[N]
 

 //version without specifying pre-allocated z
 double* Copy(int N, double* a, MList* List){return Copy(N, 0, a, List);}
 cdouble* Copy(int N, cdouble* a, MList* List){return Copy(N, 0, a, List);}
 
 //---------------------------------------------------------------------------
-/*
+/**
   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.

   delete[] rp;
   return result;
 }//det
 
 //---------------------------------------------------------------------------
-/*
+/**
   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].
 

   }
   delete[] y; return 1;
 }//EigPower
 
 //---------------------------------------------------------------------------
-/*
+/**
   function EigPowerA: EigPower with Aitken acceleration
 
   In: matrix A[N][N], initial arbitrary vector x[N].
   Out: eigenvalue l, eigenvector x[N].
 

   }
   delete[] y; return 1;
 }//EigPowerA
 
 //---------------------------------------------------------------------------
-/*
+/**
   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].

   delete[] aa[0]; delete[] aa;
   delete[] y; return 1;
 }//EigPowerI
 
 //---------------------------------------------------------------------------
-/*
+/**
   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].
 

   delete[] y;
   return 1;
 }//EigPowerS
 
 //---------------------------------------------------------------------------
-/*
+/**
   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

 }//EigPowerWielandt
 
 //---------------------------------------------------------------------------
 //NR versions of eigensystem
 
-/*
+/**
   function EigenValues: solves for eigenvalues of general system
 
   In: matrix A[N][N]
   Out: eigenvalues ev[N]
 

   BalanceSim(N, A);
   Hessenb(N, A);
   return QR(N, A, ev);
 }//EigenValues
 
-/*
+/**
   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)
 

   delete[] t;
   return result;
 }//EigSym
 
 //---------------------------------------------------------------------------
-/*
+/**
   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]
 

   delete[] rp;
   return 0;
 }//GEB
 
 //---------------------------------------------------------------------------
-/*
+/**
   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]
 

   delete[] s; delete[] rp;
   return 0;
 }//GESCP
 
 //---------------------------------------------------------------------------
-/*
+/**
   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]

     for (int m=n+1; m<N; m++) xn-=x[m]*L[m][n];
     x[n]=xn/L[n][n];
   }
 }//GExL
 
-/*
+/**
   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]
 

   GExL(N, lx, L, a);
   for (int i=0; i<N; i++) x[i]+=lx[i];
   delete[] lx;
 }//GExLAdd
 
-/*
+/**
   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]
 

     x[n]=xn/L[n][n];
     xn=0;
   }
 }//GExL1
 
-/*
+/**
   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]

   for (int i=0; i<N; i++) x[i]+=lx[i];
   delete[] lx;
 }//GExL1Add
 
 //---------------------------------------------------------------------------
-/*
+/**
   function GICP: matrix inverse using Gaussian elimination with column pivoting: inv(A)->A.
 
   In: matrix A[N][N]
   Out: matrix A[N][N]
 

 
   delete[] tm; delete[] rp;
   return result;
 }//GICP
 
-/*
+/**
   function GICP: wrapper function that does not overwrite the input matrix: inv(A)->X.
 
   In: matrix A[N][N]
   Out: matrix X[N][N]
 

   Copy(N, X, A);
   return GICP(N, X);
 }//GICP
 
 //---------------------------------------------------------------------------
-/*
+/**
   function GILT: inv(lower trangular of A)->lower trangular of A
 
   In: matrix A[N][N]
   Out: matrix A[N][N]
 

     }
   }
   return result;
 }//GILT
 
-/*
+/**
   function GIUT: inv(upper trangular of A)->upper trangular of A
 
   In: matrix A[N][N]
   Out: matrix A[N][N]
 

   }
   return result;
 }//GIUT
 
 //---------------------------------------------------------------------------
-/*
+/**
   function GISCP: matrix inverse using Gaussian elimination w. scaled column pivoting: inv(A)->A.
 
   In: matrix A[N][N]
   Out: matrix A[N][N]
 

 
   delete[] tm; delete[] s; delete[] rp;
   return result;
 }//GISCP
 
-/*
+/**
   function GISCP: wrapper function that does not overwrite input matrix A: inv(A)->X.
 
   In: matrix A[N][N]
   Out: matrix X[N][N]
 

   Copy(N, X, A);
   return GISCP(N, X);
 }//GISCP
 
 //---------------------------------------------------------------------------
-/*
+/**
   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.

   if (k>=maxiter) return 1;
   return 0;
 }//GSI
 
 //---------------------------------------------------------------------------
-/*
+/**
   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
 

     }
   }
 }//Hessenb
 
 //---------------------------------------------------------------------------
-/*
+/**
   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

     A[k][k+1]=A[k+1][k]=A[k+1][k]-v[k+1]*z[k];
   }
   delete[] u; delete[] v; delete[] z;
 }//HouseHolder
 
-/*
+/**
   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]

       MultiAdd(N-k, &Q[i][k], &Q[i][k], &v[k], -Inner(N-k, &Q[i][k], &v[k])/r2);
   }
   delete[] u; delete[] v; delete[] z;
 }//HouseHolder
 
-/*
+/**
   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

     for (int j=0; j<=l; j++) A[j][i]=A[i][j]=0.0;
   }
 }//HouseHolder
 
 //---------------------------------------------------------------------------
-/*
+/**
   function Inner: inner product z=y'x
 
   In: vectors x[N], y[N]
 
   Returns inner product of x and y.

   cfloat result=0;
   for (int i=0; i<N; i++) result+=x[i]^y[i];
   return result;
 }//Inner
 
-/*
+/**
   function Inner: inner product z=tr(Y'X)
 
   In: matrices X[M][N], Y[M][N]
 
   Returns inner product of X and Y.

   for (int m=0; m<M; m++) for (int n=0; n<N; n++) result+=X[m][n]*Y[m][n];
 	return result;
 }//Inner
 
 //---------------------------------------------------------------------------
-/*
+/**
   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.

   if (k>=maxiter) return 1;
   else return 0;
 }//JI
 
 //---------------------------------------------------------------------------
-/*
+/**
   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.

   for (int i=0; i<N; i++) {L[i][i]=1; memset(&L[i][i+1], 0, sizeof(double)*(N-1-i));}
   return 0;
 }//LDL
 
 //---------------------------------------------------------------------------
-/*
+/**
   function LQ_GS: LQ decomposition using Gram-Schmidt method
 
   In: matrix A[M][N], M<=N
   Out: matrices L[M][M], Q[M][N]
 

   }
   delete[] u;
 }//LQ_GS
 
 //---------------------------------------------------------------------------
-/*
+/**
   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].

   }
   DeAlloc2(A1);
 }//LSLinear2
 
 //---------------------------------------------------------------------------
-/*
+/**
   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.
 

   }
   delete[] diagl;
   return result;
 }//LU
 
-/*
+/**
   function LU: Solving linear system Ax=y by LU factorization
 
   In: matrix A[N][N], vector y[N]
   Out: x[N]
 

   u[N-1][N-1]=a[N-1][N-1]; for (int k=0; k<N-1; k++) u[N-1][N-1]-=l[N-1][k]*u[k][N-1]; u[N-1][N-1]/=l[N-1][N-1];
   memset(u[N-1], 0, sizeof(double)*(N-1));
 } //LU_DiagL*/
 
 //---------------------------------------------------------------------------
-/*
+/**
   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.

   }
   return 0;
 }//LU_Direct
 
 //---------------------------------------------------------------------------
-/*
+/**
   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]

   }
   for (int k=N-3; k<N-1; k++) b[0][N-1]-=b[N-1-k][k]*b[k-N+1][N-1];
   return 0;
 }//LU_PD*/
 
-/*
+/**
   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.
 

   }
   return result;
 }//LU_PD
 
 //---------------------------------------------------------------------------
-/*
+/**
   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:

   delete[] norm;
   return det*parity;
 }//LUCP
 
 //---------------------------------------------------------------------------
-/*
+/**
   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)]
 

 //function MultiplyXcy: z[M]=(x*)[M][N]y[N]
 Multiply_vect(MultiplyXcy, cdouble, cdouble*, cdouble, *x[m][n], y[n])
 Multiply_vect(MultiplyXcy, cdouble, cdouble*, cfloat, *x[m][n], y[n])
 
 //---------------------------------------------------------------------------
-/*
+/**
   function Norm1: L-1 norm of a square matrix A
 
   In: matrix A[N][N]
   Out: its L-1 norm
 

   }
   return result;
 }//Norm1
 
 //---------------------------------------------------------------------------
-/*
+/**
   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.

   }
   return 0;
 }//QL
 
 //---------------------------------------------------------------------------
-/*
+/**
   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

     } while (n>l+1);
   }
   return 0;
 }//QR
 
-/*
+/**
   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]
 

     iu=1.0/iu; for (int m=0; m<M; m++) Q[m][n]=u[m]*iu;
   }
   delete[] u;
 }//QR_GS
 
-/*
+/**
   function QR_householder: QR decomposition using householder transform
 
   In: A[M][N], M>=N
   Out: Q[M][M], R[M][N]
 

   }
   delete[] u;
 }//QR_householder
 
 //---------------------------------------------------------------------------
-/*
+/**
   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
 

   delete[] tmpi; delete[] tmpj;
   transpose(N, Q);
 }//QU
 
 //---------------------------------------------------------------------------
-/*
+/**
   function Real: extracts the real part of matrix X
 
   In: matrix x[M][N];
   Out: matrix z[M][N]
 

   return z;
 }//Real
 double** Real(int M, int N, cdouble** x, MList* List){return Real(M, N, 0, x, List);}
 
 //---------------------------------------------------------------------------
-/*
+/**
   function Roots: finds the roots of a polynomial. x^N+p[N-1]x^(N-1)+p[N-2]x^(N-2)...+p[0]
 
   In: vector p[N] of polynomial coefficients.
   Out: vector r[N] of roots.
 

   delete[] r;
   return result;
 }//Roots
 
 //---------------------------------------------------------------------------
-/*
+/**
   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

 }//SorI
 
 //---------------------------------------------------------------------------
 //Submatrix routines
 
-/*
+/**
   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.
 

 void SetSubMatrix(cdouble** z, cdouble** x, int Y1, int Y, int X1, int X)
 {
   for (int y=0; y<Y; y++) memcpy(&z[Y1+y][X1], x[y], sizeof(cdouble)*X);
 }//SetSubMatrix
 
-/*
+/**
   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].
 

 cdouble** SubMatrix(cdouble** x, int Y1, int Y, int X1, int X, MList* List)
 {
   return SubMatrix(0, x, Y1, Y, X1, X, List);
 }//SetSubMatrix
 
-/*
+/**
   function SubVector: extract a subvector of x at X1 to z[X].
 
   In: vector x no shorter than X+X1.
   Out: vector z[X].
 

 {
   return SubVector(0, x, X1, X, List);
 }//SubVector
 
 //---------------------------------------------------------------------------
-/*
+/**
   function transpose: matrix transpose: A'->A
 
   In: matrix a[N][N]
   Out: matrix a[N][N] after transpose
 

 {
   cdouble tmp;
   for (int i=1; i<N; i++) for (int j=0; j<i; j++) {tmp=a[i][j]; a[i][j]=a[j][i]; a[j][i]=tmp;}
 }//transpose
 
-/*
+/**
   function transpose: matrix transpose: A'->Z
 
   In: matrix a[M][N]
   Out: matrix z[N][M]
 

 {
   return transpose(N, M, 0, a, List);
 }//transpose
 
 //---------------------------------------------------------------------------
-/*
+/**
   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]