123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482148314841485148614871488148914901491149214931494149514961497149814991500150115021503150415051506150715081509151015111512151315141515151615171518151915201521152215231524152515261527152815291530153115321533153415351536153715381539154015411542154315441545154615471548154915501551155215531554155515561557155815591560156115621563156415651566156715681569157015711572157315741575157615771578157915801581158215831584158515861587158815891590159115921593159415951596159715981599160016011602160316041605160616071608160916101611161216131614161516161617161816191620162116221623162416251626162716281629163016311632163316341635163616371638163916401641164216431644164516461647164816491650165116521653165416551656165716581659166016611662166316641665166616671668166916701671167216731674167516761677167816791680168116821683168416851686168716881689169016911692169316941695169616971698169917001701170217031704170517061707170817091710171117121713171417151716171717181719172017211722172317241725172617271728172917301731173217331734173517361737173817391740174117421743174417451746174717481749175017511752175317541755175617571758175917601761176217631764176517661767176817691770177117721773177417751776177717781779178017811782178317841785178617871788178917901791179217931794179517961797179817991800180118021803180418051806180718081809181018111812181318141815181618171818181918201821182218231824182518261827182818291830183118321833183418351836183718381839184018411842184318441845184618471848184918501851185218531854185518561857185818591860186118621863186418651866186718681869187018711872187318741875187618771878187918801881188218831884188518861887188818891890189118921893189418951896189718981899190019011902190319041905190619071908190919101911191219131914191519161917191819191920192119221923192419251926192719281929193019311932193319341935193619371938193919401941194219431944194519461947194819491950195119521953195419551956 |
-
-
-
-
-
-
- #ifndef LOGBARRIERCG_INC
- #define LOGBARRIERCG_INC
-
- #include <Eigen/IterativeLinearSolvers>
- #include <iostream>
- #include <iomanip>
- #include <fstream>
- #include <Eigen/Eigen>
- #include "cg.h"
- #include "bicgstab.h"
- #include "lemma.h"
-
-
- namespace Lemma {
-
- template < typename Scalar >
- Scalar PhiB (const Scalar& mux, const Scalar& muy, const Scalar& minVal,
- const Scalar& maxVal, const VectorXr x) {
- Scalar phib = std::abs((x.array() - minVal).log().sum()*mux);
- phib += std::abs((maxVal - x.array()).log().sum()*muy);
- return phib;
- }
-
-
-
- template < typename Scalar >
- Scalar PhiB2 (const Scalar& mux, const Scalar& muy, const Scalar& minVal,
- const Scalar& maxVal, const VectorXr x, const int& block,
- const int &nblocks) {
- Scalar phib = std::abs((x.array() - minVal).log().sum()*mux);
-
- for (int ib=0; ib<nblocks; ++ib) {
-
- phib += Scalar(block)*std::log(maxVal - x.segment(ib*block, block).sum())*muy;
- }
- return phib;
- }
-
- template < typename Scalar >
- Scalar PhiB2 (const Scalar& minVal, const Scalar& maxVal, const VectorXr x,
- const int& block, const int &nblocks) {
- Scalar phib = std::abs((x.array() - minVal).log().sum());
-
- for (int ib=0; ib<nblocks; ++ib) {
-
- phib += Scalar(block)*std::log(maxVal - x.segment(ib*block, block).sum());
- }
- return phib;
- }
-
- template < typename Scalar >
- Scalar PhiB2_NN (const Scalar& mux, const Scalar& minVal, const VectorXr x) {
- Scalar phib = std::abs((x.array() - minVal).log().sum()*mux);
- return phib;
- }
-
-
-
-
- template < typename Scalar >
- VectorXr LogBarrierCG(const MatrixXr &A, const VectorXr &b,
- const VectorXr &x0,
- const Eigen::SparseMatrix<Scalar>& WdTWd,
- const Eigen::SparseMatrix<Scalar>& WmTWm,
- const Scalar &minVal,
- const Scalar &maxVal) {
-
-
- if (A.rows() != b.size() ) {
- std::cerr << "Matrix A is not aligned with Vector b" << "\n";
- exit(1);
- }
-
-
-
-
-
-
-
- MatrixXr ATWdTWdA = A.transpose()*WdTWd*A;
-
-
-
- Eigen::SparseMatrix<Scalar> MM (ATWdTWdA.rows(), ATWdTWdA.cols());
- for (int i=0; i<ATWdTWdA.rows(); ++i) {
- MM.insert(i,i) = 1. / ATWdTWdA(i,i);
- }
- MM.finalize();
-
- int N = A.cols();
- int M = A.rows();
- int MAXITER = 100;
- Scalar SIGMA = 1e-2;
-
-
-
-
- Scalar limsup = 1e10;
- for (int i=0; i<N; ++i) {
- VectorXr Spike = VectorXr::Zero(N);
- Spike(i) = (minVal + maxVal) / 2.;
- limsup = std::min(limsup, (ATWdTWdA*Spike).array().abs().maxCoeff());
- }
- Scalar lambda = limsup;
- Scalar mux0 = 1e-1*lambda;
- Scalar muy0 = 1e-1*lambda;
- Scalar epsilon = 1e-20;
-
-
- VectorXr x=x0;
-
-
- VectorXr b_pre = A*x;
-
- Scalar phid = (b - (b_pre)).norm();
- Scalar phim = x.norm();
- Scalar phib = PhiB(mux0, muy0, minVal, maxVal, x);
- Scalar mux = mux0;
- Scalar muy = muy0;
- Scalar tol = 1e-5*phid;
- std::fstream phireport("phimphid.dat", std::ios::out);
-
-
-
- for (int ii=0; ii<MAXITER; ++ii) {
-
-
-
-
-
-
- VectorXr Xm1 = x;
- int iter = N*M;
- mux = mux0;
- muy = muy0;
- int iloop(0);
- do {
-
-
-
- VectorXr X1 = VectorXr::Ones(N).array() / (x.array()-minVal) ;
- VectorXr Y1 = VectorXr::Ones(N).array() / (maxVal-x.array()) ;
-
- VectorXr X2 = VectorXr::Ones(N).array() / ((x.array()-minVal)*(x.array()-minVal));
- VectorXr Y2 = VectorXr::Ones(N).array() / ((maxVal-x.array())*(maxVal-x.array()));
-
-
-
-
- tol = 1e-5*phid;
- iter = N*M;
-
-
-
- VectorXr b2 = (-A.transpose()*WdTWd*(b_pre-b)).array() -
- (lambda*WmTWm*(x0-x0)).array() +
- 2.*mux*X1.array() + 2.*muy*Y1.array();
-
-
- MatrixXr A2 = ATWdTWdA;
- A2 += lambda*WmTWm;
- A2.diagonal() += mux*X2 + muy*Y2;
-
- VectorXr ztilde = CG(A2, VectorXr::Zero(N), b2, MM, iter, tol);
-
-
- VectorXr h = ztilde;
- Scalar d = std::min(1., 0.995*(x.array()/h.array().abs()).minCoeff() );
- x += d*h;
-
-
- VectorXr s1 = mux * (X2.asDiagonal()*ztilde - 2.*X1);
- VectorXr s2 = muy * (Y2.asDiagonal()*ztilde - 2.*Y1);
-
-
- mux = SIGMA/((Scalar)(N)) * std::abs( s1.dot(x) ) ;
- muy = SIGMA/((Scalar)(N)) * std::abs( s2.dot(x) ) ;
-
- b_pre = A*x;
- phid = (WdTWd*(b-b_pre)).norm() ;
- phim = (WmTWm*(x-x0)).norm();
- phib = PhiB(mux, muy, minVal, maxVal, x);
- ++iloop;
-
- } while (std::abs(phib / (phid+lambda*phim)) > epsilon);
-
-
- phireport.precision(12);
- std::cout << ii << "\t" << x.norm() << "\t" << (b-(A*x)).norm()
- << "\t" << lambda << "\t" << mux << "\t" << muy << "\t"
- << iter << "\t" << iloop << std::endl;
- phireport << ii << "\t" << x.norm() << "\t" << (b-(A*x)).norm()
- << "\t" << lambda << "\t" << mux << "\t" << muy << "\t"
- << iter << "\t" << iloop << std::endl;
-
- std::fstream modfile;
- std::string fname = "iteration" + to_string(ii) + ".dat";
- modfile.open( fname.c_str(), std::ios::out);
- modfile << ii << "\t" << x.norm() << "\t" << (b-(A*x)).norm()
- << "\t" << lambda << "\t" << mux << "\t" << muy << "\t" << iter << "\n";
- modfile << x << "\n";
- modfile.close();
-
-
-
- lambda *= .9;
-
- }
- phireport.close();
-
-
- return x;
- }
-
-
-
-
- template < typename Scalar >
- VectorXr LogBarrierCG(const MatrixXr &A, const VectorXr &b, const Scalar &minVal,
- const Scalar &maxVal) {
-
-
- if (A.rows() != b.size() ) {
- std::cerr << "Matrix A is not aligned with Vector b" << "\n";
- exit(1);
- }
-
-
- MatrixXr ATA = A.transpose()*A;
- int N = A.cols();
- int M = A.rows();
- int MAXITER = 100;
- Scalar SIGMA = 1e-2;
- Scalar delta = 1e-4;
-
-
-
-
-
-
-
- Scalar limsup = 1e10;
- for (int i=0; i<N; ++i) {
- VectorXr Spike = VectorXr::Zero(N);
- Spike(i) = (minVal + maxVal) / 2.;
- limsup = std::min(limsup, (ATA*Spike).array().abs().maxCoeff());
- }
- Scalar lambda = 1e3*limsup;
- Scalar mux0 = 1e-1*lambda;
- Scalar muy0 = 1e-1*lambda;
- Scalar epsilon = 1e-20;
-
-
-
-
-
-
-
-
-
-
-
-
- VectorXr x = VectorXr::Zero(N).array() + minVal+delta;
-
-
- VectorXr b_pre = A*x;
-
- Scalar phid = (b - (b_pre)).norm();
- Scalar phim = x.norm();
- Scalar phib = PhiB(mux0, muy0, minVal, maxVal, x);
- Scalar mux = mux0;
- Scalar muy = muy0;
- Scalar tol = 1e-5*phid;
- std::fstream phireport("phimphid.dat", std::ios::out);
-
-
-
- for (int ii=0; ii<MAXITER; ++ii) {
-
-
-
-
-
- VectorXr WmT_Wm = VectorXr::Ones(N).array()*lambda;
-
- VectorXr Xm1 = x;
- int iter = N*M;
- mux = mux0;
- muy = muy0;
- int iloop(0);
-
- do {
-
- VectorXr X1(x.size());
- VectorXr Y1(x.size());
- VectorXr X2(x.size());
- VectorXr Y2(x.size());
- VectorXr b2;
-
- MatrixXr A2;
-
-
-
- #ifdef LEMMAUSEOMP
- #pragma omp parallel sections
- {
- #endif
-
-
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- X1 = VectorXr::Ones(N).array() / (x.array()-minVal) ;
- }
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- Y1 = VectorXr::Ones(N).array() / (maxVal-x.array()) ;
- }
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- X2 = VectorXr::Ones(N).array() / ((x.array()-minVal)*(x.array()-minVal));
- }
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- Y2 = VectorXr::Ones(N).array() / ((maxVal-x.array())*(maxVal-x.array()));
- }
- #ifdef LEMMAUSEOMP
- }
- #endif
-
- #ifdef LEMMAUSEOMP
- #pragma omp parallel sections
- {
- #endif
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- b2 = -(A.transpose()*(b_pre-b)).array() -
- (WmT_Wm.array()*x.array()).array() +
- (2.*mux)*X1.array() + (2.*muy)*Y1.array();
- }
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- A2 = ATA;
- A2.diagonal().array() += WmT_Wm.array() + mux*X2.array() +
- muy*Y2.array();
- }
- #ifdef LEMMAUSEOMP
- }
- #endif
-
-
-
-
-
-
-
-
-
-
-
-
-
- tol = 1e-5*phid+mux+muy;
- iter = N*M;
-
-
-
-
-
-
-
-
-
-
-
-
-
- VectorXr ztilde = CG(A2, VectorXr::Zero(N), b2, iter, tol);
-
-
-
-
-
- VectorXr s1, s2;
-
-
-
-
-
- Scalar d = std::min(1.,0.995*(x.array()/ztilde.array().abs()).minCoeff());
- x += d*ztilde;
-
-
-
-
-
-
-
-
-
-
-
-
- b_pre = A*x;
- #ifdef LEMMAUSEOMP
- #pragma omp parallel sections
- {
- #endif
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- phib = PhiB(mux, muy, minVal, maxVal, x);
- }
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- phid = (b - (b_pre)).norm();
- }
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- phim = x.norm();
- }
- #ifdef LEMMAUSEOMP
- }
- #endif
-
-
- #ifdef LEMMAUSEOMP
- #pragma omp parallel sections
- {
- #endif
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- s1 = mux * (X2.asDiagonal()*(x+ztilde) - 2.*X1);
- mux = SIGMA/((Scalar)(N)) * std::abs( s1.dot(x) ) ;
- }
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- s2 = muy * (Y2.asDiagonal()*(x+ztilde) - 2.*Y1);
- muy = SIGMA/((Scalar)(N)) * std::abs( s2.dot(x) ) ;
- }
- #ifdef LEMMAUSEOMP
- }
- #endif
- ++iloop;
- } while (std::abs(phib / (phid+lambda*phim)) > epsilon);
-
-
- phireport.precision(12);
- std::cout << ii << "\t" << x.norm() << "\t" << (b-(A*x)).norm()
- << "\t" << lambda << "\t" << mux << "\t" << muy << "\t"
- << iter << "\t" << iloop << std::endl;
- phireport << ii << "\t" << x.norm() << "\t" << (b-(A*x)).norm()
- << "\t" << lambda << "\t" << mux << "\t" << muy << "\t"
- << iter << "\t" << iloop << std::endl;
-
-
- std::fstream modfile;
- std::string fname = "iteration" + to_string(ii) + ".dat";
- modfile.open( fname.c_str(), std::ios::out);
- modfile << ii << "\t" << x.norm() << "\t" << (b-(A*x)).norm()
- << "\t" << lambda << "\t" << mux << "\t" << muy << "\t" << iter << "\n";
- modfile << x << "\n";
- modfile.close();
-
-
- std::fstream predata;
- fname = "iteration" + to_string(ii) + "pre.dat";
- predata.open(fname.c_str(), std::ios::out);
- predata << b_pre << std::endl;
- predata.close();
-
-
-
- lambda *= .92;
-
- }
- phireport.close();
-
- return x;
- }
-
-
-
- template <typename Scalar>
- VectorXr LogBarrierCG(const MatrixXr &A, const VectorXr &b, const Scalar &minVal,
- const Scalar &maxVal, const int& block) {
-
-
- if (A.rows() != b.size() ) {
- std::cerr << "Matrix A is not aligned with Vector b" << "\n";
- exit(1);
- }
-
-
- std::fstream obsdata;
- std::string fname = "obsdata.dat";
- obsdata.open(fname.c_str(), std::ios::out);
- obsdata << b << std::endl;
- obsdata.close();
-
-
-
-
-
-
-
-
- MatrixXr ATA = A.transpose()*A;
- int N = A.cols();
- int M = A.rows();
- int MAXITER = 100;
- Scalar SIGMA = 1e-2;
- Scalar delta = 1e-4;
-
-
-
-
-
-
-
- Scalar limsup = 1e10;
- for (int i=0; i<N; ++i) {
- VectorXr Spike = VectorXr::Zero(N);
- Spike(i) = (minVal + maxVal) / 2.;
- limsup = std::min(limsup, (ATA*Spike).array().abs().maxCoeff());
- }
- Scalar lambda = 1e-6;
- Scalar mux0 = 1e-1*lambda;
- Scalar muy0 = 1e-1*lambda;
- Scalar epsilon = 1e-25;
-
-
-
-
-
-
-
-
-
-
-
-
- VectorXr x = VectorXr::Zero(N).array() + minVal+delta;
-
-
- VectorXr b_pre = A*x;
-
- int nblocks = x.size()/block;
- Scalar phid = (b - (b_pre)).norm();
- Scalar phim = x.norm();
- Scalar phib = PhiB2(mux0, muy0, minVal, maxVal, x, block, nblocks);
- Scalar mux = mux0;
- Scalar muy = muy0;
- Scalar tol = 1e-5*phid;
- std::fstream phireport("phimphid.dat", std::ios::out);
-
-
-
- for (int ii=0; ii<MAXITER; ++ii) {
-
-
-
-
-
- VectorXr WmT_Wm = VectorXr::Ones(N).array()*lambda;
-
- VectorXr Xm1 = x;
- int iter = N*M;
- mux = mux0;
- muy = muy0;
- int iloop(0);
-
-
-
-
-
-
-
-
- do {
-
- VectorXr X1(x.size());
- VectorXr Y1(x.size());
- VectorXr X2(x.size());
- VectorXr Y2(x.size());
- VectorXr b2;
-
- MatrixXr A2;
-
-
-
- #ifdef LEMMAUSEOMP
- #pragma omp parallel sections
- {
- #endif
-
-
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- X1 = VectorXr::Ones(N).array() / (x.array()-minVal) ;
- }
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- for (int ib=0; ib<nblocks; ++ib) {
-
- Y1.segment(ib*block, block) = VectorXr::Ones(block).array() /
- (maxVal - x.segment(ib*block, block).sum());
- }
-
-
-
-
- }
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- X2 = VectorXr::Ones(N).array() / ((x.array()-minVal)*(x.array()-minVal));
- }
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- for (int ib=0; ib<nblocks; ++ib) {
-
- Y2.segment(ib*block, block) = VectorXr::Ones(block).array() /
- ( (maxVal-x.segment(ib*block, block).sum()) *
- (maxVal-x.segment(ib*block, block).sum()) );
- }
-
-
-
- }
- #ifdef LEMMAUSEOMP
- }
- #endif
-
- #ifdef LEMMAUSEOMP
- #pragma omp parallel sections
- {
- #endif
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- b2 = -(A.transpose()*(b_pre-b)).array() -
- (WmT_Wm.array()*x.array()).array() +
- (2.*mux)*X1.array() + (2.*muy)*Y1.array();
- }
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- A2 = ATA;
- A2.diagonal().array() += WmT_Wm.array() + mux*X2.array() +
- muy*Y2.array();
- }
- #ifdef LEMMAUSEOMP
- }
- #endif
-
-
-
-
-
-
-
-
-
-
-
-
-
- tol = 1e-5*phid+mux+muy;
- iter = N*M;
-
-
-
-
-
-
-
-
-
-
-
-
-
- VectorXr ztilde = CG(A2, VectorXr::Zero(N), b2, iter, tol);
-
-
-
-
-
- VectorXr s1, s2;
-
-
-
-
-
- Scalar d = std::min(1.,0.995*(x.array()/ztilde.array().abs()).minCoeff());
- x += d*ztilde;
-
-
-
-
-
-
-
-
-
-
-
-
- b_pre = A*x;
-
- #ifdef LEMMAUSEOMP
- #pragma omp parallel sections
- {
- #endif
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- phib = PhiB2(mux, muy, minVal, maxVal, x, block, nblocks);
- }
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- phid = (b - (b_pre)).norm();
- }
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- phim = x.norm();
- }
- #ifdef LEMMAUSEOMP
- }
- #endif
-
-
- #ifdef LEMMAUSEOMP
- #pragma omp parallel sections
- {
- #endif
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- s1 = mux * (X2.asDiagonal()*(x+ztilde) - 2.*X1);
- mux = SIGMA/((Scalar)(N)) * std::abs( s1.dot(x) ) ;
- }
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- s2 = muy * (Y2.asDiagonal()*(x+ztilde) - 2.*Y1);
- muy = SIGMA/((Scalar)(N)) * std::abs( s2.dot(x) ) ;
- }
- #ifdef LEMMAUSEOMP
- }
- #endif
- ++iloop;
- } while (std::abs(phib / (phid+lambda*phim)) > epsilon);
-
-
- phireport.precision(12);
- std::cout << ii << "\t" << x.norm() << "\t" << (b-(A*x)).norm()
- << "\t" << lambda << "\t" << mux << "\t" << muy << "\t"
- << iter << "\t" << iloop << std::endl;
- phireport << ii << "\t" << x.norm() << "\t" << (b-(A*x)).norm()
- << "\t" << lambda << "\t" << mux << "\t" << muy << "\t"
- << iter << "\t" << iloop << std::endl;
-
- std::fstream modfile;
- std::string fname = "iteration" + to_string(ii) + ".dat";
- modfile.open( fname.c_str(), std::ios::out);
- modfile << ii << "\t" << x.norm() << "\t" << (b-(A*x)).norm()
- << "\t" << lambda << "\t" << mux << "\t" << muy << "\t" << iter << "\n";
- modfile << x << "\n";
- modfile.close();
-
-
- std::fstream predata;
- fname = "iteration" + to_string(ii) + "pre.dat";
- predata.open(fname.c_str(), std::ios::out);
- predata << b_pre << std::endl;
- predata.close();
-
-
-
- lambda *= .85;
-
- }
- phireport.close();
-
- return x;
- }
-
-
-
- template <typename Scalar>
- VectorXr LogBarrierCG(const MatrixXr &A, const VectorXr &xr,
- const VectorXr &b, const Scalar &minVal,
- const Scalar &maxVal, const int& block,
- const Eigen::SparseMatrix<Scalar>& WdTWd,
- const Eigen::SparseMatrix<Scalar>& WmTWm, Real lambda0=1e1) {
-
-
- if (A.rows() != b.size() ) {
- std::cerr << "Matrix A is not aligned with Vector b" << "\n";
- std::cerr << "A.rows() " << A.rows() << "\n";
- std::cerr << "A.cols() " << A.cols() << "\n";
- std::cerr << "b.size() " << b.size() << "\n";
- exit(1);
- }
-
-
- std::fstream obsdata;
- std::string fname = "obsdata.dat";
- obsdata.open(fname.c_str(), std::ios::out);
- obsdata << b << std::endl;
- obsdata.close();
-
-
- MatrixXr ATWdTWdA = A.transpose()*WdTWd*A;
- int N = A.cols();
- int M = A.rows();
- int MAXITER = 175;
- Scalar SIGMA = .25;
- Scalar delta = 1e-4;
-
-
-
-
-
-
-
-
-
- Scalar lambda = lambda0;
- Scalar epsilon = 1e-15;
-
-
- VectorXr x = VectorXr::Zero(N).array() + minVal+delta;
-
-
- VectorXr b_pre = A*x;
-
- int nblocks = x.size()/block;
- Scalar phid = (WdTWd*(b - b_pre)).norm();
- Scalar phim = (WmTWm*(x - xr)).norm();
- Scalar phib = PhiB2(minVal, maxVal, x, block, nblocks);
- Scalar mux = (phid + lambda*phim) / phib;
- Scalar muy = mux;
-
- std::fstream phireport("phimphid.dat", std::ios::out);
-
- Eigen::ConjugateGradient< MatrixXr > cg;
-
-
- for (int ii=0; ii<MAXITER; ++ii) {
-
- int iter = N*M;
- mux = (phid + lambda*phim) / phib;
- muy = mux;
- int iloop(0);
- int itertot(0);
- VectorXr h;
- bool cont(true);
- do {
-
-
-
- VectorXr X1(x.size());
- VectorXr Y1(x.size());
- VectorXr X2(x.size());
- VectorXr Y2(x.size());
- VectorXr b2;
- MatrixXr A2;
-
-
-
-
-
-
-
- X1 = VectorXr::Ones(N).array() / (x.array()-minVal) ;
-
- for (int ib=0; ib<nblocks; ++ib) {
- Y1.segment(ib*block, block) = VectorXr::Ones(block).array() /
- (maxVal - x.segment(ib*block, block).sum());
- }
-
- X2 = VectorXr::Ones(N).array() / ((x.array()-minVal)*(x.array()-minVal));
-
- for (int ib=0; ib<nblocks; ++ib) {
- Y2.segment(ib*block, block) = VectorXr::Ones(block).array() /
- ( (maxVal-x.segment(ib*block, block).sum()) *
- (maxVal-x.segment(ib*block, block).sum()) );
- }
-
-
-
-
-
-
-
- b2 = (A.transpose()*WdTWd*(b)).array()
-
- + (2.*mux)*X1.array() + (2.*muy)*Y1.array();
-
- A2 = ATWdTWdA;
- A2 += lambda*WmTWm;
- A2.diagonal().array() += mux*X2.array() + muy*Y2.array();
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- iter = N*M;
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- cg.compute(A2);
- VectorXr ztilde = cg.solveWithGuess(b2, x);
- iter = cg.iterations();
-
-
- ++iloop;
- itertot += iter;
-
-
-
-
-
-
- h = ztilde - x;
-
-
-
- Scalar d(1.);
- for (int ix=0; ix<x.size(); ++ix) {
- if (h[ix] < 0.) {
- d = std::min(d, (Scalar).925*(x[ix]/std::abs(h[ix])));
- }
- }
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- x += d*h;
-
-
- for (int ib=0; ib<nblocks; ++ib) {
- while (x.segment(ib*block, block).sum() >= maxVal) {
- x.segment(ib*block, block).array() *= .99;
- }
- }
-
- for (int i=0; i<x.size(); ++i) {
- if (x(i) < minVal) {
- x(i) = minVal + delta;
- }
- }
-
- b_pre = A*x;
-
- phib = PhiB2(mux, muy, minVal, maxVal, x, block, nblocks);
- phid = (WdTWd*(b-b_pre)).norm();
- phim = (WmTWm*(x-xr)).norm();
-
-
- VectorXr s1 = mux * (X2.asDiagonal()*(ztilde) - 2.*X1);
- mux = SIGMA/((Scalar)(N)) * std::abs( s1.dot(x) ) ;
- VectorXr s2 = muy * (Y2.asDiagonal()*(ztilde) - 2.*Y1);
- muy = SIGMA/((Scalar)(N)) * std::abs( s2.dot(x) ) ;
-
- if ( (std::abs(phib / (phid+lambda*phim)) < epsilon)) {
-
- cont = false;
- }
-
- } while ( cont );
-
-
-
-
- std::cout << std::setw(6) << ii << std::scientific << std::setw(18) << phim << std::setw(18) << phid
- << std::setw(18) << lambda << std::setw(18) << mux << std::setw(18) << muy
- << std::setw(12) << itertot << std::setw(12) << iloop << std::setw(18) << h.norm() << std::endl;
-
- phireport.precision(12);
- phireport << ii << "\t" << phim << "\t" << phid
- << "\t" << lambda << "\t" << mux << "\t" << muy << "\t"
- << itertot << "\t" << iloop << "\t" << h.norm() << std::endl;
-
- std::fstream modfile;
- std::string fname = "iteration" + to_string(ii) + ".dat";
- modfile.open( fname.c_str(), std::ios::out);
- modfile << ii << "\t" << phim << "\t" << phid
- << "\t" << lambda << "\t" << mux << "\t" << muy << "\t" << iter << "\n";
- modfile << x << "\n";
- modfile.close();
-
-
- std::fstream predata;
- fname = "iteration" + to_string(ii) + "pre.dat";
- predata.open(fname.c_str(), std::ios::out);
- predata << b_pre << std::endl;
- predata.close();
-
-
-
- lambda *= .9;
-
- }
- phireport.close();
-
- return x;
- }
-
-
-
-
- template <typename Scalar>
- VectorXr LogBarrierCG_NN(const MatrixXr &A, const VectorXr &xr,
- const VectorXr &b, const Scalar &minVal,
- const Eigen::SparseMatrix<Scalar>& WdTWd,
- const Eigen::SparseMatrix<Scalar>& WmTWm, Real lambda0=1e1) {
-
-
- if (A.rows() != b.size() ) {
- std::cerr << "Matrix A is not aligned with Vector b" << "\n";
- std::cerr << "A.rows() " << A.rows() << "\n";
- std::cerr << "A.cols() " << A.cols() << "\n";
- std::cerr << "b.size() " << b.size() << "\n";
- exit(1);
- }
-
-
- std::fstream obsdata;
- std::string fname = "obsdata.dat";
- obsdata.open(fname.c_str(), std::ios::out);
- obsdata << b << std::endl;
- obsdata.close();
-
-
- MatrixXr ATWdTWdA = A.transpose()*WdTWd*A;
- int N = A.cols();
- int M = A.rows();
- int MAXITER = 175;
- Scalar SIGMA = .25;
- Scalar delta = 1e-12;
-
-
-
-
-
-
-
-
-
- Scalar lambda = lambda0;
- Scalar epsilon = 1e-16;
-
-
- VectorXr x = VectorXr::Zero(N).array() + minVal+delta;
-
-
- VectorXr b_pre = A*x;
-
-
-
- Eigen::ConjugateGradient< MatrixXr, Eigen::Upper, Eigen::DiagonalPreconditioner<Real> > cg;
-
-
- Scalar phid = (WdTWd*(b - b_pre)).norm();
- Scalar phim = (WmTWm*(x - xr)).norm();
- Scalar phib = PhiB2_NN(1., minVal, x);
- Scalar mux = (phid + lambda*phim) / phib;
-
- std::fstream phireport("phimphid.dat", std::ios::out);
- mux = (phid + lambda*phim) / phib;
-
-
- for (int ii=0; ii<MAXITER; ++ii) {
-
- int iter = N*M;
- int iloop(0);
- int itertot(0);
-
- VectorXr ztilde;
- bool cont(true);
- do {
-
-
-
- VectorXr X1(x.size());
- VectorXr X2(x.size());
- VectorXr b2;
- MatrixXr A2;
-
-
-
-
-
-
-
- X1 = VectorXr::Ones(N).array() / (x.array()-minVal) ;
- X2 = VectorXr::Ones(N).array() / ((x.array()-minVal)*(x.array()-minVal));
-
-
- b2 = (A.transpose()*WdTWd*(b-b_pre)).array()
- - lambda*(WmTWm*(x-xr)).array()
- + (2.*mux)*X1.array();
-
-
- A2 = ATWdTWdA;
- A2 += lambda*WmTWm;
- A2.diagonal().array() += mux*X2.array();
-
-
-
-
-
-
- iter = N*M;
-
-
-
- cg.compute(A2);
-
- ztilde = cg.solve(b2);
-
- iter = cg.iterations();
-
-
- ++iloop;
- itertot += iter;
-
-
-
-
-
-
- Scalar d(1.);
- for (int ix=0; ix<x.size(); ++ix) {
- if (ztilde[ix] < 0.) {
- d = std::min(d, (Scalar).925*(x[ix]/std::abs(ztilde[ix])));
- }
- }
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- x += d*ztilde;
-
- for (int i=0; i<x.size(); ++i) {
- if (x(i) < minVal) {
- x(i) = minVal + delta;
- }
- }
-
- b_pre = A*x;
-
- phib = PhiB2_NN(mux, minVal, x);
- phid = (WdTWd*(b-b_pre)).norm();
- phim = (WmTWm*(x-xr)).norm();
-
-
- VectorXr s1 = mux * (X2.asDiagonal()*(ztilde) - 2.*X1);
- mux = SIGMA/((Scalar)(N)) * std::abs( s1.dot(x) ) ;
-
- if ( (std::abs(phib / (phid+lambda*phim)) < epsilon)) {
-
- cont = false;
- }
-
- } while ( cont );
-
-
-
-
- std::cout << std::setw(6) << ii << std::scientific << std::setw(18) << phim << std::setw(18) << phid
- << std::setw(18) << lambda << std::setw(18) << mux
- << std::setw(12) << itertot << std::setw(12) << iloop << std::setw(18) << ztilde.norm() << std::endl;
-
- phireport.precision(12);
- phireport << ii << "\t" << phim << "\t" << phid
- << "\t" << lambda << "\t" << mux << "\t"
- << itertot << "\t" << iloop << "\t" << ztilde.norm() << std::endl;
-
- std::fstream modfile;
- std::string fname = "iteration" + to_string(ii) + ".dat";
- modfile.open( fname.c_str(), std::ios::out);
- modfile << ii << "\t" << phim << "\t" << phid
- << "\t" << lambda << "\t" << mux << "\t" << iter << "\n";
- modfile << x << "\n";
- modfile.close();
-
-
- std::fstream predata;
- fname = "iteration" + to_string(ii) + "pre.dat";
- predata.open(fname.c_str(), std::ios::out);
- predata << b_pre << std::endl;
- predata.close();
-
-
-
- lambda *= .9;
-
- }
- phireport.close();
-
- return x;
- }
-
-
-
-
-
-
-
- template <typename Scalar>
- VectorXr Tikhonov_CG(const MatrixXr &A, const VectorXr &xr,
- const VectorXr &b,
- const Eigen::SparseMatrix<Scalar>& WdTWd,
- const Eigen::SparseMatrix<Scalar>& WmTWm, Real lambda0=1e1) {
-
-
-
- if (A.rows() != b.size() ) {
- std::cerr << "Matrix A is not aligned with Vector b" << "\n";
- std::cerr << "A.rows() " << A.rows() << "\n";
- std::cerr << "A.cols() " << A.cols() << "\n";
- std::cerr << "b.size() " << b.size() << "\n";
- exit(1);
- }
-
-
- std::fstream obsdata;
- std::string fname = "obsdata.dat";
- obsdata.open(fname.c_str(), std::ios::out);
- obsdata << b.array() << std::endl;
- obsdata.close();
-
-
-
- MatrixXr ATWdTWdA = A.transpose()*WdTWd*A;
- int N = A.cols();
- int M = A.rows();
- int MAXITER = 175;
-
-
- Scalar delta = 1e-4;
- Scalar lambda = lambda0;
- Scalar epsilon = 1e-15;
-
-
- VectorXr x = VectorXr::Zero(N).array() + delta;
-
-
- VectorXr b_pre = A*x;
- Scalar phid = (WdTWd*(b - b_pre)).norm();
- Scalar phim = (WmTWm*(x - xr)).norm();
- Scalar tol = 1e-5*phid;
- std::fstream phireport("phimphid.dat", std::ios::out);
-
-
-
- Eigen::ConjugateGradient< MatrixXr > cg;
-
-
- for (int ii=0; ii<MAXITER; ++ii) {
-
- int iter = N*M;
- int iloop(0);
- int itertot(0);
- VectorXr h;
- bool cont(true);
-
- do {
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- VectorXr b2 = A.transpose()*WdTWd*b;
-
- MatrixXr A2 = ATWdTWdA;
- A2 += lambda*WmTWm;
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- tol = 1e-5*phid;
- iter = N*M;
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- cg.compute(A2);
- x = cg.solveWithGuess(b2, x);
- b_pre = A*x;
- phid = (WdTWd*(b - b_pre)).norm();
- phim = (WmTWm*(x - xr)).norm();
-
-
- ++iloop;
- itertot += iter;
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- cont = false;
- } while ( cont );
-
-
- std::cout << std::ios::left;
- std::cout.precision(8);
- std::cout << std::setw(6) << ii << std::scientific << std::setw(18) << phim << std::setw(18) << phid
- << std::setw(18) << lambda << std::setw(18) << cg.error()
- << std::setw(12) << cg.iterations() << std::endl;
-
- phireport.precision(12);
- phireport << ii << "\t" << phim << "\t" << phid
- << "\t" << lambda << "\t" << 0.0 << "\t" << 0.0 << "\t"
- << cg.iterations() << "\t" << 0.0 << "\t" << 0.0 << std::endl;
-
- std::fstream modfile;
- std::string fname = "iteration" + to_string(ii) + ".dat";
- modfile.open( fname.c_str(), std::ios::out);
- modfile << ii << "\t" << phim << "\t" << phid
- << "\t" << lambda << "\t" << 0 << "\t" << 0 << "\t" << cg.iterations() << "\n";
- modfile << x << "\n";
- modfile.close();
-
-
- std::fstream predata;
- fname = "iteration" + to_string(ii) + "pre.dat";
- predata.open(fname.c_str(), std::ios::out);
- predata << b_pre << std::endl;
- predata.close();
-
-
-
- lambda *= .9;
-
- }
-
- phireport.close();
-
-
- return x;
- }
-
-
-
-
-
- template <typename Scalar>
- VectorXr LogBarrierCGLS(const MatrixXr &A, const VectorXr &b, const Scalar &minVal,
- const Scalar &maxVal, const int& block) {
-
-
- if (A.rows() != b.size() ) {
- std::cerr << "Matrix A is not aligned with Vector b" << "\n";
- exit(1);
- }
-
-
- std::fstream obsdata;
- std::string fname = "obsdata.dat";
- obsdata.open(fname.c_str(), std::ios::out);
- obsdata << b << std::endl;
- obsdata.close();
-
-
-
-
-
-
-
-
- MatrixXr ATA = A.transpose()*A;
- int N = A.cols();
- int M = A.rows();
- int MAXITER = 100;
- Scalar SIGMA = .25;
- Scalar delta = 1e-8;
-
-
-
-
-
-
-
- Scalar limsup = 1e10;
- for (int i=0; i<N; ++i) {
- VectorXr Spike = VectorXr::Zero(N);
- Spike(i) = (minVal + maxVal) / 2.;
- limsup = std::min(limsup, (ATA*Spike).array().abs().maxCoeff());
- }
- Scalar lambda = 1e-6;
- Scalar mux0 = 1e-1*lambda;
- Scalar muy0 = 1e-1*lambda;
- Scalar epsilon = 1e-5;
-
-
-
-
-
-
-
-
-
-
-
-
- VectorXr x = VectorXr::Zero(N).array() + minVal+delta;
-
-
- VectorXr b_pre = A*x;
-
- int nblocks = x.size()/block;
- Scalar phid = (b - (b_pre)).norm();
- Scalar phim = x.norm();
- Scalar phib = PhiB2(mux0, muy0, minVal, maxVal, x, block, nblocks);
- Scalar mux = mux0;
- Scalar muy = muy0;
- Scalar tol = 1e-5*phid;
- std::fstream phireport("phimphid.dat", std::ios::out);
-
-
-
- for (int ii=0; ii<MAXITER; ++ii) {
-
-
-
-
-
- VectorXr WmT_Wm = VectorXr::Ones(N).array()*lambda;
-
- VectorXr Xm1 = x;
- int iter = N*M;
- int iloop(0);
-
- do {
-
- VectorXr X1(x.size());
- VectorXr Y1(x.size());
- VectorXr X2(x.size());
- VectorXr Y2(x.size());
- VectorXr b2;
-
- MatrixXr A2;
-
-
-
- #ifdef LEMMAUSEOMP
- #pragma omp parallel sections
- {
- #endif
-
-
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- X1 = VectorXr::Ones(N).array() / (x.array()-minVal) ;
- }
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- for (int ib=0; ib<nblocks; ++ib) {
-
- Y1.segment(ib*block, block) = VectorXr::Ones(block).array() /
- (maxVal - x.segment(ib*block, block).sum());
- }
-
-
-
-
- }
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- X2 = VectorXr::Ones(N).array() / ((x.array()-minVal)*(x.array()-minVal));
- }
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- for (int ib=0; ib<nblocks; ++ib) {
-
- Y2.segment(ib*block, block) = VectorXr::Ones(block).array() /
- ( (maxVal-x.segment(ib*block, block).sum()) *
- (maxVal-x.segment(ib*block, block).sum()) );
- }
-
-
-
- }
- #ifdef LEMMAUSEOMP
- }
- #endif
-
- #ifdef LEMMAUSEOMP
- #pragma omp parallel sections
- {
- #endif
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- b2 = -(A.transpose()*(b_pre-b)).array() -
- (WmT_Wm.array()*x.array()).array() +
- (2.*mux)*X1.array() + (2.*muy)*Y1.array();
- }
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- A2 = ATA;
- A2.diagonal().array() += WmT_Wm.array() + mux*X2.array() +
- muy*Y2.array();
- }
- #ifdef LEMMAUSEOMP
- }
- #endif
-
-
-
-
-
-
-
-
-
-
-
-
-
- tol = 1e-5*phid+mux+muy;
- iter = N*M;
-
-
-
-
-
-
-
-
-
-
-
-
-
- VectorXr ztilde = CG(A2, VectorXr::Zero(N), b2, iter, tol);
-
-
-
-
-
- VectorXr s1, s2;
-
-
-
-
-
- Scalar d = std::min(1.,0.995*(x.array()/ztilde.array().abs()).minCoeff());
- x += d*ztilde;
-
-
-
-
-
-
-
-
-
-
-
-
- b_pre = A*x;
-
- #ifdef LEMMAUSEOMP
- #pragma omp parallel sections
- {
- #endif
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- phib = PhiB2(mux, muy, minVal, maxVal, x, block, nblocks);
- }
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- phid = (b - (b_pre)).norm();
- }
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- phim = x.norm();
- }
- #ifdef LEMMAUSEOMP
- }
- #endif
-
-
- #ifdef LEMMAUSEOMP
- #pragma omp parallel sections
- {
- #endif
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- s1 = mux * (X2.asDiagonal()*(x+ztilde) - 2.*X1);
- mux = SIGMA/((Scalar)(N)) * std::abs( s1.dot(x) ) ;
- }
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- s2 = muy * (Y2.asDiagonal()*(x+ztilde) - 2.*Y1);
- muy = SIGMA/((Scalar)(N)) * std::abs( s2.dot(x) ) ;
- }
- #ifdef LEMMAUSEOMP
- }
- #endif
- ++iloop;
-
- } while (std::abs(phib / (phid+lambda*phim)) > epsilon);
-
-
- phireport.precision(12);
- std::cout << ii << "\t" << x.norm() << "\t" << (b-(A*x)).norm()
- << "\t" << lambda << "\t" << mux << "\t" << muy << "\t"
- << iter << "\t" << iloop << std::endl;
- phireport << ii << "\t" << x.norm() << "\t" << (b-(A*x)).norm()
- << "\t" << lambda << "\t" << mux << "\t" << muy << "\t"
- << iter << "\t" << iloop << std::endl;
-
- std::fstream modfile;
- std::string fname = "iteration" + to_string(ii) + ".dat";
- modfile.open( fname.c_str(), std::ios::out);
- modfile << ii << "\t" << x.norm() << "\t" << (b-(A*x)).norm()
- << "\t" << lambda << "\t" << mux << "\t" << muy << "\t" << iter << "\n";
- modfile << x << "\n";
- modfile.close();
-
-
- std::fstream predata;
- fname = "iteration" + to_string(ii) + "pre.dat";
- predata.open(fname.c_str(), std::ios::out);
- predata << b_pre << std::endl;
- predata.close();
-
-
-
- lambda *= .85;
-
- }
- phireport.close();
-
- return x;
- }
-
-
- }
-
- #endif
|