moebiusband/CG-Bench
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Initial checkin

Browse Source
This commit is contained in:
Jan Eitzinger 2025-01-03 14:23:28 +01:00
parent fd4a93c39d
commit d355c4bbdb
25 changed files with 1659 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

176
.clang-format Normal file
View File

@ -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
...

14
.clang-tidy Normal file
View File

@ -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'

2
.gitignore vendored Normal file
View File

@ -0,0 +1,2 @@
/build
.clangd

78
Makefile Normal file
View File

@ -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)

9
config.mk Normal file
View File

@ -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

15
mk/include_CLANG.mk Normal file
View File

@ -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

14
mk/include_GCC.mk Normal file
View File

@ -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 =

14
mk/include_ICC.mk Normal file
View File

@ -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 =

38
src/allocate.c Normal file
View File

@ -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 <errno.h>
#include <stddef.h>
#include <stdio.h>
#include <stdlib.h>
#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;
}

11
src/allocate.h Normal file
View File

@ -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 <stdlib.h>
extern void* allocate(size_t alignment, size_t bytesize);
#endif

129
src/comm.c Normal file
View File

@ -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 <stdio.h>
#include <stdlib.h>
#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
}

55
src/comm.h Normal file
View File

@ -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 <mpi.h>
#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_

90
src/main.c Normal file
View File

@ -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 <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#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 <configFile>\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;
}

130
src/matrix.c Normal file
View File

@ -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 <stdbool.h>
#include <stddef.h>
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#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);
}

24
src/matrix.h Normal file
View File

@ -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_

429
src/mmio.c Normal file
View File

@ -0,0 +1,429 @@
/*
* Matrix Market I/O library for ANSI C
*
* See http://math.nist.gov/MatrixMarket for details.
*
*
*/
#include <ctype.h>
#include <stdlib.h>
#include <string.h>
#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);
}

129
src/mmio.h Normal file
View File

@ -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 <stdio.h>
#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

114
src/parameter.c Normal file
View File

@ -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 <stdio.h>
#include <stdlib.h>
#include <string.h>
#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);
}

27
src/parameter.h Normal file
View File

@ -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

51
src/progress.c Normal file
View File

@ -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 <math.h>
#include <mpi.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#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);
}

13
src/progress.h Normal file
View File

@ -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

26
src/solver.h Normal file
View File

@ -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

22
src/timing.c Normal file
View File

@ -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 <stdlib.h>
#include <time.h>
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;
}

13
src/timing.h Normal file
View File

@ -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_

36
src/util.h Normal file
View File

@ -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