diff --git a/CLAUDE.md b/CLAUDE.md index b051483..69251c3 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -14,8 +14,8 @@ decisions and plans agreed with the author, so that they survive a fresh clone. ## Build and test - Build the library with `make -j 8` (after header changes most files rebuild; run it in the background). -- Flags used: `-std=gnu++11 -O3 -DDEBUG -DNO_STRASSEN -DFORTRAN_ -DNONCLAPACK`, libs `-lcblas -lblas -llapack`. - Quick test program: `g++ -std=gnu++11 -I -O0 -DDEBUG -DNO_STRASSEN -DFORTRAN_ -DNONCLAPACK X.cc /.libs/libla.a -lcblas -lblas -llapack` +- Flags used: `-std=gnu++11 -O3 -DDEBUG -DFORTRAN_ -DNONCLAPACK`, libs `-lcblas -lblas -llapack`. + Quick test program: `g++ -std=gnu++11 -I -O0 -DDEBUG -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 "137 passed, 0 failed" (exit code = number of failures). New optimizer tests go to testoptimize.cc; @@ -84,7 +84,9 @@ decisions and plans agreed with the author, so that they survive a fresh clone. forbidden for NRMat, NRMat_from1, Mat3 (Mat3(double) works). - Real-storage gdiagonalize uses the dgeev convention (Re + i Im). - The NRMat offset view constructor is a feature. MATPTR is kept (rarely used). CSRMat is unfinished and - documented as unsupported, but kept. Strassen is experimental. CUDA (CUDALA) is deferred/experimental + documented as unsupported, but kept. NRMat::strassen (2026-11-06): native Strassen-Winograd in strassen.cc, gemm-compatible + parameters + cutoff (default LA_traits::strassen_cutoff()=2048 all types, to be tuned on HPC), explicit only (gemm stays classical), + CBLAS leaves on strided quadrants, device primitives separated for CUDA; NO_STRASSEN removed everywhere. CUDA (CUDALA) is deferred/experimental until the author has NVIDIA RTX hardware (known CUDA issues collected in suspected_bugs.md F10). - The regression module (reg.*, regsurf.*) is experimental, not much tested in practice (remarked in the manual); its NEA oscillator strength f=dE_min/dE is a placeholder for the true one. @@ -96,7 +98,7 @@ decisions and plans agreed with the author, so that they survive a fresh clone. 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. + float LM on badly scaled problems needs my_scaling=true. reg/regsurf, csrmat, CUDA stay double-only. Tensor Tucker bugs fixed 2026-11-06 (complex: factor = vt.transpose(false); unnamed tensors with inverseorder; bugtests/tensor/tucker.cc). BLAS-based NRVec members are meant for floating types only, not integers. - NRSMat: possible bool member `hermitian` toggling hermitian vs complex symmetric behaviour @@ -118,6 +120,8 @@ decisions and plans agreed with the author, so that they survive a fresh clone. all doc questions answered (2026-10-30), the manual has no open question boxes. - configure.ac and Makefile.am were changed 2026-10-30 (--enable-fortran64int, TESTS = testoptimize, targets manual/doc, doxygen; all-local builds only the manual (if pdflatex found), doxygen only on request; MANUAL_TEX lists the .tex files - add new chapters there): run autoreconf. +- configure.ac (2026-11-06): OpenBLAS (with LAPACK+CBLAS) preferred after MKL, only -lopenblas linked; --without-openblas; + not with --enable-fortran64int (its parsing moved up). Test programs then link with -lopenblas instead of -lcblas -lblas -llapack. - matexp.h: exp()/sincos() default horner=false (polynom, faster for dense); use horner=true for sparse matrices. - AuxStorage temporaries go to $TMPDIR (else cwd) and are unlinked at once; AuxStorage and DIIS are not copyable. - reg.cc/regsurf.cc cache raw pointers into their NRVec/NRMat members (set in Init): never reallocate those diff --git a/Makefile.am b/Makefile.am index 1dd3fcf..325cbf8 100644 --- a/Makefile.am +++ b/Makefile.am @@ -30,7 +30,7 @@ AM_CXXFLAGS= AM_CXXFLAGS += -g AM_CXXFLAGS += $(OPTIMIZEOPT) $(CUDAOPT) $(FORINTOPT) $(DEBUGOPT) $(PROFILEOPT) $(MATPTROPT) -AM_CXXFLAGS += -DNO_STRASSEN -DFORTRAN_ +AM_CXXFLAGS += -DFORTRAN_ AM_CXXFLAGS += $(CBLASOPT) $(CLAPACKOPT) AM_CXXFLAGS += $(MKLOPT) @@ -54,7 +54,7 @@ version.cc: always #documentation: make manual (alias doc) builds doc/manual.pdf, make doxygen the doxygen pages in doxygen-doc, #make all builds manual besides the library (each only if its tool was found by configure) MANUAL_TEX = doc/biblio.tex doc/bitvector.tex doc/building.tex doc/concepts.tex doc/fourindex.tex doc/introduction.tex doc/iterative.tex doc/lamacros.tex doc/manual.tex doc/mat.tex doc/nonclass.tex doc/optimize.tex doc/permutation.tex doc/polynomial.tex doc/quaternion.tex doc/regression.tex doc/smat.tex doc/sparse.tex doc/tensor.tex doc/utilities.tex doc/vec.tex doc/vecmat3.tex \ - doc/bib/iterative.tex doc/bib/optimize.tex doc/bib/permutation.tex doc/bib/polynomial.tex doc/bib/quaternion.tex doc/bib/utilities.tex doc/bib/vecmat3.tex + doc/bib/iterative.tex doc/bib/mat.tex doc/bib/optimize.tex doc/bib/permutation.tex doc/bib/polynomial.tex doc/bib/quaternion.tex doc/bib/utilities.tex doc/bib/vecmat3.tex EXTRA_DIST += doc/manual.pdf $(MANUAL_TEX) doc/manual.pdf: $(MANUAL_TEX) diff --git a/la_traits.h b/la_traits.h index fbbb756..8f77d4f 100644 --- a/la_traits.h +++ b/la_traits.h @@ -308,6 +308,7 @@ static inline C realpart(const std::complex &x) {return x.real();} static inline C imagpart(const std::complex &x) {return x.imag();} static inline void setrealpart(std::complex &x, const C &y) {reinterpret_cast(x)[0]=y;} static inline void setimagpart(std::complex &x, const C &y) {reinterpret_cast(x)[1]=y;} +static int strassen_cutoff() {return 2048;} //NRMat::strassen uses plain gemm for blocks with a dimension <= cutoff }; @@ -368,6 +369,7 @@ static inline C realpart(const C &x) {return x;} static inline C imagpart(const C &x) {return 0;} static inline void setrealpart(C &x, const C &y) {x=y;} static inline void setimagpart(C &x, const C &y) {} +static int strassen_cutoff() {return 2048;} //NRMat::strassen uses plain gemm for blocks with a dimension <= cutoff }; diff --git a/laerror.cc b/laerror.cc index 0b0758f..a4f84d8 100644 --- a/laerror.cc +++ b/laerror.cc @@ -149,7 +149,7 @@ throw LAerror(s1); } -//stub for f77 blas called from strassen routine +//stub for f77 blas error handler extern "C" void xerbla_(const char name[6], int *n){ char msg[1024]; strcpy(msg,"LAPACK or BLAS error in routine "); diff --git a/mat.h b/mat.h index 3009fb8..f0e6828 100644 --- a/mat.h +++ b/mat.h @@ -438,11 +438,8 @@ public: void printsorted(std::ostream &s, int direction=1, bool absvalue=false, bool stable=false, typename LA_traits::normtype thr=0, int topmax=0) const; -#ifndef NO_STRASSEN - //! Strassen's multiplication (better than \f$\mathacal{O}(n^3)\f$, analogous syntax to \see NRMat::gemm() ) - 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; -#endif + //! Strassen-Winograd multiplication (\f$\mathcal{O}(n^{2.81})\f$), the same operation as \see NRMat::gemm(); blocks with a dimension <= cutoff are multiplied by gemm + 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()); void swap(NRMat &rhs) //more efficient swap than via tmp and constructors and operator= { int tmpnn=nn; nn=rhs.nn; rhs.nn=tmpnn; diff --git a/strassen.cc b/strassen.cc index 9a350ea..2ed637e 100644 --- a/strassen.cc +++ b/strassen.cc @@ -19,32 +19,159 @@ #include "mat.h" -#ifndef NO_STRASSEN namespace LA { -/*Strassen algorithm*/ -// called routine is fortran-compatible -extern "C" void fmm(const char c_transa,const char c_transb,const int m,const int n,const int k,const double alpha, - const double *a,const int lda,const double *b,const int ldb,const double beta,double *c,const int ldc, - double *d_aux,int i_naux); -extern "C" void strassen_cutoff(int c, int c1, int c2, int c3); +/*Strassen algorithm in the Winograd variant (7 products, 15 additions), recursion on strided quadrants of the original matrices, +leaf products by BLAS gemm +V. Strassen, Gaussian elimination is not optimal, Numer. Math. 13 (1969) 354-356 +S. Winograd, On multiplication of 2x2 matrices, Linear Algebra Appl. 4 (1971) 381-388 +schedule with 2 temporaries per level: B. Boyer, J.-G. Dumas, C. Pernet, W. Zhou, Memory efficient scheduling of + Strassen-Winograd's matrix multiplication algorithm, Proc. ISSAC 2009, 55-62 +odd dimensions by dynamic peeling: S. Huss-Lederman, E. M. Jacobson, J. R. Johnson, A. Tsao, T. Turnbull, + Implementation of Strassen's algorithm for matrix multiplication, Proc. Supercomputing '96 +*/ -template<> -void NRMat::s_cutoff(const int c, const int c1, const int c2, const int c3) const -{ strassen_cutoff(c,c1,c2,c3);} +//device-dependent primitives: leaf gemm, quadrant addition, workspace (GPU version: cublas?gemm, cublas?geam, gpualloc/gpufree) +static inline CBLAS_TRANSPOSE s_trans(const char t) {return t=='n'?CblasNoTrans:(t=='c'?CblasConjTrans:CblasTrans);} -template<> -void NRMat::strassen(const double beta, const NRMat &a, const char transa, const NRMat &b, const char transb, const double alpha) +static inline void s_gemm(const char ta, const char tb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc) +{cblas_dgemm(CblasRowMajor,s_trans(ta),s_trans(tb),m,n,k,alpha,a,lda,b,ldb,beta,c,ldc);} +static inline void s_gemm(const char ta, const char tb, const int m, const int n, const int k, const float alpha, const float *a, const int lda, const float *b, const int ldb, const float beta, float *c, const int ldc) +{cblas_sgemm(CblasRowMajor,s_trans(ta),s_trans(tb),m,n,k,alpha,a,lda,b,ldb,beta,c,ldc);} +static inline void s_gemm(const char ta, const char tb, const int m, const int n, const int k, const std::complex alpha, const std::complex *a, const int lda, const std::complex *b, const int ldb, const std::complex beta, std::complex *c, const int ldc) +{cblas_zgemm(CblasRowMajor,s_trans(ta),s_trans(tb),m,n,k,&alpha,a,lda,b,ldb,&beta,c,ldc);} +static inline void s_gemm(const char ta, const char tb, const int m, const int n, const int k, const std::complex alpha, const std::complex *a, const int lda, const std::complex *b, const int ldb, const std::complex beta, std::complex *c, const int ldc) +{cblas_cgemm(CblasRowMajor,s_trans(ta),s_trans(tb),m,n,k,&alpha,a,lda,b,ldb,&beta,c,ldc);} + +//z = x+y or x-y, r x c blocks with leading dimensions; z may coincide with x or y (one pass, memory bound) +template +static void s_add(const int r, const int c, T *z, const int ldz, const T *x, const int ldx, const T *y, const int ldy, const bool minus) { -int l(transa=='n'?a.nn:a.mm); -int k(transa=='n'?a.mm:a.nn); -int kk(transb=='n'?b.nn:b.mm); -int ll(transb=='n'?b.mm:b.nn); - -if(l!=nn|| ll!=mm||k!=kk) laerror("incompatible (or undefined size) matrices in strassen"); - -copyonwrite(); -//swap transpositions and order of matrices -fmm(transb,transa,mm,nn,k,alpha,b,b.mm, a, a.mm, beta,*this, mm,NULL,0); +for(int i=0; i static T *s_alloc(const size_t n) {return n? new T[n] : NULL;} +template static void s_free(T *p) {if(p) delete[] p;} + +//operand op(M) of the product, op = identity ('n'), transposition ('t') or conjugate transposition ('c'), +//M stored row-major at p with leading dimension ld; the recursion is formulated in op() space and never transposes data: +//- at(i,j) is the address of element (i,j) of op(M), i.e. M(i,j) for 'n' and M(j,i) for 't','c'; +// a block of op(M) starting at (i,j) is then the block of M starting at at(i,j) with the same ld and flag +// (quadrant (op M)_{IJ} = op(M_{JI}), the stored quadrant being swapped for 't','c' and conjugated by gemm for 'c') +//- sums of quadrants are formed elementwise in the stored layout (S stored = M_{JI} + M_{LK} for 't','c'), which is correct +// because op(X+Y)=op(X)+op(Y), so the temporaries keep the flag of their operand and gemm applies op() at the leaves only +template +struct s_op + { + const T *p; int ld; char t; + const T *at(const int i, const int j) const {return t=='n' ? p+(size_t)i*ld+j : p+(size_t)j*ld+i;} + }; + +static inline bool s_split(const int m, const int n, const int k, const int cutoff) {return m>cutoff && n>cutoff && k>cutoff;} + +//workspace needed for all levels below (m,n,k) +static size_t s_work(const int m, const int n, const int k, const int cutoff) +{ +if(!s_split(m,n,k,cutoff)) return 0; +size_t mh=m/2, nh=n/2, kh=k/2; +return (mh*kh>mh*nh ? mh*kh : mh*nh) + kh*nh + s_work(mh,nh,kh,cutoff); +} + +//C(m x n, ldc) = alpha op(A) op(B) +template +static void s_rec(const int m, const int n, const int k, const T alpha, const s_op &A, const s_op &B, T *C, const int ldc, T *ws, const int cutoff) +{ +if(!s_split(m,n,k,cutoff)) {s_gemm(A.t,B.t,m,n,k,alpha,A.p,A.ld,B.p,B.ld,(T)0,C,ldc); return;} +const int mh=m/2, nh=n/2, kh=k/2; +//stored dimensions of the quadrants of A and B +const int ar= A.t=='n'?mh:kh, ac= A.t=='n'?kh:mh; +const int br= B.t=='n'?kh:nh, bc= B.t=='n'?nh:kh; +const T *A11=A.at(0,0), *A12=A.at(0,kh), *A21=A.at(mh,0), *A22=A.at(mh,kh); +const T *B11=B.at(0,0), *B12=B.at(0,nh), *B21=B.at(kh,0), *B22=B.at(kh,nh); +T *C11=C, *C12=C+nh, *C21=C+(size_t)mh*ldc, *C22=C21+nh; +const size_t xs= (size_t)mh*kh>(size_t)mh*nh ? (size_t)mh*kh : (size_t)mh*nh; +T *X=ws, *Y=ws+xs, *wn=Y+(size_t)kh*nh; +const s_op sx={X,ac,A.t}, sy={Y,bc,B.t}; //X, Y hold sums of quadrants of A, B in their stored layout, contiguous (ld = stored columns), same flag +const s_op a11={A11,A.ld,A.t}, a12={A12,A.ld,A.t}, a22={A22,A.ld,A.t}; +const s_op b11={B11,B.ld,B.t}, b21={B21,B.ld,B.t}, b22={B22,B.ld,B.t}; +s_add(ar,ac,X,ac,A11,A.ld,A21,A.ld,true); //S3=A11-A21 +s_add(br,bc,Y,bc,B22,B.ld,B12,B.ld,true); //T3=B22-B12 +s_rec(mh,nh,kh,alpha,sx,sy,C21,ldc,wn,cutoff); //P7 +s_add(ar,ac,X,ac,A21,A.ld,A22,A.ld,false); //S1=A21+A22 +s_add(br,bc,Y,bc,B12,B.ld,B11,B.ld,true); //T1=B12-B11 +s_rec(mh,nh,kh,alpha,sx,sy,C22,ldc,wn,cutoff); //P5 +s_add(ar,ac,X,ac,X,ac,A11,A.ld,true); //S2=S1-A11 +s_add(br,bc,Y,bc,B22,B.ld,Y,bc,true); //T2=B22-T1 +s_rec(mh,nh,kh,alpha,sx,sy,C12,ldc,wn,cutoff); //P6 +s_add(ar,ac,X,ac,A12,A.ld,X,ac,true); //S4=A12-S2 +s_rec(mh,nh,kh,alpha,sx,b22,C11,ldc,wn,cutoff); //P3 +s_rec(mh,nh,kh,alpha,a11,b11,X,nh,wn,cutoff); //P1 in X +s_add(mh,nh,C12,ldc,X,nh,C12,ldc,false); //U2=P1+P6 +s_add(mh,nh,C21,ldc,C12,ldc,C21,ldc,false); //U3=U2+P7 +s_add(mh,nh,C12,ldc,C12,ldc,C22,ldc,false); //U4=U2+P5 +s_add(mh,nh,C22,ldc,C21,ldc,C22,ldc,false); //C22=U3+P5 +s_add(mh,nh,C12,ldc,C12,ldc,C11,ldc,false); //C12=U4+P3 +s_add(br,bc,Y,bc,Y,bc,B21,B.ld,true); //T4=T2-B21 +s_rec(mh,nh,kh,alpha,a22,sy,C11,ldc,wn,cutoff); //P4 +s_add(mh,nh,C21,ldc,C21,ldc,C11,ldc,true); //C21=U3-P4 +s_rec(mh,nh,kh,alpha,a12,b21,C11,ldc,wn,cutoff); //P2 +s_add(mh,nh,C11,ldc,X,nh,C11,ldc,false); //C11=P1+P2 +//odd dimensions +const int me=2*mh, ne=2*nh, ke=2*kh; +if(k>ke) s_gemm(A.t,B.t,me,ne,1,alpha,A.at(0,ke),A.ld,B.at(ke,0),B.ld,(T)1,C,ldc); +if(n>ne) s_gemm(A.t,B.t,me,1,k,alpha,A.p,A.ld,B.at(0,ne),B.ld,(T)0,C+ne,ldc); +if(m>me) s_gemm(A.t,B.t,1,n,k,alpha,A.at(me,0),A.ld,B.p,B.ld,(T)0,C+(size_t)me*ldc,ldc); +} + +//C = beta C + alpha op(A) op(B), row-major m x n +template +static void s_top(const int m, const int n, const int k, const T alpha, const T *a, const int lda, const char ta, const T *b, const int ldb, const char tb, const T beta, T *c, const int cutoff) +{ +const s_op A={a,lda,ta}, B={b,ldb,tb}; +if(!s_split(m,n,k,cutoff)) {s_gemm(ta,tb,m,n,k,alpha,a,lda,b,ldb,beta,c,n); return;} +const size_t w=s_work(m,n,k,cutoff); +const bool acc= beta!=(T)0; +T *ws=s_alloc(w+(acc?(size_t)m*n:0)); +T *cc= acc? ws+w : c; +s_rec(m,n,k,alpha,A,B,cc,n,ws,cutoff); +if(acc) for(size_t i=0; i<(size_t)m*n; ++i) c[i]=beta*c[i]+cc[i]; +s_free(ws); +} + +static inline char s_flag(const char t, const bool cmplx) {char tt=tolower(t); return tt=='n'?'n':(tt=='c'&&cmplx?'c':'t');} + +template +static void s_strassen(NRMat &c, const T &beta, const NRMat &a, const char transa, const NRMat &b, const char transb, const T &alpha, const int cutoff) +{ +const char ta=s_flag(transa,LA_traits::is_complex()), tb=s_flag(transb,LA_traits::is_complex()); +int k(ta=='n'?a.ncols():a.nrows()); +#ifdef DEBUG +int l(ta=='n'?a.nrows():a.ncols()); +int kk(tb=='n'?b.nrows():b.ncols()); +int ll(tb=='n'?b.ncols():b.nrows()); +if (l!=c.nrows() || ll!=c.ncols() || k!=kk) laerror("incompatible matrices in NRMat::strassen(...)"); +if(b.ncols() <=0 || c.ncols()<=0) laerror("illegal matrix dimension in strassen"); #endif +SAME_LOC3(c, a, b); +NOT_GPU(c); +if(cutoff<1) laerror("illegal cutoff in strassen"); +if(alpha==(T)0 && beta==(T)1) return; +c.copyonwrite(); +s_top(c.nrows(),c.ncols(),k,alpha,(const T *)a,a.ncols(),ta,(const T *)b,b.ncols(),tb,beta,(T *)c,cutoff); +} + +#define INSTANTIZE(T) \ +template<> \ +void NRMat::strassen(const T &beta, const NRMat &a, const char transa, const NRMat &b, const char transb, const T &alpha, const int cutoff) \ +{s_strassen(*this,beta,a,transa,b,transb,alpha,cutoff);} + +INSTANTIZE(double) +INSTANTIZE(float) +INSTANTIZE(std::complex) +INSTANTIZE(std::complex) + +}//namespace diff --git a/t.cc b/t.cc index d75b9d1..65c42c4 100644 --- a/t.cc +++ b/t.cc @@ -324,16 +324,11 @@ cout <<" big1*big2 "< atest, btest,ctest; -{ -int cc,c1,c2,c3; -cin >>cc>>c1>>c2>>c3; - -atest.s_cutoff(cc,c1,c2,c3); -} +int cc; +cin >>cc; cin>>atest; cin>>btest; @@ -342,7 +337,7 @@ dtest.gemm(0., atest, 't', btest, 'n', 1.); cout << dtest; NRMat etest(atest.nrows(),btest.ncols()); -etest.strassen(0., atest, 't', btest, 'n', 1.); +etest.strassen(0., atest, 't', btest, 'n', 1., cc); cout << etest; } @@ -363,13 +358,11 @@ cout <<" classical big1*big2 "<