added lab03
This commit is contained in:
WickedJack99
2
.gitignore
vendored
2
.gitignore
vendored
@@ -0,0 +1,2 @@
|
||||
*.out*
|
||||
execs
|
||||
3
.vscode/settings.json
vendored
3
.vscode/settings.json
vendored
@@ -58,6 +58,7 @@
|
||||
"xtr1common": "cpp",
|
||||
"xutility": "cpp",
|
||||
"stop_token": "cpp",
|
||||
"thread": "cpp"
|
||||
"thread": "cpp",
|
||||
"array": "cpp"
|
||||
}
|
||||
}
|
||||
18
lab02/calc_lookup_table.py
Normal file
18
lab02/calc_lookup_table.py
Normal file
@@ -0,0 +1,18 @@
|
||||
import math
|
||||
|
||||
def calc():
|
||||
with open("output.txt", "w") as f:
|
||||
values = ""
|
||||
for i in range(0,256):
|
||||
value = fib(i)
|
||||
val = (math.sin(value) * math.tan(value) / math.sqrt(math.cos(value) + 2))
|
||||
values += f"{val}, "
|
||||
f.write(values)
|
||||
|
||||
def fib(n):
|
||||
a, b = 0, 1
|
||||
for _ in range(n):
|
||||
a, b = b, a + b
|
||||
return a
|
||||
|
||||
calc()
|
||||
File diff suppressed because one or more lines are too long
@@ -8,15 +8,17 @@ def read_csv_pandas(filename):
|
||||
return df.values.flatten().astype(float) # Konvertiert in flaches NumPy-Array (float, um Fehler zu vermeiden)
|
||||
|
||||
# Daten einlesen
|
||||
data = read_csv_pandas("outcomes-100-1000000.csv")
|
||||
data = read_csv_pandas("outcomes5.csv")
|
||||
print(str(len(data)))
|
||||
|
||||
# Plot erstellen
|
||||
plt.figure(figsize=(8, 5))
|
||||
plt.plot(data, marker="o", linestyle="-", color="b", linewidth=1, markersize=1)
|
||||
|
||||
# Styling für einen coolen Look
|
||||
plt.xlabel("n", fontsize=12)
|
||||
plt.ylabel("ms", fontsize=12)
|
||||
plt.xlabel("n/10", fontsize=12)
|
||||
plt.ylabel("p(FLOPS)", fontsize=12)
|
||||
plt.xscale('log')
|
||||
plt.grid(True, linestyle="--", alpha=0.6)
|
||||
plt.legend()
|
||||
plt.show()
|
||||
|
||||
1
lab02/outcomes2.csv
Normal file
1
lab02/outcomes2.csv
Normal file
File diff suppressed because one or more lines are too long
1
lab02/outcomes3.csv
Normal file
1
lab02/outcomes3.csv
Normal file
File diff suppressed because one or more lines are too long
1
lab02/outcomes5.csv
Normal file
1
lab02/outcomes5.csv
Normal file
File diff suppressed because one or more lines are too long
1
lab02/output.txt
Normal file
1
lab02/output.txt
Normal file
File diff suppressed because one or more lines are too long
@@ -26,13 +26,32 @@ static double ownMethod(int n) {
|
||||
return averageTime;
|
||||
}
|
||||
|
||||
static double ownMethod2(int n) {
|
||||
std::vector<double> A(n, 0);
|
||||
std::chrono::duration<double, std::milli> sumTime;
|
||||
for (int i = 0; i < 20; ++i) {
|
||||
auto startTime = std::chrono::system_clock::now();
|
||||
for (int j = 0; j < n; ++j) {
|
||||
A[j] = B[j] + C[j] * D[j];
|
||||
}
|
||||
// prevent the compiler from optimizing everything away
|
||||
volatile double dummy = A[0];
|
||||
auto endTime = std::chrono::system_clock::now();
|
||||
sumTime += endTime - startTime;
|
||||
}
|
||||
double averageTime = sumTime.count() / 20;
|
||||
double flops = (2.0 * n) / (averageTime / 1000);
|
||||
|
||||
return flops;
|
||||
}
|
||||
|
||||
static void workerThread(int start, int end) {
|
||||
std::vector<double> outcomes;
|
||||
for (int i = start; i < end; i+=10) {
|
||||
outcomes.push_back(ownMethod(i));
|
||||
outcomes.push_back(ownMethod2(i));
|
||||
}
|
||||
std::stringstream fileName;
|
||||
fileName << "outcomes-" << start << "-" << end << ".csv";
|
||||
fileName << "outcomes-" << start << "-" << end << "2.csv";
|
||||
std::ofstream MyFile(fileName.str());
|
||||
for (int j = 0; j < outcomes.size(); j++) {
|
||||
MyFile << outcomes[j] << ",";
|
||||
|
||||
71
lab03/.vscode/settings.json
vendored
Normal file
71
lab03/.vscode/settings.json
vendored
Normal file
@@ -0,0 +1,71 @@
|
||||
{
|
||||
"files.associations": {
|
||||
"chrono": "cpp",
|
||||
"cctype": "cpp",
|
||||
"clocale": "cpp",
|
||||
"cmath": "cpp",
|
||||
"cstdarg": "cpp",
|
||||
"cstddef": "cpp",
|
||||
"cstdio": "cpp",
|
||||
"cstdlib": "cpp",
|
||||
"cstring": "cpp",
|
||||
"ctime": "cpp",
|
||||
"cwchar": "cpp",
|
||||
"cwctype": "cpp",
|
||||
"array": "cpp",
|
||||
"atomic": "cpp",
|
||||
"bit": "cpp",
|
||||
"*.tcc": "cpp",
|
||||
"bitset": "cpp",
|
||||
"compare": "cpp",
|
||||
"complex": "cpp",
|
||||
"concepts": "cpp",
|
||||
"condition_variable": "cpp",
|
||||
"cstdint": "cpp",
|
||||
"deque": "cpp",
|
||||
"forward_list": "cpp",
|
||||
"list": "cpp",
|
||||
"map": "cpp",
|
||||
"set": "cpp",
|
||||
"string": "cpp",
|
||||
"unordered_map": "cpp",
|
||||
"unordered_set": "cpp",
|
||||
"vector": "cpp",
|
||||
"exception": "cpp",
|
||||
"algorithm": "cpp",
|
||||
"functional": "cpp",
|
||||
"iterator": "cpp",
|
||||
"memory": "cpp",
|
||||
"memory_resource": "cpp",
|
||||
"numeric": "cpp",
|
||||
"optional": "cpp",
|
||||
"random": "cpp",
|
||||
"ratio": "cpp",
|
||||
"string_view": "cpp",
|
||||
"system_error": "cpp",
|
||||
"tuple": "cpp",
|
||||
"type_traits": "cpp",
|
||||
"utility": "cpp",
|
||||
"fstream": "cpp",
|
||||
"initializer_list": "cpp",
|
||||
"iomanip": "cpp",
|
||||
"iosfwd": "cpp",
|
||||
"iostream": "cpp",
|
||||
"istream": "cpp",
|
||||
"limits": "cpp",
|
||||
"mutex": "cpp",
|
||||
"new": "cpp",
|
||||
"numbers": "cpp",
|
||||
"ostream": "cpp",
|
||||
"semaphore": "cpp",
|
||||
"span": "cpp",
|
||||
"sstream": "cpp",
|
||||
"stdexcept": "cpp",
|
||||
"stop_token": "cpp",
|
||||
"streambuf": "cpp",
|
||||
"thread": "cpp",
|
||||
"cinttypes": "cpp",
|
||||
"typeindex": "cpp",
|
||||
"typeinfo": "cpp"
|
||||
}
|
||||
}
|
||||
348
lab03/compute2.cpp
Normal file
348
lab03/compute2.cpp
Normal file
@@ -0,0 +1,348 @@
|
||||
#include <chrono>
|
||||
#include <cmath>
|
||||
#include <iostream>
|
||||
#include <vector>
|
||||
#include "matrix.h"
|
||||
|
||||
int fib2(int i) {
|
||||
if (i == 0) {
|
||||
return 0;
|
||||
} else if (i == 1) {
|
||||
return 1;
|
||||
} else {
|
||||
return fib2(i - 1) + fib2(i - 2);
|
||||
}
|
||||
}
|
||||
Matrix compute2(const Matrix& s, const std::vector<int>& v) {
|
||||
Matrix m(s.dim());
|
||||
const int n = v.size();
|
||||
for (int i = 0; i < n; ++i) {
|
||||
for (int j = 0; j < n; ++j) {
|
||||
double val = static_cast<double>(fib2(v[j] % 256));
|
||||
m(j, i) = s(j, i) * (sin(val) * tan(val) / sqrt(cos(val) + 2));
|
||||
}
|
||||
}
|
||||
return m;
|
||||
}
|
||||
|
||||
void print(Matrix& m) {
|
||||
std::cout << std::endl;
|
||||
|
||||
for (int i = 0; i < 3; i++) {
|
||||
for (int y = 0; y < 3; y++) {
|
||||
std::cout << m(i, y) << "\t";
|
||||
}
|
||||
std::cout << std::endl;
|
||||
}
|
||||
std::cout << std::endl;
|
||||
}
|
||||
|
||||
Matrix create_matrix(int x, int n) {
|
||||
Matrix m = Matrix(n, n);
|
||||
|
||||
for (int i = 0; i < n; i++) {
|
||||
for (int y = 0; y < n; y++) {
|
||||
m(i, y) = x;
|
||||
}
|
||||
}
|
||||
|
||||
return m;
|
||||
}
|
||||
|
||||
std::vector<int> create_vector(int x, int n) {
|
||||
std::vector v = std::vector<int>(n);
|
||||
for (int i = 0; i < n; i++) {
|
||||
v[i] = x;
|
||||
}
|
||||
return v;
|
||||
}
|
||||
|
||||
constexpr std::array<double, 256> lookupTable = {0.0,
|
||||
0.8222403186860677,
|
||||
0.8222403186860677,
|
||||
-1.5787300242524172,
|
||||
-0.020016262234217718,
|
||||
2.1451183132222993,
|
||||
-4.940042066178768,
|
||||
0.11409507657148513,
|
||||
-1.0604830926169537,
|
||||
-0.3074259803136817,
|
||||
31.766208383302768,
|
||||
0.9151534216757803,
|
||||
0.16333595403676757,
|
||||
0.16949512185916804,
|
||||
4.5543599070187406e-05,
|
||||
0.1757957713418539,
|
||||
0.18224021532829443,
|
||||
1.0206149092922863,
|
||||
-17.814542276286655,
|
||||
-0.20918076066994673,
|
||||
-1.6213971271904437,
|
||||
0.3168557058165431,
|
||||
2.298722897010151,
|
||||
0.23686197492476985,
|
||||
-2.086362829745931,
|
||||
-0.3881868942002799,
|
||||
-2.461117527398634,
|
||||
0.2808435566118653,
|
||||
1.5819909021264011,
|
||||
0.18830714860806283,
|
||||
-3.9783311640939956,
|
||||
-0.7223836412497261,
|
||||
-0.8639724628881527,
|
||||
0.0027677492169528905,
|
||||
-1.0303485832123693,
|
||||
-1.230575195462486,
|
||||
-1.45754407514109,
|
||||
0.0023462043707461705,
|
||||
-1.7453680716455173,
|
||||
-2.1290409071270853,
|
||||
-0.5251182747080108,
|
||||
0.16466950650652515,
|
||||
-0.04385364901718776,
|
||||
-0.09269399902874279,
|
||||
0.0052966445389224254,
|
||||
-0.15938094474318504,
|
||||
-0.2440476196026588,
|
||||
0.6058083982842931,
|
||||
-4.095160159163112,
|
||||
0.5228302544273172,
|
||||
-0.2944149703898431,
|
||||
-0.0914706601652273,
|
||||
0.0343847919919084,
|
||||
-0.00374146357960608,
|
||||
-0.03269530501366555,
|
||||
0.008326304813702903,
|
||||
-0.08998858152591463,
|
||||
-0.1753231596613382,
|
||||
0.35059963365479513,
|
||||
-1.582851175470022,
|
||||
2.093789103171035,
|
||||
-0.0181535233234956,
|
||||
-4.548501978245147,
|
||||
29.55052845809529,
|
||||
0.01936746685187523,
|
||||
3.1521208201441255,
|
||||
1.4657757926693418,
|
||||
-0.35010475431933086,
|
||||
-0.34162920251204093,
|
||||
2.9770086129136403e-05,
|
||||
-0.33326757588611455,
|
||||
-0.325018905514724,
|
||||
1.3297039478954402,
|
||||
4.2812466550630734,
|
||||
0.040135899514659124,
|
||||
-6.823980420549866,
|
||||
-1.896232407731212,
|
||||
-0.22067093110457817,
|
||||
0.362095020184007,
|
||||
-0.5492657903414403,
|
||||
-8.263813069231747e-05,
|
||||
0.3839613240456214,
|
||||
0.005732116510943663,
|
||||
1.0068069901920167,
|
||||
-0.06661408684607055,
|
||||
0.16505937587347172,
|
||||
-0.10618985842084362,
|
||||
-0.1548539438950448,
|
||||
0.35134836733238417,
|
||||
-1.9850782320522917,
|
||||
2.3201973898487434,
|
||||
-0.8614531626702128,
|
||||
0.6680745565561412,
|
||||
-0.00044550347885444647,
|
||||
-0.39763818064051737,
|
||||
0.7578655614877392,
|
||||
0.6072237200360099,
|
||||
0.003178329575514224,
|
||||
-1.212826954378079,
|
||||
-1.4769326433996612,
|
||||
-172.52806598178202,
|
||||
-0.5347454349381718,
|
||||
-19.190923961699053,
|
||||
-0.7917514182651304,
|
||||
0.2819302722074698,
|
||||
0.05458474508581716,
|
||||
0.019872869052966566,
|
||||
0.008222147618604445,
|
||||
0.09756055641361473,
|
||||
0.7227958766703475,
|
||||
-1.4410095077924745,
|
||||
-9.972144169933975,
|
||||
0.77699771545969,
|
||||
-4.53513322683293,
|
||||
-0.18352603012839258,
|
||||
0.0017072712007560583,
|
||||
-0.23257515516939736,
|
||||
-0.10763881722478827,
|
||||
0.884626439996002,
|
||||
-0.4810037394118251,
|
||||
-1.369797203450087,
|
||||
-0.5288521774897473,
|
||||
-2.9101872323475186,
|
||||
0.22727177754107092,
|
||||
-0.34542092073980846,
|
||||
-1.8050778597271122,
|
||||
-2.619459073093681,
|
||||
0.20132946664538756,
|
||||
0.26819837253349504,
|
||||
1.762278187625015,
|
||||
-0.011986063261174138,
|
||||
0.24469190367516808,
|
||||
1.0086060475978966,
|
||||
-3.165218395481177,
|
||||
4.189103518135543,
|
||||
-2.5760814095593942,
|
||||
-0.013859290579366626,
|
||||
1.4352697196307378,
|
||||
-1.2608389085452532,
|
||||
-4.198328999769481,
|
||||
0.1439300331481971,
|
||||
0.34262172197897917,
|
||||
0.0168015128985714,
|
||||
0.18580678551855231,
|
||||
1.4365641854331133,
|
||||
-4.79841926852741,
|
||||
-3.4019533078572475,
|
||||
-0.01590191501468279,
|
||||
1.759259558519459,
|
||||
-1.0228945308454913,
|
||||
0.29676927098386374,
|
||||
-0.09248853701134174,
|
||||
-0.8115494202333414,
|
||||
0.21410452959908136,
|
||||
-6.956781842218232,
|
||||
1.0938240947751094,
|
||||
8.314261325910216,
|
||||
0.09374743357102927,
|
||||
0.07071517812352679,
|
||||
-3.4713782052240494,
|
||||
0.051182474782728866,
|
||||
-1.2962746253140673,
|
||||
2.204810869689383,
|
||||
0.42653001839592647,
|
||||
0.15888109869643205,
|
||||
-2.2502390908963252,
|
||||
-0.5636032756994527,
|
||||
0.0012139954702342257,
|
||||
-161.44964948195522,
|
||||
0.0022234794467720603,
|
||||
-0.002981746927564041,
|
||||
-0.013589893337736567,
|
||||
-0.02176301438843992,
|
||||
0.041053147217568434,
|
||||
-0.2535229180317107,
|
||||
-0.33498174636117795,
|
||||
1.0769953736629612,
|
||||
6.899966717364552,
|
||||
0.0827438107746741,
|
||||
0.07067578121842849,
|
||||
7.848050576848035,
|
||||
-2.7861390728504825,
|
||||
0.0013635934995610396,
|
||||
-1.153162937055766,
|
||||
-6.3915800646950816,
|
||||
-0.031171002889077308,
|
||||
0.27503561575602464,
|
||||
-0.694782893937455,
|
||||
-0.2045488357018786,
|
||||
-1.3814287827861556,
|
||||
0.2639397999144297,
|
||||
-1.6283372950295938,
|
||||
0.03857424446734361,
|
||||
0.3272148908804416,
|
||||
0.21928953507794532,
|
||||
1.2939133941195287,
|
||||
-0.26723965240335223,
|
||||
0.0777139485862296,
|
||||
-1.37710795567683,
|
||||
-2.1373511044382387,
|
||||
-0.42382559169444095,
|
||||
0.21966996134691846,
|
||||
1.5086628407480521e-05,
|
||||
0.22392450662682784,
|
||||
0.12739930332047994,
|
||||
1.511939818848616,
|
||||
3.1675580990874033,
|
||||
-0.33854805156070844,
|
||||
-3.193060131852421,
|
||||
0.08157921242282763,
|
||||
0.005207776703273675,
|
||||
-1.644268309315312,
|
||||
-0.4910608493321238,
|
||||
0.12568653862630152,
|
||||
-1.1324160197145394,
|
||||
-0.042293254578116724,
|
||||
0.3565189224930081,
|
||||
0.8199298840209298,
|
||||
2.9100467396337413,
|
||||
0.14328425483708107,
|
||||
-2.7351863856444307,
|
||||
-0.3118936323791283,
|
||||
-2.3652167108089786,
|
||||
0.3364739591631545,
|
||||
0.00601840488067479,
|
||||
11.828244785427081,
|
||||
0.1124143840824127,
|
||||
1.0545971966119048,
|
||||
0.2759813103079656,
|
||||
-0.23352307906228556,
|
||||
-1.5576719903385443,
|
||||
-13.76407981937769,
|
||||
0.08754583137498555,
|
||||
1.810930803449213,
|
||||
0.5382555408863076,
|
||||
-0.002576694765818232,
|
||||
-0.8612366710787147,
|
||||
-0.07277140430056403,
|
||||
0.26460751353648865,
|
||||
0.057583427346797084,
|
||||
-6.1626683426032765,
|
||||
-0.16257003129521286,
|
||||
-5.460375467613629,
|
||||
-2.6137508557003617,
|
||||
0.7780098599544784,
|
||||
1.2899685286178642,
|
||||
-1.1290019883416802,
|
||||
-0.0969477560271896,
|
||||
2.6188550083461326,
|
||||
9.939111102377229,
|
||||
-1.8964522816616587,
|
||||
-6.9572063958158665,
|
||||
-0.21874194287663265,
|
||||
0.8151204126824054,
|
||||
-0.2744708906370598,
|
||||
-0.0006656271583815087};
|
||||
|
||||
Matrix computeNew(const Matrix& s, const std::vector<int>& v) {
|
||||
Matrix m(s.dim());
|
||||
const int n = v.size();
|
||||
for (int j = 0; j < n; ++j) {
|
||||
for (int i = 0; i < n; ++i) {
|
||||
m(j, i) = s(j, i) * lookupTable[v[j] % 256];
|
||||
}
|
||||
}
|
||||
return m;
|
||||
}
|
||||
|
||||
int main() {
|
||||
std::vector<int> v = create_vector(1, 1000);
|
||||
Matrix m = create_matrix(1, 1000);
|
||||
|
||||
std::chrono::duration<double, std::milli> sumTimeOld;
|
||||
std::chrono::duration<double, std::milli> sumTimeNew;
|
||||
auto startTimeOld = std::chrono::system_clock::now();
|
||||
volatile auto erg1 = compute2(m, v);
|
||||
auto endTimeOld = std::chrono::system_clock::now();
|
||||
sumTimeOld = endTimeOld - startTimeOld;
|
||||
|
||||
auto startTimeNew = std::chrono::system_clock::now();
|
||||
volatile auto erg2 = computeNew(m, v);
|
||||
auto endTimeNew = std::chrono::system_clock::now();
|
||||
sumTimeNew = endTimeNew - startTimeNew;
|
||||
|
||||
std::cout << "Time old: " << sumTimeOld.count()
|
||||
<< "\nTime new: " << sumTimeNew.count() << std::endl;
|
||||
|
||||
return 0;
|
||||
}
|
||||
81
lab03/jacobi_bench1.cpp
Normal file
81
lab03/jacobi_bench1.cpp
Normal file
@@ -0,0 +1,81 @@
|
||||
#include "matrix.h"
|
||||
#include <iostream>
|
||||
#include <fstream>
|
||||
#include <benchmark/benchmark.h>
|
||||
|
||||
Matrix jacobi(const Matrix& init, double eps, int maxNumIter) {
|
||||
Matrix phi(init);
|
||||
Matrix deltaPhi(phi.dim());
|
||||
const int n = phi.dim1();
|
||||
const int m = phi.dim2();
|
||||
|
||||
double dist = 1e10;
|
||||
int nIter = 0;
|
||||
while (dist > eps && nIter < maxNumIter) {
|
||||
for (int i = 1; i < n - 1; ++i) {
|
||||
for (int j = 1; j < m - 1; ++j) {
|
||||
deltaPhi(i, j) = .25 * (phi(i + 1, j) + phi(i - 1, j) + phi(i, j + 1) +
|
||||
phi(i, j - 1)) -
|
||||
phi(i, j);
|
||||
}
|
||||
}
|
||||
|
||||
dist = 0;
|
||||
for (int i = 1; i < n - 1; ++i) {
|
||||
for (int j = 1; j < n - 1; ++j) {
|
||||
dist = std::max(dist, std::abs(deltaPhi(i, j)));
|
||||
phi(i, j) += deltaPhi(i, j);
|
||||
}
|
||||
}
|
||||
|
||||
nIter++;
|
||||
}
|
||||
|
||||
std::cout << "Finished Jacobi after " << nIter
|
||||
<< " iterations with an error of " << dist << std::endl;
|
||||
return phi;
|
||||
}
|
||||
|
||||
Matrix initialCondition(int n) {
|
||||
Matrix m(n, n);
|
||||
for (int i = 0; i < n; ++i) {
|
||||
const double x = 1. * i / n;
|
||||
|
||||
for (int j = 0; j < n; ++j) {
|
||||
const double y = 1. * j / n;
|
||||
|
||||
// The lower left corner is hot
|
||||
m(i, j) += exp(-0.5*(x * x + y * y) * 10) * exp(-100*x*y);
|
||||
|
||||
// The right side is cold
|
||||
m(i, j) += -1.0 * exp(-20 * (1 - y));
|
||||
}
|
||||
}
|
||||
return m;
|
||||
}
|
||||
|
||||
void storeMatrix(const Matrix& m, const std::string& filename) {
|
||||
std::ofstream fout(filename);
|
||||
fout << m << std::endl;
|
||||
}
|
||||
|
||||
// int main() {
|
||||
// const Matrix init = initialCondition(200);
|
||||
// storeMatrix(init, "init.asc");
|
||||
// Matrix phi = jacobi(init, 1e-5, 100000);
|
||||
// storeMatrix(phi, "solution.asc");
|
||||
// }
|
||||
|
||||
void benchmarkJacobi(benchmark::State& state) {
|
||||
const Matrix init = initialCondition(100);
|
||||
|
||||
Matrix phi;
|
||||
|
||||
for (auto _ : state) {
|
||||
phi = jacobi(init, 1e-5, 100000);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
BENCHMARK(benchmarkJacobi)->Unit(benchmark::kMillisecond);
|
||||
BENCHMARK_MAIN();
|
||||
122
lab03/jacobi_bench2.cpp
Normal file
122
lab03/jacobi_bench2.cpp
Normal file
@@ -0,0 +1,122 @@
|
||||
#include <benchmark/benchmark.h>
|
||||
#include <fstream>
|
||||
#include <iostream>
|
||||
#include "matrix.h"
|
||||
|
||||
Matrix jacobi(const Matrix& init, double eps, int maxNumIter) {
|
||||
Matrix *phi = new Matrix(init);
|
||||
Matrix *tmp = new Matrix(init);
|
||||
const int n = (*phi).dim1();
|
||||
const int m = (*phi).dim2();
|
||||
|
||||
double dist = 1e10;
|
||||
int nIter = 0;
|
||||
while (dist > eps && nIter < maxNumIter) {
|
||||
dist = 0;
|
||||
|
||||
for (int i = 1; i < n - 1; ++i) {
|
||||
for (int j = 1; j < m - 1; ++j) {
|
||||
// if (nIter % 2 == 0) {
|
||||
(*tmp)(i, j) = .25 * ((*phi)(i + 1, j) + (*phi)(i - 1, j) +
|
||||
(*phi)(i, j + 1) + (*phi)(i, j - 1));
|
||||
// } else {
|
||||
// (phi)(i, j) = .25 * ((tmp)(i + 1, j) + (tmp)(i - 1, j) +
|
||||
// (tmp)(i, j + 1) + (tmp)(i, j - 1));
|
||||
// }
|
||||
dist = std::max(dist, std::abs((*tmp)(i, j) - (*phi)(i, j)));
|
||||
}
|
||||
}
|
||||
// tmp = Matrix(phi);
|
||||
std::swap(tmp, phi);
|
||||
// for (int i = 1; i < n - 1; ++i) {
|
||||
// for (int j = 1; j < n - 1; ++j) {
|
||||
// phi(i, j) = tmp(i,j);
|
||||
// }
|
||||
// }
|
||||
|
||||
// std::cout << l2phi << ", " << l2tmp << ", " << dist << std::endl;
|
||||
|
||||
nIter++;
|
||||
}
|
||||
|
||||
std::cout << "Finished Jacobi after " << nIter
|
||||
<< " iterations with an error of " << dist << std::endl;
|
||||
return (*phi);
|
||||
}
|
||||
|
||||
Matrix jacobi1(const Matrix& init, double eps, int maxNumIter) {
|
||||
Matrix phi(init);
|
||||
Matrix tmp(phi.dim());
|
||||
const int n = phi.dim1();
|
||||
const int m = phi.dim2();
|
||||
|
||||
double dist = 1e10;
|
||||
int nIter = 0;
|
||||
while (dist > eps && nIter < maxNumIter) {
|
||||
dist = 0;
|
||||
for (int i = 1; i < n - 1; ++i) {
|
||||
for (int j = 1; j < m - 1; ++j) {
|
||||
tmp(i, j) = .25 * (phi(i + 1, j) + phi(i - 1, j) + phi(i, j + 1) + phi(i, j - 1));
|
||||
dist = std::max(dist, std::abs(tmp(i, j) - phi(i, j)));
|
||||
}
|
||||
}
|
||||
|
||||
for (int i = 1; i < n - 1; ++i) {
|
||||
for (int j = 1; j < n - 1; ++j) {
|
||||
phi(i, j) = tmp(i,j);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// std::cout << l2phi << ", " << l2tmp << ", " << dist << std::endl;
|
||||
|
||||
nIter++;
|
||||
}
|
||||
|
||||
std::cout << "Finished Jacobi after " << nIter
|
||||
<< " iterations with an error of " << dist << std::endl;
|
||||
return phi;
|
||||
}
|
||||
|
||||
Matrix initialCondition(int n) {
|
||||
Matrix m(n, n);
|
||||
for (int i = 0; i < n; ++i) {
|
||||
const double x = 1. * i / n;
|
||||
|
||||
for (int j = 0; j < n; ++j) {
|
||||
const double y = 1. * j / n;
|
||||
|
||||
// The lower left corner is hot
|
||||
m(i, j) += exp(-0.5 * (x * x + y * y) * 10) * exp(-100 * x * y);
|
||||
|
||||
// The right side is cold
|
||||
m(i, j) += -1.0 * exp(-20 * (1 - y));
|
||||
}
|
||||
}
|
||||
return m;
|
||||
}
|
||||
|
||||
void storeMatrix(const Matrix& m, const std::string& filename) {
|
||||
std::ofstream fout(filename);
|
||||
fout << m << std::endl;
|
||||
}
|
||||
|
||||
// int main() {
|
||||
// const Matrix init = initialCondition(200);
|
||||
// storeMatrix(init, "init.asc");
|
||||
// Matrix phi = jacobi(init, 1e-5, 100000);
|
||||
// storeMatrix(phi, "solution.asc");
|
||||
// }
|
||||
|
||||
void benchmarkJacobi(benchmark::State& state) {
|
||||
const Matrix init = initialCondition(1000);
|
||||
|
||||
Matrix phi;
|
||||
|
||||
for (auto _ : state) {
|
||||
phi = jacobi1(init, 1e-5, 100000);
|
||||
}
|
||||
}
|
||||
|
||||
BENCHMARK(benchmarkJacobi)->Unit(benchmark::kMillisecond);
|
||||
BENCHMARK_MAIN();
|
||||
340
lab03/matrix.h
Normal file
340
lab03/matrix.h
Normal file
@@ -0,0 +1,340 @@
|
||||
/**
|
||||
* matrix.h a very simplistic class for m times n matrices.
|
||||
*/
|
||||
|
||||
#ifndef MATRIX_H
|
||||
#define MATRIX_H
|
||||
|
||||
#include <cmath>
|
||||
#include <iomanip>
|
||||
#include <iostream>
|
||||
#include <vector>
|
||||
|
||||
// A very simplistic vector class for vectors of size n
|
||||
class Vector {
|
||||
public:
|
||||
// constructors
|
||||
Vector(int n) : n_(n), data_(n_, 0) {}
|
||||
Vector(const Vector &other) = default;
|
||||
Vector(Vector &&other) = default;
|
||||
~Vector() = default;
|
||||
|
||||
// assignment operators
|
||||
Vector &operator=(const Vector &other) = default;
|
||||
Vector &operator=(Vector &&other) = default;
|
||||
|
||||
// element access
|
||||
double &operator()(int i) { return data_[i]; }
|
||||
const double &operator()(int i) const { return data_[i]; }
|
||||
|
||||
// getter functions for the dimensions
|
||||
int dim() const { return n_; }
|
||||
|
||||
// comparison operators
|
||||
bool operator==(const Vector &b) { return (data_ == b.data_); }
|
||||
bool operator!=(const Vector &b) { return (data_ != b.data_); }
|
||||
|
||||
// addition
|
||||
Vector &operator+=(const Vector &b) {
|
||||
for (int i = 0; i < n_; ++i) {
|
||||
operator()(i) += b(i);
|
||||
}
|
||||
return *this;
|
||||
}
|
||||
|
||||
// subtraction
|
||||
Vector &operator-=(const Vector &b) {
|
||||
for (int i = 0; i < n_; ++i) {
|
||||
operator()(i) -= b(i);
|
||||
}
|
||||
return *this;
|
||||
}
|
||||
|
||||
// scalar multiplication
|
||||
Vector &operator*=(double x) {
|
||||
for (int i = 0; i < n_; ++i) {
|
||||
operator()(i) *= x;
|
||||
}
|
||||
return *this;
|
||||
}
|
||||
|
||||
// dot product between two vectors
|
||||
double dot(const Vector &other) const {
|
||||
double sum = 0;
|
||||
for (int i = 0; i < n_; ++i) {
|
||||
sum += operator()(i) * other(i);
|
||||
}
|
||||
return sum;
|
||||
}
|
||||
|
||||
private:
|
||||
int n_; // vector dimension
|
||||
std::vector<double> data_; // the vectors entries
|
||||
};
|
||||
|
||||
inline double dot(const Vector &v1, const Vector &v2) { return v1.dot(v2); }
|
||||
|
||||
// Print the vector as a table
|
||||
inline std::ostream &operator<<(std::ostream &os, const Vector &a) {
|
||||
const int width = 10;
|
||||
const int precision = 4;
|
||||
|
||||
const auto originalPrecision = os.precision();
|
||||
os << std::setprecision(precision);
|
||||
|
||||
for (int i = 0; i < a.dim(); ++i) {
|
||||
os << std::setw(width) << a(i) << " ";
|
||||
}
|
||||
|
||||
os << "\n";
|
||||
|
||||
os << std::setprecision(originalPrecision);
|
||||
return os;
|
||||
}
|
||||
|
||||
// A very simple class for m times n matrices
|
||||
class Matrix {
|
||||
public:
|
||||
// constructors
|
||||
Matrix() : Matrix(0, 0) {}
|
||||
Matrix(int m, int n) : m_(m), n_(n), data_(m_ * n_, 0) {}
|
||||
Matrix(std::pair<int, int> dim) : Matrix(dim.first, dim.second) {}
|
||||
Matrix(int n) : Matrix(n, n) {}
|
||||
Matrix(const Matrix &other) = default;
|
||||
Matrix(Matrix &&other) = default;
|
||||
~Matrix() = default;
|
||||
|
||||
// assignment operators
|
||||
Matrix &operator=(const Matrix &other) = default;
|
||||
Matrix &operator=(Matrix &&other) = default;
|
||||
|
||||
// element access
|
||||
double &operator()(int i, int j) { return data_[i * n_ + j]; }
|
||||
const double &operator()(int i, int j) const { return data_[i * n_ + j]; }
|
||||
|
||||
// getter functions for the dimensions
|
||||
std::pair<int, int> dim() const { return std::pair<int, int>(m_, n_); }
|
||||
int dim1() const { return m_; }
|
||||
int dim2() const { return n_; }
|
||||
int numEntries() const { return data_.size(); }
|
||||
|
||||
// comparison operators
|
||||
bool operator==(const Matrix &b) { return (data_ == b.data_); }
|
||||
bool operator!=(const Matrix &b) { return (data_ != b.data_); }
|
||||
|
||||
// addition
|
||||
Matrix &operator+=(const Matrix &b) {
|
||||
for (int i = 0; i < m_; ++i) {
|
||||
for (int j = 0; j < n_; ++j) {
|
||||
operator()(i, j) += b(i, j);
|
||||
}
|
||||
}
|
||||
return *this;
|
||||
}
|
||||
|
||||
// subtraction
|
||||
Matrix &operator-=(const Matrix &b) {
|
||||
for (int i = 0; i < m_; ++i) {
|
||||
for (int j = 0; j < n_; ++j) {
|
||||
operator()(i, j) -= b(i, j);
|
||||
}
|
||||
}
|
||||
return *this;
|
||||
}
|
||||
|
||||
// scalar multiplication
|
||||
Matrix &operator*=(double x) {
|
||||
for (int i = 0; i < m_; ++i) {
|
||||
for (int j = 0; j < n_; ++j) {
|
||||
operator()(i, j) *= x;
|
||||
}
|
||||
}
|
||||
return *this;
|
||||
}
|
||||
|
||||
// scalar division
|
||||
Matrix &operator/=(double x) {
|
||||
for (int i = 0; i < m_; ++i) {
|
||||
for (int j = 0; j < n_; ++j) {
|
||||
operator()(i, j) /= x;
|
||||
}
|
||||
}
|
||||
return *this;
|
||||
}
|
||||
|
||||
// matrix product (only for square matrices of equal dimension)
|
||||
Matrix &operator*=(const Matrix &b) {
|
||||
if (dim1() != dim2()) {
|
||||
std::cout << "Error in matrix multiplication: no square matrix\n";
|
||||
} else if (dim1() != b.dim1() || dim2() != b.dim2()) {
|
||||
std::cout << "Error in matrix multiplication: dimensions do not match\n";
|
||||
} else {
|
||||
Matrix a = *this;
|
||||
Matrix &c = *this;
|
||||
const int m = dim1();
|
||||
for (int i = 0; i < m; ++i) {
|
||||
for (int j = 0; j < m; ++j) {
|
||||
for (int k = 0; k < m; ++k) {
|
||||
c(i, j) += a(i, k) * b(k, j);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return *this;
|
||||
}
|
||||
|
||||
public:
|
||||
int m_; // first dimension
|
||||
int n_; // second dimension
|
||||
std::vector<double> data_; // the matrix' entries
|
||||
};
|
||||
|
||||
// Print the matrix as a table
|
||||
inline std::ostream &operator<<(std::ostream &os, const Matrix &a) {
|
||||
const int width = 10;
|
||||
const int precision = 4;
|
||||
|
||||
const auto originalPrecision = os.precision();
|
||||
os << std::setprecision(precision);
|
||||
|
||||
for (int i = 0; i < a.dim1(); ++i) {
|
||||
for (int j = 0; j < a.dim2(); ++j) {
|
||||
os << std::setw(width) << a(i, j) << " ";
|
||||
}
|
||||
if (i != a.dim1() - 1)
|
||||
os << "\n";
|
||||
}
|
||||
|
||||
os << std::setprecision(originalPrecision);
|
||||
return os;
|
||||
}
|
||||
|
||||
// matrix product
|
||||
inline Matrix operator*(const Matrix &a, const Matrix &b) {
|
||||
if (a.dim2() == b.dim1()) {
|
||||
int m = a.dim1();
|
||||
int n = a.dim2();
|
||||
int p = b.dim2();
|
||||
Matrix c(m, p);
|
||||
for (int i = 0; i < m; ++i) {
|
||||
for (int j = 0; j < p; ++j) {
|
||||
for (int k = 0; k < n; ++k) {
|
||||
c(i, j) += a(i, k) * b(k, j);
|
||||
}
|
||||
}
|
||||
}
|
||||
return c;
|
||||
} else {
|
||||
return Matrix(0, 0);
|
||||
}
|
||||
}
|
||||
|
||||
inline bool equalWithinRange(const Matrix &a, const Matrix &b,
|
||||
double eps = 1e-12) {
|
||||
if (a.dim1() != b.dim1() || a.dim2() != b.dim2())
|
||||
return false;
|
||||
|
||||
int m = a.dim1();
|
||||
int n = a.dim2();
|
||||
for (int i = 0; i < m; ++i) {
|
||||
for (int j = 0; j < n; ++j) {
|
||||
if (fabs(a(i, j) - b(i, j)) > eps) {
|
||||
return false;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
// A very simple class for "3D-Matrices" (tensors) with dimension l x m x n
|
||||
class Matrix3D {
|
||||
public:
|
||||
// constructors
|
||||
Matrix3D(int l, int m, int n) : l_(l), m_(m), n_(n), data_(l) {
|
||||
for (int i = 0; i < l_; ++i) {
|
||||
data_[i] = std::vector<std::vector<double>>(m_);
|
||||
for (int j = 0; j < m_; ++j) {
|
||||
data_[i][j] = std::vector<double>(n_, 0);
|
||||
}
|
||||
}
|
||||
}
|
||||
Matrix3D(int n) : Matrix3D(n, n, n) {}
|
||||
Matrix3D(const Matrix3D &other) = default;
|
||||
Matrix3D(Matrix3D &&other) = default;
|
||||
~Matrix3D() = default;
|
||||
|
||||
// assignment operators
|
||||
Matrix3D &operator=(const Matrix3D &other) = default;
|
||||
Matrix3D &operator=(Matrix3D &&other) = default;
|
||||
|
||||
// element access
|
||||
double &operator()(int i, int j, int k) { return data_[i][j][k]; }
|
||||
const double &operator()(int i, int j, int k) const { return data_[i][j][k]; }
|
||||
|
||||
// getter functions for the dimensions
|
||||
int dim1() const { return l_; }
|
||||
int dim2() const { return m_; }
|
||||
int dim3() const { return n_; }
|
||||
|
||||
// comparison operators
|
||||
bool operator==(const Matrix3D &b) { return (data_ == b.data_); }
|
||||
bool operator!=(const Matrix3D &b) { return (data_ != b.data_); }
|
||||
|
||||
// addition
|
||||
Matrix3D &operator+=(const Matrix3D &b) {
|
||||
for (int i = 0; i < l_; ++i) {
|
||||
for (int j = 0; j < m_; ++j) {
|
||||
for (int k = 0; k < n_; ++k) {
|
||||
operator()(i, j, k) += b(i, j, k);
|
||||
}
|
||||
}
|
||||
}
|
||||
return *this;
|
||||
}
|
||||
|
||||
// substraction
|
||||
Matrix3D &operator-=(const Matrix3D &b) {
|
||||
for (int i = 0; i < l_; ++i) {
|
||||
for (int j = 0; j < m_; ++j) {
|
||||
for (int k = 0; k < n_; ++k) {
|
||||
operator()(i, j, k) -= b(i, j, k);
|
||||
}
|
||||
}
|
||||
}
|
||||
return *this;
|
||||
}
|
||||
|
||||
// scalar multiplication
|
||||
Matrix3D &operator*=(double x) {
|
||||
for (int i = 0; i < l_; ++i) {
|
||||
for (int j = 0; j < m_; ++j) {
|
||||
for (int k = 0; k < n_; ++k) {
|
||||
operator()(i, j, k) *= x;
|
||||
}
|
||||
}
|
||||
}
|
||||
return *this;
|
||||
}
|
||||
|
||||
// scalar division
|
||||
Matrix3D &operator/=(double x) {
|
||||
for (int i = 0; i < l_; ++i) {
|
||||
for (int j = 0; j < m_; ++j) {
|
||||
for (int k = 0; k < n_; ++k) {
|
||||
operator()(i, j, k) /= x;
|
||||
}
|
||||
}
|
||||
}
|
||||
return *this;
|
||||
}
|
||||
|
||||
private:
|
||||
int l_; // first dimension
|
||||
int m_; // second dimension
|
||||
int n_; // third dimension
|
||||
std::vector<std::vector<std::vector<double>>> data_; // the tensors' entries
|
||||
};
|
||||
|
||||
#endif // MATRIX_H
|
||||
Reference in New Issue
Block a user