From f4385fc3fb1f732a556e0ee6470e21cab6a30c1e Mon Sep 17 00:00:00 2001 From: Ross Brodie Date: Wed, 11 Sep 2024 10:17:26 +1000 Subject: [PATCH] Catered for deprecation of PETSC_NULL from Petsc v3.19 --- src/gaaem_version.h | 2 +- src/galeiallatonce.cpp | 1164 ++++++++--------- src/petsc_wrapper.h | 544 ++++---- .../galeiallatonce/galeiallatonce.vcxproj | 2 +- .../galeiallatonce.vcxproj.filters | 6 +- 5 files changed, 862 insertions(+), 856 deletions(-) diff --git a/src/gaaem_version.h b/src/gaaem_version.h index f95e9b6..fc883db 100644 --- a/src/gaaem_version.h +++ b/src/gaaem_version.h @@ -9,6 +9,6 @@ Author: Ross C. Brodie, Geoscience Australia. #ifndef _gaaem_version_H #define _gaaem_version_H -constexpr auto GAAEM_VERSION = "v2.0.1-Release-20240619"; +constexpr auto GAAEM_VERSION = "v2.0.2-Release-20240811"; #endif diff --git a/src/galeiallatonce.cpp b/src/galeiallatonce.cpp index 4bc553c..a863707 100644 --- a/src/galeiallatonce.cpp +++ b/src/galeiallatonce.cpp @@ -40,50 +40,50 @@ class cTDEmComponentInfo; enum eSmoothnessMethod { SM_1ST_DERIVATIVE, SM_2ND_DERIVATIVE }; -class cLineTiles{ +class cLineTiles { std::vector> L; - + public: - cLineTiles(const double linelen, const double seglen, const double overlap){ + cLineTiles(const double linelen, const double seglen, const double overlap) { calculate_tiles(linelen, seglen, overlap); }; - void calculate_tiles(const double linelen, const double seglen, const double overlap){ - int n = (int)std::round((linelen - overlap) / (seglen - overlap)); - double newseglen = (linelen + 2.0*overlap*(n-1)) / (double)n; + void calculate_tiles(const double linelen, const double seglen, const double overlap) { + int n = (int)std::round((linelen - overlap) / (seglen - overlap)); + double newseglen = (linelen + 2.0 * overlap * (n - 1)) / (double)n; L.resize(n); - for (size_t i = 0; i < (size_t)n; i++){ + for (size_t i = 0; i < (size_t)n; i++) { L[i].resize(4); - if (i == 0){ + if (i == 0) { L[i][0] = 0.0; } - else{ - L[i][0] = L[i-1][0] + newseglen - 2*overlap; + else { + L[i][0] = L[i - 1][0] + newseglen - 2 * overlap; } - + L[i][1] = L[i][0] + overlap; L[i][2] = L[i][0] + newseglen - overlap; L[i][3] = L[i][0] + newseglen; - if (i == 0){ + if (i == 0) { L[i][1] = L[i][0]; } - if (i == (size_t)(n-1)){ + if (i == (size_t)(n - 1)) { L[i][2] = L[i][3]; } } std::cout << *this; }; - + friend std::ostream& operator<<(std::ostream& os, const cLineTiles& t); }; -std::ostream& operator<<(std::ostream& os, const cLineTiles& t){ - for (size_t i = 0; i < t.L.size(); i++){ +std::ostream& operator<<(std::ostream& os, const cLineTiles& t) { + for (size_t i = 0; i < t.L.size(); i++) { os << t.L[i][0] << "\t"; os << t.L[i][1] << "\t"; os << t.L[i][2] << "\t"; @@ -92,7 +92,7 @@ std::ostream& operator<<(std::ostream& os, const cLineTiles& t){ return os; }; -class cField{ +class cField { std::vector _data; size_t _ns; @@ -108,21 +108,21 @@ class cField{ std::vector defaultvec; - double& el(const size_t& i, const size_t& j){ - return _data[i*_nb + j]; + double& el(const size_t& i, const size_t& j) { + return _data[i * _nb + j]; }; - size_t cind(){ return cn - 1; } + size_t cind() { return cn - 1; } public: - cField(){ + cField() { defined = false; cndef = false; ncndef = false; nvals = 0; }; - cField(const cBlock& b, const std::string _id, size_t _size = 1){ + cField(const cBlock& b, const std::string _id, size_t _size = 1) { defined = false; cndef = false; ncndef = false; @@ -130,28 +130,28 @@ class cField{ id = _id; def = b.getstringvalue(_id); - if (!isdefined(def)){ + if (!isdefined(def)) { nvals = 0; defined = false; return; } std::vector t = tokenize(def); - if (strcasecmp(t[0], "column") == 0){ + if (strcasecmp(t[0], "column") == 0) { cndef = true; cn = (size_t)atoi(t[1].c_str()); } - else if (strcasecmp(t[0], "-column") == 0){ + else if (strcasecmp(t[0], "-column") == 0) { ncndef = true; cn = (size_t)atoi(t[1].c_str()); } - else{ - for (size_t i = 0; i < t.size(); i++){ + else { + for (size_t i = 0; i < t.size(); i++) { defaultvec.push_back(atof(t[i].c_str())); } - if (defaultvec.size() == 1){ + if (defaultvec.size() == 1) { defaultvec.resize(nvals); - for (size_t i = 1; i < defaultvec.size(); i++){ + for (size_t i = 1; i < defaultvec.size(); i++) { defaultvec[i] = defaultvec[0]; } } @@ -159,29 +159,29 @@ class cField{ defined = true; } - void resize(const size_t ns, size_t nb = 1){ + void resize(const size_t ns, size_t nb = 1) { _ns = ns; _nb = nb; - _data.resize(ns*nb); + _data.resize(ns * nb); }; - bool parse(const std::vector& tokens, const size_t si){ - if (nvals == 0){ + bool parse(const std::vector& tokens, const size_t si) { + if (nvals == 0) { return false; } - else if (nvals == 1){ + else if (nvals == 1) { el(si, 0) = get(tokens); } - else{ + else { std::vector v = getmulti(tokens); - for (size_t bi = 0; bi < nvals; bi++){ + for (size_t bi = 0; bi < nvals; bi++) { el(si, bi) = v[bi]; } } return true; } - double get(const std::vector& tokens){ + double get(const std::vector& tokens) { double v; if (cndef) v = atof(tokens[cind()].c_str()); else if (ncndef) v = -atof(tokens[cind()].c_str()); @@ -189,27 +189,27 @@ class cField{ return v; } - std::vector getmulti(const std::vector& tokens){ + std::vector getmulti(const std::vector& tokens) { std::vector v(nvals); - if (cndef){ - for (size_t k = 0; k < nvals; k++){ + if (cndef) { + for (size_t k = 0; k < nvals; k++) { v[k] = atof(tokens[cind() + k].c_str()); } } - else if (ncndef){ - for (size_t k = 0; k < nvals; k++){ + else if (ncndef) { + for (size_t k = 0; k < nvals; k++) { v[k] = -atof(tokens[cind() + k].c_str()); } } - else{ - for (size_t k = 0; k < nvals; k++){ + else { + for (size_t k = 0; k < nvals; k++) { v[k] = defaultvec[k]; } } return v; } - double& operator()(const size_t& i, const size_t& j = 0){ + double& operator()(const size_t& i, const size_t& j = 0) { return el(i, j); }; }; @@ -217,7 +217,7 @@ class cField{ class cTDEmComponentInfo { public: - + bool Use = false; bool InvertTotalField = false; bool EstimateNoiseFromModel = false; @@ -231,11 +231,11 @@ class cTDEmComponentInfo { cTDEmComponentInfo() {}; - cTDEmComponentInfo(const cBlock& b, size_t nwindows, bool inverttotalfield) + cTDEmComponentInfo(const cBlock& b, size_t nwindows, bool inverttotalfield) { InvertTotalField = inverttotalfield; nw = nwindows; - if (b.Entries.size() == 0){ + if (b.Entries.size() == 0) { Use = false; return; } @@ -248,22 +248,22 @@ class cTDEmComponentInfo { EstimateNoiseFromModel = b.getboolvalue("EstimateNoiseFromModel"); - if (EstimateNoiseFromModel){ + if (EstimateNoiseFromModel) { fdmn = cField(b, "MultiplicativeNoise", nwindows); fdan = cField(b, "AdditiveNoise", nwindows); } - else{ + else { fdn = cField(b, "Noise", nwindows); } } - size_t ndata(){ + size_t ndata() { if (Use)return nw; else return 0; }; - bool allocate_data_arrays(const size_t nlocalsamples){ + bool allocate_data_arrays(const size_t nlocalsamples) { if (Use == false)return true; fds.resize(nlocalsamples, nw); fdp.resize(nlocalsamples, 1); @@ -273,7 +273,7 @@ class cTDEmComponentInfo { return true; }; - bool parse(const std::vector tokens, const size_t localsampleindex){ + bool parse(const std::vector tokens, const size_t localsampleindex) { if (Use == false)return true; fds.parse(tokens, localsampleindex); fdp.parse(tokens, localsampleindex); @@ -283,35 +283,35 @@ class cTDEmComponentInfo { return true; }; - std::vector data(const size_t localsampleindex){ + std::vector data(const size_t localsampleindex) { std::vector v; if (Use == false)return v; v.resize(nw); - for (size_t wi = 0; wi < nw; wi++){ + for (size_t wi = 0; wi < nw; wi++) { v[wi] = fds(localsampleindex, wi); - if (InvertTotalField){ + if (InvertTotalField) { v[wi] += fdp(localsampleindex); } } return v; }; - std::vector noise(const size_t localsampleindex){ + std::vector noise(const size_t localsampleindex) { std::vector v; if (Use == false)return v; v.resize(nw); - if (EstimateNoiseFromModel){ - for (size_t wi = 0; wi < nw; wi++){ + if (EstimateNoiseFromModel) { + for (size_t wi = 0; wi < nw; wi++) { double an = fdan(localsampleindex, wi); double pmn = fdmn(localsampleindex, wi); double s = fds(localsampleindex, wi); - double mn = 0.01*pmn*s; - v[wi] = sqrt(an*an + mn*mn); + double mn = 0.01 * pmn * s; + v[wi] = sqrt(an * an + mn * mn); } } - else{ - for (size_t wi = 0; wi < nw; wi++){ + else { + for (size_t wi = 0; wi < nw; wi++) { v[wi] = fdn(localsampleindex, wi); } } @@ -320,24 +320,24 @@ class cTDEmComponentInfo { }; -class cSystemInfo{ +class cSystemInfo { private: size_t nw; - + public: cTDEmSystem T; - bool InvertTotalField; + bool InvertTotalField; std::vector Comp; - cSystemInfo(){ }; - bool initialise(const cBlock& b){ + cSystemInfo() { }; + bool initialise(const cBlock& b) { std::string stm = b.getstringvalue("SystemFile"); T.readsystemdescriptorfile(stm); nw = T.NumberOfWindows; bool status; status = b.getvalue("InvertTotalField", InvertTotalField); - if (status == false){ + if (status == false) { InvertTotalField = false; } @@ -345,59 +345,59 @@ class cSystemInfo{ Comp[0] = cTDEmComponentInfo(b.findblock("XComponent"), nw, InvertTotalField); Comp[1] = cTDEmComponentInfo(b.findblock("YComponent"), nw, InvertTotalField); Comp[2] = cTDEmComponentInfo(b.findblock("ZComponent"), nw, InvertTotalField); - + Comp[0].basedindex = 0; Comp[1].basedindex = 0; Comp[2].basedindex = 0; - if (Comp[0].Use) Comp[1].basedindex += nw; + if (Comp[0].Use) Comp[1].basedindex += nw; if (Comp[0].Use) Comp[2].basedindex += nw; - if (Comp[1].Use) Comp[2].basedindex += nw; + if (Comp[1].Use) Comp[2].basedindex += nw; return true; }; - size_t ndata(){ + size_t ndata() { return Comp[0].ndata() + Comp[1].ndata() + Comp[2].ndata(); } - inline size_t dindex(const size_t& component, const size_t& window){ + inline size_t dindex(const size_t& component, const size_t& window) { return Comp[component].basedindex + window; } - bool allocate_data_arrays(const size_t nlocalsamples){ - for (size_t k = 0; k < Comp.size(); k++){ + bool allocate_data_arrays(const size_t nlocalsamples) { + for (size_t k = 0; k < Comp.size(); k++) { Comp[k].allocate_data_arrays(nlocalsamples); } return true; } - bool parse(const std::vector tokens, const size_t localsampleindex){ - for (size_t ci = 0; ci < Comp.size(); ci++){ + bool parse(const std::vector tokens, const size_t localsampleindex) { + for (size_t ci = 0; ci < Comp.size(); ci++) { Comp[ci].parse(tokens, localsampleindex); } return true; } - std::vector data(const size_t localsampleindex){ - std::vector v; - for (size_t ci = 0; ci < Comp.size(); ci++){ - if (Comp[ci].Use){ + std::vector data(const size_t localsampleindex) { + std::vector v; + for (size_t ci = 0; ci < Comp.size(); ci++) { + if (Comp[ci].Use) { append(v, Comp[ci].data(localsampleindex)); } } return v; } - std::vector noise(const size_t localsampleindex){ - std::vector v; - for (size_t ci = 0; ci < Comp.size(); ci++){ - if (Comp[ci].Use){ + std::vector noise(const size_t localsampleindex) { + std::vector v; + for (size_t ci = 0; ci < Comp.size(); ci++) { + if (Comp[ci].Use) { append(v, Comp[ci].noise(localsampleindex)); } } return v; } - bool forward_model(const std::vector& conductivity, const std::vector& thickness, const cTDEmGeometry& geometry){ + bool forward_model(const std::vector& conductivity, const std::vector& thickness, const cTDEmGeometry& geometry) { T.setconductivitythickness(conductivity, thickness); T.setgeometry(geometry); T.LEM.calculation_type = cLEM::CalculationType::FORWARDMODEL; @@ -408,7 +408,7 @@ class cSystemInfo{ return true; } - bool forward_model_and_derivatives(const std::vector& conductivity, const std::vector& thickness, const cTDEmGeometry& geometry, std::vector& predicted, std::vector>& derivatives, const bool computederivatives, const std::vector UGI){ + bool forward_model_and_derivatives(const std::vector& conductivity, const std::vector& thickness, const cTDEmGeometry& geometry, std::vector& predicted, std::vector>& derivatives, const bool computederivatives, const std::vector UGI) { size_t nlayers = conductivity.size(); T.setconductivitythickness(conductivity, thickness); T.setgeometry(geometry); @@ -418,90 +418,90 @@ class cSystemInfo{ T.setprimaryfields(); T.setsecondaryfields(); - + //Save for later derivative calculations std::vector X = T.X; std::vector Y = T.Y; std::vector Z = T.Z; - if (InvertTotalField){ + if (InvertTotalField) { X += T.PrimaryX; Y += T.PrimaryY; Z += T.PrimaryZ; } - + predicted.resize(ndata()); - for (size_t ci = 0; ci < Comp.size(); ci++){ + for (size_t ci = 0; ci < Comp.size(); ci++) { if (Comp[ci].Use == false)continue; - for (size_t wi = 0; wi < T.NumberOfWindows; wi++){ - predicted[dindex(ci,wi)] = T.secondary(ci, wi); - if (Comp[ci].InvertTotalField){ + for (size_t wi = 0; wi < T.NumberOfWindows; wi++) { + predicted[dindex(ci, wi)] = T.secondary(ci, wi); + if (Comp[ci].InvertTotalField) { predicted[dindex(ci, wi)] += T.primary(ci); - } + } } } - if (computederivatives == true){ - + if (computederivatives == true) { + derivatives.resize(ndata()); - for (size_t di = 0; di < ndata(); di++){ + for (size_t di = 0; di < ndata(); di++) { derivatives[di].resize(nlayers + UGI.size()); } - for (size_t li = 0; li < nlayers; li++){ + for (size_t li = 0; li < nlayers; li++) { T.LEM.calculation_type = cLEM::CalculationType::CONDUCTIVITYDERIVATIVE; T.LEM.derivative_layer = li; T.setupcomputations(); T.setprimaryfields(); T.setsecondaryfields(); - - for (size_t ci = 0; ci < Comp.size(); ci++){ + + for (size_t ci = 0; ci < Comp.size(); ci++) { if (Comp[ci].Use == false)continue; - for (size_t wi = 0; wi < T.NumberOfWindows; wi++){ + for (size_t wi = 0; wi < T.NumberOfWindows; wi++) { derivatives[dindex(ci, wi)][li] = T.secondary(ci, wi); - if (Comp[ci].InvertTotalField){ + if (Comp[ci].InvertTotalField) { derivatives[dindex(ci, wi)][li] += T.primary(ci); - } + } } } - } + } - for (size_t gi = 0; gi < UGI.size(); gi++){ - if (cTDEmGeometry::elementtype(UGI[gi]) == cTDEmGeometry::ElementType::rx_pitch){ + for (size_t gi = 0; gi < UGI.size(); gi++) { + if (cTDEmGeometry::elementtype(UGI[gi]) == cTDEmGeometry::ElementType::rx_pitch) { std::vector dxbdp; std::vector dzbdp; - T.drx_pitch(X, Z, geometry.rx_pitch, dxbdp, dzbdp); - for (size_t ci = 0; ci < Comp.size(); ci++){ + T.drx_pitch(X, Z, geometry.rx_pitch, dxbdp, dzbdp); + for (size_t ci = 0; ci < Comp.size(); ci++) { if (Comp[ci].Use == false)continue; - for (size_t wi = 0; wi < T.NumberOfWindows; wi++){ + for (size_t wi = 0; wi < T.NumberOfWindows; wi++) { if (ci == 0) derivatives[dindex(ci, wi)][gi + nlayers] = dxbdp[wi]; else if (ci == 1) derivatives[dindex(ci, wi)][gi + nlayers] = 0.0; - else derivatives[dindex(ci, wi)][gi + nlayers] = dzbdp[wi]; + else derivatives[dindex(ci, wi)][gi + nlayers] = dzbdp[wi]; } } } - else{ + else { T.LEM.calculation_type = cTDEmGeometry::derivativetype(UGI[gi]); T.LEM.derivative_layer = INT_MAX; T.setupcomputations(); T.setprimaryfields(); - T.setsecondaryfields(); - for (size_t ci = 0; ci < Comp.size(); ci++){ + T.setsecondaryfields(); + for (size_t ci = 0; ci < Comp.size(); ci++) { if (Comp[ci].Use == false) continue; - for (size_t wi = 0; wi < T.NumberOfWindows; wi++){ + for (size_t wi = 0; wi < T.NumberOfWindows; wi++) { derivatives[dindex(ci, wi)][gi + nlayers] = T.secondary(ci, wi); - if (Comp[ci].InvertTotalField){ + if (Comp[ci].InvertTotalField) { derivatives[dindex(ci, wi)][gi + nlayers] += T.primary(ci); - } + } } } - } + } } } return true; } }; -class cEarthInfo{ +class cEarthInfo { private: size_t nlayers; @@ -512,8 +512,8 @@ class cEarthInfo{ cField cstd; cField tstd; - cEarthInfo(){}; - cEarthInfo(const cBlock& b){ + cEarthInfo() {}; + cEarthInfo(const cBlock& b) { nlayers = b.getsizetvalue("NumberOfLayers"); cref = cField(b, "ReferenceModel.Conductivity", nlayers); tref = cField(b, "ReferenceModel.Thickness", nlayers - 1); @@ -522,9 +522,9 @@ class cEarthInfo{ tstd = cField(b, "StdDevReferenceModel.Thickness", nlayers - 1); } - size_t numlayers(){ return nlayers; } + size_t numlayers() { return nlayers; } - bool allocate_data_arrays(const size_t nlocalsamples){ + bool allocate_data_arrays(const size_t nlocalsamples) { cref.resize(nlocalsamples, nlayers); tref.resize(nlocalsamples, nlayers - 1); cstd.resize(nlocalsamples, nlayers); @@ -532,7 +532,7 @@ class cEarthInfo{ return true; } - bool parse(const std::vector tokens, const size_t localsampleindex){ + bool parse(const std::vector tokens, const size_t localsampleindex) { cref.parse(tokens, localsampleindex); tref.parse(tokens, localsampleindex); cstd.parse(tokens, localsampleindex); @@ -541,7 +541,7 @@ class cEarthInfo{ } }; -class cInversionOptions{ +class cInversionOptions { public: double CorrelationRadius; @@ -555,9 +555,9 @@ class cInversionOptions{ double AlphaR = 1; eSmoothnessMethod VerticalSmoothnessMethod = SM_2ND_DERIVATIVE; - cInversionOptions(){}; + cInversionOptions() {}; - cInversionOptions(const cBlock& b){ + cInversionOptions(const cBlock& b) { CorrelationRadius = b.getdoublevalue("CorrelationRadius"); InverseDistancePower = b.getdoublevalue("InverseDistancePower"); @@ -570,30 +570,30 @@ class cInversionOptions{ AlphaR = b.getdoublevalue("AlphaReferenceModel"); std::string sm = b.getstringvalue("VerticalSmoothnessMethod"); - if (!isdefined(sm)){ + if (!isdefined(sm)) { VerticalSmoothnessMethod = SM_2ND_DERIVATIVE; } - else if (strcasecmp(sm, "Minimise1stDerivatives") == 0){ + else if (strcasecmp(sm, "Minimise1stDerivatives") == 0) { VerticalSmoothnessMethod = SM_1ST_DERIVATIVE; } - else if (strcasecmp(sm, "Minimize1stDerivatives") == 0){ + else if (strcasecmp(sm, "Minimize1stDerivatives") == 0) { VerticalSmoothnessMethod = SM_1ST_DERIVATIVE; } - else if (strcasecmp(sm, "Minimise2ndDerivatives") == 0){ + else if (strcasecmp(sm, "Minimise2ndDerivatives") == 0) { VerticalSmoothnessMethod = SM_2ND_DERIVATIVE; } - else if (strcasecmp(sm, "Minimize2ndDerivatives") == 0){ + else if (strcasecmp(sm, "Minimize2ndDerivatives") == 0) { VerticalSmoothnessMethod = SM_2ND_DERIVATIVE; } - else{ - glog.logmsg(0,"Unknown SmoothnessMethod %s\n", sm.c_str()); + else { + glog.logmsg(0, "Unknown SmoothnessMethod %s\n", sm.c_str()); std::string e = strprint("Error: exception thrown from %s (%d) %s\n", __FILE__, __LINE__, __FUNCTION__); throw e; } } }; -class cInputOptions{ +class cInputOptions { public: std::string DataFile; @@ -604,63 +604,63 @@ class cInputOptions{ std::vector IncludeLines; std::vector> IncludeLineRanges; - cInputOptions(){}; + cInputOptions() {}; - cInputOptions(const cBlock& b){ + cInputOptions(const cBlock& b) { bool status; - + status = b.getvalue("HeaderLines", HeaderLines); - if (status == false){ + if (status == false) { HeaderLines = 0; } status = b.getvalue("Subsample", SubSample); - if (status == false){ + if (status == false) { SubSample = 1; } - if (b.getvalue("DataFile", DataFile) == false){ - glog.logmsg(0, "Input DataFile was not specified\n"); + if (b.getvalue("DataFile", DataFile) == false) { + glog.logmsg(0, "Input DataFile was not specified\n"); std::string e = strprint("Error: exception thrown from %s (%d) %s\n", __FILE__, __LINE__, __FUNCTION__); throw e; } - else if (exists(DataFile) == false){ - glog.logmsg(0, "Input DataFile %s not found\n", DataFile.c_str()); + else if (exists(DataFile) == false) { + glog.logmsg(0, "Input DataFile %s not found\n", DataFile.c_str()); std::string e = strprint("Error: exception thrown from %s (%d) %s\n", __FILE__, __LINE__, __FUNCTION__); throw e; } std::string pfile; - if (b.getvalue("IncludePolygon", pfile)){ + if (b.getvalue("IncludePolygon", pfile)) { IncludePolygons.push_back(cPolygon(pfile)); } std::string s; - if (b.getvalue("IncludeLines", s)){ + if (b.getvalue("IncludeLines", s)) { parse_include_lines(s); } } - void parse_include_lines(const std::string& s){ + void parse_include_lines(const std::string& s) { std::pair r; std::vector t = tokenize(s); - for (size_t i = 0; i < t.size(); i++){ + for (size_t i = 0; i < t.size(); i++) { - if (t[i][0] == ':' || t[i][t[i].size() - 1] == ':'){ - glog.logmsg(0, "Error bad token (%s) when parsing Input.IncludeLines\n", t[i].c_str()); + if (t[i][0] == ':' || t[i][t[i].size() - 1] == ':') { + glog.logmsg(0, "Error bad token (%s) when parsing Input.IncludeLines\n", t[i].c_str()); std::string e = strprint("Error: exception thrown from %s (%d) %s\n", __FILE__, __LINE__, __FUNCTION__); throw e; } - else if (std::sscanf(t[i].c_str(), "%d:%d", &r.first, &r.second) == 2){ + else if (std::sscanf(t[i].c_str(), "%d:%d", &r.first, &r.second) == 2) { IncludeLineRanges.push_back(r); } - else if (std::sscanf(t[i].c_str(), "%d", &r.first) == 1){ + else if (std::sscanf(t[i].c_str(), "%d", &r.first) == 1) { IncludeLines.push_back(r.first); } - else{ - glog.logmsg(0, "Error bad token (%s) when parsing Input.IncludeLines\n", t[i].c_str()); + else { + glog.logmsg(0, "Error bad token (%s) when parsing Input.IncludeLines\n", t[i].c_str()); std::string e = strprint("Error: exception thrown from %s (%d) %s\n", __FILE__, __LINE__, __FUNCTION__); throw e; } @@ -668,20 +668,20 @@ class cInputOptions{ } - bool is_line_included(const int& line){ + bool is_line_included(const int& line) { - if (IncludeLines.size() == 0 && IncludeLineRanges.size() == 0) return true; + if (IncludeLines.size() == 0 && IncludeLineRanges.size() == 0) return true; auto rit = std::find_if( IncludeLineRanges.begin(), IncludeLineRanges.end(), [&line](const std::pair& r) - { - return (line >= r.first && line <= r.second); - } + { + return (line >= r.first && line <= r.second); + } ); if (rit != IncludeLineRanges.end())return true; - + auto lit = std::find(IncludeLines.begin(), IncludeLines.end(), line); if (lit != IncludeLines.end())return true; @@ -691,7 +691,7 @@ class cInputOptions{ bool is_point_included(const cPoint& p) { bool insidestatus = true; - for (size_t pi = 0; pi < IncludePolygons.size(); pi++){ + for (size_t pi = 0; pi < IncludePolygons.size(); pi++) { insidestatus = IncludePolygons[pi].isinside(p); if (insidestatus) break; } @@ -700,29 +700,29 @@ class cInputOptions{ }; -class cOutputOptions{ +class cOutputOptions { public: std::string LogFile; std::string DataFile; bool PredictedData = true; - bool ObservedData = true; - bool Noise = true; + bool ObservedData = true; + bool Noise = true; bool PositiveLayerBottomDepths = false; bool NegativeLayerBottomDepths = false; bool InterfaceElevations = false; - cOutputOptions(){}; - cOutputOptions(const cBlock& b){ + cOutputOptions() {}; + cOutputOptions(const cBlock& b) { - if (b.getvalue("LogFile", LogFile) == false){ - glog.logmsg(0, "Output LogFile was not specified\n"); + if (b.getvalue("LogFile", LogFile) == false) { + glog.logmsg(0, "Output LogFile was not specified\n"); std::string e = strprint("Error: exception thrown from %s (%d) %s\n", __FILE__, __LINE__, __FUNCTION__); throw e; } - if (b.getvalue("DataFile", DataFile) == false){ - glog.logmsg(0, "Output DataFile was not specified\n"); + if (b.getvalue("DataFile", DataFile) == false) { + glog.logmsg(0, "Output DataFile was not specified\n"); std::string e = strprint("Error: exception thrown from %s (%d) %s\n", __FILE__, __LINE__, __FUNCTION__); throw e; } @@ -743,24 +743,24 @@ class cGeometryInfo { bool solve = false; cField ref; cField std; - - cGeometryInfo(){ }; - cGeometryInfo(const cBlock& b){ - ref = cField(b, "Reference"); + cGeometryInfo() { }; + + cGeometryInfo(const cBlock& b) { + ref = cField(b, "Reference"); b.getvalue("Solve", solve); - if (solve){ + if (solve) { std = cField(b, "Uncertainty"); } } }; -class cAllAtOnceInverter{ +class cAllAtOnceInverter { private: cBlock Control; - + cMpiEnv mpienv; cMpiComm mpicomm; cRadiusSearcher RS; @@ -772,7 +772,7 @@ class cAllAtOnceInverter{ cField fdfiducial; cField fdx; cField fdy; - cField fdelevation; + cField fdelevation; std::vector G; std::vector UGI; @@ -793,7 +793,7 @@ class cAllAtOnceInverter{ size_t nparampersample; size_t nparam; - + size_t nlocaldata; size_t nlocalparam; cOwnership sown; @@ -827,14 +827,14 @@ class cAllAtOnceInverter{ cAllAtOnceInverter(int argc, char** argv) { - if(argc < 2){ + if (argc < 2) { glog.logmsg("%s\n", commandlinestring(argc, argv).c_str()); glog.logmsg("%s\n", versionstring(GAAEM_VERSION, __TIME__, __DATE__).c_str()); glog.logmsg("Usage: %s control_file_name\n", argv[0]); glog.logmsg("Too few command line arguments\n"); exit(1); } - else if(argc > 2){ + else if (argc > 2) { glog.logmsg("%s\n", commandlinestring(argc, argv).c_str()); glog.logmsg("%s\n", versionstring(GAAEM_VERSION, __TIME__, __DATE__).c_str()); glog.logmsg("Usage: %s control_file_name\n", argv[0]); @@ -849,7 +849,7 @@ class cAllAtOnceInverter{ mpirank = mpicomm.rank(); std::string ControlFile = std::string(argv[1]); - if(exists(ControlFile) == false){ + if (exists(ControlFile) == false) { glog.logmsg(0, "%s\n", commandlinestring(argc, argv).c_str()); glog.logmsg(0, "%s\n", versionstring(GAAEM_VERSION, __TIME__, __DATE__).c_str()); glog.logmsg(0, "Controlfile %s was not found\n", ControlFile.c_str()); @@ -860,16 +860,16 @@ class cAllAtOnceInverter{ Control = cBlock(ControlFile); OutputOp = cOutputOptions(Control.findblock("Output")); std::string s = strprint(".%04d", mpirank); - + OutputOp.LogFile = insert_after_filename(OutputOp.LogFile, s); glog.open(OutputOp.LogFile); - glog.logmsg(0,"Opening log file %s\n", OutputOp.LogFile.c_str()); - glog.logmsg(0,"Logfile opened on %s\n", timestamp().c_str()); - glog.logmsg(0,"Control file %s\n", Control.Filename.c_str()); - glog.logmsg(0,"Version %s Compiled at %s on %s\n", GAAEM_VERSION, __TIME__, __DATE__); - glog.logmsg(0,"Working directory %s\n", getcurrentdirectory().c_str()); - glog.logmsg(0,"Processes=%lu\tRank=%lu\n", mpisize, mpirank); - glog.logmsg(0,"Processor name = %s\n", mpipname.c_str()); + glog.logmsg(0, "Opening log file %s\n", OutputOp.LogFile.c_str()); + glog.logmsg(0, "Logfile opened on %s\n", timestamp().c_str()); + glog.logmsg(0, "Control file %s\n", Control.Filename.c_str()); + glog.logmsg(0, "Version %s Compiled at %s on %s\n", GAAEM_VERSION, __TIME__, __DATE__); + glog.logmsg(0, "Working directory %s\n", getcurrentdirectory().c_str()); + glog.logmsg(0, "Processes=%lu\tRank=%lu\n", mpisize, mpirank); + glog.logmsg(0, "Processor name = %s\n", mpipname.c_str()); if (mpirank == 0) Control.print(); glog.log(Control.get_as_string()); @@ -882,28 +882,28 @@ class cAllAtOnceInverter{ std::vector bv = Control.findblocks("EMSystem"); T.resize(bv.size()); - for (size_t i = 0; i < bv.size(); i++){ - T[i].initialise(bv[i]); + for (size_t i = 0; i < bv.size(); i++) { + T[i].initialise(bv[i]); std::string stmfile = bv[i].getstringvalue("SystemFile"); std::string str = T[i].T.STM.get_as_string(); - glog.logmsg(0, "==============System file %s\n", stmfile.c_str()); - glog.logmsg(0, str.c_str()); + glog.logmsg(0, "==============System file %s\n", stmfile.c_str()); + glog.logmsg(0, str.c_str()); } - glog.logmsg(0, "==========================================================================\n"); + glog.logmsg(0, "==========================================================================\n"); nchan = calculate_nchan(); - glog.logmsg(0, "\nStarting setup\n"); + glog.logmsg(0, "\nStarting setup\n"); setup(); - glog.logmsg(0, "\nStarting iterations\n"); + glog.logmsg(0, "\nStarting iterations\n"); iterate(); - glog.logmsg(0, "\nFinishing at at %s\n", timestamp().c_str()); - glog.logmsg(0, "Elapsed time = %.2lf\n", stopwatch.etimenow()); + glog.logmsg(0, "\nFinishing at at %s\n", timestamp().c_str()); + glog.logmsg(0, "Elapsed time = %.2lf\n", stopwatch.etimenow()); glog.close(); }; - ~cAllAtOnceInverter(){} + ~cAllAtOnceInverter() {} - bool get_columns(){ + bool get_columns() { cBlock b = Control.findblock("Input.Columns"); fdsurvey = cField(b, "SurveyNumber"); @@ -913,30 +913,30 @@ class cAllAtOnceInverter{ fdfiducial = cField(b, "FidNumber"); fdx = cField(b, "Easting"); fdy = cField(b, "Northing"); - fdelevation = cField(b, "GroundElevation"); + fdelevation = cField(b, "GroundElevation"); return true; } - bool get_geometry_columns(){ + bool get_geometry_columns() { G.resize(10); cBlock g = Control.findblock("Input.Geometry"); - for (size_t i = 0; i < G.size(); i++){ + for (size_t i = 0; i < G.size(); i++) { std::string fname = cTDEmGeometry::element_name(i); - cBlock b = g.findblock(fname); - if (b.Name.size() == 0){ + cBlock b = g.findblock(fname); + if (b.Name.size() == 0) { glog.logmsg(0, "Could not find block for geometry parameter %s\n", fname.c_str()); - std::string e = strprint("Error: exception thrown from %s (%d) %s\n", __FILE__, __LINE__, __FUNCTION__); + std::string e = strprint("Error: exception thrown from %s (%d) %s\n", __FILE__, __LINE__, __FUNCTION__); throw(e); } G[i] = cGeometryInfo(b); - } + } UGI = unknown_geometry_indices(); return true; } - size_t count_closer_than(const std::vector distance, const double& value){ + size_t count_closer_than(const std::vector distance, const double& value) { size_t nn = std::count_if( distance.begin(), distance.end(), [&value](const double& item) { return item <= value; } @@ -944,58 +944,58 @@ class cAllAtOnceInverter{ return nn; } - void read_conductivity_logs(){ + void read_conductivity_logs() { cBlock b = Control.findblock("ConductivityLogs"); bool use = b.getboolvalue("Use"); - if (use){ + if (use) { ConductivityLogMaximumDistance = b.getdoublevalue("MaximumDistance"); ConductivityLogPercentError = b.getdoublevalue("PercentError"); - glog.logmsg(0, "Reading conductivity logs\n"); + glog.logmsg(0, "Reading conductivity logs\n"); std::string ldir = b.getstringvalue("Directory"); auto flist = getfilelist(ldir, "con"); - for (size_t k = 0; k < flist.size(); k++){ + for (size_t k = 0; k < flist.size(); k++) { cConductivityLog clog(flist[k], true); cPoint p(clog.x, clog.y); - if (InputOp.is_point_included(p)){ + if (InputOp.is_point_included(p)) { std::vector ndistances; std::vector neighbours = RS.findneighbourstopoint(clog.x, clog.y, ndistances, ConductivityLogMaximumDistance); - if (neighbours.size()>0){ + if (neighbours.size() > 0) { clog = cConductivityLog(flist[k], false); ConductivityLogs.push_back(clog); } } } - glog.logmsg(0, "There are %lu conductivity logs available\n", flist.size()); - glog.logmsg(0, "There are %lu conductivity logs close enough to be included in the inversion\n", ConductivityLogs.size()); - for (size_t i = 0; i < ConductivityLogs.size(); i++){ - glog.logmsg(0, "%s\n", ConductivityLogs[i].infostring().c_str()); + glog.logmsg(0, "There are %lu conductivity logs available\n", flist.size()); + glog.logmsg(0, "There are %lu conductivity logs close enough to be included in the inversion\n", ConductivityLogs.size()); + for (size_t i = 0; i < ConductivityLogs.size(); i++) { + glog.logmsg(0, "%s\n", ConductivityLogs[i].infostring().c_str()); } } } - - bool count_samples(){ - if (mpirank == 0){ + bool count_samples() { + + if (mpirank == 0) { FILE* fp = fileopen(InputOp.DataFile, "r"); - if (fp == NULL){ - glog.logmsg(0, "Unable to open input DataFile %s\n", InputOp.DataFile.c_str()); + if (fp == NULL) { + glog.logmsg(0, "Unable to open input DataFile %s\n", InputOp.DataFile.c_str()); std::string e = strprint("Error: exception thrown from %s (%d) %s\n", __FILE__, __LINE__, __FUNCTION__); throw(e); } std::string s; - for (size_t i = 0; i < InputOp.HeaderLines; i++){ + for (size_t i = 0; i < InputOp.HeaderLines; i++) { filegetline(fp, s); } size_t k = 0; - while (filegetline(fp, s)){ - if (k % InputOp.SubSample == 0){ + while (filegetline(fp, s)) { + if (k % InputOp.SubSample == 0) { std::vector tokens = tokenize(s); int line = (int)fdline.get(tokens); - if (InputOp.is_line_included(line)){ + if (InputOp.is_line_included(line)) { cPoint p(fdx.get(tokens), fdy.get(tokens)); - if (InputOp.is_point_included(p)){ + if (InputOp.is_point_included(p)) { filerecordindex.push_back(k); } } @@ -1007,27 +1007,27 @@ class cAllAtOnceInverter{ mpicomm.bcast(filerecordindex); nsamples = filerecordindex.size(); - if (nsamples == 0){ - glog.logmsg(0, "There were no samples in the included lines and/or line ranges and/or polygon\n"); + if (nsamples == 0) { + glog.logmsg(0, "There were no samples in the included lines and/or line ranges and/or polygon\n"); std::string e = strprint("Error: exception thrown from %s (%d) %s\n", __FILE__, __LINE__, __FUNCTION__); throw e; } ndata = calculate_ndata(); - nparampersample = calculate_nparampersample(); + nparampersample = calculate_nparampersample(); nparam = calculate_nparam(); sown.set_petsc_default(mpisize, mpirank, (PetscInt)nsamples); //Make sure the ownserships are intergral number of samples - cOwnership down((PetscInt)(sown.start*nchan), (PetscInt)(sown.end*nchan)); - cOwnership pown((PetscInt)(sown.start*nparampersample), (PetscInt)(sown.end*nparampersample)); - nlocaldata = down.nlocal(); + cOwnership down((PetscInt)(sown.start * nchan), (PetscInt)(sown.end * nchan)); + cOwnership pown((PetscInt)(sown.start * nparampersample), (PetscInt)(sown.end * nparampersample)); + nlocaldata = down.nlocal(); nlocalparam = pown.nlocal(); return true; } - bool allocate_data_arrays(){ + bool allocate_data_arrays() { size_t nl = sown.nlocal(); fdsurvey.resize(nl); @@ -1039,20 +1039,20 @@ class cAllAtOnceInverter{ fdy.resize(nl); fdelevation.resize(nl); - for (size_t gi = 0; gi < G.size(); gi++){ + for (size_t gi = 0; gi < G.size(); gi++) { G[gi].ref.resize(nl); - if (G[gi].solve){ + if (G[gi].solve) { G[gi].std.resize(nl); } } - for (size_t k = 0; k < T.size(); k++){ + for (size_t k = 0; k < T.size(); k++) { T[k].allocate_data_arrays(nl); } E.allocate_data_arrays(nl); return true; }; - bool parse(const std::vector tokens, const size_t localindex){ + bool parse(const std::vector tokens, const size_t localindex) { fdsurvey.parse(tokens, localindex); fddate.parse(tokens, localindex); fdflight.parse(tokens, localindex); @@ -1062,14 +1062,14 @@ class cAllAtOnceInverter{ fdy.parse(tokens, localindex); fdelevation.parse(tokens, localindex); - for (size_t gi = 0; gi < G.size(); gi++){ + for (size_t gi = 0; gi < G.size(); gi++) { G[gi].ref.parse(tokens, localindex); - if (G[gi].solve){ + if (G[gi].solve) { G[gi].std.parse(tokens, localindex); } } - for (size_t ti = 0; ti < T.size(); ti++){ + for (size_t ti = 0; ti < T.size(); ti++) { T[ti].parse(tokens, localindex); } @@ -1077,7 +1077,7 @@ class cAllAtOnceInverter{ return true; } - bool read_data(){ + bool read_data() { RS.x.resize(nsamples); RS.y.resize(nsamples); RS.elevation.resize(nsamples); @@ -1085,8 +1085,8 @@ class cAllAtOnceInverter{ allocate_data_arrays(); FILE* fp = fileopen(InputOp.DataFile, "r"); - if (fp == NULL){ - glog.logmsg(0, "Unable to open input DataFile %s\n", InputOp.DataFile.c_str()); + if (fp == NULL) { + glog.logmsg(0, "Unable to open input DataFile %s\n", InputOp.DataFile.c_str()); std::string e = strprint("Error: exception thrown from %s (%d) %s\n", __FILE__, __LINE__, __FUNCTION__); throw e; } @@ -1094,13 +1094,13 @@ class cAllAtOnceInverter{ std::string s; size_t rec = 0; size_t gsi = 0; - while (filegetline(fp, s)){ - if (rec == filerecordindex[gsi]){ + while (filegetline(fp, s)) { + if (rec == filerecordindex[gsi]) { std::vector tokens = tokenize(s); RS.x[gsi] = fdx.get(tokens); RS.y[gsi] = fdy.get(tokens); RS.elevation[gsi] = fdelevation.get(tokens); - if (sown.owns((PetscInt)gsi)){ + if (sown.owns((PetscInt)gsi)) { size_t lsi = sown.localind((PetscInt)gsi); parse(tokens, lsi); } @@ -1116,19 +1116,19 @@ class cAllAtOnceInverter{ } - void test_write_neighbours(){ + void test_write_neighbours() { size_t k = 15; std::vector ndistances; std::vector n = RS.findneighbourstopoint(RS.x[k], RS.y[k], ndistances); //std::vector n = RS.findneighbourstopoint(630325.2,6405732.5,ndistances,500.0); std::vector xn(n.size()); std::vector yn(n.size()); - for (size_t i = 0; i < n.size(); i++){ + for (size_t i = 0; i < n.size(); i++) { xn[i] = RS.x[n[i]]; yn[i] = RS.y[n[i]]; } - if (mpirank == 0){ + if (mpirank == 0) { write_xy("samples.txt", RS.x, RS.y); write_xy("neighbours.txt", xn, yn); } @@ -1138,83 +1138,83 @@ class cAllAtOnceInverter{ void write_xy(const std::string& filename, const std::vector& x, const std::vector& y) { FILE* fp = fileopen(filename, "w"); - for (size_t i = 0; i < x.size(); i++){ + for (size_t i = 0; i < x.size(); i++) { fprintf(fp, "%lf,%lf\n", x[i], y[i]); } fclose(fp); } - size_t calculate_nchan(){ + size_t calculate_nchan() { size_t n = 0; - for (size_t i = 0; i < T.size(); i++){ + for (size_t i = 0; i < T.size(); i++) { n += T[i].ndata(); } return n; } - inline size_t calculate_ndata(){ - return nsamples*nchan; + inline size_t calculate_ndata() { + return nsamples * nchan; } - inline size_t calculate_nparampersample(){ - return nlayers + UGI.size(); + inline size_t calculate_nparampersample() { + return nlayers + UGI.size(); } - size_t calculate_nparam(){ - return nsamples*nparampersample; + size_t calculate_nparam() { + return nsamples * nparampersample; } - size_t dindex(const size_t& iglobalsample, const size_t& ichan){ - return iglobalsample*nchan + ichan; + size_t dindex(const size_t& iglobalsample, const size_t& ichan) { + return iglobalsample * nchan + ichan; } - size_t gpindex_c(const size_t& iglobalsample, const size_t& ilayer){ + size_t gpindex_c(const size_t& iglobalsample, const size_t& ilayer) { return iglobalsample * nparampersample + ilayer; } - size_t gpindex_g(const size_t& iglobalsample, const size_t& igparam){ - return iglobalsample*nparampersample + nlayers + igparam; + size_t gpindex_g(const size_t& iglobalsample, const size_t& igparam) { + return iglobalsample * nparampersample + nlayers + igparam; } - size_t localsampleindex(const size_t& _globalsampleindex){ + size_t localsampleindex(const size_t& _globalsampleindex) { return _globalsampleindex - (size_t)sown.start; } - size_t globalsampleindex(const size_t& _localsampleindex){ + size_t globalsampleindex(const size_t& _localsampleindex) { return _localsampleindex + (size_t)sown.start; } - std::vector local_data(){ + std::vector local_data() { std::vector v; v.reserve(nlocaldata); - for (size_t si = 0; si < (size_t)sown.nlocal(); si++){ - for (size_t ti = 0; ti < T.size(); ti++){ - append(v,T[ti].data(si)); + for (size_t si = 0; si < (size_t)sown.nlocal(); si++) { + for (size_t ti = 0; ti < T.size(); ti++) { + append(v, T[ti].data(si)); } } return v; } - std::vector local_noise(){ + std::vector local_noise() { std::vector v; - v.reserve(nlocaldata); - for (size_t si = 0; si < (size_t)sown.nlocal(); si++){ - for (size_t ti = 0; ti < T.size(); ti++){ - append(v,T[ti].noise(si)); + v.reserve(nlocaldata); + for (size_t si = 0; si < (size_t)sown.nlocal(); si++) { + for (size_t ti = 0; ti < T.size(); ti++) { + append(v, T[ti].noise(si)); } } return v; } - std::vector local_mref(){ + std::vector local_mref() { std::vector v(nlocalparam); size_t i = 0; - for (size_t si = 0; si < (size_t)sown.nlocal(); si++){ - for (size_t li = 0; li < nlayers; li++){ + for (size_t si = 0; si < (size_t)sown.nlocal(); si++) { + for (size_t li = 0; li < nlayers; li++) { v[i] = std::log10(E.cref(si, li)); i++; } - for (size_t gi = 0; gi < UGI.size(); gi++){ + for (size_t gi = 0; gi < UGI.size(); gi++) { v[i] = G[UGI[gi]].ref(si); i++; } @@ -1222,45 +1222,45 @@ class cAllAtOnceInverter{ return v; } - std::vector local_mstd(){ + std::vector local_mstd() { std::vector v(nlocalparam); size_t i = 0; - for (size_t si = 0; si < (size_t)sown.nlocal(); si++){ - for (size_t li = 0; li < nlayers; li++){ + for (size_t si = 0; si < (size_t)sown.nlocal(); si++) { + for (size_t li = 0; li < nlayers; li++) { v[i] = E.cstd(si, li); i++; } - for (size_t gi = 0; gi < UGI.size(); gi++){ + for (size_t gi = 0; gi < UGI.size(); gi++) { v[i] = G[UGI[gi]].std(si); - i++; + i++; } } return v; } - void J_create(){ + void J_create() { - glog.logmsg(0, "Creating matrix J\n"); + glog.logmsg(0, "Creating matrix J\n"); J.create_sparse("J", mpicomm, (PetscInt)nlocaldata, (PetscInt)nlocalparam, (PetscInt)ndata, (PetscInt)nparam); std::vector d_nnz(J.nlocalrows(), 0); std::vector o_nnz(J.nlocalrows(), 0); - for (size_t si = (size_t)sown.start; si < (size_t)sown.end; si++){ - for (size_t ci = 0; ci < nchan; ci++){ + for (size_t si = (size_t)sown.start; si < (size_t)sown.end; si++) { + for (size_t ci = 0; ci < nchan; ci++) { PetscInt di = (PetscInt)dindex(si, ci); PetscInt lri = J.lri((PetscInt)di); - for (size_t li = 0; li < nlayers; li++){ + for (size_t li = 0; li < nlayers; li++) { PetscInt pi = (PetscInt)gpindex_c(si, li); - if (J.inownerdiagonalblock(di, pi)){ + if (J.inownerdiagonalblock(di, pi)) { d_nnz[lri]++; } else o_nnz[lri]++;; } - for (size_t gi = 0; gi < UGI.size(); gi++){ + for (size_t gi = 0; gi < UGI.size(); gi++) { PetscInt pi = (PetscInt)gpindex_g(si, gi); - if (J.inownerdiagonalblock(di, pi)){ + if (J.inownerdiagonalblock(di, pi)) { d_nnz[lri]++; } - else o_nnz[lri]++; + else o_nnz[lri]++; } } @@ -1268,55 +1268,55 @@ class cAllAtOnceInverter{ J.preallocate(d_nnz, o_nnz); } - void J_set_nonzero_pattern(){ - glog.logmsg(0, "Assembling matrix G\n"); + void J_set_nonzero_pattern() { + glog.logmsg(0, "Assembling matrix G\n"); PetscErrorCode ierr; - for (size_t si = (size_t)sown.start; si < (size_t)sown.end; si++){ - for (size_t ci = 0; ci < nchan; ci++){ + for (size_t si = (size_t)sown.start; si < (size_t)sown.end; si++) { + for (size_t ci = 0; ci < nchan; ci++) { PetscInt di = (PetscInt)dindex(si, ci); - for (size_t li = 0; li < nlayers; li++){ + for (size_t li = 0; li < nlayers; li++) { PetscInt pi = (PetscInt)gpindex_c(si, li); ierr = MatSetValue(J.mat(), di, pi, 0.0, INSERT_VALUES); CHKERR(ierr); } - for (size_t gi = 0; gi < UGI.size(); gi++){ + for (size_t gi = 0; gi < UGI.size(); gi++) { PetscInt pi = (PetscInt)gpindex_g(si, gi); ierr = MatSetValue(J.mat(), di, pi, 0.0, INSERT_VALUES); CHKERR(ierr); } } } J.assemble(); - glog.logmsg(0, "Finished assembling matrix G\n"); + glog.logmsg(0, "Finished assembling matrix G\n"); return; } - void Wd_create(){ - glog.logmsg(0, "Creating matrix Wd\n"); + void Wd_create() { + glog.logmsg(0, "Creating matrix Wd\n"); Wd.create_diagonal_to_power("Wd", dstd, -2.0); Wd *= (1.0 / Wd.nglobalrows()); - glog.logmsg(0, "Finished creating matrix Wd\n"); + glog.logmsg(0, "Finished creating matrix Wd\n"); return; }; - void Wr_create(){ - glog.logmsg(0, "Creating matrix Wr\n"); + void Wr_create() { + glog.logmsg(0, "Creating matrix Wr\n"); Wr.create_diagonal_to_power("Wr", mstd, -2.0); Wr *= (InversionOp.AlphaR / Wr.nglobalrows()); - glog.logmsg(0, "Finished creating matrix Wr\n"); + glog.logmsg(0, "Finished creating matrix Wr\n"); return; }; - void V_create_1st_derivative(){ + void V_create_1st_derivative() { - glog.logmsg(0, "Creating matrix V\n"); - - size_t nglobalconstraints = (nlayers - 1)*nsamples; - size_t nlocalconstraints = (nlayers - 1)*sown.nlocal(); + glog.logmsg(0, "Creating matrix V\n"); + + size_t nglobalconstraints = (nlayers - 1) * nsamples; + size_t nlocalconstraints = (nlayers - 1) * sown.nlocal(); V.create_sparse("V", mpicomm, (PetscInt)nlocalconstraints, (PetscInt)nlocalparam, (PetscInt)nglobalconstraints, (PetscInt)nparam); V.preallocate(2, 0); PetscInt gri = V.gri(0); - for (size_t si = (size_t)sown.start; si < (size_t)sown.end; si++){ - for (size_t li = 0; li < nlayers - 1; li++){ + for (size_t si = (size_t)sown.start; si < (size_t)sown.end; si++) { + for (size_t li = 0; li < nlayers - 1; li++) { PetscInt pa = (PetscInt)gpindex_c(si, li); PetscInt pb = (PetscInt)gpindex_c(si, li + 1); V.set(gri, pa, 1.0); @@ -1326,20 +1326,20 @@ class cAllAtOnceInverter{ } V.assemble(); V *= std::sqrt(InversionOp.AlphaV / (double)V.nglobalrows()); - glog.logmsg(0, "Finished creating matrix V\n"); + glog.logmsg(0, "Finished creating matrix V\n"); }; - void V_create_2nd_derivative(){ + void V_create_2nd_derivative() { - glog.logmsg(0, "Creating matrix V\n"); - size_t nglobalconstraints = (nlayers - 2)*nsamples; - size_t nlocalconstraints = (nlayers - 2)*sown.nlocal(); + glog.logmsg(0, "Creating matrix V\n"); + size_t nglobalconstraints = (nlayers - 2) * nsamples; + size_t nlocalconstraints = (nlayers - 2) * sown.nlocal(); V.create_sparse("V", mpicomm, (PetscInt)nlocalconstraints, (PetscInt)nlocalparam, (PetscInt)nglobalconstraints, (PetscInt)nparam); V.preallocate(3, 0); PetscInt gri = V.gri(0); - for (size_t si = (size_t)sown.start; si < (size_t)sown.end; si++){ - for (size_t li = 1; li < nlayers - 1; li++){ + for (size_t si = (size_t)sown.start; si < (size_t)sown.end; si++) { + for (size_t li = 1; li < nlayers - 1; li++) { PetscInt pa = (PetscInt)gpindex_c(si, li - 1); PetscInt pb = (PetscInt)gpindex_c(si, li); PetscInt pc = (PetscInt)gpindex_c(si, li + 1); @@ -1351,12 +1351,12 @@ class cAllAtOnceInverter{ } V.assemble(); V *= std::sqrt(InversionOp.AlphaV / (double)V.nglobalrows()); - glog.logmsg(0, "Finished creating matrix V\n"); + glog.logmsg(0, "Finished creating matrix V\n"); }; - void H_create_elevation(){ + void H_create_elevation() { - glog.logmsg(0, "Creating matrix H\n"); + glog.logmsg(0, "Creating matrix H\n"); std::vector t = get_thicknesses_ref(0); std::vector d = get_interface_depths(t); d.push_back(d.back() + t.back()); @@ -1365,23 +1365,23 @@ class cAllAtOnceInverter{ std::vector d_nnz(H.nlocalrows(), 0); std::vector o_nnz(H.nlocalrows(), 0); size_t nsum = 0; - for (size_t gsi = (size_t)sown.start; gsi < (size_t)sown.end; gsi++){ + for (size_t gsi = (size_t)sown.start; gsi < (size_t)sown.end; gsi++) { double elev = RS.elevation[gsi]; std::vector ndistance; std::vector neighbours = RS.findneighbours(gsi, ndistance); nsum += neighbours.size(); - for (size_t li = 0; li < nlayers; li++){ + for (size_t li = 0; li < nlayers; li++) { PetscInt gri = (PetscInt)gpindex_c(gsi, li); PetscInt lri = (PetscInt)H.lri(gri); d_nnz[lri]++;//non-zero for this samples layer - for (size_t ni = 0; ni < neighbours.size(); ni++){ + for (size_t ni = 0; ni < neighbours.size(); ni++) { size_t ngsi = neighbours[ni]; double nelev = RS.elevation[ngsi]; std::vector nd = d + (elev - nelev); std::vector fo = fractionaloverlaps(d[li], d[li + 1], nd); - for (size_t nli = 0; nli < nlayers; nli++){ - if (fo[nli]>0){ + for (size_t nli = 0; nli < nlayers; nli++) { + if (fo[nli] > 0) { PetscInt gci = (PetscInt)gpindex_c(ngsi, nli); if (H.inownerdiagonalblock(gri, gci)) d_nnz[lri]++; else o_nnz[lri]++; @@ -1392,36 +1392,36 @@ class cAllAtOnceInverter{ } H.preallocate(d_nnz, o_nnz); nsum = mpicomm.sum(nsum); - glog.logmsg(0, "Average number of neighbours per sample = %lu\n", nsum / nsamples); + glog.logmsg(0, "Average number of neighbours per sample = %lu\n", nsum / nsamples); //Set entries - glog.logmsg(0, "Assembling matrix H\n"); + glog.logmsg(0, "Assembling matrix H\n"); //loop over each sample int count = 0; - for (size_t gsi = (size_t)sown.start; gsi < (size_t)sown.end; gsi++){ + for (size_t gsi = (size_t)sown.start; gsi < (size_t)sown.end; gsi++) { double elev = RS.elevation[gsi]; //Get neighbours and set weights std::vector ndistances; std::vector neighbours = RS.findneighbours(gsi, ndistances); std::vector nweights(neighbours.size(), 0.0); - for (size_t ni = 0; ni < neighbours.size(); ni++){ + for (size_t ni = 0; ni < neighbours.size(); ni++) { //In case there are coincident samples make minimum distance 10 m if (ndistances[ni] < 10.0)ndistances[ni] = 10.0; - nweights[ni] = std::pow(ndistances[ni], -1.0*InversionOp.InverseDistancePower); + nweights[ni] = std::pow(ndistances[ni], -1.0 * InversionOp.InverseDistancePower); } nweights /= sum(nweights); //loop over this sample's layers - for (size_t li = 0; li < nlayers; li++){ + for (size_t li = 0; li < nlayers; li++) { PetscInt gri = (PetscInt)gpindex_c(gsi, li); //Loop over each neighbour/layer double wsum = 0.0; std::vector> values(neighbours.size()); - for (size_t ni = 0; ni < neighbours.size(); ni++){ + for (size_t ni = 0; ni < neighbours.size(); ni++) { size_t ngsi = neighbours[ni]; double nelev = RS.elevation[ngsi]; std::vector nd = d + (elev - nelev); @@ -1436,9 +1436,9 @@ class cAllAtOnceInverter{ PetscInt gci = (PetscInt)gpindex_c(gsi, li); PetscErrorCode ierr = MatSetValue(H.mat(), gri, gci, -wsum, INSERT_VALUES); CHKERR(ierr); //Set values for neighbours - for (size_t ni = 0; ni < neighbours.size(); ni++){ + for (size_t ni = 0; ni < neighbours.size(); ni++) { size_t ngsi = neighbours[ni]; - for (size_t nli = 0; nli < nlayers; nli++){ + for (size_t nli = 0; nli < nlayers; nli++) { if (values[ni][nli] == 0.0) continue; gci = (PetscInt)gpindex_c(ngsi, nli); PetscErrorCode ierr = MatSetValue(H.mat(), gri, gci, values[ni][nli], INSERT_VALUES); CHKERR(ierr); @@ -1449,35 +1449,35 @@ class cAllAtOnceInverter{ }//loop over samples H.assemble(); H *= std::sqrt(InversionOp.AlphaH / (double)H.nglobalrows()); - glog.logmsg(0, "Finished assembling matrix H\n"); + glog.logmsg(0, "Finished assembling matrix H\n"); return; } - void B_create_elevation_new(){ + void B_create_elevation_new() { - double clogerr = 0.5*(std::log10(100.0 + ConductivityLogPercentError) - std::log10(100.0 - ConductivityLogPercentError)); - glog.logmsg(0, "Creating matrix B\n"); + double clogerr = 0.5 * (std::log10(100.0 + ConductivityLogPercentError) - std::log10(100.0 - ConductivityLogPercentError)); + glog.logmsg(0, "Creating matrix B\n"); std::vector lthick = get_thicknesses_ref(0); std::vector ldepth = get_interface_depths(lthick); ldepth.push_back(ldepth.back() + lthick.back()); PetscInt nconstraints = 0; - for (size_t bi = 0; bi < ConductivityLogs.size(); bi++){ + for (size_t bi = 0; bi < ConductivityLogs.size(); bi++) { cConductivityLog& clog = ConductivityLogs[bi]; std::vector ndistances; std::vector neighbours = RS.findneighbourstopoint(clog.x, clog.y, ndistances, ConductivityLogMaximumDistance); - for (size_t ni = 0; ni < neighbours.size(); ni++){ + for (size_t ni = 0; ni < neighbours.size(); ni++) { size_t ngsi = neighbours[ni]; double nelev = RS.elevation[ngsi]; std::vector nd = ldepth + (clog.z - nelev); - for (size_t li = 0; li < nlayers; li++){ - if (clog.interval_has_overlap(nd[li], nd[li + 1])){ + for (size_t li = 0; li < nlayers; li++) { + if (clog.interval_has_overlap(nd[li], nd[li + 1])) { nconstraints += 1; } } } } - glog.logmsg(0, "nconductivitylogconstraints=%d\n", nconstraints); + glog.logmsg(0, "nconductivitylogconstraints=%d\n", nconstraints); B.create_sparse("B", mpicomm, PETSC_DECIDE, (PetscInt)nlocalparam, (PetscInt)nconstraints, (PetscInt)nparam); PetscInt nlocalconstraints = (PetscInt)B.nlocalrows(); @@ -1490,19 +1490,19 @@ class cAllAtOnceInverter{ std::vector d_nnz(B.nlocalrows(), 0); std::vector o_nnz(B.nlocalrows(), 0); PetscInt gri = 0; - for (size_t bi = 0; bi < ConductivityLogs.size(); bi++){ + for (size_t bi = 0; bi < ConductivityLogs.size(); bi++) { cConductivityLog& clog = ConductivityLogs[bi]; //Get neighbours std::vector ndistances; std::vector neighbours = RS.findneighbourstopoint(clog.x, clog.y, ndistances, ConductivityLogMaximumDistance); - for (size_t ni = 0; ni < neighbours.size(); ni++){ + for (size_t ni = 0; ni < neighbours.size(); ni++) { size_t ngsi = neighbours[ni]; double nelev = RS.elevation[ngsi]; std::vector nd = ldepth + (clog.z - nelev); - for (size_t li = 0; li < nlayers; li++){ - if (clog.interval_has_overlap(nd[li], nd[li + 1])){ - if (B.ownsrow(gri)){ + for (size_t li = 0; li < nlayers; li++) { + if (clog.interval_has_overlap(nd[li], nd[li + 1])) { + if (B.ownsrow(gri)) { PetscInt lri = B.lri(gri); PetscInt gci = (PetscInt)gpindex_c(ngsi, li); if (B.inownerdiagonalblock(gri, gci)) d_nnz[lri]++; @@ -1517,37 +1517,37 @@ class cAllAtOnceInverter{ //Set values double* clogref_loc = clogref.getlocalarray(); - double* clogstd_loc = clogstd.getlocalarray(); + double* clogstd_loc = clogstd.getlocalarray(); gri = 0; - for (size_t bi = 0; bi < ConductivityLogs.size(); bi++){ + for (size_t bi = 0; bi < ConductivityLogs.size(); bi++) { cConductivityLog& clog = ConductivityLogs[bi]; std::vector ndistances; std::vector neighbours = RS.findneighbourstopoint(clog.x, clog.y, ndistances, ConductivityLogMaximumDistance); - + std::vector nweights(neighbours.size()); - for (size_t ni = 0; ni < neighbours.size(); ni++){ + for (size_t ni = 0; ni < neighbours.size(); ni++) { size_t ngsi = neighbours[ni]; double nelev = RS.elevation[ngsi]; std::vector nd = ldepth + (clog.z - nelev); - + //In case there are coincident positions make minimum distance 10 m if (ndistances[ni] < 10.0) ndistances[ni] = 10.0; - nweights[ni] = std::pow(ndistances[ni], -1.0*InversionOp.InverseDistancePower); + nweights[ni] = std::pow(ndistances[ni], -1.0 * InversionOp.InverseDistancePower); - for (size_t li = 0; li < nlayers; li++){ + for (size_t li = 0; li < nlayers; li++) { PetscInt gci = (PetscInt)gpindex_c(ngsi, li); - if (clog.interval_has_overlap(nd[li], nd[li + 1])){ - if (B.ownsrow(gri)){ + if (clog.interval_has_overlap(nd[li], nd[li + 1])) { + if (B.ownsrow(gri)) { PetscInt lri = (PetscInt)B.lri(gri); - + //Set rhs vector's value size_t npoints; double linear_mean, log10_mean; clog.interval_means(nd[li], nd[li + 1], npoints, linear_mean, log10_mean); clogref_loc[lri] = std::log10(log10_mean); clogstd_loc[lri] = clogerr;//Converted to log10 decades - + //Set matrix entry B.set(gri, gci, 1.0); }//is owner @@ -1561,28 +1561,28 @@ class cAllAtOnceInverter{ B.assemble(); Wb.create_diagonal_to_power("Wb", clogstd, -2.0); Wb *= (InversionOp.AlphaB / Wb.nglobalrows()); - glog.logmsg(0, "Finished creating matrix B\n"); + glog.logmsg(0, "Finished creating matrix B\n"); } - void B_create_elevation(){ + void B_create_elevation() { - double clogerr = 0.5*(std::log10(100.0 + ConductivityLogPercentError) - std::log10(100.0 - ConductivityLogPercentError)); - glog.logmsg(0, "Creating matrix B\n"); + double clogerr = 0.5 * (std::log10(100.0 + ConductivityLogPercentError) - std::log10(100.0 - ConductivityLogPercentError)); + glog.logmsg(0, "Creating matrix B\n"); std::vector t = get_thicknesses_ref(0); std::vector d = get_interface_depths(t); d.push_back(d.back() + t.back()); PetscInt nconstraints = 0; - for (size_t bi = 0; bi < ConductivityLogs.size(); bi++){ + for (size_t bi = 0; bi < ConductivityLogs.size(); bi++) { cConductivityLog& clog = ConductivityLogs[bi]; std::vector ndistances; std::vector neighbours = RS.findneighbourstopoint(clog.x, clog.y, ndistances, ConductivityLogMaximumDistance); if (neighbours.size() == 0)continue; - for (size_t li = 0; li < nlayers; li++){ + for (size_t li = 0; li < nlayers; li++) { bool hasoverlap = clog.interval_has_overlap(d[li], d[li + 1]); if (hasoverlap) nconstraints += 1; } } - glog.logmsg(0, "nconductivitylogconstraints=%d\n", nconstraints); + glog.logmsg(0, "nconductivitylogconstraints=%d\n", nconstraints); B.create_sparse("B", mpicomm, PETSC_DECIDE, (PetscInt)nlocalparam, (PetscInt)nconstraints, (PetscInt)nparam); PetscInt nlocalconstraints = B.nlocalrows(); @@ -1592,7 +1592,7 @@ class cAllAtOnceInverter{ std::vector d_nnz(B.nlocalrows(), 0); std::vector o_nnz(B.nlocalrows(), 0); PetscInt gri = 0; - for (size_t bi = 0; bi < ConductivityLogs.size(); bi++){ + for (size_t bi = 0; bi < ConductivityLogs.size(); bi++) { cConductivityLog& clog = ConductivityLogs[bi]; //Get neighbours @@ -1600,18 +1600,18 @@ class cAllAtOnceInverter{ std::vector neighbours = RS.findneighbourstopoint(clog.x, clog.y, ndistances, ConductivityLogMaximumDistance); if (neighbours.size() == 0)continue; - for (size_t li = 0; li < nlayers; li++){ + for (size_t li = 0; li < nlayers; li++) { bool hasoverlap = clog.interval_has_overlap(d[li], d[li + 1]); - if (hasoverlap){ - if (B.ownsrow(gri)){ + if (hasoverlap) { + if (B.ownsrow(gri)) { PetscInt lri = B.lri(gri); - for (size_t ni = 0; ni < neighbours.size(); ni++){ + for (size_t ni = 0; ni < neighbours.size(); ni++) { size_t ngsi = neighbours[ni]; double nelev = RS.elevation[ngsi]; std::vector nd = d + (clog.z - nelev); std::vector fo = fractionaloverlaps(d[li], d[li + 1], nd); - for (size_t nli = 0; nli < fo.size(); nli++){ - if (fo[nli]>0){ + for (size_t nli = 0; nli < fo.size(); nli++) { + if (fo[nli] > 0) { PetscInt gci = (PetscInt)gpindex_c(ngsi, nli); if (B.inownerdiagonalblock(gri, gci)) d_nnz[lri]++; else o_nnz[lri]++; @@ -1631,7 +1631,7 @@ class cAllAtOnceInverter{ double* clogref_loc = clogref.getlocalarray(); double* clogstd_loc = clogstd.getlocalarray(); gri = -1; - for (size_t bi = 0; bi < ConductivityLogs.size(); bi++){ + for (size_t bi = 0; bi < ConductivityLogs.size(); bi++) { cConductivityLog& clog = ConductivityLogs[bi]; std::vector ndistances; std::vector neighbours = RS.findneighbourstopoint(clog.x, clog.y, ndistances, ConductivityLogMaximumDistance); @@ -1639,15 +1639,15 @@ class cAllAtOnceInverter{ //Set neighbour weights std::vector nweights(neighbours.size()); - for (size_t ni = 0; ni < neighbours.size(); ni++){ + for (size_t ni = 0; ni < neighbours.size(); ni++) { //In case there are coincident positions make minimum distance 10 m if (ndistances[ni] < 10.0) ndistances[ni] = 10.0; - nweights[ni] = std::pow(ndistances[ni], -1.0*InversionOp.InverseDistancePower); + nweights[ni] = std::pow(ndistances[ni], -1.0 * InversionOp.InverseDistancePower); } nweights /= sum(nweights); - for (size_t li = 0; li < nlayers; li++){ + for (size_t li = 0; li < nlayers; li++) { size_t npoints; double linear_mean, log10_mean; bool hasoverlap = clog.interval_means(d[li], d[li + 1], npoints, linear_mean, log10_mean); @@ -1659,7 +1659,7 @@ class cAllAtOnceInverter{ //Loop over each neighbour/layer double wsum = 0.0; std::vector> values(neighbours.size()); - for (size_t ni = 0; ni < neighbours.size(); ni++){ + for (size_t ni = 0; ni < neighbours.size(); ni++) { size_t ngsi = neighbours[ni]; double nelev = RS.elevation[ngsi]; std::vector nd = d + (clog.z - nelev); @@ -1671,7 +1671,7 @@ class cAllAtOnceInverter{ } //if wsum is tiny there is very little overlap with any neighbour - if (wsum < 0.001){ + if (wsum < 0.001) { clogref_loc[lri] = 0.0; clogstd_loc[lri] = clogerr;//Converted to log10 decades continue; @@ -1682,9 +1682,9 @@ class cAllAtOnceInverter{ clogstd_loc[lri] = clogerr;//Converted to log10 decades //Set values for neighbours - for (size_t ni = 0; ni < neighbours.size(); ni++){ + for (size_t ni = 0; ni < neighbours.size(); ni++) { size_t ngsi = neighbours[ni]; - for (size_t nli = 0; nli < nlayers; nli++){ + for (size_t nli = 0; nli < nlayers; nli++) { if (values[ni][nli] == 0.0) continue; values[ni][nli] /= wsum; PetscInt gci = (PetscInt)gpindex_c(ngsi, nli); @@ -1697,55 +1697,55 @@ class cAllAtOnceInverter{ clogstd.restorelocalarray(clogstd_loc); B.assemble(); Wb.create_diagonal_to_power("Wb", clogstd, -2.0); - if (Wb.nglobalrows() > 0){ + if (Wb.nglobalrows() > 0) { Wb *= (InversionOp.AlphaB / Wb.nglobalrows()); } - glog.logmsg(0, "Finished creating matrix B\n"); + glog.logmsg(0, "Finished creating matrix B\n"); } - void report_matrix_memory_usage(){ - glog.logmsg(0, "J matrix global memory %.3lf MiB\n", J.globalmemory() / 1e6); - glog.logmsg(0, "Wd matrix global memory %.3lf MiB\n", Wd.globalmemory() / 1e6); - glog.logmsg(0, "V matrix global memory %.3lf MiB\n", V.globalmemory() / 1e6); - glog.logmsg(0, "H matrix global memory %.3lf MiB\n", H.globalmemory() / 1e6); - glog.logmsg(0, "B matrix global memory %.3lf MiB\n", B.globalmemory() / 1e6); - glog.logmsg(0, "Wb matrix global memory %.3lf MiB\n", Wb.globalmemory() / 1e6); - glog.logmsg(0, "Wr matrix global memory %.3lf MiB\n", Wr.globalmemory() / 1e6); - glog.logmsg(0, "P matrix global memory %.3lf MiB\n", P.globalmemory() / 1e6); + void report_matrix_memory_usage() { + glog.logmsg(0, "J matrix global memory %.3lf MiB\n", J.globalmemory() / 1e6); + glog.logmsg(0, "Wd matrix global memory %.3lf MiB\n", Wd.globalmemory() / 1e6); + glog.logmsg(0, "V matrix global memory %.3lf MiB\n", V.globalmemory() / 1e6); + glog.logmsg(0, "H matrix global memory %.3lf MiB\n", H.globalmemory() / 1e6); + glog.logmsg(0, "B matrix global memory %.3lf MiB\n", B.globalmemory() / 1e6); + glog.logmsg(0, "Wb matrix global memory %.3lf MiB\n", Wb.globalmemory() / 1e6); + glog.logmsg(0, "Wr matrix global memory %.3lf MiB\n", Wr.globalmemory() / 1e6); + glog.logmsg(0, "P matrix global memory %.3lf MiB\n", P.globalmemory() / 1e6); double total = J.globalmemory() + Wd.globalmemory() + V.globalmemory() + H.globalmemory() + B.globalmemory() + Wb.globalmemory() + Wr.globalmemory() + P.globalmemory(); - glog.logmsg(0, "Total matrix global memory %.3lf MiB\n", total / 1e6); + glog.logmsg(0, "Total matrix global memory %.3lf MiB\n", total / 1e6); } - std::vector get_conductivity_model(const size_t& localsampleindex, const double* mlocal){ + std::vector get_conductivity_model(const size_t& localsampleindex, const double* mlocal) { std::vector c(nlayers); size_t pi = nparampersample * localsampleindex; - for (size_t li = 0; li < nlayers; li++){ - c[li] = std::pow(10.0, mlocal[pi+li]); + for (size_t li = 0; li < nlayers; li++) { + c[li] = std::pow(10.0, mlocal[pi + li]); } return c; } - std::vector get_thicknesses_ref(const size_t localsampleindex){ + std::vector get_thicknesses_ref(const size_t localsampleindex) { std::vector thickness(nlayers - 1); - for (size_t i = 0; i < nlayers - 1; i++){ + for (size_t i = 0; i < nlayers - 1; i++) { thickness[i] = E.tref(localsampleindex, i); } return thickness; } - std::vector get_interface_depths(const std::vector& thickness){ + std::vector get_interface_depths(const std::vector& thickness) { std::vector d(thickness.size() + 1); d[0] = 0; - for (size_t i = 0; i < thickness.size(); i++){ + for (size_t i = 0; i < thickness.size(); i++) { d[i + 1] = d[i] + thickness[i]; } return d; } - cTDEmGeometry get_geometry_ref(const size_t localsampleindex){ + cTDEmGeometry get_geometry_ref(const size_t localsampleindex) { cTDEmGeometry geom; - for (size_t gi = 0; gi < G.size(); gi++){ + for (size_t gi = 0; gi < G.size(); gi++) { geom[gi] = G[gi].ref(localsampleindex); } return geom; @@ -1753,25 +1753,25 @@ class cAllAtOnceInverter{ cTDEmGeometry get_geometry_model(const size_t& localsampleindex, const double* mlocal) { - cTDEmGeometry geom = get_geometry_ref(localsampleindex); + cTDEmGeometry geom = get_geometry_ref(localsampleindex); size_t pi = nparampersample * localsampleindex; - for (size_t gi = 0; gi < UGI.size(); gi++){ - geom[UGI[gi]] = mlocal[pi + nlayers + gi]; + for (size_t gi = 0; gi < UGI.size(); gi++) { + geom[UGI[gi]] = mlocal[pi + nlayers + gi]; } - return geom; + return geom; } - std::vector unknown_geometry_indices(){ + std::vector unknown_geometry_indices() { std::vector indices; - for (size_t gi = 0; gi < G.size(); gi++){ + for (size_t gi = 0; gi < G.size(); gi++) { if (G[gi].solve) indices.push_back(gi); } return indices; } - bool forwardmodel_and_jacobian(const cPetscDistVector& m, cPetscDistVector& g, const bool computejacobian){ - - static double natlog10 = std::log(10.0); + bool forwardmodel_and_jacobian(const cPetscDistVector& m, cPetscDistVector& g, const bool computejacobian) { + + static double natlog10 = std::log(10.0); std::vector predicted; std::vector> derivatives; @@ -1781,35 +1781,35 @@ class cAllAtOnceInverter{ double* glocal = g.getlocalarray(); cOwnership gdist = g.ownership(); //glog.logmsg(0, "Starting forward modelling\n"); - for (size_t si = (size_t)sown.start; si < (size_t)sown.end; si++){ + for (size_t si = (size_t)sown.start; si < (size_t)sown.end; si++) { size_t lsi = sown.localind((PetscInt)si); //size_t gpi = gpindex_c(si, 0); //size_t lpi = mdist.localind(gpi); - - std::vector conductivity = get_conductivity_model(lsi,mlocal); - std::vector thickness = get_thicknesses_ref(lsi); - cTDEmGeometry geometry = get_geometry_model(lsi,mlocal); - + + std::vector conductivity = get_conductivity_model(lsi, mlocal); + std::vector thickness = get_thicknesses_ref(lsi); + cTDEmGeometry geometry = get_geometry_model(lsi, mlocal); + size_t gdi = dindex(si, 0); size_t ldi = gdist.localind((PetscInt)gdi); - for (size_t ti = 0; ti < T.size(); ti++){ + for (size_t ti = 0; ti < T.size(); ti++) { T[ti].forward_model_and_derivatives(conductivity, thickness, geometry, predicted, derivatives, computejacobian, UGI); - for (size_t k = 0; k < predicted.size(); k++){ + for (size_t k = 0; k < predicted.size(); k++) { glocal[ldi + k] = predicted[k]; } - if (computejacobian){ - for (size_t k = 0; k < derivatives.size(); k++){ - for (size_t li = 0; li < nlayers; li++){ + if (computejacobian) { + for (size_t k = 0; k < derivatives.size(); k++) { + for (size_t li = 0; li < nlayers; li++) { //multiply by natural log(10) as parameters are in log base 10 units const double val = natlog10 * conductivity[li] * derivatives[k][li]; - J.set((PetscInt)(gdi + k), (PetscInt)gpindex_c(si,li), val); + J.set((PetscInt)(gdi + k), (PetscInt)gpindex_c(si, li), val); } - for (size_t gi = 0; gi < UGI.size(); gi++){ - const double val = derivatives[k][nlayers+gi]; - J.set((PetscInt)(gdi + k), (PetscInt)gpindex_g(si,gi), val); + for (size_t gi = 0; gi < UGI.size(); gi++) { + const double val = derivatives[k][nlayers + gi]; + J.set((PetscInt)(gdi + k), (PetscInt)gpindex_g(si, gi), val); } } } @@ -1820,20 +1820,20 @@ class cAllAtOnceInverter{ m.restorelocalreadonlyarray(mlocal); g.restorelocalarray(glocal); if (computejacobian) J.assemble(); - + //glog.logmsg(0, "Finished forward modelling\n"); return true; } - bool setup(){ + bool setup() { count_samples(); - glog.logmsg(0, "nsamples=%lu\n", nsamples); - glog.logmsg(0, "nchannels=%lu\n", nchan); - glog.logmsg(0, "ndata=%lu\n", ndata); - glog.logmsg(0, "nlayers=%lu\n", nlayers); - glog.logmsg(0, "nparampersample=%lu\n", nparampersample); - glog.logmsg(0, "nparam=%lu\n", nparam); + glog.logmsg(0, "nsamples=%lu\n", nsamples); + glog.logmsg(0, "nchannels=%lu\n", nchan); + glog.logmsg(0, "ndata=%lu\n", ndata); + glog.logmsg(0, "nlayers=%lu\n", nlayers); + glog.logmsg(0, "nparampersample=%lu\n", nparampersample); + glog.logmsg(0, "nparam=%lu\n", nparam); read_data(); read_conductivity_logs(); @@ -1866,10 +1866,10 @@ class cAllAtOnceInverter{ //clogref.writetextfile("clogref.vec"); //clogstd.writetextfile("clogstd.vec"); - if (InversionOp.VerticalSmoothnessMethod == SM_1ST_DERIVATIVE){ + if (InversionOp.VerticalSmoothnessMethod == SM_1ST_DERIVATIVE) { V_create_1st_derivative(); } - else if (InversionOp.VerticalSmoothnessMethod == SM_2ND_DERIVATIVE){ + else if (InversionOp.VerticalSmoothnessMethod == SM_2ND_DERIVATIVE) { V_create_2nd_derivative(); } //V.writetextfile("V.smat"); @@ -1884,13 +1884,13 @@ class cAllAtOnceInverter{ J_set_nonzero_pattern(); //J.writetextfile("J.smat"); - glog.logmsg(0, "Creating preconditioner\n"); + glog.logmsg(0, "Creating preconditioner\n"); P.create_identity("P", mpicomm, (PetscInt)nlocalparam, (PetscInt)nparam); //P.writetextfile("P.smat"); - glog.logmsg(0, "Finished creating preconditioner\n"); + glog.logmsg(0, "Finished creating preconditioner\n"); report_matrix_memory_usage(); - + return true; }; @@ -1908,27 +1908,27 @@ class cAllAtOnceInverter{ // Ax=b cPetscDistVector x(xVec); cPetscDistVector b(bVec); - b = (J ^ (Wd*(J*x))); - if (InversionOp.AlphaV > 0.0)b += ((V ^ (V*x)) *= lambda); - if (InversionOp.AlphaH > 0.0)b += ((H ^ (H*x)) *= lambda); - if (InversionOp.AlphaB > 0.0)b += ((B ^ (Wb*(B*x))) *= lambda); - if (InversionOp.AlphaR > 0.0)b += ((Wr*x) *= lambda); + b = (J ^ (Wd * (J * x))); + if (InversionOp.AlphaV > 0.0)b += ((V ^ (V * x)) *= lambda); + if (InversionOp.AlphaH > 0.0)b += ((H ^ (H * x)) *= lambda); + if (InversionOp.AlphaB > 0.0)b += ((B ^ (Wb * (B * x))) *= lambda); + if (InversionOp.AlphaR > 0.0)b += ((Wr * x) *= lambda); return 0; } - cPetscDistVector cg_solve(cPetscDistShellMatrix& A, const cPetscDistVector& m, const cPetscDistVector& g){ + cPetscDistVector cg_solve(cPetscDistShellMatrix& A, const cPetscDistVector& m, const cPetscDistVector& g) { cPetscDistVector b("b", mpicomm, (PetscInt)nlocalparam, (PetscInt)nparam); - glog.logmsg(0, "Starting CG solve\n"); + glog.logmsg(0, "Starting CG solve\n"); cStopWatch sw; - b = J ^ (Wd*(dobs - g + J*m)); - if (InversionOp.AlphaB > 0) b += lambda*((B ^ (Wb^clogref))); - if (InversionOp.AlphaR > 0) b += lambda*((Wr^mref)); + b = J ^ (Wd * (dobs - g + J * m)); + if (InversionOp.AlphaB > 0) b += lambda * ((B ^ (Wb ^ clogref))); + if (InversionOp.AlphaR > 0) b += lambda * ((Wr ^ mref)); cPetscDistVector mtrial = A.solve_CG(P, b, m); - glog.logmsg(0, "Finished CG solve time=%lf\n", sw.etimenow()); - glog.logmsg(0, A.convergence_summary().c_str()); + glog.logmsg(0, "Finished CG solve time=%lf\n", sw.etimenow()); + glog.logmsg(0, A.convergence_summary().c_str()); return (mtrial - m); }; @@ -1937,82 +1937,82 @@ class cAllAtOnceInverter{ // Ax=b cPetscDistVector x(xVec); cPetscDistVector b(bVec); - b = (J ^ (Wd*(J*x))); - if (InversionOp.AlphaV > 0.0)b += ((V ^ (V*x)) *= lambda); - if (InversionOp.AlphaH > 0.0)b += ((H ^ (H*x)) *= lambda); - if (InversionOp.AlphaB > 0.0)b += ((B ^ (Wb*(B*x))) *= lambda); - if (InversionOp.AlphaR > 0.0)b += ((Wr*x) *= lambda); + b = (J ^ (Wd * (J * x))); + if (InversionOp.AlphaV > 0.0)b += ((V ^ (V * x)) *= lambda); + if (InversionOp.AlphaH > 0.0)b += ((H ^ (H * x)) *= lambda); + if (InversionOp.AlphaB > 0.0)b += ((B ^ (Wb * (B * x))) *= lambda); + if (InversionOp.AlphaR > 0.0)b += ((Wr * x) *= lambda); return 0; } - - cPetscDistVector cg_solve_dm(cPetscDistShellMatrix& A, const cPetscDistVector& m, const cPetscDistVector& g){ + + cPetscDistVector cg_solve_dm(cPetscDistShellMatrix& A, const cPetscDistVector& m, const cPetscDistVector& g) { cPetscDistVector b("b", mpicomm, (PetscInt)nlocalparam, (PetscInt)nparam); - glog.logmsg(0, "Starting CG solve\n"); + glog.logmsg(0, "Starting CG solve\n"); cStopWatch sw; - b = J ^ (Wd*(dobs - g)); - if (InversionOp.AlphaV > 0) b -= lambda*((V ^ (V*m))); - if (InversionOp.AlphaH > 0) b -= lambda*((H ^ (H*m))); - if (InversionOp.AlphaB > 0) b -= lambda*((B ^ (B*m))); - if (InversionOp.AlphaB > 0) b += lambda*((B ^ (Wb^clogref))); - if (InversionOp.AlphaR > 0) b += lambda*((Wr^ (mref-m))); + b = J ^ (Wd * (dobs - g)); + if (InversionOp.AlphaV > 0) b -= lambda * ((V ^ (V * m))); + if (InversionOp.AlphaH > 0) b -= lambda * ((H ^ (H * m))); + if (InversionOp.AlphaB > 0) b -= lambda * ((B ^ (B * m))); + if (InversionOp.AlphaB > 0) b += lambda * ((B ^ (Wb ^ clogref))); + if (InversionOp.AlphaR > 0) b += lambda * ((Wr ^ (mref - m))); cPetscDistVector mtrial = A.solve_CG(P, b, m); - glog.logmsg(0, "Finished CG solve time=%lf\n", sw.etimenow()); - glog.logmsg(0, A.convergence_summary().c_str()); + glog.logmsg(0, "Finished CG solve time=%lf\n", sw.etimenow()); + glog.logmsg(0, A.convergence_summary().c_str()); return (mtrial - m); }; - double PhiD(const cPetscDistVector& g){ + double PhiD(const cPetscDistVector& g) { return Wd.vtAv(dobs - g); } - double PhiV(const cPetscDistVector& m){ - return (V*m).l2_norm(); + double PhiV(const cPetscDistVector& m) { + return (V * m).l2_norm(); } - double PhiH(const cPetscDistVector& m){ - return (H*m).l2_norm(); + double PhiH(const cPetscDistVector& m) { + return (H * m).l2_norm(); } - double PhiB(const cPetscDistVector& m){ - return Wb.vtAv(B*m - clogref); + double PhiB(const cPetscDistVector& m) { + return Wb.vtAv(B * m - clogref); } - double PhiR(const cPetscDistVector& m){ + double PhiR(const cPetscDistVector& m) { return Wr.vtAv(m - mref); } - void find_stepfactor(const cPetscDistVector& m, const cPetscDistVector& dm, const double& currentphid, const double& targetphid, double& bestsf, double& bestphid, double& improvement){ + void find_stepfactor(const cPetscDistVector& m, const cPetscDistVector& dm, const double& currentphid, const double& targetphid, double& bestsf, double& bestphid, double& improvement) { - glog.logmsg(0, "Finding step factor\n"); + glog.logmsg(0, "Finding step factor\n"); cStopWatch sw; cPetscDistVector gtrial = dobs; cPetscDistVector mtrial = m; cInversionLineSearcher LS(targetphid); double sf; - while (LS.next_x(sf)){ - mtrial = m + sf*dm; + while (LS.next_x(sf)) { + mtrial = m + sf * dm; forwardmodel_and_jacobian(mtrial, gtrial, false); double phid = PhiD(gtrial); LS.add_pair(sf, phid); } //size_t index = LS.nearest_index(bestsf, bestphid); - auto best = LS.nearest(); - glog.logmsg(0, "Find stepfactor time=%lf\n", sw.etimenow()); - improvement = 100.0*(currentphid - bestphid) / currentphid; - glog.logmsg(0, "Step factor = %.5lf\n", best.first); - glog.logmsg(0, "Found PhiD = %.5lf\n", best.second); - glog.logmsg(0, "Improvement = %.5lf%%\n", improvement); - + auto best = LS.nearest(); + glog.logmsg(0, "Find stepfactor time=%lf\n", sw.etimenow()); + improvement = 100.0 * (currentphid - bestphid) / currentphid; + glog.logmsg(0, "Step factor = %.5lf\n", best.first); + glog.logmsg(0, "Found PhiD = %.5lf\n", best.second); + glog.logmsg(0, "Improvement = %.5lf%%\n", improvement); + //if (mpirank == 0){ // std::string stepsfile = strprint("output//steps//steps_%02llu.txt", mLastIteration); // LS.writetextfile(stepsfile); //} - + if (mpirank == 0) { std::string stepsfile = strprint("output//steps//steps_%02llu.txt", mLastIteration); std::ofstream of(stepsfile); @@ -2022,8 +2022,8 @@ class cAllAtOnceInverter{ return; } - void iterate(){ - + void iterate() { + cPetscDistVector m = mref; m.setname("m"); @@ -2034,12 +2034,12 @@ class cAllAtOnceInverter{ bool keepgoing = true; size_t iteration = 1; - while (keepgoing){ - glog.logmsg(0, "\n\nItaration = %lu\n", iteration); + while (keepgoing) { + glog.logmsg(0, "\n\nItaration = %lu\n", iteration); cStopWatch sw; cPetscDistVector g("g", mpicomm, (PetscInt)nlocaldata, (PetscInt)ndata); forwardmodel_and_jacobian(m, g, true); - glog.logmsg(0, "Forward modelling time=%lf\n", sw.etimenow()); + glog.logmsg(0, "Forward modelling time=%lf\n", sw.etimenow()); double phiv = PhiV(m); double phih = PhiH(m); double phib = PhiB(m); double phir = PhiR(m); @@ -2047,53 +2047,53 @@ class cAllAtOnceInverter{ double phi = phid + phiv + phih + phib + phir; double targetphid = phid * 0.7; if (targetphid < InversionOp.MinimumPhiD) targetphid = InversionOp.MinimumPhiD; - + log_iteration_msg(lambda, phi, phiv, phih, phib, phir, phid, targetphid); - - cPetscDistVector dm = cg_solve(A, m, g); + + cPetscDistVector dm = cg_solve(A, m, g); double bestsf, bestphid, improvement; find_stepfactor(m, dm, phid, targetphid, bestsf, bestphid, improvement); - + iteration++; - if (improvement > 0.0){ - m += bestsf*dm; + if (improvement > 0.0) { + m += bestsf * dm; forwardmodel_and_jacobian(m, g, false); mLastPhiD = bestphid; mLastLambda = lambda; mLastIteration = iteration; - write_results(OutputOp.DataFile, m, g); + write_results(OutputOp.DataFile, m, g); } - if (improvement < InversionOp.MinimumPercentageImprovement){ + if (improvement < InversionOp.MinimumPercentageImprovement) { keepgoing = false; //lambda = lambda * 0.7; } - if (iteration > InversionOp.MaximumIterations){ + if (iteration > InversionOp.MaximumIterations) { keepgoing = false; } - if (bestphid < InversionOp.MinimumPhiD){ + if (bestphid < InversionOp.MinimumPhiD) { keepgoing = false; } } }; - void log_iteration_msg(const double& lam, const double& phi, const double& phiv, const double& phih, const double& phib, const double& phir, const double& phid, const double& targetphid){ - glog.logmsg(0, "Current Lambda = %lf\n", lam); - glog.logmsg(0, "Current Phi = %lf\n", phi); - glog.logmsg(0, "Current PhiV = %lf\n", phiv); - glog.logmsg(0, "Current PhiH = %lf\n", phih); - glog.logmsg(0, "Current PhiB = %lf\n", phib); - glog.logmsg(0, "Current PhiR = %lf\n", phir); - glog.logmsg(0, "Current PhiD = %lf\n", phid); - glog.logmsg(0, "Target PhiD = %lf\n", targetphid); - } - - void write_results(const std::string& filename, const cPetscDistVector& m, const cPetscDistVector& g){ - - for (int p = 0; p < mpisize; p++){ - if (p == mpirank){ - if (mpirank == 0){ + void log_iteration_msg(const double& lam, const double& phi, const double& phiv, const double& phih, const double& phib, const double& phir, const double& phid, const double& targetphid) { + glog.logmsg(0, "Current Lambda = %lf\n", lam); + glog.logmsg(0, "Current Phi = %lf\n", phi); + glog.logmsg(0, "Current PhiV = %lf\n", phiv); + glog.logmsg(0, "Current PhiH = %lf\n", phih); + glog.logmsg(0, "Current PhiB = %lf\n", phib); + glog.logmsg(0, "Current PhiR = %lf\n", phir); + glog.logmsg(0, "Current PhiD = %lf\n", phid); + glog.logmsg(0, "Target PhiD = %lf\n", targetphid); + } + + void write_results(const std::string& filename, const cPetscDistVector& m, const cPetscDistVector& g) { + + for (int p = 0; p < mpisize; p++) { + if (p == mpirank) { + if (mpirank == 0) { FILE* fp = fileopen(filename, "w"); fclose(fp); } @@ -2104,8 +2104,8 @@ class cAllAtOnceInverter{ } - void append_my_results(const std::string& filename, const cPetscDistVector& m, const cPetscDistVector& g){ - + void append_my_results(const std::string& filename, const cPetscDistVector& m, const cPetscDistVector& g) { + const double* mlocal = m.getlocalreadonlyarray(); const double* glocal = g.getlocalreadonlyarray(); //cOwnership mdist = m.ownership(); @@ -2115,11 +2115,11 @@ class cAllAtOnceInverter{ cOutputFileInfo OI; std::string buf; FILE* fp = fileopen(filename, "a"); - for (size_t lsi = 0; lsi < (size_t)sown.nlocal(); lsi++){ + for (size_t lsi = 0; lsi < (size_t)sown.nlocal(); lsi++) { size_t gsi = sown.globalind((PetscInt)lsi); //size_t lpi = mdist.localind(gpindex_c(gsi, 0)); std::vector conductivity = get_conductivity_model(lsi, mlocal); - std::vector thickness = get_thicknesses_ref(lsi); + std::vector thickness = get_thicknesses_ref(lsi); thickness.push_back(thickness.back()); cTDEmGeometry gref = get_geometry_ref(lsi); cTDEmGeometry ginv = get_geometry_model(lsi, mlocal); @@ -2157,16 +2157,16 @@ class cAllAtOnceInverter{ OI.addfield("elevation", 'F', 10, 2); OI.setunits("m"); OI.setdescription("Ground elevation relative to sea-level"); buf += strprint("%10.2lf", fdelevation(lsi)); - - for (size_t gi = 0; gi < G.size(); gi++){ + + for (size_t gi = 0; gi < G.size(); gi++) { OI.addfield(gref.element_name(gi), 'F', 9, 2); OI.setunits(gref.units(gi)); OI.setdescription(gref.description(gi)); buf += strprint("%9.2lf", gref[gi]); } - - for (size_t gi = 0; gi < UGI.size(); gi++){ - OI.addfield("inverted_"+ginv.element_name(UGI[gi]), 'F', 9, 2); + + for (size_t gi = 0; gi < UGI.size(); gi++) { + OI.addfield("inverted_" + ginv.element_name(UGI[gi]), 'F', 9, 2); OI.setunits(ginv.units(UGI[gi])); OI.setdescription("Inverted " + ginv.description(UGI[gi])); buf += strprint("%9.2lf", ginv[UGI[gi]]); @@ -2178,95 +2178,95 @@ class cAllAtOnceInverter{ OI.addfield("conductivity", 'E', 15, 6, nlayers); OI.setunits("S/m"); OI.setdescription("Layer conductivity"); - for (size_t li = 0; li < nlayers; li++){ + for (size_t li = 0; li < nlayers; li++) { buf += strprint("%15.6le", conductivity[li]); } OI.addfield("thickness", 'F', 9, 2, nlayers); OI.setunits("m"); OI.setdescription("Layer thickness"); - for (size_t li = 0; li < nlayers; li++){ + for (size_t li = 0; li < nlayers; li++) { buf += strprint("%9.2lf", thickness[li]); } - if (OutputOp.PositiveLayerBottomDepths){ + if (OutputOp.PositiveLayerBottomDepths) { OI.addfield("depth_bottom", 'F', 9, 2, nlayers); OI.setunits("m"); OI.setdescription("Depth to bottom of layer"); double tsum = 0.0; - for (size_t i = 0; i < nlayers; i++){ + for (size_t i = 0; i < nlayers; i++) { buf += strprint("%9.2lf", tsum); tsum += thickness[i]; } } - if (OutputOp.NegativeLayerBottomDepths){ + if (OutputOp.NegativeLayerBottomDepths) { OI.addfield("depth_bottom_negative", 'F', 9, 2, nlayers); OI.setunits("m"); OI.setdescription("Negative of depth to bottom of layer"); double tsum = 0.0; - for (size_t i = 0; i < nlayers; i++){ + for (size_t i = 0; i < nlayers; i++) { tsum += thickness[i]; buf += strprint("%9.2lf", -tsum); } } - if (OutputOp.InterfaceElevations){ + if (OutputOp.InterfaceElevations) { OI.addfield("elevation_interfaces", 'F', 9, 2, nlayers); OI.setunits("m"); OI.setdescription("Elevation of interfaces"); double etop = fdelevation(lsi); - for (size_t i = 0; i < nlayers; i++){ + for (size_t i = 0; i < nlayers; i++) { buf += strprint("%9.2lf", etop); etop -= thickness[i]; } } char cid[3] = { 'X', 'Y', 'Z' }; - if (OutputOp.ObservedData){ - for (size_t si = 0; si < T.size(); si++){ - cSystemInfo& S = T[si]; + if (OutputOp.ObservedData) { + for (size_t si = 0; si < T.size(); si++) { + cSystemInfo& S = T[si]; std::string sys = strprint("EMSystem_%lu_", si + 1); - for (size_t ci = 0; ci < S.Comp.size(); ci++){ + for (size_t ci = 0; ci < S.Comp.size(); ci++) { cTDEmComponentInfo& C = S.Comp[ci]; if (C.Use == false)continue; - if (S.InvertTotalField){ + if (S.InvertTotalField) { OI.addfield("observed_" + sys + cid[ci] + "P", 'E', 15, 6); OI.setdescription("Observed " + sys + cid[ci] + "-component primary field"); buf += strprint("%15.6le", C.fdp(lsi, 0)); - } + } OI.addfield("observed_" + sys + cid[ci] + "S", 'E', 15, 6, C.nw); OI.setdescription("Observed " + sys + cid[ci] + "-component secondary field windows"); - for (size_t w = 0; w < C.nw; w++){ - buf += strprint("%15.6le", C.fds(lsi,w)); - } + for (size_t w = 0; w < C.nw; w++) { + buf += strprint("%15.6le", C.fds(lsi, w)); + } } } } - if (OutputOp.Noise){ - for (size_t si = 0; si < T.size(); si++){ + if (OutputOp.Noise) { + for (size_t si = 0; si < T.size(); si++) { cSystemInfo& S = T[si]; std::string sys = strprint("EMSystem_%lu_", si + 1); - for (size_t ci = 0; ci < S.Comp.size(); ci++){ + for (size_t ci = 0; ci < S.Comp.size(); ci++) { cTDEmComponentInfo& C = S.Comp[ci]; - if (C.Use == false)continue; + if (C.Use == false)continue; OI.addfield("noise_" + sys + cid[ci] + "S", 'E', 15, 6, C.nw); OI.setdescription("Estimated noise " + sys + cid[ci] + "-component secondary field windows"); - for (size_t w = 0; w < C.nw; w++){ + for (size_t w = 0; w < C.nw; w++) { buf += strprint("%15.6le", C.fdn(lsi, w)); } } } } - if (OutputOp.PredictedData){ + if (OutputOp.PredictedData) { size_t ldi = gdist.localind((PetscInt)dindex(gsi, 0)); - for (size_t si = 0; si < T.size(); si++){ + for (size_t si = 0; si < T.size(); si++) { cSystemInfo& S = T[si]; S.forward_model(conductivity, thickness, ginv); std::string sys = strprint("EMSystem_%lu_", si + 1); - for (size_t ci = 0; ci < S.Comp.size(); ci++){ + for (size_t ci = 0; ci < S.Comp.size(); ci++) { cTDEmComponentInfo& C = S.Comp[ci]; if (C.Use == false)continue; - if (S.InvertTotalField){ - OI.addfield("predicted_" + sys + cid[ci] + "P", 'E', 15, 6); + if (S.InvertTotalField) { + OI.addfield("predicted_" + sys + cid[ci] + "P", 'E', 15, 6); OI.setdescription("Predicted " + sys + cid[ci] + "-component primary field"); if (ci == 0) buf += strprint("%15.6le", S.T.PrimaryX); else if (ci == 1) buf += strprint("%15.6le", S.T.PrimaryY); @@ -2274,16 +2274,16 @@ class cAllAtOnceInverter{ OI.addfield("predicted_" + sys + cid[ci] + "S", 'E', 15, 6, C.nw); OI.setdescription("Predicted " + sys + cid[ci] + "-component secondary field windows"); - for (size_t w = 0; w < C.nw; w++){ - if (ci == 0) buf += strprint("%15.6le", S.T.X[w]); + for (size_t w = 0; w < C.nw; w++) { + if (ci == 0) buf += strprint("%15.6le", S.T.X[w]); else if (ci == 1) buf += strprint("%15.6le", S.T.Y[w]); else buf += strprint("%15.6le", S.T.Z[w]); } } - else{ + else { OI.addfield("predicted_" + sys + cid[ci] + "S", 'E', 15, 6, C.nw); OI.setdescription("Predicted " + sys + cid[ci] + "-component secondary field windows"); - for (size_t w = 0; w < C.nw; w++){ + for (size_t w = 0; w < C.nw; w++) { buf += strprint("%15.6le", glocal[ldi]); ldi++; } @@ -2291,7 +2291,7 @@ class cAllAtOnceInverter{ } } } - + //Inversion parameters OI.addfield("SamplePhiD", 'E', 15, 6); @@ -2312,14 +2312,14 @@ class cAllAtOnceInverter{ //Carriage return buf += strprint("\n"); - if ((int)lsi == sown.nlocal() - 1 || buf.size() >= 2048){ + if ((int)lsi == sown.nlocal() - 1 || buf.size() >= 2048) { fprintf(fp, buf.c_str()); fflush(fp); buf.resize(0); } OI.lockfields(); - if (lsi == 0 && mpirank == 0){ + if (lsi == 0 && mpirank == 0) { sFilePathParts fpp = getfilepathparts(filename); std::string hdrfile = fpp.directory + fpp.prefix + ".hdr"; OI.write_simple_header(hdrfile); @@ -2333,7 +2333,7 @@ class cAllAtOnceInverter{ fclose(fp); }; - std::vector get_sample_phid(const cPetscDistVector& g){ + std::vector get_sample_phid(const cPetscDistVector& g) { cPetscDistVector nr2 = (g - dobs) / dstd; nr2.pow(2.0); @@ -2341,11 +2341,11 @@ class cAllAtOnceInverter{ const double* nr2local = nr2.getlocalreadonlyarray(); cOwnership nr2dist = nr2.ownership(); std::vector samplephid((size_t)sown.nlocal()); - for (size_t lsi = 0; lsi < (size_t)sown.nlocal(); lsi++){ + for (size_t lsi = 0; lsi < (size_t)sown.nlocal(); lsi++) { size_t gsi = sown.globalind((PetscInt)lsi); size_t ldi = nr2dist.localind((PetscInt)dindex(gsi, 0)); double sum = 0.0; - for (size_t i = 0; i < nchan; i++){ + for (size_t i = 0; i < nchan; i++) { sum += nr2local[ldi + i]; } samplephid[lsi] = sum / nchan; @@ -2357,38 +2357,38 @@ class cAllAtOnceInverter{ }; int main(int argc, char** argv) -{ +{ PetscErrorCode ierr; - try{ - ierr = PetscInitialize(&argc, &argv, PETSC_NULL, PETSC_NULL); + try { + ierr = PetscInitialize(&argc, &argv, (char*)NULL, (char*)NULL); - if(argc < 2){ + if (argc < 2) { glog.logmsg(0, "%s\n", commandlinestring(argc, argv).c_str()); glog.logmsg(0, "%s\n", versionstring(GAAEM_VERSION, __TIME__, __DATE__).c_str()); glog.logmsg(0, "Usage: %s control_file_name\n", argv[0]); glog.logmsg(0, "Too few command line arguments\n"); } - else if(argc > 2){ + else if (argc > 2) { glog.logmsg(0, "%s\n", commandlinestring(argc, argv).c_str()); glog.logmsg(0, "%s\n", versionstring(GAAEM_VERSION, __TIME__, __DATE__).c_str()); glog.logmsg(0, "Usage: %s control_file_name\n", argv[0]); glog.logmsg(0, "Too many command line arguments\n"); } - else{ + else { cAllAtOnceInverter(argc, argv); } ierr = PetscFinalize(); CHKERRQ(ierr); } - catch (const std::string msg){ - glog.logmsg(0,"%s", msg.c_str()); + catch (const std::string msg) { + glog.logmsg(0, "%s", msg.c_str()); ierr = PetscFinalize(); CHKERRQ(ierr); } - catch (const std::runtime_error& e){ - glog.logmsg(0,"%s", e.what()); + catch (const std::runtime_error& e) { + glog.logmsg(0, "%s", e.what()); ierr = PetscFinalize(); CHKERRQ(ierr); } - catch (const std::exception& e){ - glog.logmsg(0,"%s", e.what()); + catch (const std::exception& e) { + glog.logmsg(0, "%s", e.what()); ierr = PetscFinalize(); CHKERRQ(ierr); } diff --git a/src/petsc_wrapper.h b/src/petsc_wrapper.h index f7a9f35..ce06c84 100644 --- a/src/petsc_wrapper.h +++ b/src/petsc_wrapper.h @@ -26,51 +26,56 @@ class cPetscDistShellMatrix;//forward declaration only #define __SDIR__ "src\\" #endif +#if PETSC_VERSION_LT(3,19,0) +// In Petsc Version 3.19.0 PETSC_NULL was deprecated and replaced with PETSC_NULL +#define PETSC_NULLPTR NULL +#endif + #define CHKERR(ierr) cPetscObject::chkerrabort(ierr, __LINE__, __FUNCTION__,__FILE__,__SDIR__) #define dbgprint cPetscObject::debugprint(__LINE__, __FUNCTION__,__FILE__,__SDIR__) -class cSparseMatrixEntries{ +class cSparseMatrixEntries { public: std::vector rowind; std::vector colind; std::vector value; - cSparseMatrixEntries(){ resize(0); } - cSparseMatrixEntries(const PetscInt nnz){ resize(nnz); } + cSparseMatrixEntries() { resize(0); } + cSparseMatrixEntries(const PetscInt nnz) { resize(nnz); } - void resize(const PetscInt nnz){ + void resize(const PetscInt nnz) { rowind.resize(nnz); colind.resize(nnz); value.resize(nnz); - }; + }; }; -class cOwnership{ - +class cOwnership { + public: PetscInt start;//index of first PetscInt end;//index of last + 1 - cOwnership(){ + cOwnership() { start = -1; - end = -1; + end = -1; }; - cOwnership(const PetscInt _start, const PetscInt _end){ + cOwnership(const PetscInt _start, const PetscInt _end) { start = _start; - end = _end; + end = _end; } - cOwnership(const int size, const int rank, const PetscInt nglobal){ + cOwnership(const int size, const int rank, const PetscInt nglobal) { set_petsc_default(size, rank, nglobal); } void set_petsc_default(const int size, const int rank, const PetscInt nglobal) - { + { end = 0; - for (PetscInt i = 0; i <= rank; i++){ + for (PetscInt i = 0; i <= rank; i++) { PetscInt n = nglobal / (PetscInt)size + ((nglobal % (PetscInt)size) > i); end += n; start = end - n; @@ -81,15 +86,15 @@ class cOwnership{ { return end - start; }; //number of local items - - PetscInt globalind(const PetscInt& ilocal) const + + PetscInt globalind(const PetscInt& ilocal) const { - return start + ilocal; + return start + ilocal; } PetscInt localind(const PetscInt& iglobal) const { - return iglobal - start; + return iglobal - start; } bool owns(const PetscInt& iglobal) const @@ -100,25 +105,25 @@ class cOwnership{ }; -class cPetscObject{ - +class cPetscObject { + private: - -public: - - cPetscObject(){ }; + +public: + + cPetscObject() { }; virtual const PetscObject* pobj() const = 0; void setname(const std::string& inname) { - if (pobj()){ + if (pobj()) { PetscErrorCode ierr = PetscObjectSetName(*pobj(), inname.c_str()); CHKERR(ierr); } } const char* cname() const { - if (pobj()){ + if (pobj()) { const char* name; PetscErrorCode ierr = PetscObjectGetName(*pobj(), &name); CHKERR(ierr); @@ -129,7 +134,7 @@ class cPetscObject{ const std::string name() const { - if (cname()){ + if (cname()) { return std::string(cname()); } return std::string("unnamed"); @@ -143,7 +148,7 @@ class cPetscObject{ return comm; } - int mpisize() const + int mpisize() const { int size; PetscErrorCode ierr = MPI_Comm_size(mpicomm(), &size); @@ -151,20 +156,20 @@ class cPetscObject{ return size; } - int mpirank() const + int mpirank() const { int rank; PetscErrorCode ierr = MPI_Comm_rank(mpicomm(), &rank); CHKERR(ierr); return rank; } - PetscErrorCode mpibarrier() const + PetscErrorCode mpibarrier() const { PetscErrorCode ierr = MPI_Barrier(mpicomm()); CHKERR(ierr); return ierr; } - std::string mpiprocname() const + std::string mpiprocname() const { char name[MPI_MAX_PROCESSOR_NAME + 1]; int len; @@ -191,20 +196,20 @@ class cPetscObject{ printf("sizeof(PetscMPIInt) = %zu\n", sizeof(PetscMPIInt)); }; - static void chkerrabort(PetscErrorCode ierr, const int linenumber, const char* functionname, const char* srcfile, const char* srcdirectory){ + static void chkerrabort(PetscErrorCode ierr, const int linenumber, const char* functionname, const char* srcfile, const char* srcdirectory) { do { - if (PetscUnlikely(ierr)){ + if (PetscUnlikely(ierr)) { #if PETSC_VERSION_LT(3,5,0) - PetscError(PETSC_COMM_SELF,linenumber,functionname,srcfile,srcdirectory,ierr,PETSC_ERROR_REPEAT," "); + PetscError(PETSC_COMM_SELF, linenumber, functionname, srcfile, srcdirectory, ierr, PETSC_ERROR_REPEAT, " "); #else - PetscError(PETSC_COMM_SELF,linenumber,functionname,srcfile,ierr,PETSC_ERROR_REPEAT," "); + PetscError(PETSC_COMM_SELF, linenumber, functionname, srcfile, ierr, PETSC_ERROR_REPEAT, " "); #endif MPI_Abort(PETSC_COMM_SELF, ierr); } - } while (0); + } while (0); } - static void debugprint(const int linenumber, const char* functionname, const char* srcfile, const char* srcdirectory){ + static void debugprint(const int linenumber, const char* functionname, const char* srcfile, const char* srcdirectory) { PetscPrintf(PETSC_COMM_WORLD, "**Debug: at line %d : function %s : file %s : directory %s\n", linenumber, functionname, srcfile, srcdirectory); } @@ -219,20 +224,20 @@ class cPetscDistVector : public cPetscObject { private: bool isowner = true; - Vec* pVec = PETSC_NULL; + Vec* pVec = PETSC_NULLPTR; - const PetscObject* pobj() const{ + const PetscObject* pobj() const { if (pVec) return (PetscObject*)pVec; - else return (PetscObject*)PETSC_NULL; + else return (PetscObject*)PETSC_NULLPTR; } void deallocate() { - if (allocated()){ - if (isowner){ + if (allocated()) { + if (isowner) { //printf("Destroying Vector %s\n", name().c_str()); PetscErrorCode ierr = VecDestroy(pVec); CHKERR(ierr); - pVec = PETSC_NULL; + pVec = PETSC_NULLPTR; } } } @@ -240,32 +245,32 @@ class cPetscDistVector : public cPetscObject { void resize(const MPI_Comm comm, const PetscInt _localsize, const PetscInt _globalsize) { std::string myname = name(); - if (allocated()){ - if (localsize() == _localsize && globalsize() == _globalsize){ + if (allocated()) { + if (localsize() == _localsize && globalsize() == _globalsize) { return; } - else deallocate(); + else deallocate(); } create(myname, comm, _localsize, _globalsize); } public: - + cPetscDistVector() {}; cPetscDistVector(const std::string inname, const MPI_Comm comm, const PetscInt localsize, const PetscInt globalsize, const double value = 0.0) - { + { create(inname, comm, localsize, globalsize, value); } - + cPetscDistVector(const MPI_Comm comm, const PetscInt localsize, const PetscInt globalsize) - { - create(comm, localsize, globalsize); + { + create(comm, localsize, globalsize); } cPetscDistVector(const cPetscDistVector& rhs) - { - PetscErrorCode ierr; + { + PetscErrorCode ierr; pVec = new Vec; ierr = VecDuplicate(rhs.vec(), ncvecptr()); CHKERR(ierr); ierr = VecCopy(rhs.vec(), vec()); CHKERR(ierr); @@ -273,13 +278,13 @@ class cPetscDistVector : public cPetscObject { } cPetscDistVector(Vec& vec) - { + { isowner = false; pVec = &vec; } cPetscDistVector(cPetscDistVector&& rhs) - { + { //std::swap(isowner, rhs.isowner); std::swap(pVec, rhs.pVec); } @@ -311,7 +316,7 @@ class cPetscDistVector : public cPetscObject { } void create(const MPI_Comm comm, const PetscInt _localsize, const PetscInt _globalsize) - { + { pVec = new Vec; PetscErrorCode ierr = VecCreateMPI(comm, _localsize, _globalsize, ncvecptr()); CHKERR(ierr); } @@ -319,14 +324,14 @@ class cPetscDistVector : public cPetscObject { void create(const std::string& inname, const MPI_Comm comm, const PetscInt _localsize, const PetscInt _globalsize) { create(comm, _localsize, _globalsize); - setname(inname); + setname(inname); } void create(const std::string& inname, const MPI_Comm comm, const PetscInt _localsize, const PetscInt _globalsize, const double value) { - create(comm, _localsize, _globalsize); - setname(inname); - set(value); + create(comm, _localsize, _globalsize); + setname(inname); + set(value); } cOwnership ownership() const @@ -339,7 +344,7 @@ class cPetscDistVector : public cPetscObject { PetscInt globalsize() const { - if (vec()){ + if (vec()) { PetscInt N; PetscErrorCode ierr = VecGetSize(vec(), &N); CHKERR(ierr); return N; @@ -349,7 +354,7 @@ class cPetscDistVector : public cPetscObject { PetscInt localsize() const { - if (pVec){ + if (pVec) { PetscInt n; PetscErrorCode ierr = VecGetLocalSize(vec(), &n); CHKERR(ierr); return n; @@ -360,7 +365,7 @@ class cPetscDistVector : public cPetscObject { const double* getlocalreadonlyarray() const { const double* localarray; - if (allocated()==false){ + if (allocated() == false) { PetscPrintf(mpicomm(), "Error cPetscDistVector::getlocalarray() (%s) attempting to access NULL vec() pointer\n", name().c_str()); throw(strprint("Error: exception throw from %s (%d) %s\n", __FILE__, __LINE__, __FUNCTION__)); } @@ -376,7 +381,7 @@ class cPetscDistVector : public cPetscObject { double* getlocalarray() { double* localarray; - if (vec() == 0){ + if (vec() == 0) { PetscPrintf(mpicomm(), "Error cPetscDistVector::getlocalarray() (%s) attempting to access NULL vec() pointer\n", cname()); throw(strprint("Error: exception throw from %s (%d) %s\n", __FILE__, __LINE__, __FUNCTION__)); } @@ -399,7 +404,7 @@ class cPetscDistVector : public cPetscObject { void set_local(const std::vector& v) { double* a = getlocalarray(); - for (PetscInt i = 0; igetlocalarray(); } - ~cPetscDistVectorLocalView(){ + ~cPetscDistVectorLocalView() { pv->restorelocalarray(data); } - double& operator[](const size_t& i){ + double& operator[](const size_t& i) { return data[i]; } @@ -692,15 +697,15 @@ class cPetscDistMatrix : public cPetscObject { friend class cPetscDistShellMatrix; private: - - Mat* pMat = PETSC_NULL; + + Mat* pMat = PETSC_NULLPTR; bool isowner = true; std::vector mRowInd; std::vector mColInd; - const PetscObject* pobj() const{ + const PetscObject* pobj() const { if (pMat) return (PetscObject*)pMat; - else return (PetscObject*)PETSC_NULL; + else return (PetscObject*)PETSC_NULLPTR; } bool allocated() const { @@ -710,16 +715,16 @@ class cPetscDistMatrix : public cPetscObject { void deallocate() { - if (allocated()){ - if (isowner){ + if (allocated()) { + if (isowner) { //printf("Destroying Matrix %s\n", name().c_str()); PetscErrorCode ierr = MatDestroy(pMat); CHKERR(ierr); - pMat = PETSC_NULL; + pMat = PETSC_NULLPTR; } } } - bool init_sparse(const MPI_Comm comm, const PetscInt& _nlocalrows, const PetscInt& _nlocalcols, const PetscInt& _nglobalrows, const PetscInt& _nglobalcols){ + bool init_sparse(const MPI_Comm comm, const PetscInt& _nlocalrows, const PetscInt& _nlocalcols, const PetscInt& _nglobalrows, const PetscInt& _nglobalcols) { deallocate(); pMat = new Mat; PetscErrorCode ierr; @@ -731,12 +736,12 @@ class cPetscDistMatrix : public cPetscObject { } public: - + //Constructors & Destructors cPetscDistMatrix() {}; - + cPetscDistMatrix(const cPetscDistMatrix& rhs) - { + { PetscErrorCode ierr; pMat = new Mat; ierr = MatDuplicate(rhs.mat(), MAT_COPY_VALUES, pMat); CHKERR(ierr); @@ -758,7 +763,7 @@ class cPetscDistMatrix : public cPetscObject { Mat& ncmat() const { return *pMat; } const Mat* matptr() const { return pMat; } Mat* ncmatptr() const { return pMat; } - + PetscInt nglobalrows() const { PetscErrorCode ierr; @@ -771,7 +776,7 @@ class cPetscDistMatrix : public cPetscObject { { PetscErrorCode ierr; PetscInt M, N; -ierr = MatGetSize(mat(), &M, &N); CHKERR(ierr); + ierr = MatGetSize(mat(), &M, &N); CHKERR(ierr); return N; }; @@ -779,7 +784,7 @@ ierr = MatGetSize(mat(), &M, &N); CHKERR(ierr); { PetscErrorCode ierr; PetscInt m, n; -ierr = MatGetLocalSize(mat(), &m, &n); CHKERR(ierr); + ierr = MatGetLocalSize(mat(), &m, &n); CHKERR(ierr); return m; }; @@ -790,7 +795,7 @@ ierr = MatGetLocalSize(mat(), &m, &n); CHKERR(ierr); //see http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatMPIAIJSetPreallocation.html PetscErrorCode ierr; PetscInt m, n; -ierr = MatGetLocalSize(mat(), &m, &n); CHKERR(ierr); + ierr = MatGetLocalSize(mat(), &m, &n); CHKERR(ierr); return n; }; @@ -820,21 +825,21 @@ ierr = MatGetLocalSize(mat(), &m, &n); CHKERR(ierr); return rowownership().localind(globalrow); }; - void create_sparse(const std::string& name, const MPI_Comm comm, const PetscInt& _nlocalrows, const PetscInt& _nlocalcols, const PetscInt& _nglobalrows, const PetscInt& _nglobalcols){ + void create_sparse(const std::string& name, const MPI_Comm comm, const PetscInt& _nlocalrows, const PetscInt& _nlocalcols, const PetscInt& _nglobalrows, const PetscInt& _nglobalcols) { init_sparse(comm, _nlocalrows, _nlocalcols, _nglobalrows, _nglobalcols); - setname(name); + setname(name); } bool create_diagonal_structure(const std::string& name, const MPI_Comm comm, const PetscInt& localsize, const PetscInt& globalsize) { create_sparse(name, comm, localsize, localsize, globalsize, globalsize); PetscErrorCode ierr; - ierr = MatMPIAIJSetPreallocation(mat(), 1, PETSC_NULL, 0, PETSC_NULL); CHKERR(ierr); - const cOwnership r = rowownership(); - for (PetscInt gi = r.start; gi& localdiagonal) - { - create_diagonal_structure(name,comm,(PetscInt)localdiagonal.size(),PETSC_DETERMINE); - set_diagonal_local(localdiagonal); + { + create_diagonal_structure(name, comm, (PetscInt)localdiagonal.size(), PETSC_DETERMINE); + set_diagonal_local(localdiagonal); return true; } @@ -892,7 +897,7 @@ ierr = MatGetLocalSize(mat(), &m, &n); CHKERR(ierr); //} //Set the matrix type - if (mattype != "mpidense" && mattype != "elemental"){ + if (mattype != "mpidense" && mattype != "elemental") { PetscPrintf(comm, "Error: cPetscDistMatrix::create_dense() the matrix type should be either be 'mpidense' or 'elemental' \n"); throw(strprint("Error: exception throw from %s (%d) %s\n", __FILE__, __LINE__, __FUNCTION__)); } @@ -922,15 +927,15 @@ ierr = MatGetLocalSize(mat(), &m, &n); CHKERR(ierr); if (_nglobalcols <= 0) _nglobalcols = cmax + 1; PetscInt nnz = (PetscInt)rowind.size(); - if ((int)colind.size() != nnz || (int)values.size() != nnz){ + if ((int)colind.size() != nnz || (int)values.size() != nnz) { printf("cPetscDistMatrix::create_globalindices(...) rowind, colind, and vals must be the same size\n"); throw(strprint("Error: exception throw from %s (%d) %s\n", __FILE__, __LINE__, __FUNCTION__)); } - if (rmin<0 || rmax >= _nglobalrows){ + if (rmin < 0 || rmax >= _nglobalrows) { printf("cPetscDistMatrix::create_globalindices(...) row index out of range\n"); throw(strprint("Error: exception throw from %s (%d) %s\n", __FILE__, __LINE__, __FUNCTION__)); } - if (cmin<0 || cmax >= _nglobalcols){ + if (cmin < 0 || cmax >= _nglobalcols) { printf("cPetscDistMatrix::create_globalindices(...) column index out of range\n"); throw(strprint("Error: exception throw from %s (%d) %s\n", __FILE__, __LINE__, __FUNCTION__)); } @@ -945,13 +950,13 @@ ierr = MatGetLocalSize(mat(), &m, &n); CHKERR(ierr); cOwnership cown = columnownership(); //Preallocate - for (PetscInt i = 0; i nlc)d_nnz[k] = nlc; if (o_nnz[k] > _nglobalcols - nlc)o_nnz[k] = _nglobalcols - nlc; } ierr = MatMPIAIJSetPreallocation(mat(), 0, d_nnz.data(), 0, o_nnz.data()); CHKERR(ierr); //Now set the values - for (PetscInt i = 0; i < nnz; ++i){ - if (rown.owns(rowind[i])){ + for (PetscInt i = 0; i < nnz; ++i) { + if (rown.owns(rowind[i])) { ierr = MatSetValue(mat(), rowind[i], colind[i], values[i], ADD_VALUES); CHKERR(ierr); } } @@ -978,19 +983,20 @@ ierr = MatGetLocalSize(mat(), &m, &n); CHKERR(ierr); bool preallocate(const PetscInt d_nz, const PetscInt o_nz) const { - PetscErrorCode ierr = MatMPIAIJSetPreallocation(mat(), d_nz, PETSC_NULL, o_nz, PETSC_NULL); CHKERR(ierr); + PetscErrorCode ierr = MatMPIAIJSetPreallocation(mat(), d_nz, PETSC_NULLPTR, o_nz, PETSC_NULLPTR); CHKERR(ierr); return true; } bool preallocate(std::vector& d_nnz, std::vector& o_nnz) const { - PetscErrorCode ierr = MatMPIAIJSetPreallocation(mat(), (PetscInt)PETSC_NULL, d_nnz.data(), (PetscInt)PETSC_NULL, o_nnz.data()); + const PetscInt ignore = 0; + PetscErrorCode ierr = MatMPIAIJSetPreallocation(mat(), ignore, d_nnz.data(), ignore, o_nnz.data()); CHKERR(ierr); return true; } bool inownerdiagonalblock(const PetscInt& ri, const PetscInt& ci) const - { + { if (rowownership().owns(ri) && columnownership().owns(ci)) return true; else return false; } @@ -1001,13 +1007,13 @@ ierr = MatGetLocalSize(mat(), &m, &n); CHKERR(ierr); } bool isassembled() const { - if (mat() == PETSC_NULL)return false; + if (mat() == PETSC_NULLPTR)return false; PetscBool status; PetscErrorCode ierr = MatAssembled(mat(), &status); CHKERR(ierr); if (status == PETSC_TRUE)return true; return false; } - + void set(const PetscInt globalrow, const PetscInt globalcol, double value) { PetscErrorCode ierr; @@ -1015,47 +1021,47 @@ ierr = MatGetLocalSize(mat(), &m, &n); CHKERR(ierr); CHKERR(ierr); }; - bool set_diagonal_local(const std::vector& localdiagonal){ - - if (nglobalrows() != nglobalcols()){ + bool set_diagonal_local(const std::vector& localdiagonal) { + + if (nglobalrows() != nglobalcols()) { printf("cPetscDistMatrix::set_diagonal_local() not a square matrix\n"); throw(strprint("Error: exception throw from %s (%d) %s\n", __FILE__, __LINE__, __FUNCTION__)); } - else if((int)localdiagonal.size() != nlocalrows()){ + else if ((int)localdiagonal.size() != nlocalrows()) { printf("cPetscDistMatrix::set_diagonal_local() incompatible size\n"); throw(strprint("Error: exception throw from %s (%d) %s\n", __FILE__, __LINE__, __FUNCTION__)); } - + PetscErrorCode ierr; const cOwnership r = rowownership(); PetscInt k = 0; - for (PetscInt gi = r.start; gi < r.end; gi++){ + for (PetscInt gi = r.start; gi < r.end; gi++) { ierr = MatSetValue(mat(), gi, gi, localdiagonal[k], INSERT_VALUES); CHKERR(ierr); k++; } - assemble(); + assemble(); return true; } - PetscErrorCode setrow(const PetscInt globalrow, const std::vector& globalcols, const std::vector& values) - { + PetscErrorCode setrow(const PetscInt globalrow, const std::vector& globalcols, const std::vector& values) + { PetscErrorCode ierr; - ierr = MatSetValues(mat(),1,&globalrow, (PetscInt)globalcols.size(),globalcols.data(),values.data(),INSERT_VALUES); CHKERR(ierr); + ierr = MatSetValues(mat(), 1, &globalrow, (PetscInt)globalcols.size(), globalcols.data(), values.data(), INSERT_VALUES); CHKERR(ierr); return ierr; } - PetscErrorCode setrowdense(const PetscInt globalrow, const std::vector& globalcols, const std::vector& rowvalues) - { - PetscErrorCode ierr; + PetscErrorCode setrowdense(const PetscInt globalrow, const std::vector& globalcols, const std::vector& rowvalues) + { + PetscErrorCode ierr; PetscInt nlr = nlocalrows(); - PetscInt nc = nglobalcols(); - double* a=NULL; - Mat lMat; + PetscInt nc = nglobalcols(); + double* a = NULL; + Mat lMat; PetscInt localrow = lri(globalrow); ierr = MatDenseGetLocalMatrix(mat(), &lMat); CHKERR(ierr); ierr = MatDenseGetArray(lMat, &a); CHKERR(ierr); - for(PetscInt j=0; j1000)k = 1000; - else if (n>100)k = 100; - else if (n>10)k = 10; + if (n > 1000)k = 1000; + else if (n > 100)k = 100; + else if (n > 10)k = 10; else k = 1; - if (n%k == 0){ + if (n % k == 0) { cPetscDistShellMatrix* pA = (cPetscDistShellMatrix*)context; - if (pA->mpirank() == 0){ + if (pA->mpirank() == 0) { //printf("CG Iteration %6d residual norm = %8.6le\n", n, rnorm); } if (n == 0) pA->_rnorm_start = rnorm; @@ -1611,13 +1617,13 @@ class cPetscDistShellMatrix : public cPetscObject { } cPetscDistVector solve_CG(const cPetscDistMatrix& P, const cPetscDistVector& b, const cPetscDistVector& initialguess = cPetscDistVector()) - { + { //A x = b cPetscDistVector x;//solution vector that will be returned - if (initialguess.globalsize() == 0){ - x.create("x",mpicomm(), nlocalcols(), nglobalcols(), 0.0); + if (initialguess.globalsize() == 0) { + x.create("x", mpicomm(), nlocalcols(), nglobalcols(), 0.0); } - else{ + else { x = initialguess; } @@ -1652,38 +1658,38 @@ class cPetscDistShellMatrix : public cPetscObject { //ierr = KSPSetTolerances(ksp, PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT, maxits); CHKERR(ierr); //Point to function that monitors/outputs convergence progress at each CG iteration - ierr = KSPMonitorSet(ksp, kspmonitor, this, PETSC_NULL); CHKERR(ierr); + ierr = KSPMonitorSet(ksp, kspmonitor, this, PETSC_NULLPTR); CHKERR(ierr); //See also PETSC functions KSPMonitorDefault, KSPMonitorTrueResidualNorm, KSPMonitorTrueResidualMaxNorm - //ierr = KSPMonitorSet(ksp, KSPMonitorDefault, PETSC_NULL, PETSC_NULL); CHKERR(ierr); - + //ierr = KSPMonitorSet(ksp, KSPMonitorDefault, PETSC_NULLPTR, PETSC_NULLPTR); CHKERR(ierr); + //Setup and run the Conjugate Gradients algorithm ierr = KSPSetUp(ksp); CHKERR(ierr); double t1 = gettime(); ierr = KSPSolve(ksp, b.vec(), x.vec()); CHKERR(ierr); double t2 = gettime(); - + ierr = KSPGetIterationNumber(ksp, &_niterations); CHKERR(ierr); ierr = KSPGetResidualNorm(ksp, &_rnorm_end); CHKERR(ierr); ierr = KSPGetConvergedReason(ksp, &_conv_reason_code); CHKERR(ierr); - + _solve_time = t2 - t1; _conv_reason_str = std::string(KSPConvergedReasons[_conv_reason_code]); - if (_conv_reason_code <= 0){ + if (_conv_reason_code <= 0) { _converged = false; } - else{ + else { _converged = true; } - ierr = KSPDestroy(&ksp); + ierr = KSPDestroy(&ksp); return x; } - - std::string convergence_summary(){ + + std::string convergence_summary() { std::string s; if (_converged) s += strprint("CG Converged\n"); - else s += strprint("CG Diverged\n"); + else s += strprint("CG Diverged\n"); s += strprint("\tCG Iterations=%d\n", _niterations); s += strprint("\tReason = %d (%s)\n", _conv_reason_code, _conv_reason_str.c_str()); s += strprint("\tResidual norm at start = %8.6le\n", _rnorm_start); @@ -1695,11 +1701,11 @@ class cPetscDistShellMatrix : public cPetscObject { } //Operators - cPetscDistVector operator*(const cPetscDistVector& x) const - { - cPetscDistVector b(mpicomm(), nlocalcols(), nglobalcols()); - vecmult(x, b); - return b; + cPetscDistVector operator*(const cPetscDistVector& x) const + { + cPetscDistVector b(mpicomm(), nlocalcols(), nglobalcols()); + vecmult(x, b); + return b; } }; diff --git a/visualstudio/galeiallatonce/galeiallatonce.vcxproj b/visualstudio/galeiallatonce/galeiallatonce.vcxproj index fd272f0..7103b6f 100644 --- a/visualstudio/galeiallatonce/galeiallatonce.vcxproj +++ b/visualstudio/galeiallatonce/galeiallatonce.vcxproj @@ -97,6 +97,7 @@ + @@ -104,7 +105,6 @@ - diff --git a/visualstudio/galeiallatonce/galeiallatonce.vcxproj.filters b/visualstudio/galeiallatonce/galeiallatonce.vcxproj.filters index 87fe14c..177c6e6 100644 --- a/visualstudio/galeiallatonce/galeiallatonce.vcxproj.filters +++ b/visualstudio/galeiallatonce/galeiallatonce.vcxproj.filters @@ -34,9 +34,6 @@ Header Files - - Header Files - Header Files @@ -52,5 +49,8 @@ Header Files + + Header Files + \ No newline at end of file