+++ /dev/null
-/* atlc - arbitrary transmission line calculator, for the analysis of
-transmission lines are directional couplers.
-
-Copyright (C) 2002. Dr. David Kirkby, PhD (G8WRB).
-
-This program is free software; you can redistribute it and/or
-modify it under the terms of the GNU General Public License
-as published by the Free Software Foundation; either package_version 2
-of the License, or (at your option) any later package_version.
-
-This program is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-GNU General Public License for more details.
-
-You should have received a copy of the GNU General Public License
-along with this program; if not, write to the Free Software
-Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
-USA.
-
-Dr. David Kirkby, e-mail drkirkby at ntlworld.com
-
-*/
-
-#include "config.h"
-
-#define SLOW
-
-#ifndef ENABLE_POSIX_THREADS
-
-/* If we are not comiling for multi-threaded use, a function
-is defined, but it is not used. This keeps the linker happy */
-
-double finite_difference_multi_threaded()
-{
- return(0.0);
-}
-#endif
-
-
-#ifdef ENABLE_POSIX_THREADS
-
-#ifdef HAVE_SYS_TYPES_H
-#include <sys/types.h>
-#endif
-
-#ifdef HAVE_STDLIB_H
-#include <stdlib.h>
-#endif
-
-#include <pthread.h>
-
-#ifdef HAVE_SCHED_H
-#include <sched.h>
-#endif
-
-#ifdef HAVE_STRING_H
-#include <string.h>
-#endif
-
-#include "definitions.h"
-#include "exit_codes.h"
-
-extern int coupler;
-extern int width, height;
-extern double **Vij, **Er;
-extern signed char **oddity;
-extern int number_of_workers;
-
-
-/* The algorithm for this is based on one in the book 'Foundations of Multithraded,
-Parallel and Distributed Programming' by G. A. Andrews, (2000). Especially
-chapter 11, Fig 11.6, "Red/Black Gauss-Seidel using shared vairables'.
-
-Basically the array is a considered a chess board. Since we only use in the computation the
-values above, below, to the left and the right of the current pixel, these will all be the
-same colour on the chess board (black or white). So we compute all the white squares, which
-can all be done in parallel, since they don't depend on each other. Once the white squares
-are done, the black ones can be done. Again, these don't depend on the white squares, but
-only the black ones. So the algorithm does
-
-1) Covers the array with a large chess board with black and white squares.
-
-
-B W B W B W B W B W B W B W B W B W B W B W B W B W B W
-W B W B W B W B W B W B W B W B W B W B W B W B W B W B
-B W B W B W B W B W B W B W B W B W B W B W B W B W B W
-W B W B W B W B W B W B W B W B W B W B W B W B W B W B
-B W B W B W B W B W B W B W B W B W B W B W B W B W B W
-W B W B W B W B W B W B W B W B W B W B W B W B W B W B
-B W B W B W B W B W B W B W B W B W B W B W B W B W B W
-W B W B W B W B W B W B W B W B W B W B W B W B W B W B
-B W B W B W B W B W B W B W B W B W B W B W B W B W B W
-W B W B W B W B W B W B W B W B W B W B W B W B W B W B
-B W B W B W B W B W B W B W B W B W B W B W B W B W B W
-W B W B W B W B W B W B W B W B W B W B W B W B W B W B
-
-
-2) Split the array into a number of columns, one for each CPU.
-
-B W B W B W B W B W B W B W B W B W B W B W B W B W B W
-W B W B W B W B W B W B W B W B W B W B W B W B W B W B
-B W B W B W B W B W B W B W B W B W B W B W B W B W B W
-W B W B W B W B W B W B W B W B W B W B W B W B W B W B
-B W B W B W B W B W B W B W B W B W B W B W B W B W B W
-W B W B W B W B W B W B W B W B W B W B W B W B W B W B
-B W B W B W B W B W B W B W B W B W B W B W B W B W B W
-W B W B W B W B W B W B W B W B W B W B W B W B W B W B
-B W B W B W B W B W B W B W B W B W B W B W B W B W B W
-W B W B W B W B W B W B W B W B W B W B W B W B W B W B
-B W B W B W B W B W B W B W B W B W B W B W B W B W B W
-W B W B W B W B W B W B W B W B W B W B W B W B W B W B
-
-3) Do all columns in parallel on the black squares.
-
-B B B B B B B B B B B B B B
- B B B B B B B B B B B B B B
-B B B B B B B B B B B B B B
- B B B B B B B B B B B B B B
-B B B B B B B B B B B B B B
- B B B B B B B B B B B B B B
-B B B B B B B B B B B B B B
- B B B B B B B B B B B B B B
-B B B B B B B B B B B B B B
- B B B B B B B B B B B B B B
-B B B B B B B B B B B B B B
- B B B B B B B B B B B B B B
-
-Thread 0 Thread 1 Thread 2 Thread 3
-
-4) Wait until all the black squares are done.
-
-5) Compute all the white squares in parallel.
-
- W W W W W W W W W W W W W W
-W W W W W W W W W W W W W W
- W W W W W W W W W W W W W W
-W W W W W W W W W W W W W W
- W W W W W W W W W W W W W W
-W W W W W W W W W W W W W W
- W W W W W W W W W W W W W W
-W W W W W W W W W W W W W W
- W W W W W W W W W W W W W W
-W W W W W W W W W W W W W W
- W W W W W W W W W W W W W W
-W W W W W W W W W W W W W W
-Thread 0 Thread 1 Thread 2 Thread 3
-
-
-6) Repeat a large number (typically 100) times.
-
-7) Check for convergence.
-
-8) Stop if the solution has converged.
-*/
-
-extern double r;
-
-pthread_mutex_t barrier; /* mutex semaphore for the barrier */
-pthread_cond_t go; /* condition variable for leaving */
-int numArrived = 0; /* count of the number who have arrived */
-
-#define BARRIER2
-
-#ifdef BARRIER1
-void Barrier() {
-#ifndef DEBUG
- pthread_mutex_lock(&barrier);
-#else
- if( pthread_mutex_lock(&barrier) != 0 )
- exit_with_msg_and_exit_code("pthread_mutex_lock failed",PTHREAD_MUTEX_LOCK_FAILED);
-#endif
-
- numArrived++;
- if (numArrived == number_of_workers) {
- numArrived = 0;
- pthread_cond_broadcast(&go);
- } else
- pthread_cond_wait(&go, &barrier);
-#ifndef DEBUG
- pthread_mutex_unlock(&barrier);
-#else
- if( pthread_mutex_unlock(&barrier) != 0 )
- exit_with_msg_and_exit_code("pthread_mutex_unlock failed",PTHREAD_MUTEX_UNLOCK_FAILED);
-#endif
-}
-#endif
-
-#ifdef BARRIER2
-int numArrived2[2] = {0, 0};
-int barrierNdx = 0;
-
-void Barrier() {
- int localNdx;
-#ifndef DEBUG
- pthread_mutex_lock(&barrier);
-#else
- if( pthread_mutex_lock(&barrier) != 0 )
- exit_with_msg_and_exit_code("pthread_mutex_lock failed in finite_difference_multi_threaded.c",PTHREAD_MUTEX_LOCK_FAILED);
-#endif
-
- numArrived2[barrierNdx]++;
- if (numArrived2[barrierNdx] == number_of_workers) {
- barrierNdx = (barrierNdx + 1)%2; /* toggle */
- numArrived2[barrierNdx] = 0; /* reset other count */
-
-#ifndef DEBUG
- pthread_cond_broadcast(&go);
-#else
- if( pthread_cond_broadcast(&go) != 0)
- exit_with_msg_and_exit_code("pthread_cond_broadcast failed in finite_difference_multi_threaded.c",PTHREAD_COND_BROADCAST_FAILED);
-#endif
- }
- else
- {
- localNdx = barrierNdx; /* wait on "current" numArrived. */
- while (numArrived2[localNdx] != number_of_workers)
-#ifndef DEBUG
- pthread_cond_wait(&go, &barrier);
-#else
- if( pthread_cond_wait(&go, &barrier) != 0)
- exit_with_msg_and_exit_code("pthread_cond_wait failed finite_difference_multi_threaded.c",PTHREAD_COND_WAIT_FAILED);
-#endif
- }
-#ifndef DEBUG
- pthread_mutex_unlock(&barrier);
-#else
- if( pthread_mutex_unlock(&barrier) != 0 )
- exit_with_msg_and_exit_code("pthread_mutex_unlock failed finite_difference_multi_threaded.c",PTHREAD_MUTEX_UNLOCK_FAILED);
-#endif
-}
-#endif
-
-
-
-/* Each Worker computes values in one strip of the grids.
- The main worker loop does two computations to avoid copying from
- one grid to the other. */
-
-void *worker(void *arg) {
- int myid = (int) arg;
- double r_over_4,a,b,c,d,e,f,g,h;
- double one_minus_r;
- int i, j, iters, jstart;
- int firstcol, lastcol;
- int strip_width=width/number_of_workers;
-
- firstcol = myid*strip_width+1;
- lastcol = firstcol + strip_width - 1;
- r_over_4=r/4.0;
- one_minus_r=1-r;
- if(myid == number_of_workers -1)
- lastcol=width-2;
- Barrier();
- for (iters = 1; iters <= ITERATIONS; iters++)
- {
- for(i= firstcol ; i <= lastcol; ++i){
- if(i%2 ==1 )
- jstart=1;
- else
- jstart=2;
- for(j=jstart ; j < height-1 ; j+=2){
- if(oddity[i][j] == ORDINARY_INTERIOR_POINT) /* Same dielectric all around */
- Vij[i][j]=r_over_4*(Vij[i][j+1]+Vij[i+1][j]+Vij[i][j-1]+Vij[i-1][j])+one_minus_r*Vij[i][j];
- else if(oddity[i][j] > ORDINARY_INTERIOR_POINT) /* only update dielectrics, not conductors */
- {
- a=(Er[i][j] * Er[i][j-1] * Vij[i][j-1])/(Er[i][j] + Er[i][j-1]);
- b=(Er[i][j] * Er[i][j+1] * Vij[i][j+1])/(Er[i][j] + Er[i][j+1]);
- c=(Er[i][j] * Er[i-1][j] * Vij[i-1][j])/(Er[i][j] + Er[i-1][j]);
- d=(Er[i][j] * Er[i+1][j] * Vij[i+1][j])/(Er[i][j] + Er[i+1][j]);
-
- e=(Er[i][j] * Er[i][j-1])/(Er[i][j]+Er[i][j-1]);
- f=(Er[i][j] * Er[i][j+1])/(Er[i][j]+Er[i][j+1]);
- g=(Er[i][j] * Er[i-1][j])/(Er[i][j]+Er[i-1][j]);
- h=(Er[i][j] * Er[i+1][j])/(Er[i][j]+Er[i+1][j]);
-
- Vij[i][j]=r*(a+b+c+d)/(e+f+g+h) + (1-r)*Vij[i][j];
- }
- }
- }
- Barrier();
- for(i= firstcol ; i <= lastcol; ++i){
- if(i%2 ==1 )
- jstart=2;
- else
- jstart=1;
- for(j=jstart ; j < height -1; j+=2){
- if(oddity[i][j] == ORDINARY_INTERIOR_POINT) /* Same dielectric all around */
- Vij[i][j]=r_over_4*(Vij[i][j+1]+Vij[i+1][j]+Vij[i][j-1]+Vij[i-1][j])+one_minus_r*Vij[i][j];
- else if(oddity[i][j] > ORDINARY_INTERIOR_POINT) /* only update dielectrics, not conductors */
- {
- a=(Er[i][j] * Er[i][j-1] * Vij[i][j-1])/(Er[i][j] + Er[i][j-1]);
- b=(Er[i][j] * Er[i][j+1] * Vij[i][j+1])/(Er[i][j] + Er[i][j+1]);
- c=(Er[i][j] * Er[i-1][j] * Vij[i-1][j])/(Er[i][j] + Er[i-1][j]);
- d=(Er[i][j] * Er[i+1][j] * Vij[i+1][j])/(Er[i][j] + Er[i+1][j]);
-
- e=(Er[i][j] * Er[i][j-1])/(Er[i][j]+Er[i][j-1]);
- f=(Er[i][j] * Er[i][j+1])/(Er[i][j]+Er[i][j+1]);
- g=(Er[i][j] * Er[i-1][j])/(Er[i][j]+Er[i-1][j]);
- h=(Er[i][j] * Er[i+1][j])/(Er[i][j]+Er[i+1][j]);
-
- Vij[i][j]=r*(a+b+c+d)/(e+f+g+h) + (1-r)*Vij[i][j];
- }
- }
- }
- Barrier();
- }
- Barrier();
- return(0);
-}
-
-double finite_difference_multi_threaded()
-{
- int i, j, ret, thread_number;
- pthread_t thread_id[MAXIMUM_PROCESSING_DEVICES];
-
-
- double capacitance_per_metre, energy_per_metre;
-
-
- /* initialize mutex and condition variable */
- pthread_mutex_init(&barrier, NULL);
- pthread_cond_init(&go, NULL);
- for(thread_number=0;thread_number<number_of_workers;thread_number++){
- ret=pthread_create(&thread_id[thread_number],NULL, worker,(void *)thread_number);
- if(ret != 0)
- exit_with_msg_and_exit_code("failed to create thread in finite_difference_multi_threaded.c",THREAD_CREATION_FAILED);
- }
- /* Wait for each thread to join - i.e. once they are all finished. */
- for(thread_number=0;thread_number<number_of_workers;thread_number++)
- {
- ret=pthread_join(thread_id[thread_number], NULL);
- if(ret != 0)
- exit_with_msg_and_exit_code("failed to join thread in finite_difference_multi_threaded.c",THREAD_FAILED_TO_JOIN);
- }
- energy_per_metre=0.0;
- /*
- The following code is wrong, but worked on all systems apart from AIX !!
-
- for(i=0;i<width;++i)
- for(j=0;j<height;++j)
- energy_per_metre+=find_energy_per_metre(i,j);
- */
-
- for(i=1;i<width-1;++i)
- for(j=1;j<height-1;++j)
- energy_per_metre+=find_energy_per_metre(i,j);
-
- if(coupler==FALSE)
- capacitance_per_metre=2*energy_per_metre;
- else
- capacitance_per_metre=energy_per_metre;
- return(capacitance_per_metre);
-}
-#endif
projects / debian / atlc / blobdiff