.gitignore

@@ -6,6 +6,5 @@
 *.pdf
 *.eps
 cmake-build-*
-build
 *.idea
 .DS_Store

.travis.yml

@@ -0,0 +1,27 @@
+matrix:
+  include:
+  - language: cpp
+    sudo: enabled
+    before_install:
+      - sudo apt-get update
+      - sudo apt-get install -y gcc liblapack-dev cmake libfftw3-dev python-dev
+        #- sudo apt-get install -y gcc liblapack-dev cmake libfftw3-dev libgtest-dev python-dev
+      - pip install numpy --user
+      - pip install scipy --user
+      - sudo apt-get install -y libopenmpi-dev openmpi-bin
+        #- install -d gtest_build
+        #- pushd gtest_build
+        #- cmake -DCMAKE_CXX_COMPILER=g++ /usr/src/gtest/
+        #- make
+        #- sudo cp *.a /usr/lib
+        #- popd 
+
+    install:
+      - mkdir build
+      - cd build
+      - cmake -DTesting=ON ../
+      - VERBOSE=1 make 
+    script:
+      - export OMP_NUM_THREADS=1
+      - make test ARGS='-V'
+

CMakeLists.txt

@@ -3,7 +3,7 @@ cmake_minimum_required(VERSION 2.8.4 FATAL_ERROR)
 project(sparse CXX C)
 
 # add common modules from ../common/cmake
-list(APPEND CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/common/cmake)
+list(APPEND CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake)
 
 if(NOT CMAKE_BUILD_TYPE)
   set(CMAKE_BUILD_TYPE "Release" CACHE STRING "Type of build" FORCE)
@@ -13,21 +13,48 @@ message(STATUS "Build type: " ${CMAKE_BUILD_TYPE})
 
 set(CMAKE_CXX_FLAGS_DEBUG "-std=c++11 -D_DEBUG -g -Wall -Wextra -pedantic -Wcast-align -Wcast-qual -Wconversion -Wdisabled-optimization -Wendif-labels -Wfloat-equal -Winit-self -Winline -Wlogical-op -Wmissing-include-dirs -Wnon-virtual-dtor -Wold-style-cast -Woverloaded-virtual -Wpacked -Wpointer-arith -Wredundant-decls -Wshadow -Wsign-promo -Wswitch-default -Wswitch-enum -Wunsafe-loop-optimizations -Wvariadic-macros -Wwrite-strings ")
 
-set(CMAKE_CXX_FLAGS_RELEASE "-std=c++11 -O3 -fopenmp")
+set(CMAKE_CXX_FLAGS_RELEASE "-std=c++11 -O3")
 
 #include directories
-include_directories(${CMAKE_SOURCE_DIR}/include ${CMAKE_SOURCE_DIR}/thirdparty/cpplapack/include)
-
-find_package(OpenMP)
-if(OPENMP_FOUND)
-  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
-endif(OPENMP_FOUND)
-
-#Find BLAS
-find_package(BLAS REQUIRED)
-
-#Find LAPACK
-find_package(LAPACK REQUIRED)
+include_directories(${CMAKE_SOURCE_DIR}/thirdparty/cpplapack/include)
+
+#find_package(OpenMP)
+#if(OPENMP_FOUND)
+#  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
+#endif(OPENMP_FOUND)
+
+if(USE_MKL)
+  if(NOT MKLROOT)
+    set(MKLROOT $ENV{MKLROOT})
+  endif()
+  if(NOT MKLROOT)
+    message(FATAL "MKLROOT is not set")
+  endif()
+
+  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -mkl")
+  include_directories(${MKLROOT}/include/fftw)
+else()
+  #Find BLAS
+  find_package(BLAS REQUIRED)
+
+  #Find LAPACK
+  find_package(LAPACK REQUIRED)
+
+  #Find FFTW3
+  find_package(FFTW REQUIRED)
+  include_directories(${FFTW_INCLUDE_DIRS})
+endif()
+
+# Build and enable tests
+# testing setup
+# enable_testing() must be called in the top-level CMakeLists.txt before any add_subdirectory() is called.
+option(Testing "Enable testing" ON)
+
+if (Testing)
+  enable_testing(test)
+  add_subdirectory(test/c++)
+  add_subdirectory(test/python)
+endif()
 
 #Find SWIG
 #find_package(SWIG REQUIRED)
@@ -52,6 +79,4 @@ find_package(LAPACK REQUIRED)
 # CPPLAPACK
 #include_directories(${CPPLAPACK_INCLUDE_DIR})
 
-add_subdirectory(src)
-
-#add_subdirectory(python)
+add_subdirectory(c++)

README.md

@@ -1,8 +1,6 @@
+# SpM
 
-SPM
-====
-
-Sparse Modeling tool for analytical continuation.
+Sparse Modeling tool for analytical continuation from the imaginary-time Green's function to real-requency spectral function.
 
 The algorithm is presented in the article
 
@@ -14,55 +12,65 @@ The algorithm is presented in the article
 
 This package is distributed under GNU General Public License version 3 ([GPL v3](http://www.gnu.org/licenses/gpl-3.0.en.html)).
 
-We kindly ask you to cite the article above
+We kindly ask you to cite the below articles
 in publications that include results obtained using this package.
 
-## Author
-Junya Otsuki, Kazuyoshi Yoshimi, Hiroshi Shinaoka, Masayuki Ohzeki
+* The article about original algorithm
+    * J. Otsuki, M. Ohzeki, H. Shinaoka, K. Yoshimi,  
+    "Sparse modeling approach to analytical continuation of imaginary-time quantum Monte Carlo data"  
+    [Phys. Rev. E 95, 061302(R) (2017).](https://doi.org/10.1103/PhysRevE.95.061302)
+* The article about this package
+    * K. Yoshimi, J. Otsuki, Y. Motoyama, M. Ohzeki, and H. Shinaoka, "SpM: Sparse modeling tool for analytic continuation of imaginary-time Green's function" [Comput. Phys. Commun. 244, 319-323 (2019)](https://www.sciencedirect.com/science/article/pii/S0010465519302103).
+* The article about the SpM-Pade method (please cite if you use)
+    * Y. Motoyama, K. Yoshimi, and J. Otsuki, "Robust analytic continuation combining the advantages of the sparse modeling approach and Pade approximation" [arXiv:2109.08370](https://arxiv.org/abs/2109.08370)
+
+## Authors
+
+Junya Otsuki, Kazuyoshi Yoshimi, Yuichi Motoyama, Hiroshi Shinaoka, Masayuki Ohzeki
 
 ## Requirement
 
 * LAPACK, BLAS
+* FFTW3
 * cpplapack (included in this package)
 
-
 ## How to build
 
 ### Getting the source codes
 
 Download the latest source codes by
 
-	$ git clone https://github.com/j-otsuki/SpM.git spm.src
+    $ git clone https://github.com/j-otsuki/SpM.git spm.src
 
-Then, the source codes are downloaded in the directory ``spm.src``.
+Then, the source codes are downloaded in the directory `spm.src`.
 
 ### Using Cmake
 
-Build with cmake command is done in a separate directory, e.g. ``spm.build``.
+Build with cmake command is done in a separate directory, e.g. `spm.build`.
 Type the following commands:
 
-	$ mkdir spm.build && cd spm.build
-	$ cmake ../spm.src
-	$ make
+    $ mkdir spm.build && cd spm.build
+    $ cmake ../spm.src
+    $ make
 
-Then, the executable file ``SpM.out`` is created in directory ``spm.build/src``.
+Then, the executable file `SpM.out` is created in directory `spm.build/src`.
 
 ## Sample scripts
 
-Some sample data are provided in ``samples`` directory:
+Some sample data are provided in `samples` directory:
 
-* ``samples/fermion``  # sample for fermionic spectrum (data in the article)
-* ``samples/boson``  # sample for bosonic spectrum
+* `samples/fermion`  # sample for fermionic spectrum (data in the article)
+* `samples/boson`  # sample for bosonic spectrum
 
 A script file, `run.sh`, is provided to run through the program.
 Enter into the directory `samples/fermion`, and execute the script by
 
     $ ./run.sh
 
-You may need to change the path to ``SpM.out`` in the script.
-If succeeded, results including graphs in eps format are created in ``output`` directory.
+You may need to change the path to `SpM.out` in the script.
+If succeeded, results including graphs in eps format are created in `output` directory.
 For details, see the document linked below.
 
-
 ## Official page
+
 The official page of the SpM is [here](https://spm-lab.github.io/SpM/manual/build/html/index.html).  

c++/CMakeLists.txt

@@ -0,0 +1,3 @@
+#include directories
+include_directories(include)
+add_subdirectory(src)

c++/include/Gf.h

@@ -0,0 +1,55 @@
+#ifndef _GF_H
+#define _GF_H
+
+#include <vector>
+#include <complex>
+#include "pade.h"
+
+class Gf{
+public:
+    Gf(): flag_pade(false) {};
+
+    enum statistics{ FERMION, BOSON };
+
+    // initialize by G(tau). tau=[0:beta), N points (N=2^integer is recommended)
+    void set_Gtau(const std::vector<double> &G_tau, const statistics stat, const double beta, const double tail=1.);
+
+    // initialize by G(iw_n). n=[0:N/2)
+    void set_Giw(const std::vector< std::complex<double> > &G_iw, const statistics stat, const double beta, const double tail=1.);
+
+    int get_size() { return N; };
+    double get_tail() { return tail; };
+    double get_beta() { return beta; };
+    statistics get_statistics() { return stat; };
+
+    // get G(tau)
+    void get_Gtau(std::vector<double> &G_tau);
+
+    // get G(iw)
+    void get_Giw(std::vector< std::complex<double> > &G_iw);
+
+    // get G(w). w=[omega_min:omega_max], n-points
+    // void get_Gomega(std::vector< std::complex<double> > &G_omega, const double omega_min, const double omega_max, const int n);
+
+    // return G(omega)
+    std::complex<double> Gomega(const double omega);
+
+    // return rho(omega)
+    double rho(const double omega);
+
+private:
+    int N; // == beta/dt
+    double beta;
+    double tail;
+    statistics stat;
+    std::vector<double> gtau;  // size N
+    std::vector< std::complex<double> > giw;  // size N/2
+
+    bool flag_pade;
+    Pade pade;
+
+    void compute_tau2iw();  // G(tau) --> G(iw)
+    void compute_iw2tau();  // G(iw) --> G(tau)
+};
+
+#endif // _GF_H

include/SVD_matrix.h → c++/include/SVD_matrix.h

@@ -40,10 +40,10 @@ class SVD_matrix {
     void rearrange_col(std::vector<int> &);  // rearrange column of A (actaully VT)
     void print_basis(std::string _file_U, std::string _file_V, int);
 
-    CPPL::dcovector transform_x2sv(const CPPL::dcovector &);  // V^t * x  (to original basis)
-    CPPL::dcovector transform_y2sv(const CPPL::dcovector &);  // U^t * y  (to original basis)
-    CPPL::dcovector transform_sv2x(const CPPL::dcovector &);  // V * x'  (to SV basis)
-    CPPL::dcovector transform_sv2y(const CPPL::dcovector &);  // U * y'  (to SV basis)
+    CPPL::dcovector transform_x2sv(const CPPL::dcovector &);  // V^t * x  (to SV basis)
+    CPPL::dcovector transform_y2sv(const CPPL::dcovector &);  // U^t * y  (to SV basis)
+    CPPL::dcovector transform_sv2x(const CPPL::dcovector &);  // V * x'  (to original basis)
+    CPPL::dcovector transform_sv2y(const CPPL::dcovector &);  // U * y'  (to original basis)
     void OutputSVD(std::string _file_SVD);
 private:
     CPPL::dcovector S_temp;

include/admm_svd.h → c++/include/admm_svd.h

@@ -35,12 +35,14 @@ class SVD_matrix;
 struct admm_result{
 	// w/ sv  : results in SV basis
 	// w/o sv : results in omega-tau basis
-	std::vector<double> x, xsv, z1, z1sv, z2, z2sv;
+	std::vector<double> x, xsv, z1, z1sv, z2, z2sv, z3, z3sv;
 	std::vector<double> y, ysv, y_recovered_x, y_recovered_z1, ysv_recovered_x, ysv_recovered_z1;
 };
 struct admm_info{
 	double res1_pri, res1_dual, res2_pri, res2_dual;  // residual errors
+	double res3_pri, res3_dual; // spmpade
 	double mse, mse_full;
+	int l0_norm;
 	double l1_norm, sum_x_calc, negative_weight;
 	int iter;
 };
@@ -52,13 +54,16 @@ class admm_svd{
 
 	// [optional]
 	void set_sumrule(double sum_x);
+	void set_rho_ref(std::vector<double>const &rho_ref, std::vector<double> const & rho_ref_weight);
+	void set_omega_coeff(std::vector<double> const& omega_coeff);
 	void set_nonnegative(bool _flag);
 	void set_fileout_iter(const std::string filename);  // output convergence in a file.  unset if filename=""
 	void set_print_level(int);  // 0: none,  1: results,  2: verbose
 
 	// [required]
-	void set_coef(double lambda, double penalty1=1.0, double penalty2=1.0, bool flag_penalty_auto=false);
+	void set_coef(double lambda, double penalty1=1.0, double penalty2=1.0, double penalty3=1.0, bool flag_penalty_auto=false);
 	void set_y(const std::vector<double> &y);
+	int get_svd_num(){return svd_num;}
 
 	// [optional]
 	void clear_x();
@@ -77,28 +82,42 @@ class admm_svd{
 
 	bool flag_nonnegative;
 	bool flag_sumrule;
+	bool flag_rho_ref;
 	double sum_x;
+	std::vector<double> rho_ref;
+	std::vector<double> rho_ref_weight;
 
-	CPPL::dcovector x, Vx, z1, u1, z2, u2;  // 1 for L1 norm, 2 for non-negativity
+	CPPL::dcovector x, Vx, WVx, z1, u1, z2, u2;  // 1 for L1 norm, 2 for non-negativity
+	CPPL::dcovector z3, u3; // SpM-Pade
 	CPPL::dcovector y, y_sv;
 
 	double regulariz, penalty1, penalty2;
+	double penalty3; // SpMPade
 	bool flag_penalty_auto;
 	static const int PENALTY_UPDATE_INTERVAL;
 
+	std::vector<double> omega_coeff; // omega
+
 	int print_level;
 	bool flag_fileout_iter;
 	std::string file_iter;
 
+	int svd_num;
 	void pre_update();  // This function must be called before calling update_x, and must be recalled when one of values of lambda, penalty1, penalty2 are changed
 	void update_x();
 
+	CPPL::dcovector transform_x2sv(CPPL::dcovector const& v);
+	CPPL::dcovector transform_sv2x(CPPL::dcovector const& v);
+
 	// quantities used in functions update_x and set_y (set in function pre_update)
 	struct quantities_for_update{
-		CPPL::dgbmatrix Y, B;  // diagonal matrix
-		CPPL::dgematrix C;
+		CPPL::dgematrix Y, B;  // diagonal matrix
+		CPPL::dgbmatrix Pade_B3;
+		CPPL::dgematrix Pade_B0, C;
 		CPPL::dcovector Yy, w;
+		CPPL::dcovector Pade_y;
 		CPPL::drovector v_row;
+        CPPL::dcovector rho_w;
 		double sum_Vw;
 	} pre;
 };

include/common.h → c++/include/common.h

@@ -24,6 +24,7 @@ CPPL::dcovector vec2cppl_col(const std::vector<double> &);  // vector -> CPPL
 std::vector<double> cppl2vec(const CPPL::dcovector &);  // CPPL -> vector
 std::vector<double> cppl2vec(const CPPL::drovector &);  // CPPL -> vector
 
+int norm_l0(const CPPL::dcovector &);  // L0 norm
 double norm_l1(const CPPL::dcovector &);  // L1 norm
 double norm_l2(const CPPL::dcovector &);  // L2 norm
 double norm_l2_sq(const CPPL::dcovector &);  // square of L2 norm
@@ -39,6 +40,11 @@ CPPL::drovector drovector_all1(int n);  // return a rovector with elements all b
 CPPL::dcovector positive(const CPPL::dcovector &);  // set negative elements at zero
 CPPL::dgematrix positive(const CPPL::dgematrix &);
 
+// return (v + coeff*ref) / (1+coeff)
+CPPL::dcovector project_ref(const CPPL::dcovector &v,
+                            const std::vector<double> &ref,
+                            const std::vector<double> &coeff); 
+
 // return shrinked matrix
 CPPL::dgematrix low_rank_matrix(CPPL::dgematrix &mat, int m, int n);