1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
|
/***********************************************/
/**
* @file normalsScale.cpp
*
* @brief Scales rows and columns of a normal equation system.
*
* @author Torsten Mayer-Guerr
* @date 2011-02-23
*/
/***********************************************/
// Latex documentation
#define DOCSTRING docstring
static const char *docstring = R"(
Scales rows and columns of a system of \configFile{inputfileNormalEquation}{normalEquation}
given by a diagonal matrix \configFile{inputfileFactorVector}{matrix} $\M S$
\begin{equation}
\bar{\M N} := \M S \M N \M S \qquad\text{and}\qquad \bar{\M n} := \M S \M n.
\end{equation}
The estimated solution is now
\begin{equation}
\bar{\M x} := \M S^{-1} \M x.
\end{equation}
This is effectively the same as rescaling columns of the design matrix.
This program is useful when combining normal equations from different sources,
for example in case the units of certain parameters don't match.
)";
/***********************************************/
#include "programs/program.h"
#include "parallel/matrixDistributed.h"
#include "files/fileMatrix.h"
#include "files/fileNormalEquation.h"
/***** CLASS ***********************************/
/** @brief Scales rows and columns of a normal equation system.
* @ingroup programsGroup */
class NormalsScale
{
public:
void run(Config &config, Parallel::CommunicatorPtr comm);
};
GROOPS_REGISTER_PROGRAM(NormalsScale, PARALLEL, "Scales rows and columns of a normal equation system", NormalEquation)
/***********************************************/
void NormalsScale::run(Config &config, Parallel::CommunicatorPtr comm)
{
try
{
FileName outName, normalsName, factorName;
renameDeprecatedConfig(config, "outputfileNormalequation", "outputfileNormalEquation", date2time(2020, 6, 3));
renameDeprecatedConfig(config, "inputfileNormalequation", "inputfileNormalEquation", date2time(2020, 6, 3));
readConfig(config, "outputfileNormalEquation", outName, Config::MUSTSET, "", "");
readConfig(config, "inputfileNormalEquation", normalsName, Config::MUSTSET, "", "");
readConfig(config, "inputfileFactorVector", factorName, Config::MUSTSET, "", "Vector containing the factors");
if(isCreateSchema(config)) return;
// ==================================
logStatus<<"read normal equations <"<<outName<<">"<<Log::endl;
MatrixDistributed normal;
Matrix n;
NormalEquationInfo info;
readFileNormalEquation(normalsName, info, normal, n, comm);
logInfo<<" number of parameters: "<<normal.parameterCount()<<Log::endl;
logInfo<<" number of right hand sides: "<<info.lPl.rows()<<Log::endl;
logStatus<<"read factor vector <"<<factorName<<">"<<Log::endl;
Vector factor;
readFileMatrix(factorName, factor);
if(Parallel::isMaster(comm) && (factor.rows() != n.rows()))
throw(Exception("Dimension error factor("+factor.rows()%"%i) != n("s+n.rows()%"%i)"s));
// ==================================
logStatus<<"scale normal matrix"<<Log::endl;
for(UInt i=0; i<normal.blockCount(); i++)
for(UInt k=i; k<normal.blockCount(); k++)
if(normal.isMyRank(i,k))
{
Matrix &N = normal.N(i,k);
for(UInt z=0; z<N.rows(); z++)
N.row(z) *= factor(z+normal.blockIndex(i));
for(UInt s=0; s<N.columns(); s++)
N.column(s) *= factor(s+normal.blockIndex(k));
}
if(Parallel::isMaster(comm))
for(UInt z=0; z<n.rows(); z++)
n.row(z) *= factor(z)*factor(z);
// ==================================
// write normal equations
// ----------------------
logStatus<<"write normal equations to <"<<outName<<">"<<Log::endl;
writeFileNormalEquation(outName, info, normal, n);
}
catch(std::exception &e)
{
GROOPS_RETHROW(e)
}
}
/***********************************************/
|
