GNU Gama, vyrovnání měření zprostředkujících s diagonální maticí vah

Součástí projektu GNU Gama je C++ knihovna tříd a funkcí. Třída Adj implementuje úlohu vyrovnání měření zprostředkujících

$\mathbf {A} \mathbf {x} =\mathbf {l} ,\qquad \mathbf {P}$

kde A je první matice plánu, x vektor určovaných parametrů, l vektor absolutních členů a P matice vah.

Pro řídkou matici A a diagonální váhovou matici P může zapsat vstupní data v následujícím formátu

  počet_neznámých   počet_měření

následuje popis jednotlivých řádků řídké soustavy

  N  index_1  index_2 ... index_N
  váha  pravá_strana  koef_1  koef_2  ... koef_N

Například soustavu

$\mathbf {A} ={\begin{pmatrix}1.0&0&0&0\\0.9&0&0.1&0\\2.0&0&3.0&0.1\\0.5&0.6&0.2&0.3\\-0.1&0&0.1&3.0\\0.6&-0.1&-0.9&0.4\end{pmatrix}},\qquad \mathbf {l} ={\begin{pmatrix}1.03\\1.18\\11.38\\3.51\\12.21\\-0.66\end{pmatrix}},\qquad \mathbf {P} =\mathrm {diag} {\begin{pmatrix}1\\1.1\\1.2\\1.3\\1.4\\1.5\end{pmatrix}}$

zapíšeme jako

   4  6

   1          1
   1   1.03   1
   2          1    3
  1.1  1.18  0.9  0.1
   3          1    3    4
  1.2 11.38  2.0  3.0  0.1   
   4          1    2    3    4
  1.3  3.51  0.5  0.6  0.2  0.3
   3          1    3    4
  1.4 12.21 -0.1  0.1  3.0
   4          1    2    3    4
  1.5 -0.66  0.6 -0.1 -0.9  0.4

počet mezer není významný, zde jsou vícenásobné mezery použity pouze pro formátování pro lepší čitelnost. Tento vstupní formát čte metoda read_gama_local_old_format(cin) třídy AdjInputData.

/*******************************************************************

Vstupni data:
------------

Zapisuji se pouze nenulove koeficienty rovnic oprav. Prvni zaznam
obsahuhue pouze dve cela cisla

   pocet_neznamych  pocet_mereni

dale pak nasleduji dva radky pro kazdou rovnici

   N index_1 index_2 ... index_N
   vaha prava_strana koef_1 koef_2 ... koef_N

kde N je pocet nenulovych prvku v danem radku


Volani programu:
---------------

   kombinace-01 [envelope|gso|svd|cholesky] < vstup.txt > vystup.xml

Implicitni algoritmus je reseni ridke matice v pametovem modelu
envelope.

********************************************************************/

#include <iostream>
#include <sstream>
#include <cmath>
#include <vector>
#include <gnu_gama/adj/adj.h>

using std::cin;
using std::cout;
using std::endl;
using std::cerr;
using std::sqrt;
using std::vector;
using std::string;


int main(int argc, char* argv[])
{
  using namespace GNU_gama;

  // objekt typu AdjInputData musi byt vytvoren dynamicky, protoze
  // pointer prejima objekt Adj, ktery pro vtupni data vola destruktor

  AdjInputData data;
  data.read_gama_local_old_format(cin);

  Adj adj;
  adj.set(&data);

  string alg = "envelope";
  if (argc == 2)
    {
      if      (alg == string("envelope")) adj.set_algorithm(Adj::envelope);
      else if (alg == string("gso")     ) adj.set_algorithm(Adj::gso);
      else if (alg == string("svd")     ) adj.set_algorithm(Adj::svd);
      else if (alg == string("cholesky")) adj.set_algorithm(Adj::cholesky);
      else
        {
          cerr <<
            "volani : kombinace-01 [envelope|gso|svd|cholesky] "
            "< vstup > vystup\n\n";
        }
    }

  Vec<> x = adj.x();
  Vec<> r = adj.r();

  cout << "# algoritmus = " << alg << endl;
  cout << "# nezname    = " << x.dim() << endl;
  cout << "# opravy     = " << r.dim() << endl;
  cout << "# defekt     = " << adj.defect() << endl;
  cout << "# [pvv]      = " << adj.rtr() << endl;
  cout << endl;

  double m02 = r.dim() > x.dim() ? adj.rtr()/(r.dim() - x.dim()) : 0.0;

  for (int i=1; i<=x.dim(); i++)
    {
      cout << i << "\t" << x(i) << "\t" << sqrt(m02*adj.q_xx(i,i)) << endl;
    }

  // cout << endl;
  //
  // for (int i=1; i<=r.dim(); i++)
  //   {
  //     cout << i << "\t" << r(i) << "\t" << sqrt(m02*adj.q_bb(i,i)) << endl;
  //   }
}