Line data Source code
1 : //----------------------------------------------------------------------
2 : // This file is part of PARSEC (http://physik.rwth-aachen.de/parsec)
3 : // a parametrized simulation engine for cosmic rays.
4 : //
5 : // Copyright (C) 2011 Martin Erdmann, Peter Schiffer, Tobias Winchen
6 : // RWTH Aachen University, Germany
7 : // Contact: winchen@physik.rwth-aachen.de
8 : //
9 : // This program is free software: you can redistribute it and/or
10 : // modify it under the terms of the GNU General Public License as
11 : // published by the Free Software Foundation, either version 3 of
12 : // the License, or (at your option) any later version.
13 : //
14 : // This program is distributed in the hope that it will be useful,
15 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : // GNU General Public License for more details.
18 : //
19 : // You should have received a copy of the GNU General Public License
20 : // along with this program. If not, see <http://www.gnu.org/licenses/>.
21 : //----------------------------------------------------------------------
22 :
23 : #include "crpropa/magneticLens/ModelMatrix.h"
24 : #include <ctime>
25 :
26 : #include <Eigen/Core>
27 : namespace crpropa
28 : {
29 :
30 0 : void serialize(const string &filename, const ModelMatrixType& matrix)
31 : {
32 0 : ofstream outfile(filename.c_str(), ios::binary);
33 0 : if (!outfile)
34 : {
35 0 : throw runtime_error("Can't write file: " + filename);
36 : }
37 :
38 0 : uint32_t C = 0;
39 : double val;
40 :
41 0 : C = (uint32_t) (matrix.nonZeros());
42 0 : outfile.write((char*) &C, sizeof(uint32_t));
43 0 : C = (uint32_t) (matrix.rows());
44 0 : outfile.write((char*) &C, sizeof(uint32_t));
45 0 : C = (uint32_t) (matrix.cols());
46 0 : outfile.write((char*) &C, sizeof(uint32_t));
47 :
48 : // serialize non zero elements
49 0 : for (size_t col_idx = 0; col_idx < matrix.cols(); ++col_idx)
50 : {
51 0 : for (ModelMatrixType::InnerIterator it(matrix,col_idx); it; ++it)
52 : {
53 : it.value();
54 0 : C = (uint32_t) it.row();
55 0 : outfile.write((char*) &C, sizeof(uint32_t));
56 :
57 0 : C = (uint32_t) it.col();
58 0 : outfile.write((char*) &C, sizeof(uint32_t));
59 :
60 0 : val = it.value();
61 0 : outfile.write((char*) &val, sizeof(double));
62 0 : if (outfile.fail())
63 : {
64 0 : throw runtime_error("Error writing file: " + filename);
65 : }
66 : }
67 : }
68 :
69 0 : outfile.close();
70 0 : }
71 :
72 :
73 0 : void deserialize(const string &filename, ModelMatrixType& matrix)
74 : {
75 0 : ifstream infile(filename.c_str(), ios::binary);
76 0 : if (!infile)
77 : {
78 0 : throw runtime_error("Can't read file: " + filename);
79 : }
80 :
81 : uint32_t nnz, nRows, nColumns;
82 0 : infile.read((char*) &nnz, sizeof(uint32_t));
83 0 : infile.read((char*) &nRows, sizeof(uint32_t));
84 0 : infile.read((char*) &nColumns, sizeof(uint32_t));
85 0 : matrix.resize(nRows, nColumns);
86 0 : matrix.reserve(nnz);
87 :
88 : uint32_t row, column;
89 : double val;
90 : std::vector< Eigen::Triplet<double> > triplets;
91 0 : triplets.resize(nnz);
92 0 : for (size_t i = 0; i < nnz; i++)
93 : {
94 0 : infile.read((char*) &row, sizeof(uint32_t));
95 0 : infile.read((char*) &column, sizeof(uint32_t));
96 0 : infile.read((char*) &val, sizeof(double));
97 : //M(size1,size2) = val;
98 0 : triplets[i] = Eigen::Triplet<double>(row, column, val);
99 : }
100 0 : matrix.setFromTriplets(triplets.begin(), triplets.end());
101 0 : matrix.makeCompressed();
102 0 : }
103 :
104 :
105 0 : double norm_1(const ModelVectorType &v)
106 : {
107 0 : return v.cwiseAbs().sum();
108 : }
109 :
110 :
111 0 : void normalizeColumns(ModelMatrixType &matrix){
112 0 : for (size_t i=0; i< matrix.cols(); i++)
113 : {
114 0 : ModelVectorType v = matrix.col(i);
115 0 : double rn = norm_1(v);
116 0 : matrix.col(i) = v/rn;
117 : }
118 0 : }
119 :
120 :
121 0 : double maximumOfSumsOfColumns(const ModelMatrixType &matrix)
122 : {
123 : double summax = 0;
124 0 : for (size_t i = 0; i < matrix.cols(); i++)
125 : {
126 0 : double sum = matrix.col(i).sum();
127 0 : if (sum > summax)
128 : summax = sum;
129 : }
130 0 : return summax;
131 : }
132 :
133 0 : void normalizeMatrix(ModelMatrixType& matrix, double norm)
134 : {
135 0 : matrix /= norm;
136 0 : }
137 :
138 0 : void prod_up(const ModelMatrixType& matrix, double* model)
139 : {
140 :
141 : // copy storage of model, as matrix vector product cannot be done
142 : // in place
143 :
144 0 : const size_t mSize = matrix.cols();
145 0 : double *origVectorStorage = new double[mSize];
146 : memcpy(origVectorStorage, model, mSize * sizeof(double));
147 :
148 :
149 :
150 : Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, 1> > origVectorAdaptor(origVectorStorage, mSize);
151 : Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, 1> > modelVectorAdaptor(model, mSize);
152 :
153 : // perform the optimized product
154 0 : modelVectorAdaptor = matrix * origVectorAdaptor;
155 :
156 : // clean up
157 0 : delete[] origVectorStorage;
158 0 : }
159 :
160 :
161 : } // namespace parsec
|
