diff matrix.cpp @ 5:5f3c32dc6e17

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

--- a/matrix.cpp	Tue Oct 05 17:03:27 2010 +0100
+++ b/matrix.cpp	Wed Oct 06 15:19:49 2010 +0100
@@ -2,8 +2,12 @@
 #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.
 
@@ -60,7 +64,7 @@
 }//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.
@@ -93,7 +97,7 @@
 //---------------------------------------------------------------------------
 //matrix duplication routines
 
-/*
+/**
   function Copy: duplicate the matrix A as matrix Z.
 
   In: matrix A[M][N]
@@ -128,7 +132,7 @@
 //---------------------------------------------------------------------------
 //vector duplication routines
 
-/*
+/**
   function Copy: duplicating vector a as vector z
 
   In: vector a[N]
@@ -153,7 +157,7 @@
 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]
@@ -235,7 +239,7 @@
 }//det
 
 //---------------------------------------------------------------------------
-/*
+/**
   function EigPower: power method for solving dominant eigenvalue and eigenvector
 
   In: matrix A[N][N], initial arbitrary vector x[N].
@@ -268,7 +272,7 @@
 }//EigPower
 
 //---------------------------------------------------------------------------
-/*
+/**
   function EigPowerA: EigPower with Aitken acceleration
 
   In: matrix A[N][N], initial arbitrary vector x[N].
@@ -301,7 +305,7 @@
 }//EigPowerA
 
 //---------------------------------------------------------------------------
-/*
+/**
   function EigPowerI: Inverse power method for solving the eigenvalue given an approximate non-zero
   eigenvector.
 
@@ -342,7 +346,7 @@
 }//EigPowerI
 
 //---------------------------------------------------------------------------
-/*
+/**
   function EigPowerS: symmetric power method for solving the dominant eigenvalue with its eigenvector
 
   In: matrix A[N][N], initial arbitrary vector x[N].
@@ -375,7 +379,7 @@
 }//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).
 
@@ -415,7 +419,7 @@
 //---------------------------------------------------------------------------
 //NR versions of eigensystem
 
-/*
+/**
   function EigenValues: solves for eigenvalues of general system
 
   In: matrix A[N][N]
@@ -430,7 +434,7 @@
   return QR(N, A, ev);
 }//EigenValues
 
-/*
+/**
   function EigSym: Solves real symmetric eigensystem A
 
   In: matrix A[N][N]
@@ -449,7 +453,7 @@
 }//EigSym
 
 //---------------------------------------------------------------------------
-/*
+/**
   function GEB: Gaussian elimination with backward substitution for solving linear system Ax=b.
 
   In: coefficient matrix A[N][N], vector b[N]
@@ -491,7 +495,7 @@
 }//GEB
 
 //---------------------------------------------------------------------------
-/*
+/**
   function GESCP: Gaussian elimination with scaled column pivoting for solving linear system Ax=b
 
   In: matrix A[N][N], vector b[N]
@@ -537,7 +541,7 @@
 }//GESCP
 
 //---------------------------------------------------------------------------
-/*
+/**
   function GExL: solves linear system xL=a, L being lower-triangular. This is used in LU factorization
   for solving linear systems.
 
@@ -556,7 +560,7 @@
   }
 }//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]
@@ -572,7 +576,7 @@
   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
@@ -591,7 +595,7 @@
   }
 }//GExL1
 
-/*
+/**
   function GExL1Add: solves linear system *L=(0, 0, ..., 0, a)', L being lower-triangular, and add the
   solution * to x[].
 
@@ -609,7 +613,7 @@
 }//GExL1Add
 
 //---------------------------------------------------------------------------
-/*
+/**
   function GICP: matrix inverse using Gaussian elimination with column pivoting: inv(A)->A.
 
   In: matrix A[N][N]
@@ -701,7 +705,7 @@
   return result;
 }//GICP
 
-/*
+/**
   function GICP: wrapper function that does not overwrite the input matrix: inv(A)->X.
 
   In: matrix A[N][N]
@@ -716,7 +720,7 @@
 }//GICP
 
 //---------------------------------------------------------------------------
-/*
+/**
   function GILT: inv(lower trangular of A)->lower trangular of A
 
   In: matrix A[N][N]
@@ -742,7 +746,7 @@
   return result;
 }//GILT
 
-/*
+/**
   function GIUT: inv(upper trangular of A)->upper trangular of A
 
   In: matrix A[N][N]
@@ -769,7 +773,7 @@
 }//GIUT
 
 //---------------------------------------------------------------------------
-/*
+/**
   function GISCP: matrix inverse using Gaussian elimination w. scaled column pivoting: inv(A)->A.
 
   In: matrix A[N][N]
@@ -827,7 +831,7 @@
   return result;
 }//GISCP
 
-/*
+/**
   function GISCP: wrapper function that does not overwrite input matrix A: inv(A)->X.
 
   In: matrix A[N][N]
@@ -842,7 +846,7 @@
 }//GISCP
 
 //---------------------------------------------------------------------------
-/*
+/**
   function GSI: Gaussian-Seidel iterative algorithm for solving linear system Ax=b. Breaks down if any
   Aii=0, like the Jocobi method JI(...).
 
@@ -878,7 +882,7 @@
 }//GSI
 
 //---------------------------------------------------------------------------
-/*
+/**
   function Hessenb: reducing a square matrix A to upper Hessenberg form
 
   In: matrix A[N][N]
@@ -933,7 +937,7 @@
 }//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.
 
@@ -978,7 +982,7 @@
   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.
 
@@ -1029,7 +1033,7 @@
   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.
 
@@ -1107,7 +1111,7 @@
 }//HouseHolder
 
 //---------------------------------------------------------------------------
-/*
+/**
   function Inner: inner product z=y'x
 
   In: vectors x[N], y[N]
@@ -1146,7 +1150,7 @@
   return result;
 }//Inner
 
-/*
+/**
   function Inner: inner product z=tr(Y'X)
 
   In: matrices X[M][N], Y[M][N]
@@ -1161,7 +1165,7 @@
 }//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.
 
@@ -1194,7 +1198,7 @@
 }//JI
 
 //---------------------------------------------------------------------------
-/*
+/**
   function LDL: LDL' decomposition A=LDL', where L is lower triangular and D is diagonal identical l and
   a allowed.
 
@@ -1233,7 +1237,7 @@
 }//LDL
 
 //---------------------------------------------------------------------------
-/*
+/**
   function LQ_GS: LQ decomposition using Gram-Schmidt method
 
   In: matrix A[M][N], M<=N
@@ -1262,7 +1266,7 @@
 }//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.
 
@@ -1296,7 +1300,7 @@
 }//LSLinear2
 
 //---------------------------------------------------------------------------
-/*
+/**
   function LU: LU decomposition A=LU, where L is lower triangular with diagonal entries 1 and U is upper
   triangular.
 
@@ -1333,7 +1337,7 @@
   return result;
 }//LU
 
-/*
+/**
   function LU: Solving linear system Ax=y by LU factorization
 
   In: matrix A[N][N], vector y[N]
@@ -1388,7 +1392,7 @@
 } //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,
@@ -1437,7 +1441,7 @@
 }//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]
@@ -1510,7 +1514,7 @@
   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]
@@ -1537,7 +1541,7 @@
 }//LU_PD
 
 //---------------------------------------------------------------------------
-/*
+/**
   function LUCP: LU decomposition A=LU with column pivoting
 
   In: matrix A[N][N]
@@ -1626,7 +1630,7 @@
 }//LUCP
 
 //---------------------------------------------------------------------------
-/*
+/**
   function maxind: returns the index of the maximal value of data[from:(to-1)].
 
   In: vector data containing at least $to entries.
@@ -1686,7 +1690,7 @@
 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]
@@ -1706,7 +1710,7 @@
 }//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
@@ -1776,7 +1780,7 @@
 }//QL
 
 //---------------------------------------------------------------------------
-/*
+/**
   function QR: nr version of QR method for solving upper Hessenberg system A. This is compatible with
   Hessenb method.
 
@@ -1908,7 +1912,7 @@
   return 0;
 }//QR
 
-/*
+/**
   function QR_GS: QR decomposition A=QR using Gram-Schmidt method
 
   In: matrix A[M][N], M>=N
@@ -1936,7 +1940,7 @@
   delete[] u;
 }//QR_GS
 
-/*
+/**
   function QR_householder: QR decomposition using householder transform
 
   In: A[M][N], M>=N
@@ -1974,7 +1978,7 @@
 }//QR_householder
 
 //---------------------------------------------------------------------------
-/*
+/**
   function QU: Unitary decomposition A=QU, where Q is unitary and U is upper triangular
 
   In: matrix A[N][N]
@@ -2004,7 +2008,7 @@
 }//QU
 
 //---------------------------------------------------------------------------
-/*
+/**
   function Real: extracts the real part of matrix X
 
   In: matrix x[M][N];
@@ -2021,7 +2025,7 @@
 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.
@@ -2051,7 +2055,7 @@
 }//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
@@ -2090,7 +2094,7 @@
 //---------------------------------------------------------------------------
 //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]
@@ -2108,7 +2112,7 @@
   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]
@@ -2128,7 +2132,7 @@
   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.
@@ -2149,7 +2153,7 @@
 }//SubVector
 
 //---------------------------------------------------------------------------
-/*
+/**
   function transpose: matrix transpose: A'->A
 
   In: matrix a[N][N]
@@ -2169,7 +2173,7 @@
   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]
@@ -2190,7 +2194,7 @@
 }//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