Commit f5449681 authored by Tandri Gauksson's avatar Tandri Gauksson

FastMatMult synced with CodeExpert

parent f5afa20c
# Code Expert Generic Project Configuration.
#
# In case this file is omitted, the execution is done with the default values.
# You can either define a cmd or a script file to be executed.
# Default to script file setting a cmd overrides it.
cmd:
# compile: echo "COMPILE"
# run: echo "RUN"
# test: echo "TEST"
# submit: echo "SUBMIT"
# All scripts are set by default to "/scripts/${ACTION}.sh"
script:
compile: /scripts/compile.sh
run: /scripts/run.sh
test: /scripts/test.sh
submit: /scripts/submit.sh
\ No newline at end of file
descriptions:
default:
name: Default
default: true
# Optionally add mulitple versions here.
# Use the file name without the .md extension as key and add a `name:` to it.
\ No newline at end of file
# Problem 2-4: Fast matrix multiplication with EIGEN
> For the task description of this exercise, please refer to [NCSE19_Problems.pdf](
https://www.sam.math.ethz.ch/~grsam/NCSE19/HOMEWORK/NCSE19_Problems.pdf).
> Open "strassen.hpp" and fill in the missing code in between the delimiters `// START` and `// END` according to the instructions preceded by `// TO DO:`.
# Tests of gramschmidt.hpp
> There are three functions in strassen.hpp and main.cpp has tests for the first two. The third function is a timing exercise, which should not be tested automatically.
***
> `strassenMatMult(const MatrixXd& A, const MatrixXd& B)`: We test this function by calculating norm(Strassen(A,Id) - A) and also norm(Strassen(Id,A) - A) for a random 128*128 matrix A. The tolerance is 10^{-9}.
***
> `test_strassen()`: This function is a test written by the student. It should create two random matrices A and B, and return err = norm(Strassen(A,B) - A*B).
The test is passed if err < 10^{-9}.
\ No newline at end of file
#include <Eigen/Dense>
#include <iostream>
#include "strassen.hpp"
using namespace Eigen;
int main() {
double tolerance = 1e-9;
int n = 128;
std::srand(5);
MatrixXd A = MatrixXd::Random(n,n);
std::cout << "Works for identity matrix: "
<< ((strassenMatMult(A, MatrixXd::Identity(n,n)) - A).norm()<tolerance)
<< " "
<< ((strassenMatMult(MatrixXd::Identity(n,n),A) - A).norm()<tolerance)
<< std::endl;
std::cout << "Works for two random matrices: "
<<(test_strassen() < tolerance) << std::endl;
std::cout << "Output of test_strassen()" << test_strassen() << std::endl;
time_strassen();
}
#!/bin/bash
echo "Compiling ..."
# compile c code
#mkdir -p bin
#g++ test_strassen.cpp -I/usr/include/eigen3/ -fdiagnostics-color=always --pedantic -Wextra -Wall -std=c++14 \
# -o bin/test_strassen.out
#g++ timing.cpp -I/usr/include/eigen3/ -fdiagnostics-color=always --pedantic -Wextra -Wall -std=c++14 \
# -o bin/timing.out
mkdir -p bin
find . -iname '*.cpp' | sort | xargs \
g++ -fdiagnostics-color=always -std=c++11 \
-I/usr/local/include/python3.7m \
-I/usr/local/lib/python3.7/site-packages/numpy/core/include \
-I/usr/include/eigen3/ \
-lpython3.7m \
-lpthread -lutil -ldl \
-Xlinker -export-dynamic \
-o bin/a.out
echo "Compilation successful"
\ No newline at end of file
#!/bin/bash
export PYTHONIOENCODING="UTF-8"
# compile (call compile script)
bash "${WORKDIR}/scripts/compile.sh"
echo "Running script ... "
# run
bin/a.out
\ No newline at end of file
#!/bin/bash
# I would like to delete this file, but I cannot.
export PYTHONIOENCODING="UTF-8"
# compile (call compile script)
bash "${WORKDIR}/scripts/compile.sh"
echo "Running code ..."
echo " "
echo " "
# run
echo "Task 1-4.b) 'Testing the correctness of the Strassen Alorithm' "
echo " "
bin/test_strassen.out
echo " "
echo "Task 1-4.c) 'Comparing the timings of A*B and the Strassen Algorithm' "
echo " "
bin/timing.out
\ No newline at end of file
#!/bin/bash
echo "Not implemented"
\ No newline at end of file
#!/bin/bash
# run (call run script)
bash "${WORKDIR}/scripts/run.sh"
\ No newline at end of file
////
//// Copyright (C) 2016 SAM (D-MATH) @ ETH Zurich
//// Author(s): lfilippo <filippo.leonardi@sam.math.ethz.ch>
//// Contributors: tille, jgacon, dcasati
//// This file is part of the NumCSE repository.
////
#include <Eigen/Dense>
#include <iomanip>
#include <iostream>
#include <vector>
#include "timer.h"
using namespace Eigen;
/* \brief Compute the Matrix product $A \times B$ using Strassen's algorithm.
* \param[in] A Matrix $2^k \times 2^k$
* \param[in] B Matrix $2^k \times 2^k$
* \param[out] Matrix product of A and B of dim $2^k \times 2^k$
*/
/* SAM_LISTING_BEGIN_0 */
MatrixXd strassenMatMult(const MatrixXd& A, const MatrixXd& B) {
// Ensure square matrix
assert(A.rows() == A.cols() && "Matrix A must be square");
assert(B.rows() == B.cols() && "Matrix B must be square");
// Matrix dimension must be a power of 2
assert(A.rows() % 2 == 0 && "Matrix dimensions must be a power of two.");
const unsigned n = A.rows();
// The function is recursive and acts on blocks of size $n/2 \times n/2$
// i.e. exploits fast product of 2x2 block matrix
if ( n==2 ) { // End of recursion
MatrixXd C(2, 2);
C << A(0,0)*B(0,0) + A(0,1)*B(1,0),
A(0,0)*B(0,1) + A(0,1)*B(1,1),
A(1,0)*B(0,0) + A(1,1)*B(1,0),
A(1,0)*B(0,1) + A(1,1)*B(1,1);
return C;
}
MatrixXd Q0(n/2, n/2), Q1(n/2, n/2), Q2(n/2, n/2), Q3(n/2, n/2),
Q4(n/2, n/2), Q5(n/2, n/2), Q6(n/2, n/2);
MatrixXd C(n,n);
// TO DO: (2-4.a) Finish Strassen's algorithm.
// Hint: Use strassenMatMult() to fill in the matrices Q0,..., Q6,
// and then use Q0,..., Q6 to fill in the matrix C.
// Note that comma-initialization works for block matrices.
// START:
MatrixXd A11 = A.topLeftCorner(n/2, n/2);
MatrixXd A12 = A.topRightCorner(n/2, n/2);
MatrixXd A21 = A.bottomLeftCorner(n/2, n/2);
MatrixXd A22 = A.bottomRightCorner(n/2, n/2);
MatrixXd B11 = B.topLeftCorner(n/2, n/2);
MatrixXd B12 = B.topRightCorner(n/2, n/2);
MatrixXd B21 = B.bottomLeftCorner(n/2, n/2);
MatrixXd B22 = B.bottomRightCorner(n/2, n/2);
Q0 = strassenMatMult(A11 + A22, B11 + B22);
Q1 = strassenMatMult(A21 + A22, B11);
Q2 = strassenMatMult(A11, B12 - B22);
Q3 = strassenMatMult(A22, B21 - B11);
Q4 = strassenMatMult(A11 + A12, B22);
Q5 = strassenMatMult(A21 - A11, B11 + B12);
Q6 = strassenMatMult(A12 - A22, B21 + B22);
C << Q0 + Q3 - Q4 + Q6,
Q2 + Q4,
Q1 + Q3,
Q0 + Q2 - Q1 + Q5;
// END
return C;
}
/* SAM_LISTING_END_0 */
/* SAM_LISTING_BEGIN_! */
double test_strassen() {
// Check algorithm for correctness
// Size of the matrix is a power of 2.
int k = 4;
int n = std::pow(2, k);
double err=1;
// TO DO: (2-4.b) Generate two random matrices using MatrixXd::Random(n,n).
// and compare the result of strassenMatMult() and the built-in matrix multiplication.
// Save the difference in the double err.
// START
MatrixXd A = MatrixXd::Random(n,n);
MatrixXd B = MatrixXd::Random(n,n);
// Testing matrix multiplication
MatrixXd AB = strassenMatMult(A,B);
// Eigen Matrix multiplication
MatrixXd AxB = A*B;
err = (AB-AxB).norm();
// END
return err;
}
/* SAM_LISTING_END_1 */
/* SAM_LISTING_END_2 */
void time_strassen() {
// Unfortunately, the time limit set by Code Expert may not be sufficiently
// long to show the advantage of Strassen's algorithm for large matrices.
// Minimum number of repetitions
unsigned int repetitions = 5;
// Display header column
std::cout << std::setw(4) << "k"
<< std::setw(15) << "A*B"
<< std::setw(15) << "Strassen" << std::endl;
for(unsigned k = 3; k <= 6; k++) {
unsigned int n = std::pow(2, k);
// Initialize random input matricies
MatrixXd A = MatrixXd::Random(n, n);
MatrixXd B = MatrixXd::Random(n, n);
// Timer to collect runtime of each individual run
Timer timer, timer_own;
// TO DO: Use the .start(), and .stop() methods of timer and timer_own to
// measure the runtimes of Eigen's built-in multiplication and Strassen's algorithm.
// Perform a few trials for each k (use the variable repetitions).
// START
for(unsigned int r = 0; r < repetitions; ++r) {
// initialize memory for result matrix
MatrixXd AxB, AxB2;
// Benchmark Eigen's matrix multiplication
timer.start(); // start timer
AxB=(A*B); // do the multiplication
timer.stop(); // stop timer
// Benchmark Strassen's matrix multiplication
timer_own.start(); // start timer
AxB2=strassenMatMult(A, B); // do the multiplication
timer_own.stop(); // stop timer
}
// END
// Print runtimes
std::cout << std::setw(4) << k // power
<< std::setprecision(3) << std::setw(15) << std::scientific
<< timer.min() // eigen timing
<< std::setprecision(3) << std::setw(15) << std::scientific
<< timer_own.min() // strassing timing
<< std::endl;
}
}
/* SAM_LISTING_END_2 */
\ No newline at end of file
// "testname", "input", "expected output", "points", "execution time limit (ms)"
"Identity test", "(any input)", "Works for identity matrix: 1 1", "1", "10000",
"Test Strassen", "(any input)", "Works for two random matrices: 1", "1", "10000",
# ifndef TIMER_HPP
# define TIMER_HPP
# include <iostream>
# include <chrono>
# include <vector>
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* USAGE: Timer t; *
* t.start(); *
* for (bla) { ... stuff happening ...; t.lap(); } *
* double min = t.min(), *
* mean = t.mean(); *
* *
* OR *
* *
* Timer t; *
* t.start(); *
* ... stuff happening ... *
* t.stop(); *
* *
* NOTE: stop() and lap() are equivalent! *
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
class Timer {
typedef std::chrono::nanoseconds prec;
typedef std::chrono::high_resolution_clock clock;
typedef std::chrono::duration<double> duration_t;
private:
static const unsigned divisor = 1e9;
clock::time_point t_start, t_end;
duration_t t_min;
std::vector<duration_t> t_laps;
public:
Timer();
void start();
void stop();
void lap();
void reset();
double duration() const;
double mean() const;
double min() const;
};
// start the timer
void Timer::start(){
t_start = clock::now();
t_end = t_start;
}
// stop the timer (equivalent to lap)
void Timer::stop(){
// stop is just another lap
lap();
}
// new lap
void Timer::lap(){
clock::time_point tmp = clock::now();
// get laptime
duration_t laptime;
laptime = tmp - t_end;
// check if this lap was faster
if (t_min > laptime || t_laps.size() == 0) {
t_min = laptime;
}
// save time of this lap
t_laps.push_back(laptime);
// save total time
t_end = tmp;
}
// idle constructor
Timer::Timer() {
#ifndef NDEBUG
static bool runonce = true;
if (runonce) {
std::cerr << "Warning: Timer was build as DEBUG." << std::endl;
runonce = false;
}
#endif
}
// resets all values
void Timer::reset(){
t_laps = std::vector<duration_t>();
start();
}
// returns total duration timer has been running
double Timer::duration() const {
if (t_laps.size() > 0){
// returning time in seconds! thats what the divisor is for
auto dur = std::chrono::duration_cast<prec>(t_end - t_start);
return double(dur.count())/divisor;
}
else {
std::cerr << "Before calling Timer::duration() you need to call Timer::lap() or Timer::stop()!\n";
return 0.;
}
}
// returns mean of all laps
double Timer::mean() const {
if (t_laps.size() > 0){
// save total time in std::chrono units
auto total_time = t_laps[0];
for (unsigned int i = 1; i < t_laps.size(); ++i) {
total_time += t_laps[i];
}
// convert time to double
auto total_dur = std::chrono::duration_cast<prec>(total_time);
double avg = double(total_dur.count())/divisor;
return avg;
}
else {
std::cerr << "Before calling Timer::mean() you need to call Timer::lap() or Timer::stop()!\n";
return 0.;
}
}
// returns minimum of all laps
double Timer::min() const {
auto min = std::chrono::duration_cast<prec>(t_min);
return double(min.count())/divisor;
}
# endif
This diff is collapsed.
# Code Expert Generic Project Configuration.
#
# In case this file is omitted, the execution is done with the default values.
# You can either define a cmd or a script file to be executed.
# Default to script file setting a cmd overrides it.
cmd:
# compile: echo "COMPILE"
# run: echo "RUN"
# test: echo "TEST"
# submit: echo "SUBMIT"
# All scripts are set by default to "/scripts/${ACTION}.sh"
script:
compile: /scripts/compile.sh
run: /scripts/run.sh
test: /scripts/test.sh
submit: /scripts/submit.sh
\ No newline at end of file
descriptions:
default:
name: Default
default: true
# Optionally add mulitple versions here.
# Use the file name without the .md extension as key and add a `name:` to it.
\ No newline at end of file
# Problem 2-4: Fast matrix multiplication with EIGEN
> For the task description of this exercise, please refer to [NCSE19_Problems.pdf](
https://www.sam.math.ethz.ch/~grsam/NCSE19/HOMEWORK/NCSE19_Problems.pdf).
> Open "strassen.hpp" and fill in the missing code in between the delimiters `// START` and `// END` according to the instructions preceded by `// TO DO:`.
# Tests of gramschmidt.hpp
> There are three functions in strassen.hpp and main.cpp has tests for the first two. The third function is a timing exercise, which should not be tested automatically.
***
> `strassenMatMult(const MatrixXd& A, const MatrixXd& B)`: We test this function by calculating norm(Strassen(A,Id) - A) and also norm(Strassen(Id,A) - A) for a random 128*128 matrix A. The tolerance is 10^{-9}.
***
> `test_strassen()`: This function is a test written by the student. It should create two random matrices A and B, and return err = norm(Strassen(A,B) - A*B).
The test is passed if err < 10^{-9}.
\ No newline at end of file
#include <Eigen/Dense>
#include <iostream>
#include "strassen.hpp"
using namespace Eigen;
int main() {
double tolerance = 1e-9;
int n = 128;
std::srand(5);
MatrixXd A = MatrixXd::Random(n,n);
std::cout << "Works for identity matrix: "
<< ((strassenMatMult(A, MatrixXd::Identity(n,n)) - A).norm()<tolerance)
<< " "
<< ((strassenMatMult(MatrixXd::Identity(n,n),A) - A).norm()<tolerance)
<< std::endl;
std::cout << "Works for two random matrices: "
<<(test_strassen() < tolerance) << std::endl;
std::cout << "Output of test_strassen()" << test_strassen() << std::endl;
time_strassen();
}
#!/bin/bash
echo "Compiling ..."
# compile c code
#mkdir -p bin
#g++ test_strassen.cpp -I/usr/include/eigen3/ -fdiagnostics-color=always --pedantic -Wextra -Wall -std=c++14 \
# -o bin/test_strassen.out
#g++ timing.cpp -I/usr/include/eigen3/ -fdiagnostics-color=always --pedantic -Wextra -Wall -std=c++14 \
# -o bin/timing.out
mkdir -p bin
find . -iname '*.cpp' | sort | xargs \
g++ -fdiagnostics-color=always -std=c++11 \
-I/usr/local/include/python3.7m \
-I/usr/local/lib/python3.7/site-packages/numpy/core/include \
-I/usr/include/eigen3/ \
-lpython3.7m \
-lpthread -lutil -ldl \
-Xlinker -export-dynamic \
-o bin/a.out
echo "Compilation successful"
\ No newline at end of file
#!/bin/bash
export PYTHONIOENCODING="UTF-8"
# compile (call compile script)
bash "${WORKDIR}/scripts/compile.sh"
echo "Running script ... "
# run
bin/a.out
\ No newline at end of file
#!/bin/bash
# I would like to delete this file, but I cannot.
export PYTHONIOENCODING="UTF-8"
# compile (call compile script)
bash "${WORKDIR}/scripts/compile.sh"
echo "Running code ..."
echo " "
echo " "
# run
echo "Task 1-4.b) 'Testing the correctness of the Strassen Alorithm' "
echo " "
bin/test_strassen.out
echo " "
echo "Task 1-4.c) 'Comparing the timings of A*B and the Strassen Algorithm' "
echo " "
bin/timing.out
\ No newline at end of file
#!/bin/bash
echo "Not implemented"
\ No newline at end of file
#!/bin/bash
# run (call run script)
bash "${WORKDIR}/scripts/run.sh"
\ No newline at end of file
////
//// Copyright (C) 2016 SAM (D-MATH) @ ETH Zurich
//// Author(s): lfilippo <filippo.leonardi@sam.math.ethz.ch>
//// Contributors: tille, jgacon, dcasati
//// This file is part of the NumCSE repository.
////
#include <Eigen/Dense>
#include <iomanip>
#include <iostream>
#include <vector>
#include "timer.h"
using namespace Eigen;
/* \brief Compute the Matrix product $A \times B$ using Strassen's algorithm.
* \param[in] A Matrix $2^k \times 2^k$
* \param[in] B Matrix $2^k \times 2^k$
* \param[out] Matrix product of A and B of dim $2^k \times 2^k$
*/
/* SAM_LISTING_BEGIN_0 */
MatrixXd strassenMatMult(const MatrixXd& A, const MatrixXd& B) {
// Ensure square matrix
assert(A.rows() == A.cols() && "Matrix A must be square");
assert(B.rows() == B.cols() && "Matrix B must be square");
// Matrix dimension must be a power of 2
assert(A.rows() % 2 == 0 && "Matrix dimensions must be a power of two.");
const unsigned n = A.rows();
// The function is recursive and acts on blocks of size $n/2 \times n/2$
// i.e. exploits fast product of 2x2 block matrix
if ( n==2 ) { // End of recursion
MatrixXd C(2, 2);
C << A(0,0)*B(0,0) + A(0,1)*B(1,0),
A(0,0)*B(0,1) + A(0,1)*B(1,1),
A(1,0)*B(0,0) + A(1,1)*B(1,0),
A(1,0)*B(0,1) + A(1,1)*B(1,1);
return C;
}
MatrixXd Q0(n/2, n/2), Q1(n/2, n/2), Q2(n/2, n/2), Q3(n/2, n/2),
Q4(n/2, n/2), Q5(n/2, n/2), Q6(n/2, n/2);
MatrixXd C(n,n);
// TO DO: (2-4.a) Finish Strassen's algorithm.
// Hint: Use strassenMatMult() to fill in the matrices Q0,..., Q6,
// and then use Q0,..., Q6 to fill in the matrix C.
// Note that comma-initialization works for block matrices.
// START:
// END
return C;
}
/* SAM_LISTING_END_0 */
/* SAM_LISTING_BEGIN_! */
double test_strassen() {
// Check algorithm for correctness
// Size of the matrix is a power of 2.
int k = 4;
int n = std::pow(2, k);
double err=1;
// TO DO: (2-4.b) Generate two random matrices using MatrixXd::Random(n,n).
// and compare the result of strassenMatMult() and the built-in matrix multiplication.
// Save the difference in the double err.
// START
// END
return err;