text:
// g++ -O3 -DNDEBUG -I.. -L /usr/lib64/atlas/ benchBlasGemm.cpp -o benchBlasGemm -lrt -lcblas
// possible options:
//    -DEIGEN_DONT_VECTORIZE
//    -msse2

// #define EIGEN_DEFAULT_TO_ROW_MAJOR
#define _FLOAT

#include <iostream>

#include <Eigen/Core>
#include "BenchTimer.h"

// include the BLAS headers
extern "C" {
#include <cblas.h>
}
#include <string>

#ifdef _FLOAT
typedef float Scalar;
#define CBLAS_GEMM cblas_sgemm
#else
typedef double Scalar;
#define CBLAS_GEMM cblas_dgemm
#endif

typedef Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic> MyMatrix;
void bench_eigengemm(MyMatrix& mc, const MyMatrix& ma, const MyMatrix& mb, int nbloops);
void check_product(int M, int N, int K);
void check_product(void);

int main(int argc, char *argv[])
{
  // disable SSE exceptions
  #ifdef __GNUC__
  {
    int aux;
    asm(
    "stmxcsr   %[aux]           \n\t"
    "orl       $32832, %[aux]   \n\t"
    "ldmxcsr   %[aux]           \n\t"
    : : [aux] "m" (aux));
  }
  #endif

  int nbtries=1, nbloops=1, M, N, K;

  if (argc==2)
  {
    if (std::string(argv[1])=="check")
      check_product();
    else
      M = N = K = atoi(argv[1]);
  }
  else if ((argc==3) && (std::string(argv[1])=="auto"))
  {
    M = N = K = atoi(argv[2]);
    nbloops = 1000000000/(M*M*M);
    if (nbloops<1)
      nbloops = 1;
    nbtries = 6;
  }
  else if (argc==4)
  {
    M = N = K = atoi(argv[1]);
    nbloops = atoi(argv[2]);
    nbtries = atoi(argv[3]);
  }
  else if (argc==6)
  {
    M = atoi(argv[1]);
    N = atoi(argv[2]);
    K = atoi(argv[3]);
    nbloops = atoi(argv[4]);
    nbtries = atoi(argv[5]);
  }
  else
  {
    std::cout << "Usage: " << argv[0] << " size  \n";
    std::cout << "Usage: " << argv[0] << " auto size\n";
    std::cout << "Usage: " << argv[0] << " size nbloops nbtries\n";
    std::cout << "Usage: " << argv[0] << " M N K nbloops nbtries\n";
    std::cout << "Usage: " << argv[0] << " check\n";
    std::cout << "Options:\n";
    std::cout << "    size       unique size of the 2 matrices (integer)\n";
    std::cout << "    auto       automatically set the number of repetitions and tries\n";
    std::cout << "    nbloops    number of times the GEMM routines is executed\n";
    std::cout << "    nbtries    number of times the loop is benched (return the best try)\n";
    std::cout << "    M N K      sizes of the matrices: MxN  =  MxK * KxN (integers)\n";
    std::cout << "    check      check eigen product using cblas as a reference\n";
    exit(1);
  }

  double nbmad = double(M) * double(N) * double(K) * double(nbloops);

  if (!(std::string(argv[1])=="auto"))
    std::cout << M << " x " << N << " x " << K << "\n";

  Scalar alpha, beta;
  MyMatrix ma(M,K), mb(K,N), mc(M,N);
  ma = MyMatrix::Random(M,K);
  mb = MyMatrix::Random(K,N);
  mc = MyMatrix::Random(M,N);

  Eigen::BenchTimer timer;

  // we simply compute c += a*b, so:
  alpha = 1;
  beta = 1;

  // bench cblas
  // ROWS_A, COLS_B, COLS_A, 1.0,  A, COLS_A, B, COLS_B, 0.0, C, COLS_B);
  if (!(std::string(argv[1])=="auto"))
  {
    timer.reset();
    for (uint k=0 ; k<nbtries ; ++k)
    {
        timer.start();
        for (uint j=0 ; j<nbloops ; ++j)
              #ifdef EIGEN_DEFAULT_TO_ROW_MAJOR
              CBLAS_GEMM(CblasRowMajor, CblasNoTrans, CblasNoTrans, M, N, K, alpha, ma.data(), K, mb.data(), N, beta, mc.data(), N);
              #else
              CBLAS_GEMM(CblasColMajor, CblasNoTrans, CblasNoTrans, M, N, K, alpha, ma.data(), M, mb.data(), K, beta, mc.data(), M);
              #endif
        timer.stop();
    }
    if (!(std::string(argv[1])=="auto"))
      std::cout << "cblas: " << timer.value() << " (" << 1e-3*floor(1e-6*nbmad/timer.value()) << " GFlops/s)\n";
    else
        std::cout << M << " : " << timer.value() << " ; " << 1e-3*floor(1e-6*nbmad/timer.value()) << "\n";
  }

  // clear
  ma = MyMatrix::Random(M,K);
  mb = MyMatrix::Random(K,N);
  mc = MyMatrix::Random(M,N);

  // eigen
//   if (!(std::string(argv[1])=="auto"))
  {
      timer.reset();
      for (uint k=0 ; k<nbtries ; ++k)
      {
          timer.start();
          bench_eigengemm(mc, ma, mb, nbloops);
          timer.stop();
      }
      if (!(std::string(argv[1])=="auto"))
        std::cout << "eigen : " << timer.value() << " (" << 1e-3*floor(1e-6*nbmad/timer.value()) << " GFlops/s)\n";
      else
        std::cout << M << " : " << timer.value() << " ; " << 1e-3*floor(1e-6*nbmad/timer.value()) << "\n";
  }

  std::cout << "l1: " << Eigen::l1CacheSize() << std::endl;
  std::cout << "l2: " << Eigen::l2CacheSize() << std::endl;

  return 0;
}

using namespace Eigen;

void bench_eigengemm(MyMatrix& mc, const MyMatrix& ma, const MyMatrix& mb, int nbloops)
{
  for (uint j=0 ; j<nbloops ; ++j)
      mc.noalias() += ma * mb;
}

#define MYVERIFY(A,M) if (!(A)) { \
    std::cout << "FAIL: " << M << "\n"; \
  }
void check_product(int M, int N, int K)
{
  MyMatrix ma(M,K), mb(K,N), mc(M,N), maT(K,M), mbT(N,K), meigen(M,N), mref(M,N);
  ma = MyMatrix::Random(M,K);
  mb = MyMatrix::Random(K,N);
  maT = ma.transpose();
  mbT = mb.transpose();
  mc = MyMatrix::Random(M,N);

  MyMatrix::Scalar eps = 1e-4;

  meigen = mref = mc;
  CBLAS_GEMM(CblasColMajor, CblasNoTrans, CblasNoTrans, M, N, K, 1, ma.data(), M, mb.data(), K, 1, mref.data(), M);
  meigen += ma * mb;
  MYVERIFY(meigen.isApprox(mref, eps),". * .");

  meigen = mref = mc;
  CBLAS_GEMM(CblasColMajor, CblasTrans, CblasNoTrans, M, N, K, 1, maT.data(), K, mb.data(), K, 1, mref.data(), M);
  meigen += maT.transpose() * mb;
  MYVERIFY(meigen.isApprox(mref, eps),"T * .");

  meigen = mref = mc;
  CBLAS_GEMM(CblasColMajor, CblasTrans, CblasTrans, M, N, K, 1, maT.data(), K, mbT.data(), N, 1, mref.data(), M);
  meigen += (maT.transpose()) * (mbT.transpose());
  MYVERIFY(meigen.isApprox(mref, eps),"T * T");

  meigen = mref = mc;
  CBLAS_GEMM(CblasColMajor, CblasNoTrans, CblasTrans, M, N, K, 1, ma.data(), M, mbT.data(), N, 1, mref.data(), M);
  meigen += ma * mbT.transpose();
  MYVERIFY(meigen.isApprox(mref, eps),". * T");
}

void check_product(void)
{
  int M, N, K;
  for (uint i=0; i<1000; ++i)
  {
    M = internal::random<int>(1,64);
    N = internal::random<int>(1,768);
    K = internal::random<int>(1,768);
    M = (0 + M) * 1;
    std::cout << M << " x " << N << " x " << K << "\n";
    check_product(M, N, K);
  }
}

max_stars_repo_name: eundersander/bps-nav

max_issues_repo_name: eundersander/bps-nav

max_forks_repo_name: eundersander/bps-nav


text: #include <boost/mpl/aux_/test/test_case.hpp>

ext: hpp

max_stars_repo_name: miathedev/BoostForArduino

max_issues_repo_name: miathedev/BoostForArduino

max_forks_repo_name: miathedev/BoostForArduino


text:
// Software License for MTL
//
// Copyright (c) 2007 The Trustees of Indiana University.
//               2008 Dresden University of Technology and the Trustees of Indiana University.
//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
// All rights reserved.
// Authors: Peter Gottschling and Andrew Lumsdaine
//
// This file is part of the Matrix Template Library
//
// See also license.mtl.txt in the distribution.

#include <iostream>
#include <boost/numeric/mtl/mtl.hpp>
#include <boost/numeric/itl/itl.hpp>

template <typename PC>
void test(const PC& pc, bool expected)
{
    if (is_identity(pc) != expected) {
	std::cerr << "is_identity should be " << (expected ? "true" : "false") << '\n';
	throw "wrong result";
    }
}

int main()
{
    // For a more realistic example set sz to 1000 or larger
    const int size = 2, N = size * size;

    typedef mtl::compressed2D<double>  matrix_type;
    mtl::compressed2D<double>          A(N, N), dia(N, N);
    laplacian_setup(A, size, size);

    itl::pc::identity<matrix_type>     identity(A);
    itl::pc::ic_0<matrix_type>         ic_0(A);
    itl::pc::ilu_0<matrix_type>        ilu_0(A);
    itl::pc::diagonal<matrix_type>     diagonal(A);

    test(identity, true);
    test(ic_0, false);
    test(ilu_0, false);
    test(diagonal, false);

    return 0;
}

ext: cpp


text:
// COLMAP - Structure-from-Motion and Multi-View Stereo.
// Copyright (C) 2017  Johannes L. Schoenberger <jsch at inf.ethz.ch>
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program 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 General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see <http://www.gnu.org/licenses/>.

#include "controllers/incremental_mapper.h"
#include "controllers/bundle_adjustment.h"
#include "util/logging.h"
#include "util/misc.h"
#include "util/option_manager.h"
#include <Eigen/Geometry>

using namespace colmap;

namespace {
    typedef Eigen::Matrix3d R_t;
    typedef Eigen::Vector3d tvec_t;
    typedef Eigen::Vector4d qvec_t;
    typedef Eigen::Matrix<double, 9, 9> cov_t;
}

void ReadCameraMeasurements(const std::string& path,
                         std::vector<std::string>* image_names,
                         std::vector<tvec_t>* ricI_vec,
                         std::vector<tvec_t>* rciC_vec,
                         std::vector<R_t>* RIC_vec,
                         std::vector<cov_t>* cov_vec) {
  std::vector<std::string> lines = ReadTextFileLines(path);
  for (const auto line : lines) {
    std::stringstream line_parser(line);

    // Buffers to store data in
    std::string image_name = "";
    tvec_t ricI;
    tvec_t rciC;
    R_t RIC;
    cov_t cov;

    // Read data into buffers
    // Translation
    line_parser >>
        image_name >>
        ricI(0) >>
        ricI(1) >>
        ricI(2) >>
        rciC(0) >>
        rciC(1) >>
        rciC(2);

    // Rotation
    for(int i = 0; i < 3; i ++) {
        for(int j = 0; j < 3; j++) {
            line_parser >> RIC(j, i);
        }
    }

    // Covariance
    for(int i = 0; i < 9; i++) {
        for(int j = 0; j < 9; j++) {
            line_parser >> cov(j, i);
        }
    }

    // Push data into vectors
    image_names->push_back(image_name);
    ricI_vec->push_back(ricI);
    rciC_vec->push_back(rciC);
    RIC_vec->push_back(RIC);
    cov_vec->push_back(cov);
  }
}

int main(int argc, char** argv) {
  InitializeGlog(argv);

  std::string import_path;
  std::string priors_path;
  std::string aligned_path;
  std::string metadata_path;

  OptionManager options;
  options.AddRequiredOption("import_path", &import_path);
  options.AddRequiredOption("priors_path", &priors_path);
  options.AddRequiredOption("aligned_path", &aligned_path);
  options.AddRequiredOption("metadata_path", &metadata_path);
  options.AddBundleAdjustmentOptions();
  options.Parse(argc, argv);

  Reconstruction reconstruction;
  reconstruction.Read(import_path);

  /* GAME PLAN
  * 1) Read in the measurements
  * 2) Apply similarity transform to data. Bring reconstruction into global frame.
  * 3) Insert the measurements into Colmap
  * 4) Run global BA
  */

  /* 1) Read in the camera measurements */
  std::vector<std::string> image_names;
  std::vector<tvec_t> ricI_vec;
  std::vector<tvec_t> rciC_vec;
  std::vector<R_t> RIC_vec;
  std::vector<cov_t> cov_vec;
  ReadCameraMeasurements(
          metadata_path,
          &image_names,
          &ricI_vec,
          &rciC_vec,
          &RIC_vec,
          &cov_vec);

  /* 2) Align reconstruction with ECEF measurements to bring it into global frame. */
  RANSACOptions ransac_options;
  ransac_options.max_error = 0.01;
  bool alignment_success = reconstruction.AlignRobust(
          image_names,
          ricI_vec,
          3,
          ransac_options);

  if (alignment_success) {
    std::cout << " => Alignment succeeded" << std::endl;
    reconstruction.Write(aligned_path);
    reconstruction.WriteText(aligned_path);
  } else {
    std::cout << " => Alignment failed" << std::endl;
    return EXIT_FAILURE;
  }

  /* 3) Insert the measurements into reconstruction */
  std::unordered_map< std::string, std::tuple<tvec_t, qvec_t, Eigen::Matrix<double, 6, 6> > > image_priors;
  for(size_t i = 0; i < rciC_vec.size(); i++) {
    image_priors.insert({
                image_names[i],
                std::make_tuple(ricI_vec[i],
                                RotationMatrixToQuaternion(RIC_vec[i]),
                                cov_vec[i].bottomRightCorner<6,6>())      // Bottom right corner is [eIC, ricI]
    });
  }

  reconstruction.AddPriors(image_priors);

  /* 4) Run global BA */
  // Configure BA
  options.bundle_adjustment->cov.compute = false;
  options.bundle_adjustment->priors = true;
  // options.bundle_adjustment->refine_focal_length = true;
  // options.bundle_adjustment->refine_extra_params = true;
  // options.bundle_adjustment->refine_principal_point = true;
  // options.bundle_adjustment->loss_function_type =
  //     BundleAdjustmentOptions::LossFunctionType::CAUCHY;
  // options.bundle_adjustment->loss_function_scale = 5;
  options.bundle_adjustment->solver_options.max_num_iterations = 100;

  // Run BA
  BundleAdjustmentController ba_controller(options, &reconstruction);
  ba_controller.Start();
  ba_controller.Wait();

  // Save output
  reconstruction.Write(priors_path);
  reconstruction.WriteText(priors_path);

  std::cout << "Success!" << std::endl;

  return EXIT_SUCCESS;
}

ext: cc

max_stars_repo_name: radionavlab/colmap

max_issues_repo_name: radionavlab/colmap

max_forks_repo_name: radionavlab/colmap