add jacobi solver as example program
This commit is contained in:
Stefan Lankes
parent
3fea08710d
commit
6e255fe27e
2 changed files with 208 additions and 1 deletions
|
|
@ -3,6 +3,7 @@ NEWLIB = ../x86/i586-metalsvm-elf32
|
|||
MAKE = make
|
||||
STRIP_DEBUG = --strip-debug
|
||||
KEEP_DEBUG = --only-keep-debug
|
||||
LDFLAGS =
|
||||
|
||||
# other implicit rules
|
||||
%.o : %.c
|
||||
|
|
@ -10,7 +11,13 @@ KEEP_DEBUG = --only-keep-debug
|
|||
|
||||
default: all
|
||||
|
||||
all: hello tests
|
||||
all: hello tests jacobi
|
||||
|
||||
jacobi: jacobi.o
|
||||
$(CC_FOR_TARGET) -T link.ld -o $@ $(LDFLAGS) $< -lm
|
||||
$(OBJCOPY_FOR_TARGET) $(KEEP_DEBUG) $@ $@.sym
|
||||
$(OBJCOPY_FOR_TARGET) $(STRIP_DEBUG) $@
|
||||
chmod a-x $@.sym
|
||||
|
||||
tests: tests.o
|
||||
$(CC_FOR_TARGET) -T link.ld -o $@ $(LDFLAGS) $<
|
||||
|
|
|
|||
200
newlib/examples/jacobi.c
Normal file
200
newlib/examples/jacobi.c
Normal file
|
|
@ -0,0 +1,200 @@
|
|||
/*
|
||||
* Copyright 2011 Stefan Lankes, Alexander Pilz, Maximilian Marx, Michael Ober,
|
||||
* Chair for Operating Systems, RWTH Aachen University
|
||||
*
|
||||
* Licensed under the Apache License, Version 2.0 (the "License");
|
||||
* you may not use this file except in compliance with the License.
|
||||
* You may obtain a copy of the License at
|
||||
*
|
||||
* http://www.apache.org/licenses/LICENSE-2.0
|
||||
*
|
||||
* Unless required by applicable law or agreed to in writing, software
|
||||
* distributed under the License is distributed on an "AS IS" BASIS,
|
||||
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
* See the License for the specific language governing permissions and
|
||||
* limitations under the License.
|
||||
*
|
||||
*/
|
||||
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include <math.h>
|
||||
#include <time.h>
|
||||
#include <unistd.h>
|
||||
#include <errno.h>
|
||||
#undef errno
|
||||
extern int errno;
|
||||
|
||||
#define MATRIX_SIZE 256
|
||||
#define MAXVALUE 1337
|
||||
#define PAGE_SIZE 4096
|
||||
#define CACHE_SIZE (256*1024)
|
||||
#define ALIGN(x,a) (((x)+(a)-1)&~((a)-1))
|
||||
|
||||
static int generate_empty_matrix(double*** A , unsigned int N) {
|
||||
unsigned int iCnt;
|
||||
int i,j;
|
||||
|
||||
*A = (double**) malloc((N+1)*sizeof(double*));
|
||||
|
||||
if (*A == NULL)
|
||||
return -2; /* Error */
|
||||
|
||||
(*A)[0] = (double*) malloc((N+1)*N*sizeof(double));
|
||||
|
||||
if (**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+1)]);
|
||||
}
|
||||
|
||||
memset(**A, 0, (N+1)*N*sizeof(double)); /* Fill matrix values with 0 */
|
||||
|
||||
srand( 42 /*(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;
|
||||
|
||||
for (j = 0; j < N; j++)
|
||||
{
|
||||
if (i != j)
|
||||
{
|
||||
double c = ((double)rand()) / ((double)RAND_MAX) * MAXVALUE;
|
||||
|
||||
sum += fabs(c);
|
||||
(*A)[i][j] = c;
|
||||
(*A)[i][N] += 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;
|
||||
(*A)[i][N] += sum + 2.0;
|
||||
}
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
int main(int argc, char **argv)
|
||||
{
|
||||
double* temp;
|
||||
unsigned int i, j, iter_start, iter_end;
|
||||
unsigned int iterations = 0;
|
||||
double error, norm, norm_res, max = 0.0;
|
||||
double** A=0;
|
||||
double* X;
|
||||
double* X_old, xi;
|
||||
double start,stop;
|
||||
|
||||
if (generate_empty_matrix(&A,MATRIX_SIZE) < 0)
|
||||
{
|
||||
printf("generate_empty_matrix() failed...\n");
|
||||
fflush(stdout);
|
||||
exit(-1);
|
||||
|
||||
}
|
||||
|
||||
printf("generate_empty_matrix() done...\n");
|
||||
fflush(stdout);
|
||||
|
||||
X=(double*) malloc(MATRIX_SIZE*sizeof(double));
|
||||
X_old=(double*) malloc(MATRIX_SIZE*sizeof(double));
|
||||
if(X == NULL || X_old == NULL)
|
||||
{
|
||||
printf("X or X_old is NULL...\n");
|
||||
exit(-1);
|
||||
}
|
||||
|
||||
for(i=0; i<MATRIX_SIZE; i++)
|
||||
{
|
||||
X[i] = ((double)rand()) / ((double)RAND_MAX) * 10.0;
|
||||
X_old[i] = 0.0;
|
||||
}
|
||||
|
||||
printf("start calculation...\n");
|
||||
fflush(stdout);
|
||||
|
||||
iter_start = 0;
|
||||
iter_end = MATRIX_SIZE;
|
||||
|
||||
//start = RCCE_wtime();
|
||||
|
||||
while(1)
|
||||
{
|
||||
iterations++;
|
||||
|
||||
temp = X_old;
|
||||
X_old = X;
|
||||
X = temp;
|
||||
|
||||
for (i=iter_start; i<iter_end; i++)
|
||||
{
|
||||
for(j=0, xi=0.0; j<i; j++)
|
||||
xi += A[i][j]* X_old[j];
|
||||
|
||||
for(j=i+1; j<MATRIX_SIZE; j++)
|
||||
xi += A[i][j] * X_old[j];
|
||||
X[i] = (A[i][MATRIX_SIZE] - xi) / A[i][i];
|
||||
}
|
||||
|
||||
if (iterations % 5000 == 0 ) {/* calculate the Euclidean norm between X_old and X*/
|
||||
norm_res = norm = 0.0;
|
||||
for (i=iter_start; i<iter_end; i++)
|
||||
norm += (X_old[i] - X[i]) * (X_old[i] - X[i]);
|
||||
|
||||
/* check the break condition */
|
||||
norm_res /= (double) MATRIX_SIZE;
|
||||
|
||||
if (norm_res < 0.0000001)
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
//stop = RCCE_wtime();
|
||||
|
||||
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], A[i][MATRIX_SIZE]);
|
||||
}
|
||||
}
|
||||
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("calculation time: %f s\n", stop-start);
|
||||
|
||||
free((void*) X_old);
|
||||
free((void*) X);
|
||||
|
||||
return 0;
|
||||
}
|
||||
Loading…
Add table
Reference in a new issue