main()

int main	(	int	iargc,
		char *	argv[]
	)

Definition at line 87 of file lsqexample.cc.

References linear::op::dotNxM(), linear::lapack::dposv(), printmatrix(), TF_LSQEXAMPLE_CC_VERSION, and linear::op::transposeNxM().

{

  // size of matrix for test
  const int n=3;
  const int m=4;

  // set up random number generator
  tfxx::numeric::RNGgaussian rng;

  // create matrix of partial derivatives
  Tmatrix G(m,n);
  // create data vector
  Tmatrix d(m,1);
  // create vector for given model parameters
  Tmatrix mknown(n,1);

  // fill matrix with random numbers
  // enclose code in block to make declarations local
  {
    aff::Iterator<Tmatrix> I(G);
    while (I.valid())
    {
      *I = rng();
      ++I;
    }
  }

  // fill vector with random numbers
  // enclose code in block to make declarations local
  {
    aff::Iterator<Tmatrix> I(mknown);
    while (I.valid())
    {
      *I = rng();
      ++I;
    }
  }

  cout << TF_LSQEXAMPLE_CC_VERSION "\n" << endl;
  cout << "Solve least squares problem" << endl;
  cout << "===========================" << endl;

  cout << "\n"
          "Residuals:\n"
          "r = G * m - d\n"
          "Find solution such that |r|^2 becomes minimal" << endl;

  // setup data vector
  d=linear::op::dotNxM(G,mknown);

  // report system
  cout << "\nMatrix of partial derivatives (G):\n";
  printmatrix(G);
  cout << "\nData vector (d):\n";
  printmatrix(d);
 
  // setup system of linear equations
  // --------------------------------
  // system matrix
  Tmatrix M=linear::op::dotNxM(linear::op::transposeNxM(G),G);
  // right hand side
  Tmatrix rhs=linear::op::dotNxM(linear::op::transposeNxM(G),d);

  cout << "\nSystem of linear equations to be solved:\n"
          "M * m = rhs\n" << endl;

  cout << "With matrix M=G^T*G:\n";
  printmatrix(M);
  cout << "and vector rhs=G^T*d:\n";
  printmatrix(rhs);

  // solve system
  // ------------
  // solution vetcor (optimal model parameters)
  Tmatrix mopt=linear::lapack::dposv(M,rhs,'U',false);

  cout << "\nsolution vector mopt:\n";
  printmatrix(mopt);
  cout << "\nexpected solution (mknown):\n";
  printmatrix(mknown);
  cout << "\nresidual (mopt-mknown):\n";
  printmatrix(mopt-mknown);
}

Here is the call graph for this function: