@@ -59,16 +59,16 @@ class UserOperator : public VirtualGenerator<double> {
5959int main (int argc, char *argv[]) {
6060
6161 // Check the number of parameters
62- if (argc != 2 ) {
62+ if (argc > 2 ) {
6363 // Tell the user how to run the program
64- cerr << " Usage: " << argv[0 ] << " outputpath " << endl;
64+ cerr << " Usage: " << argv[0 ] << " output_folder " << endl;
6565 /* "Usage messages" are a conventional way of telling the user
6666 * how to run a program if they enter the command incorrectly.
6767 */
6868 return 1 ;
6969 }
7070
71- std::string outputpath = argv[1 ];
71+ std::string output_folder = argc == 2 ? argv[1 ] : " ./ " ;
7272
7373 // Execution policy
7474 auto &policy = exec_compat::par;
@@ -77,15 +77,16 @@ int main(int argc, char *argv[]) {
7777 const int number_points = 10000 ;
7878 const int spatial_dimension = 3 ;
7979 vector<double > coordinates (spatial_dimension * number_points);
80- create_sphere ( number_points, coordinates.data ());
80+ create_rotated_ellipse ( 3 , 4 ., 1 ., 0 ., 0 ., number_points, coordinates.data ()); // 2d ellipse in 3d
8181
8282 // Cluster tree builder
8383 ClusterTreeBuilder<double > recursive_build_strategy;
8484 recursive_build_strategy.set_maximal_leaf_size (500 );
85+ recursive_build_strategy.set_partitioning_strategy (std::make_shared<Partitioning_N<double , ComputeLargestExtent<double >, RegularSplitting<double >>>());
8586
8687 // HMatrix parameters
8788 const double epsilon = 0.01 ;
88- const double eta = 10 ;
89+ const double eta = 200 ;
8990 char symmetry = ' S' ;
9091 char UPLO = ' L' ;
9192
@@ -97,18 +98,18 @@ int main(int argc, char *argv[]) {
9798 HMatrix<double > hmatrix = hmatrix_builder.build (policy, A, htool::HMatrixTreeBuilder<double >(epsilon, eta, symmetry, UPLO));
9899
99100 // Output
100- // save_leaves_with_rank(hmatrix, outputpath + "hmatrix");
101+ save_leaves_with_rank (hmatrix, output_folder + " / hmatrix" );
101102 print_tree_parameters (hmatrix, std::cout);
102103 print_hmatrix_information (hmatrix, std::cout);
103104
104- // sequential y= A*x
105+ // y= A*x
105106 std::vector<double > x (number_points, 1 ), y (number_points, 0 ), ref (number_points, 0 );
106107 add_hmatrix_vector_product (policy, ' N' , double (1 ), hmatrix, x.data (), double (0 ), y.data ());
107108 ref = A * x;
108109 std::cout << " relative error on matrix vector product : " ;
109110 std::cout << norm2 (ref - y) / norm2 (ref) << " \n " ;
110111
111- // sequential z = A^-1 y
112+ // z = A^-1 y
112113 std::vector<double > z (number_points, 0 );
113114 z = y;
114115 MatrixView<double > z_view (number_points, 1 , z.data ());
0 commit comments