diff --git a/.gitignore b/.gitignore index 36b2056..d738ded 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,3 @@ /build .clangd +exe-* diff --git a/data/matrix_band_klein.mtx b/data/matrix_band_klein.mtx new file mode 100644 index 0000000..027c4c4 --- /dev/null +++ b/data/matrix_band_klein.mtx @@ -0,0 +1,300 @@ +%%MatrixMarket matrix coordinate real general +100 100 298 + 1 1 2 + 2 2 3 + 3 3 3 + 4 4 3 + 5 5 3 + 6 6 3 + 7 7 3 + 8 8 3 + 9 9 3 + 10 10 3 + 11 11 3 + 12 12 3 + 13 13 3 + 14 14 3 + 15 15 3 + 16 16 3 + 17 17 3 + 18 18 3 + 19 19 3 + 20 20 3 + 21 21 3 + 22 22 3 + 23 23 3 + 24 24 3 + 25 25 3 + 26 26 3 + 27 27 3 + 28 28 3 + 29 29 3 + 30 30 3 + 31 31 3 + 32 32 3 + 33 33 3 + 34 34 3 + 35 35 3 + 36 36 3 + 37 37 3 + 38 38 3 + 39 39 3 + 40 40 3 + 41 41 3 + 42 42 3 + 43 43 3 + 44 44 3 + 45 45 3 + 46 46 3 + 47 47 3 + 48 48 3 + 49 49 3 + 50 50 3 + 51 51 3 + 52 52 3 + 53 53 3 + 54 54 3 + 55 55 3 + 56 56 3 + 57 57 3 + 58 58 3 + 59 59 3 + 60 60 3 + 61 61 3 + 62 62 3 + 63 63 3 + 64 64 3 + 65 65 3 + 66 66 3 + 67 67 3 + 68 68 3 + 69 69 3 + 70 70 3 + 71 71 3 + 72 72 3 + 73 73 3 + 74 74 3 + 75 75 3 + 76 76 3 + 77 77 3 + 78 78 3 + 79 79 3 + 80 80 3 + 81 81 3 + 82 82 3 + 83 83 3 + 84 84 3 + 85 85 3 + 86 86 3 + 87 87 3 + 88 88 3 + 89 89 3 + 90 90 3 + 91 91 3 + 92 92 3 + 93 93 3 + 94 94 3 + 95 95 3 + 96 96 3 + 97 97 3 + 98 98 3 + 99 99 3 +100 100 2 + 1 2 -1 + 2 3 -1 + 3 4 -1 + 4 5 -1 + 5 6 -1 + 6 7 -1 + 7 8 -1 + 8 9 -1 + 9 10 -1 + 10 11 -1 + 11 12 -1 + 12 13 -1 + 13 14 -1 + 14 15 -1 + 15 16 -1 + 16 17 -1 + 17 18 -1 + 18 19 -1 + 19 20 -1 + 20 21 -1 + 21 22 -1 + 22 23 -1 + 23 24 -1 + 24 25 -1 + 25 26 -1 + 26 27 -1 + 27 28 -1 + 28 29 -1 + 29 30 -1 + 30 31 -1 + 31 32 -1 + 32 33 -1 + 33 34 -1 + 34 35 -1 + 35 36 -1 + 36 37 -1 + 37 38 -1 + 38 39 -1 + 39 40 -1 + 40 41 -1 + 41 42 -1 + 42 43 -1 + 43 44 -1 + 44 45 -1 + 45 46 -1 + 46 47 -1 + 47 48 -1 + 48 49 -1 + 49 50 -1 + 50 51 -1 + 51 52 -1 + 52 53 -1 + 53 54 -1 + 54 55 -1 + 55 56 -1 + 56 57 -1 + 57 58 -1 + 58 59 -1 + 59 60 -1 + 60 61 -1 + 61 62 -1 + 62 63 -1 + 63 64 -1 + 64 65 -1 + 65 66 -1 + 66 67 -1 + 67 68 -1 + 68 69 -1 + 69 70 -1 + 70 71 -1 + 71 72 -1 + 72 73 -1 + 73 74 -1 + 74 75 -1 + 75 76 -1 + 76 77 -1 + 77 78 -1 + 78 79 -1 + 79 80 -1 + 80 81 -1 + 81 82 -1 + 82 83 -1 + 83 84 -1 + 84 85 -1 + 85 86 -1 + 86 87 -1 + 87 88 -1 + 88 89 -1 + 89 90 -1 + 90 91 -1 + 91 92 -1 + 92 93 -1 + 93 94 -1 + 94 95 -1 + 95 96 -1 + 96 97 -1 + 97 98 -1 + 98 99 -1 + 99 100 -1 + 2 1 -1 + 3 2 -1 + 4 3 -1 + 5 4 -1 + 6 5 -1 + 7 6 -1 + 8 7 -1 + 9 8 -1 + 10 9 -1 + 11 10 -1 + 12 11 -1 + 13 12 -1 + 14 13 -1 + 15 14 -1 + 16 15 -1 + 17 16 -1 + 18 17 -1 + 19 18 -1 + 20 19 -1 + 21 20 -1 + 22 21 -1 + 23 22 -1 + 24 23 -1 + 25 24 -1 + 26 25 -1 + 27 26 -1 + 28 27 -1 + 29 28 -1 + 30 29 -1 + 31 30 -1 + 32 31 -1 + 33 32 -1 + 34 33 -1 + 35 34 -1 + 36 35 -1 + 37 36 -1 + 38 37 -1 + 39 38 -1 + 40 39 -1 + 41 40 -1 + 42 41 -1 + 43 42 -1 + 44 43 -1 + 45 44 -1 + 46 45 -1 + 47 46 -1 + 48 47 -1 + 49 48 -1 + 50 49 -1 + 51 50 -1 + 52 51 -1 + 53 52 -1 + 54 53 -1 + 55 54 -1 + 56 55 -1 + 57 56 -1 + 58 57 -1 + 59 58 -1 + 60 59 -1 + 61 60 -1 + 62 61 -1 + 63 62 -1 + 64 63 -1 + 65 64 -1 + 66 65 -1 + 67 66 -1 + 68 67 -1 + 69 68 -1 + 70 69 -1 + 71 70 -1 + 72 71 -1 + 73 72 -1 + 74 73 -1 + 75 74 -1 + 76 75 -1 + 77 76 -1 + 78 77 -1 + 79 78 -1 + 80 79 -1 + 81 80 -1 + 82 81 -1 + 83 82 -1 + 84 83 -1 + 85 84 -1 + 86 85 -1 + 87 86 -1 + 88 87 -1 + 89 88 -1 + 90 89 -1 + 91 90 -1 + 92 91 -1 + 93 92 -1 + 94 93 -1 + 95 94 -1 + 96 95 -1 + 97 96 -1 + 98 97 -1 + 99 98 -1 +100 99 -1 diff --git a/src/main.c b/src/main.c index abdc431..77a8fee 100644 --- a/src/main.c +++ b/src/main.c @@ -7,67 +7,42 @@ #include #include "comm.h" +#include "matrix.h" #include "progress.h" #include "solver.h" #include "timing.h" -static FILE* initResidualWriter() -{ - FILE* fp; - fp = fopen("residual.dat", "w"); - - if (fp == NULL) { - printf("Error!\n"); - exit(EXIT_FAILURE); - } - - return fp; -} - -static void writeResidual(FILE* fp, double ts, double res) -{ - fprintf(fp, "%f, %f\n", ts, res); -} - int main(int argc, char** argv) { int rank; double timeStart, timeStop; - Parameter p; + Matrix m; Solver s; - commInit(s.comm, argc, argv); - initParameter(&p); + // commInit(s.comm, argc, argv); - FILE* fp; - if (commIsMaster(s.comm)) fp = initResidualWriter(); + // FILE* fp; + // if (commIsMaster(s.comm)) fp = initResidualWriter(); if (argc != 2) { printf("Usage: %s \n", argv[0]); exit(EXIT_SUCCESS); } - readParameter(&p, argv[1]); - commPartition(s.comm, p.jmax, p.imax); - if (commIsMaster(s.comm)) { - printParameter(&p); - } + matrixRead(&m, argv[1]); + matrixDump(&m); + exit(EXIT_SUCCESS); - initSolver(&s, &p); - // #ifndef VERBOSE - // initProgress(s.te); - // #endif - - double te = p.te; + double te = 1000; double t = 0.0; double res = 0.0; timeStart = getTimeStamp(); while (t <= te) { - if (commIsMaster(s.comm)) writeResidual(fp, t, res); + // if (commIsMaster(s.comm)) writeResidual(fp, t, res); - t += p.dt; + t += 0.5; #ifdef VERBOSE if (commIsMaster(s.comm)) { @@ -78,13 +53,10 @@ int main(int argc, char** argv) #endif } timeStop = getTimeStamp(); -#ifndef VERBOSE - stopProgress(); -#endif if (commIsMaster(s.comm)) { printf("Solution took %.2fs\n", timeStop - timeStart); } - if (commIsMaster(s.comm)) fclose(fp); + // if (commIsMaster(s.comm)) fclose(fp); commFinalize(s.comm); return EXIT_SUCCESS; } diff --git a/src/matrix.c b/src/matrix.c index 1e19394..2ae408a 100644 --- a/src/matrix.c +++ b/src/matrix.c @@ -31,7 +31,6 @@ static inline int compareRow(const void* a, const void* b) void matrixRead(Matrix* m, char* filename) { - int ret_code; MM_typecode matcode; FILE* f; int M, N, nz; @@ -81,6 +80,8 @@ void matrixRead(Matrix* m, char* filename) exit(EXIT_FAILURE); } + printf("Read matrix %s with %d non zeroes and %d rows", filename, nz, M); + m->nr = M; m->nnz = nz; Entry* mm; @@ -121,10 +122,59 @@ void matrixRead(Matrix* m, char* filename) } fclose(f); + size_t mms = cursor; // sort by column - qsort(mm, cursor, sizeof(Entry), compareColumn); - + qsort(mm, mms, sizeof(Entry), compareColumn); // sort by row - qsort(mm, cursor, sizeof(Entry), compareRow); + qsort(mm, mms, sizeof(Entry), compareRow); + + m->rowPtr = (CG_UINT*)allocate(64, (m->nr + 1) * sizeof(CG_UINT)); + m->colInd = (CG_UINT*)allocate(64, m->nnz * sizeof(CG_UINT)); + m->val = (CG_FLOAT*)allocate(64, m->nnz * sizeof(CG_FLOAT)); + + int* valsPerRow = (int*)allocate(64, m->nr * sizeof(int)); + + for (int i = 0; i < m->nr; i++) { + valsPerRow[i] = 0; + } + + for (int i = 0; i < mms; i++) { + valsPerRow[mm[i].row]++; + } + + int* offsets = (int*)allocate(64, (m->nr + 1) * sizeof(int)); + m->rowPtr[0] = 0; + + // convert to CRS format + for (int rowID = 0; rowID < m->nr; rowID++) { + + m->rowPtr[rowID + 1] = m->rowPtr[rowID] + valsPerRow[rowID]; + + // loop over all elements in Row + for (int id = m->rowPtr[rowID]; id < m->rowPtr[rowID + 1]; id++) { + m->val[id] = (CG_UINT)mm[id].val; + m->colInd[id] = (CG_FLOAT)mm[id].col; + } + } +} + +void matrixDump(Matrix* m) +{ + CG_UINT numRows = m->nr; + CG_UINT* rowPtr = m->rowPtr; + CG_UINT* colInd = m->colInd; + CG_FLOAT* val = m->val; + + printf("Matrix: %lld non zeroes, number of rows %lld\n", numRows, m->nnz); + + for (int rowID = 0; rowID < numRows; rowID++) { + printf("Row [%d]: ", rowID); + + for (size_t rowEntry = rowPtr[rowID]; rowEntry < rowPtr[rowID + 1]; rowEntry++) { + printf("[%lld]:%.2f ", colInd[rowEntry], val[rowEntry]); + } + + printf("\n"); + } } diff --git a/src/matrix.h b/src/matrix.h index 0a912dd..7f3a941 100644 --- a/src/matrix.h +++ b/src/matrix.h @@ -20,5 +20,6 @@ typedef struct { } Matrix; extern void matrixRead(Matrix* m, char* filename); +extern void matrixDump(Matrix* m); #endif // __MATRIX_H_ diff --git a/src/solver.c b/src/solver.c new file mode 100644 index 0000000..556ad68 --- /dev/null +++ b/src/solver.c @@ -0,0 +1,28 @@ +/* Copyright (C) NHR@FAU, University Erlangen-Nuremberg. + * All rights reserved. + * Use of this source code is governed by a MIT-style + * license that can be found in the LICENSE file.*/ +#include + +#include "matrix.h" +#include "solver.h" +#include "util.h" + +void spMVM(Matrix* m, double* x, double* y) +{ + CG_UINT numRows = m->nr; + CG_UINT* rowPtr = m->rowPtr; + CG_UINT* colInd = m->colInd; + CG_FLOAT* val = m->val; + + for (size_t rowID = 0; rowID < numRows; rowID++) { + CG_FLOAT tmp = y[rowID]; + + // loop over all elements in row + for (size_t rowEntry = rowPtr[rowID]; rowEntry < rowPtr[rowID + 1]; rowEntry++) { + tmp += val[rowEntry] * x[colInd[rowEntry]]; + } + + y[rowID] = tmp; + } +} diff --git a/src/solver.h b/src/solver.h index 4a0fc20..be35da6 100644 --- a/src/solver.h +++ b/src/solver.h @@ -1,26 +1,24 @@ -/* - * Copyright (C) NHR@FAU, University Erlangen-Nuremberg. +/* Copyright (C) NHR@FAU, University Erlangen-Nuremberg. * All rights reserved. This file is part of nusif-solver. * Use of this source code is governed by a MIT style - * license that can be found in the LICENSE file. - */ + * license that can be found in the LICENSE file. */ #ifndef __SOLVER_H_ #define __SOLVER_H_ #include "comm.h" -#include "mpi.h" +#include "matrix.h" #include "parameter.h" +#include "util.h" typedef struct { /* geometry and grid information */ /* parameters */ double eps, omega; int itermax; - int levels, presmooth, postsmooth; - double **r, **e; + /* communication */ Comm* comm; } Solver; void initSolver(Solver*, Parameter*); -double solve(Solver*, double*, double*); -#endif +void spMVM(Matrix* m, CG_FLOAT* x, CG_FLOAT* y); +#endif // __SOLVER_H_