From d355c4bbdb1856ffdaeeecc4619035c93121d828 Mon Sep 17 00:00:00 2001 From: Jan Eitzinger Date: Fri, 3 Jan 2025 14:23:28 +0100 Subject: [PATCH] Initial checkin --- .clang-format | 176 ++++++++++++++++++ .clang-tidy | 14 ++ .gitignore | 2 + Makefile | 78 ++++++++ config.mk | 9 + mk/include_CLANG.mk | 15 ++ mk/include_GCC.mk | 14 ++ mk/include_ICC.mk | 14 ++ src/allocate.c | 38 ++++ src/allocate.h | 11 ++ src/comm.c | 129 +++++++++++++ src/comm.h | 55 ++++++ src/main.c | 90 ++++++++++ src/matrix.c | 130 ++++++++++++++ src/matrix.h | 24 +++ src/mmio.c | 429 ++++++++++++++++++++++++++++++++++++++++++++ src/mmio.h | 129 +++++++++++++ src/parameter.c | 114 ++++++++++++ src/parameter.h | 27 +++ src/progress.c | 51 ++++++ src/progress.h | 13 ++ src/solver.h | 26 +++ src/timing.c | 22 +++ src/timing.h | 13 ++ src/util.h | 36 ++++ 25 files changed, 1659 insertions(+) create mode 100644 .clang-format create mode 100644 .clang-tidy create mode 100644 .gitignore create mode 100644 Makefile create mode 100644 config.mk create mode 100644 mk/include_CLANG.mk create mode 100644 mk/include_GCC.mk create mode 100644 mk/include_ICC.mk create mode 100644 src/allocate.c create mode 100644 src/allocate.h create mode 100644 src/comm.c create mode 100644 src/comm.h create mode 100644 src/main.c create mode 100644 src/matrix.c create mode 100644 src/matrix.h create mode 100644 src/mmio.c create mode 100644 src/mmio.h create mode 100644 src/parameter.c create mode 100644 src/parameter.h create mode 100644 src/progress.c create mode 100644 src/progress.h create mode 100644 src/solver.h create mode 100644 src/timing.c create mode 100644 src/timing.h create mode 100644 src/util.h diff --git a/.clang-format b/.clang-format new file mode 100644 index 0000000..38404b9 --- /dev/null +++ b/.clang-format @@ -0,0 +1,176 @@ +--- +Language: Cpp +# BasedOnStyle: WebKit +AccessModifierOffset: -4 +AlignAfterOpenBracket: DontAlign +AlignArrayOfStructures: None +AlignConsecutiveAssignments: Consecutive +AlignConsecutiveBitFields: None +AlignConsecutiveDeclarations: None +AlignConsecutiveMacros: Consecutive +AlignEscapedNewlines: Right +AlignOperands: Align +AlignTrailingComments: true +AllowAllArgumentsOnNextLine: false +AllowAllParametersOfDeclarationOnNextLine: true +AllowShortEnumsOnASingleLine: true +AllowShortBlocksOnASingleLine: Never +AllowShortCaseLabelsOnASingleLine: false +AllowShortFunctionsOnASingleLine: All +AllowShortLambdasOnASingleLine: All +AllowShortIfStatementsOnASingleLine: OnlyFirstIf +AllowShortLoopsOnASingleLine: false +AlwaysBreakAfterDefinitionReturnType: None +AlwaysBreakAfterReturnType: None +AlwaysBreakBeforeMultilineStrings: false +AlwaysBreakTemplateDeclarations: MultiLine +AttributeMacros: + - __capability +BinPackArguments: false +BinPackParameters: false +BraceWrapping: + AfterCaseLabel: false + AfterClass: false + AfterControlStatement: Never + AfterEnum: false + AfterFunction: true + AfterNamespace: false + AfterObjCDeclaration: false + AfterStruct: false + AfterUnion: false + AfterExternBlock: false + BeforeCatch: false + BeforeElse: false + BeforeLambdaBody: false + BeforeWhile: false + IndentBraces: false + SplitEmptyFunction: true + SplitEmptyRecord: true + SplitEmptyNamespace: true +BreakBeforeBinaryOperators: None +BreakBeforeBraces: WebKit +BreakBeforeInheritanceComma: false +BreakInheritanceList: BeforeColon +BreakBeforeTernaryOperators: true +BreakConstructorInitializersBeforeComma: false +BreakConstructorInitializers: BeforeComma +BreakAfterJavaFieldAnnotations: false +BreakStringLiterals: true +ColumnLimit: 90 +CommentPragmas: '^ IWYU pragma:' +CompactNamespaces: false +ConstructorInitializerIndentWidth: 4 +ContinuationIndentWidth: 4 +Cpp11BracedListStyle: false +DeriveLineEnding: true +DerivePointerAlignment: false +DisableFormat: false +EmptyLineAfterAccessModifier: Never +EmptyLineBeforeAccessModifier: LogicalBlock +ExperimentalAutoDetectBinPacking: false +BasedOnStyle: '' +ConstructorInitializerAllOnOneLineOrOnePerLine: false +AllowAllConstructorInitializersOnNextLine: true +FixNamespaceComments: false +ForEachMacros: + - foreach + - Q_FOREACH + - BOOST_FOREACH +IfMacros: + - KJ_IF_MAYBE +IncludeBlocks: Preserve +IncludeCategories: + - Regex: '^"(llvm|llvm-c|clang|clang-c)/' + Priority: 2 + SortPriority: 0 + CaseSensitive: false + - Regex: '^(<|"(gtest|gmock|isl|json)/)' + Priority: 3 + SortPriority: 0 + CaseSensitive: false + - Regex: '.*' + Priority: 1 + SortPriority: 0 + CaseSensitive: false +IncludeIsMainRegex: '(Test)?$' +IncludeIsMainSourceRegex: '' +IndentAccessModifiers: false +IndentCaseLabels: false +IndentCaseBlocks: false +IndentGotoLabels: true +IndentPPDirectives: None +IndentExternBlock: AfterExternBlock +IndentWidth: 4 +IndentWrappedFunctionNames: false +InsertTrailingCommas: None +JavaScriptQuotes: Leave +JavaScriptWrapImports: true +KeepEmptyLinesAtTheStartOfBlocks: true +LambdaBodyIndentation: Signature +MacroBlockBegin: '' +MacroBlockEnd: '' +MaxEmptyLinesToKeep: 1 +NamespaceIndentation: Inner +ObjCBinPackProtocolList: Auto +ObjCBlockIndentWidth: 4 +ObjCBreakBeforeNestedBlockParam: true +ObjCSpaceAfterProperty: true +ObjCSpaceBeforeProtocolList: true +PenaltyBreakAssignment: 200 +PenaltyBreakBeforeFirstCallParameter: 19 +PenaltyBreakComment: 300 +PenaltyBreakFirstLessLess: 120 +PenaltyBreakString: 1000 +PenaltyBreakTemplateDeclaration: 10 +PenaltyExcessCharacter: 1000000 +PenaltyReturnTypeOnItsOwnLine: 60 +PenaltyIndentedWhitespace: 0 +PointerAlignment: Left +PPIndentWidth: -1 +ReferenceAlignment: Pointer +ReflowComments: true +ShortNamespaceLines: 1 +SortIncludes: CaseSensitive +SortJavaStaticImport: Before +SortUsingDeclarations: true +SpaceAfterCStyleCast: false +SpaceAfterLogicalNot: false +SpaceAfterTemplateKeyword: true +SpaceBeforeAssignmentOperators: true +SpaceBeforeCaseColon: false +SpaceBeforeCpp11BracedList: true +SpaceBeforeCtorInitializerColon: true +SpaceBeforeInheritanceColon: true +SpaceBeforeParens: ControlStatements +SpaceAroundPointerQualifiers: Default +SpaceBeforeRangeBasedForLoopColon: true +SpaceInEmptyBlock: false +SpaceInEmptyParentheses: false +SpacesBeforeTrailingComments: 1 +SpacesInAngles: Never +SpacesInConditionalStatement: false +SpacesInContainerLiterals: true +SpacesInCStyleCastParentheses: false +SpacesInLineCommentPrefix: + Minimum: 1 + Maximum: -1 +SpacesInParentheses: false +SpacesInSquareBrackets: false +SpaceBeforeSquareBrackets: false +BitFieldColonSpacing: Both +Standard: Latest +StatementAttributeLikeMacros: + - Q_EMIT +StatementMacros: + - Q_UNUSED + - QT_REQUIRE_VERSION +TabWidth: 8 +UseCRLF: false +UseTab: Never +WhitespaceSensitiveMacros: + - STRINGIZE + - PP_STRINGIZE + - BOOST_PP_STRINGIZE + - NS_SWIFT_NAME + - CF_SWIFT_NAME +... diff --git a/.clang-tidy b/.clang-tidy new file mode 100644 index 0000000..acdff1e --- /dev/null +++ b/.clang-tidy @@ -0,0 +1,14 @@ +--- +Checks: 'clang-diagnostic-*,clang-analyzer-*,clang-bugprone-*,readability-identifier-naming' +WarningsAsErrors: true +HeaderFilterRegex: '.*' +AnalyzeTemporaryDtors: false +CheckOptions: + - key: readability-identifier-naming.StructCase + value: 'CamelCase' + - key: readability-identifier-naming.FunctionCase + value: 'camelBack' + - key: readability-identifier-naming.VariableCase + value: 'camelBack' + - key: readability-identifier-naming.GlobalConstantCase + value: 'UPPER_CASE' diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..36b2056 --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +/build +.clangd diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..8b8c60a --- /dev/null +++ b/Makefile @@ -0,0 +1,78 @@ +#======================================================================================= +# 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. +#======================================================================================= + +#CONFIGURE BUILD SYSTEM +TARGET = exe-$(TOOLCHAIN) +BUILD_DIR = ./build/$(TOOLCHAIN) +SRC_DIR = ./src +MAKE_DIR = ./mk +Q ?= @ + +#DO NOT EDIT BELOW +include config.mk +include $(MAKE_DIR)/include_$(TOOLCHAIN).mk +INCLUDES += -I$(SRC_DIR)/includes -I$(BUILD_DIR) + +VPATH = $(SRC_DIR) +ASM = $(patsubst $(SRC_DIR)/%.c, $(BUILD_DIR)/%.s,$(wildcard $(SRC_DIR)/*.c)) +OBJ = $(patsubst $(SRC_DIR)/%.c, $(BUILD_DIR)/%.o,$(wildcard $(SRC_DIR)/*.c)) +SRC = $(wildcard $(SRC_DIR)/*.h $(SRC_DIR)/*.c) +CPPFLAGS := $(CPPFLAGS) $(DEFINES) $(OPTIONS) $(INCLUDES) +c := , +clist = $(subst $(eval) ,$c,$(strip $1)) + +define CLANGD_TEMPLATE +CompileFlags: + Add: [$(call clist,$(INCLUDES)), $(call clist,$(CPPFLAGS))] + Compiler: clang +endef + +${TARGET}: $(BUILD_DIR) .clangd $(OBJ) + $(info ===> LINKING $(TARGET)) + $(Q)${LD} ${LFLAGS} -o $(TARGET) $(OBJ) $(LIBS) + +$(BUILD_DIR)/%.o: %.c $(MAKE_DIR)/include_$(TOOLCHAIN).mk + $(info ===> COMPILE $@) + $(CC) -c $(CPPFLAGS) $(CFLAGS) $< -o $@ + $(Q)$(CC) $(CPPFLAGS) -MT $(@:.d=.o) -MM $< > $(BUILD_DIR)/$*.d + +$(BUILD_DIR)/%.s: %.c + $(info ===> GENERATE ASM $@) + $(CC) -S $(CPPFLAGS) $(CFLAGS) $< -o $@ + +.PHONY: clean distclean info asm format + +clean: + $(info ===> CLEAN) + @rm -rf $(BUILD_DIR) + +distclean: + $(info ===> DIST CLEAN) + @rm -rf build + @rm -f $(TARGET) + @rm -f tags .clangd + +info: + $(info $(CFLAGS)) + $(Q)$(CC) $(VERSION) + +asm: $(BUILD_DIR) $(ASM) + +format: + @for src in $(SRC) ; do \ + echo "Formatting $$src" ; \ + clang-format -i $$src ; \ + done + @echo "Done" + +$(BUILD_DIR): + @mkdir -p $(BUILD_DIR) + +.clangd: + $(file > .clangd,$(CLANGD_TEMPLATE)) + +-include $(OBJ:.o=.d) diff --git a/config.mk b/config.mk new file mode 100644 index 0000000..7d53436 --- /dev/null +++ b/config.mk @@ -0,0 +1,9 @@ +# Supported: GCC, CLANG, ICC +TOOLCHAIN ?= CLANG +# ENABLE_OPENMP ?= false + +#Feature options +OPTIONS += -DARRAY_ALIGNMENT=64 +#OPTIONS += -DVERBOSE_AFFINITY +#OPTIONS += -DVERBOSE_DATASIZE +#OPTIONS += -DVERBOSE_TIMER diff --git a/mk/include_CLANG.mk b/mk/include_CLANG.mk new file mode 100644 index 0000000..429894d --- /dev/null +++ b/mk/include_CLANG.mk @@ -0,0 +1,15 @@ +CC = mpicc +LD = $(CC) + +ifeq ($(ENABLE_OPENMP),true) +OPENMP = -fopenmp +#OPENMP = -Xpreprocessor -fopenmp #required on Macos with homebrew libomp +LIBS = # -lomp +endif + +VERSION = --version +CFLAGS = -Ofast -std=c99 $(OPENMP) +#CFLAGS = -Ofast -fnt-store=aggressive -std=c99 $(OPENMP) #AMD CLANG +LFLAGS = $(OPENMP) +DEFINES = -D_GNU_SOURCE +INCLUDES = -I~/.local/include diff --git a/mk/include_GCC.mk b/mk/include_GCC.mk new file mode 100644 index 0000000..427e798 --- /dev/null +++ b/mk/include_GCC.mk @@ -0,0 +1,14 @@ +CC = gcc +GCC = gcc +LINKER = $(CC) + +ifeq ($(ENABLE_OPENMP),true) +OPENMP = -fopenmp +endif + +VERSION = --version +CFLAGS = -Ofast -ffreestanding -std=c99 $(OPENMP) +LFLAGS = $(OPENMP) +DEFINES = -D_GNU_SOURCE +INCLUDES = +LIBS = diff --git a/mk/include_ICC.mk b/mk/include_ICC.mk new file mode 100644 index 0000000..e172b15 --- /dev/null +++ b/mk/include_ICC.mk @@ -0,0 +1,14 @@ +CC = icc +GCC = gcc +LINKER = $(CC) + +ifeq ($(ENABLE_OPENMP),true) +OPENMP = -qopenmp +endif + +VERSION = --version +CFLAGS = -fast -xHost -qopt-streaming-stores=always -std=c99 -ffreestanding $(OPENMP) +LFLAGS = $(OPENMP) +DEFINES = -D_GNU_SOURCE +INCLUDES = +LIBS = diff --git a/src/allocate.c b/src/allocate.c new file mode 100644 index 0000000..cf2efd6 --- /dev/null +++ b/src/allocate.c @@ -0,0 +1,38 @@ +/* + * 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 +#include +#include + +#include "allocate.h" + +void* allocate(size_t alignment, size_t bytesize) +{ + int errorCode; + void* ptr; + + errorCode = posix_memalign(&ptr, alignment, bytesize); + + if (errorCode) { + if (errorCode == EINVAL) { + fprintf(stderr, "Error: Alignment parameter is not a power of two\n"); + exit(EXIT_FAILURE); + } + if (errorCode == ENOMEM) { + fprintf(stderr, "Error: Insufficient memory to fulfill the request\n"); + exit(EXIT_FAILURE); + } + } + + if (ptr == NULL) { + fprintf(stderr, "Error: posix_memalign failed!\n"); + exit(EXIT_FAILURE); + } + + return ptr; +} diff --git a/src/allocate.h b/src/allocate.h new file mode 100644 index 0000000..48d8f68 --- /dev/null +++ b/src/allocate.h @@ -0,0 +1,11 @@ +/* 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. */ +#ifndef __ALLOCATE_H_ +#define __ALLOCATE_H_ +#include + +extern void* allocate(size_t alignment, size_t bytesize); + +#endif diff --git a/src/comm.c b/src/comm.c new file mode 100644 index 0000000..fab8f98 --- /dev/null +++ b/src/comm.c @@ -0,0 +1,129 @@ +/* + * 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. + */ +#include +#include + +#include "comm.h" + +// subroutines local to this module +int sizeOfRank(int rank, int size, int N) +{ + return N / size + ((N % size > rank) ? 1 : 0); +} + +void commReduction(double* v, int op) +{ +#ifdef _MPI + if (op == MAX) { + MPI_Allreduce(MPI_IN_PLACE, v, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); + } else if (op == SUM) { + MPI_Allreduce(MPI_IN_PLACE, v, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + } +#endif +} + +void commPrintConfig(Comm* c) +{ +#ifdef _MPI + fflush(stdout); + MPI_Barrier(MPI_COMM_WORLD); + if (commIsMaster(c)) { + printf("Communication setup:\n"); + } + + for (int i = 0; i < c->size; i++) { + if (i == c->rank) { + printf("\tRank %d of %d\n", c->rank, c->size); + printf("\tNeighbours (bottom, top, left, right): %d %d, %d, %d\n", + c->neighbours[BOTTOM], + c->neighbours[TOP], + c->neighbours[LEFT], + c->neighbours[RIGHT]); + printf("\tIs boundary:\n"); + printf("\t\tLEFT: %d\n", commIsBoundary(c, LEFT)); + printf("\t\tRIGHT: %d\n", commIsBoundary(c, RIGHT)); + printf("\t\tBOTTOM: %d\n", commIsBoundary(c, BOTTOM)); + printf("\t\tTOP: %d\n", commIsBoundary(c, TOP)); + printf("\tCoordinates (i,j) %d %d\n", c->coords[IDIM], c->coords[JDIM]); + printf("\tDims (i,j) %d %d\n", c->dims[IDIM], c->dims[JDIM]); + printf("\tLocal domain size (i,j) %dx%d\n", c->imaxLocal, c->jmaxLocal); + fflush(stdout); + } + MPI_Barrier(MPI_COMM_WORLD); + } +#endif +} + +void commInit(Comm* c, int argc, char** argv) +{ +#ifdef _MPI + MPI_Init(&argc, &argv); + MPI_Comm_rank(MPI_COMM_WORLD, &(c->rank)); + MPI_Comm_size(MPI_COMM_WORLD, &(c->size)); +#else + c->rank = 0; + c->size = 1; +#endif +} + +void commTestInit(Comm* c, double* p, double* f, double* g) +{ + int imax = c->imaxLocal; + int jmax = c->jmaxLocal; + int rank = c->rank; + + for (int j = 0; j < jmax + 2; j++) { + for (int i = 0; i < imax + 2; i++) { + p[j * (imax + 2) + i] = rank; + f[j * (imax + 2) + i] = rank; + g[j * (imax + 2) + i] = rank; + } + } +} + +static void testWriteFile(char* filename, double* grid, int imax, int jmax) +{ + FILE* fp = fopen(filename, "w"); + + if (fp == NULL) { + printf("Error!\n"); + exit(EXIT_FAILURE); + } + + for (int j = 0; j < jmax + 2; j++) { + for (int i = 0; i < imax + 2; i++) { + fprintf(fp, "%.2f ", grid[j * (imax + 2) + i]); + } + fprintf(fp, "\n"); + } + + fclose(fp); +} + +void commTestWrite(Comm* c, double* p, double* f, double* g) +{ + int imax = c->imaxLocal; + int jmax = c->jmaxLocal; + int rank = c->rank; + + char filename[30]; + snprintf(filename, 30, "ptest-%d.dat", rank); + testWriteFile(filename, p, imax, jmax); + + snprintf(filename, 30, "ftest-%d.dat", rank); + testWriteFile(filename, f, imax, jmax); + + snprintf(filename, 30, "gtest-%d.dat", rank); + testWriteFile(filename, g, imax, jmax); +} + +void commFinalize(Comm* c) +{ +#ifdef _MPI + MPI_Finalize(); +#endif +} diff --git a/src/comm.h b/src/comm.h new file mode 100644 index 0000000..c1bf93f --- /dev/null +++ b/src/comm.h @@ -0,0 +1,55 @@ +/* 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.*/ +#ifndef __COMM_H_ +#define __COMM_H_ +#if defined(_MPI) +#include +#endif + +enum direction { LEFT = 0, RIGHT, BOTTOM, TOP, NDIRS }; +enum dimension { IDIM = 0, JDIM, NDIMS }; +enum cdimension { CJDIM = 0, CIDIM }; +enum layer { HALO = 0, BULK }; +enum op { MAX = 0, SUM }; + +typedef struct { + int rank; + int size; +#if defined(_MPI) + MPI_Comm comm; + MPI_Datatype bufferTypes[NDIRS]; + MPI_Aint sdispls[NDIRS]; + MPI_Aint rdispls[NDIRS]; +#endif + int neighbours[NDIRS]; + int coords[NDIMS], dims[NDIMS]; + int imaxLocal, jmaxLocal; +} Comm; + +extern int sizeOfRank(int rank, int size, int N); +extern void commInit(Comm* c, int argc, char** argv); +extern void commTestInit(Comm* c, double* p, double* f, double* g); +extern void commTestWrite(Comm* c, double* p, double* f, double* g); +extern void commFinalize(Comm* c); +extern void commPartition(Comm* c, int jmax, int imax); +extern void commPrintConfig(Comm*); +extern void commExchange(Comm*, double*); +extern void commShift(Comm* c, double* f, double* g); +extern void commReduction(double* v, int op); +extern int commIsBoundary(Comm* c, int direction); +extern void commUpdateDatatypes(Comm*, Comm*, int, int); +extern void commFreeCommunicator(Comm*); +extern void commCollectResult(Comm* c, + double* ug, + double* vg, + double* pg, + double* u, + double* v, + double* p, + int jmax, + int imax); + +static inline int commIsMaster(Comm* c) { return c->rank == 0; } +#endif // __COMM_H_ diff --git a/src/main.c b/src/main.c new file mode 100644 index 0000000..abdc431 --- /dev/null +++ b/src/main.c @@ -0,0 +1,90 @@ +/* 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 +#include + +#include "comm.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; + Solver s; + + commInit(s.comm, argc, argv); + initParameter(&p); + + 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); + } + + initSolver(&s, &p); + // #ifndef VERBOSE + // initProgress(s.te); + // #endif + + double te = p.te; + double t = 0.0; + double res = 0.0; + + timeStart = getTimeStamp(); + while (t <= te) { + + if (commIsMaster(s.comm)) writeResidual(fp, t, res); + + t += p.dt; + +#ifdef VERBOSE + if (commIsMaster(s.comm)) { + printf("TIME %f , TIMESTEP %f\n", t, d.dt); + } +#else + printProgress(t); +#endif + } + timeStop = getTimeStamp(); +#ifndef VERBOSE + stopProgress(); +#endif + if (commIsMaster(s.comm)) { + printf("Solution took %.2fs\n", timeStop - timeStart); + } + if (commIsMaster(s.comm)) fclose(fp); + commFinalize(s.comm); + return EXIT_SUCCESS; +} diff --git a/src/matrix.c b/src/matrix.c new file mode 100644 index 0000000..1e19394 --- /dev/null +++ b/src/matrix.c @@ -0,0 +1,130 @@ +/* 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 +#include +#include +#include + +#include "allocate.h" +#include "matrix.h" +#include "mmio.h" +#include "util.h" + +static inline int compareColumn(const void* a, const void* b) +{ + const Entry* a_ = (const Entry*)a; + const Entry* b_ = (const Entry*)b; + + return (a_->col > b_->col) - (a_->col < b_->col); +} + +static inline int compareRow(const void* a, const void* b) +{ + const Entry* a_ = (const Entry*)a; + const Entry* b_ = (const Entry*)b; + + return (a_->row > b_->row) - (a_->row < b_->row); +} + +void matrixRead(Matrix* m, char* filename) +{ + int ret_code; + MM_typecode matcode; + FILE* f; + int M, N, nz; + + if ((f = fopen(filename, "r")) == NULL) { + printf("Unable to open file"); + } + + if (mm_read_banner(f, &matcode) != 0) { + printf("Could not process Matrix Market banner.\n"); + exit(EXIT_FAILURE); + } + + if (!((mm_is_real(matcode) || mm_is_pattern(matcode) || mm_is_integer(matcode)) && + mm_is_matrix(matcode) && mm_is_sparse(matcode))) { + fprintf(stderr, "Sorry, this application does not support "); + fprintf(stderr, "Market Market type: [%s]\n", mm_typecode_to_str(matcode)); + exit(EXIT_FAILURE); + } + + bool compatible_flag = (mm_is_sparse(matcode) && + (mm_is_real(matcode) || mm_is_pattern(matcode) || + mm_is_integer(matcode))) && + (mm_is_symmetric(matcode) || mm_is_general(matcode)); + bool sym_flag = mm_is_symmetric(matcode); + bool pattern_flag = mm_is_pattern(matcode); + bool complex_flag = mm_is_complex(matcode); + + if (!compatible_flag) { + printf("The matrix market file provided is not supported.\n Reason :\n"); + if (!mm_is_sparse(matcode)) { + printf(" * matrix has to be sparse\n"); + } + + if (!mm_is_real(matcode) && !(mm_is_pattern(matcode))) { + printf(" * matrix has to be real or pattern\n"); + } + + if (!mm_is_symmetric(matcode) && !mm_is_general(matcode)) { + printf(" * matrix has to be either general or symmetric\n"); + } + + exit(EXIT_FAILURE); + } + + if (mm_read_mtx_crd_size(f, &M, &N, &nz) != 0) { + exit(EXIT_FAILURE); + } + + m->nr = M; + m->nnz = nz; + Entry* mm; + + if (sym_flag) { + mm = (Entry*)allocate(64, nz * 2 * sizeof(Entry)); + } else { + mm = (Entry*)allocate(64, nz * sizeof(Entry)); + } + + size_t cursor = 0; + int row, col; + double v; + + for (size_t i = 0; i < nz; i++) { + + if (pattern_flag) { + fscanf(f, "%d %d\n", &row, &col); + v = 1.; + } else if (complex_flag) { + fscanf(f, "%d %d %lg %*g\n", &row, &col, &v); + } else { + fscanf(f, "%d %d %lg\n", &row, &col, &v); + } + + row--; /* adjust from 1-based to 0-based */ + col--; + + mm[cursor].row = row; + mm[cursor].col = col; + mm[cursor++].val = v; + + if (sym_flag && (row != col)) { + mm[cursor].row = col; + mm[cursor].col = row; + mm[cursor++].val = v; + } + } + + fclose(f); + + // sort by column + qsort(mm, cursor, sizeof(Entry), compareColumn); + + // sort by row + qsort(mm, cursor, sizeof(Entry), compareRow); +} diff --git a/src/matrix.h b/src/matrix.h new file mode 100644 index 0000000..0a912dd --- /dev/null +++ b/src/matrix.h @@ -0,0 +1,24 @@ +/* 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.*/ +#ifndef __MATRIX_H_ +#define __MATRIX_H_ + +#include "util.h" + +typedef struct { + int row; + int col; + double val; +} Entry; + +typedef struct { + CG_UINT nr, nnz; // number of rows and non zeros + CG_UINT *colInd, *rowPtr; // colum Indices, row Pointer + CG_FLOAT* val; +} Matrix; + +extern void matrixRead(Matrix* m, char* filename); + +#endif // __MATRIX_H_ diff --git a/src/mmio.c b/src/mmio.c new file mode 100644 index 0000000..080181b --- /dev/null +++ b/src/mmio.c @@ -0,0 +1,429 @@ +/* + * Matrix Market I/O library for ANSI C + * + * See http://math.nist.gov/MatrixMarket for details. + * + * + */ +#include +#include +#include + +#include "mmio.h" + +int mm_read_unsymmetric_sparse( + const char* fname, int* M_, int* N_, int* nz_, double** val_, int** I_, int** J_) +{ + FILE* f; + MM_typecode matcode; + int M, N, nz; + int i; + double* val; + int *I, *J; + + if ((f = fopen(fname, "r")) == NULL) return -1; + + if (mm_read_banner(f, &matcode) != 0) { + printf("mm_read_unsymetric: Could not process Matrix Market banner "); + printf(" in file [%s]\n", fname); + return -1; + } + + if (!(mm_is_real(matcode) && mm_is_matrix(matcode) && mm_is_sparse(matcode))) { + fprintf(stderr, "Sorry, this application does not support "); + fprintf(stderr, "Market Market type: [%s]\n", mm_typecode_to_str(matcode)); + return -1; + } + + /* find out size of sparse matrix: M, N, nz .... */ + + if (mm_read_mtx_crd_size(f, &M, &N, &nz) != 0) { + fprintf(stderr, "read_unsymmetric_sparse(): could not parse matrix size.\n"); + return -1; + } + + *M_ = M; + *N_ = N; + *nz_ = nz; + + /* reseve memory for matrices */ + + I = (int*)malloc(nz * sizeof(int)); + J = (int*)malloc(nz * sizeof(int)); + val = (double*)malloc(nz * sizeof(double)); + + *val_ = val; + *I_ = I; + *J_ = J; + + /* NOTE: when reading in doubles, ANSI C requires the use of the "l" */ + /* specifier as in "%lg", "%lf", "%le", otherwise errors will occur */ + /* (ANSI C X3.159-1989, Sec. 4.9.6.2, p. 136 lines 13-15) */ + + for (i = 0; i < nz; i++) { + fscanf(f, "%d %d %lg\n", &I[i], &J[i], &val[i]); + I[i]--; /* adjust from 1-based to 0-based */ + J[i]--; + } + fclose(f); + + return 0; +} + +int mm_is_valid(MM_typecode matcode) +{ + if (!mm_is_matrix(matcode)) return 0; + if (mm_is_dense(matcode) && mm_is_pattern(matcode)) return 0; + if (mm_is_real(matcode) && mm_is_hermitian(matcode)) return 0; + if (mm_is_pattern(matcode) && (mm_is_hermitian(matcode) || mm_is_skew(matcode))) + return 0; + return 1; +} + +int mm_read_banner(FILE* f, MM_typecode* matcode) +{ + char line[MM_MAX_LINE_LENGTH]; + char banner[MM_MAX_TOKEN_LENGTH]; + char mtx[MM_MAX_TOKEN_LENGTH]; + char crd[MM_MAX_TOKEN_LENGTH]; + char data_type[MM_MAX_TOKEN_LENGTH]; + char storage_scheme[MM_MAX_TOKEN_LENGTH]; + char* p; + + mm_clear_typecode(matcode); + + if (fgets(line, MM_MAX_LINE_LENGTH, f) == NULL) return MM_PREMATURE_EOF; + + if (sscanf(line, "%s %s %s %s %s", banner, mtx, crd, data_type, storage_scheme) != 5) + return MM_PREMATURE_EOF; + + for (p = mtx; *p != '\0'; *p = tolower(*p), p++) + ; /* convert to lower case */ + for (p = crd; *p != '\0'; *p = tolower(*p), p++) + ; + for (p = data_type; *p != '\0'; *p = tolower(*p), p++) + ; + for (p = storage_scheme; *p != '\0'; *p = tolower(*p), p++) + ; + + /* check for banner */ + if (strncmp(banner, MatrixMarketBanner, strlen(MatrixMarketBanner)) != 0) + return MM_NO_HEADER; + + /* first field should be "mtx" */ + if (strcmp(mtx, MM_MTX_STR) != 0) return MM_UNSUPPORTED_TYPE; + mm_set_matrix(matcode); + + /* second field describes whether this is a sparse matrix (in coordinate + storgae) or a dense array */ + + if (strcmp(crd, MM_SPARSE_STR) == 0) mm_set_sparse(matcode); + else if (strcmp(crd, MM_DENSE_STR) == 0) + mm_set_dense(matcode); + else + return MM_UNSUPPORTED_TYPE; + + /* third field */ + + if (strcmp(data_type, MM_REAL_STR) == 0) mm_set_real(matcode); + else if (strcmp(data_type, MM_COMPLEX_STR) == 0) + mm_set_complex(matcode); + else if (strcmp(data_type, MM_PATTERN_STR) == 0) + mm_set_pattern(matcode); + else if (strcmp(data_type, MM_INT_STR) == 0) + mm_set_integer(matcode); + else + return MM_UNSUPPORTED_TYPE; + + /* fourth field */ + + if (strcmp(storage_scheme, MM_GENERAL_STR) == 0) mm_set_general(matcode); + else if (strcmp(storage_scheme, MM_SYMM_STR) == 0) + mm_set_symmetric(matcode); + else if (strcmp(storage_scheme, MM_HERM_STR) == 0) + mm_set_hermitian(matcode); + else if (strcmp(storage_scheme, MM_SKEW_STR) == 0) + mm_set_skew(matcode); + else + return MM_UNSUPPORTED_TYPE; + + return 0; +} + +int mm_write_mtx_crd_size(FILE* f, int M, int N, int nz) +{ + if (fprintf(f, "%d %d %d\n", M, N, nz) != 3) return MM_COULD_NOT_WRITE_FILE; + else + return 0; +} + +int mm_read_mtx_crd_size(FILE* f, int* M, int* N, int* nz) +{ + char line[MM_MAX_LINE_LENGTH]; + int num_items_read; + + /* set return null parameter values, in case we exit with errors */ + *M = *N = *nz = 0; + + /* now continue scanning until you reach the end-of-comments */ + do { + if (fgets(line, MM_MAX_LINE_LENGTH, f) == NULL) return MM_PREMATURE_EOF; + } while (line[0] == '%'); + + /* line[] is either blank or has M,N, nz */ + if (sscanf(line, "%d %d %d", M, N, nz) == 3) return 0; + + else + do { + num_items_read = fscanf(f, "%d %d %d", M, N, nz); + if (num_items_read == EOF) return MM_PREMATURE_EOF; + } while (num_items_read != 3); + + return 0; +} + +int mm_read_mtx_array_size(FILE* f, int* M, int* N) +{ + char line[MM_MAX_LINE_LENGTH]; + int num_items_read; + /* set return null parameter values, in case we exit with errors */ + *M = *N = 0; + + /* now continue scanning until you reach the end-of-comments */ + do { + if (fgets(line, MM_MAX_LINE_LENGTH, f) == NULL) return MM_PREMATURE_EOF; + } while (line[0] == '%'); + + /* line[] is either blank or has M,N, nz */ + if (sscanf(line, "%d %d", M, N) == 2) return 0; + + else /* we have a blank line */ + do { + num_items_read = fscanf(f, "%d %d", M, N); + if (num_items_read == EOF) return MM_PREMATURE_EOF; + } while (num_items_read != 2); + + return 0; +} + +int mm_write_mtx_array_size(FILE* f, int M, int N) +{ + if (fprintf(f, "%d %d\n", M, N) != 2) return MM_COULD_NOT_WRITE_FILE; + else + return 0; +} + +/*-------------------------------------------------------------------------*/ + +/******************************************************************/ +/* use when I[], J[], and val[]J, and val[] are already allocated */ +/******************************************************************/ + +int mm_read_mtx_crd_data( + FILE* f, int M, int N, int nz, int I[], int J[], double val[], MM_typecode matcode) +{ + int i; + if (mm_is_complex(matcode)) { + for (i = 0; i < nz; i++) + if (fscanf(f, "%d %d %lg %lg", &I[i], &J[i], &val[2 * i], &val[2 * i + 1]) != + 4) + return MM_PREMATURE_EOF; + } else if (mm_is_real(matcode)) { + for (i = 0; i < nz; i++) { + if (fscanf(f, "%d %d %lg\n", &I[i], &J[i], &val[i]) != 3) + return MM_PREMATURE_EOF; + } + } + + else if (mm_is_pattern(matcode)) { + for (i = 0; i < nz; i++) + if (fscanf(f, "%d %d", &I[i], &J[i]) != 2) return MM_PREMATURE_EOF; + } else + return MM_UNSUPPORTED_TYPE; + + return 0; +} + +int mm_read_mtx_crd_entry( + FILE* f, int* I, int* J, double* real, double* imag, MM_typecode matcode) +{ + if (mm_is_complex(matcode)) { + if (fscanf(f, "%d %d %lg %lg", I, J, real, imag) != 4) return MM_PREMATURE_EOF; + } else if (mm_is_real(matcode)) { + if (fscanf(f, "%d %d %lg\n", I, J, real) != 3) return MM_PREMATURE_EOF; + + } + + else if (mm_is_pattern(matcode)) { + if (fscanf(f, "%d %d", I, J) != 2) return MM_PREMATURE_EOF; + } else + return MM_UNSUPPORTED_TYPE; + + return 0; +} + +/************************************************************************ + mm_read_mtx_crd() fills M, N, nz, array of values, and return + type code, e.g. 'MCRS' + + if matrix is complex, values[] is of size 2*nz, + (nz pairs of real/imaginary values) +************************************************************************/ + +int mm_read_mtx_crd(char* fname, + int* M, + int* N, + int* nz, + int** I, + int** J, + double** val, + MM_typecode* matcode) +{ + int ret_code; + FILE* f; + + if (strcmp(fname, "stdin") == 0) f = stdin; + else if ((f = fopen(fname, "r")) == NULL) + return MM_COULD_NOT_READ_FILE; + + if ((ret_code = mm_read_banner(f, matcode)) != 0) return ret_code; + + if (!(mm_is_valid(*matcode) && mm_is_sparse(*matcode) && mm_is_matrix(*matcode))) + return MM_UNSUPPORTED_TYPE; + + if ((ret_code = mm_read_mtx_crd_size(f, M, N, nz)) != 0) return ret_code; + + *I = (int*)malloc(*nz * sizeof(int)); + *J = (int*)malloc(*nz * sizeof(int)); + *val = NULL; + + if (mm_is_complex(*matcode)) { + *val = (double*)malloc(*nz * 2 * sizeof(double)); + ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *I, *J, *val, *matcode); + if (ret_code != 0) return ret_code; + } else if (mm_is_real(*matcode)) { + *val = (double*)malloc(*nz * sizeof(double)); + ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *I, *J, *val, *matcode); + if (ret_code != 0) return ret_code; + } + + else if (mm_is_pattern(*matcode)) { + ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *I, *J, *val, *matcode); + if (ret_code != 0) return ret_code; + } + + if (f != stdin) fclose(f); + return 0; +} + +int mm_write_banner(FILE* f, MM_typecode matcode) +{ + char* str = mm_typecode_to_str(matcode); + int ret_code; + + ret_code = fprintf(f, "%s %s\n", MatrixMarketBanner, str); + free(str); + if (ret_code != 2) return MM_COULD_NOT_WRITE_FILE; + else + return 0; +} + +int mm_write_mtx_crd(char fname[], + int M, + int N, + int nz, + int I[], + int J[], + double val[], + MM_typecode matcode) +{ + FILE* f; + int i; + + if (strcmp(fname, "stdout") == 0) f = stdout; + else if ((f = fopen(fname, "w")) == NULL) + return MM_COULD_NOT_WRITE_FILE; + + /* print banner followed by typecode */ + fprintf(f, "%s ", MatrixMarketBanner); + fprintf(f, "%s\n", mm_typecode_to_str(matcode)); + + /* print matrix sizes and nonzeros */ + fprintf(f, "%d %d %d\n", M, N, nz); + + /* print values */ + if (mm_is_pattern(matcode)) + for (i = 0; i < nz; i++) + fprintf(f, "%d %d\n", I[i], J[i]); + else if (mm_is_real(matcode)) + for (i = 0; i < nz; i++) + fprintf(f, "%d %d %20.16g\n", I[i], J[i], val[i]); + else if (mm_is_complex(matcode)) + for (i = 0; i < nz; i++) + fprintf(f, "%d %d %20.16g %20.16g\n", I[i], J[i], val[2 * i], val[2 * i + 1]); + else { + if (f != stdout) fclose(f); + return MM_UNSUPPORTED_TYPE; + } + + if (f != stdout) fclose(f); + + return 0; +} + +/** + * Create a new copy of a string s. mm_strdup() is a common routine, but + * not part of ANSI C, so it is included here. Used by mm_typecode_to_str(). + * + */ +char* mm_strdup(const char* s) +{ + int len = strlen(s); + char* s2 = (char*)malloc((len + 1) * sizeof(char)); + return strcpy(s2, s); +} + +char* mm_typecode_to_str(MM_typecode matcode) +{ + char buffer[MM_MAX_LINE_LENGTH]; + char* types[4]; + char* mm_strdup(const char*); + int error = 0; + + /* check for MTX type */ + if (mm_is_matrix(matcode)) types[0] = MM_MTX_STR; + else + error = 1; + + /* check for CRD or ARR matrix */ + if (mm_is_sparse(matcode)) types[1] = MM_SPARSE_STR; + else if (mm_is_dense(matcode)) + types[1] = MM_DENSE_STR; + else + return NULL; + + /* check for element data type */ + if (mm_is_real(matcode)) types[2] = MM_REAL_STR; + else if (mm_is_complex(matcode)) + types[2] = MM_COMPLEX_STR; + else if (mm_is_pattern(matcode)) + types[2] = MM_PATTERN_STR; + else if (mm_is_integer(matcode)) + types[2] = MM_INT_STR; + else + return NULL; + + /* check for symmetry type */ + if (mm_is_general(matcode)) types[3] = MM_GENERAL_STR; + else if (mm_is_symmetric(matcode)) + types[3] = MM_SYMM_STR; + else if (mm_is_hermitian(matcode)) + types[3] = MM_HERM_STR; + else if (mm_is_skew(matcode)) + types[3] = MM_SKEW_STR; + else + return NULL; + + sprintf(buffer, "%s %s %s %s", types[0], types[1], types[2], types[3]); + return mm_strdup(buffer); +} diff --git a/src/mmio.h b/src/mmio.h new file mode 100644 index 0000000..db660fa --- /dev/null +++ b/src/mmio.h @@ -0,0 +1,129 @@ +/* + * Matrix Market I/O library for ANSI C + * + * See http://math.nist.gov/MatrixMarket for details. + * + * + */ +#ifndef MM_IO_H +#define MM_IO_H +#include + +#define MM_MAX_LINE_LENGTH 1025 +#define MatrixMarketBanner "%%MatrixMarket" +#define MM_MAX_TOKEN_LENGTH 64 + +typedef char MM_typecode[4]; + +char* mm_typecode_to_str(MM_typecode matcode); + +int mm_read_banner(FILE* f, MM_typecode* matcode); +int mm_read_mtx_crd_size(FILE* f, int* M, int* N, int* nz); +int mm_read_mtx_array_size(FILE* f, int* M, int* N); + +int mm_write_banner(FILE* f, MM_typecode matcode); +int mm_write_mtx_crd_size(FILE* f, int M, int N, int nz); +int mm_write_mtx_array_size(FILE* f, int M, int N); + +/********************* MM_typecode query fucntions ***************************/ + +#define mm_is_matrix(typecode) ((typecode)[0] == 'M') + +#define mm_is_sparse(typecode) ((typecode)[1] == 'C') +#define mm_is_coordinate(typecode) ((typecode)[1] == 'C') +#define mm_is_dense(typecode) ((typecode)[1] == 'A') +#define mm_is_array(typecode) ((typecode)[1] == 'A') + +#define mm_is_complex(typecode) ((typecode)[2] == 'C') +#define mm_is_real(typecode) ((typecode)[2] == 'R') +#define mm_is_pattern(typecode) ((typecode)[2] == 'P') +#define mm_is_integer(typecode) ((typecode)[2] == 'I') + +#define mm_is_symmetric(typecode) ((typecode)[3] == 'S') +#define mm_is_general(typecode) ((typecode)[3] == 'G') +#define mm_is_skew(typecode) ((typecode)[3] == 'K') +#define mm_is_hermitian(typecode) ((typecode)[3] == 'H') + +int mm_is_valid(MM_typecode matcode); /* too complex for a macro */ + +/********************* MM_typecode modify fucntions ***************************/ + +#define mm_set_matrix(typecode) ((*typecode)[0] = 'M') +#define mm_set_coordinate(typecode) ((*typecode)[1] = 'C') +#define mm_set_array(typecode) ((*typecode)[1] = 'A') +#define mm_set_dense(typecode) mm_set_array(typecode) +#define mm_set_sparse(typecode) mm_set_coordinate(typecode) + +#define mm_set_complex(typecode) ((*typecode)[2] = 'C') +#define mm_set_real(typecode) ((*typecode)[2] = 'R') +#define mm_set_pattern(typecode) ((*typecode)[2] = 'P') +#define mm_set_integer(typecode) ((*typecode)[2] = 'I') + +#define mm_set_symmetric(typecode) ((*typecode)[3] = 'S') +#define mm_set_general(typecode) ((*typecode)[3] = 'G') +#define mm_set_skew(typecode) ((*typecode)[3] = 'K') +#define mm_set_hermitian(typecode) ((*typecode)[3] = 'H') + +#define mm_clear_typecode(typecode) \ + ((*typecode)[0] = (*typecode)[1] = (*typecode)[2] = ' ', (*typecode)[3] = 'G') + +#define mm_initialize_typecode(typecode) mm_clear_typecode(typecode) + +/********************* Matrix Market error codes ***************************/ + +#define MM_COULD_NOT_READ_FILE 11 +#define MM_PREMATURE_EOF 12 +#define MM_NOT_MTX 13 +#define MM_NO_HEADER 14 +#define MM_UNSUPPORTED_TYPE 15 +#define MM_LINE_TOO_LONG 16 +#define MM_COULD_NOT_WRITE_FILE 17 + +/******************** Matrix Market internal definitions ******************** + + MM_matrix_typecode: 4-character sequence + + ojbect sparse/ data storage dense + type scheme + + string position: [0] [1] [2] [3] + + Matrix typecode: M(atrix) C(oord) R(eal) G(eneral) + A(array) C(omplex) + H(ermitian) P(attern) S(ymmetric) I(nteger) K(kew) + + ***********************************************************************/ + +#define MM_MTX_STR "matrix" +#define MM_ARRAY_STR "array" +#define MM_DENSE_STR "array" +#define MM_COORDINATE_STR "coordinate" +#define MM_SPARSE_STR "coordinate" +#define MM_COMPLEX_STR "complex" +#define MM_REAL_STR "real" +#define MM_INT_STR "integer" +#define MM_GENERAL_STR "general" +#define MM_SYMM_STR "symmetric" +#define MM_HERM_STR "hermitian" +#define MM_SKEW_STR "skew-symmetric" +#define MM_PATTERN_STR "pattern" + +/* high level routines */ + +int mm_write_mtx_crd(char fname[], + int M, + int N, + int nz, + int I[], + int J[], + double val[], + MM_typecode matcode); +int mm_read_mtx_crd_data( + FILE* f, int M, int N, int nz, int I[], int J[], double val[], MM_typecode matcode); +int mm_read_mtx_crd_entry( + FILE* f, int* I, int* J, double* real, double* img, MM_typecode matcode); + +int mm_read_unsymmetric_sparse( + const char* fname, int* M_, int* N_, int* nz_, double** val_, int** I_, int** J_); + +#endif diff --git a/src/parameter.c b/src/parameter.c new file mode 100644 index 0000000..e5e7029 --- /dev/null +++ b/src/parameter.c @@ -0,0 +1,114 @@ +/* + * 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. + */ +#include +#include +#include + +#include "parameter.h" +#include "util.h" +#define MAXLINE 4096 + +void initParameter(Parameter* param) +{ + param->xlength = 1.0; + param->ylength = 1.0; + param->imax = 100; + param->jmax = 100; + param->itermax = 1000; + param->eps = 0.0001; + param->omg = 1.8; + param->levels = 5; + param->presmooth = 5; + param->postsmooth = 5; +} + +void readParameter(Parameter* param, const char* filename) +{ + FILE* fp = fopen(filename, "r"); + char line[MAXLINE]; + int i; + + if (!fp) { + fprintf(stderr, "Could not open parameter file: %s\n", filename); + exit(EXIT_FAILURE); + } + + while (!feof(fp)) { + line[0] = '\0'; + fgets(line, MAXLINE, fp); + for (i = 0; line[i] != '\0' && line[i] != '#'; i++) + ; + line[i] = '\0'; + + char* tok = strtok(line, " "); + char* val = strtok(NULL, " "); + +#define PARSE_PARAM(p, f) \ + if (strncmp(tok, #p, sizeof(#p) / sizeof(#p[0]) - 1) == 0) { \ + param->p = f(val); \ + } +#define PARSE_STRING(p) PARSE_PARAM(p, strdup) +#define PARSE_INT(p) PARSE_PARAM(p, atoi) +#define PARSE_REAL(p) PARSE_PARAM(p, atof) + + if (tok != NULL && val != NULL) { + PARSE_REAL(xlength); + PARSE_REAL(ylength); + PARSE_INT(imax); + PARSE_INT(jmax); + PARSE_INT(itermax); + PARSE_REAL(eps); + PARSE_REAL(omg); + PARSE_REAL(re); + PARSE_REAL(tau); + PARSE_REAL(gamma); + PARSE_REAL(dt); + PARSE_REAL(te); + PARSE_REAL(gx); + PARSE_REAL(gy); + PARSE_STRING(name); + PARSE_INT(bcLeft); + PARSE_INT(bcRight); + PARSE_INT(bcBottom); + PARSE_INT(bcTop); + PARSE_INT(levels); + PARSE_INT(presmooth); + PARSE_INT(postsmooth); + PARSE_REAL(u_init); + PARSE_REAL(v_init); + PARSE_REAL(p_init); + } + } + + fclose(fp); +} + +void printParameter(Parameter* param) +{ + printf("Parameters for %s\n", param->name); + printf("Boundary conditions Left:%d Right:%d Bottom:%d Top:%d\n", + param->bcLeft, + param->bcRight, + param->bcBottom, + param->bcTop); + printf("\tReynolds number: %.2f\n", param->re); + printf("\tInit arrays: U:%.2f V:%.2f P:%.2f\n", + param->u_init, + param->v_init, + param->p_init); + printf("Geometry data:\n"); + printf("\tDomain box size (x, y): %.2f, %.2f\n", param->xlength, param->ylength); + printf("\tCells (x, y): %d, %d\n", param->imax, param->jmax); + printf("Timestep parameters:\n"); + printf("\tDefault stepsize: %.2f, Final time %.2f\n", param->dt, param->te); + printf("\tTau factor: %.2f\n", param->tau); + printf("Iterative solver parameters:\n"); + printf("\tMax iterations: %d\n", param->itermax); + printf("\tepsilon (stopping tolerance) : %f\n", param->eps); + printf("\tgamma (stopping tolerance) : %f\n", param->gamma); + printf("\tomega (SOR relaxation): %f\n", param->omg); +} diff --git a/src/parameter.h b/src/parameter.h new file mode 100644 index 0000000..ed765eb --- /dev/null +++ b/src/parameter.h @@ -0,0 +1,27 @@ +/* + * 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. + */ +#ifndef __PARAMETER_H_ +#define __PARAMETER_H_ + +typedef struct { + double xlength, ylength; + int imax, jmax; + int itermax; + double eps, omg; + double re, tau, gamma; + double te, dt; + double gx, gy; + char* name; + int bcLeft, bcRight, bcBottom, bcTop; + double u_init, v_init, p_init; + int levels, presmooth, postsmooth; +} Parameter; + +void initParameter(Parameter*); +void readParameter(Parameter*, const char*); +void printParameter(Parameter*); +#endif diff --git a/src/progress.c b/src/progress.c new file mode 100644 index 0000000..2358e43 --- /dev/null +++ b/src/progress.c @@ -0,0 +1,51 @@ +/* + * 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. + */ +#include +#include +#include +#include +#include +#include "progress.h" + +static double _end; +static int _current; + +void initProgress(double end) +{ + _end = end; + _current = 0; + + printf("[ ]"); + fflush(stdout); +} + +void printProgress(double current) +{ + int new = (int)rint((current / _end) * 10.0); + + if (new > _current) { + char progress[11]; + _current = new; + progress[0] = 0; + + for (int i = 0; i < 10; i++) { + if (i < _current) { + sprintf(progress + strlen(progress), "#"); + } else { + sprintf(progress + strlen(progress), " "); + } + } + printf("\r[%s]", progress); + } + fflush(stdout); +} + +void stopProgress() +{ + printf("\n"); + fflush(stdout); +} diff --git a/src/progress.h b/src/progress.h new file mode 100644 index 0000000..b1937cc --- /dev/null +++ b/src/progress.h @@ -0,0 +1,13 @@ +/* + * 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. + */ +#ifndef __PROGRESS_H_ +#define __PROGRESS_H_ + +extern void initProgress(double); +extern void printProgress(double); +extern void stopProgress(); +#endif diff --git a/src/solver.h b/src/solver.h new file mode 100644 index 0000000..4a0fc20 --- /dev/null +++ b/src/solver.h @@ -0,0 +1,26 @@ +/* + * 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. + */ +#ifndef __SOLVER_H_ +#define __SOLVER_H_ +#include "comm.h" +#include "mpi.h" +#include "parameter.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 diff --git a/src/timing.c b/src/timing.c new file mode 100644 index 0000000..78b01c4 --- /dev/null +++ b/src/timing.c @@ -0,0 +1,22 @@ +/* + * 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 + +double getTimeStamp(void) +{ + struct timespec ts; + clock_gettime(CLOCK_MONOTONIC, &ts); + return (double)ts.tv_sec + (double)ts.tv_nsec * 1.e-9; +} + +double getTimeResolution(void) +{ + struct timespec ts; + clock_getres(CLOCK_MONOTONIC, &ts); + return (double)ts.tv_sec + (double)ts.tv_nsec * 1.e-9; +} diff --git a/src/timing.h b/src/timing.h new file mode 100644 index 0000000..ed05a8c --- /dev/null +++ b/src/timing.h @@ -0,0 +1,13 @@ +/* + * 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. + */ +#ifndef __TIMING_H_ +#define __TIMING_H_ + +extern double getTimeStamp(void); +extern double getTimeResolution(void); + +#endif // __TIMING_H_ diff --git a/src/util.h b/src/util.h new file mode 100644 index 0000000..7218da6 --- /dev/null +++ b/src/util.h @@ -0,0 +1,36 @@ +/* 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. */ +#ifndef __UTIL_H_ +#define __UTIL_H_ + +#ifndef MIN +#define MIN(x, y) ((x) < (y) ? (x) : (y)) +#endif + +#ifndef MAX +#define MAX(x, y) ((x) > (y) ? (x) : (y)) +#endif + +#ifndef ABS +#define ABS(a) ((a) >= 0 ? (a) : -(a)) +#endif + +#define DEBUG_MESSAGE debug_printf + +#ifndef MAXLINE +#define MAXLINE 4096 +#endif + +#define CG_UINT unsigned long long int + +#if PRECISION == 1 +#define CG_FLOAT float +#define PRECISION_STRING "single" +#else +#define CG_FLOAT double +#define PRECISION_STRING "double" +#endif + +#endif