1/*
2 * Copyright (c) 2004, Bull S.A..  All rights reserved.
3 * Created by: Sebastien Decugis
4
5 * This program is free software; you can redistribute it and/or modify it
6 * under the terms of version 2 of the GNU General Public License as
7 * published by the Free Software Foundation.
8 *
9 * This program is distributed in the hope that it would be useful, but
10 * WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
12 *
13 * You should have received a copy of the GNU General Public License along
14 * with this program; if not, write the Free Software Foundation, Inc.,
15 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
16 *
17
18 * This file is a scalability test for the pthread_cond_timedwait function.
19 *
20 * It aims to measure the time between end of timeout and actual wakeup
21 * with different factors.
22
23 * The steps are:
24 * -*- Number of threads waiting on the conditionnal variable
25 *     -> for an increaing number of threads,
26 *        -> create the other threads which will do a pthread_cond_wait on the same cond/mutex
27 *        -> When the threads are waiting, create one thread which will measure the time
28 *        -> once the timeout has expired and the measure is done, broadcast the condition.
29 *     -> do each measure 10 times (with different attributes i.e.)
30 *
31 * -*- other possible influencial parameters
32 *     -> To be defined.
33 */
34
35 /* We are testing conformance to IEEE Std 1003.1, 2003 Edition */
36#define _POSIX_C_SOURCE 200112L
37
38#ifndef WITHOUT_XOPEN
39#define _XOPEN_SOURCE 600
40#endif
41
42/********************************************************************************************/
43/****************************** standard includes *****************************************/
44/********************************************************************************************/
45#include <pthread.h>
46#include <stdarg.h>
47#include <stdio.h>
48#include <stdlib.h>
49#include <unistd.h>
50
51#include <time.h>
52#include <errno.h>
53#include <math.h>
54
55/********************************************************************************************/
56/******************************   Test framework   *****************************************/
57/********************************************************************************************/
58#include "testfrmw.h"
59#include "testfrmw.c"
60 /* This header is responsible for defining the following macros:
61  * UNRESOLVED(ret, descr);
62  *    where descr is a description of the error and ret is an int (error code for example)
63  * FAILED(descr);
64  *    where descr is a short text saying why the test has failed.
65  * PASSED();
66  *    No parameter.
67  *
68  * Both three macros shall terminate the calling process.
69  * The testcase shall not terminate in any other maneer.
70  *
71  * The other file defines the functions
72  * void output_init()
73  * void output(char * string, ...)
74  *
75  * Those may be used to output information.
76  */
77
78/********************************************************************************************/
79/********************************** Configuration ******************************************/
80/********************************************************************************************/
81#ifndef SCALABILITY_FACTOR
82#define SCALABILITY_FACTOR 1
83#endif
84#ifndef VERBOSE
85#define VERBOSE 1
86#endif
87
88#ifndef WITHOUT_ALTCLK
89#define USE_ALTCLK		/* make tests with MONOTONIC CLOCK if supported */
90#endif
91
92#define MES_TIMEOUT  (1000000)	/* ns, offset for the pthread_cond_timedwait call */
93
94#ifdef PLOT_OUTPUT
95#undef VERBOSE
96#define VERBOSE 0
97#endif
98
99// #define USE_CANCEL  /* Will cancel the threads instead of signaling the cond */#define
100
101/********************************************************************************************/
102/***********************************    Test case   *****************************************/
103/********************************************************************************************/
104
105long altclk_ok, pshared_ok;
106
107typedef struct {
108	pthread_cond_t *cnd;
109	pthread_mutex_t *mtx;
110	int *predicate;
111	int *tnum;
112} test_t;
113
114struct {
115	int mutex_type;
116	int pshared;
117	clockid_t cid;
118	char *desc;
119} test_scenar[] = {
120	{
121	PTHREAD_MUTEX_DEFAULT, PTHREAD_PROCESS_PRIVATE, CLOCK_REALTIME,
122		    "Default"}
123	, {
124	PTHREAD_MUTEX_DEFAULT, PTHREAD_PROCESS_SHARED, CLOCK_REALTIME,
125		    "PShared"}
126#ifndef WITHOUT_XOPEN
127	, {
128	PTHREAD_MUTEX_NORMAL, PTHREAD_PROCESS_PRIVATE, CLOCK_REALTIME,
129		    "Normal"}
130	, {
131	PTHREAD_MUTEX_NORMAL, PTHREAD_PROCESS_SHARED, CLOCK_REALTIME,
132		    "Normal+PShared"}
133	, {
134	PTHREAD_MUTEX_ERRORCHECK, PTHREAD_PROCESS_PRIVATE,
135		    CLOCK_REALTIME, "Errorcheck"}
136	, {
137	PTHREAD_MUTEX_ERRORCHECK, PTHREAD_PROCESS_SHARED,
138		    CLOCK_REALTIME, "Errorcheck+PShared"}
139	, {
140	PTHREAD_MUTEX_RECURSIVE, PTHREAD_PROCESS_PRIVATE,
141		    CLOCK_REALTIME, "Recursive"}
142	, {
143	PTHREAD_MUTEX_RECURSIVE, PTHREAD_PROCESS_SHARED, CLOCK_REALTIME,
144		    "Recursive+PShared"}
145#endif
146#ifdef USE_ALTCLK
147	, {
148	PTHREAD_MUTEX_DEFAULT, PTHREAD_PROCESS_PRIVATE, CLOCK_MONOTONIC,
149		    "Monotonic"}
150	, {
151	PTHREAD_MUTEX_DEFAULT, PTHREAD_PROCESS_SHARED, CLOCK_MONOTONIC,
152		    "PShared+Monotonic"}
153#ifndef WITHOUT_XOPEN
154	, {
155	PTHREAD_MUTEX_NORMAL, PTHREAD_PROCESS_PRIVATE, CLOCK_MONOTONIC,
156		    "Normal+Monotonic"}
157	, {
158	PTHREAD_MUTEX_NORMAL, PTHREAD_PROCESS_SHARED, CLOCK_MONOTONIC,
159		    "Normal+PShared+Monotonic"}
160	, {
161	PTHREAD_MUTEX_ERRORCHECK, PTHREAD_PROCESS_PRIVATE,
162		    CLOCK_MONOTONIC, "Errorcheck+Monotonic"}
163	, {
164	PTHREAD_MUTEX_ERRORCHECK, PTHREAD_PROCESS_SHARED,
165		    CLOCK_MONOTONIC, "Errorcheck+PShared+Monotonic"}
166	, {
167	PTHREAD_MUTEX_RECURSIVE, PTHREAD_PROCESS_PRIVATE,
168		    CLOCK_MONOTONIC, "Recursive+Monotonic"}
169	, {
170	PTHREAD_MUTEX_RECURSIVE, PTHREAD_PROCESS_SHARED,
171		    CLOCK_MONOTONIC, "Recursive+PShared+Monotonic"}
172#endif
173#endif
174};
175
176#define NSCENAR (sizeof(test_scenar) / sizeof(test_scenar[0]))
177
178pthread_attr_t ta;
179
180/* The next structure is used to save the tests measures */
181typedef struct __mes_t {
182	int nthreads;
183	long _data[NSCENAR];	/* As we store µsec values, a long type should be amply enough. */
184	struct __mes_t *next;
185} mes_t;
186
187/**** do_measure
188 * This function will do a timedwait on the cond cnd after locking mtx.
189 * Once the timedwait times out, it will read the clock cid then
190 * compute the difference and put it into ts.
191 * This function must be called once test is ready, as the timeout will be very short. */
192void do_measure(pthread_mutex_t * mtx, pthread_cond_t * cnd, clockid_t cid,
193		struct timespec *ts)
194{
195	int ret, rc;
196
197	struct timespec ts_cnd, ts_clk;
198
199	/* lock the mutex */
200	ret = pthread_mutex_lock(mtx);
201	if (ret != 0) {
202		UNRESOLVED(ret, "Unable to lock mutex");
203	}
204
205	/* Prepare the timeout parameter */
206	ret = clock_gettime(cid, &ts_cnd);
207	if (ret != 0) {
208		UNRESOLVED(errno, "Unable to read clock");
209	}
210
211	ts_cnd.tv_nsec += MES_TIMEOUT;
212	while (ts_cnd.tv_nsec >= 1000000000) {
213		ts_cnd.tv_nsec -= 1000000000;
214		ts_cnd.tv_sec++;
215	}
216
217	/* Do the timedwait */
218	do {
219		rc = pthread_cond_timedwait(cnd, mtx, &ts_cnd);
220		/* Re-read the clock as soon as possible */
221		ret = clock_gettime(cid, &ts_clk);
222		if (ret != 0) {
223			UNRESOLVED(errno, "Unable to read clock");
224		}
225	}
226	while (rc == 0);
227	if (rc != ETIMEDOUT) {
228		UNRESOLVED(rc,
229			   "Timedwait returned an unexpected error (expected ETIMEDOUT)");
230	}
231
232	/* Compute the difference time */
233	ts->tv_sec = ts_clk.tv_sec - ts_cnd.tv_sec;
234	ts->tv_nsec = ts_clk.tv_nsec - ts_cnd.tv_nsec;
235	if (ts->tv_nsec < 0) {
236		ts->tv_nsec += 1000000000;
237		ts->tv_sec -= 1;
238	}
239
240	if (ts->tv_sec < 0) {
241		FAILED
242		    ("The function returned from wait with a timeout before the time has passed\n");
243	}
244
245	/* unlock the mutex */
246	ret = pthread_mutex_unlock(mtx);
247	if (ret != 0) {
248		UNRESOLVED(ret, "Unable to unlock mutex");
249	}
250
251	return;
252}
253
254void *waiter(void *arg)
255{
256	test_t *dt = (test_t *) arg;
257
258	int ret;
259
260	ret = pthread_mutex_lock(dt->mtx);
261	if (ret != 0) {
262		UNRESOLVED(ret, "Mutex lock failed in waiter");
263	}
264#ifdef USE_CANCEL
265	pthread_cleanup_push((void *)pthread_mutex_unlock, (void *)(dt->mtx));
266#endif
267
268	/* This thread is ready to wait */
269	*(dt->tnum) += 1;
270
271	do {
272		ret = pthread_cond_wait(dt->cnd, dt->mtx);
273	}
274	while ((ret == 0) && (*(dt->predicate) == 0));
275	if (ret != 0) {
276		UNRESOLVED(ret, "pthread_cond_wait failed in waiter");
277	}
278#ifdef USE_CANCEL
279	pthread_cleanup_pop(0);	/* We could put 1 and avoid the next line, but we would miss the return code */
280#endif
281
282	ret = pthread_mutex_unlock(dt->mtx);
283	if (ret != 0) {
284		UNRESOLVED(ret, "Mutex unlock failed in waiter");
285	}
286
287	return NULL;
288}
289
290/**** do_threads_test
291 * This function is responsible for all the testing with a given # of threads
292 *  nthreads is the amount of threads to create.
293 * the return value is:
294 *  < 0 if function was not able to create enough threads.
295 *  cumulated # of nanoseconds otherwise.
296 */
297long do_threads_test(int nthreads, mes_t * measure)
298{
299	int ret;
300
301	int scal, i, j;
302
303	pthread_t *th;
304
305	int s;
306	pthread_mutexattr_t ma;
307	pthread_condattr_t ca;
308
309	pthread_cond_t cnd;
310	pthread_mutex_t mtx;
311	int predicate;
312	int tnum;
313
314	test_t td;
315
316	struct timespec ts, ts_cumul;
317
318	td.mtx = &mtx;
319	td.cnd = &cnd;
320	td.predicate = &predicate;
321	td.tnum = &tnum;
322
323	/* Allocate space for the threads structures */
324	th = (pthread_t *) calloc(nthreads, sizeof(pthread_t));
325	if (th == NULL) {
326		UNRESOLVED(errno, "Not enough memory for thread storage");
327	}
328#ifdef PLOT_OUTPUT
329	output("%d", nthreads);
330#endif
331	/* For each test scenario (mutex and cond attributes) */
332	for (s = 0; s < NSCENAR; s++) {
333		/* Initialize the attributes */
334		ret = pthread_mutexattr_init(&ma);
335		if (ret != 0) {
336			UNRESOLVED(ret,
337				   "Unable to initialize mutex attribute object");
338		}
339
340		ret = pthread_condattr_init(&ca);
341		if (ret != 0) {
342			UNRESOLVED(ret,
343				   "Unable to initialize cond attribute object");
344		}
345
346		/* Set the attributes according to the scenar and the impl capabilities */
347		ret = pthread_mutexattr_settype(&ma, test_scenar[s].mutex_type);
348		if (ret != 0) {
349			UNRESOLVED(ret, "Unable to set mutex type");
350		}
351
352		if (pshared_ok > 0) {
353			ret =
354			    pthread_mutexattr_setpshared(&ma,
355							 test_scenar[s].
356							 pshared);
357			if (ret != 0) {
358				UNRESOLVED(ret,
359					   "Unable to set mutex process-shared");
360			}
361
362			ret =
363			    pthread_condattr_setpshared(&ca,
364							test_scenar[s].pshared);
365			if (ret != 0) {
366				UNRESOLVED(ret,
367					   "Unable to set cond process-shared");
368			}
369		}
370#ifdef USE_ALTCLK
371		if (altclk_ok > 0) {
372			ret =
373			    pthread_condattr_setclock(&ca, test_scenar[s].cid);
374			if (ret != 0) {
375				UNRESOLVED(ret, "Unable to set clock for cond");
376			}
377		}
378#endif
379
380		ts_cumul.tv_sec = 0;
381		ts_cumul.tv_nsec = 0;
382
383#if VERBOSE > 1
384		output("Starting case %s\n", test_scenar[s].desc);
385#endif
386
387		for (scal = 0; scal < 5 * SCALABILITY_FACTOR; scal++) {
388			/* Initialize the mutex, the cond, and other data */
389			ret = pthread_mutex_init(&mtx, &ma);
390			if (ret != 0) {
391				UNRESOLVED(ret, "Mutex init failed");
392			}
393
394			ret = pthread_cond_init(&cnd, &ca);
395			if (ret != 0) {
396				UNRESOLVED(ret, "Cond init failed");
397			}
398
399			predicate = 0;
400			tnum = 0;
401
402			/* Create the waiter threads */
403			for (i = 0; i < nthreads; i++) {
404				ret =
405				    pthread_create(&th[i], &ta, waiter,
406						   (void *)&td);
407				if (ret != 0) {	/* We reached the limits */
408#if VERBOSE > 1
409					output
410					    ("Limit reached with %i threads\n",
411					     i);
412#endif
413#ifdef USE_CANCEL
414					for (j = i - 1; j >= 0; j--) {
415						ret = pthread_cancel(th[j]);
416						if (ret != 0) {
417							UNRESOLVED(ret,
418								   "Unable to cancel a thread");
419						}
420					}
421#else
422					predicate = 1;
423					ret = pthread_cond_broadcast(&cnd);
424					if (ret != 0) {
425						UNRESOLVED(ret,
426							   "Unable to broadcast the condition");
427					}
428#endif
429					for (j = i - 1; j >= 0; j--) {
430						ret = pthread_join(th[j], NULL);
431						if (ret != 0) {
432							UNRESOLVED(ret,
433								   "Unable to join a canceled thread");
434						}
435					}
436					free(th);
437					return -1;
438				}
439			}
440			/* All waiter threads are created now */
441#if VERBOSE > 5
442			output("%i waiter threads created successfully\n", i);
443#endif
444
445			ret = pthread_mutex_lock(&mtx);
446			if (ret != 0) {
447				UNRESOLVED(ret, "Unable to lock mutex");
448			}
449
450			/* Wait for all threads be waiting */
451			while (tnum < nthreads) {
452				ret = pthread_mutex_unlock(&mtx);
453				if (ret != 0) {
454					UNRESOLVED(ret, "Mutex unlock failed");
455				}
456
457				sched_yield();
458
459				ret = pthread_mutex_lock(&mtx);
460				if (ret != 0) {
461					UNRESOLVED(ret, "Mutex lock failed");
462				}
463			}
464
465			/* All threads are now waiting - we do the measure */
466
467			ret = pthread_mutex_unlock(&mtx);
468			if (ret != 0) {
469				UNRESOLVED(ret, "Mutex unlock failed");
470			}
471#if VERBOSE > 5
472			output("%i waiter threads are waiting; start measure\n",
473			       tnum);
474#endif
475
476			do_measure(&mtx, &cnd, test_scenar[s].cid, &ts);
477
478#if VERBOSE > 5
479			output("Measure for %s returned %d.%09d\n",
480			       test_scenar[s].desc, ts.tv_sec, ts.tv_nsec);
481#endif
482
483			ts_cumul.tv_sec += ts.tv_sec;
484			ts_cumul.tv_nsec += ts.tv_nsec;
485			if (ts_cumul.tv_nsec >= 1000000000) {
486				ts_cumul.tv_nsec -= 1000000000;
487				ts_cumul.tv_sec += 1;
488			}
489
490			/* We can release the threads */
491			predicate = 1;
492			ret = pthread_cond_broadcast(&cnd);
493			if (ret != 0) {
494				UNRESOLVED(ret,
495					   "Unable to broadcast the condition");
496			}
497#if VERBOSE > 5
498			output("Joining the waiters...\n");
499#endif
500
501			/* We will join every threads */
502			for (i = 0; i < nthreads; i++) {
503				ret = pthread_join(th[i], NULL);
504				if (ret != 0) {
505					UNRESOLVED(ret,
506						   "Unable to join a thread");
507				}
508
509			}
510
511			/* Destroy everything */
512			ret = pthread_mutex_destroy(&mtx);
513			if (ret != 0) {
514				UNRESOLVED(ret, "Unable to destroy mutex");
515			}
516
517			ret = pthread_cond_destroy(&cnd);
518			if (ret != 0) {
519				UNRESOLVED(ret, "Unable to destroy cond");
520			}
521		}
522
523#ifdef PLOT_OUTPUT
524		output(" %d.%09d", ts_cumul.tv_sec, ts_cumul.tv_nsec);
525#endif
526
527		measure->_data[s] = ts_cumul.tv_sec * 1000000 + (ts_cumul.tv_nsec / 1000);	/* We reduce precision */
528
529		/* Destroy the mutex attributes */
530		ret = pthread_mutexattr_destroy(&ma);
531		if (ret != 0) {
532			UNRESOLVED(ret, "Unable to destroy mutex attribute");
533		}
534
535		ret = pthread_condattr_destroy(&ca);
536		if (ret != 0) {
537			UNRESOLVED(ret, "Unable to destroy cond attribute");
538		}
539	}
540
541	/* Free the memory */
542	free(th);
543
544#if VERBOSE > 2
545	output("%5d threads; %d.%09d s (%i loops)\n", nthreads, ts_cumul.tv_sec,
546	       ts_cumul.tv_nsec, scal);
547#endif
548
549#ifdef PLOT_OUTPUT
550	output("\n");
551#endif
552
553	return ts_cumul.tv_sec * 1000000000 + ts_cumul.tv_nsec;
554}
555
556/* Forward declaration */
557int parse_measure(mes_t * measures);
558
559/* Main
560 */
561int main(int argc, char *argv[])
562{
563	int ret, nth;
564	long dur;
565
566	/* Initialize the measure list */
567	mes_t sentinel;
568	mes_t *m_cur, *m_tmp;
569	m_cur = &sentinel;
570	m_cur->next = NULL;
571
572	/* Initialize the output */
573	output_init();
574
575	/* Test machine capabilities */
576	/* -> clockid_t; pshared; ... */
577	altclk_ok = sysconf(_SC_CLOCK_SELECTION);
578	if (altclk_ok > 0)
579		altclk_ok = sysconf(_SC_MONOTONIC_CLOCK);
580#ifndef USE_ALTCLK
581	if (altclk_ok > 0)
582		output
583		    ("Implementation supports the MONOTONIC CLOCK but option is disabled in test.\n");
584#endif
585
586	pshared_ok = sysconf(_SC_THREAD_PROCESS_SHARED);
587
588#if VERBOSE > 0
589	output("Test starting\n");
590	output(" Process-shared primitive %s be tested\n",
591	       (pshared_ok > 0) ? "will" : "won't");
592	output(" Alternative clock for cond %s be tested\n",
593	       (altclk_ok > 0) ? "will" : "won't");
594#endif
595
596	/* Prepare thread attribute */
597	ret = pthread_attr_init(&ta);
598	if (ret != 0) {
599		UNRESOLVED(ret, "Unable to initialize thread attributes");
600	}
601
602	ret = pthread_attr_setstacksize(&ta, sysconf(_SC_THREAD_STACK_MIN));
603	if (ret != 0) {
604		UNRESOLVED(ret, "Unable to set stack size to minimum value");
605	}
606#ifdef PLOT_OUTPUT
607	output("# COLUMNS %d #threads", NSCENAR + 1);
608	for (nth = 0; nth < NSCENAR; nth++)
609		output(" %s", test_scenar[nth].desc);
610	output("\n");
611#endif
612
613	/* Do the testing */
614	nth = 0;
615	do {
616		nth += 100 * SCALABILITY_FACTOR;
617
618		/* Create a new measure item */
619		m_tmp = malloc(sizeof(mes_t));
620		if (m_tmp == NULL) {
621			UNRESOLVED(errno,
622				   "Unable to alloc memory for measure saving");
623		}
624		m_tmp->nthreads = nth;
625		m_tmp->next = NULL;
626
627		/* Run the test */
628		dur = do_threads_test(nth, m_tmp);
629
630		/* If test was success, add this measure to the list. Otherwise, free the mem */
631		if (dur >= 0) {
632			m_cur->next = m_tmp;
633			m_cur = m_tmp;
634		} else {
635			free(m_tmp);
636		}
637	}
638	while (dur >= 0);
639
640	/* We will now parse the results to determine if the measure is ~ constant or is growing. */
641
642	ret = parse_measure(&sentinel);
643
644	/* Free the memory from the list */
645	m_cur = sentinel.next;
646	while (m_cur != NULL) {
647		m_tmp = m_cur;
648		m_cur = m_tmp->next;
649		free(m_tmp);
650	}
651
652	if (ret != 0) {
653		FAILED("This function is not scalable");
654	}
655#if VERBOSE > 0
656	output("The function is scalable\n");
657#endif
658
659	PASSED;
660}
661
662/***
663 * The next function will seek for the better model for each series of measurements.
664 *
665 * The tested models are: -- X = # threads; Y = latency
666 * -> Y = a;      -- Error is r1 = avg((Y - Yavg)²);
667 * -> Y = aX + b; -- Error is r2 = avg((Y -aX -b)²);
668 *                -- where a = avg ((X - Xavg)(Y - Yavg)) / avg((X - Xavg)²)
669 *                --         Note: We will call _q = sum((X - Xavg) * (Y - Yavg));
670 *                --                       and  _d = sum((X - Xavg)²);
671 *                -- and   b = Yavg - a * Xavg
672 * -> Y = c * X^a;-- Same as previous, but with log(Y) = a log(X) + b; and b = log(c). Error is r3
673 * -> Y = exp(aX + b); -- log(Y) = aX + b. Error is r4
674 *
675 * We compute each error factor (r1, r2, r3, r4) then search which is the smallest (with ponderation).
676 * The function returns 0 when r1 is the best for all cases (latency is constant) and !0 otherwise.
677 */
678
679struct row {
680	long X;			/* the X values -- copied from function argument */
681	long Y[NSCENAR];	/* the Y values -- copied from function argument */
682	long _x;		/* Value X - Xavg */
683	long _y[NSCENAR];	/* Value Y - Yavg */
684	double LnX;		/* Natural logarithm of X values */
685	double LnY[NSCENAR];	/* Natural logarithm of Y values */
686	double _lnx;		/* Value LnX - LnXavg */
687	double _lny[NSCENAR];	/* Value LnY - LnYavg */
688};
689
690int parse_measure(mes_t * measures)
691{
692	int ret, i, r;
693
694	mes_t *cur;
695
696	double Xavg, Yavg[NSCENAR];
697	double LnXavg, LnYavg[NSCENAR];
698
699	int N;
700
701	double r1[NSCENAR], r2[NSCENAR], r3[NSCENAR], r4[NSCENAR];
702
703	/* Some more intermediate vars */
704	long double _q[3][NSCENAR];
705	long double _d[3][NSCENAR];
706
707	long double t;		/* temp value */
708
709	struct row *Table = NULL;
710
711	/* Initialize the datas */
712	Xavg = 0.0;
713	LnXavg = 0.0;
714	for (i = 0; i < NSCENAR; i++) {
715		Yavg[i] = 0.0;
716		LnYavg[i] = 0.0;
717		r1[i] = 0.0;
718		r2[i] = 0.0;
719		r3[i] = 0.0;
720		r4[i] = 0.0;
721		_q[0][i] = 0.0;
722		_q[1][i] = 0.0;
723		_q[2][i] = 0.0;
724		_d[0][i] = 0.0;
725		_d[1][i] = 0.0;
726		_d[2][i] = 0.0;
727	}
728	N = 0;
729	cur = measures;
730
731#if VERBOSE > 1
732	output("Data analysis starting\n");
733#endif
734
735	/* We start with reading the list to find:
736	 * -> number of elements, to assign an array
737	 * -> average values
738	 */
739	while (cur->next != NULL) {
740		cur = cur->next;
741
742		N++;
743
744		Xavg += (double)cur->nthreads;
745		LnXavg += log((double)cur->nthreads);
746
747		for (i = 0; i < NSCENAR; i++) {
748			Yavg[i] += (double)cur->_data[i];
749			LnYavg[i] += log((double)cur->_data[i]);
750		}
751	}
752
753	/* We have the sum; we can divide to obtain the average values */
754	Xavg /= N;
755	LnXavg /= N;
756
757	for (i = 0; i < NSCENAR; i++) {
758		Yavg[i] /= N;
759		LnYavg[i] /= N;
760	}
761
762#if VERBOSE > 1
763	output(" Found %d rows and %d columns\n", N, NSCENAR);
764#endif
765
766	/* We will now alloc the array ... */
767	Table = calloc(N, sizeof(struct row));
768	if (Table == NULL) {
769		UNRESOLVED(errno, "Unable to alloc space for results parsing");
770	}
771
772	/* ... and fill it */
773	N = 0;
774	cur = measures;
775
776	while (cur->next != NULL) {
777		cur = cur->next;
778
779		Table[N].X = (long)cur->nthreads;
780		Table[N]._x = Table[N].X - Xavg;
781		Table[N].LnX = log((double)cur->nthreads);
782		Table[N]._lnx = Table[N].LnX - LnXavg;
783		for (i = 0; i < NSCENAR; i++) {
784			Table[N].Y[i] = cur->_data[i];
785			Table[N]._y[i] = Table[N].Y[i] - Yavg[i];
786			Table[N].LnY[i] = log((double)cur->_data[i]);
787			Table[N]._lny[i] = Table[N].LnY[i] - LnYavg[i];
788		}
789
790		N++;
791	}
792
793	/* We won't need the list anymore -- we'll work with the array which should be faster. */
794#if VERBOSE > 1
795	output(" Data was stored in an array.\n");
796#endif
797
798	/* We need to read the full array at least twice to compute all the error factors */
799
800	/* In the first pass, we'll compute:
801	 * -> r1 for each scenar.
802	 * -> "a" factor for linear (0), power (1) and exponential (2) approximations -- with using the _d and _q vars.
803	 */
804#if VERBOSE > 1
805	output("Starting first pass...\n");
806#endif
807	for (r = 0; r < N; r++) {
808		for (i = 0; i < NSCENAR; i++) {
809			r1[i] +=
810			    ((double)Table[r]._y[i] / N) *
811			    (double)Table[r]._y[i];
812
813			_q[0][i] += Table[r]._y[i] * Table[r]._x;
814			_d[0][i] += Table[r]._x * Table[r]._x;
815
816			_q[1][i] += Table[r]._lny[i] * Table[r]._lnx;
817			_d[1][i] += Table[r]._lnx * Table[r]._lnx;
818
819			_q[2][i] += Table[r]._lny[i] * Table[r]._x;
820			_d[2][i] += Table[r]._x * Table[r]._x;
821		}
822	}
823
824	/* First pass is terminated; a2 = _q[0]/_d[0]; a3 = _q[1]/_d[1]; a4 = _q[2]/_d[2] */
825
826	/* In the first pass, we'll compute:
827	 * -> r2, r3, r4 for each scenar.
828	 */
829
830#if VERBOSE > 1
831	output("Starting second pass...\n");
832#endif
833	for (r = 0; r < N; r++) {
834		for (i = 0; i < NSCENAR; i++) {
835			/* r2 = avg((y - ax -b)²);  t = (y - ax - b) = (y - yavg) - a (x - xavg); */
836			t = (Table[r]._y[i] -
837			     ((_q[0][i] * Table[r]._x) / _d[0][i]));
838			r2[i] += t * t / N;
839
840			/* r3 = avg((y - c.x^a) ²);
841			   t = y - c * x ^ a
842			   = y - log (LnYavg - (_q[1]/_d[1]) * LnXavg) * x ^ (_q[1]/_d[1])
843			 */
844			t = (Table[r].Y[i]
845			     - (logl(LnYavg[i] - (_q[1][i] / _d[1][i]) * LnXavg)
846				* powl(Table[r].X, (_q[1][i] / _d[1][i]))
847			     ));
848			r3[i] += t * t / N;
849
850			/* r4 = avg((y - exp(ax+b))²);
851			   t = y - exp(ax+b)
852			   = y - exp(_q[2]/_d[2] * x + (LnYavg - (_q[2]/_d[2] * Xavg)));
853			   = y - exp(_q[2]/_d[2] * (x - Xavg) + LnYavg);
854			 */
855			t = (Table[r].Y[i]
856			     - expl((_q[2][i] / _d[2][i]) * Table[r]._x +
857				    LnYavg[i]));
858			r4[i] += t * t / N;
859
860		}
861	}
862
863#if VERBOSE > 1
864	output("All computing terminated.\n");
865#endif
866	ret = 0;
867	for (i = 0; i < NSCENAR; i++) {
868#if VERBOSE > 1
869		output("\nScenario: %s\n", test_scenar[i].desc);
870
871		output("  Model: Y = k\n");
872		output("       k = %g\n", Yavg[i]);
873		output("    Divergence %g\n", r1[i]);
874
875		output("  Model: Y = a * X + b\n");
876		output("       a = %Lg\n", _q[0][i] / _d[0][i]);
877		output("       b = %Lg\n",
878		       Yavg[i] - ((_q[0][i] / _d[0][i]) * Xavg));
879		output("    Divergence %g\n", r2[i]);
880
881		output("  Model: Y = c * X ^ a\n");
882		output("       a = %Lg\n", _q[1][i] / _d[1][i]);
883		output("       c = %Lg\n",
884		       logl(LnYavg[i] - (_q[1][i] / _d[1][i]) * LnXavg));
885		output("    Divergence %g\n", r2[i]);
886
887		output("  Model: Y = exp(a * X + b)\n");
888		output("       a = %Lg\n", _q[2][i] / _d[2][i]);
889		output("       b = %Lg\n",
890		       LnYavg[i] - ((_q[2][i] / _d[2][i]) * Xavg));
891		output("    Divergence %g\n", r2[i]);
892#endif
893		/* Compare r1 to other values, with some ponderations */
894		if ((r1[i] > 1.1 * r2[i]) || (r1[i] > 1.2 * r3[i])
895		    || (r1[i] > 1.3 * r4[i]))
896			ret++;
897#if VERBOSE > 1
898		else
899			output(" Sanction: OK\n");
900#endif
901	}
902
903	/* We need to free the array */
904	free(Table);
905
906	/* We're done */
907	return ret;
908}
909