1 // -*- C++ -*-
2 /**
3 * @file HMatrixFactory.cxx
4 * @brief This file supplies support for HMat
5 *
6 * Copyright 2005-2021 Airbus-EDF-IMACS-ONERA-Phimeca
7 *
8 * This library is free software: you can redistribute it and/or modify
9 * it under the terms of the GNU Lesser General Public License as published by
10 * the Free Software Foundation, either version 3 of the License, or
11 * (at your option) any later version.
12 *
13 * This library is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU Lesser General Public License for more details.
17 *
18 * You should have received a copy of the GNU Lesser General Public License
19 * along with this library. If not, see <http://www.gnu.org/licenses/>.
20 */
21 #include "openturns/HMatrixFactory.hxx"
22 #include "openturns/HMatrix.hxx"
23 #include "openturns/HMatrixImplementation.hxx"
24 #include "openturns/Sample.hxx"
25 #include "openturns/Log.hxx"
26
27 #ifdef OPENTURNS_HAVE_HMAT
28 #include <hmat/config.h>
29 #ifdef HMAT_HAVE_STARPU
30 # include <hmat/hmat_parallel.h>
31 #else
32 # include <hmat/hmat.h>
33 #endif
34 #endif /* OPENTURNS_HAVE_HMAT */
35
36 BEGIN_NAMESPACE_OPENTURNS
37
CLASSNAMEINIT(HMatrixFactory)38 CLASSNAMEINIT(HMatrixFactory)
39
40 /* Default constructor */
41 HMatrixFactory::HMatrixFactory()
42 : PersistentObject()
43 {
44 // Nothing to do
45 }
46
clone() const47 HMatrixFactory * HMatrixFactory::clone() const
48 {
49 return new HMatrixFactory( *this );
50 }
51
52
53 /** Tell whether HMat support is available */
54 Bool
IsAvailable()55 HMatrixFactory::IsAvailable()
56 {
57 #ifdef OPENTURNS_HAVE_HMAT
58 return true;
59 #else
60 return false;
61 #endif
62 }
63
64
65 HMatrix
build(const Sample & sample,UnsignedInteger outputDimension,Bool symmetric,const HMatrixParameters & parameters)66 HMatrixFactory::build(const Sample & sample, UnsignedInteger outputDimension, Bool symmetric, const HMatrixParameters & parameters)
67 {
68 #ifndef OPENTURNS_HAVE_HMAT
69 throw NotYetImplementedException(HERE) << "OpenTURNS has been built without HMat support";
70 #else
71 hmat_interface_t *hmatInterface = (hmat_interface_t*) calloc(1, sizeof(hmat_interface_t));
72 hmat_settings_t settings;
73
74 #ifdef HMAT_HAVE_STARPU
75 if (!ResourceMap::GetAsBool("HMatrix-ForceSequential"))
76 hmat_init_starpu_interface(hmatInterface, HMAT_DOUBLE_PRECISION);
77 else
78 #endif
79 hmat_init_default_interface(hmatInterface, HMAT_DOUBLE_PRECISION);
80
81 hmat_get_parameters(&settings);
82
83 settings.maxLeafSize = ResourceMap::GetAsUnsignedInteger("HMatrix-MaxLeafSize");
84 settings.validationErrorThreshold = ResourceMap::GetAsScalar("HMatrix-ValidationError");
85 settings.validateCompression = settings.validationErrorThreshold > 0;
86 settings.validationReRun = ResourceMap::GetAsUnsignedInteger("HMatrix-ValidationRerun");
87 settings.validationDump = ResourceMap::GetAsUnsignedInteger("HMatrix-ValidationDump");
88
89 hmat_set_parameters(&settings);
90
91 if (0 != hmatInterface->init())
92 {
93 throw NotYetImplementedException(HERE) << "Unable to initialize HMat library";
94 }
95
96 const UnsignedInteger size = sample.getSize();
97 const UnsignedInteger inputDimension = sample.getDimension();
98 double* points = new double[inputDimension * outputDimension * size];
99 for (UnsignedInteger i = 0; i < size; ++i)
100 {
101 for (UnsignedInteger j = 0; j < outputDimension; ++j)
102 std::copy(&sample(i, 0), &sample(i, 0) + inputDimension, points + i * inputDimension * outputDimension + j * inputDimension);
103 }
104
105 hmat_clustering_algorithm_t* algo;
106 const String clusteringAlgorithm = parameters.getClusteringAlgorithm();
107 if (clusteringAlgorithm == "median")
108 algo = hmat_create_clustering_median();
109 else if (clusteringAlgorithm == "geometric")
110 algo = hmat_create_clustering_geometric();
111 else if (clusteringAlgorithm == "hybrid")
112 algo = hmat_create_clustering_hybrid();
113 else
114 throw InvalidArgumentException(HERE) << "Unknown clustering method: " << clusteringAlgorithm << ", valid choices are: median, geometric or hybrid";
115
116 hmat_cluster_tree_t* ct = hmat_create_cluster_tree(points, inputDimension, outputDimension * size, algo);
117 hmat_delete_clustering(algo);
118 delete[] points;
119
120 Scalar eta = parameters.getAdmissibilityFactor();
121 hmat_admissibility_t* admissibility = hmat_create_admissibility_standard(eta);
122 hmat_matrix_t* ptrHMat = hmatInterface->create_empty_hmatrix_admissibility(ct, ct, symmetric, admissibility);
123 hmat_delete_admissibility(admissibility);
124 return HMatrix(new HMatrixImplementation(hmatInterface, ct, outputDimension * size, ptrHMat));
125 #endif /* OPENTURNS_HAVE_HMAT */
126 }
127
__repr__() const128 String HMatrixFactory::__repr__() const
129 {
130 return OSS() << "class=" << getClassName();
131 }
132
133 /* Method save() stores the object through the StorageManager */
save(Advocate & adv) const134 void HMatrixFactory::save(Advocate & adv) const
135 {
136 PersistentObject::save(adv);
137 }
138
139
140 /* Method load() reloads the object from the StorageManager */
load(Advocate & adv)141 void HMatrixFactory::load(Advocate & adv)
142 {
143 PersistentObject::load(adv);
144 }
145
146
147 END_NAMESPACE_OPENTURNS
148