text:
/*******************************************************************
FUNCTION NAME:   interpolate - performs the interpolation to calculate
                 the pixel value e.g. during resampling

INPUT: interpolation  - interpolation method (nearest neighbor, bilinear, sinc)
       inbuf          - input image buffer
       nLines         - number of lines (image buffer)
       nSamples       - number of samples per line (image buffer)
       xLine          - line of interest
       xSample        - sample position within line
       weighting      - weighting function to be applied (Kaiser, Hamming)
       nKernel        - kernel size for sinc function
*******************************************************************/
#include <assert.h>
#include <math.h>

#include <glib.h>
#include <gsl/gsl_sf_bessel.h>
#include "asf.h"
#include "asf_raster.h"
#include "float_image.h"

#define SQR(X) ((X)*(X))

// Some parameters and weights for bicubic interpolation
#define CUBE(X) ((X)*(X)*(X))
#define P(X) ( X > 0  ? X : 0 )

// For sinc interpolation, the goal is to weight source values using a
// sinc function.  We want to position the peak of the sinc at the
// point of interest.  The sinc filter itself will interpolate a small
// number of point, perhaps 8 or 16, for example.  In order to achieve
// the desired location of the peak, we want to introduce a shift in
// the sinc function computed.  However, we don't want to recompute
// the sinc function for every point to be interpolated!  The solution
// is to precompute an array of slightly shifted sinc functions, and
// then choose the appropriate one for the point we are interpolating.
// For example, if we have known values for 1 and 2, and we want to
// interpolate a value for 1.2, we choose the (1.2 - 1.0) * (2 - 1) *
// NUM_SINCS sinc function to evaluate the result.
#define NUM_SINCS 512 // Number of sinc functions to compute.

// For convenience using buffered image functions
#define GET_PIXEL(x, y) float_image_get_pixel (inbuf, x, y)
#define SET_PIXEL(x, y, value) float_image_set_pixel (inbuf, x, y, value)

// Fetch pixel ii, jj from inbuf, where inbuf is taken to be an image
// nSamples wide and nLines high.  If ii, jj is outside the image, a
// "reflected" value is returned.  For example, these two calls will
// return the same pixel:
//
//     get_pixel_with_reflection (inbuf, sample_count, line_count, -1, -2)
//     get_pixel_with_reflection (inbuf, sample_count, line_count, 1, 2);
//
// If the reflected indicies are still outside the image, an exception
// is triggered.

//static float
//get_pixel_with_reflection (float *inbuf, int nSamples, int nLines, int sample,
//			   int line)
//{
//  int ii = sample;		// Convenience alias.
//  int jj = line;		// Convenience alias.
//  // Handle reflection at image edges.
//  if ( G_UNLIKELY (ii < 0) ) {
//    ii = -ii;
//  }
//  else if ( G_UNLIKELY (ii >= nSamples) ) {
//    ii = ii - (ii - nSamples + 1) - 1;
//  }
//  if ( G_UNLIKELY (jj < 0) ) {
//    jj = -jj;
//  }
//  else if ( G_UNLIKELY (jj >= nLines) ) {
//    jj = jj - (jj - nLines + 1) - 1;
//  }
//  // Now we better be in the image.
//  assert (ii >= 0 && ii < nSamples);
//  assert (jj >= 0 && jj < nLines);
//
//  return inbuf[nSamples * jj + ii];
//}

void samples2coefficients(FloatImage *inbuf, char dimension)
{
  int ii, kk;
  float z = sqrt(3.0) - 2.0;
  float lambda = (1.0 - z) * (1.0 - 1.0/z);
  float value;

  // Apply overall gain
  for (ii=0; ii<inbuf->size_x; ii++)
    for (kk=0; kk<inbuf->size_y; kk++) {
      value = GET_PIXEL(ii,kk) * lambda;
      SET_PIXEL(ii,kk,value);
    }

  for (kk=0; kk<inbuf->size_y; kk++) {
    // Causal initialization (first pixel of line)

  }

  // Causal recursion

  // Anticausal initialization

  // Anticausal recursion

}

float interpolate(interpolate_type_t interpolation, FloatImage *inbuf, float yLine,
		  float xSample, weighting_type_t weighting, int sinc_points)
{
  int ix, iy, ii, kk, ll, xI[4], yI[4];
  float a00, a10, a01, a11;
  float value, dx, dy, wx[4], wy[4], w;

  /* Get closed pixel position to start from */

  switch ( interpolation )
    {
    case NEAREST:
      ix = (int) (xSample + 0.5);
      iy = (int) (yLine + 0.5);
      value = GET_PIXEL(ix,iy);
      break;

    case BILINEAR:
      assert (xSample >= 0.0);
      assert (yLine >= 0.0);
      assert (xSample <= inbuf->size_x - 1);
      assert (yLine <= inbuf->size_y - 1);
      ix = floor(xSample);
      iy = floor(yLine);

      a00 = GET_PIXEL(ix,iy);
      a10 = GET_PIXEL(ix+1,iy) - GET_PIXEL(ix,iy);
      a01 = GET_PIXEL(ix,iy+1) - GET_PIXEL(ix,iy);
      a11 = GET_PIXEL(ix,iy) - GET_PIXEL(ix+1,iy) - GET_PIXEL(ix,iy+1)
	+ GET_PIXEL(ix+1,iy+1);
      value = (a00 + a10 * (xSample - ix) + a01 * (yLine - iy)
	       + a11 * (xSample - ix) * (yLine - iy));
      break;

    case BICUBIC:
      assert (xSample >= 0.0);
      assert (yLine >= 0.0);
      assert (xSample <= inbuf->size_x - 1);
      assert (yLine <= inbuf->size_y - 1);
      ix = floor(xSample);
      iy = floor(yLine);
      dx = xSample - ix;
      dy = yLine - iy;

      // Calculating weights
      for (ii=-1; ii<=2; ii++) {
	wx[ii+1] = 1.0/6.0 * (CUBE(P(ii-dx+2)) - 4*CUBE(P(ii-dx+1))
			      + 6*CUBE(P(ii-dx)) - 4*CUBE(P(ii-dx-1)));
	wy[ii+1] = 1.0/6.0 * (CUBE(P(dy-ii+2)) - 4*CUBE(P(dy-ii+1))
			      + 6*CUBE(P(dy-ii)) - 4*CUBE(P(dy-ii-1)));
      }

      // Get the interpolated pixel value
      value = 0.0;
      for (ii=-1; ii<=2; ii++)
	for (kk=-1; kk<=2; kk++)
	  value += GET_PIXEL(ix+ii,iy+kk) * wx[ii+1] * wy[kk+1];

      break;

    case SPLINES:
      // Determine B-Spline coefficients.
      // This is a essentially a pre-filtering applied to the entire image
      // that converts the image from a sample representation into a
      // representation based on B-spline coefficients.
      // Without this step the result would look blurry at the end.

      // Convert samples to coefficients in x direction
      samples2coefficients(inbuf, 'x');

      // Convert samples to coefficients in y direction
      samples2coefficients(inbuf, 'y');

      // Calculate interpolation indices
      ii = (int)floor(xSample) - 1;
      kk = (int)floor(yLine) - 1;
      for (ll=0; ll<=3; ll++) {
	xI[ll] = ii++;
	yI[ll] = kk++;
      }

      // Calculate interpolation weights
      w = xSample - (float)xI[1];
      wx[3] = (1.0/6.0)*CUBE(xSample);
      wx[0] = (1.0/6.0) + (1.0/2.0)*w*(w-1.0) - wx[3];
      wx[2] = w + wx[0] - 2.0*wx[3];
      wx[1] = 1.0 - wx[0] - wx[2] - wx[3];
      w = yLine - (float)yI[1];
      wy[3] = (1.0/6.0)*CUBE(xSample);
      wy[0] = (1.0/6.0) + (1.0/2.0)*w*(w-1.0) - wy[3];
      wy[2] = w + wy[0] - 2.0*wy[3];
      wy[1] = 1.0 - wy[0] - wy[2] - wy[3];

      // Get interpolation value
      value = 0.0;
      for (kk=0; kk<=3; kk++) {
	w = 0.0;
	for (ii=0; ii<=3; ii++)
	  w += wx[ii] * 1;
	value += wy[kk] * w;
      }
      break;
    case SINC:
      break;
    default:
      assert (FALSE);
      break;
  }

  return value;
}

max_stars_repo_name: glshort/MapReady

max_issues_repo_name: glshort/MapReady

max_forks_repo_name: glshort/MapReady


text:
#pragma once

#include "zerrors.h"

#include <string>
#include <vector>
#include <numeric>

#include <cwchar>
#include <mk_wcswidth/mk_wcwidth.h>

#include <gsl/span>

namespace terminal_editor {

#if !defined(USE_NATIVE_WCWIDTH)

/// Returns number of columns required to represent given character.
/// Returns 0 for zero width caracters.
/// Throws when non-printable character is passed as argument.
/// @param ucs      UTF-32 code point.
inline int wcwidth(uint32_t ucs) {
    auto width = mk::wcwidth(ucs);
    ZASSERT(width >= 0) << "Non-printable character encountered: " << ucs;
    return width;
}

#else

inline int wcwidth(uint32_t ucs) {
    auto width = ::wcwidth(ucs);
    ZASSERT(width >= 0) << "Non-printable character encountered: " << ucs;
    return width;
}

#endif

/// Returns number of columns required to represent given string.
/// Returns 0 for zero width strings.
/// Throws when non-printable characters are encountered.
/// @param ucs      UTF-32 code point.
inline int wcswidth(gsl::span<const uint32_t> text) {
    return std::accumulate(text.begin(), text.end(), 0, [](int width, const uint32_t ucs) { return width + wcwidth(ucs); });
}

} // namespace terminal_editor

ext: h

max_stars_repo_name: Zbyl/terminal-editor

max_issues_repo_name: Zbyl/terminal-editor

max_forks_repo_name: Zbyl/terminal-editor


text:
// Copyright (c) Microsoft Corporation.
// Licensed under the MIT License.

#pragma once

#include "BaseWorker.h"

#include "MageSlam.h"

#include <arcana/scheduling/state_machine.h>
#include <gsl/gsl>
#include <memory>

namespace mage
{
    struct MageContext;
    struct MageSlamSettings;
    struct FrameData;

    class Fuser;

    class Runtime
    {
    public:
        Runtime(const MageSlamSettings& settings, MageContext& context, Fuser& fuser, mira::state_machine_driver& driver);
        ~Runtime();

        void Run(gsl::span<const MAGESlam::CameraConfiguration> cameras);

        void TrackMono(std::shared_ptr<FrameData> frame);
        void TrackStereo(std::shared_ptr<FrameData> one, std::shared_ptr<FrameData> two);

        void AddSample(const mage::SensorSample& sample);
    private:
        struct Impl;
        const std::unique_ptr<Impl> m_impl;
    };
}

max_stars_repo_name: syntheticmagus/mageslam

max_issues_repo_name: syntheticmagus/mageslam

max_forks_repo_name: syntheticmagus/mageslam


text:
#pragma once

#include "Cesium3DTilesSelection/Library.h"
#include <glm/mat4x4.hpp>
#include <glm/vec2.hpp>
#include <gsl/span>

namespace CesiumGeometry {
struct Rectangle;
}

namespace CesiumGltf {
struct ImageCesium;
struct Model;
} // namespace CesiumGltf

namespace Cesium3DTilesSelection {

class Tile;
class RasterOverlayTile;

/**
 * @brief When implemented for a rendering engine, allows renderer resources to
 * be created and destroyed under the control of a {@link Tileset}.
 *
 * It is not supposed to be used directly by clients. It is implemented
 * for specific rendering engines to provide an infrastructure for preparing the
 * data of a {@link Tile} so that it can be used for rendering.
 *
 * Instances of this class are associated with a {@link Tileset}, in the
 * {@link TilesetExternals} structure that can be obtained
 * via {@link Tileset::getExternals}.
 */
class CESIUM3DTILESSELECTION_API IPrepareRendererResources {
public:
  virtual ~IPrepareRendererResources() = default;

  /**
   * @brief Prepares renderer resources for the given tile.
   *
   * This method is invoked in the load thread, and it may not modify the tile.
   *
   * @param model The glTF model to prepare.
   * @param transform The tile's transformation.
   * @returns Arbitrary data representing the result of the load process. This
   * data is passed to {@link prepareInMainThread} as the `pLoadThreadResult`
   * parameter.
   */
  virtual void* prepareInLoadThread(
      const CesiumGltf::Model& model,
      const glm::dmat4& transform) = 0;

  /**
   * @brief Further prepares renderer resources.
   *
   * This is called after {@link prepareInLoadThread}, and unlike that method,
   * this one is called from the same thread that called
   * {@link Tileset::updateView}.
   *
   * @param tile The tile to prepare.
   * @param pLoadThreadResult The value returned from
   * {@link prepareInLoadThread}.
   * @returns Arbitrary data representing the result of the load process.
   * Note that the value returned by {@link prepareInLoadThread} will _not_ be
   * automatically preserved and passed to {@link free}. If you need to free
   * that value, do it in this method before returning. If you need that value
   * later, add it to the object returned from this method.
   */
  virtual void* prepareInMainThread(Tile& tile, void* pLoadThreadResult) = 0;

  /**
   * @brief Frees previously-prepared renderer resources.
   *
   * This method is always called from the thread that called
   * {@link Tileset::updateView} or deleted the tileset.
   *
   * @param tile The tile for which to free renderer resources.
   * @param pLoadThreadResult The result returned by
   * {@link prepareInLoadThread}. If {@link prepareInMainThread} has
   * already been called, this parameter will be `nullptr`.
   * @param pMainThreadResult The result returned by
   * {@link prepareInMainThread}. If {@link prepareInMainThread} has
   * not yet been called, this parameter will be `nullptr`.
   */
  virtual void free(
      Tile& tile,
      void* pLoadThreadResult,
      void* pMainThreadResult) noexcept = 0;

  /**
   * @brief Prepares a raster overlay tile.
   *
   * This method is invoked in the load thread, and it may not modify the tile.
   *
   * @param image The raster tile image to prepare.
   * @returns Arbitrary data representing the result of the load process. This
   * data is passed to {@link prepareRasterInMainThread} as the
   * `pLoadThreadResult` parameter.
   */
  virtual void*
  prepareRasterInLoadThread(const CesiumGltf::ImageCesium& image) = 0;

  /**
   * @brief Further preprares a raster overlay tile.
   *
   * This is called after {@link prepareRasterInLoadThread}, and unlike that
   * method, this one is called from the same thread that called
   * {@link Tileset::updateView}.
   *
   * @param rasterTile The raster tile to prepare.
   * @param pLoadThreadResult The value returned from
   * {@link prepareInLoadThread}.
   * @returns Arbitrary data representing the result of the load process. Note
   * that the value returned by {@link prepareRasterInLoadThread} will _not_ be
   * automatically preserved and passed to {@link free}. If you need to free
   * that value, do it in this method before returning. If you need that value
   * later, add it to the object returned from this method.
   */
  virtual void* prepareRasterInMainThread(
      const RasterOverlayTile& rasterTile,
      void* pLoadThreadResult) = 0;

  /**
   * @brief Frees previously-prepared renderer resources for a raster tile.
   *
   * This method is always called from the thread that called
   * {@link Tileset::updateView} or deleted the tileset.
   *
   * @param rasterTile The tile for which to free renderer resources.
   * @param pLoadThreadResult The result returned by
   * {@link prepareRasterInLoadThread}. If {@link prepareRasterInMainThread}
   * has already been called, this parameter will be `nullptr`.
   * @param pMainThreadResult The result returned by
   * {@link prepareRasterInMainThread}. If {@link prepareRasterInMainThread}
   * has not yet been called, this parameter will be `nullptr`.
   */
  virtual void freeRaster(
      const RasterOverlayTile& rasterTile,
      void* pLoadThreadResult,
      void* pMainThreadResult) noexcept = 0;

  /**
   * @brief Attaches a raster overlay tile to a geometry tile.
   *
   * @param tile The geometry tile.
   * @param overlayTextureCoordinateID The ID of the overlay texture coordinate
   * set to use.
   * @param rasterTile The raster overlay tile to add. The raster tile will have
   * been previously prepared with a call to {@link prepareRasterInLoadThread}
   * followed by {@link prepareRasterInMainThread}.
   * @param pMainThreadRendererResources The renderer resources for this raster
   * tile, as created and returned by {@link prepareRasterInMainThread}.
   * @param textureCoordinateRectangle Defines the range of texture coordinates
   * in which this raster tile should be applied, in the order west, south, east
   * north. Each coordinate is in the range 0.0 (southwest corner) to 1.0
   * (northeast corner).
   * @param translation The translation to apply to the texture coordinates
   * identified by `overlayTextureCoordinateID`. The texture coordinates to use
   * to sample the raster image are computed as `overlayTextureCoordinates *
   * scale + translation`.
   * @param scale The scale to apply to the texture coordinates identified by
   * `overlayTextureCoordinateID`. The texture coordinates to use to sample the
   * raster image are computed as `overlayTextureCoordinates * scale +
   * translation`.
   */
  virtual void attachRasterInMainThread(
      const Tile& tile,
      uint32_t overlayTextureCoordinateID,
      const RasterOverlayTile& rasterTile,
      void* pMainThreadRendererResources,
      const CesiumGeometry::Rectangle& textureCoordinateRectangle,
      const glm::dvec2& translation,
      const glm::dvec2& scale) = 0;

  /**
   * @brief Detaches a raster overlay tile from a geometry tile.
   *
   * @param tile The geometry tile.
   * @param overlayTextureCoordinateID The ID of the overlay texture coordinate
   * set to which the raster tile was previously attached.
   * @param rasterTile The raster overlay tile to remove.
   * @param pMainThreadRendererResources The renderer resources for this raster
   * tile, as created and returned by {@link prepareRasterInMainThread}.
   * @param textureCoordinateRectangle Defines the range of texture coordinates
   * in which this raster tile should be applied, in the order west, south, east
   * north. Each coordinate is in the range 0.0 (southwest corner) to 1.0
   * (northeast corner).
   */
  virtual void detachRasterInMainThread(
      const Tile& tile,
      uint32_t overlayTextureCoordinateID,
      const RasterOverlayTile& rasterTile,
      void* pMainThreadRendererResources,
      const CesiumGeometry::Rectangle& textureCoordinateRectangle) noexcept = 0;
};

} // namespace Cesium3DTilesSelection

max_stars_repo_name: zrkcode/cesium-native

max_issues_repo_name: zrkcode/cesium-native

max_forks_repo_name: zrkcode/cesium-native


text:
/* -*- c++ -*- */
/*
 * Copyright 2015 Free Software Foundation, Inc.
 *
 * This 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, or (at your
 * option) any later version.
 *
 * This software 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 software; see the file COPYING.  If not, write to
 * the Free Software Foundation, Inc., 51 Franklin Street,
 * Boston, MA 02110-1301, USA.
 */

#ifndef INCLUDED_ldpc_G_matrix_impl_H
#define INCLUDED_ldpc_G_matrix_impl_H

#include <gsl/gsl_randist.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>

#include "fec_mtrx_impl.h"
#include <gnuradio/fec/ldpc_G_matrix.h>
#include <gnuradio/fec/ldpc_H_matrix.h>
#include <gnuradio/logger.h>

namespace gr {
  namespace fec {
    namespace code {
      /*!
       * \brief Class for storing H or G matrix
       * \ingroup error_coding_blk
       *
       * \details
       * This class stores a GSL matrix variable, specifically
       * either a:
       *
       * 1) Generator matrix, G, in the standard format G = [I P],
       *    where I is an identity matrix and P is the parity
       *    submatrix.
       *
       * or
       *
       * 2) Parity matrix, H, in the standard format H = [P' I],
       *    where P' is the transpose of the parity submatrix and I
       *    is an identity matrix.
       *
       * This variable can used by the ldpc_gen_mtrx_encoder and
       * ldpc_bit_flip_decoder classes.
       */
      class ldpc_G_matrix_impl
        : public fec_mtrx_impl, public ldpc_G_matrix
      {
      private:
        // GSL matrix structure for transpose of G
        gsl_matrix *d_G_transp_ptr;

        gsl_matrix *d_H_obj;

        //! Get the generator matrix (used during encoding)
        const gsl_matrix *G_transpose() const;

        gr::logger_ptr d_logger;
        gr::logger_ptr d_debug_logger;

      public:
        ldpc_G_matrix_impl(const std::string filename);

        void encode(unsigned char *outbuffer,
                    const unsigned char *inbuffer) const;

        void decode(unsigned char *outbuffer,
                    const float *inbuffer,
                    unsigned int frame_size,
                    unsigned int max_iterations) const;

        unsigned int n() const { return fec_mtrx_impl::n(); }

        unsigned int k() const { return fec_mtrx_impl::k(); }

        gsl_matrix* generate_H();

        gr::fec::code::fec_mtrx_sptr get_base_sptr();

        /*!
         * \brief Destructor
         * \details
         * Calls the gsl_matrix_free function to free memory.
         */
        virtual ~ldpc_G_matrix_impl();
      };

    }
  }
}

#endif /* INCLUDED_ldpc_G_matrix_impl_H */


text:
#include <stdio.h>
#include <string.h>
#include <time.h>
#include <MaxSLiCInterface.h>
#include "Maxfiles.h"
#include "dgemm.h"

#define TILE_SIZE DGEMM_tileSize
#define TILE_SIZE_2D (TILE_SIZE * TILE_SIZE)
#define FREQUENCY DGEMM_frequency

static size_t sizeInTiles(size_t n) {
	return (n + TILE_SIZE - 1) / TILE_SIZE;
}

#ifdef USE_BLAS
#include <cblas.h>
#endif

void dgemm_model(
		const char* trans_a, const char* trans_b, size_t m, size_t n, size_t k, double alpha,
		const double* A, size_t lda, const double* B, size_t ldb, double beta, double* C, size_t ldc)
{
	// TODO support transform

#ifdef USE_BLAS
	cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc);
#else
	for (size_t mm = 0; mm < m; ++mm) {
		for (size_t nn = 0; nn < n; ++nn) {
			C[mm*ldc+nn] *= beta;
			for (size_t kk = 0; kk < k; ++kk) {
				C[mm*ldc+nn] += alpha * A[mm*lda+kk] * B[kk*ldb+nn];
			}
		}
	}
#endif
}

static struct timespec diff(struct timespec start, struct timespec finish) {
	struct timespec diff;
	if ((finish.tv_nsec < start.tv_nsec)) {
		diff.tv_sec  = finish.tv_sec  - start.tv_sec  - 1;
		diff.tv_nsec = finish.tv_nsec - start.tv_nsec + 1000000000;
	} else {
		diff.tv_sec  = finish.tv_sec  - start.tv_sec;
		diff.tv_nsec = finish.tv_nsec - start.tv_nsec;
	}

	return diff;
}

static int check_result(int m, int n, const double* expected, const double* actual) {
	for (int i = 0; i < m; ++i) {
		for (int j = 0; j < n; ++j) {
			int idx = i*n + j;
			if (expected[idx] != actual[idx]) {
				return 1;
			}
		}
	}
	return 0;
}

int main(int argc, char** argv) {
	unsigned seed = time(NULL);
	printf("Random seed: %u\n", seed);
	srand(seed);

	size_t m;
	size_t n;
	size_t k;

	if (argc == 1) {
		m = (random() % (5 * TILE_SIZE)) + TILE_SIZE;
		n = (random() % (5 * TILE_SIZE)) + TILE_SIZE;
		k = (random() % (5 * TILE_SIZE)) + TILE_SIZE;
	} else if (argc == 2) {
		m = n = k = atoi(argv[1]);
	} else {
		printf("Usage: %s [size]\n", argv[0]);
		exit(1);
	}

	double* A   = calloc(m * k, sizeof(double));
	double* B   = calloc(k * n, sizeof(double));
	double* Csw = calloc(m * n, sizeof(double));
	double* Chw = calloc(m * n, sizeof(double));

	printf("Matrix dimensions: m = %lu, n = %lu, k = %lu\n", m, n, k);

	dgemm_init(m, n, k);

	size_t mTiles = sizeInTiles(m);
	size_t nTiles = sizeInTiles(n);
	size_t kTiles = sizeInTiles(k);

	printf("DFE tile size: %d\n", TILE_SIZE);
	printf("DFE compute dimensions: m = %zu, n = %zu, k = %zu\n",
			mTiles * TILE_SIZE, nTiles * TILE_SIZE, kTiles * TILE_SIZE);

	double cpuPoints = ((double) m) * n * k;
	double dfePoints = ((double) mTiles) * nTiles * kTiles * TILE_SIZE_2D;
	printf("DFE compute efficiency: %f\n", cpuPoints / (dfePoints * TILE_SIZE));
	printf("DFE frequency: %d MHz\n", FREQUENCY);
	// add one tile to account for stream offset TODO measure pipeline depth properly?
	printf("DFE predicted compute time: %f s\n", (dfePoints + TILE_SIZE_2D) / (((double) FREQUENCY) * 1e6));

	for (size_t i = 0; i < m*k; ++i) {
		A[i] = random() % 100;
	}

	for (size_t i = 0; i < k*n; ++i) {
		B[i] = random() % 100;
	}

	for (size_t i = 0; i < m*n; ++i) {
		Csw[i] = random() % 100;
		Chw[i] = Csw[i];
	}

	// TODO random
	double alpha = 1;
	double beta  = 0;

	struct timespec start;
	struct timespec finish;

	for (int xx = 0; xx < 5; ++xx) {
		printf("Repeat %d\n", xx);
		printf("Running HW... "); fflush(stdout);
		clock_gettime(CLOCK_REALTIME, &start);
		dgemm("n", "n", m, n, k, alpha, A, k, B, n, beta, Chw, n);
		clock_gettime(CLOCK_REALTIME, &finish);

		struct timespec hw_time  = diff(start, finish);

		printf("took: %ld.%09ld s\n", hw_time.tv_sec, hw_time.tv_nsec);
	}

	printf("Running SW... "); fflush(stdout);
	clock_gettime(CLOCK_REALTIME, &start);
	dgemm_model("n", "n", m, n, k, alpha, A, k, B, n, beta, Csw, n);
	clock_gettime(CLOCK_REALTIME, &finish);

	struct timespec sw_time = diff(start, finish);

	printf("took: %ld.%09ld s\n", sw_time.tv_sec, sw_time.tv_nsec);

	printf("Comparing results... ");
	int failed = check_result(m, n, Csw, Chw);
	printf("%s\n", failed ? "WRONG" : "CORRECT");
	printf("Done.\n");
	return failed;
}

ext: c

max_stars_repo_name: nilsv/Dense-Matrix-Multiplication


text:
/**
 *
 * @file core_ztrmdm.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.4.5
 * @author Dulceneia Becker
 * @date 2011-1-18
 **/
/*
 * @precisions normal z -> c d s
 */
#include <lapacke.h>
#include "dplasma_cores.h"
#include "dplasma_zcores.h"

#if defined(PARSEC_HAVE_STRING_H)
#include <string.h>
#endif  /* defined(PARSEC_HAVE_STRING_H) */
#if defined(PARSEC_HAVE_STDARG_H)
#include <stdarg.h>
#endif  /* defined(PARSEC_HAVE_STDARG_H) */
#include <stdio.h>
#ifdef PARSEC_HAVE_LIMITS_H
#include <limits.h>
#endif

#include <cblas.h>
#include <core_blas.h>

#define max(a, b) ((a) > (b) ? (a) : (b))
#define min(a, b) ((a) < (b) ? (a) : (b))

int CORE_ztrmdm(int uplo, int N, PLASMA_Complex64_t *A, int LDA);

/***************************************************************************//**
 *
 * @ingroup CORE_PLASMA_Complex64_t
 *
 * CORE_ztrmdm scales the strictly upper or strictly lower triangular part of a
 * square matrix by the inverse of a diagonal matrix, ie performs either
 *
 *    A := L / D  or  A := U / D (only for triangular part above or below diagonal)
 *
 * where:
 *
 *    L is a strictly lower triangular matrix stored as the strictly lower triangle in A
 *    U is a strictly upper triangular matrix stored as the strictly upper triangle in A
 *    D is a diagonal matrix stored as the diagonal in A
 *
 * The diagonal elements of A are not changed.
 *
*******************************************************************************
 *
 * @param[in] uplo
 *         INTEGER
 *         @arg PlasmaLower: Lower triangle of A is stored and scaled.
 *         @arg PlasmaUpper: Upper triangle of A is stored and scaled.
 *
 * @param[in] N
 *         INTEGER
 *         The number of rows and columns of A.  N >= 0.
 *
 * @param[in,out] A
 *         PLASMA_Complex64_t array, dimension (LDA,N)
 *
 *         On entry, the triangular matrix A. If uplo = 'U', the leading
 *         N-by-N upper triangular part of A contains the upper
 *         triangular part of the matrix A, and the strictly lower
 *         triangular part of A is not referenced.  If uplo = 'L', the
 *         leading N-by-N lower triangular part of A contains the lower
 *         triangular part of the matrix A, and the strictly upper
 *         triangular part of A is not referenced.
 *
 *         On exit, the strictly lower or upper triangular part of A
 *         scaled by the diagonal elements of A.
 *
 * @param[in] LDA
 *         INTEGER
 *         The leading dimension of the array A.  LDA >= max(1,N).
 *
*******************************************************************************
 *
 * @return
 *          \retval PLASMA_SUCCESS successful exit
 *          \retval <0 if -i, the i-th argument had an illegal value
 *
 ******************************************************************************/
#if defined(PLASMA_PARSEC_HAVE_WEAK)
#pragma weak CORE_ztrmdm = PCORE_ztrmdm
#define CORE_ztrmdm PCORE_ztrmdm
#endif
int CORE_ztrmdm(int uplo, int N, PLASMA_Complex64_t *A, int LDA)
{
    static PLASMA_Complex64_t zone = 1.0;

    PLASMA_Complex64_t alpha;
    int j;

    /* Check input arguments */
    if (uplo != PlasmaUpper && uplo != PlasmaLower) {
        coreblas_error(1, "Illegal value of UPLO");
        return -1;
    }
    if (N < 0) {
        coreblas_error(2, "Illegal value of N");
        return -2;
    }
    if (LDA < max(1, N)) {
        coreblas_error(1, "Illegal value of LDA");
        return -4;
    }

    /* Quick return */
    if (max(N, 0) == 0)
        return PLASMA_SUCCESS;

    /**/

    if (uplo==PlasmaLower) {
        for (j=0; j<(N-1); j++) {
                alpha = zone / A[LDA*j+j];
                cblas_zscal(N-j-1, CBLAS_SADDR(alpha), &A[LDA*j+j+1], 1);
        }

    } else if (uplo==PlasmaUpper) {
        for (j=1; j<N; j++) {
                alpha = zone / A[LDA*j+j];
                cblas_zscal(j, CBLAS_SADDR(alpha), &A[LDA*j], 1);
        }

    }

    return PLASMA_SUCCESS;

}

max_stars_repo_name: NLAFET/ABFT

max_issues_repo_name: NLAFET/ABFT

max_forks_repo_name: NLAFET/ABFT


text:
/**
 *
 * @file core_zlaset2.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Hatem Ltaief
 * @date 2010-11-15
 * @precisions normal z -> c d s
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 * @ingroup CORE_PLASMA_Complex64_t
 *
 *  CORE_zlaset2 - Sets the elements of the matrix A to alpha.
 *  Not LAPACK compliant! Read below.
 *
 *******************************************************************************
 *
 * @param[in] uplo
 *          Specifies which elements of the matrix are to be set
 *          = PlasmaUpper: STRICT Upper part of A is set to alpha;
 *          = PlasmaLower: STRICT Lower part of A is set to alpha;
 *          = PlasmaUpperLower: ALL elements of A are set to alpha.
 *          Not LAPACK Compliant.
 *
 * @param[in] M
 *          The number of rows of the matrix A.  M >= 0.
 *
 * @param[in] N
 *         The number of columns of the matrix A.  N >= 0.
 *
 * @param[in] alpha
 *         The constant to which the elements are to be set.
 *
 * @param[in,out] A
 *         On entry, the M-by-N tile A.
 *         On exit, A has been set to alpha accordingly.
 *
 * @param[in] LDA
 *         The leading dimension of the array A.  LDA >= max(1,M).
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_zlaset2 = PCORE_zlaset2
#define CORE_zlaset2 PCORE_zlaset2
#endif
void CORE_zlaset2(PLASMA_enum uplo, int M, int N,
                  PLASMA_Complex64_t alpha, PLASMA_Complex64_t *A, int LDA)
{
    if (uplo == PlasmaUpper) {
        LAPACKE_zlaset_work(
            LAPACK_COL_MAJOR,
            lapack_const(uplo),
            M, N-1, alpha, alpha, A+LDA, LDA);
    }
    else if (uplo == PlasmaLower) {
        LAPACKE_zlaset_work(
            LAPACK_COL_MAJOR,
            lapack_const(uplo),
            M-1, N, alpha, alpha, A+1, LDA);
    }
    else {
        LAPACKE_zlaset_work(
            LAPACK_COL_MAJOR,
            lapack_const(uplo),
            M, N, alpha, alpha, A, LDA);
    }
}