Commit e1e56c68 authored by Stefan Lankes's avatar Stefan Lankes

add template

parent d30be97f
MAKE = make
CC = gcc
CXX = g++
# optimzation flags
#CFLAGS = -Wall -O3
# flags for debugging
CFLAGS = -Wall -O0 -g -pthread
CXXFLAGS = $(CFLAGS)
LDFLAGS =
RM = rm -rf
ASM = nasm
ASMFLAGS = -f elf32 -O0 -g -F dwarf
NAME = jacobi_v1
C_source = main.c init_matrix.c
CPP_source =
ASM_source =
# extend this for other object files
OBJS += $(patsubst %.cpp, %.o, $(filter %.cpp, $(CPP_source)))
OBJS += $(patsubst %.c, %.o, $(filter %.c, $(C_source)))
OBJS += $(patsubst %.asm, %.o, $(filter %.asm, $(ASM_source)))
# other implicit rules
%.o : %.c
$(CC) -c $(CFLAGS) -o $@ $<
%.o : %.cpp
$(CXX) -c $(CXXFLAGS) -o $@ $<
%.o : %.asm
$(ASM) $(ASMFLAGS) -o $@ $<
default:
$(MAKE) $(NAME)
all:
$(MAKE) $(NAME)
$(NAME): $(OBJS)
$(CC) $(CFLAGS) -o $(NAME) $(OBJS) $(LDFLAGS)
clean:
$(RM) *.o $(NAME)
depend:
$(CC) -MM $(CFLAGS) $(C_source) $(CPP_source) > Makefile.dep
-include Makefile.dep
# DO NOT DELETE
#include "init_matrix.h"
#include <stdlib.h>
#include <time.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#define MAXVALUE 2.0
int generate_empty_matrix(double ***A, double **b, unsigned int N)
{
unsigned int iCnt;
if ((A == NULL) || (b == NULL)) /* No memory allocated */
return -1; /* Error */
*A = (double **)malloc(N * sizeof(double *));
if (*A == NULL)
return -2; /* Error */
**A = (double *)malloc(N * N * sizeof(double));
if (**A == NULL) {
free(*A); /* Clean up */
*A = NULL;
return -2; /* Error */
}
for (iCnt = 1; iCnt < N; iCnt++) { /* Assign pointers in the first "real index"; Value from 1 to N (0 yet set, value N means N+1) */
(*A)[iCnt] = &((*A)[0][iCnt * N]);
}
memset(**A, 0, N * N * sizeof(double)); /* Fill matrix values with 0 */
*b = (double *)malloc(N * sizeof(double));
if (*b == NULL) {
free(**A);
free(*A);
*A = NULL;
*b = NULL;
return -2; /* Error */
}
memset(*b, 0, N * sizeof(double)); /* Fill vector values with 0 */
return 0;
}
void clean_matrix(double ***A)
{
if (A != NULL) {
if (*A != NULL) {
if (**A != NULL)
free(**A);
free(*A);
*A = NULL;
}
}
}
/* frees the (result) vector */
void clean_vector(double **b)
{
if ((b != NULL) && (*b != NULL)) {
free(*b);
*b = NULL;
}
}
int generate_base_matrix(double ***A, double **b, unsigned int N)
{
int errorcode;
errorcode = generate_empty_matrix(A, b, N);
if (errorcode < 0) { /* Forward the errorcode */
return errorcode;
} else {
unsigned int i;
for (i = 0; i < N; i++) {
(*b)[i] = 0.0;
(*A)[i][i] = 1.0; /* Fill diag of identity matrix */
}
return 0;
}
}
int init_matrix(double ***A, double **b, unsigned int N)
{
unsigned int j, i, errorcode;
errorcode = generate_empty_matrix(A, b, N);
if (errorcode < 0) /* forward errorcode */
return errorcode;
srand((unsigned)time(NULL)); /* init random number generator */
/*
* initialize the system of linear equations
* the result vector is one
*/
for (i = 0; i < N; i++) {
double sum = 0.0;
(*b)[i] = 0.0f;
for (j = 0; j < N; j++) {
if (i != j) {
double c =
((double)rand()) / ((double)RAND_MAX) *
MAXVALUE;
sum += fabs(c);
(*A)[i][j] = c;
(*b)[i] += c;
}
}
/*
* The Jacobi method will always converge if the matrix A is strictly or irreducibly diagonally dominant.
* Strict row diagonal dominance means that for each row, the absolute value of the diagonal term is
* greater than the sum of absolute values of other terms: |A[i][i]| > Sum |A[i][j]| with (i != j)
*/
(*A)[i][i] = sum + 2.0;
(*b)[i] += sum + 2.0;
}
return 0;
}
/*
* Errorcodes in all integer functions:
*
* -1: Given pointer is 0 (can't work on it)
* -2: Allocation error
*
*/
#ifndef __INIT_MATRIX_H__
#define __INIT_MATRIX_H__
#ifdef __cplusplus
extern "C" {
#endif
/*
* Creates identity matrix (NxN) and an empty result vector
* simply solveable
*
* *********
* * 1 0 *
* * 1 *
* * 0 1 *
* *********
*/
int generate_base_matrix(double ***A, double **b, unsigned int N);
/* Creates empty NxN matrix and empty result vector */
int generate_empty_matrix(double ***A, double **b, unsigned int N);
/* creates a solveable system of linear equations */
int init_matrix(double ***A, double **b, unsigned int N);
/* frees all pointers below upper level of A */
void clean_matrix(double ***A);
/* frees the (result) vector */
void clean_vector(double **b);
#ifdef __cplusplus
}
#endif
#endif
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <string.h>
#include <math.h>
#include <sys/time.h>
#include "init_matrix.h"
#define MATRIX_SIZE (1024)
double **A;
double *b;
double *X;
double *X_old;
double *temp;
int main(int argc, char **argv)
{
unsigned int i, j;
unsigned int iterations = 0;
double error, xi, norm, max = 0.0;
struct timeval start, end;
printf("\nInitialize system of linear equations...\n");
/* allocate memory for the system of linear equations */
init_matrix(&A, &b, MATRIX_SIZE);
X = (double *)malloc(sizeof(double) * MATRIX_SIZE);
X_old = (double *)malloc(sizeof(double) * MATRIX_SIZE);
/* a "random" solution vector */
for (i = 0; i < MATRIX_SIZE; i++) {
X[i] = ((double)rand()) / ((double)RAND_MAX) * 10.0;
X_old[i] = 0.0;
}
printf("Start Jacobi method...\n");
gettimeofday(&start, NULL);
/* TODO: Hier muss die Aufgabe geloest werden */
gettimeofday(&end, NULL);
if (MATRIX_SIZE < 16) {
printf("Print the solution...\n");
/* print solution */
for (i = 0; i < MATRIX_SIZE; i++) {
for (j = 0; j < MATRIX_SIZE; j++)
printf("%8.2f\t", A[i][j]);
printf("*\t%8.2f\t=\t%8.2f\n", X[i], b[i]);
}
}
printf("Check the result...\n");
/*
* check the result
* X[i] have to be 1
*/
for (i = 0; i < MATRIX_SIZE; i++) {
error = fabs(X[i] - 1.0f);
if (max < error)
max = error;
if (error > 0.01f)
printf("Result is on position %d wrong (%f != 1.0)\n",
i, X[i]);
}
printf("maximal error is %f\n", max);
printf("\nmatrix size: %d x %d\n", MATRIX_SIZE, MATRIX_SIZE);
printf("number of iterations: %d\n", iterations);
printf("Time : %lf sec\n",
(double)(end.tv_sec - start.tv_sec) + (double)(end.tv_usec -
start.tv_usec) /
1000000.0);
/* frees the allocated memory */
free(X_old);
free(X);
clean_matrix(&A);
clean_vector(&b);
return 0;
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment