/*

Copyright (C) 2008-2021 Michele Martone

This file is part of librsb.

librsb is free software; you can redistribute it and/or modify it
under the terms of the GNU Lesser General Public License as published
by the Free Software Foundation; either version 3 of the License, or
(at your option) any later version.

librsb 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 Lesser General Public
License for more details.

You should have received a copy of the GNU Lesser General Public
License along with librsb; see the file COPYING.
If not, see <http://www.gnu.org/licenses/>.

*/
/* @cond INNERDOC  */
/*!
 * @file
 * @author Michele Martone
 * @brief Auto tuning mechanisms.
 * */

#include <unistd.h>	/* sysconf */
#include "rsb_internals.h"
#include "rsb_lock.h"
RSB_INTERNALS_COMMON_HEAD_DECLS

#define RSB_TUNE_WITH_LIKWID 0 /* TODO: for future reference */

#if RSB_TUNE_WITH_LIKWID
#include <likwid.h>
#define RSB_LIKWID_MARKER_INIT	likwid_markerInit()
#define RSB_LIKWID_MARKER_EXIT likwid_markerClose()
#define RSB_LIKWID_MARKER_R_START(R) likwid_markerStartRegion(R)
#define RSB_LIKWID_MARKER_R_STOP(R) likwid_markerStopRegion(R)
#else /* RSB_TUNE_WITH_LIKWID */
#define RSB_LIKWID_MARKER_INIT
#define RSB_LIKWID_MARKER_EXIT
#define RSB_LIKWID_MARKER_R_START(R)
#define RSB_LIKWID_MARKER_R_STOP(R)
#endif /* RSB_TUNE_WITH_LIKWID */

#if RSB_TUNE_WITH_LIKWID
#define RSB_TM_LIKWID_MARKER_R_START(R) if(want_likwid == RSB_BOOL_TRUE)RSB_LIKWID_MARKER_R_START(R)
#define RSB_TM_LIKWID_MARKER_R_STOP(R)  if(want_likwid == RSB_BOOL_TRUE)RSB_LIKWID_MARKER_R_STOP(R)
#else
#define RSB_TM_LIKWID_MARKER_R_START(R)
#define RSB_TM_LIKWID_MARKER_R_STOP(R)
#endif /* RSB_TUNE_WITH_LIKWID */
#if RSB_TUNE_WITH_LIKWID
//#define want_likwid RSB_BOOL_TRUE
#endif /* RSB_TUNE_WITH_LIKWID */

#define rsb__calloca_vector rsb__calloc_vector
#define RSB_CONDITIONAL_FREE_ALLOCA RSB_CONDITIONAL_FREE

#if 0
#define RSB_THREAD_STATS RSB_PRINT_THREAD_STATS
#else
#define RSB_THREAD_STATS
#endif

#define RSB_AT_NO_VERBOSE(VERBOSE) ((VERBOSE)<1)
#define RSB_VERBOSE_FOR_MSG(VERBOSE) ((VERBOSE)>2)
#define RSB_VERBOSE_FOR_SAVING(VERBOSE) ((VERBOSE)>1)

#define RSB_ATR_SAVENAME(FN,ERRVAL,RFPFX,BASEFILENAME,TYPECODE,TRANSA,NRHS,SVAL,ROP,RVAL,TVAL,MTX,IOT,W_SPMM)	 { \
       	char ename[RSB_MAX_FILENAME_LENGTH]; \
	if(TVAL>0) \
	{	\
	       	sprintf(ename,"--%s-tuned_%s%d_%0lgx%0dth.eps",(W_SPMM?"mv":"sv"),ROP,RVAL,SVAL,TVAL); \
       	} \
       	else \
       	{ \
		sprintf(ename,"--base.eps"); \
	} \
	sprintf(FN,"%s%s--%c-%c-%ld%s",RFPFX,BASEFILENAME,TYPECODE,RSB_TRANSPOSITION_AS_CHAR(TRANSA),(long int)(NRHS),ename); \
	}

#if 0
#define RSB_ATR_SAVE(ERRVAL,RFPFX,BASEFILENAME,TYPECODE,TRANSA,NRHS,SVAL,ROP,TVAL,MTX,IOT,W_SPMM) { RSB_NULL_STATEMENT_FOR_COMPILER_HAPPINESS; }
#else
#define RSB_ATR_SAVE(ERRVAL,RFPFX,BASEFILENAME,TYPECODE,TRANSA,NRHS,SVAL,ROP,RVAL,TVAL,MTX,IOT,W_SPMM) { \
	char fname[RSB_MAX_FILENAME_LENGTH+32]; \
	(IOT) -= rsb_time(); \
	RSB_ATR_SAVENAME(fname,ERRVAL,RFPFX,BASEFILENAME,TYPECODE,TRANSA,NRHS,SVAL,ROP,RVAL,TVAL,MTX,IOT,W_SPMM)	  \
       	(ERRVAL) = rsb_mtx_rndr(fname,MTX,RSB_DEFAULT_MATRIX_RENDERING_ROWS,RSB_DEFAULT_MATRIX_RENDERING_COLS,/*RSB_MARF_EPS_B*/RSB_MARF_EPS_L); \
       	(IOT) += rsb_time(); \
	RSB_STDOUT("Saved plot to  %s\n",fname); \
}
#endif

#define RSB_AP "" /* FIXME */
#define RSB_INVALID_SVAL -1.0 /* */
#define RSB_PERF_ZERO RSB_REAL_ZERO /* TODO: to rsb_tune.h */

#define RSB_AT_SAMPLE_MINDT 0.0

#if RSB_WANT_MKL
static double rsb__estimate_mflops_per_op_spmv(size_t Ec, rsb_type_t typecode, rsb_flags_t flags)
/* generated by the RSB_M4_ESTIMATE_MFLOPS_PER_MOP_FUNCTION macro */
{
	/*!
	 * \ingroup gr_internals
	 * A function which returns the approximate count of floating point operations
	 * needed to perform the "spmv_uaua" matrix operation.
	 * In the case of symmetric/hermitian, the number of operations is multiplied by two.
	 * In the case of a complex type, the number of operations is multiplied by six:
	 * (a +bi)*(c+di) = (ac-bd)+(ad+bc)i
	 * accounts for four real multiplications and two real additions.
	 * FIXME: complexity is NOT taken in consideration for non-SPMV/SPSV operations
	 */

	const double M_  = 1000000.0;
	double Me = Ec;

	if(RSB_IS_MATRIX_TYPE_COMPLEX(typecode)) { Me=8*Ec; } else { Me=2*Ec; }
	if(( flags & RSB_FLAG_ANY_SYMMETRY )){ Me*=2; }/* slightly optimistic : should subtract the diagonal elements count */
	Me /= M_;
	return Me;
}
#endif /* RSB_WANT_MKL */

void rsb__tattr_init(struct rsb_tattr_t* TTRP, const struct rsb_mtx_t*MTXAP, rsb_coo_idx_t nrA, rsb_nnz_idx_t nnz, rsb_type_t typecode, rsb_flags_t flags, rsb_int_t nrhs)
{
	if(!TTRP)
		goto err;
#if !RSB_WANT_MKL
	if(!MTXAP)
		goto err;
#endif /* RSB_WANT_MKL */

	RSB_BZERO_P(TTRP);
	(TTRP)->mint = RSB_TRACE_MAX_THREADS_P1;
	(TTRP)->maxt = 0;
	if(MTXAP)
	{
		(TTRP)->mtxAc = *(MTXAP); /*(TTRP)->mtxAp = (MTXAP);*/
		rsb__init_blank_pointers(&(TTRP)->mtxAc);
       		(TTRP)->bpn = ((rsb_real_t)rsb__get_index_storage_amount(MTXAP))/((MTXAP)->nnz);
       		(TTRP)->ofe = rsb__estimate_mflops_per_op_spmv_uaua(MTXAP);
	}
#if RSB_WANT_MKL
	else
	{
		/* used by rsbench for int-indexed CSR */
		size_t idx = sizeof(int);
		idx *= nrA + nnz ;
       		(TTRP)->bpn = idx;
       		(TTRP)->bpn /= nnz;
       		(TTRP)->ofe = rsb__estimate_mflops_per_op_spmv(nnz, typecode, flags);
	}
#endif /* RSB_WANT_MKL */
       	(TTRP)->ttt = -rsb_time(); /* FIXME: if no rsb__tattr_finalize, whi shall update this ? */
       	(TTRP)->ofe *= nrhs;
	(TTRP)->btpo = RSB_CONST_IMPOSSIBLY_BIG_TIME;
	(TTRP)->dtpo = RSB_TIME_ZERO;
	(TTRP)->vl = 0;
err:
	return;
}

void rsb__tattr_sets(struct rsb_tattr_t* ttrp, rsb_int_t dnt, rsb_int_t nt, rsb_time_t tpo, rsb_int_t bnt, rsb_int_t nits)
{
	/* TODO: verbose to a stream */

	if(!ttrp)
		goto ret;

	if(nt == 0)
	{
		nt = dnt = bnt = rsb_get_num_threads();
	}

	if(nt>0 && nt < RSB_TRACE_MAX_THREADS_P1)
	{
       		(ttrp)->tpo[(nt)] = (tpo),
       		(ttrp)->nit[(nt)] = (nits);
	}
	else
	{
		RSB_ERROR("Threads specification (%d) is wrong!",nt);
		RSB_ASSERT(0);/* TODO: write error case here */
	}
	(ttrp)->mint = RSB_MIN((ttrp)->mint,( nt));
	(ttrp)->maxt = RSB_MAX((ttrp)->maxt,( nt));
	(ttrp)->optt = (bnt);
	(ttrp)->deft = (dnt);
	(ttrp)->btpo = RSB_MIN((ttrp)->btpo,tpo);
	if((dnt)==(nt))
	       	(ttrp)->dtpo = (tpo);

	if(ttrp->vl>0)
		RSB_STDOUT("# for %d threads, %d best threads, %d default threads, %0.4lgs p.op, %0.4lg Mflops\n",nt,bnt,dnt,tpo,(ttrp->ofe)/tpo);
ret:
	return;
}

// #define RSB_ITS "%d" /* operations */
#define RSB_ITS "%6.3lf" /* bytes per non zero */
#define RSB_BPS "%6.3lf" /* bytes per non zero */
#define RSB_PPS "%6.2le" /* performance (mflops) */
#define RSB_SPS "%6.3lf" /* speedup */
#define RSB_CPS "%6.1lf" /* percentage */
#define RSB_TPS "%8.2le" /* time */

#define RSB_TTRP_OK(TTRP) ( (TTRP)->mint <= (TTRP)->maxt && (TTRP)->maxt > 0 )

static void rsb__tattr_dump(FILE*FP, struct rsb_tattr_t* TTRP)
{
	rsb_int_t curt; /* FIXME: this shall go to a pair of .dat/.gp files */ /* FIXME: bpnz printout ? */

	if(!TTRP)
		goto ret;

       	if( RSB_TTRP_OK (TTRP) )
	{
		RSB_FPRINTF(FP,"%s",RSB_AP);
		if(!RSB_MATRIX_UNSUPPORTED_TYPE(TTRP->mtxAc.typecode))
			RSB_FPRINTF(FP,"# " RSB_PRINTF_MATRIX_AT_SUMMARY_ARGS( &((TTRP)->mtxAc))), RSB_FPRINTF(FP,"\n");
		RSB_FPRINTF(FP,RSB_AP"# [mint:%d ... optt:%d ... maxt:%d], btpo:"RSB_PPS" bpn:"RSB_BPS" dtn:%d\n",(TTRP)->mint,(TTRP)->optt,(TTRP)->maxt,(TTRP)->btpo,(TTRP)->bpn,(TTRP)->deft);
		RSB_FPRINTF(FP,RSB_AP"# threads perf bpn tpo subm times\n");
		for(curt=(TTRP)->mint;curt<=(TTRP)->maxt;curt++)
		{
			RSB_FPRINTF(FP,RSB_AP" %d "RSB_TPS" "RSB_BPS" "RSB_TPS" %d %d\n",curt,(TTRP)->ofe/(TTRP)->tpo[curt],(TTRP)->bpn,(TTRP)->tpo[curt],
					RSB_MAX(1,(TTRP)->mtxAc.all_leaf_matrices_n),(TTRP)->nit[curt]
					);
		}
		RSB_FPRINTF(FP,RSB_AP"\n"RSB_AP"\n");
	}
ret:
	return;
}

static void rsb__attr_init(struct rsb_attr_t* TTRP, rsb_bool_t DT, const struct rsb_mtx_t*MTXAP, const char*BASEFILENAME, rsb_trans_t TRANSA, rsb_int_t NRHS, enum rsb_op_t op)
{
	if(!MTXAP || ! TTRP)
		goto ret;

	RSB_BZERO_P(TTRP);
	(TTRP)->lfp = stdout;
	(TTRP)->pfp = stdout;
	(TTRP)->dtr = (DT);
	(TTRP)->nrhs = (NRHS);
	(TTRP)->transA = (TRANSA);

	if((TTRP)->dtr)
	{
       		/*char ename[RSB_MAX_FILENAME_LENGTH];*/
       		char fname[RSB_MAX_FILENAME_LENGTH];
		char * ocs = op == rsb_op_spmv ? "mv" : "sv";

		(TTRP)->typecode = (MTXAP)->typecode;
		sprintf(fname,"%s--%c-%c-%d--%s%s",BASEFILENAME,(MTXAP)->typecode,RSB_TRANSPOSITION_AS_CHAR(TRANSA),NRHS,ocs,/*ename*/"-tuning_trace.dat");
		(TTRP)->lfp = rsb__util_fopen(fname,"w");
		sprintf(fname,"%s--%c-%c-%d--%s%s",BASEFILENAME,(MTXAP)->typecode,RSB_TRANSPOSITION_AS_CHAR(TRANSA),NRHS,ocs,/*ename*/"-tuning_trace.gnu");
		sprintf((TTRP)->bname,"%s--%c-%c-%d--%s%s",BASEFILENAME,(MTXAP)->typecode,RSB_TRANSPOSITION_AS_CHAR(TRANSA),NRHS,ocs,/*ename*/"-tuning_trace");
		(TTRP)->pfp = rsb__util_fopen(fname,"w");
		sprintf((TTRP)->mname,"%s",BASEFILENAME);
	}
ret:
	return;
}

void rsb__attr_dump(struct rsb_attr_t*TTRP)
{
	if(TTRP && TTRP->dtr)
	{
	       	rsb_int_t tri; /* FIME: this shall go to a pair of .dat/.gp files */
	       	rsb_int_t hmkl = ( RSB_TTRP_OK ( &((TTRP)->clattr)) );
		rsb_int_t smaxt = rsb__set_num_threads(/*RSB_THREADS_GET_MAX*/RSB_THREADS_GET_MAX_SYS),msmaxth = smaxt;
		rsb_time_t bptpo = RSB_CONST_IMPOSSIBLY_BIG_TIME, fdtpo = RSB_CONST_IMPOSSIBLY_BIG_TIME;
		rsb_perf_t ofe = RSB_PERF_ZERO;
		rsb_int_t br = 0;
		rsb_time_t ttt = RSB_TIME_ZERO, stt = RSB_TIME_ZERO;
		struct rsb_mtx_t mtxAc = (TTRP)->tattra[0].mtxAc ;
      	 	RSB_FPRINTF((TTRP)->lfp,RSB_AP"# BEGIN TRACE RECORD [%d entries]\n",(TTRP)->trc);

		rsb__tattr_dump((TTRP)->lfp, (& ((TTRP)->clattr)) );

       		for(tri=0;tri<(TTRP)->trc;tri++) {
       			RSB_FPRINTF((TTRP)->lfp,RSB_AP"# RECORD %d:\n",tri+1);
			rsb__tattr_dump((TTRP)->lfp, (&((TTRP)->tattra[tri])) );
			if( (TTRP)->tattra[ br ].btpo > (TTRP)->tattra[ tri ].btpo ) br = tri;
		}
		((TTRP)->br) = br;
       		RSB_FPRINTF((TTRP)->lfp,RSB_AP"# best path (threads perf bpn tpo subm spdp spdppcnt): \n");

       		for(tri=0;tri<(TTRP)->trc;tri++)
		{
       			if( tri == 0 )
				fdtpo = (TTRP)->tattra[0].dtpo;
			if( bptpo > (TTRP)->tattra[tri].btpo )
			{
	       			ofe = (TTRP)->tattra[tri].ofe;
				bptpo = (TTRP)->tattra[tri].btpo;
				RSB_FPRINTF((TTRP)->lfp,RSB_AP" %d "RSB_TPS" "RSB_BPS" "RSB_TPS" %d "RSB_SPS" "RSB_CPS,(TTRP)->tattra[tri].optt,ofe/(TTRP)->tattra[tri].btpo,(TTRP)->tattra[tri].bpn,(TTRP)->tattra[tri].btpo,(TTRP)->tattra[tri].mtxAc.all_leaf_matrices_n,fdtpo/(TTRP)->tattra[tri].btpo, RSB_SPEEDUP_TO_PCT(fdtpo/(TTRP)->tattra[tri].btpo)
						);
				/* FIXME: may put here speedup over MKL */
				RSB_FPRINTF((TTRP)->lfp,"\n");
			}
		}
		ttt = ((TTRP)->tattra[0 ].ttt);
       		for(tri=1;tri<(TTRP)->trc;tri++)
			stt += ((TTRP)->tattra[tri].ttt);
		RSB_FPRINTF((TTRP)->lfp,RSB_AP"\n"RSB_AP"\n");
       		RSB_FPRINTF((TTRP)->lfp,RSB_AP"# SUMMARY AUTOTUNING ttt="RSB_TPS" stt="RSB_TPS" \n",ttt,stt);
       		RSB_FPRINTF((TTRP)->lfp,RSB_AP"# SUMMARY: mtx ");
       		RSB_FPRINTF((TTRP)->lfp," sym transa typecode nr nc nnz nrhs");
       		if( hmkl )
       			RSB_FPRINTF((TTRP)->lfp," mkl-d mkl-t rsb-vs-mkl-d rsb-vs-mkl-t rsb-vs-mkl-s"),
       			RSB_FPRINTF((TTRP)->lfp," k-mkl-rsb-t k-mkl-vs-rsb-s k-rsb-rsb-t k-rsb-vs-rsb-s");
       		RSB_FPRINTF((TTRP)->lfp," rsb-d rsb-t rsb-s rsb-vs-rsb-t rsb-vs-rsb-s");
       		RSB_FPRINTF((TTRP)->lfp," bpn-d bpn-s");
       		RSB_FPRINTF((TTRP)->lfp," subm-d subm-s");
       		RSB_FPRINTF((TTRP)->lfp," speedup-d speedup-s");
       		RSB_FPRINTF((TTRP)->lfp,"\n");
       		RSB_FPRINTF((TTRP)->lfp,RSB_AP"# SUMMARY: %s ",(TTRP)->mname);
       		RSB_FPRINTF((TTRP)->lfp," %c %c %c %ld %ld %ld %ld",
				rsb__do_get_symmetry_char(&mtxAc), RSB_TRANSPOSITION_AS_CHAR((TTRP)->transA), mtxAc.typecode,
				(long int)mtxAc.nr,
				(long int)mtxAc.nc,
				(long int)mtxAc.nnz,
				(long int)(TTRP)->nrhs
				);
       		if( hmkl )
       			RSB_FPRINTF((TTRP)->lfp," "RSB_PPS " "RSB_PPS" "RSB_SPS" "RSB_SPS" "RSB_SPS" "
				,ofe/(TTRP)->clattr.dtpo, ofe/(TTRP)->clattr.btpo
				,(TTRP)->clattr.dtpo/(TTRP)->tattra[0 ].dtpo
				,(TTRP)->clattr.btpo/(TTRP)->tattra[0 ].btpo
				,(TTRP)->clattr.btpo/(TTRP)->tattra[br].btpo
					),
       			RSB_FPRINTF((TTRP)->lfp," "RSB_ITS" "RSB_ITS" "RSB_ITS" "RSB_ITS" "
				,( ttt / ((TTRP)->clattr.dtpo-(TTRP)->tattra[0 ].btpo) )
				,( stt / ((TTRP)->clattr.dtpo-(TTRP)->tattra[br].btpo) )
				,( ttt / ((TTRP)->tattra[0 ].dtpo-(TTRP)->tattra[0 ].btpo) )
				,( stt / ((TTRP)->tattra[0 ].dtpo-(TTRP)->tattra[br].btpo) )
				//,(TTRP)->clattr.dtpo-(TTRP)->tattra[0 ].dtpo
				   );
       		RSB_FPRINTF((TTRP)->lfp," "RSB_PPS" "RSB_PPS" "RSB_PPS" "RSB_SPS" "RSB_SPS,
				ofe/(TTRP)->tattra[0].dtpo,ofe/(TTRP)->tattra[0].btpo,ofe/(TTRP)->tattra[(TTRP)->br].btpo
				,(TTRP)->tattra[0].dtpo/(TTRP)->tattra[0 ].btpo
				,(TTRP)->tattra[0].dtpo/(TTRP)->tattra[br].btpo
				);
       		RSB_FPRINTF((TTRP)->lfp," "RSB_BPS" "RSB_BPS,(TTRP)->tattra[0].bpn,(TTRP)->tattra[br].bpn);
       		RSB_FPRINTF((TTRP)->lfp," %d %d",(TTRP)->tattra[0 ].mtxAc.all_leaf_matrices_n,(TTRP)->tattra[br].mtxAc.all_leaf_matrices_n);
       		RSB_FPRINTF((TTRP)->lfp," "RSB_BPS" "RSB_BPS,
				(TTRP)->tattra[0 ].tpo[1]/(TTRP)->tattra[0 ].tpo[(TTRP)->tattra[0 ].optt],
				(TTRP)->tattra[br].tpo[1]/(TTRP)->tattra[br].tpo[(TTRP)->tattra[br].optt]
				);
       		RSB_FPRINTF((TTRP)->lfp,"\n");
       		RSB_FPRINTF((TTRP)->lfp,RSB_AP"\n"RSB_AP"# END TRACE RECORD\n");

  	 	RSB_FPRINTF((TTRP)->pfp,
"set terminal postscript eps color\n\
set output '%s.eps'\n\
set xlabel 'threads'\n\
set ylabel 'Mflops'\n\
set title  'Trace of empirical auto-tuning of %s in the RSB format'\n\
#set title  'Trace of empirical auto-tuning of matrix %s in the RSB format'\n\
set xrange [ 1 : %d ] noreverse nowriteback\n\
plot ",(TTRP)->bname,(TTRP)->bname,(TTRP)->mname,smaxt);

 if(hmkl)
 	RSB_FPRINTF((TTRP)->pfp,"'%s.dat' index 0   u 1:2    with linespoints lw 2 lc rgb 'blue'   title 'MKL', ", (TTRP)->bname),
 	RSB_FPRINTF((TTRP)->pfp,"'%s.dat' index 0   u 1:2:((%lg*$3)) with circles  lc rgb 'blue'  fs border 0.15 notitle, ",(TTRP)->bname,0.2); // 0.2 -> 0.1
 if((TTRP)->trc > 0)
  RSB_FPRINTF((TTRP)->pfp,"'%s.dat' index %d   u 1:2    with linespoints lw 2 lc rgb 'red' title 'RSB-t', ", (TTRP)->bname,hmkl?1:0),
  RSB_FPRINTF((TTRP)->pfp,"'%s.dat' index %d   u 1:2:((%lg*$3)) with circles  lc rgb 'red'  fs border 0.15 notitle, ",(TTRP)->bname,hmkl?1:0,0.2); // 0.2 -> 0.1
 if((TTRP)->trc > 1)
  RSB_FPRINTF((TTRP)->pfp,"'%s.dat' index %d:%d:1 u 1:2 with linespoints lw 2 lc    rgb 'green' title 'RSB-s', ",
		  (TTRP)->bname, hmkl?2:1, (TTRP)->trc-(hmkl?0:1) ),
  RSB_FPRINTF((TTRP)->pfp,"'%s.dat' index %d:%d:1 u 1:2:((%lg*$3)) with circles  lc rgb 'green' fs border 0.15 notitle, ",
		  (TTRP)->bname, hmkl?2:1, (TTRP)->trc-(hmkl?0:1), 0.2 ),
  RSB_FPRINTF((TTRP)->pfp,"'%s.dat' index %d u 1:2 with linespoints lw 2 lc    rgb 'black' title 'top RSB', ",
		  (TTRP)->bname, (TTRP)->trc+(hmkl?1:0) );
 // '%s.dat' index 2   u 1:2:3:4 with vectors lw    2 lc rgb 'pink'  title 'pink', 
					//, bname
  for(tri=0;tri<(TTRP)->trc;tri++) {
   msmaxth = RSB_MIN(msmaxth,((TTRP)->tattra[tri].maxt-(TTRP)->tattra[tri].mint)/2);
  }
  RSB_FPRINTF((TTRP)->pfp,"'%s.dat' index %d:%d using 1:2:(sprintf(\"%%d subm\",$5)) every %d::%d with labels left offset 1 notitle",
		  (TTRP)->bname,0,(TTRP)->trc+(hmkl?0:-1), RSB_TRACE_MAX_THREADS_P1, msmaxth);
  RSB_FPRINTF((TTRP)->pfp,",'%s.dat' index %d using 1:2:(sprintf(\"%%6.3f x\",$6)) "/*every %d::%d*/" with labels left offset 0 notitle",
		  (TTRP)->bname,  (TTRP)->trc+(hmkl?1: 0)/*, RSB_TRACE_MAX_THREADS_P1, 1*/);

		if((TTRP)->dtr)
			fclose((TTRP)->lfp),
			fclose((TTRP)->pfp);
	}
}

#define RSB_ATTR_PREPARE(TTRP,ATTRP,ERRVAL) { if(ATTRP && ATTRP->dtr){  (TTRP)=&((ATTRP)->tattra[(ATTRP)->trc]); /*(TTRP)->ttt=-rsb_time();*//* (ATTRP)->trc++;*/ if ((ATTRP)->trc >= RSB_CONST_MAX_TUNING_SAMPLES )ERRVAL=RSB_ERR_INTERNAL_ERROR;} else {(TTRP)=NULL;} }
#define RSB_ATTR_ADVANCE(TTRP,ATTRP,ET) { if(ATTRP && ATTRP->dtr){/*(TTRP)=&((ATTRP)->tattra[(ATTRP)->trc]);*/ (TTRP)->ttt += rsb_time(); (TTRP)->ttt += (ET); (ATTRP)->trc++;  } else {(TTRP)=NULL;} /* printf("ADVANCING TO %d/%d  (%zd bytes)\n",(ATTRP)->trc,RSB_CONST_MAX_TUNING_SAMPLES,sizeof(struct rsb_attr_t));*/}

rsb_err_t rsb__do_bench_spxm(rsb_time_t *tpop, rsb_int_t *timesp, rsb_trans_t transA, const void * alphap, const struct rsb_mtx_t * mtxAp, rsb_coo_idx_t nrhs, rsb_flags_t order, const void * Bp, rsb_nnz_idx_t ldB, const void * betap, void * Cp, rsb_nnz_idx_t ldC, rsb_time_t maxdt, rsb_int_t mintimes, enum rsb_op_t op, rsb_int_t maxtimes, int verbose, rsb_int_t *tnp, struct rsb_ts_t * tstp)
{
	/*!
	       	Iterates calling the computing kernel.
		Does not modify matrix structure.
		If RSB_DT_THREADS_TUNE_TNP(tnp) it modifes temporarily executing threads count.
		Unless RSB_SOME_ERROR(errval), shall give in output *tpop > RSB_TIME_ZERO.

	 	TODO: can consider: const rsb_bool_t want_inner_flush = RSB_BOOL_FALSE; if(want_inner_flush == RSB_BOOL_TRUE) RSB_DO_ERROR_CUMULATE(errval,rsb__flush_cache(0));
	*/
	rsb_err_t errval = RSB_ERR_NO_ERROR;
	rsb_err_t kerrval = RSB_ERR_NO_ERROR;
	rsb_time_t it = RSB_TIME_ZERO, ct = RSB_TIME_ZERO;	/* initial/current time */
	rsb_time_t dt = it, tt = RSB_TIME_ZERO; /* elapsed (delta) / total  time */
	rsb_time_t bt = RSB_CONST_IMPOSSIBLY_BIG_TIME, wt = RSB_TIME_ZERO; /* best / worst  time */
	rsb_time_t ss = RSB_TIME_ZERO; /* sum of squares */
	const rsb_time_t jf = RSB_CACHED_TIMER_GRANULARITY; /* jiffie */
	const rsb_time_t mindt = RSB_AT_SAMPLE_MINDT; /* minimal time */
	rsb_int_t times = 0;
	rsb_int_t otn = 0;

	RSB_DO_REINIT_SINGLE_VALUE_GET(RSB_IO_WANT_EXECUTING_THREADS,&otn,errval);

	if(tnp)
	{
		RSB_DO_REINIT_SINGLE_VALUE_SET(RSB_IO_WANT_EXECUTING_THREADS, tnp, errval);
		if(RSB_SOME_ERROR(errval))
		{
		       	RSB_PERR_GOTO(err,RSB_ERRM_ES)
		}
	}

	dt = it = rsb_time();

	switch(op)
	{
		case(rsb_op_spmv):
		do
		{
			kerrval |= rsb__do_spmm(transA, alphap, mtxAp, nrhs, order, Bp, ldB, betap, Cp, ldC, RSB_OP_FLAG_DEFAULT);
			RSB_SAMPLE_STAT(it,ct,dt,tt,bt,wt,ss,jf,times);
		}
		while(RSB_REPEAT(ct-it,times,mindt,mintimes,maxdt,maxtimes));
		break;
		case(rsb_op_spsvlt):
		do
		{
			kerrval |= rsb__do_spsm(transA, alphap, mtxAp, nrhs, order, betap, Bp, ldB, Cp, ldC);
			RSB_SAMPLE_STAT(it,ct,dt,tt,bt,wt,ss,jf,times);
		}
		while(RSB_REPEAT(ct-it,times,mindt,mintimes,maxdt,maxtimes));
		break;
		case(rsb_op_nop):
		/* no operation, no error */
		break;
		default:
		/* no operation, error */
	      	RSB_PERR_GOTO(done,RSB_ERRM_ES)
	}

done:
	if(RSB_SOME_ERROR(kerrval))
	{
		RSB_ERROR("Benchmarked kernel error (propagating it further): %s\n", rsb__get_errstr_ptr(kerrval));
	}

	if(tnp)
	{
		RSB_DO_REINIT_SINGLE_VALUE_SET(RSB_IO_WANT_EXECUTING_THREADS,&otn,errval);
		if(RSB_SOME_ERROR(errval))
		{
		       	RSB_PERR_GOTO(err,RSB_ERRM_ES)
		}
	}
err:

#if RSB_AT_WANT_BEST_TIME
	RSB_ASSIGN_IF(tpop,bt)
#else
	RSB_ASSIGN_IF(tpop,( tt ) / times)
#endif
	RSB_ASSIGN_IF(timesp,times)

#if RSB_ALLOW_INTERNAL_GETENVS
	if(getenv("RSB_VERBOSE_TUNING_INNER"))
		verbose = rsb__util_atoi(getenv("RSB_VERBOSE_TUNING_INNER"));
#endif /* RSB_ALLOW_INTERNAL_GETENVS*/

	if( RSB_UNLIKELY( verbose >= RSB_AUT0_TUNING_VERBOSE ) )
		RSB_STAT_DUMP(it,tnp?*tnp:otn,ct,dt,tt,bt,wt,ss,times);
	RSB_STAT_TAKE(it,tnp?*tnp:otn,ct,dt,tt,bt,wt,ss,times,tstp);

	RSB_DO_ERROR_CUMULATE(errval,kerrval);
	return errval;
}

struct rsb_oas_t /* spmv/spsv/spmm/spsm op args structure, for rsb__dbf_t */
{
       	enum rsb_op_t op; /* e.g. may expand to combined N+T transposition benchmarking ... */
       	rsb_trans_t transA;
       	const void * alphap;
       	rsb_coo_idx_t nrhs;
       	rsb_flags_t order;
       	const void * Bp;
       	rsb_nnz_idx_t ldB;
       	const void * betap;
       	void * Cp;
       	rsb_nnz_idx_t ldC;
};

struct rsb_aoa_t /* autotuning output args */
{
	rsb_time_t otpo, btpo; /* original / best time per operation */
	rsb_int_t cc; /* clones count (total) */
	rsb_time_t ct; /* clones time */
	struct rsb_ts_t otpos, btpos; /* original/best time per operation stats */
};

struct rsb_aia_t /* autotuning input args */
{
	/* non RSB specific */
	rsb_time_t maxt; /* max operation iterations time */
	rsb_int_t mintimes; /* min operation iterations */
	rsb_int_t maxtimes; /* max operation iterations */
	const char * rfpfx; /* rendering filename prefix */
	const char * mtxns; /* matrix name string */
	rsb_perf_t ofe; /* operation flops estimate */
	int verbose; /* TODO: indicate valid values */
	rsb_int_t *tnp; /* threads number */
	/* RSB specific */
	const struct rsb_mtx_t * mtxAp;
	rsb_int_t maxr; /* max rounds */
	rsb_int_t maxms; /* max merge steps */
	rsb_int_t maxss; /* max split steps */
	rsb_int_t continue_partitioning; /* continue partitioning and pretend faster: mostly here for debug */
 	/* rsb_time_t *iotp, rsb_time_t *tttp, rsb_time_t *rctp */
};

#if 0
struct rsb_ats_t /* autotuning times struct (FIXME: still unused) */
{
 	rsb_time_t *iotp;
	rsb_time_t *tttp;
	rsb_time_t *rctp;
};
#endif

struct rsb_bos_t /* benchmarkable operation structure, for rsb__dbf_t */
{
	struct rsb_oas_t oas;
       	rsb_int_t *timesp;
       	rsb_time_t maxdt;
       	rsb_int_t mintimes;
       	rsb_int_t maxtimes;
       	int verbose;
       	rsb_int_t *tnp;
	struct rsb_ts_t*otsp;
};

typedef rsb_err_t (*rsb_bf_t) (rsb_time_t *tpop, const struct rsb_mtx_t * mtxAp, const void*p) ;

static rsb_err_t rsb__dbf_t(rsb_time_t *tpop, const struct rsb_mtx_t * mtxAp, const void*p)
{
	rsb_err_t errval = RSB_ERR_NO_ERROR;
	struct rsb_bos_t bos;

	if(!p)
	{
		errval = RSB_ERR_BADARGS; 
		RSB_PERR_GOTO(err,RSB_ERRM_ES)
	}

	bos = *(const struct rsb_bos_t*) p;

	errval = rsb__do_bench_spxm(tpop, bos.timesp, bos.oas.transA, bos.oas.alphap, mtxAp, bos.oas.nrhs, bos.oas.order, bos.oas.Bp, bos.oas.ldB, bos.oas.betap, bos.oas.Cp, bos.oas.ldC, bos.maxdt, bos.mintimes, bos.oas.op, bos.maxtimes, bos.verbose, bos.tnp, bos.otsp);
err:
	return errval;
}

static rsb_err_t rsb__do_bench_spxm_t(rsb_time_t *otpop, rsb_time_t *tpop, rsb_int_t *tnp, const struct rsb_mtx_t * mtxAp, rsb_time_t maxt, rsb_int_t mintimes, rsb_int_t maxtimes, int verbose, struct rsb_oas_t*oasp, struct rsb_tattr_t*ttrp, struct rsb_ts_t * tstp)
{
	/*!
		Thread tuning.
	       	No matrix structure modification is attempted.

		FIXME: Unfinished: errors may leave library status undefined.
		FIXME: expand usage of rsb__dbf_t.
	 */
	rsb_int_t otn = 0;
	rsb_err_t errval = RSB_ERR_NO_ERROR;
	rsb_int_t btimes = 0, times = 0;
	rsb_time_t ctpo = RSB_TIME_ZERO, btpo = RSB_TIME_ZERO, otpo = RSB_TIME_ZERO;
	rsb_int_t mt = rsb__set_num_threads(/*RSB_THREADS_GET_MAX*/RSB_THREADS_GET_MAX_SYS);
	rsb_perf_t ofe = RSB_PERF_ZERO; /* operation flops estimate */
	struct rsb_ts_t ots;
	struct rsb_ts_t btpos;

	RSB_BZERO_P(&ots);

	if(!tpop || !oasp )
	{
		errval = RSB_ERR_BADARGS; 
		RSB_PERR_GOTO(err,RSB_ERRM_ES)
	}

	if( !oasp->alphap || !oasp->betap || !mtxAp || !oasp->Bp || !oasp->Cp )
	{
		errval = RSB_ERR_BADARGS; 
		RSB_PERR_GOTO(err,RSB_ERRM_ES)
	}

	ofe = rsb__estimate_mflops_per_op_spmv_uaua(mtxAp) * oasp->nrhs;

	RSB_THREAD_STATS
	RSB_DO_REINIT_SINGLE_VALUE_GET(RSB_IO_WANT_EXECUTING_THREADS,&otn,errval);
	if(RSB_SOME_ERROR(errval))
	{
	       	RSB_PERR_GOTO(err,RSB_ERRM_ES)
	}

	rsb__tattr_init(ttrp,mtxAp,mtxAp->nr,mtxAp->nnz,mtxAp->typecode,mtxAp->flags,oasp->nrhs);

	if( RSB_DT_THREADS_TUNE_TNP(tnp) )
	{
		const rsb_int_t tinc = -1; /* FIXME: experimental */
		rsb_int_t mit, mat = mt, cut;
		rsb_time_t wtpo = RSB_TIME_ZERO ; /* worst time per operation */
		rsb_time_t ltpo = RSB_TIME_ZERO ; /* last time per operation */
		const rsb_time_t tst = -rsb_time();
		rsb_int_t diar = 0; /* degradations in a row */
		rsb_int_t mdiar = 0; /* maximal degradations in a row */

		if(tinc > 1)
			mit = 1, mat = mt, mdiar = mt+1;
		else
			mit = mt, mat = 1, mdiar = 2;

		if(verbose > 1)
			RSB_STDOUT("Will test in range (%d to %d) threads, with increase of %d and tolerate max %d consecutive degradations.\n",mit,mat,tinc,mdiar);

		for(cut=mit;tinc>0?(cut<=mat):(cut>=1);cut+=tinc)
		{
#if 0
			errval = rsb__do_bench_spxm(&ctpo, &times, oasp->transA, oasp->alphap, mtxAp, oasp->nrhs, oasp->order, oasp->Bp, oasp->ldB, oasp->betap, oasp->Cp, oasp->ldC, maxt, mintimes, oasp->op, maxtimes, verbose, &cut, NULL);
#else
			struct rsb_bos_t bos = { *oasp, &times, maxt, mintimes, maxtimes, verbose, &cut, &ots };
			const rsb_bf_t bf = &rsb__dbf_t; /* FIXME: in the future, bf shall be an argument */

			bf(&ctpo,mtxAp,&bos);
#endif

			if( ltpo == RSB_TIME_ZERO ) 
				ltpo = ctpo;

			if( ctpo > wtpo || wtpo == RSB_TIME_ZERO ) 
				wtpo = ctpo;

			if( ctpo < btpo || btpo == RSB_TIME_ZERO ) 
			{
				*tnp = cut;
				btpo = ctpo;
				btpos = ots;
				btimes = times;
			}

			/* if( ctpo > ltpo ) */
			if( ctpo > btpo )
				diar ++;
			else
				diar = 0;

			rsb__tattr_sets(ttrp,otn,cut,ctpo,*tnp,times);

			if(verbose)
			RSB_STDOUT("%c %d threads: %0.4lgs  (%0.2lg Mflops) (%d/%d degradations so far)  %c%c\n",cut==otn?'~':' ',cut,ctpo,ofe/ctpo,
					diar,mdiar,
					ctpo<otpo ? '+' : (ctpo==otpo ?'=':'-'), ctpo<btpo ? '+' : ' ' );
			ltpo = ctpo;

			if( diar >= mdiar )
				break;
		} /* cut */

		if(verbose)
			RSB_STDOUT("Best threads choice is %d; starting threads were %d; max speed gap is %0.2lgx; search took %0.2lgs.\n",*tnp,otn,wtpo/btpo,tst+rsb_time());
		goto done;
	} /* RSB_DT_THREADS_TUNE_TNP(tnp) */

	/* FIXME: need a better tuning + exploration mechanism. */

	if(1)
	{
#if 0
		errval = rsb__do_bench_spxm(&btpo, &btimes, oasp->transA, oasp->alphap, mtxAp, oasp->nrhs, oasp->order, oasp->Bp, oasp->ldB, oasp->betap, oasp->Cp, oasp->ldC, maxt, mintimes, oasp->op, maxtimes, verbose, NULL,NULL);
#else
		const struct rsb_bos_t bos = { *oasp, &btimes, maxt, mintimes, maxtimes, verbose, NULL, &ots };
		const rsb_bf_t bf = &rsb__dbf_t; /* FIXME: in the future, bf shall be an argument */

		RSB_DO_ERROR_CUMULATE(errval, bf(&btpo,mtxAp,&bos)); /* replaces rsb__do_bench_spxm */
#endif
		if(RSB_SOME_ERROR(errval))
		{
		       	RSB_PERR_GOTO(err,RSB_ERRM_ES)
		}
		otpo = btpo;
		btpos = ots;
		rsb__tattr_sets(ttrp,otn,otn,btpo,otn,btimes);

		if(verbose)
			RSB_STDOUT("Reference operation time is %lg s (%0.4lg Mflops) with %d threads.\n",otpo,ofe/otpo,otn);
		goto done;
	}
done:
	*tpop = btpo;
	RSB_ASSIGN_IF(otpop,otpo)
	RSB_ASSIGN_IF(tstp,btpos)
	goto err;
err:
	return errval;
}

static rsb_err_t rsb__do_tune_spxm_round( rsb_real_t*sfp, struct rsb_mtx_t ** mtxOpp, rsb_real_t*osvalp, const struct rsb_aia_t*aiap, rsb_time_t *iotp, rsb_time_t *tttp, rsb_time_t *rctp, struct rsb_aoa_t*aoap, struct rsb_oas_t*oasp, struct rsb_attr_t *attrp )
{
	/*
	 * A round of auto-tuning.
	 *
	 * Can check for different blockings (by re-invoking the constructor).
	 *
	 * See RSB_DT_SAME_THREADS_TNP RSB_DT_THREADS_TUNE_TNP RSB_DT_SPEC_THREADS_TNP.
	 *
	 * if (mtxOpp && *mtxOpp): its poined area may be freed.
	 * otherwise wiil use aiap->mtxAp and update *mtxOpp.
	 * If mtxOpp==NULL, only thread tuning will be possible.
	 */
	rsb_err_t errval = RSB_ERR_NO_ERROR;
	rsb_real_t rsf = RSB_REAL_ZERO; /* result speedup factor */
	double bsval = RSB_REAL_ZERO, /* best subdivision [multiplier] value */
	       sval = RSB_INVALID_SVAL, ospv = RSB_REAL_ZERO;
	const double svalr = 4.0, svalm = 2.0;
	const rsb_time_t omaxt = aiap->maxt/(svalr*svalr);
	rsb_time_t rtpo = RSB_TIME_ZERO;
	struct rsb_mtx_t * mtxOp = NULL, *mtxCp = NULL; /* Output, Clone */
	const struct rsb_mtx_t * mtxBp = NULL; /* Best */
	rsb_flags_t sflags = RSB_FLAG_QUAD_PARTITIONING | RSB_FLAG_USE_HALFWORD_INDICES; /* structural flags */
	rsb_int_t otn = 0,mt/*,tries=2*rsb_get_num_threads()*/;
	rsb_perf_t ofe = RSB_PERF_ZERO; /* operation flops estimate */
	rsb_time_t iot = RSB_TIME_ZERO,	/* I/O time */
                   att = rsb_time(), 	/* auto tuning time */
		   rct = RSB_TIME_ZERO; /* reference cloning  time */
	rsb_thread_t btn = 0; /* best      threads number */
	rsb_thread_t rtn = 0; /* requested threads number */
	rsb_thread_t ctn = 0; /* current threads number */
	rsb_bool_t only_one_sample = (mtxOpp == NULL) && RSB_DT_SAME_THREADS_TNP(aiap->tnp);
	struct rsb_tattr_t * ttrp = NULL; /* thread tuning record pointer  */
	struct rsb_aoa_t aoa; /* autotuning output args */
	struct rsb_ts_t rtpos; /* running time per operation structure */

	RSB_BZERO_P(&aoa);

	mt = rsb__set_num_threads(/*RSB_THREADS_GET_MAX*/RSB_THREADS_GET_MAX_SYS); 

	if(mtxOpp && *mtxOpp)
		mtxBp = *mtxOpp;
	else
		mtxBp = aiap->mtxAp;

	ofe = rsb__estimate_mflops_per_op_spmv_uaua(mtxBp) * oasp->nrhs; /* FIXME: oasp dep. */

	if( !oasp || /*!mtxAp ||*/ /* *mtxOpp || */ (!RSB_IS_VALID_TRANS(oasp->transA)) || RSB_INVALID_COO_COUNT(oasp->nrhs)  /*|| !RSB_IS_VALID_THREAD_SPEC(ltn) || !RSB_IS_VALID_THREAD_SPEC(utn)*/)
	{
		errval = RSB_ERR_BADARGS;
		RSB_PERR_GOTO(ret,RSB_ERRM_ES);
	}
	
	RSB_DO_REINIT_SINGLE_VALUE_GET(RSB_IO_WANT_EXECUTING_THREADS,&otn,errval);
	
	if(RSB_DT_SPEC_THREADS_TNP(aiap->tnp))
	       	ctn = rtn = *aiap->tnp;

	if(aiap->verbose)
	{
		if( mtxOpp == NULL )
		{
			if( RSB_DT_SAME_THREADS_TNP(aiap->tnp) )
				; /* RSB_STDOUT("Will sample performance at default threads.\n"); */
			else
				RSB_STDOUT("Started tuning inner round on given matrix instance.\n");
		}
		else
		{
			RSB_STDOUT("Started tuning inner round: will search for an optimal matrix instance."); /* with or without threads ?? */
			RSB_STDOUT("\n""Starting with requested %d threads ; current default %d ; at most %d.\n",rtn,otn,mt);
		}
	}

	//do
	{
		RSB_ATTR_PREPARE(ttrp,attrp,errval);
		if(RSB_SOME_ERROR(errval))
		{
		       	RSB_PERR_GOTO(err,RSB_ERRM_EM); 
		}
		errval = rsb__do_bench_spxm_t(&aoa.otpo,&aoa.btpo, /*(mtxOpp) ? NULL :*/ aiap->tnp, mtxBp, omaxt, aiap->mintimes, aiap->maxtimes, aiap->verbose, oasp, ttrp, &aoa.otpos);
		rtpos = aoa.btpos = aoa.otpos; /* first sample is also best so far */
		RSB_ATTR_ADVANCE(ttrp,attrp,0.0);

		if(RSB_SOME_ERROR(errval))
		{
		       	RSB_PERR_GOTO(err,RSB_ERRM_EM); 
		}
	}

	if( mtxOpp == NULL )
	{
		if(aiap->verbose)
		{
			if(aiap->tnp && *aiap->tnp!=0)
				RSB_STDOUT("Seems like best threads choice with this matrix is %d threads, with speedup %lg x.\n",*aiap->tnp,aoa.otpo/aoa.btpo);
		}
		RSB_DO_REINIT_SINGLE_VALUE_GET(RSB_IO_WANT_SUBDIVISION_MULTIPLIER,&bsval,errval);
		goto ann; /* job done */
	}

	btn = otn;

	rct -= rsb_time();
	/* Make a pure COO clone beforehand, so to make sure not to get 1 to 1 cloned copies after. */
	errval = rsb__mtx_clone(&mtxCp,mtxBp->typecode,RSB_TRANSPOSITION_N,NULL,mtxBp,RSB_DO_FLAG_FILTEROUT(mtxBp->flags,sflags));
	rct += rsb_time();
	/*if(RSB_VERBOSE_FOR_SAVING(aiap->verbose))
		RSB_ATR_SAVE(errval,rfpfx,mtxns,mtxBp->typecode,transA,nrhs,"r",0,bsval,mtxCp,iot,op);
	*/

	if(RSB_SOME_ERROR(errval))
	{
		RSB_PERR_GOTO(ret,RSB_ERRM_NL);
	}

	if(RSB_VERBOSE_FOR_MSG(aiap->verbose))
	{
		RSB_STDOUT("Copied to a disposable COO clone: ");
		RSB_STDOUT(RSB_PRINTF_MATRIX_AT_SUMMARY_ARGS(mtxCp));
		RSB_STDOUT("\n");
	}
	/* Save subdivision multiplier values. */
	if(osvalp)
		ospv = *osvalp;
	else
		RSB_DO_REINIT_SINGLE_VALUE_GET(RSB_IO_WANT_SUBDIVISION_MULTIPLIER,&ospv,errval);
	
	if(RSB_SOME_ERROR(errval))
	{
		RSB_PERR_GOTO(ret,RSB_ERRM_NL);
       	}
	/* TODO: shall make sure that RSB_REAL_ZERO is never returned */
	bsval = ospv;

	aoa.otpo = aoa.btpo; /* best time with this matrix */

	if(aiap->verbose)
		RSB_STDOUT("Starting autotuning stage, with subdivision of %lg (current threads=%d, requested threads=%d, max threads = %d).\n",ospv,otn,rtn,mt);

	/* For different thread counts */
	/* ... */
	// for(tni=...;tni<...;++tni)
	/* For different cache blocking solutions */
	for(sval=(ospv/svalr);sval<=(ospv*svalr);sval*=svalm)
	{
		rsb_time_t cct = - rsb_time();
		/* Set cache blocking */
		RSB_DO_REINIT_SINGLE_VALUE_SET(RSB_IO_WANT_SUBDIVISION_MULTIPLIER,&sval,errval);
		if(RSB_SOME_ERROR(errval)){ RSB_PERR_GOTO(err,RSB_ERRM_EM); }

		if(otn != ctn)
			RSB_DO_REINIT_SINGLE_VALUE_SET(RSB_IO_WANT_EXECUTING_THREADS,&ctn,errval);
		rct -= rsb_time();
		errval = rsb__mtx_clone(&mtxOp,RSB_NUMERICAL_TYPE_SAME_TYPE,RSB_TRANSPOSITION_N,NULL,mtxCp,mtxBp->flags|sflags);
		rct += rsb_time();
		if(RSB_SOME_ERROR(errval))
			{ RSB_PERR_GOTO(err,RSB_ERRM_EM); }
		if(otn != ctn)
			RSB_DO_REINIT_SINGLE_VALUE_SET(RSB_IO_WANT_EXECUTING_THREADS,&otn,errval);
		if(RSB_VERBOSE_FOR_MSG(aiap->verbose))
		{
			RSB_STDOUT("Considering a candidate clone: ");
			RSB_STDOUT(RSB_PRINTF_MATRIX_AT_SUMMARY_ARGS(mtxOp));
			RSB_STDOUT("\n");
		}

		if(RSB_VERBOSE_FOR_SAVING(aiap->verbose))
			RSB_ATR_SAVE(errval,aiap->rfpfx,aiap->mtxns,mtxBp->typecode,oasp->transA,oasp->nrhs,sval,"r",0,ctn,mtxOp,iot,oasp->op);

		/* FIXME: need to check whether matrix is effectively more/less subdivided; if not then skipping */
		RSB_ATTR_PREPARE(ttrp,attrp,errval);
		if(RSB_SOME_ERROR(errval))
			{ RSB_PERR_GOTO(err,RSB_ERRM_IE); }
		errval = rsb__do_bench_spxm_t( NULL, &rtpo, aiap->tnp ? &ctn : NULL, mtxOp, omaxt, aiap->mintimes, aiap->maxtimes, aiap->verbose, oasp, ttrp, &rtpos);
		cct += rsb_time();
		RSB_ATTR_ADVANCE(ttrp,attrp,cct);
		if(RSB_SOME_ERROR(errval)){ RSB_PERR_GOTO(err,RSB_ERRM_EM); }

		/* Free the "slower" matrix, keep the "optimized" one. */
		if(aiap->verbose)
		{
			rsb_real_t bpn = 0;
			rsb_blk_idx_t lfc = 0;
			errval = rsb__do_mtx_get_info(mtxOp, RSB_MIF_INDEX_STORAGE_IN_BYTES_PER_NNZ__TO__RSB_REAL_T, &bpn);
			if(RSB_SOME_ERROR(errval)){ RSB_PERR_GOTO(err,RSB_ERRM_EM) }
			errval = rsb__do_mtx_get_info(mtxOp, RSB_MIF_LEAVES_COUNT__TO__RSB_BLK_INDEX_T, &lfc);
			if(RSB_SOME_ERROR(errval)){ RSB_PERR_GOTO(err,RSB_ERRM_EM) }
			// RSB_STDOUT("For subdivision %lg (%d leaves, %5.4lg bytes/nz, %d threads), challenging %lg s with %lg s (speedup %lg x)\n",sval,(int)lfc,bpn,(int)(tnp?*tnp:otn),rbtpo,rtpo,rbtpo/rtpo);
			RSB_STDOUT("Challenging best inner round reference (%lg s/%d threads) with: subdivision %lg, %d leaves, %5.4lg bytes/nz, %lg s/%d threads (speedup %lg x), same?%c.\n", aoa.btpo,(int)btn, sval, (int)lfc, bpn, rtpo,(int)ctn, aoa.btpo/rtpo,((mtxOp==*mtxOpp)?'y':'n'));
		}

		if(rtpo < aoa.btpo)
		{
			if( mtxOp != *mtxOpp )
			{
				if(/*RSB_VERBOSE_FOR_MSG */(aiap->verbose) && *mtxOpp)
				{
					RSB_STDOUT("New good candidate clone found: discarding old candidate clone: ");
					RSB_STDOUT(RSB_PRINTF_MATRIX_AT_SUMMARY_ARGS(*mtxOpp));
					RSB_STDOUT("\n");
				}
#if RSB_AT_DESTROYS_MTX
				RSB_MTX_FREE(*mtxOpp);
#endif /* RSB_AT_DESTROYS_MTX */
			}
			*mtxOpp = mtxOp;
			mtxBp = mtxOp;
			mtxOp = NULL;
			aoa.btpo = rtpo;
			aoa.btpos = rtpos;
			bsval = sval;
			btn = ctn;
		}
		else
		{
			if(/*RSB_VERBOSE_FOR_MSG*/(aiap->verbose))
			{
				RSB_STDOUT("New candidate clone performs slowly; discarding it: ");
				RSB_STDOUT(RSB_PRINTF_MATRIX_AT_SUMMARY_ARGS(mtxOp));
				RSB_STDOUT("\n");
			}
			RSB_MTX_FREE(mtxOp);
		}
		
		if( aiap->tnp && *aiap->tnp > mt )
		{
			/*  */
			break;
		}
	}

	if(aiap->tnp)
	{
		*aiap->tnp = btn;
	}
ann:
	if( !only_one_sample )
	if(aiap->verbose)
	{
		RSB_STDOUT("%s %s performance",(mtxOpp==NULL && RSB_DT_SAME_THREADS_TNP(aiap->tnp))?"Measured":"Best",( oasp->op == rsb_op_spmv ?"sparse multiply":"sparse solve"));
		if(mtxOpp!=NULL)
			RSB_STDOUT(" with subdivision multiplier of %lg",bsval);
		if(aiap->tnp && *aiap->tnp!=0)
			RSB_STDOUT(" using %d threads",(int)btn );
		RSB_STDOUT(": %lg Mflops.\n", (ofe)/(aoa.btpo) );
	}
	/* Pick up the best solution for return. */
err:
	/* FIXME: in some error cases the RSB_IO_WANT_SUBDIVISION_MULTIPLIER seems not to be be resetted.  */
	if(RSB_SOME_ERROR(errval))
	{
	       	RSB_MTX_FREE(mtxOp);
	}
	/* Restore subdivision multiplier values. */
	if( ospv != RSB_REAL_ZERO )
		RSB_DO_REINIT_SINGLE_VALUE_SET(RSB_IO_WANT_SUBDIVISION_MULTIPLIER,&ospv,errval);
	RSB_ASSIGN_IF(osvalp,bsval);
ret:
	rsf = aoa.otpo / aoa.btpo;
	att = ( rsb_time() - att ) - iot;
	if( aiap->verbose && RSB_DT_THREADS_TUNE_TNP(aiap->tnp) )
	{
		RSB_STDOUT("Auto-tuning inner round complete. Gained a total speedup of %lg x (= %lg : %lg)).\n",rsf,aoa.otpo,aoa.btpo);
	}
	RSB_ASSIGN_IF(sfp,rsf) /* speedup factor */
	RSB_ASSIGN_IF(iotp,iot)
	RSB_ASSIGN_IF(tttp,att)
	RSB_ASSIGN_IF(rctp,rct)
	RSB_ASSIGN_IF(aoap,aoa)
	RSB_MTX_FREE(mtxCp);
	RSB_DO_ERR_RETURN(errval)
}

#define RSB_TIME_BETTER_THRESHOLD 0.01
/* #define RSB_TIME_BETTER_THAN(TNEW,TOLD) ( (TNEW) < ((TOLD)*(1.0+RSB_TIME_BETTER_THRESHOLD)) ) */
#define RSB_TIME_BETTER_THAN(TNEW,TOLD) ( (TOLD) > ((TNEW)*(1.0+RSB_TIME_BETTER_THRESHOLD)) )
#define RSB_TIME_DEFINITELY_WORSE_THAN(TNEW,TOLD) ( (TNEW) > ((TOLD)*(1.0+0.02)) )
#define RSB_TIME_MUCH_WORSE_THAN(TNEW,TOLD) ( (TNEW) > ((TOLD)*(1.0+0.20)) )

static rsb_err_t rsb__rblk_tune_inner(struct rsb_mtx_t **mtxOpp, rsb_submatrix_idx_t ms, const struct rsb_aia_t*aiap, rsb_int_t*epsp, struct rsb_ts_t*otposp, struct rsb_ts_t*btposp, rsb_time_t*otpopp, rsb_int_t *tnp,  rsb_time_t*btpopp, rsb_int_t *ccp, rsb_time_t *ctp, rsb_bool_t want_split, struct rsb_oas_t*oasp, struct rsb_attr_t *attrp )
{
	/**
	 * Re-blocking (split/merge) based autotuning.
	 * This function makes a clone of the original matrix and works on that.
	 * Works by:
	 *  - taking a first clone of a matrix, benchmarking it
	 *  - taking a second clone and reblocking it
	 *  - keeping the best and iterating
	 * So in principle, three instances may be at a given time.
	 * Additionally, one should keep in mind that  rsb__mtx_split/rsb__leaves_merge  use temporary memory.
	 * 
	 * FIXME: in case of an error (e.g. insufficient memory), it is still possible that a new (repartitioned) matrix will be returned, and the error return code will be reset.
	 * TODO: ideally, one should have reversible merge and split operators and no cloning, but reversal only would be necessary (maybe in librsb-1.3); one may create such "reversal" structures: an array of transition actions.
	 * TODO: support for truly (no clone) in place tuning, via the "reversal" based algorithms.
	 */
	rsb_err_t errval = RSB_ERR_NO_ERROR;
	rsb_time_t tmt = RSB_TIME_ZERO, tst = RSB_TIME_ZERO, /* total [merge/split,sort] time */
		   tlt = RSB_TIME_ZERO, tat = RSB_TIME_ZERO, /* total [elapsed,analysis] */
		   tpt = RSB_TIME_ZERO, tct = RSB_TIME_ZERO, /* total [partitioning(=splitting/merging),cloning] time */
		   st = RSB_TIME_ZERO, at = RSB_TIME_ZERO, lt = RSB_TIME_ZERO; /* sort, analysis, elapsed, */
	rsb_submatrix_idx_t eps = 0; /* effective (executed) partitioning steps */
	rsb_submatrix_idx_t pps = 0; /* profitable (executed) partitioning steps */
	rsb_submatrix_idx_t nsbp,nsap,nsdp; /* number of submatrices [before/after/during] partition */
	rsb_submatrix_idx_t cc = 0; /* cloning count */
#ifndef NDEBUG
	rsb_submatrix_idx_t nt = rsb_get_num_threads(); /*  */
#endif
	rsb_submatrix_idx_t manp = 0; /* maximum number of passes (merges/splits) */
	int wom = 0, db = /* tnp == NULL ? 1 : */ 0; /* want only merging, do bench */
#if 0
	int wom = da = 0; /* do autotuning */
#endif
	rsb_time_t tpop = RSB_TIME_ZERO, otpop = RSB_TIME_ZERO, btpop = RSB_TIME_ZERO; /* time per operation (t.p.op.), old t.p.op., best t.p.op. */
	struct rsb_mtx_t *mtxWp = NULL, *mtxBp = NULL; /* matrices: work, best  */
	rsb_int_t diar = 0; /* degradations in a row */
	rsb_int_t niar = 0; /* neutral in a row */
	const rsb_int_t mdiar = RSB_MIN(3,ms); /* max [tolerated] degradations in a row [before skipping tuning] */
	rsb_int_t itn = 0, otn = 0; /* initial/optimal threads number */
	const rsb_char_t*const pos = (want_split) ? "split" : "merge"; /* partition operation string  */
	const int wv = /*1*/aiap->verbose; /* want verbose */
	double mpoe = 0.0; /* estimated megaflops per operation */
	struct rsb_tattr_t * ttrp = NULL; /* thread tuning record pointer  */
	rsb_time_t iot = RSB_TIME_ZERO;	/* I/O time */
        const struct rsb_mtx_t*mtxAp = ( aiap->mtxAp ? aiap->mtxAp : *mtxOpp );
	struct rsb_ts_t otpos, btpos;
	struct rsb_ts_t tpos;

	RSB_BZERO_P(& tpos);
	RSB_BZERO_P(&otpos);
	RSB_BZERO_P(&btpos);

        RSB_DEBUG_ASSERT(mtxAp);
	RSB_ASSERT(nt>0);

#if RSB_ALLOW_INTERNAL_GETENVS
	if(getenv("RSB_MERGE_WOM")) wom = rsb__util_atoi(getenv("RSB_MERGE_WOM"));
	if(getenv("RSB_MERGE_DB")) db = rsb__util_atoi(getenv("RSB_MERGE_DB"));
	/*if(getenv("RSB_MERGE_DA")) da = rsb__util_atoi(getenv("RSB_MERGE_DA")); */
#endif /* RSB_ALLOW_INTERNAL_GETENVS */

	if(!mtxAp)
	{
		errval = RSB_ERR_BADARGS;
		RSB_PERR_GOTO(ret,RSB_ERRM_ES);
	}

	RSB_ASSERT( ! ( aiap->mtxAp != NULL && mtxOpp != NULL ) );
	RSB_ASSERT( mtxAp->all_leaf_matrices_n > 0 );

	if( aiap->mtxAp != NULL && mtxOpp != NULL )
	{
		errval = RSB_ERR_BADARGS;
		RSB_PERR_GOTO(ret,RSB_ERRM_ES);
	}

	nsap = nsbp = mtxAp->all_leaf_matrices_n;
#if 0
	{
		const rsb_int_t maxpasses = 20; /* max number of merging passes (FIXME: this is an arbitrary constant ...) */
		manp = RSB_MAX( nt, nsbp / maxpasses );
	}
#endif
	mpoe = aiap->ofe;

	tct -= rsb_time();
	mtxWp = rsb__mtx_clone_simple(mtxAp);
	tct += rsb_time();
	cc++;

	if(!mtxWp)
	{
		errval = RSB_ERR_ENOMEM;
		RSB_PERR_GOTO(ret, RSB_ERRM_E_MTXAP);
	}

	RSB_ATTR_PREPARE(ttrp,attrp,errval);

	if(RSB_SOME_ERROR(errval))
	{
	       	RSB_PERR_GOTO(err,RSB_ERRM_EM); 
	}

	if(db)
		errval = rsb__do_bench_spxm(&tpop,NULL,oasp->transA,oasp->alphap,mtxWp,oasp->nrhs,oasp->order,oasp->Bp,oasp->ldB,oasp->betap,oasp->Cp,oasp->ldC,aiap->maxt,aiap->mintimes,oasp->op,aiap->maxtimes,wv,NULL,&tpos);
	else
	{
		RSB_ASSIGN_IF_SP(otn,tnp);
		errval = rsb__do_bench_spxm_t(NULL,&tpop,&otn,mtxAp,aiap->maxt,aiap->mintimes,aiap->maxtimes,wv,oasp,ttrp,&tpos);
		itn = otn;
	}
	RSB_ATTR_ADVANCE(ttrp,attrp,0.0);
	RSB_ASSERT( ( tpop != RSB_TIME_ZERO ) || RSB_SOME_ERROR(errval) );

	if(RSB_SOME_ERROR(errval))
	{
	       	RSB_PERR_GOTO(err,RSB_ERRM_EM); 
	}

	btpop = otpop = tpop;
	btpos = otpos = tpos;

	if(wv)
		RSB_STDOUT("Starting %s (%s threads) based auto-tuning procedure (transA=%c, nrhs=%ld, order=%s) (max %ld steps, inclusive %ld grace steps) on: ",pos,RSB_DT_SAME_THREADS_TNP(tnp)?"same":(RSB_DT_THREADS_TUNE_TNP(tnp)?"and":"user-supplied"),RSB_TRANSPOSITION_AS_CHAR(oasp->transA),(long int)oasp->nrhs,RSB__ORDER_AS_STRING(oasp->order),(long int)ms,(long int)mdiar),
		RSB_STDOUT(RSB_PRINTF_MATRIX_AT_SUMMARY_ARGS(mtxWp)),
	       	RSB_STDOUT(" (tpop: %0.4lg  Mflops: %3.3lf)\n",tpop,mpoe/tpop);

	tpt -= rsb_time();

	for(eps=0;eps<ms && ( want_split || ( mtxWp->all_leaf_matrices_n>manp ) ) ;++eps)
	{
		rsb_int_t ctn = 0;
		rsb_time_t cct = - rsb_time();

		tpop = RSB_TIME_ZERO;
		tmt -= rsb_time();
		nsdp = mtxWp->all_leaf_matrices_n;

		if(want_split)
		{
			errval = rsb__mtx_realloc_with_spare_leaves(&mtxWp, RSB_TMP_OVERALLOC_MTX*rsb__submatrices(mtxWp));
			if(RSB_SOME_ERROR(errval))
			{
				RSB_PERR_GOTO(err,RSB_ERRM_ES);
			}

			errval = rsb__mtx_split(mtxWp,manp,&st,&at,&lt,wv,0); /* FIXME: manp is not used yet by rsb__mtx_split */

			if(RSB_SOME_ERROR(errval))
			{
				RSB_PERR_GOTO(err,RSB_ERRM_ES);
			}		
		}
		else
		{
			errval = rsb__leaves_merge(mtxWp,manp,&st,&at,&lt,wv,0);

			if(RSB_SOME_ERROR(errval))
			{
				/* FIXME: this is a critical error and mtxWp shall be deallocated */
				RSB_PERR_GOTO(err,RSB_ERRM_ES);
			}
		}

		tst += st;
		tlt += lt;
		tat += at;
		tmt += rsb_time(); 
		nsap = mtxWp->all_leaf_matrices_n;
		if(wom)
		       	continue; /* want only merging/splitting */
		/* full tuning: */
#if 0
		if(da)
		/*otn = - manp;*/ errval = rsb__tune_spxx(NULL, NULL, &otn, mintimes, 0, 0, RSB_CONST_AT_OP_SAMPLES_MIN, RSB_CONST_AT_OP_SAMPLES_MAX, RSB_AT_TIME_AUTO, transA, alphap, mtxWp, nrhs, order, Bp, ldB, betap, Cp, ldC, op, /*RSB_AUT0_TUNING_SILENT*//*RSB_AUT0_TUNING_VERBOSE*/wv, NULL, NULL, NULL, NULL, NULL );
#endif

		if(RSB_SOME_ERROR(errval))
		{
			RSB_PERR_GOTO(err,RSB_ERRM_ES);
		}

		RSB_ATTR_PREPARE(ttrp,attrp,errval);
		if(RSB_SOME_ERROR(errval))
		{
		       	RSB_PERR_GOTO(err,RSB_ERRM_EM); 
		}

		/* just benchmark: */
		if(db)
			errval = rsb__do_bench_spxm(&tpop,NULL,oasp->transA,oasp->alphap,mtxWp,oasp->nrhs,oasp->order,oasp->Bp,oasp->ldB,oasp->betap,oasp->Cp,oasp->ldC,aiap->maxt,aiap->mintimes,oasp->op,aiap->maxtimes,wv,NULL,&tpos);
		else
		{
#if RSB_TUNE_WITH_LIKWID
			static rsb_bool_t want_likwid = RSB_BOOL_TRUE;
			char likwid_marker_string[2*RSB_MAX_FILENAME_LENGTH];
			RSB_ATR_SAVENAME(likwid_marker_string,errval,aiap->rfpfx,aiap->mtxns,mtxWp->typecode,oasp->transA,oasp->nrhs,/*sval*/1.0,pos,pps,/*ctn*/otn,mtxWp,iot,oasp->op);
			RSB_TM_LIKWID_MARKER_R_START(likwid_marker_string);
#endif /* RSB_TUNE_WITH_LIKWID */
			ctn = 0;
			RSB_ASSIGN_IF_SP(ctn,tnp);
			errval = rsb__do_bench_spxm_t(NULL,&tpop,&ctn,mtxWp,aiap->maxt,aiap->mintimes,aiap->maxtimes,wv,oasp,ttrp,&tpos); /* FIXME: need first argument */
#if RSB_TUNE_WITH_LIKWID
			RSB_TM_LIKWID_MARKER_R_STOP(likwid_marker_string);
#endif /* RSB_TUNE_WITH_LIKWID */
		}
		RSB_ASSERT( ( tpop != RSB_TIME_ZERO ) || RSB_SOME_ERROR(errval) );
		cct += rsb_time();
		RSB_ATTR_ADVANCE(ttrp,attrp,cct); /* FIXME: need to adjust cct time */

		if(RSB_SOME_ERROR(errval))
		{
			RSB_PERR_GOTO(err,RSB_ERRM_ES);
		}

		if(ctn == 0)
		{
			RSB_DO_REINIT_SINGLE_VALUE_GET(RSB_IO_WANT_EXECUTING_THREADS,&ctn,errval);
		}

		if(tpop)
		if(wv)
		{
			RSB_STDOUT("After %s step %d: tpop: %0.4lg s   ~Mflops: %3.3lf   nsubm:%d otn:%d\n",pos,eps+1,tpop,mpoe/tpop,nsap,ctn);
		}

		if( ! ( aiap->continue_partitioning < 0 ) )
		if( mtxOpp )
		{
			if( RSB_TIME_BETTER_THAN(tpop,btpop) || ( aiap->continue_partitioning > 0 ) )
			{
				tct -= rsb_time(); 
				RSB_MTX_FREE(mtxBp);
				mtxBp = mtxWp;
				mtxWp = rsb__mtx_clone_simple(mtxWp);
				otn = ctn;
				if( mtxWp == NULL )
				{
					RSB_STDOUT("Cloning failed: probably out of memory !\n"); /* FIXME: temporary. */
					errval = RSB_ERR_ENOMEM;
					/* A non critical error: execution can continue. */
					RSB_PERR_GOTO(ret,RSB_ERRM_EM);
				}
				tct += rsb_time(); 
				cc ++;
				++pps;
				pps += diar; /* due to the present speedup, we count the degradations and neutral as benign */
				pps += niar;
				diar = 0;
				niar = 0;
				if(wv)
				RSB_STDOUT("Applying %s (%d -> %d leaves, %d th.) yielded SPEEDUP of %6.3lfx: %0.4lgs -> %0.4lgs, so taking this instance.\n",pos,nsdp,nsap,ctn,btpop/tpop,btpop,tpop);
				btpop = tpop;
				btpos = tpos;

				if(RSB_VERBOSE_FOR_SAVING(wv))
					RSB_ATR_SAVE(errval,aiap->rfpfx,aiap->mtxns,mtxWp->typecode,oasp->transA,oasp->nrhs,/*sval*/1.0,pos,pps,/*ctn*/otn,mtxWp,iot,oasp->op);
				/* FIXME: subdivision parameter sval ignored here  ... */
			}
			else /* either niar++ or diar++ */
#if 1
			if( ! RSB_TIME_DEFINITELY_WORSE_THAN(tpop,btpop) )
			{
				++niar;
				if(wv)
					RSB_STDOUT("Applying %s (%d -> %d leaves, %d th.) yielded NEGLIGIBLE change (%dth in a row) (old/new=%.5lfx): %0.4lgs -> %0.4lgs, so IGNORING this instance.\n",pos,nsdp,nsap,ctn,niar,btpop/tpop,btpop,tpop);

			}
#else
			if ( eps == ms-1 ) /* better, minimal threshold if last  */
			{
				diar = 0; // FIXME: here also makes sense to retain last matrix, even if no clone is necessary anymore
				if(wv)
				       	RSB_STDOUT("Negligible improvement: not taking a clone yet.\n");
			}
#endif
		}

		if( RSB_TIME_DEFINITELY_WORSE_THAN(tpop,btpop) || ( aiap->continue_partitioning < 0 ) )
		{
			++diar;
			if(wv) RSB_STDOUT("Applying %s (%d -> %d leaves, %d th.) yielded SLOWDOWN (%dth of %d tolerable) of %6.3lfx: %0.4lgs -> %0.4lgs.\n",pos,nsdp,nsap,ctn,diar,mdiar,tpop/btpop,btpop,tpop);
		}

		if( diar > mdiar || RSB_TIME_MUCH_WORSE_THAN(tpop,btpop) )
		{
			if(wv)
			RSB_STDOUT("Skipping further %s based tests after %d definite performance degradations in a row%s.\n",pos,diar,(RSB_TIME_MUCH_WORSE_THAN(tpop,btpop))?" (and last exceeding limit)":"");
			/* to stop process when no more matrices to merge/split */
			eps ++; /* effective steps */
			goto noloop;
		}

		if ( (!want_split) && nsap == 1 )
		{
			if(wv)
			RSB_STDOUT("Merged all the matrix leaves: no reason to continue merging.\n");
			eps ++; /* effective steps */
			goto noloop;
		}
	} /* eps */
noloop:
	tpt += rsb_time();

	if(RSB_AT_NO_VERBOSE(wv))
		goto ret;

	RSB_STDOUT("A total of %d %s steps (of max %d) (%d -> %d subms) took %0.4lgs (of which %0.4lgs partitioning, %0.4lgs I/O); computing times: %0.4lgs in par. loops, %0.4lgs sorting, %0.4lgs analyzing)\n",eps,pos,ms,nsbp,nsap,tpt,tmt,iot,tlt,tst,tat);

	if(tpop)
	{
		RSB_STDOUT("Total %s + benchmarking process took %0.4lgs, equivalent to %3.1lf/%3.1lf new/old ops (%0.4lgs for %d clones -- as %3.1lf/%3.1lf ops, or %3.1lf/%3.1lf ops per clone)",pos,tpt,tpt/btpop,tpt/otpop,tct,cc,tct/btpop,tct/otpop,(tct/cc)/btpop,(tct/cc)/otpop);
		RSB_STDOUT(", SPEEDUP of %6.3lfx",otpop/btpop);
		if(otpop>btpop)
			RSB_STDOUT("\n");
		else
			RSB_STDOUT(" (NO SPEEDUP)\n");
	}

	if(otpop>btpop)
	{
		RSB_STDOUT("Applying multi-%s (%d -> %d leaves, %d steps, %d -> %d th.sp.) yielded SPEEDUP of %6.3lfx (%0.4lgs -> %0.4lgs), will amortize in %10.1lf ops by saving %0.4lgs per op.\n",pos,mtxAp->all_leaf_matrices_n,mtxBp->all_leaf_matrices_n,pps,itn,otn,otpop/btpop,otpop,btpop,((tpt)/(otpop-btpop)),otpop-btpop);
	}
	goto ret;
err:
	RSB_ERROR("A recoverable error happened during tuning!\n");
ret:
	RSB_MTX_FREE(mtxWp);
        
        if( mtxBp && mtxOpp && *mtxOpp != mtxBp )
	{
	        /* RSB_MTX_FREE(*mtxOpp); */
		if( btpop >= otpop )
		{
			if( aiap->continue_partitioning > 0 )
				btpop = otpop * 0.8; /* for not breaking consistency in callers .. */
			if( aiap->verbose && aiap->continue_partitioning > 0 )
				RSB_STDOUT("Taking clone by force and pretending its time is %5.3lgx the original (%2.3es instead of %2.3es).\n",0.8,btpop,otpop);
		}

	        RSB_ASSIGN_IF(mtxOpp,mtxBp);
	}

	btpop = RSB_MIN(btpop,otpop); /* this is a parachute to signalize autotuner's failure without breaking consistency. */
	RSB_ASSIGN_IF(btpopp,btpop)
	RSB_ASSIGN_IF(btposp,btpos)
	RSB_ASSIGN_IF(otpopp,otpop)
	RSB_ASSIGN_IF(otposp,otpos)
	RSB_ASSIGN_IF(ccp,*ccp+cc)
	RSB_ASSIGN_IF(ctp,*ctp+tct)
	otn = ( otn == 0 ) ? rsb__set_num_threads(RSB_THREADS_GET) /*rsb_want_executing_threads() */ /* rsb_get_num_threads()*/ : otn;
	RSB_ASSIGN_IF(tnp,otn)
	RSB_ASSIGN_IF(epsp,pps)

	return errval;
}

static rsb_err_t rsb__tune_spxx_blk(struct rsb_mtx_t **mtxOpp, rsb_int_t*epsp, rsb_int_t *otcp, const struct rsb_aia_t*aiap, struct rsb_aoa_t*aoap, struct rsb_oas_t*oasp, struct rsb_attr_t *attrp )
{
	/*
	 * Try to obtain a block autotuned clone of a matrix.
	 *
	 * otcp: optimal thread count (pointer)
	 * ...
	 * */
	rsb_err_t errval = RSB_ERR_NO_ERROR;
	const int wv = 1; /* want verbose */
	rsb_int_t ootc = otcp ? *otcp : 0;

        RSB_DEBUG_ASSERT(aiap->mtxAp ? aiap->mtxAp : *mtxOpp);

	if(! aiap->maxms )
	{
		goto wos;
	}

	errval = rsb__rblk_tune_inner(mtxOpp,aiap->maxms,aiap,epsp,&aoap->otpos,&aoap->btpos,&aoap->otpo,otcp,&aoap->btpo,&aoap->cc,&aoap->ct,RSB_BOOL_FALSE,oasp,attrp);

	if(RSB_SOME_ERROR(errval))
	{
		RSB_PERR_GOTO(ret,RSB_ERRM_ES);
	}

	/* if( ( aiap->mtxAp && *mtxOpp ) && ( aiap->mtxAp != *mtxOpp ) ) */
	if( aoap->otpo != aoap->btpo )
	{
		/* autotuning (probably) succeeded */
		goto ret;
	}

	if(wv && aiap->maxss && aiap->verbose)
	{
		RSB_STDOUT("Merging based autotuning FAILED (=NO SPEEDUP); let's try splitting then...\n");
	}

	if( otcp && *otcp != ootc ) /* if 'polluted' *otcp contents... */
		*otcp = ootc; /* ...revert them. TODO: do this within rsb__rblk_tune_inner. */
wos:	/* want only split */
	if(! aiap->maxss )
	{
		goto ret;
	}

	errval = rsb__rblk_tune_inner(mtxOpp,aiap->maxss,aiap,epsp,&aoap->otpos,&aoap->btpos,&aoap->otpo,otcp,&aoap->btpo,&aoap->cc,&aoap->ct,RSB_BOOL_TRUE,oasp,attrp);

	if(RSB_SOME_ERROR(errval))
	{
		RSB_PERR_GOTO(ret,RSB_ERRM_ES);
	}
ret:
	return errval;
}

rsb_err_t rsb__mtx_ms_check(struct rsb_mtx_t ** mtxOpp)
{
	rsb_err_t errval = RSB_ERR_NO_ERROR;
	const rsb_submatrix_idx_t manp = 0;
	rsb_time_t st, at, lt;
	const int wv = 1;
	const int mm = 1;
	const int kc = 1;
	struct rsb_mtx_t *mtxWp = *mtxOpp;
	rsb_long_t nsubm = -1;
	int i;

	for(nsubm=-1, i=0;i<mm && nsubm != rsb__submatrices(mtxWp);++i)
	{
		nsubm = rsb__submatrices(mtxWp);
		errval = rsb__mtx_realloc_with_spare_leaves(&mtxWp, RSB_TMP_OVERALLOC_MTX*rsb__submatrices(mtxWp));
		assert(!errval);
		errval = rsb__mtx_split(mtxWp,manp,&st,&at,&lt,wv,kc); /* FIXME: manp is not used yet by rsb__mtx_split */
		assert(!errval);
		assert(RSB_BOOL_TRUE==rsb__mtx_chk(mtxWp));
	}

	for(nsubm=-1, i=0;i<mm && nsubm != rsb__submatrices(mtxWp);++i)
	{
		nsubm = rsb__submatrices(mtxWp);
		errval = rsb__leaves_merge(mtxWp,manp,&st,&at,&lt,wv,kc);
		assert(!errval);
		assert(RSB_BOOL_TRUE==rsb__mtx_chk(mtxWp));
	}

	for(nsubm=-1, i=0;i<mm && nsubm != rsb__submatrices(mtxWp);++i)
	{
		nsubm = rsb__submatrices(mtxWp);
		errval = rsb__mtx_realloc_with_spare_leaves(&mtxWp, RSB_TMP_OVERALLOC_MTX*rsb__submatrices(mtxWp));
		assert(!errval);
		assert(RSB_BOOL_TRUE==rsb__mtx_chk(mtxWp));
		errval = rsb__mtx_split(mtxWp,manp,&st,&at,&lt,wv,kc); /* FIXME: manp is not used yet by rsb__mtx_split */
		assert(!errval);
		assert(RSB_BOOL_TRUE==rsb__mtx_chk(mtxWp));
	}


	*mtxOpp=mtxWp;
	goto err;
err:
	assert(!errval);
	return errval;
}

static rsb_err_t rsb__tune_spxx_bos(struct rsb_mtx_t ** mtxOpp, rsb_real_t *tsfp, rsb_int_t*epsp, const struct rsb_aia_t*aiap, struct rsb_aoa_t*aoap, struct rsb_oas_t*oasp,  struct rsb_attr_t *attrp)
{
	/*
	 * Original back-end function to librsb-1.1 auto tuning.
	 * Enriched with librsb-1.2 auto tuning functionality.
	 * It uses the auto tuning technique following user's request.
	 * It allocates/deallocates the necessary vectors if these are not provided.
	 */
	rsb_err_t errval = RSB_ERR_NO_ERROR;
	const rsb_time_t tg = RSB_CACHED_TIMER_GRANULARITY;
	rsb_aligned_t alpha[RSB_CONST_ENOUGH_ALIGNED_FOR_ANY_TYPE];
	rsb_aligned_t  beta[RSB_CONST_ENOUGH_ALIGNED_FOR_ANY_TYPE];
	rsb_real_t tsf = 1.0, sfm = 1.0; /* total speed-up factor / speed up factor multiplier */
	rsb_int_t ri, gr = RSB_MIN(2,aiap->maxr); /* round index / grace rounds */
	const rsb_int_t rc = aiap->maxr; /* round count */
	double ospv = RSB_REAL_ZERO; /* optimal subdivision parameter value */
	void *tBp = NULL, *tCp = NULL; /* temporary  B / C  pointers */
	const struct rsb_mtx_t * mtxBp = NULL; /* Base matrix */
	rsb_time_t tatt = RSB_TIME_ZERO, tiot = RSB_TIME_ZERO; /* total  auto tuning / i.o. time */
	rsb_time_t tct = RSB_TIME_ZERO; /* total constructor time */
	rsb_time_t fotpo = RSB_TIME_ZERO; /* first original time per operation */
	rsb_int_t mt = rsb__set_num_threads(RSB_THREADS_GET_MAX_SYS);
	rsb_int_t tn = RSB_DT_SPEC_THREADS_TNP(aiap->tnp) ? *aiap->tnp : 0;
	const rsb_bool_t only_one_sample = ( (mtxOpp == NULL) && RSB_DT_SAME_THREADS_TNP(aiap->tnp) ) || ( aiap->maxr == 0 );
	struct rsb_aoa_t aoa; /* autotuning output args */
#if RSB_ALLOW_INTERNAL_GETENVS
	if(getenv("RSB_AT_GR" )) gr  = rsb__util_atoi(getenv("RSB_AT_GR" )); /* override: grace rounds */
#endif /* RSB_ALLOW_INTERNAL_GETENVS */

	RSB_BZERO_P(&aoa);
	RSB_THREAD_STATS

	rsb__attr_init(attrp,aiap->verbose>1?RSB_BOOL_TRUE:RSB_BOOL_FALSE,aiap->mtxAp,aiap->mtxns,oasp->transA,oasp->nrhs,oasp->op);

	if( mtxOpp && *mtxOpp )
		mtxBp = *mtxOpp;
	else
		mtxBp = aiap->mtxAp;

	if( RSB_DT_SPEC_THREADS_TNP(aiap->tnp) && *aiap->tnp > mt )
	{
		/* if max thread count is exceeded, the code will crash (this shall be fixed for resilience purposes) */
		errval = RSB_ERR_BADARGS;
	       	RSB_PERR_GOTO(err,"Requested an excessive threads count !");
	}
	
	if( ( mtxOpp != NULL ) && ( *mtxOpp != NULL ) && ( aiap->mtxAp != NULL ) )
	{
		errval = RSB_ERR_BADARGS;
	       	RSB_PERR_GOTO(err,RSB_ERRM_ES);
	}
	
	if( mtxBp == NULL )
	{
		errval = RSB_ERR_BADARGS;
	       	RSB_PERR_GOTO(err,RSB_ERRM_ES);
	}
	
	if( ( mtxOpp == NULL ) && ( aiap->mtxAp != NULL ) && ( aiap->tnp == NULL) )
	{
		errval = RSB_ERR_BADARGS;
	}
	
/*
	if( tnp == NULL && mtxOpp == NULL )
	{
		errval = RSB_ERR_BADARGS;
	       	RSB_PERR_GOTO(err,RSB_ERRM_ES);
	}
*/
	if(/*RSB_VERBOSE_FOR_MSG*/(aiap->verbose))
	{
		if(only_one_sample)
			RSB_STDOUT("Will use autotuning routine to sample matrix: ");
		else
			RSB_STDOUT("Will autotune matrix: ");
		RSB_STDOUT(RSB_PRINTF_MATRIX_AT_SUMMARY_ARGS(mtxBp));
		RSB_STDOUT(".\n");
		RSB_STDOUT("Parameters: verbosity:%d mintimes:%d maxtimes:%d mindt:%.4lg maxdt:%.4lg\n",(int)aiap->verbose,(int)aiap->mintimes,(int)aiap->maxtimes,(double)0 /*aiap->mindt*/,(double)aiap->maxt);
	}

	if(aiap->tnp && mtxOpp == NULL && mtxBp->all_leaf_matrices_n<2 )
	{
		/* FIXME: special case */
		if(RSB_VERBOSE_FOR_MSG(aiap->verbose))
			RSB_STDOUT("Selected threads tuning but matrix has only %d leaves. Skipping.\n",mtxBp->all_leaf_matrices_n);
	       	/* RSB_PERR_GOTO(err,RSB_ERRM_ES); */
	}

	if( oasp->alphap == NULL )
		oasp->alphap = &alpha[0];
	if( oasp->betap == NULL )
		oasp->betap = &beta[0];
	if( oasp->nrhs < 1 )
		oasp->nrhs = 1;

	if( oasp->Bp == NULL )
	{
		const rsb_nnz_idx_t ndB = ( oasp->order == RSB_FLAG_WANT_COLUMN_MAJOR_ORDER ) ? oasp->nrhs : RSB_MTX_TRANSPOSED_COLS(mtxBp,oasp->transA);

		oasp->Bp = tBp = rsb__calloca_vector(oasp->ldB*ndB,mtxBp->typecode);
		if( !tBp )
		{
			errval = RSB_ERR_INTERNAL_ERROR;
		       	RSB_PERR_GOTO(err,RSB_ERRM_ES);
		}
	}

	if( oasp->Cp == NULL )
	{
		const rsb_nnz_idx_t ndC = ( oasp->order == RSB_FLAG_WANT_COLUMN_MAJOR_ORDER ) ? oasp->nrhs : RSB_MTX_TRANSPOSED_ROWS(mtxBp,oasp->transA);

		oasp->Cp = tCp = rsb__calloca_vector(oasp->ldC*ndC,mtxBp->typecode);
		if( !tCp )
		{
			errval = RSB_ERR_INTERNAL_ERROR;
			RSB_PERR_GOTO(err,RSB_ERRM_ES);
		}
	}

	if(aiap->verbose)
	if(aiap->maxr > 1) /* FIXME: this is for the cases the routine is invoked from the merge/split based one. */
	{
		rsb_int_t ctn = 0.0;
		RSB_DO_REINIT_SINGLE_VALUE_GET(RSB_IO_WANT_EXECUTING_THREADS,&ctn,errval);
		if( only_one_sample )
			RSB_STDOUT("Sampling");
		else
			RSB_STDOUT("Starting autotuning");
		RSB_STDOUT(" (%zd x %lg s stages, transA=%c, nrhs=%zd, timer gran.=%lg), %zd suggested as starting thread count%s.\n",(rsb_printf_int_t)aiap->maxr,aiap->maxt,RSB_TRANSPOSITION_AS_CHAR(oasp->transA),(rsb_printf_int_t)oasp->nrhs,tg,(rsb_printf_int_t)ctn,(tn==0)?"(default)":"");
	}

	if( mtxOpp )
	if(RSB_VERBOSE_FOR_SAVING(aiap->verbose))
		RSB_ATR_SAVE(errval,aiap->rfpfx,aiap->mtxns,mtxBp->typecode,oasp->transA,oasp->nrhs,ospv,"r",0,0,mtxBp,tiot,oasp->op);

	if(RSB_SOME_ERROR(errval = rsb__fill_with_ones(beta ,mtxBp->typecode,1,1))){ RSB_PERR_GOTO(err,RSB_ERRM_ES)}
	if(RSB_SOME_ERROR(errval = rsb__fill_with_ones(alpha,mtxBp->typecode,1,1))){ RSB_PERR_GOTO(err,RSB_ERRM_ES)}

	if(aiap->maxr == 0)
	{
		tsf = 1.0; /* use this routine as benchmark only: no speedup */
		if(/*RSB_VERBOSE_FOR_MSG*/(aiap->verbose))
			RSB_STDOUT("Taking a single sample of performance (zero autotuning rounds requested).\n");
		errval = rsb__do_bench_spxm_t(NULL,&fotpo,NULL,mtxBp,aiap->maxt,aiap->mintimes,aiap->maxtimes,aiap->verbose,oasp,NULL/* TODO: attrp/tattrp */,&aoa.otpos);
		aoa.otpo = aoa.btpo = fotpo;
		aoa.btpos = aoa.otpos;
		RSB_ASSIGN_IF(aoap,aoa)
		goto err;
	}

	RSB_DO_REINIT_SINGLE_VALUE_GET(RSB_IO_WANT_SUBDIVISION_MULTIPLIER,&ospv,errval);
	for(ri=0;ri<rc;++ri)
	{
		rsb_time_t ratt = RSB_TIME_ZERO, riot = RSB_TIME_ZERO, rct = RSB_TIME_ZERO; /* round  auto tuning / i.o. / constructor  time */

		if( aiap->maxms > 0 || aiap->maxss > 0 )
		{
			struct rsb_mtx_t * mtxOp = mtxOpp ? *mtxOpp : NULL;
			rsb_int_t eps = 0; /* effective partitiong steps */

			tatt = - rsb_time();
		        RSB_DEBUG_ASSERT(
				( (aiap->mtxAp != NULL) && (mtxOpp && *mtxOpp == NULL ) )||
				( (aiap->mtxAp == NULL) && (mtxOpp && *mtxOpp != NULL ) )
			);
			errval = rsb__tune_spxx_blk(mtxOpp,&eps,aiap->tnp,aiap,&aoa,oasp,attrp);
			fotpo = aoa.otpo;
#if RSB_AT_DESTROYS_MTX
			if(mtxOpp && *mtxOpp != mtxOp)
			{
		 		RSB_MTX_FREE(mtxOp);
			}
#endif /* RSB_AT_DESTROYS_MTX */
			RSB_ASSIGN_IF(epsp,eps)
			RSB_ASSIGN_IF(aoap,aoa)
			tatt += rsb_time();
			tct += aoa.ct;
			goto done; /* no reason to continue */
		}

		/* Old style tuning. */

		errval = rsb__do_tune_spxm_round(&sfm,mtxOpp,&ospv,aiap,&riot,&ratt,&rct,&aoa,oasp,attrp);
		tatt += ratt;
		tiot += riot;
		tct += rct;

		if( fotpo == RSB_TIME_ZERO ) /* first iteration's baseline */
		{
			aoap->otpos = aoa.otpos;
			fotpo = aoa.otpo;
		}
		aoap->btpos = aoa.btpos;

		if(RSB_SOME_ERROR(errval)){ RSB_PERR_GOTO(err,RSB_ERRM_EM); }

		tsf *= sfm;

		if( mtxOpp && ( mtxBp != *mtxOpp ) )
		{
			if( *mtxOpp )
			if( RSB_VERBOSE_FOR_MSG(aiap->verbose) )
			{
				RSB_STDOUT("Found a better candidate clone matrix (original is gone): ");
				RSB_STDOUT(RSB_PRINTF_MATRIX_AT_SUMMARY_ARGS(*mtxOpp));
				RSB_STDOUT("\n ");
			}
			mtxBp = *mtxOpp;
		}

		/* if(rc>1) */
		if( !only_one_sample )
		if(aiap->verbose)
		{
			RSB_STDOUT("Last tuner inner round (%d of %d) took %lg s (eq. to %6.1lg/%6.1lg old/new op.times), gained local/global speedup %lg x (%lg : %lg) / %lg x (%lg : %lg). ",
					ri+1,rc,ratt,ratt/fotpo,ratt/aoa.btpo,/*tsf*/ aoa.otpo/aoa.btpo,aoa.otpo,aoa.btpo, fotpo/aoa.btpo,fotpo,aoa.btpo);
			if(aoa.btpo < aoa.otpo)
				RSB_STDOUT("This is amortizable in %zd op.times.\n", (size_t)ceil((ratt)/(aoa.otpo-aoa.btpo)));
			else
				RSB_STDOUT("This is not amortizable !\n");
		}

		if( ( ( ( sfm <= 1.0) && mtxOpp != NULL && !( sfm > 1.0) ) || ( mtxOpp == NULL && ( aiap->tnp == NULL || tn == *aiap->tnp ) ) ) )
		{
			if( !only_one_sample )
			if(aiap->verbose)
			{
				RSB_STDOUT("Auto tuning inner round %d did not find a configuration better than the original.\n",ri+1);
			}
			if(gr > 1)
			{
				gr--;
				if( !only_one_sample )
				if(aiap->verbose)
					RSB_STDOUT("Will attempt %d more 'grace' rounds.\n",gr);
			}
			else
			{
				if( !only_one_sample )
				if(aiap->verbose)
				if(ri+1<rc)
					RSB_STDOUT("Skipping further auto tuning.\n");
				goto done;
			}
		}
		/* if(!mtxOpp)
			tnpp = *tnp; */
	} /* ri */
done:
	tsf = fotpo/aoa.btpo; /* FIXME */

	if( ! only_one_sample )
	if(aiap->verbose)
	{
		RSB_STDOUT("In %d tuning rounds (tot. %0.2lgs, %0.2lgs for constructor, %d clones) obtained ",RSB_MIN(ri+1,rc),tatt,tct,aoa.cc);
		if(fotpo==aoa.btpo)
	       		RSB_STDOUT("NO speedup (best stays %0.4lg Mflops).",aiap->ofe/fotpo);
		else
		{
	      		RSB_STDOUT("%s %6.1lf%% (%0.4lgx) (from %0.4lg to %0.4lg Mflops).",
				(fotpo < aoa.btpo) ? "a dubious SLOWDOWN of " : "a SPEEDUP of",
				RSB_SPEEDUP_TO_PCT(fotpo/aoa.btpo),fotpo/aoa.btpo,aiap->ofe/fotpo,aiap->ofe/aoa.btpo);
		}
		if(tiot)
			RSB_STDOUT(" Employed %0.2lgs for I/O of matrix plots.",tiot);
		RSB_STDOUT("\n");
	}
err:
	RSB_THREAD_STATS
	if(!RSB_SOME_ERROR(errval))
		RSB_ASSIGN_IF(tsfp,tsf)
	RSB_CONDITIONAL_FREE_ALLOCA(tBp);
	RSB_CONDITIONAL_FREE_ALLOCA(tCp);

	return errval;
}

rsb_err_t rsb__tune_spxx( struct rsb_mtx_t ** mtxOpp, rsb_real_t *tsfp, rsb_int_t *tnp, rsb_int_t maxr, rsb_int_t maxms, rsb_int_t maxss, rsb_int_t mintimes, rsb_int_t maxtimes, rsb_time_t maxt, rsb_trans_t transA, const void * alphap, const struct rsb_mtx_t * mtxAp, rsb_coo_idx_t nrhs, rsb_flags_t order, const void * Bp, rsb_nnz_idx_t ldB, const void * betap, void * Cp, rsb_nnz_idx_t ldC, enum rsb_op_t op, rsb_int_t*epsp, rsb_time_t*otpopp, rsb_time_t*btpopp, int verbose, const char * fprfn, const char * mtxns, struct rsb_attr_t *attrp, struct rsb_ts_t*otposp, struct rsb_ts_t*btposp )
{
	/*
		otpopp = original time per operation pointer
		btpopp =   best   time per operation pointer
	  
	 	Note: it is possible that *mtxOpp will be overwritten and an error code returned.
	 */
	rsb_err_t errval = RSB_ERR_BADARGS;
        
	if( mtxAp != NULL && mtxOpp != NULL && *mtxOpp != NULL )
	{
		errval = RSB_ERR_BADARGS;
		RSB_PERR_GOTO(ret,RSB_ERRM_ES);
	}
	RSB_ASSERT( ! ( mtxAp != NULL && mtxOpp != NULL && *mtxOpp != NULL ) );

        if(!mtxOpp || !*mtxOpp)
        {
		if(verbose && ( maxms || maxss ) )
			RSB_STDOUT("Setting old tuning style.");
                maxms = maxss = 0; /* merge / split cannot apply in this case */
        }

#if RSB_ALLOW_INTERNAL_GETENVS
	if(getenv("RSB_MAXMS"))
		maxms = rsb__util_atoi(getenv("RSB_MAXMS")); /* override: max merge steps */
	if(getenv("RSB_MAXSS"))
		maxss = rsb__util_atoi(getenv("RSB_MAXSS")); /* override: max subdivide steps */
	if(getenv("RSB_MAXR" ))
		maxr  = rsb__util_atoi(getenv("RSB_MAXR" )); /* override: max rounds */
#endif /* RSB_ALLOW_INTERNAL_GETENVS */

        if( mtxAp || ( mtxOpp && *mtxOpp ) )
        {
		/*rsb_int_t mintimes = RSB_CONST_AT_OP_SAMPLES_MIN ;*/ /* min iterations for taking operation time samples */
		/*rsb_int_t maxtimes = RSB_CONST_AT_OP_SAMPLES_MAX ;*/ /* max iterations for taking operation time samples */
		const rsb_time_t tg = RSB_CACHED_TIMER_GRANULARITY;
        	struct rsb_oas_t oas = { op, transA, alphap, nrhs, order, Bp, ldB, betap, Cp, ldC }; /* FIXME: these may be overwritten */
        	struct rsb_aia_t aia = { maxt, mintimes, maxtimes, fprfn, mtxns, RSB_PERF_ZERO, verbose, tnp, mtxAp, maxr, maxms, maxss, 0 /* =cont.part. */ };
        	struct rsb_aoa_t aoa;
		const struct rsb_mtx_t * mtxp = mtxAp ? mtxAp : *mtxOpp /*aia.mtxAp*/;

		if(RSB_SOME_ERROR(errval = rsb__set_ldX_for_spmm(transA, mtxp, nrhs, order, &oas.ldB, &oas.ldC)))
			RSB_PERR_GOTO(ret,RSB_ERRM_EM);

		RSB_BZERO_P(&aoa);
		aia.ofe = rsb__estimate_mflops_per_op_spmv_uaua(mtxp) * nrhs;
	        /* maxt = ( maxt <= RSB_TIME_ZERO ) ? RSB_AUT0_TUNING_DEFAULT_TIME : maxt;*/

	        if ( aia.maxt < RSB_TIME_ZERO )
			aia.mintimes = - ceil( aia.maxt );
        	aia.maxt = ( aia.maxt <= RSB_TIME_ZERO ) ? tg : aia.maxt; /* FIXME: this code seems obscure */

        	aia.mtxns = ( aia.mtxns == NULL ) ? "rsb" : aia.mtxns; /* matrix name string */
        	aia.verbose = RSB_MAX(rsb_global_session_handle.verbose_tuning,aia.verbose);
#if RSB_ALLOW_INTERNAL_GETENVS
		if(getenv("RSB_VERBOSE_TUNING"))
			aia.verbose = rsb__util_atoi(getenv("RSB_VERBOSE_TUNING"));
#endif /* RSB_ALLOW_INTERNAL_GETENVS*/
		if( aia.maxr > RSB_CONST_MAX_TUNING_ROUNDS )
		{
#if 1
			if(RSB_VERBOSE_FOR_MSG(aia.verbose))
				RSB_STDOUT("The specified %d tuning rounds is exaggerated: forcing down to %d ;)\n",aia.maxr,RSB_CONST_MAX_TUNING_ROUNDS);
			aia.maxr = RSB_MIN(aia.maxr,RSB_CONST_MAX_TUNING_ROUNDS);
#else
			errval = RSB_ERR_INTERNAL_ERROR;
		       	RSB_PERR_GOTO(err,RSB_ERRM_IE);
#endif
		}
#if RSB_ALLOW_INTERNAL_GETENVS
		if(getenv("RSB_MINTIMES"))
			aia.mintimes = rsb__util_atoi(getenv("RSB_MINTIMES"));
		if(getenv("RSB_CONTINUE_PARTITIONING"))
			aia.continue_partitioning = rsb__util_atoi(getenv("RSB_CONTINUE_PARTITIONING"));
		if(getenv("RSB_MAXTIMES"))
			aia.maxtimes = rsb__util_atoi(getenv("RSB_MAXTIMES"));
#endif /* RSB_ALLOW_INTERNAL_GETENVS*/

#if RSB_ALLOW_INTERNAL_GETENVS
		if(rsb__getenv_int_t("RSB_AT_CHECKS",0)==1)
        	if( mtxOpp )
		{
        		errval = rsb__mtx_ms_check( mtxOpp );
			if(RSB_SOME_ERROR(errval))
				RSB_PERR_GOTO(ret,RSB_ERRM_EM);
			goto ret;
		}
#endif /* RSB_ALLOW_INTERNAL_GETENVS*/
        	errval = rsb__tune_spxx_bos( mtxOpp, tsfp, epsp, &aia, &aoa, &oas, attrp );

#if 0
		//if(RSB_VERBOSE_FOR_MSG(aia.verbose))
		if(aia.verbose)
		{
			RSB_STAT_DUMP_TS(aoa.btpos);
			RSB_STAT_DUMP_TS(aoa.otpos);
		}
#endif
        	RSB_ASSIGN_IF(btposp,aoa.btpos)
        	RSB_ASSIGN_IF(otposp,aoa.otpos)
        	RSB_ASSIGN_IF(btpopp,aoa.btpo)
        	RSB_ASSIGN_IF(otpopp,aoa.otpo)
        }
        else
                goto ret;
ret:
	return errval;
}

#define RSB__ATTRP_FOR_SPXM_TUNE ( rsb_global_session_handle.verbose_tuning > 2 ? rsb__calloc(sizeof(struct rsb_attr_t)) : NULL )

rsb_err_t rsb__do_tune_spmm(struct rsb_mtx_t ** mtxOpp, rsb_real_t *sfp, rsb_int_t *tnp, rsb_int_t maxr, rsb_time_t maxt, rsb_trans_t transA, const void * alphap, const struct rsb_mtx_t * mtxAp, rsb_coo_idx_t nrhs, rsb_flags_t order, const void * Bp, rsb_nnz_idx_t ldB, const void * betap, void * Cp, rsb_nnz_idx_t ldC)
{
	rsb_err_t errval = RSB_ERR_BADARGS;
	struct rsb_attr_t * attrp = RSB__ATTRP_FOR_SPXM_TUNE;

	errval = rsb__tune_spxx( mtxOpp, sfp, tnp, RSB_MAX(maxr,1), RSB_CONST_DEF_MS_AT_AUTO_STEPS, RSB_CONST_DEF_MS_AT_AUTO_STEPS, RSB_CONST_AT_OP_SAMPLES_MIN, RSB_CONST_AT_OP_SAMPLES_MAX, maxt, transA, alphap, mtxAp, nrhs, order, Bp, ldB, betap, Cp, ldC, rsb_op_spmv, NULL, NULL, NULL, RSB_AUT0_TUNING_SILENT, NULL, NULL, attrp, NULL, NULL);
	rsb__attr_dump(attrp);
	RSB_CONDITIONAL_FREE(attrp);
	return errval;
}

rsb_err_t rsb__do_tune_spsm(struct rsb_mtx_t ** mtxOpp, rsb_real_t *sfp, rsb_int_t *tnp, rsb_int_t maxr, rsb_time_t maxt, rsb_trans_t transA, const void * alphap, const struct rsb_mtx_t * mtxAp, rsb_coo_idx_t nrhs, rsb_flags_t order, const void * Bp, rsb_nnz_idx_t ldB, const void * betap, void * Cp, rsb_nnz_idx_t ldC)
{
	rsb_err_t errval = RSB_ERR_BADARGS;
	struct rsb_attr_t * attrp = RSB__ATTRP_FOR_SPXM_TUNE;

	errval = rsb__tune_spxx( mtxOpp, sfp, tnp, RSB_MAX(maxr,1), RSB_CONST_DEF_MS_AT_AUTO_STEPS, RSB_CONST_DEF_MS_AT_AUTO_STEPS, RSB_CONST_AT_OP_SAMPLES_MIN, RSB_CONST_AT_OP_SAMPLES_MAX , maxt, transA, alphap, mtxAp, nrhs, order, Bp, ldB, betap, Cp, ldC, rsb_op_spsvlt, NULL, NULL, NULL, RSB_AUT0_TUNING_SILENT, NULL, NULL, attrp, NULL, NULL);
	rsb__attr_dump(attrp);
	RSB_CONDITIONAL_FREE(attrp);
	return errval;
}

/* @endcond */