diff --git a/CLAUDE.md b/CLAUDE.md index 21cd71a..761aa38 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -18,7 +18,7 @@ decisions and plans agreed with the author, so that they survive a fresh clone. Quick test program: `g++ -std=gnu++11 -I -O0 -DDEBUG -DNO_STRASSEN -DFORTRAN_ -DNONCLAPACK X.cc /.libs/libla.a -lcblas -lblas -llapack` (add `-fsanitize=address` for memory checks). - Regression test of the optimizers: `make check` builds `testoptimize`; `./testoptimize` must report - "123 passed, 0 failed" (exit code = number of failures). New optimizer tests go to testoptimize.cc; + "137 passed, 0 failed" (exit code = number of failures). New optimizer tests go to testoptimize.cc; test problems live in testoptimize.h. - `t` (t.cc) is a collection of if(0)/if(1) test blocks, most read n from stdin (`echo 10 | ./t`). - NONCBLAS check: copy the repo with `cp -a` (keeps timestamps, otherwise automake reruns and fails), @@ -90,10 +90,15 @@ decisions and plans agreed with the author, so that they survive a fresh clone. its NEA oscillator strength f=dE_min/dE is a placeholder for the true one. ## Planned (only when the author asks) -- float version: float and complex must ALWAYS use BLAS/LAPACK (s*/c* routines, noncblas wrappers) - like double/complex, never generic loops. INSTANTIZE(float) in optimize.cc currently needs float - versions of diagonalize, linear_solve, singular_decomposition. BLAS-based NRVec members are meant for - floating types only, not integers. +- float version (IN PROGRESS since 2026-11-05): float and complex must ALWAYS use BLAS/LAPACK (s*/c* routines, + noncblas wrappers) like double/complex, never generic loops. Done: noncblas, vec, mat, smat, nonclass (+ complexify + BLAS for NRSMat/NRMat), sparsemat, sparsesmat, tensor, fourindex, polynomial, header templates (matexp, + conjgrad, gmres, davidson, lanczos, diis: no changes needed beyond the defaults). Float twins are generated from the double specializations (placed right after them). Precision- + dependent defaults live in LA_precision (la_traits.h, only double and float; double values unchanged). optimize done + 2026-11-05: own trait opt_precision in optimize.h, INSTANTIZE(float), float tests in testoptimize.cc (test_float); + float LM on badly scaled problems needs my_scaling=true. reg/regsurf, csrmat, strassen, CUDA stay double-only. + Open: Tensor Tucker bugs (complex roundtrip, inverseorder=true abort, also double; bugtests/tensor_tucker_bug.cc). + BLAS-based NRVec members are meant for floating types only, not integers. - NRSMat: possible bool member `hermitian` toggling hermitian vs complex symmetric behaviour (Tensor(const NRSMat&) must then respect it). - SparseSMat: possible member `symmetry_type` (0 none, 1 symmetric, 2 hermitian) for NRSMat conversion, diff --git a/doc/concepts.tex b/doc/concepts.tex index 6cb6ee7..d4f7df2 100644 --- a/doc/concepts.tex +++ b/doc/concepts.tex @@ -153,12 +153,36 @@ the mixed real--complex functions such as \code{gemv} of a complex vector with a The functions \code{LA::real(c)} and \code{LA::imag(c)} return references to the real and imaginary part of a \code{std::complex}, which can be used as l-values. +Default thresholds of the numerical methods depend on the floating point precision and are collected in +\code{LA_precision} (\header{la_traits.h}), which is defined for \code{T = double} and \code{float} only (any other type +gives a compile-time error); functions with an argument of type \code{T} use \code{LA_precision::normtype>}, +so that complex types take the values of their real component. The double values are those used by the library before. +\begin{center} +\begin{tabular}{llll} +\toprule +member & double & float & used by \\ +\midrule +\code{svd_threshold()} & $10^{-12}$ & $10^{-5}$ & \code{svd_solve}, \code{Tensor::Tucker} \\ +\code{taylor_precision()} & $10^{-14}$ & $10^{-7}$ & \code{exp}, \code{sincos}, \code{exptimes} (\header{matexp.h}) \\ +\code{conjgrad_tol()} & $10^{-8}$ & $10^{-5}$ & \code{conjgrad} \\ +\code{gmres_eps()} & $10^{-7}$ & $10^{-5}$ & \code{gmres} \\ +\code{eigensolver_eps()} & $10^{-6}$ & $10^{-4}$ & \code{davidson}, \code{lanczos} \\ +\code{newton_threshold()} & $10^{-14}$ & $10^{-6}$ & \code{Polynomial::newton} \\ +\code{det_tiny()} & $10^{-300}$ & $10^{-37}$ & underflow guard of the determinant in \code{linear_solve} \\ +\code{sparse_epsilon()} & $10^{-15}$ & $10^{-7}$ & zero elements of \code{SparseMat}, \code{SparseSMat} (\code{SPARSEEPSILON}) \\ +\code{cholesky_epsilon()} & $10^{-16}$ & $10^{-8}$ & pivots skipped by \code{SparseSMat::cholesky} (\code{CHOLESKYEPSILON}) \\ +\bottomrule +\end{tabular} +\end{center} + \subsection{Element types and instantiation} The class templates are explicitly instantiated in the library for \code{double}, \code{std::complex} and the integer types \code{char}, \code{short}, \code{int}, \code{long}, \code{long long} and their unsigned -variants. The BLAS and LAPACK backed functions are implemented for \code{double} and \code{std::complex} -(some also for mixed real--complex arguments); calling them for other types results in a run-time +variants; \code{NRVec}, \code{NRMat}, \code{NRSMat}, \code{SparseMat}, \code{SparseSMat}, \code{Tensor}, \code{fourindex} (with \code{unsigned char} indices) and \code{Polynomial} also for \code{float} and +\code{std::complex} (single precision support, 2026, in progress for the other classes). The BLAS and LAPACK backed +functions are implemented for \code{double} and \code{std::complex} (some also for mixed real--complex arguments) +and for their single precision counterparts, which always use the single precision BLAS/LAPACK routines (\code{s}/\code{c}); calling them for other types results in a run-time \code{laerror} (``not implemented for this type'') or a link error. Functions implemented in the headers work for any type with the necessary arithmetic operations. diff --git a/doc/fourindex.tex b/doc/fourindex.tex index 06e8b36..bf2c165 100644 --- a/doc/fourindex.tex +++ b/doc/fourindex.tex @@ -24,9 +24,9 @@ Indices are opaque to the sparse classes (they are never compared with the dimen 1-based indexing are possible there; the dense classes use indices from 1. Only \code{fourindex}, \code{fourindex_ext} and -\code{fourindex_dense} are explicitly instantiated in fourindex.cc; +\code{fourindex_dense} and the same with \code{float} are explicitly instantiated in fourindex.cc; everything else is implemented in the header and instantiated on use, except \code{fscanf()} and -\code{fprintf()} which exist only for \code{fourindex}. +\code{fprintf()} which exist only for \code{fourindex} and \code{fourindex}. \section{Symmetry types}\label{sec:fourindex:symmetry} @@ -131,7 +131,7 @@ The symmetry must be set (not \code{undefined_symmetry}) for \code{piterator} an \lafunc{unsigned long get(int fd, bool withattr=true)}{read the same format until end of file, \emph{adding} the elements to the list; returns their number.} \lafunc{template std::ostream& operator<<(std::ostream &s, const fourindex &x)}{text output: the dimension, lines \code{i j k l value}, and a final line with four terminators and 0.} \lafunc{template std::istream& operator>>(std::istream &s, fourindex &x)}{read the text format; \code{x} is resized to the read dimension (the symmetry is not set) and reading stops at the first line containing the terminator of \code{x} as an index.} -\lafunc{void fscanf(FILE *f); void fprintf(FILE *f, char *format) const}{C-style text I/O, only for \code{fourindex}; line format \code{value i j k l}. \code{fscanf} adds elements until a line in which any index equals the terminator (compared after conversion to the index type, so $-1$ in the file matches the default terminator 255); \code{fprintf} applies \code{format} (which must consume a \code{double} and four \code{int}s) to each element and ends with a line \code{0 t t t t} with the terminator $t$. With terminator 0 this is the FCIDUMP-like convention (blocks separated by zero lines, read by consecutive \code{fscanf} calls; a one-electron line $i\,j\,0\,0$ also ends a block).} +\lafunc{void fscanf(FILE *f); void fprintf(FILE *f, char *format) const}{C-style text I/O, only for \code{fourindex} and \code{}; line format \code{value i j k l}. \code{fscanf} adds elements until a line in which any index equals the terminator (compared after conversion to the index type, so $-1$ in the file matches the default terminator 255); \code{fprintf} applies \code{format} (which must consume a \code{double} and four \code{int}s) to each element and ends with a line \code{0 t t t t} with the terminator $t$. With terminator 0 this is the FCIDUMP-like convention (blocks separated by zero lines, read by consecutive \code{fscanf} calls; a one-electron line $i\,j\,0\,0$ also ends a block).} \end{lareference} \section{Disc storage: \texttt{fourindex\_ext}}\label{sec:fourindex:ext} diff --git a/doc/iterative.tex b/doc/iterative.tex index 2507a61..04fba89 100644 --- a/doc/iterative.tex +++ b/doc/iterative.tex @@ -14,6 +14,9 @@ or of the normal equations $A\Tr Ax=A\Tr b$ for rectangular $A$ (least squares); \item the class \code{DIIS} (\header{diis.h}): Pulay's DIIS convergence acceleration of a general fixed-point iteration; \item \code{bisection_find()} and \code{interpolation_find()} (\header{bisection.h}): search in a sorted array. \end{itemize} +All of them work for \code{T = float} and \code{std::complex} as well (where the double versions support the +corresponding real or complex type); their default thresholds then come from \code{LA_precision} +(section~\ref{sec:concepts:traits}). Dense (LAPACK based) diagonalization, linear solvers and matrix functions are in chapter~\ref{ch:nonclass}; the small subspace problems of \code{davidson} and \code{lanczos} are solved by the routines described there (\code{diagonalize}, \code{gdiagonalize}, \code{diagonalize3}, \code{linear_solve}). The quasi-Newton GDIIS optimizer built on @@ -96,7 +99,7 @@ problem. \laheader{davidson.h} \begin{lareference} -\lafunc{template void davidson(const Matrix &bigmat, NRVec &eivals, NRVec *eivecs, const char *eivecsfile, int nroots=1, const bool verbose=0, const double eps=1e-6, const bool incore=1, int maxit=100, const int maxkrylov = 500, void (*initguess)(NRVec &)=NULL, const typename LA_traits::normtype *target=NULL)}{ +\lafunc{template void davidson(const Matrix &bigmat, NRVec &eivals, NRVec *eivecs, const char *eivecsfile, int nroots=1, const bool verbose=0, const double eps=LA_precision::eigensolver_eps(), const bool incore=1, int maxit=100, const int maxkrylov = 500, void (*initguess)(NRVec &)=NULL, const typename LA_traits::normtype *target=NULL)}{ Davidson's method \cite{iterative:davidson1975} for the \code{nroots} lowest eigenvalues of \code{bigmat}. \begin{itemize} \item \code{eivals}: on return the eigenvalues in its first \code{nroots} elements (ascending, or by increasing distance from @@ -167,7 +170,7 @@ implicitly given one, built column by column by \code{gemv} on unit vectors (for \laheader{lanczos.h} \begin{lareference} -\lafunc{template void lanczos(const Matrix &bigmat, NRVec &eivals, NRVec *eivecs, const char *eivecsfile, int nroots=1, const bool verbose=0, const double eps=1e-6, const bool incore=1, int maxit=100, const int maxkrylov = 1000, void (*initguess)(NRVec &)=NULL, const typename LA_traits::normtype *target=NULL)}{ +\lafunc{template void lanczos(const Matrix &bigmat, NRVec &eivals, NRVec *eivecs, const char *eivecsfile, int nroots=1, const bool verbose=0, const double eps=LA_precision::eigensolver_eps(), const bool incore=1, int maxit=100, const int maxkrylov = 1000, void (*initguess)(NRVec &)=NULL, const typename LA_traits::normtype *target=NULL)}{ Lanczos method for the \code{nroots} lowest (or closest to \code{*target}) eigenvalues of a hermitian matrix. The parameters have the same meaning as for \code{davidson}, with these differences: \code{bigmat.issymmetric()} must be \code{true} (\code{laerror} otherwise); the start vector is normalized by the function; \code{incore=false} keeps only the Lanczos vectors @@ -194,7 +197,7 @@ the method fails with \code{laerror} when more than \code{maxit} iterations are \laheader{conjgrad.h} \begin{lareference} -\lafunc{template bool conjgrad(const Matrix &bigmat, const NRVec &b, NRVec &x, const bool doguess=true, const double tol=1e-8, const int itmax=1000, const bool verbose=true, bool issquare=1, const bool precondition=1)}{ +\lafunc{template bool conjgrad(const Matrix &bigmat, const NRVec &b, NRVec &x, const bool doguess=true, const double tol=LA_precision::conjgrad_tol(), const int itmax=1000, const bool verbose=true, bool issquare=1, const bool precondition=1)}{ Solves $Ax=b$ by the (Jacobi preconditioned) conjugate gradient method (the source refers to the book of Stoer and Bulirsch). Returns \code{true} if converged, \code{false} after \code{itmax}+1 iterations without convergence; \code{x} then holds the last iterate. @@ -227,7 +230,7 @@ $p\cdot Ap$ (or $|Ap|^2$) or the residual product $r\cdot z$ vanishes. Only real \laheader{gmres.h} \begin{lareference} -\lafunc{template bool gmres(const Matrix &bigmat, const NRVec &b, NRVec &x, const bool doguess=1, const double eps=1e-7, const int MAXIT=50, const bool verbose=1, bool square=1, const bool precondition=1, int neustart=0, const int incore=1)}{ +\lafunc{template bool gmres(const Matrix &bigmat, const NRVec &b, NRVec &x, const bool doguess=1, const double eps=LA_precision::gmres_eps(), const int MAXIT=50, const bool verbose=1, bool square=1, const bool precondition=1, int neustart=0, const int incore=1)}{ Solves $Ax=b$ by restarted GMRES with Givens rotations (after H.~R.~Schwarz, original implementation by M.~Warken). Returns \code{true} if converged; otherwise \code{false}, with \code{x} holding the last iterate. \begin{itemize} diff --git a/doc/nonclass.tex b/doc/nonclass.tex index af641b6..344d7cd 100644 --- a/doc/nonclass.tex +++ b/doc/nonclass.tex @@ -119,7 +119,7 @@ scaled). Returns \code{info}.} \lafunc{const NRMat calcsvdinverse(NRMat a, double 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, double thres=1e-12, int *nullity=NULL); NRMat svd_solve(NRMat &a, const NRMat &b, double thres=1e-12, int *nullity=NULL)}{Least-squares +\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 solution $x=V\,\mathrm{diag}(1/s_i)\,U\Hc b$ discarding singular values $\le$\code{thres} (no contamination by the null space); the number of discarded ones is stored in \code{*nullity}. Here (unlike \code{linear_solve}) the right-hand sides are the \emph{columns} of the matrix \code{b} ($m\times n_{\mathrm{rhs}}$); \code{a} ($m\times n$, $m\ge n$) is destroyed.} @@ -223,8 +223,8 @@ so they can be used with implicitly defined (``direct'') matrices as in chapter~ needs \code{size()}, \code{operator*=}, \code{axpy} and \code{resize()}. The exponential uses scaling and squaring: the power $p$ is chosen so that $\|A\|/2^p\lesssim 0.75$ -(\code{mynextpow2}), the Taylor series of the scaled argument is truncated when the next term would be below $10^{-14}$ -relative to the norm, and the result is squared $p$ times. The scaling is put into the Taylor coefficients, +(\code{mynextpow2}), the Taylor series of the scaled argument is truncated when the next term would be below +\code{LA_precision::taylor_precision()} ($10^{-14}$ for double, $10^{-7}$ for float) relative to the norm, and the result is squared $p$ times. The scaling is put into the Taylor coefficients, so the argument itself is not copied. \code{maxpower} (if ${}\ge0$) limits $p$, \code{maxtaylor} (if ${}\ge0$) limits the order of the Taylor polynomial, which is useful and exact for nilpotent matrices (e.g.\ coupled-cluster excitation operators). diff --git a/doc/optimize.tex b/doc/optimize.tex index 1cf4b3f..1d9a226 100644 --- a/doc/optimize.tex +++ b/doc/optimize.tex @@ -23,12 +23,13 @@ Methods are written on top of a small infrastructure (evaluation counter, iterat common drivers (line search, first order, trust region), so that a new method usually consists of a short class supplying the search direction, the update rule or the trust-region step (section~\ref{sec:optimize:drivers}). -All templates are parameterized by the real type \code{T}, but optimize.cc instantiates them for \code{T = double} only. -The capability interfaces, adapters and the hyperparameter classes are header-only templates; the optimizers, the adapters with -member functions in optimize.cc (\code{FrozenOptProblem}, \code{ParameterShift}, \code{FiniteDifference}, -\code{OptLeastSquares}) and the auxiliary functions exist in the library only for \code{double}. +All templates are parameterized by the real type \code{T}; optimize.cc instantiates them for \code{T = double} and +\code{float}. The capability interfaces, adapters and the hyperparameter classes are header-only templates; the optimizers, the +adapters with member functions in optimize.cc (\code{FrozenOptProblem}, \code{ParameterShift}, \code{FiniteDifference}, +\code{OptLeastSquares}) and the auxiliary functions exist in the library for \code{double} and \code{float}. The common drivers \code{linesearch_minimize}, \code{firstorder_minimize} and \code{trustregion_minimize}, the extension -points for user-written methods, are explicitly instantiated for \code{double} as well. +points for user-written methods, are explicitly instantiated for both types as well. The defaults depending on the precision +are taken from \code{opt_precision} (section~\ref{sec:optimize:hyper}). The data members of all classes of this chapter are small, except vectors held by \code{FrozenOptProblem} and \code{optevaluator}; these are ordinary \code{NRVec} objects following the reference counting and copy-on-demand scheme of @@ -254,6 +255,29 @@ hyperparameters maxstep(200) maxfunc(0) f_thr(0) x_thr(0) sad +- Cobyla_hyperparameters rhobeg(0.5) rhoend(1e-6) flexible_radius(true); maxstep=100000 \end{lstlisting} +The defaults given above and below are those for \code{double}. The precision-dependent ones are taken from the trait +\code{opt_precision} (\header{optimize.h}, analogous to \code{LA_precision} of section~\ref{sec:concepts:traits}), defined +for \code{T = double} and \code{float} only; the double values are those used by the library before. In float the gradient norm +reachable by a line search is limited by the resolution of $f$ (about $\epsilon|f|$ divided by the step), hence the larger +\code{g_thr}. +\begin{center} +\begin{tabular}{llll} +\toprule +member & double & float & default of \\ +\midrule +\code{fd_central_step()} & $10^{-6}$ & $5\cdot10^{-3}$ & \code{hh} of \code{FiniteDifference}, \code{ParameterShift} \\ +\code{fd_forward_step()} & $10^{-7}$ & $3\cdot10^{-4}$ & \code{h} of \code{OptLeastSquares}, \code{my_fd_step} (LM) \\ +\code{g_thr()} & $10^{-6}$ & $10^{-3}$ & \code{my_g_thr} \\ +\code{sy_thr()} & $10^{-12}$ & $10^{-6}$ & \code{my_sy_thr} (quasi-Newton) \\ +\code{sr1_thr()} & $10^{-8}$ & $10^{-4}$ & \code{my_sr1_thr} (SR1) \\ +\code{update_thr()} & $10^{-8}$ & $10^{-4}$ & \code{my_update_thr} (trust region) \\ +\code{nm_f_thr()}, \code{nm_x_thr()} & $10^{-8}$, $10^{-6}$ & $10^{-5}$, $10^{-4}$ & \code{my_f_thr}, \code{my_x_thr} (Nelder-Mead) \\ +\code{powell_f_thr()}, \code{powell_x_thr()} & $10^{-10}$, $10^{-6}$ & $10^{-6}$, $10^{-4}$ & \code{my_f_thr}, \code{my_x_thr} (Powell) \\ +\code{rhoend()} & $10^{-6}$ & $10^{-4}$ & \code{my_rhoend} (COBYLA) \\ +\bottomrule +\end{tabular} +\end{center} + \paragraph{Base classes.} \begin{lareference} \lafunc{template class hyperparameters}{Parameters common to all methods.} @@ -671,16 +695,19 @@ saddle order.} For least-squares problems providing \code{OptResiduals} (typically derived from \code{OptLeastSquares}): the step $h$ solves $(J\Tr J+\mu D)h=-J\Tr r$ (\code{linear_solve}), $D=I$ or $D=\mathrm{diag}(J\Tr J)$ (Marquardt scaling); the damping starts at -$\mu=\tau\max_i(J\Tr J)_{ii}$ and is updated by the gain ratio $\varrho$ \cite{optimize:madsen2004} (alg.~3.16): for $\varrho>0$ +$\mu=\tau\max_i(J\Tr J)_{ii}$ (with scaling $\mu=\tau$, the same relative to $D$) and is updated by the gain ratio $\varrho$ \cite{optimize:madsen2004} (alg.~3.16): for $\varrho>0$ the step is accepted and $\mu\leftarrow\mu\max(1/3,1-(2\varrho-1)^3)$, $\nu=2$, otherwise $\mu\leftarrow\mu\nu$, $\nu\leftarrow2\nu$. The Jacobian comes from \code{jacobian()} of the problem or from forward differences (relative step \code{my_fd_step}, $n$ residual evaluations counted in \code{my_maxfunc}). Convergence: $|J\Tr r|<$ \code{my_g_thr}, the general $f$/$x$ criterion, or -when a rejected step becomes negligible at machine precision ($|h|\le10\epsilon(|x|+\epsilon)$). The function value returned is +when a rejected step becomes negligible at machine precision ($h_i^2\le(10\epsilon)^2(x_i^2+\epsilon^2)$ for all $i$, +tested per component to be scale invariant, cf.\ \cite{optimize:madsen2004} eq.~3.16). The function value returned is $f=\frac12|r|^2$. The problem's \code{funcgrad()} is not used. \begin{lareference} \lafunc{template class LM_hyperparameters : public gradient_hyperparameters}{} \lafunc{T my_tau}{Initial damping factor $\tau$ (\code{1e-3}).} -\lafunc{bool my_scaling}{Marquardt scaling $D=\mathrm{diag}(J\Tr J)$ (default \code{false}, $D=I$).} +\lafunc{bool my_scaling}{Marquardt scaling $D=\mathrm{diag}(J\Tr J)$ (default \code{false}, $D=I$). Recommended for +badly scaled parameters, in float usually necessary (otherwise the steps in the weakly curved directions fall below the +resolution of $f$).} \lafunc{T my_fd_step}{Relative step of the forward-difference Jacobian (\code{1e-7}).} \lafunc{template class LM_Optimizer : public Optimizer}{\code{laerror} if the problem does not provide \code{OptResiduals}. Verbose name Levenberg-Marquardt.} diff --git a/doc/polynomial.tex b/doc/polynomial.tex index 25ec789..e0c6fb2 100644 --- a/doc/polynomial.tex +++ b/doc/polynomial.tex @@ -10,7 +10,7 @@ This chapter describes the polynomial class \code{Polynomial} and the exact r \code{template class Polynomial : public NRVec} represents $p(x)=\sum_{i=0}^{d}c_ix^i$; the coefficient $c_i$ is the vector element \code{p[i]} (indexed from 0, i.e.\ by the power) and the vector has $d+1$ elements. The member \code{size()} is made private to prevent confusion; use \code{degree()}. The leading coefficient may be zero (e.g.\ after a subtraction); \code{simplify()} removes such terms. The class is explicitly -instantiated for \code{int}, \code{unsigned int}, \code{double}, \code{std::complex}; root finding and SVD-based gcd require a floating point type. +instantiated for \code{int}, \code{unsigned int}, \code{double}, \code{std::complex}, \code{float}, \code{std::complex}; root finding and SVD-based gcd require a floating point type. All inherited \code{NRVec} functionality (norms, \code{randomize}, \code{clear}, I/O operators printing the number of coefficients followed by the coefficients, \ldots) is available. \code{Polynomial} follows the copy-on-demand scheme of section~\ref{sec:concepts:refcount}: its modifying members call @@ -80,7 +80,7 @@ vector, hence also a \code{Polynomial} (section~\ref{sec:nonclass:matexp}). \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{T newton(const T x0, const typename LA_traits::normtype thr=1e-14, const int maxit=1000) const}{Newton iteration from \code{x0} until $|p(x)|<$ \code{thr}; no convergence report.} +\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$.} \lafunc{template NRMat Sylvester(const Polynomial &p, const Polynomial &q)}{Sylvester matrix of dimension $\deg p+\deg q$ (its determinant is the resultant).} diff --git a/doc/sparse.tex b/doc/sparse.tex index e79208d..63ece73 100644 --- a/doc/sparse.tex +++ b/doc/sparse.tex @@ -23,7 +23,7 @@ Common definitions in \code{sparsemat.h}: \begin{lareference} \lafunc{typedef unsigned int SPMatindex; typedef int SPMatindexdiff}{Index type and index difference type. All row and column indices count from 0.} -\lafunc{const double SPARSEEPSILON=1e-15}{Threshold: elements with absolute value below (or equal to) it are dropped by \code{simplify()} and by most constructors and accumulating operations.} +\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}. @@ -43,7 +43,7 @@ which no longer exists, is not included by any file and is not installed; it is The nonzero elements form a singly linked list of \code{matel} (pointer \code{list}); new elements are prepended. The list is in no particular order and may contain several entries with the same $(i,j)$, which represent their sum; \code{simplify()} merges such duplicates, removes elements with -$|a_{ij}|\le$ \code{SPARSEEPSILON} and leaves the list sorted by rows. Operations needing an ordering +$|a_{ij}|\le$ $\epsilon_{\rm sp}$ and leaves the list sorted by rows. Operations needing an ordering (\code{simplify}, sparse products, \code{dot}) build private arrays of pointers sorted by rows or by columns (member \code{sort()}); these arrays and the cached element count \code{nonzero} are kept until the next operation which calls the private \code{unsort()}: \code{add()} (and everything built on it), \code{simplify}, @@ -86,7 +86,7 @@ Most modifying members call \code{copyonwrite()}: \code{simplify}, \code{operato \lafunc{SparseMat()}{Undefined empty matrix ($0\times0$, no reference count).} \lafunc{SparseMat(const SPMatindex n, const SPMatindex m)}{Defined empty (all zero) $n\times m$ matrix, not symmetric.} \lafunc{SparseMat(const SparseMat &rhs)}{Shallow copy; an empty matrix gets its own reference count.} -\lafunc{explicit SparseMat(const NRMat &rhs)}{From a dense matrix, elements with $|a_{ij}|>$ \code{SPARSEEPSILON}.} +\lafunc{explicit SparseMat(const NRMat &rhs)}{From a dense matrix, elements with $|a_{ij}|>$ $\epsilon_{\rm sp}$.} \lafunc{explicit SparseMat(const NRSMat &rhs)}{From a packed symmetric matrix; symmetric flag set, lower triangle stored.} \lafunc{~SparseMat()}{Destructor.} \lafunc{SparseMat & operator=(const SparseMat &rhs)}{Shallow copy (also copies the symmetric flag).} @@ -107,8 +107,8 @@ duplicates summed), \code{explicit NRSMat(const SparseMat&)} (requires the \subsection{Element insertion and access} \begin{lareference} \lafunc{inline void add(const SPMatindex n, const SPMatindex m, const T elem)}{Prepend element $(n,m)$; it is summed with existing entries at the same position. No checks, no \code{copyonwrite()}, see section~\ref{sec:sparse:refcount}. The matrix must have been dimensioned (constructor with dimensions or \code{resize}).} -\lafunc{void addsafe(const SPMatindex n, const SPMatindex m, const T elem)}{As \code{add}, but ignores $|$\code{elem}$|\le$ \code{SPARSEEPSILON}, calls \code{copyonwrite()} and creates the reference count of an undefined matrix (range check only with the lowercase \code{debug} macro).} -\lafunc{void simplify(const double sparseepsilon=SPARSEEPSILON)}{Sum duplicates, remove elements with $|a_{ij}|\le$ \code{sparseepsilon}, sort by rows; symmetric matrices are first moved to the upper triangle ($i\le j$).} +\lafunc{void addsafe(const SPMatindex n, const SPMatindex m, const T elem)}{As \code{add}, but ignores $|$\code{elem}$|\le\epsilon_{\rm sp}$, calls \code{copyonwrite()} and creates the reference count of an undefined matrix (range check only with the lowercase \code{debug} macro).} +\lafunc{void simplify(const typename LA_traits::normtype sparseepsilon=LA_precision::sparse_epsilon())}{Sum duplicates, remove elements with $|a_{ij}|\le$ \code{sparseepsilon}, sort by rows; symmetric matrices are first moved to the upper triangle ($i\le j$).} \lafunc{unsigned int length() const}{Number of list elements (cached until the next modification, see section~\ref{sec:sparse:storage}).} \lafunc{unsigned int sort(int type) const}{Build the internal array sorted by rows (\code{type=0}) or columns (\code{type=1}) and return the number of elements. Logically const, but unfolds a symmetric matrix (\code{setunsymmetric}). Normally called internally only.} \lafunc{inline matel *getlist() const; void setlist(matel *l)}{Raw access to the list head (no reference count handling; \code{setlist} does not free the old list).} @@ -245,7 +245,7 @@ copies). \lafunc{SparseSMat()}{Undefined empty matrix.} \lafunc{explicit SparseSMat(const SPMatindex n); explicit SparseSMat(const SPMatindex n, const SPMatindex m)}{Empty (zero) $n\times n$ or $n\times m$ matrix.} \lafunc{SparseSMat(const SparseSMat &rhs)}{Shallow copy.} -\lafunc{explicit SparseSMat(const NRSMat &rhs)}{From a packed symmetric (complex: hermitian) matrix (both triangles inserted, the upper one conjugated, $|a_{ij}|>$ \code{SPARSEEPSILON}).} +\lafunc{explicit SparseSMat(const NRSMat &rhs)}{From a packed symmetric (complex: hermitian) matrix (both triangles inserted, the upper one conjugated, $|a_{ij}|>$ $\epsilon_{\rm sp}$).} \lafunc{explicit SparseSMat(const NRMat &rhs)}{From a dense $n\times m$ matrix, nonzero elements copied as they are (no symmetrization).} \lafunc{explicit SparseSMat(const SparseMat &rhs)}{From a \code{SparseMat} (same dimensions); duplicates are summed, a symmetric \code{SparseMat} is unfolded without conjugation (complex symmetric, as \code{SparseMat} defines it).} \lafunc{explicit SparseSMat(const CSRMat &rhs)}{Declared but not implemented (\code{CSRMat} is not supported).} @@ -266,7 +266,7 @@ Dense from sparse: \code{explicit NRMat(const SparseSMat&)} copies all sto \lafunc{inline void add(const SPMatindex n, const SPMatindex m, const T elem, const bool both=true)}{Add \code{elem} to element $(n,m)$ and, if \code{both} and $n\ne m$, its complex conjugate to $(m,n)$. Creates the row map if needed; no \code{copyonwrite()}; range check in DEBUG only. Resulting zeros remain until \code{simplify()}.} \lafunc{inline const std::map *line(const SPMatindex n) const}{Read-only access to row $n$ (\code{NULL} if empty).} \lafunc{inline ROWTYPE & operator[](const SPMatindex i)}{Row $i$ as a reference to the pointer \code{std::map*} (\code{typedef ROWTYPE}), possibly \code{NULL}. Non-const only.} -\lafunc{unsigned long long simplify()}{Remove elements with $|a_{ij}|<$ \code{SPARSEEPSILON} and empty rows; returns the number of remaining elements.} +\lafunc{unsigned long long simplify()}{Remove elements with $|a_{ij}|<$ $\epsilon_{\rm sp}$ and empty rows; returns the number of remaining elements.} \lafunc{inline unsigned long long length()}{\code{simplify()} (non-const).} \lafunc{int nrows() const; int ncols() const; inline bool issymmetric() const}{Dimensions; \code{issymmetric()} is always true.} \lafunc{iterator begin() const; iterator end() const}{Iterator over all stored elements, row by row in increasing column order. \code{*it} is a \code{matel} \emph{copy} held in the iterator (writing to it does not change the matrix). It is meant for comparison with \code{end()} or via \code{it.notend()}; \code{iterator(const SparseSMat&)} starts at the first element.} @@ -288,7 +288,7 @@ Dense from sparse: \code{explicit NRMat(const SparseSMat&)} copies all sto \lafunc{const T* diagonalof(NRVec &r, const bool divide=0, bool cache=false) const}{As for \code{SparseMat} (square matrices only).} \lafunc{SparseSMat transpose(bool conj=false) const}{Transposed (optionally conjugated) copy.} \lafunc{void transposeme() const}{Raises an error (not implemented).} -\lafunc{SparseSMat cholesky(void) const}{Pivoted Cholesky decomposition of a positive semidefinite symmetric/hermitian matrix (not checked): returns $R$ with $R\Hc R = A$ (\code{result.transpose(true)*result} reproduces the input). The pivoting is done a priori by sorting the diagonal; pivots with real part below \code{CHOLESKYEPSILON}$=10^{-16}$ are skipped (rank-deficient case), elements below \code{SPARSEEPSILON} are dropped. $R$ is upper triangular only in the permuted basis.} +\lafunc{SparseSMat cholesky(void) const}{Pivoted Cholesky decomposition of a positive semidefinite symmetric/hermitian matrix (not checked): returns $R$ with $R\Hc R = A$ (\code{result.transpose(true)*result} reproduces the input). The pivoting is done a priori by sorting the diagonal; pivots with real part below \code{LA_precision::cholesky_epsilon()} ($10^{-16}$ for double, the old \code{CHOLESKYEPSILON}, $10^{-8}$ for float) are skipped (rank-deficient case), elements below $\epsilon_{\rm sp}$ are dropped. $R$ is upper triangular only in the permuted basis.} \lafunc{SparseSMat submatrix(const int fromrow, const int torow, const int fromcol, const int tocol) const}{Block (inclusive bounds).} \lafunc{void storesubmatrix(const int fromrow, const int fromcol, const SparseSMat &rhs)}{Store \code{rhs} into the block starting at $(\mathtt{fromrow},\mathtt{fromcol})$: the old elements of the block are removed first (as \code{NRMat::storesubmatrix}); calls \code{copyonwrite()}.} \end{lareference} @@ -302,7 +302,7 @@ Dense from sparse: \code{explicit NRMat(const SparseSMat&)} copies all sto \lafunc{SparseSMat otimes_sparse(const NRVec &lhs, const NRVec &rhs, const bool conjugate=false, const T &scale=1)}{\emph{Symmetrized} outer product, symmetric (complex: hermitian) regardless of the input vectors (unlike \code{SparseMat}): with $M=\mathtt{scale}\cdot l\,r\Tr$ (or $l\,r\Hc$ if \code{conjugate}) the result is $S=M+M\Hc$ ($=l\,r\Tr+r\,l\Tr$ for real types), including the doubled diagonal $S_{ii}=M_{ii}+M_{ii}^*$. Note that \code{otimes_sparse(x,x,true)} is $2\,x x\Hc$; use \code{scale=0.5} for $x x\Hc$. -Elements with $|S_{ij}|\le$ \code{SPARSEEPSILON} are omitted. Both vectors must have the same length (\code{laerror} otherwise).} +Elements with $|S_{ij}|\le$ $\epsilon_{\rm sp}$ are omitted. Both vectors must have the same length (\code{laerror} otherwise).} \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 diff --git a/doc/tensor.tex b/doc/tensor.tex index d870de5..a9c41fa 100644 --- a/doc/tensor.tex +++ b/doc/tensor.tex @@ -18,8 +18,8 @@ Typical uses are many-body quantities (amplitudes, integrals, density matrices) in groups of indices should be exploited for storage. Limitations stated in the source: \begin{itemize} -\item the implementation is ``rudimentary''; only \code{Tensor} and \code{Tensor >} - are instantiated (tensor.cc); +\item the implementation is ``rudimentary''; only \code{Tensor}, \code{Tensor >} and their single precision + counterparts \code{Tensor}, \code{Tensor >} are instantiated (tensor.cc); \item index types \code{LA_index} and \code{LA_largeindex} are \code{int}, so the total size is limited to $2^{31}$ elements (\code{NRVec} limitation); \item symmetry \emph{between} index groups (e.g.\ $(ij|kl)=(kl|ij)$) is not supported; one would have to @@ -291,7 +291,7 @@ Names are permuted together with the indices unless stated otherwise. \begin{lareference} \lafunc{void apply_permutation_algebra(const Tensor &rhs, const PermutationAlgebra &pa, bool inverse=false, T alpha=1, T beta=0, bool conjugate_by_parity=false)}{$\mathtt{*this} \leftarrow \beta\,\mathtt{*this}$, then for all stored elements $I$ of \code{*this}: $\mathtt{this}(I) \mathrel{+}= \alpha \sum_P c_P\, \mathtt{rhs}(P(I))$, with $P$ acting on the flat index list (inverse permutations if \code{inverse}); with \code{conjugate_by_parity} odd permutations take the complex conjugate. Can represent general (anti)symmetrizers, e.g.\ \code{general_antisymmetrizer} (chapter~\ref{ch:permutation}). \code{*this} may have a symmetry which \code{rhs} lacks. The names of \code{*this} are taken from (or checked against) the names of \code{rhs} permuted by the first permutation. Element-by-element, not optimally efficient.} \lafunc{void apply_permutation_algebra(const NRVec &rhsvec, const PermutationAlgebra &pa, bool inverse=false, T alpha=1, T beta=0)}{the same with \code{rhs} being the outer product of the tensors in \code{rhsvec}, which is not formed explicitly (e.g.\ Grassmann products building a density matrix from cumulants); rather inefficient.} -\lafunc{NRVec > Tucker(typename LA_traits::normtype thr=1e-12, bool inverseorder=false)}{HOSVD/Tucker decomposition: \code{*this} is replaced by the flat core tensor and the factor matrices (one per index, columns orthonormal) are returned; singular values below \code{thr} are dropped, which reduces the corresponding range of the core.} +\lafunc{NRVec > Tucker(typename LA_traits::normtype thr=LA_precision::svd_threshold(), bool inverseorder=false)}{HOSVD/Tucker decomposition: \code{*this} is replaced by the flat core tensor and the factor matrices (one per index, columns orthonormal) are returned; singular values below \code{thr} are dropped, which reduces the corresponding range of the core.} \lafunc{Tensor inverseTucker(const NRVec > &x, bool inverseorder=false) const}{rebuild the tensor from the core (\code{*this}) and the matrices returned by \code{Tucker} with the same \code{inverseorder}.} \lafunc{Tensor linear_transform(const NRVec > &x) const}{transform every index of group $g$ by the matrix \code{x[g]} (\code{x[g].ncols()} must equal the range; the new range is \code{x[g].nrows()}), preserving the group symmetries: $t'_{\dots a \dots} = \sum_i x_{ai}\, t_{\dots i\dots}$.} \lafunc{Tensor linear_transform(const int g, const Tensor &x) const; Tensor linear_transform(const int g, const NRMat &x) const}{transform only group \code{g}, keeping its position and symmetry; \code{x} must be a flat rank-2 tensor (useful for raising/lowering indices) or a matrix.}