ExampleLapack.cpp

/* This file is part of the BIAS library (Basic ImageAlgorithmS).

Copyright (C) 2003-2009    (see file CONTACT for details)
  Multimediale Systeme der Informationsverarbeitung
  Institut fuer Informatik
  Christian-Albrechts-Universitaet Kiel


BIAS is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version.

BIAS is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License
along with BIAS; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/


/** @example ExampleLapack.cpp
    @relates SVD, Lapack
    @ingroup g_examples    
    @brief  example for Lapack MahalanobisDistance usage
    @author  koeser 11/2007 
*/

#include <MathAlgo/SVD.hh>
#include <Base/Math/Random.hh>
#include <Base/Debug/TimeMeasure.hh>
#include <MathAlgo/Lapack.hh>


using namespace BIAS;
using namespace std;


int main()
{
  Random rand;
  unsigned int matdim = 10;
  Matrix<double> Cov(matdim,matdim,MatrixIdentity), 
    factor(matdim,matdim,MatrixIdentity);
  BIAS::Vector<double> dif(matdim);
  for (unsigned int i=0; i<matdim; i++) {
    for (unsigned int j=0; j<matdim; j++) {
      factor[i][j] = rand.GetUniformDistributed(0.01,10);
    }
    dif[i] = rand.GetUniformDistributed(-10,10);
    factor[i][i] += rand.GetUniformDistributed(10, 50);
  }
  Cov = factor.Transpose()*factor;
  Cov.MakeSymmetric();
  cout<<"Cov is "<< Cov<<endl;
  
  SVD mytestsvd;


  double svddist=0,lapackdist=0;
  TimeMeasure tsvd, tsvd3;
  for (unsigned int i=0; i<1000; i++) {
    tsvd.Start();
    svddist = BIAS::Vector<double>(mytestsvd.Invert(Cov)*dif).ScalarProduct(dif);
    tsvd.Stop();
    tsvd3.Start();
    lapackdist = MahalanobisDistance(Cov,dif);
    tsvd3.Stop();
  }

  cerr << "svd: ";
  tsvd.Print();
  cerr << "user: "<<(double)tsvd.GetUserTime() << " ms";
  cerr << "real: "<<(double)tsvd.GetRealTime() << " us";
  cerr << "lapack: ";
  tsvd3.Print();
  cerr << "user: "<<(double)tsvd3.GetUserTime() << " ms";
  cerr << "real: "<<(double)tsvd3.GetRealTime() << " us"<<endl;

  cout<<"distances are "<<svddist<<" "<<lapackdist <<endl;
  return 0;


}