/* Compilation: gcc -o exam1 exam1.c -lpthread */

// Параллельное умножение квадратной матрицы А размерностью NxN
// на вектор B, результат - вектор C.
// Создается nt потоков, каждый из которых умножает ленту матрицы
// A размерностью (N/nt)xN на вектор B.
// Умножение повторяется nc раз при неизменной матрице A и
// различных векторах B.

#include <stdlib.h>
#include <unistd.h>
#include <stdio.h>
#include <string.h>
#include <signal.h>
#define _REENTRANT
#include <pthread.h>
#include <sched.h>

typedef struct {
  pthread_t tid;
  int first;
  int last;
  pthread_mutex_t mutx;
  pthread_cond_t condx;
  } ThreadRecord;

int N;	// Размер матрицы A и векторов B, C
double *A, *B, *C;

int done=0;
ThreadRecord *threads; // Массив мьютексов и условн. переменных

void *
mysolver (void *arg_p) {
  ThreadRecord *thr;
  int i, j;

  thr = (ThreadRecord *)arg_p;
  pthread_mutex_lock( &(thr->mutx) );

  while ( 1 ) {
fprintf (stderr, "\t%d\n", thr->tid);
    for (i=thr->first; i<=thr->last; i++) {
      C[i] = 0;
      for (j=0; j<=N; j++)
        C[i] += A[N*i+j]*B[j];
      };
    pthread_cond_wait(&(thr->condx), &(thr->mutx));
    };
  }

int
main (int argc, char* argv[]) {
  int nt;	// Кол-во расчетных потоков
  int nc;	// Кол-во циклов
  int i, j, k;

  pthread_attr_t pattr;

  if (argc != 4) {
    fprintf(stderr, "Исп: %s размерность #потоков #циклов\n", argv[0]);
    exit (1);
    };

  N = atoi(argv[1]);
  nt = atoi(argv[2]);
  nc = atoi(argv[3]);

  if ( N%nt ) {
    fprintf(stderr, "Размерность д.б. кратна кол-ву потоков\n", argv[0]);
    exit (2);
    };

  // Выделение памяти
  A = (double *) malloc( sizeof(double)*N*N );
  B = (double *) malloc( sizeof(double)*N );
  C = (double *) malloc( sizeof(double)*N );

  // Инициализация A и B случ. значениями
  srand(time(NULL));
  for (i=0; i<N; i++) {
    B[i] = (double) rand();
    for (j=0; j<N; j++)
      A[i*N+j] = (double) rand();
    };

  pthread_attr_init (&pattr);
  pthread_attr_setscope (&pattr, PTHREAD_SCOPE_SYSTEM);
  pthread_attr_setdetachstate (&pattr,PTHREAD_CREATE_JOINABLE);

  threads = (ThreadRecord *) calloc (nt, sizeof(ThreadRecord));

  j = N/nt;
  for (i=0; i<nt; i++) {
    threads[i].first = j*i; // Индекс нач. строки А для потока
    threads[i].last = j*(i+1)-1; // Индекс конечн. строки А для потока
    pthread_mutex_init (&(threads[i].mutx), NULL);
    pthread_cond_init (&(threads[i].condx), NULL);
    if ( pthread_create (&(threads[i].tid), &pattr, mysolver, (void *) &(threads[i])) )
      perror("pthread_create");
    sched_yield(); // Попытка стартовать поток как можно быстрее
    };

  for (k=1; k<nc; k++) {
    usleep(100000); // Попытка дать шанс расчетному потоку
                    // вернуть себе мьютекс - весьма критичная величина!!!!!!!!!!!!!!!!
fprintf(stderr, "---------\n");
    for (i=0; i<nt; i++) {
      pthread_mutex_lock( &(threads[i].mutx) );
      };
    for (j=0; j<N; j++)
        B[j] = (double) rand();

    for (i=0; i<nt; i++) {
      pthread_mutex_unlock( &(threads[i].mutx) );
      pthread_cond_signal( &(threads[i].condx) );
      sched_yield();
      };
    };

  exit (0);
  }