From cb28835fe90529a5160968fa50cf10baa9ad4c08 Mon Sep 17 00:00:00 2001 From: Jiri Pittner Date: Thu, 5 Nov 2026 21:14:13 +0100 Subject: [PATCH] documentation update --- doc/building.tex | 1 - doc/concepts.tex | 1 + doc/mat.tex | 55 ++++++++++++++++++++++++------------------ doc/nonclass.tex | 59 +++++++++++++++++++++++----------------------- doc/polynomial.tex | 6 ++--- doc/smat.tex | 19 +++++++-------- doc/sparse.tex | 14 +++++------ 7 files changed, 82 insertions(+), 73 deletions(-) diff --git a/doc/building.tex b/doc/building.tex index 5456064..be42018 100644 --- a/doc/building.tex +++ b/doc/building.tex @@ -77,7 +77,6 @@ functions do not match. \code{LONG_FORTRAN_INT} & Fortran integers are 64-bit (\code{long}) \\ \code{MATPTR} & row pointer representation of \code{NRMat} \\ \code{CUDALA} & GPU support via CUBLAS (experimental) \\ -\code{NO_STRASSEN} & do not compile the experimental Strassen matrix multiplication interface (\header{strassen.cc}, needs an external library); always set in \header{Makefile.am} \\ \bottomrule \end{longtable} diff --git a/doc/concepts.tex b/doc/concepts.tex index d4f7df2..38a7e04 100644 --- a/doc/concepts.tex +++ b/doc/concepts.tex @@ -145,6 +145,7 @@ Its main members are: \lafunc{static void put(int fd, const T &x, ...); static void get(int fd, T &x, ...); multiput; multiget}{binary I/O, recursive for nested types} \lafunc{static T conjugate(const T &x); realpart(x); imagpart(x); setrealpart(x,y); setimagpart(x,y)}{complex helpers, trivial for real types} \lafunc{static bool bigger(const T &x, const T &y); smaller(x,y)}{comparison used for sorting; for complex types it is an error} +\lafunc{static int strassen_cutoff()}{scalar types only: default \code{cutoff} of \code{NRMat::strassen} (section~\ref{sec:mat:strassen}), currently 2048 for all types} \end{lareference} The auxiliary template \code{isscalar} distinguishes scalars from classes, \code{LA_traits_complex} gives for a complex type the component type and the corresponding real vector and matrix types (used by diff --git a/doc/mat.tex b/doc/mat.tex index 3462cfd..6b16b89 100644 --- a/doc/mat.tex +++ b/doc/mat.tex @@ -30,9 +30,9 @@ reference count is larger than one (controlled by \code{_LA_count_check}); read possibly shared matrix through a const reference. \paragraph{Element types.} \code{mat.cc} explicitly instantiates \code{NRMat} for -\code{double}, \code{std::complex}, and the integer types \code{char}, \code{short}, \code{int}, +\code{double}, \code{float}, \code{std::complex}, \code{std::complex}, and the integer types \code{char}, \code{short}, \code{int}, \code{long}, \code{long long} and their unsigned variants. Many numerical members are only -specialized (BLAS-backed) for \code{double} and \code{complex} and have no generic +specialized (BLAS-backed) for the four floating point types and have no generic definition, so using them with an integer type results in a link error. Table~\ref{tab:mat:types} summarizes this. Nested element types (e.g.\ \code{NRMat >}) are supported by the reference counting (type-nested \code{copyonwrite}, sec.~\ref{sec:concepts:traits}), but several @@ -47,10 +47,10 @@ Members & Available for \\ \midrule \texttt{gemm}, \texttt{norm}, \texttt{axpy(T,const NRMat\&)}, \texttt{dot}, \texttt{transpose}, \texttt{diagmultl}, \texttt{diagmultr}, \texttt{operator*(const NRSMat\&)}, \texttt{diagonalset}, \texttt{inverse}, \texttt{svdinverse}, \texttt{absadjoint}, -\texttt{randomize}, \texttt{amax}, \texttt{amin} & \texttt{double}, \texttt{complex} \\ -\texttt{operator*(const NRMat\&)} & \texttt{double}, \texttt{complex}, \texttt{int} \\ -\texttt{diagonalof}, \texttt{orthonormalize}, \texttt{strassen}, \texttt{gemm} with \texttt{SparseMat} & \texttt{double} only \\ -\texttt{nonhermiticity} & \texttt{complex} only \\ +\texttt{randomize}, \texttt{amax}, \texttt{amin}, \texttt{strassen} & \texttt{double}, \texttt{float}, \texttt{complex}, \texttt{complex} \\ +\texttt{operator*(const NRMat\&)} & the same and \texttt{int} \\ +\texttt{diagonalof}, \texttt{orthonormalize}, \texttt{gemm} with \texttt{SparseMat} & \texttt{double}, \texttt{float} \\ +\texttt{nonhermiticity} & \texttt{complex}, \texttt{complex} \\ everything else & all instantiated types \\ \bottomrule \end{tabular} @@ -73,14 +73,23 @@ MATPTR is rarely used and less tested than the default representation; the libra regression test (\code{testoptimize}). \paragraph{Strassen multiplication.}\label{sec:mat:strassen} -Unless \code{NO_STRASSEN} is defined, \code{NRMat} declares \code{strassen()} and \code{s_cutoff()} -(specialized for \code{double} only, in \code{strassen.cc}). They call the external Fortran-compatible -routines \code{fmm} and \code{strassen_cutoff}, which are \emph{not} part of \LA{} -(the test code in \code{t.cc} mentions linking \code{-lstrassen}). Strassen multiplication is never -used implicitly: \code{operator*} and \code{gemm} always call BLAS \code{gemm}. Note that -\code{Makefile.am} unconditionally adds \code{-DNO_STRASSEN}, so in the standard build these members do not exist. -\lawarning{The Strassen interface is \emph{experimental}: it has not been used or tested for a long time. It is kept for possible future work, -in particular if fast (sub-$O(n^3)$) matrix multiplication algorithms become practically useful.} +\code{strassen()} (in \code{strassen.cc}, for \code{double}, \code{float}, \code{complex} and \code{complex}) computes the +same $C=\beta C+\alpha\,\mathrm{op}(A)\,\mathrm{op}(B)$ as \code{gemm} with the same parameters, plus an optional \code{cutoff}, by the +algorithm of Strassen \cite{mat:strassen1969} in the variant of Winograd (7 products and 15 additions of quadrants per level) +\cite{mat:winograd1971}, which needs $O(n^{\log_2 7})=O(n^{2.81})$ operations. The recursion works on strided quadrants of the original +matrices (no copies of $A$ and $B$, transpositions are passed to the products), the products of the blocks are computed by BLAS +\code{gemm}, the order of operations follows the memory-saving schedule with two temporaries per level of \cite{mat:boyer2009} and odd +dimensions are treated by dynamic peeling \cite{mat:husslederman1996}. A level of recursion is done only while all of $m$, $n$, $k$ are larger +than \code{cutoff} (default \code{LA_traits::strassen_cutoff()}, 2048 for all types), so products with a dimension $\le$ +\code{cutoff} are plain \code{gemm} calls. Workspace of about $\frac23 mn$ elements is allocated (plus $mn$ if $\beta\ne0$). +Strassen multiplication is never used implicitly: \code{operator*} and \code{gemm} always use the classical algorithm. + +Whether it pays off depends on the BLAS and the machine: on top of an optimized multithreaded BLAS the saving of one product in eight +is largely eaten by the memory-bound additions, so a gain is to be expected mainly for large matrices ($n\gtrsim4000$) and single-threaded +BLAS; the cutoff should be tuned by a benchmark. CPU only (\code{laerror} for matrices on the GPU); the device-dependent operations +are separated in \code{strassen.cc} for a future CUDA version. +\lawarning{The error of Strassen multiplication is bounded only normwise, not componentwise as for the classical algorithm, and the bound +grows by roughly a factor of 12 per level of recursion. Small elements of the result may therefore be less accurate than with \code{gemm}.} \section{Construction, assignment, conversion}\label{sec:mat:constr} \begin{lareference} @@ -93,7 +102,7 @@ in particular if fast (sub-$O(n^3)$) matrix multiplication algorithms become pra \lafunc{explicit NRMat(const int n) = delete}{A single integer argument is forbidden (it would otherwise silently create a $3\times3$ matrix through \code{Mat3}); use \code{NRMat(n,m)}.} \lafunc{NRMat(const NRMat &rhs)}{Shallow copy (shares data).} -\lafunc{NRMat(const typename LA_traits_complex::NRMat_Noncomplex_type &rhs, bool imagpart = false)}{Complexifying (non-explicit) constructor, defined for \code{T = complex}: \code{rhs} (\code{NRMat}) becomes the real part (or the imaginary part if \code{imagpart}), the other part is zero.} +\lafunc{NRMat(const typename LA_traits_complex::NRMat_Noncomplex_type &rhs, bool imagpart = false)}{Complexifying (non-explicit) constructor, defined for \code{T = complex} and \code{complex}: \code{rhs} (\code{NRMat} resp.\ \code{NRMat}) becomes the real part (or the imaginary part if \code{imagpart}), the other part is zero.} \lafunc{explicit NRMat(const NRMat > &rhs, bool imagpart = false)}{Decomplexifying constructor, defined for \code{T = double}: real part (or imaginary part if \code{imagpart}) of \code{rhs}.} \lafunc{explicit NRMat(const NRSMat &rhs)}{Full square matrix from a packed symmetric one ($A_{ij}=A_{ji}$; for complex types hermitian, $A_{ij}=A_{ji}^*$).} \lafunc{explicit NRMat(const NRVec &rhs, const int n, const int m, const int offset = 0)}{View of the vector data from position \code{offset} as an $n\times m$ matrix. Without MATPTR the data are \emph{shared} (reference count of the vector is incremented); with MATPTR they are copied. A nonzero offset gives a view into the middle of the vector without copying (a deliberate feature); @@ -140,8 +149,8 @@ Note that scalar addition and subtraction act on the \emph{diagonal} only, i.e.\ \lafunc{NRMat & operator^=(const NRMat &rhs)}{Hadamard (element-wise) product $A_{ij}\leftarrow A_{ij}B_{ij}$.} \lafunc{void axpy(const T alpha, const NRMat &x)}{$A\leftarrow A+\alpha X$ (BLAS \code{axpy}).} \lafunc{const NRMat operator*(const NRMat &rhs) const}{Matrix product $AB$ (BLAS \code{gemm}; simple loops for \code{int}).} -\lafunc{void gemm(const T &beta, const NRMat &a, const char transa, const NRMat &b, const char transb, const T &alpha)}{$C\leftarrow\alpha\,\mathrm{op}(A)\,\mathrm{op}(B)+\beta C$ with $C$ = \code{*this}, which must already have the correct dimensions (it is not resized; if $\beta=0$ its contents may be uninitialized). \code{trans}: \code{'n'} no transposition; for \code{double} any other character means transposition; for \code{complex} \code{'t'} transposition and \code{'c'} conjugate transposition (case insensitive). Returns immediately when $\alpha=0$ and $\beta=1$. \code{*this} must not share storage with \code{a} or \code{b}.} -\lafunc{void gemm(const T &beta, const SparseMat &a, const char transa, const NRMat &b, const char transb, const T &alpha)}{The same with a sparse first factor; \code{double} only (the complex version calls \code{laerror("not implemented yet")}), see chapter~\ref{ch:sparse}.} +\lafunc{void gemm(const T &beta, const NRMat &a, const char transa, const NRMat &b, const char transb, const T &alpha)}{$C\leftarrow\alpha\,\mathrm{op}(A)\,\mathrm{op}(B)+\beta C$ with $C$ = \code{*this}, which must already have the correct dimensions (it is not resized; if $\beta=0$ its contents may be uninitialized). \code{trans}: \code{'n'} no transposition; for real types any other character means transposition; for complex types \code{'t'} transposition and \code{'c'} conjugate transposition (case insensitive). Returns immediately when $\alpha=0$ and $\beta=1$. \code{*this} must not share storage with \code{a} or \code{b}.} +\lafunc{void gemm(const T &beta, const SparseMat &a, const char transa, const NRMat &b, const char transb, const T &alpha)}{The same with a sparse first factor; \code{double} and \code{float} only (the complex versions call \code{laerror("not implemented yet")}), see chapter~\ref{ch:sparse}.} \lafunc{const NRMat operator*(const NRSMat &rhs) const}{$AS$ with packed symmetric (real) or hermitian (complex) $S$ (BLAS \code{dspmv}/\code{zhpmv} per row).} \lafunc{SparseSMat operator*(const SparseSMat &rhs) const}{Dense times sparse, the result is sparse (chapter~\ref{ch:sparse}).} \lafunc{const NRVec operator*(const NRVec &rhs) const; const NRVec > operator*(const NRVec > &rhs) const}{Matrix--vector product via \code{NRVec::gemv} (the second form: real matrix times complex vector).} @@ -153,7 +162,7 @@ Note that scalar addition and subtraction act on the \emph{diagonal} only, i.e.\ \lafunc{const NRMat operator&(const NRMat &rhs) const; const NRMat oplus(const NRMat &rhs) const}{Direct sum $A\oplus B$ (block-diagonal matrix). \code{oplus} returns the other operand when one is empty; \code{operator&} is plain-data only.} \lafunc{const NRMat operator|(const NRMat &rhs) const}{Direct (Kronecker) product $A\otimes B$, $(A\otimes B)_{i n_B+k,\,j m_B+l}=A_{ij}B_{kl}$.} \lafunc{const NRMat otimes(const NRMat &rhs, bool reversecolumns = false) const}{Kronecker product as \code{operator|}; with \code{reversecolumns} the column index is ordered the other way, $R_{i n_B+k,\,l m_A+j}=A_{ij}B_{kl}$. Returns the empty operand when one is empty.} -\lafunc{void strassen(const T beta, const NRMat &a, const char transa, const NRMat &b, const char transb, const T alpha); void s_cutoff(const int, const int, const int, const int) const}{Only without \code{NO_STRASSEN}, \code{double} only: Strassen multiplication with \code{gemm} semantics (dimensions always checked) and setting of the cutoffs of the external Strassen library (section~\ref{sec:mat:strassen}).} +\lafunc{void strassen(const T &beta, const NRMat &a, const char transa, const NRMat &b, const char transb, const T &alpha, const int cutoff=LA_traits::strassen_cutoff())}{Strassen-Winograd multiplication with the parameters and semantics of \code{gemm}; blocks with a dimension $\le$ \code{cutoff} are multiplied by \code{gemm} (section~\ref{sec:mat:strassen}).} \end{lareference} \section{Other member functions}\label{sec:mat:other} @@ -167,13 +176,13 @@ Note that scalar addition and subtraction act on the \emph{diagonal} only, i.e.\ \lafunc{const NRMat submatrix(const NRVec &rows, const NRVec &cols) const}{Copy of the elements $A_{r_i c_j}$ for given (arbitrary) lists of row and column indices; indices are always range-checked.} \lafunc{void storesubmatrix(const int fromrow, const int fromcol, const NRMat &rhs)}{Store \code{rhs} as a block with the upper left corner at (\code{fromrow},\code{fromcol}). Plain data only.} \lafunc{void storesubmatrix(const NRVec &rows, const NRVec &cols, const NRMat &rhs)}{Scatter: $A_{r_i c_j}\leftarrow R_{ij}$.} -\lafunc{const T* diagonalof(NRVec &r, const bool divide = 0, bool cache = false) const}{\code{double} only. For a square matrix copy the diagonal into \code{r}, or with \code{divide} divide \code{r} element-wise by the (nonzero) diagonal elements. For a non-square matrix the diagonal of $A\Tr A$ (squared column norms) is used instead; \code{r} must have \code{ncols()} elements. Returns pointer to \code{r}'s data, or NULL if \code{divide}. \code{cache} is unused. Same interface as in other matrix classes, so that it can be used by iterative solvers as a preconditioner (chapter~\ref{ch:iterative}).} +\lafunc{const T* diagonalof(NRVec &r, const bool divide = 0, bool cache = false) const}{\code{double} and \code{float} only. For a square matrix copy the diagonal into \code{r}, or with \code{divide} divide \code{r} element-wise by the (nonzero) diagonal elements. For a non-square matrix the diagonal of $A\Tr A$ (squared column norms) is used instead; \code{r} must have \code{ncols()} elements. Returns pointer to \code{r}'s data, or NULL if \code{divide}. \code{cache} is unused. Same interface as in other matrix classes, so that it can be used by iterative solvers as a preconditioner (chapter~\ref{ch:iterative}).} \lafunc{void diagonalset(const NRVec &r)}{Set the diagonal; \code{r} has $\min(n,m)$ elements.} \end{lareference} \subsection{Transposition and conjugation} \begin{lareference} -\lafunc{const NRMat transpose(bool conj = false) const}{$A\Tr$, or $A\Hc$ if \code{conj} (complex). \code{double} and \code{complex} only.} +\lafunc{const NRMat transpose(bool conj = false) const}{$A\Tr$, or $A\Hc$ if \code{conj} (complex). Floating point types only.} \lafunc{NRMat& transposeme(const int n = 0)}{In-place transposition of a square matrix (no conjugation); with \code{n>0} only the leading $n\times n$ block is transposed (this also works for a non-square matrix with $n\le\min(\text{rows},\text{cols})$). Dimensions are not changed, so a non-square matrix cannot be transposed in place.} \lafunc{NRMat& conjugateme(); const NRMat conjugate() const}{Complex conjugation in place resp.\ by value; no-op for real types.} \end{lareference} @@ -198,7 +207,7 @@ Note that scalar addition and subtraction act on the \emph{diagonal} only, i.e.\ \lafunc{const T sum() const; const NRVec rsum() const; const NRVec csum() const}{Sum of all elements; vector of row sums (length \code{nrows()}, sum over the column index); vector of column sums (length \code{ncols()}).} \lafunc{inline const T amax() const; inline const T amin() const}{The element (not its index) with the largest resp.\ smallest absolute value; for complex the measure is $|\mathrm{Re}|+|\mathrm{Im}|$ (BLAS \code{i?amax} convention).} \lafunc{const typename LA_traits::normtype nonsymmetry() const}{$\sum_{i>j}|A_{ij}-A_{ji}|^2$ (square of the norm, no square root).} -\lafunc{const typename LA_traits::normtype nonhermiticity() const}{\code{complex} only: $\sum_{i\ge j}|A_{ij}-A_{ji}^*|^2$.} +\lafunc{const typename LA_traits::normtype nonhermiticity() const}{Complex types only: $\sum_{i\ge j}|A_{ij}-A_{ji}^*|^2$.} \lafunc{NRMat::normtype> abs() const}{Matrix of element-wise absolute values.} \lafunc{NRMat::normtype> diffabs(const NRMat &rhs) const}{Element-wise $|A_{ij}|-|B_{ij}|$.} \lafunc{const bool operator==(const NRMat &rhs) const; const bool operator!=(const NRMat &rhs) const}{Equality of dimensions and contents; for plain types a bitwise (\code{memcmp}) comparison, element-wise for nested types.} @@ -207,7 +216,7 @@ Note that scalar addition and subtraction act on the \emph{diagonal} only, i.e.\ \subsection{Orthonormalization, inverse, random fill} \begin{lareference} -\lafunc{void orthonormalize(const bool rowcol, const NRSMat *metric = NULL)}{\code{double} only. Modified Gram--Schmidt orthonormalization of the rows (\code{rowcol==true}) or columns (\code{false}) of the matrix in place, optionally with respect to the inner product $x\Tr S y$ given by \code{metric}. A zero vector (or non-positive metric) raises \code{laerror}.} +\lafunc{void orthonormalize(const bool rowcol, const NRSMat *metric = NULL)}{\code{double} and \code{float} only. Modified Gram--Schmidt orthonormalization of the rows (\code{rowcol==true}) or columns (\code{false}) of the matrix in place, optionally with respect to the inner product $x\Tr S y$ given by \code{metric}. A zero vector (or non-positive metric) raises \code{laerror}.} \lafunc{NRMat inverse() const}{$A^{-1}$ via \code{calcinverse} (LAPACK linear solve, chapter~\ref{ch:nonclass}).} \lafunc{NRMat svdinverse(const typename LA_traits::normtype thr = 0) const}{Pseudo-inverse via SVD, singular values below \code{thr} are dropped.} \lafunc{NRMat absadjoint(const typename LA_traits::normtype thr = 0) const}{Adjugate (transposed cofactor) matrix computed via SVD up to the factor $\det U\det V$ (a sign, or a phase in the complex case), cf.\ \code{calcsvdinverse} in chapter~\ref{ch:nonclass}. \code{thr} is not used in this mode.} @@ -322,6 +331,6 @@ NRMat ds = re.oplus(im); // direct sum, 4x4 block diagonal \lawarning{Writing through \code{operator[]}, \code{operator()} or \code{operator T*()} into a matrix obtained by copy construction or \code{=} modifies all sharing matrices unless \code{copyonwrite()} is called first. In DEBUG builds even \emph{reading} through the non-const operators of a shared matrix raises an error; use a const reference.} \lawarning{\code{NRMat(const NRVec&, n, m, offset)} shares the vector data (without MATPTR). With a nonzero offset, if the matrix happens to be the last owner of the data, it will \code{delete[]} a pointer that is not the beginning of the allocated block; the constructor therefore prints a warning. This is by design (a view without copying, avoiding the cost of a copy); the condition of use is: keep the vector alive longer than the matrix, or make an independent copy (\code{|=}) if the matrix has to outlive it.} -\lawarning{Many numerical members (\code{gemm}, \code{transpose}, \code{norm}, \ldots) exist only for \code{double} and \code{complex} (Table~\ref{tab:mat:types}); with other element types the program fails to link.} +\lawarning{Many numerical members (\code{gemm}, \code{transpose}, \code{norm}, \ldots) exist only for the floating point types \code{double}, \code{float}, \code{complex} and \code{complex} (Table~\ref{tab:mat:types}); with other element types the program fails to link.} \lawarning{\code{transposeme()} does not conjugate and works only on square matrices (or square leading blocks); use \code{transpose(true)} for $A\Hc$ or for a non-square matrix.} \lawarning{In the non-DEBUG build no dimension checks are done in most operations (\code{operator+=}, \code{gemm}, \code{operator=(const T&)}, \ldots); mismatches lead to silent memory corruption.} diff --git a/doc/nonclass.tex b/doc/nonclass.tex index 344d7cd..a3fedd8 100644 --- a/doc/nonclass.tex +++ b/doc/nonclass.tex @@ -6,8 +6,8 @@ This chapter describes the free (non-member) functions operating on the dense cl \code{NRVec}, \code{NRMat} and \code{NRSMat}: linear systems, determinants, inverses, Cholesky decomposition, symmetric/hermitian and general eigenvalue problems, singular value decomposition, matrix functions, and a few auxiliary products and conversions. Most of them are thin wrappers of -LAPACK routines (\header{nonclass.cc}); they exist for \code{T}${}={}$\code{double} and -\code{std::complex} only. The templates in \header{matexp.h} (matrix polynomials, matrix +LAPACK routines (\header{nonclass.cc}); they exist for \code{T}${}={}$\code{double}, \code{float}, +\code{std::complex} and \code{std::complex} only. The templates in \header{matexp.h} (matrix polynomials, matrix exponential, $\exp(A)\,x$, BCH expansion) are generic: they work for every class that provides the few operations listed in section~\ref{sec:nonclass:matexp}, including sparse matrices (chapter~\ref{ch:sparse}) and user-defined matrix classes which only implement \code{gemv}. @@ -44,8 +44,8 @@ growing leading block. (sec.~\ref{sec:concepts:errors}), e.g.\ ``singular matrix'' or ``convergence problem''. Exceptions: \code{linear_solve_x} and \code{multiply_by_inverse} return \code{info}. -\paragraph{Element types and availability.} The functions are declared for \code{double} and -\code{std::complex} by the macro \code{declare_la}, but not all declared complex variants are +\paragraph{Element types and availability.} The functions are declared for \code{double}, \code{float}, +\code{std::complex} and \code{std::complex} by the macro \code{declare_la}, but not all declared complex variants are defined in \header{nonclass.cc}. Calling an undefined one gives a link error (table~\ref{tab:nonclass:avail}). \begin{table}[htbp] @@ -54,15 +54,15 @@ defined in \header{nonclass.cc}. Calling an undefined one gives a link error \toprule Function & Available for \\ \midrule -\texttt{linear\_solve(NRMat\&, NRMat*, ...)} & \texttt{double}, \texttt{complex} \\ -\texttt{linear\_solve(NRMat\&, NRVec\&, ...)}, \texttt{linear\_solve(NRSMat\&, ...)} & \texttt{double}, \texttt{complex} (complex \texttt{NRSMat}: hermitian) \\ +\texttt{linear\_solve(NRMat\&, NRMat*, ...)} & \texttt{double}, \texttt{float}, \texttt{complex}, \texttt{complex} \\ +\texttt{linear\_solve(NRMat\&, NRVec\&, ...)}, \texttt{linear\_solve(NRSMat\&, ...)} & \texttt{double}, \texttt{float}, \texttt{complex}, \texttt{complex} (complex \texttt{NRSMat}: hermitian) \\ \texttt{determinant}, \texttt{determinant\_destroy} & \texttt{NRMat}, \texttt{NRSMat}, real and complex \\ \texttt{diagonalize}, \texttt{singular\_decomposition}, \texttt{cholesky}, \texttt{diagofproduct}, -\texttt{trace2(NRMat,NRMat)}, \texttt{linear\_solve\_x}, \texttt{multiply\_by\_inverse}, \texttt{MatrixNorm}, \texttt{CondNumber} & \texttt{double}, \texttt{complex} \\ +\texttt{trace2(NRMat,NRMat)}, \texttt{linear\_solve\_x}, \texttt{multiply\_by\_inverse}, \texttt{MatrixNorm}, \texttt{CondNumber} & \texttt{double}, \texttt{float}, \texttt{complex}, \texttt{complex} \\ \texttt{gdiagonalize} & real matrix (real or complex eigenvector storage), complex matrix \\ -\texttt{trace2} with \texttt{NRSMat} & \texttt{double}, \texttt{complex} (hermitian) \\ +\texttt{trace2} with \texttt{NRSMat} & \texttt{double}, \texttt{float}, \texttt{complex}, \texttt{complex} (hermitian) \\ \texttt{diagonalize3}, \texttt{realmatrixfunction}, \texttt{complexmatrixfunction}, -\texttt{adjustphases}, \texttt{xcopy}, \texttt{xaxpy}, \texttt{xdot} & \texttt{double} only \\ +\texttt{adjustphases}, \texttt{xcopy}, \texttt{xaxpy}, \texttt{xdot} & \texttt{double}, \texttt{float} \\ templates of \texttt{matexp.h}, \texttt{power}, \texttt{commutator}, \texttt{twoside\_transform}, \dots & any type providing the needed operations \\ \bottomrule \end{tabular} @@ -106,17 +106,17 @@ $n_{\mathrm{rhs}}=0$).} \lafunc{void linear_solve(NRSMat &a, NRMat *b, T *det=0, int n=0); void linear_solve(NRSMat &a, NRVec &b, T *det=0, int n=0)}{Symmetric (real) or hermitian (complex), possibly indefinite $A$ in packed storage, \code{dspsv}/\code{zhpsv} (Bunch--Kaufman factorization; the determinant is computed from its $1\times1$ and $2\times2$ blocks and is real in the hermitian case). \code{n} selects a leading block, which is contiguous in the packed storage.} -\lafunc{int linear_solve_x(NRMat &A, T *B, const int rhsCount, const int eqCount, const bool eq, const bool saveA, double *rcond)}{Expert solver +\lafunc{int linear_solve_x(NRMat &A, T *B, const int rhsCount, const int eqCount, const bool eq, const bool saveA, typename LA_traits::normtype *rcond)}{Expert solver (\code{dgesvx}/\code{zgesvx}): solves $A x=b$ for \code{rhsCount} right-hand sides stored consecutively (each of length \code{eqCount}) in the raw array \code{B}, using the leading \code{eqCount}${}\times{}$\code{eqCount} block of \code{A}. \code{eq}: equilibrate $A$ (for badly conditioned matrices); then \code{A} is overwritten unless \code{saveA}. The reciprocal condition number is stored in \code{*rcond} if non-NULL. Returns LAPACK \code{info} (0 success, $>0$ singular or ill-conditioned). The solution replaces \code{B}.} -\lafunc{int multiply_by_inverse(NRMat &A, NRMat &B, bool useEq=false, double *rcond=NULL)}{$B\leftarrow A B^{-1}$ for square $A$, $B$ via +\lafunc{int multiply_by_inverse(NRMat &A, NRMat &B, bool useEq=false, typename LA_traits::normtype *rcond=NULL)}{$B\leftarrow A B^{-1}$ for square $A$, $B$ via \code{dgesvx}/\code{zgesvx}; \code{useEq} as \code{eq} above. \code{A} is not modified (with \code{useEq} a temporary copy is scaled). Returns \code{info}.} \lafunc{const NRMat calcinverse(NRMat a, T *det=NULL)}{$A^{-1}$ by \code{linear_solve} with unit right-hand sides; optionally the determinant. \code{NRMat::inverse()} (chapter~\ref{ch:mat}) calls it.} -\lafunc{const NRMat calcsvdinverse(NRMat a, double thr=0, bool do_absadjoint=false)}{Inverse via SVD, $V\,\mathrm{diag}(1/s_i)\,U\Hc$, +\lafunc{const NRMat calcsvdinverse(NRMat a, typename LA_traits::normtype thr=0, bool do_absadjoint=false)}{Inverse via SVD, $V\,\mathrm{diag}(1/s_i)\,U\Hc$, singular values $s_i<$\code{thr} are treated as zero (pseudoinverse). With \code{do_absadjoint} returns $V\,\mathrm{diag}(\prod_{j\ne i}s_j)\,U\Hc$, the adjugate matrix up to the phase $\det U\det V$. Square matrices only.} \lafunc{NRVec svd_solve(NRMat &a, const NRVec &b, normtype thres=LA_precision::svd_threshold(), int *nullity=NULL); NRMat svd_solve(NRMat &a, const NRMat &b, normtype thres=LA_precision::svd_threshold(), int *nullity=NULL)}{Least-squares @@ -126,7 +126,7 @@ the number of discarded ones is stored in \code{*nullity}. Here (unlike \code{li \lafunc{const LA_traits::elementtype determinant(MAT a)}{Determinant by LU (\code{NRMat}) or Bunch--Kaufman (\code{NRSMat}) factorization of a copy.} \lafunc{const LA_traits::elementtype determinant_destroy(MAT &a)}{The same, overwriting \code{a}. Unless \code{HAS_MKL} is defined, a dummy right-hand side is passed to \code{linear_solve} (work-around for OpenBLAS, which skips the factorization without right-hand sides).} -\lafunc{void cholesky(NRMat &a, bool upper=1); void cholesky(NRMat > &a, bool upper=1)}{Cholesky factorization of a +\lafunc{void cholesky(NRMat &a, bool upper=1)}{Cholesky factorization of a positive definite matrix in place (\code{dpotrf}/\code{zpotrf}). \code{upper=1}: \code{a} becomes the upper triangular $U$ with $A=U\Hc U$; \code{upper=0}: lower triangular $L$ with $A=L L\Hc$. The other triangle is zeroed. Raises \code{laerror} if $A$ is not positive definite.} \end{lareference} @@ -138,15 +138,15 @@ Eigenvalues are always returned in \textbf{ascending order} (LAPACK). Only the \ (row index ${}\ge{}$ column index) of an \code{NRMat} is referenced (verified); \code{NRSMat} stores only that triangle and is interpreted as hermitian (sec.~\ref{sec:smat:hermitian}). \begin{lareference} -\lafunc{void diagonalize(NRMat &a, NRVec &w, const bool eivec=1, const bool corder=1, int n=0, NRMat *b=NULL, const int itype=1)}{Eigenvalues +\lafunc{void diagonalize(NRMat &a, NRVec::normtype> &w, const bool eivec=1, const bool corder=1, int n=0, NRMat *b=NULL, const int itype=1)}{Eigenvalues \code{w} (size must equal \code{a.nrows()}) and, if \code{eivec}, eigenvectors stored in \code{a} (columns for \code{corder=1}, see table~\ref{tab:nonclass:eivec}) of a real symmetric or complex hermitian matrix (\code{dsyev}/\code{zheev}). If \code{b} is given, the generalized problem with positive definite $B$ is solved (\code{dsygv}/\code{zhegv}): \code{itype=1}: $A x=\lambda B x$, \code{itype=2}: $A B x=\lambda x$, \code{itype=3}: $B A x=\lambda x$; eigenvectors are $B$-normalized ($x\Hc B x=1$ for \code{itype} 1, 2), \code{b} is overwritten by its Cholesky factor.} -\lafunc{void diagonalize(NRSMat &a, NRVec &w, NRMat *v, const bool corder=1, int n=0, NRSMat *b=NULL, const int itype=1)}{The same for packed +\lafunc{void diagonalize(NRSMat &a, NRVec::normtype> &w, NRMat *v, const bool corder=1, int n=0, NRSMat *b=NULL, const int itype=1)}{The same for packed storage (\code{dspev}/\code{zhpev}, \code{dspgv}/\code{zhpgv}); eigenvectors go to \code{*v} (square, at least $n\times n$) or are not computed if \code{v==NULL}. \code{a} (and \code{b}) are destroyed.} -\lafunc{void diagonalize3(NRVec &d, NRVec &d1, NRMat *v, const bool corder=1, int n0=0)}{Real symmetric tridiagonal matrix +\lafunc{void diagonalize3(NRVec &d, NRVec &d1, NRMat *v, const bool corder=1, int n0=0)}{(also for \code{float}) Real symmetric tridiagonal matrix (\code{dsteqr}) with diagonal \code{d} and sub-diagonal \code{d1} (size ${}\ge n-1$). On return \code{d} holds the eigenvalues (ascending), \code{d1} is destroyed, \code{*v} (exactly $n\times n$) the eigenvectors. \code{n0>0} uses only the first \code{n0} elements. Used by \code{lanczos} (chapter~\ref{ch:iterative}).} \end{lareference} @@ -158,7 +158,7 @@ may contain zeros for infinite eigenvalues). Right (\code{vr}) and left (\code{v pointers are non-NULL; storage according to table~\ref{tab:nonclass:eivec}. Without sorting, the eigenvalues come in LAPACK order (unsorted, complex conjugate pairs adjacent with positive imaginary part first). \begin{lareference} -\lafunc{void gdiagonalize(NRMat &a, NRVec &wr, NRVec &wi, NRMat *vl=NULL, NRMat *vr=NULL, const bool corder=1, int n=0, const int sorttype=0, const int biorthonormalize=0, NRMat *b=NULL, NRVec *beta=NULL)}{Real +\lafunc{void gdiagonalize(NRMat &a, NRVec &wr, NRVec &wi, NRMat *vl=NULL, NRMat *vr=NULL, const bool corder=1, int n=0, const int sorttype=0, const int biorthonormalize=0, NRMat *b=NULL, NRVec *beta=NULL)}{(also for \code{float}) Real matrix, eigenvalues $w_r+\mathrm{i}\,w_i$, eigenvectors in the real storage of \code{dgeev}: for a real eigenvalue the vector is real; for a complex pair $(j,j{+}1)$ the vectors $j$ and $j{+}1$ hold the real and imaginary part. Verified: the right eigenvector belonging to $w_{r,j}+\mathrm{i}w_{i,j}$ ($w_{i,j}>0$) is $\mathrm{Re}-\mathrm{i}\,\mathrm{Im}$, and $\mathrm{Re}+\mathrm{i}\,\mathrm{Im}$ of \code{vl} satisfies $y\Tr A=\lambda y\Tr$ -- @@ -167,15 +167,15 @@ note that this sign is opposite to the \code{dgeev} documentation. 2 real eigenvalues first, then complex ones, each group ordered as for 1. Eigenvectors are permuted accordingly. \code{biorthonormalize} (needs \code{vl} and \code{vr}, not with \code{b}): 1 scales the left, 2 the right eigenvectors so that $y_k\Tr x_k=1$; complex pairs are rotated to biorthonormality. \code{a} is destroyed (and used as scratch space).} -\lafunc{void gdiagonalize(NRMat &a, NRVec > &w, NRMat > *vl=NULL, NRMat > *vr=NULL, const bool corder=1, int n=0, const int sorttype=0, const int biorthonormalize=0, NRMat *b=NULL, NRVec *beta=NULL)}{Template -specialized for \code{T=double} (calls the previous routine and converts to complex eigenvalues and eigenvectors) and \code{T=complex} (\code{zgeev}/\code{zggev}). +\lafunc{void gdiagonalize(NRMat &a, NRVec::complextype> &w, NRMat::complextype> *vl=NULL, NRMat::complextype> *vr=NULL, const bool corder=1, int n=0, const int sorttype=0, const int biorthonormalize=0, NRMat *b=NULL, NRVec *beta=NULL)}{Template +specialized for \code{T=double}, \code{float} (calls the previous routine and converts to complex eigenvalues and eigenvectors) and \code{T=complex}, \code{complex} (\code{zgeev}/\code{zggev}, \code{cgeev}/\code{cggev}). For the complex matrix, \code{sorttype>0} raises \code{laerror} (not implemented) and \code{biorthonormalize} scales so that $\bar y_k\Tr x_k=1$ in the row storage.} -\lafunc{void gdiagonalize(NRMat > &a, NRVec &wr, NRVec &wi, NRMat > *vl=NULL, NRMat > *vr=NULL, const bool corder=1, int n=0, const int sorttype=0, const int biorthonormalize=0, NRMat > *b=NULL, NRVec > *beta=NULL)}{Complex +\lafunc{void gdiagonalize(NRMat > &a, NRVec &wr, NRVec &wi, NRMat > *vl=NULL, NRMat > *vr=NULL, const bool corder=1, int n=0, const int sorttype=0, const int biorthonormalize=0, NRMat > *b=NULL, NRVec > *beta=NULL)}{(also for \code{complex}) Complex matrix with real and imaginary parts of the eigenvalues in separate vectors; exists to make \code{davidson} compile for complex \code{T}.} \end{lareference} \section{Singular value decomposition}\label{sec:nonclass:svd} \begin{lareference} -\lafunc{void singular_decomposition(NRMat &a, NRMat *u, NRVec &s, NRMat *v, const bool vnotdagger=0, int m=0, int n=0)}{SVD +\lafunc{void singular_decomposition(NRMat &a, NRMat *u, NRVec::normtype> &s, NRMat *v, const bool vnotdagger=0, int m=0, int n=0)}{SVD $A=U\,\mathrm{diag}(s)\,V\Hc$ of the $m\times n$ matrix \code{a} (\code{dgesvd}/\code{zgesvd}; \code{m}, \code{n} select a leading block, default whole matrix). \code{s} (size ${}\ge\min(m,n)$) receives the singular values in descending order. \code{u}: \code{NULL} (not computed), $m\times m$ (all columns) or $m\times\min(m,n)$ (economy size, chosen by \code{u} being non-square). \code{v}: \code{NULL}, $n\times n$ or $\min(m,n)\times n$, receives $V\Hc$ @@ -186,11 +186,12 @@ $m\times\min(m,n)$ (economy size, chosen by \code{u} being non-square). \code{v} Functions of matrices computed via the eigendecomposition $f(A)=V f(\Lambda) V^{-1}$. For general matrices the eigenvalue problem is solved in complex arithmetic and the result converted back to the element type of the argument; for a real argument the imaginary part is \emph{silently dropped} (its size is not checked). The matrix must be -diagonalizable; nothing is done about defective or ill-conditioned eigenvector matrices. +diagonalizable; nothing is done about defective or ill-conditioned eigenvector matrices. The functions listed below with +\code{double} in the signature exist also with \code{float} (and \code{complex} in place of \code{complex}). \begin{lareference} -\lafunc{NRMat matrixfunction(NRSMat a, double (*f)(double))}{$f(A)$ for a real symmetric or complex hermitian packed matrix, real function \code{f} +\lafunc{NRMat matrixfunction(NRSMat a, typename LA_traits::normtype (*f)(typename LA_traits::normtype))}{$f(A)$ for a real symmetric or complex hermitian packed matrix, real function \code{f} of the (real) eigenvalues.} -\lafunc{NRMat matrixfunction(NRMat a, std::complex (*f)(const std::complex &))}{$f(A)$ for a general real or complex +\lafunc{NRMat matrixfunction(NRMat a, typename LA_traits::complextype (*f)(const typename LA_traits::complextype &))}{$f(A)$ for a general real or complex matrix by \code{gdiagonalize}, $f(A)=\sum_k f(\lambda_k)\,x_k y_k\Tr/(y_k\Tr x_k)$.} \lafunc{NRMat realmatrixfunction(NRMat a, double (*f)(double))}{$f(A)$ for a real symmetric matrix in full storage (only the lower triangle is used).} \lafunc{NRMat > complexmatrixfunction(NRMat a, double (*fre)(double), double (*fim)(double))}{Complex function @@ -285,7 +286,7 @@ with two \code{gemm} calls; \code{trx}, \code{tryy} transpose the factors (no co \section{Norms and condition number}\label{sec:nonclass:norms} \begin{lareference} \lafunc{typename LA_traits::normtype MatrixNorm(const MAT &A, const char norm)}{Matrix norm by \code{dlange}/\code{zlange}, implemented for -\code{NRMat} and \code{NRMat >} (explicit specializations of the template, which can be added for other types). +\code{NRMat} of \code{double}, \code{float}, \code{std::complex} and \code{std::complex} (explicit specializations of the template, which can be added for other types). \code{'1'} or \code{'O'}: 1-norm (maximum column sum); \code{'I'}: $\infty$-norm (maximum row sum); \code{'F'} or \code{'E'}: Frobenius norm; \code{'M'}: largest absolute value of an element (lower case accepted). Works for non-square matrices.} \lafunc{typename LA_traits::normtype CondNumber(const MAT &A, const char norm)}{Condition number $\|A\|\,\|A^{-1}\|$ of a square matrix in the 1-norm @@ -296,10 +297,10 @@ infinity for a singular matrix. It is LAPACK's estimate, a lower bound that is u \section{Real and imaginary parts}\label{sec:nonclass:complex} \begin{lareference} \lafunc{const LA_traits::realtype realpart(const T &a); const LA_traits::realtype imagpart(const T &a)}{Real or imaginary part as a real object; -specialized for \code{NRVec}, \code{NRMat}, \code{NRSMat} and \code{SparseMat} over \code{complex}.} +specialized for \code{NRVec}, \code{NRMat}, \code{NRSMat} and \code{SparseMat} over \code{complex} and \code{complex}.} \lafunc{const LA_traits::complextype realmatrix(const T &a); const LA_traits::complextype imagmatrix(const T &a)}{Complex object $A$ resp.\ $\mathrm{i}A$ from a real one; -for \code{NRMat} and \code{SparseSMat}.} -\lafunc{const LA_traits::complextype complexmatrix(const T &re, const T &im)}{Complex object $\mathrm{re}+\mathrm{i}\,\mathrm{im}$; for \code{NRMat} and \code{SparseSMat}.} +for \code{NRMat} and \code{SparseSMat} of \code{double} and \code{float}.} +\lafunc{const LA_traits::complextype complexmatrix(const T &re, const T &im)}{Complex object $\mathrm{re}+\mathrm{i}\,\mathrm{im}$; for \code{NRMat} and \code{SparseSMat} of \code{double} and \code{float}.} \end{lareference} For the same purpose the classes provide \code{complexify()} (chapters~\ref{ch:vec}, \ref{ch:mat}). The real and imaginary part functions are CPU only (sec.~\ref{sec:concepts:gpu}). @@ -314,7 +315,7 @@ Application from the source: overlap of Slater determinants with signed spin-orb \lafunc{NRMat reconstructmatrix(const M &implicitmat)}{Explicit matrix of an implicit one (applying \code{gemv} to unit vectors), for debugging; both template arguments must be given. Equivalent to \code{explicit_matrix} of \header{davidson.h} (chapter~\ref{ch:iterative}).} \lafunc{void xcopy(int n, const T *x, int incx, T *y, int incy); void xaxpy(int n, const T &a, const T *x, int incx, T *y, int incy); T xdot(int n, const T *x, int incx, const T *y, int incy)}{Type-generic -names of the CBLAS routines (to ease switching between \code{float}/\code{double}/complex in user code); at present specialized only for \code{double}.} +names of the CBLAS routines (to ease switching between \code{float}/\code{double}/complex in user code); specialized for \code{double} and \code{float}.} \end{lareference} \section{Examples}\label{sec:nonclass:examples} diff --git a/doc/polynomial.tex b/doc/polynomial.tex index e0c6fb2..ecef977 100644 --- a/doc/polynomial.tex +++ b/doc/polynomial.tex @@ -77,9 +77,9 @@ vector, hence also a \code{Polynomial} (section~\ref{sec:nonclass:matexp}). \subsection{Roots, interpolation, gcd} \begin{lareference} \lafunc{NRMat companion() const}{companion matrix ($d\times d$, ones on the superdiagonal, last row $-c_i/c_d$) whose characteristic polynomial is $p/c_d$.} -\lafunc{NRVec::complextype> roots() const}{all roots as eigenvalues of the companion matrix (\code{gdiagonalize}, chapter~\ref{ch:nonclass}); for \code{double} - and \code{std::complex} only (laerror for integer types).} -\lafunc{NRVec realroots(const typename LA_traits::normtype thr) const}{roots with $|\mathrm{Im}|<$ \code{thr}, real parts sorted ascending (\code{T = double}).} +\lafunc{NRVec::complextype> roots() const}{all roots as eigenvalues of the companion matrix (\code{gdiagonalize}, chapter~\ref{ch:nonclass}); for the floating point types + \code{double}, \code{float}, \code{std::complex} and \code{std::complex} only (laerror for integer types).} +\lafunc{NRVec realroots(const typename LA_traits::normtype thr) const}{roots with $|\mathrm{Im}|<$ \code{thr}, real parts sorted ascending (floating point types).} \lafunc{T newton(const T x0, const typename LA_traits::normtype thr=LA_precision::newton_threshold(), const int maxit=1000) const}{Newton iteration from \code{x0} until $|p(x)|<$ \code{thr}; no convergence report.} \lafunc{template Polynomial polyfromroots(const NRVec &roots, const int skip= -1)}{monic $\prod_{i\ne\mathrm{skip}}(x-r_i)$.} \lafunc{template Polynomial lagrange_interpolation(const NRVec &x, const NRVec &y)}{interpolating polynomial of degree $n-1$ through the $n$ points; laerror for repeated $x$.} diff --git a/doc/smat.tex b/doc/smat.tex index ebfeda5..f0b873c 100644 --- a/doc/smat.tex +++ b/doc/smat.tex @@ -14,11 +14,10 @@ matrices) are returned as \code{NRMat}. \code{NRSMat_from1} (section~\ref{sec:smat:from1}) is the same class with the two-index \code{operator()} counting from 1. -The file \code{smat.cc} contains explicit instantiations for the element types \code{double}, -\code{std::complex} and the integer types \code{char}, \code{short}, \code{int}, \code{long}, \code{long long} and -their unsigned variants (\code{float} is not instantiated). -Most numerical members have BLAS-based specializations only for the two types \code{double} and -\code{std::complex}. Some members are defined only for these two types: +The file \code{smat.cc} contains explicit instantiations for the element types \code{double}, \code{float}, +\code{std::complex}, \code{std::complex} and the integer types \code{char}, \code{short}, \code{int}, \code{long}, \code{long long} and +their unsigned variants. +Most numerical members have BLAS-based specializations only for the four floating point types. Some members are defined only for these types: \code{amax}, \code{amin}, \code{norm}, \code{dot}, \code{axpy}, \code{randomize}, \code{inverse} and the products with \code{NRMat} and \code{NRSMat}. @@ -105,7 +104,7 @@ are CPU only (\code{NOT_GPU}). \lafunc{NRSMat(const T &a, const int n)}{Order $n$, all $n(n+1)/2$ stored elements (i.e.\ also the off-diagonal ones) set to \code{a}.} \lafunc{NRSMat(const T *a, const int n)}{Copy of $n(n+1)/2$ packed elements (lower triangle by rows) from \code{a}.} \lafunc{NRSMat(const NRSMat &rhs)}{Shallow copy.} -\lafunc{NRSMat(const typename LA_traits_complex::NRSMat_Noncomplex_type &rhs, bool imagpart = false)}{Complex matrix from a real one: \code{rhs} becomes the real part (\code{imagpart=false}) or the imaginary part (\code{imagpart=true}), the other part is zero. Defined only for \code{T=std::complex}.} +\lafunc{NRSMat(const typename LA_traits_complex::NRSMat_Noncomplex_type &rhs, bool imagpart = false)}{Complex matrix from a real one: \code{rhs} becomes the real part (\code{imagpart=false}) or the imaginary part (\code{imagpart=true}), the other part is zero. Defined only for \code{T=std::complex} and \code{std::complex}.} \lafunc{explicit NRSMat(const NRMat &rhs)}{Symmetric part $(A+A\Tr)/2$ of a square dense matrix (deep copy; nonsquare input is checked only in DEBUG).} \lafunc{explicit NRSMat(const NRVec &rhs, const int n)}{Interprets a vector of length $n(n+1)/2$ as packed matrix of order $n$. \emph{Shares} the data with \code{rhs} (no copy); length checked only in DEBUG.} \lafunc{explicit NRSMat(const SparseMat &rhs); explicit NRSMat(const SparseSMat &rhs)}{Dense from sparse, see \ref{ch:sparse}.} @@ -161,20 +160,20 @@ generated from the in-place ones. \lafunc{void gemv(const T beta, NRVec > &r, const char trans, const T alpha, const NRVec > &x) const}{The same, real matrix and complex vectors.} \lafunc{void axpy(const T alpha, const NRSMat &x)}{$S \leftarrow S + \alpha x$. For hermitian complex matrices the result is hermitian only for real $\alpha$ (not checked).} -\lafunc{void rank2update(const T alpha, const NRVec &x, const NRVec &y)}{Rank-2 update $S \leftarrow S + \alpha x y\Hc + \bar\alpha y x\Hc$ (real: $\alpha(xy\Tr+yx\Tr)$); \code{dspr2}/\code{zhpr2} for \code{double}/\code{complex}, generic loop otherwise. A symmetric rank-1 update $S \leftarrow S + \beta x x\Hc$ (real $\beta$) is obtained as \code{rank2update(beta/2,x,x)}; there is no separate rank-1 member.} +\lafunc{void rank2update(const T alpha, const NRVec &x, const NRVec &y)}{Rank-2 update $S \leftarrow S + \alpha x y\Hc + \bar\alpha y x\Hc$ (real: $\alpha(xy\Tr+yx\Tr)$); \code{dspr2}/\code{zhpr2}/\code{sspr2}/\code{chpr2} for \code{double}/\code{complex}/\code{float}/\code{complex}, generic loop otherwise. A symmetric rank-1 update $S \leftarrow S + \beta x x\Hc$ (real $\beta$) is obtained as \code{rank2update(beta/2,x,x)}; there is no separate rank-1 member.} \lafunc{const T dot(const NRSMat &rhs) const}{$\sum_{i\ge j} \bar A_{ij} B_{ij}$, i.e.\ the dot product of the \emph{packed} vectors (off-diagonal elements counted once; \code{ddot}/\code{zdotc}).} \lafunc{const T dot(const NRVec &rhs) const}{Dot product of the packed data with a vector of length $n(n+1)/2$.} \lafunc{const typename LA_traits::normtype norm(const T scalar = (T)0) const}{Euclidean norm of the packed data of $S - \mathtt{scalar}\cdot 1$; off-diagonal elements are counted \emph{once}, so this is not the Frobenius norm of the full matrix.} \lafunc{const T trace() const}{Trace.} \lafunc{const T sum(bool onlytriangle=false) const}{Sum of all elements of the full matrix, or of the stored triangle only.} \lafunc{const NRVec rsum() const; const NRVec csum() const}{Row sums and column sums of the full matrix (equal for real, complex conjugate of each other for hermitian matrices).} -\lafunc{inline const T amax() const; inline const T amin() const}{Stored element with the largest/smallest absolute value (complex: $|\Re|+|\Im|$); \code{double} and \code{complex} only.} +\lafunc{inline const T amax() const; inline const T amin() const}{Stored element with the largest/smallest absolute value (complex: $|\Re|+|\Im|$); floating point types only.} \lafunc{const T* diagonalof(NRVec &r, const bool divide = 0, bool cache = false) const}{\code{divide=false}: store the diagonal into \code{r} and return \code{&r[0]}; \code{divide=true}: divide \code{r} elementwise by the nonzero diagonal elements and return \code{NULL}. \code{cache} is unused. Interface used by iterative solvers (\ref{ch:iterative}).} -\lafunc{NRSMat inverse()}{Inverse via \code{NRMat(*this).inverse()}, symmetrized back. \code{double} and \code{complex} only.} +\lafunc{NRSMat inverse()}{Inverse via \code{NRMat(*this).inverse()}, symmetrized back. Floating point types only.} \lafunc{NRSMat& conjugateme(); const NRSMat conjugate() const}{Complex conjugation of all stored elements (no-op for real types).} \lafunc{NRSMat::normtype> abs() const}{Elementwise absolute values.} \lafunc{NRSMat::normtype> diffabs(const NRSMat &rhs) const}{Elementwise $|A_{ij}|-|B_{ij}|$.} -\lafunc{void randomize(const typename LA_traits::normtype &x)}{Uniformly distributed elements in $(-x,x)$; complex: real and imaginary parts independent, imaginary part of the diagonal set to zero. \code{double} and \code{complex} only.} +\lafunc{void randomize(const typename LA_traits::normtype &x)}{Uniformly distributed elements in $(-x,x)$; complex: real and imaginary parts independent, imaginary part of the diagonal set to zero. Floating point types only.} \lafunc{const bool operator==(const NRSMat &rhs) const; const bool operator!=(const NRSMat &rhs) const}{Exact elementwise comparison.} \end{lareference} diff --git a/doc/sparse.tex b/doc/sparse.tex index 63ece73..e184ecf 100644 --- a/doc/sparse.tex +++ b/doc/sparse.tex @@ -26,8 +26,8 @@ Common definitions in \code{sparsemat.h}: \lafunc{const double SPARSEEPSILON=1e-15}{Kept for user code. The library itself uses the precision-dependent threshold \code{LA_precision::sparse_epsilon()} (section~\ref{sec:concepts:traits}; $10^{-15}$ for double, i.e.\ $\epsilon_{\rm sp}$, and $10^{-7}$ for float): elements with absolute value below (or equal to) it are dropped by \code{simplify()} and by most constructors and accumulating operations; it is denoted $\epsilon_{\rm sp}$ below.} \lafunc{template struct matel {T elem; SPMatindex row; SPMatindex col; matel *next;}}{List element of \code{SparseMat}; also the value type returned by the iterators of both classes.} \end{lareference} -Both classes are explicitly instantiated for \code{double} and \code{std::complex}. -Several operations are implemented for \code{double} only. \code{SparseMat} nowhere conjugates when +Both classes are explicitly instantiated for \code{double}, \code{float}, \code{std::complex} and \code{std::complex}. +Several operations are implemented for the real types \code{double} and \code{float} only. \code{SparseMat} nowhere conjugates when it mirrors elements of a symmetric matrix (source comments: ``not OK for complex hermitean''), so its complex symmetric matrices are complex symmetric, not hermitian. \code{SparseSMat} stores both triangles and its mirroring operations follow the hermitian convention of \code{NRSMat} (section~\ref{sec:smat:hermitian}). @@ -132,7 +132,7 @@ Binary operators with a matrix operand must not combine a symmetric left operand \lafunc{const T trace() const}{Trace.} \lafunc{const typename LA_traits::normtype norm(const T scalar=(T)0) const}{Frobenius norm of $A-\mathtt{scalar}\cdot1$ of the full matrix (symmetric off-diagonal elements counted twice). Calls \code{simplify()}.} \lafunc{const T dot(const SparseMat &rhs) const}{$\sum_{ij} a_{ij} b_{ij}$ without conjugation (also for complex types), general (not symmetric) matrices only; merge of the two lists sorted by rows.} -\lafunc{const T dot(const NRMat &rhs) const}{$\sum_{ij} a_{ij} b_{ij}$ with a dense matrix; \code{double} only.} +\lafunc{const T dot(const NRMat &rhs) const}{$\sum_{ij} a_{ij} b_{ij}$ with a dense matrix; \code{double} and \code{float} only.} \lafunc{const T* diagonalof(NRVec &r, const bool divide=0, bool cache=false) const}{Square matrix: diagonal into \code{r} (\code{divide=false}, returns \code{&r[0]}) or divide \code{r} by the nonzero diagonal elements (\code{divide=true}, returns \code{NULL}). Rectangular matrix: the diagonal of $A\Tr A$ (assuming a simplified matrix), for preconditioning. \code{cache} is unused.} \end{lareference} @@ -140,7 +140,7 @@ Binary operators with a matrix operand must not combine a symmetric left operand \begin{lareference} \lafunc{inline const NRVec operator*(const NRVec &rhs) const}{$A x$.} \lafunc{void gemv(const T beta, NRVec &r, const char trans, const T alpha, const NRVec &x, bool treat_as_symmetric=false) const}{$r \leftarrow \beta r + \alpha\,\mathrm{op}(A)x$; forwards to \code{NRVec::gemv}. Any \code{trans} other than \code{'n'} means transposition (no conjugation). With \code{treat_as_symmetric} a general matrix is used as if it had the symmetric flag. The sizes of \code{x} and (in DEBUG) of the result are checked.} -\lafunc{inline const NRMat operator*(const NRMat &rhs) const}{$A B$ via \code{NRMat::gemm} with a sparse first operand; \code{double} only (\code{complex} raises ``not implemented yet'').} +\lafunc{inline const NRMat operator*(const NRMat &rhs) const}{$A B$ via \code{NRMat::gemm} with a sparse first operand; \code{double} and \code{float} only (the complex types raise ``not implemented yet'').} \lafunc{const SparseMat operator*(const SparseMat &rhs) const}{Sparse product of two general matrices (error if one is symmetric); the result is simplified.} \lafunc{void gemm(const T beta, const SparseMat &a, const char transa, const SparseMat &b, const char transb, const T alpha)}{$C \leftarrow \beta C + \alpha\,\mathrm{op}(A)\,\mathrm{op}(B)$ for general \code{a}, \code{b}; if \code{*this} is symmetric (possible only with $\beta\neq0$) only elements with $i\le j$ are added. Result simplified. \code{a} and \code{b} must not be \code{*this} (checked in DEBUG).} \end{lareference} @@ -148,7 +148,7 @@ Related functions in other classes: \code{NRVec::gemv(const T beta, const SparseMat &a, const char trans, const T alpha, const NRVec &x, const bool treat_as_symmetric=false)}, \code{NRVec::operator*(const SparseMat&)} ($x\Tr A$), \code{NRMat::gemm(const T &beta, const SparseMat &a, const char transa, const NRMat &b, const char transb, const T &alpha)} -(\code{double} only) and \code{NRMat::operator+=}, \code{operator-=(const SparseMat&)}. +(\code{double} and \code{float} only) and \code{NRMat::operator+=}, \code{operator-=(const SparseMat&)}. \subsection{Transposition, permutations, direct sums} \begin{lareference} @@ -180,7 +180,7 @@ e.g.\ \code{NRVec>} can be written and read; the ORCA format is def \lafunc{inline const SparseMat commutator(const SparseMat &x, const SparseMat &y, const bool trx=0, const bool tryy=0)}{$[\mathrm{op}(x),\mathrm{op}(y)] = xy-yx$, via two \code{gemm} calls (general matrices only).} \lafunc{inline const SparseMat anticommutator(const SparseMat &x, const SparseMat &y, const bool trx=0, const bool tryy=0)}{$xy+yx$.} \end{lareference} -In \ref{ch:nonclass}: \code{realpart()} and \code{imagpart()} of a \code{SparseMat>}. The generic +In \ref{ch:nonclass}: \code{realpart()} and \code{imagpart()} of a \code{SparseMat>} or \code{SparseMat>}. The generic matrix functions in \code{matexp.h} (e.g.\ \code{exp()}) work with \code{SparseMat}, and the iterative solvers of \ref{ch:iterative} use its \code{gemv}, \code{diagonalof}, \code{issymmetric}, \code{nrows} interface. @@ -306,7 +306,7 @@ Elements with $|S_{ij}|\le$ $\epsilon_{\rm sp}$ are omitted. Both vectors must h \lafunc{SparseSMat NRMat::operator*(const SparseSMat &rhs) const}{Declared in mat.h: dense times sparse, general product $AS$ with a sparse result (row $i$ accumulated as $\sum_k A_{ik}\,(\text{row }k\text{ of }S)$, zero $A_{ik}$ skipped; simplified). Instantiated for \code{double} and \code{std::complex}.} \end{lareference} In \ref{ch:nonclass}: \code{complexmatrix()}, \code{realmatrix()} and \code{imagmatrix()} build a -\code{SparseSMat>} from real ones. The generic \code{exp()} of \code{matexp.h} works with +\code{SparseSMat>} (\code{SparseSMat>}) from real ones. The generic \code{exp()} of \code{matexp.h} works with \code{SparseSMat}. \section{CSRMat}\label{sec:sparse:csr}