diff --git a/Makefile b/Makefile
index 803ab991aa7dd7d172da46c1bb10650b484cc58c..3a7f943df664ce6aa46ab5052fcfd89e73a56649 100644
--- a/Makefile
+++ b/Makefile
@@ -3,7 +3,7 @@ CFLAGS = -Iinclude -Wall -g3
LDFLAGS = -lm
LIBHQR=src/libhqr.a
-default: testing_pivgen testing_treedraw
+default: testing_pivgen testing_treedraw testing_tileinit
src/queue.o: include/queue.h include/common.h
src/libhqr.o: include/common.h include/libhqr.h
@@ -24,6 +24,9 @@ testing_pivgen : testings/testing_pivgen.o ${LIBHQR}
testing_treedraw : testings/testing_treedraw.o ${LIBHQR}
$(CC) -o $@ $^ $(CFLAGS) $(LDFLAGS)
+testing_tileinit : testings/testing_tileinit.o ${LIBHQR}
+ $(CC) -o $@ $^ $(CFLAGS) $(LDFLAGS)
+
clean :
rm -f include/*~
rm -f src/*.o src/*~
@@ -31,4 +34,4 @@ clean :
cleanall: clean
rm -f ${LIBHQR}
- rm -f testing_pivgen testing_treedraw tree.svg
+ rm -f testing_pivgen testing_treedraw testing_tileinit tree.svg
diff --git a/include/libhqr.h b/include/libhqr.h
index 9e8baa03fdacdcab59a0fc56461a0bf2bbed2a8c..5be6d282152952ceeda0d5b6f1433d3aff2fd8f6 100644
--- a/include/libhqr.h
+++ b/include/libhqr.h
@@ -89,6 +89,19 @@ typedef struct libhqr_tiledesc_s{
int p; /**< The number of nodes per column in the data distribution */
} libhqr_tiledesc_t;
+
+/**
+ * @brief Minimal structure to stock the information for each tile
+ */
+typedef struct libhqr_tileinfo_s{
+ int type; /**< The type of the tile */
+ int currpiv; /**< Number of the time annihilating the current tile */
+ int nextpiv; /**< The next tile killed by the tile */
+ int prevpiv; /**< The previous tile killed by the tile */
+ int first_nextpiv; /**< The first next tile killed */
+ int first_prevpiv; /**< The first previous tile killed */
+} libhqr_tileinfo_t;
+
struct libhqr_tree_s;
typedef struct libhqr_tree_s libhqr_tree_t;
@@ -172,6 +185,7 @@ struct libhqr_tree_s {
int a;
/** Size of highest level tree (distributed one) */
int p;
+ int domino;
void *args;
};
@@ -190,6 +204,14 @@ int libhqr_hqr_init( libhqr_tree_t *qrtree,
int a, int p, int domino, int tsrr );
void libhqr_hqr_finalize( libhqr_tree_t *qrtree );
+
+void libhqr_matrix_init(libhqr_tree_t *qrtree, const libhqr_tree_t *qrtree_init);
+int rdmtx_gettype(const libhqr_tree_t *qrtree, int k, int m);
+int rdmtx_currpiv(const libhqr_tree_t *qrtree, int k, int m);
+int rdmtx_nextpiv(const libhqr_tree_t *qrtree, int k, int p, int m);
+int rdmtx_prevpiv(const libhqr_tree_t *qrtree, int k, int p, int m);
+void libhqr_matrix_finalize(libhqr_tree_t *qrtree);
+
/*
* function for drawing the tree
*/
diff --git a/src/libhqr.c b/src/libhqr.c
index 434bf2fb4e9cdca5237def9ce60e89086a119cdc..2750739839116e72aefb4d3cf31563e6cb0344a4 100644
--- a/src/libhqr.c
+++ b/src/libhqr.c
@@ -170,6 +170,73 @@ static void hqr_low_greedy_init( hqr_subpiv_t *arg, int minMN);
static void hqr_low_binary_init( hqr_subpiv_t *arg);
static void hqr_low_fibonacci_init(hqr_subpiv_t *arg, int minMN);
+
+/****************************************************
+ * Reading functions
+ **************************************************
+ *
+ * Common parameters for the following functions
+ * qrtree - tree used for the factorization
+ * k - Step k of the QR factorization
+ * m - line anhilated
+ */
+
+int rdmtx_gettype(const libhqr_tree_t *qrtree, int k, int m){
+ int id;
+ libhqr_tileinfo_t *arg = (libhqr_tileinfo_t*)(qrtree->args);
+ id = (k * qrtree->mt) + m;
+ return arg[id].type;
+}
+
+int rdmtx_currpiv(const libhqr_tree_t *qrtree, int k, int m){
+ int id, perm_m, p, a;
+ libhqr_tileinfo_t *arg = (libhqr_tileinfo_t*)(qrtree->args);
+ perm_m = m;
+ p = qrtree->p;
+ a = qrtree->a;
+ myassert( (p==1) || (perm_m / (p*a)) == (m / (p*a)) );
+ myassert( (p==1) || (perm_m % p) == (m % p) );
+ id = (k * qrtree->mt) + m;
+ return arg[id].currpiv;
+}
+
+/*
+ * Extra parameter:
+ * p - line used as pivot
+ */
+
+int rdmtx_nextpiv(const libhqr_tree_t *qrtree, int k, int p, int m){
+ int id;
+ libhqr_tileinfo_t *arg = (libhqr_tileinfo_t*)(qrtree->args);
+ int gmt = qrtree->mt;
+ myassert( m > p && p >= k );
+ myassert( m == gmt || p == rdmtx_currpiv( qrtree, k, m ) );
+ if(m == qrtree->mt){
+ id = (k * qrtree->mt) + p;
+ return arg[id].first_nextpiv;
+ }
+ else{
+ id = (k * qrtree->mt) + m;
+ return arg[id].nextpiv;
+ }
+}
+
+int rdmtx_prevpiv(const libhqr_tree_t *qrtree, int k, int p, int m){
+ int id;
+ libhqr_tileinfo_t *arg = (libhqr_tileinfo_t*)(qrtree->args);
+ int gmt = qrtree->mt;
+ myassert( m >= p && p >= k && m < gmt );
+ myassert( m == p || p == rdmtx_currpiv( qrtree, k, m ) );
+
+ if(m == p){
+ id = (k * qrtree->mt) + p;
+ return arg[id].first_prevpiv;
+ }
+ else{
+ id = (k * qrtree->mt) + m;
+ return arg[id].prevpiv;
+ }
+}
/****************************************************
* Common ipiv
***************************************************
@@ -189,11 +256,10 @@ static void hqr_low_fibonacci_init(hqr_subpiv_t *arg, int minMN);
* The number of geqrt to execute in the panel k
*/
static int hqr_getnbgeqrf( const libhqr_tree_t *qrtree, int k ) {
- hqr_args_t *arg = (hqr_args_t*)(qrtree->args);
int a = qrtree->a;
int p = qrtree->p;
int gmt = qrtree->mt;
- int domino = arg->domino;
+ int domino = qrtree->domino;
int pa = p * a;
int nb_1, nb_2, nb_3;
int nb_11, nb_12;
@@ -236,12 +302,9 @@ static int hqr_getnbgeqrf( const libhqr_tree_t *qrtree, int k ) {
*/
static int hqr_getm( const libhqr_tree_t *qrtree, int k, int i )
{
- hqr_args_t *arg = (hqr_args_t*)(qrtree->args);
int a = qrtree->a;
int p = qrtree->p;
- int gmt = qrtree->mt + 1;
- int domino = arg->domino;
- int *perm = arg->perm + gmt * k;
+ int domino = qrtree->domino;
int pos1, j, pa = p * a;
int nbextra1 = nbextra1_formula;
@@ -258,7 +321,7 @@ static int hqr_getm( const libhqr_tree_t *qrtree, int k, int i )
pos1 = ( ( (p * (k+1)) + pa - 1 ) / pa ) * pa;
else
pos1 = ( ( (p + k ) + pa - 1 ) / pa ) * pa;
- return perm[ pos1 + (j/p) * pa + j%p ];
+ return pos1 + (j/p) * pa + j%p;
}
}
@@ -270,11 +333,10 @@ static int hqr_getm( const libhqr_tree_t *qrtree, int k, int i )
*/
static int hqr_geti( const libhqr_tree_t *qrtree, int k, int m )
{
- hqr_args_t *arg = (hqr_args_t*)(qrtree->args);
int a = qrtree->a;
int p = qrtree->p;
- int domino = arg->domino;
- int lm = hqr_getinvperm( qrtree, k, m );
+ int domino = qrtree->domino;
+ int lm = m;
int pos1, j, pa = p * a;
int nbextra1 = nbextra1_formula;
@@ -306,12 +368,11 @@ static int hqr_geti( const libhqr_tree_t *qrtree, int k, int m )
* 3 - if m is reduced in distributed
*/
static int hqr_gettype( const libhqr_tree_t *qrtree, int k, int m ) {
- hqr_args_t *arg = (hqr_args_t*)(qrtree->args);
int a = qrtree->a;
int p = qrtree->p;
- int domino = arg->domino;
+ int domino = qrtree->domino;
- int lm = hqr_getinvperm( qrtree, k, m );
+ int lm = m;
myassert( lm >= k );
/* Element to be reduce in distributed */
@@ -1267,20 +1328,15 @@ static void hqr_high_greedy_init(hqr_subpiv_t *arg, int minMN){
static int hqr_currpiv(const libhqr_tree_t *qrtree, int k, int m)
{
hqr_args_t *arg = (hqr_args_t*)(qrtree->args);
- int tmp, tmpk, perm_m;
- int lm, rank, *perm;
+ int tmp, tmpk;
+ int lm, rank;
int a = qrtree->a;
int p = qrtree->p;
- int domino = arg->domino;
+ int domino = qrtree->domino;
int gmt = qrtree->mt;
- perm_m = hqr_getinvperm( qrtree, k, m );
- lm = perm_m / p; /* Local index in the distribution over p domains */
- rank = perm_m % p; /* Staring index in this distribution */
- perm = arg->perm + (gmt+1) * k;
-
- myassert( (p==1) || (perm_m / (p*a)) == (m / (p*a)) );
- myassert( (p==1) || (perm_m % p) == (m % p) );
+ lm = m / p; /* Local index in the distribution over p domains */
+ rank = m % p; /* Staring index in this distribution */
/* TS level common to every case */
if ( domino ) {
@@ -1289,20 +1345,20 @@ static int hqr_currpiv(const libhqr_tree_t *qrtree, int k, int m)
case 0:
tmp = lm / a;
if ( tmp == k / a )
- return perm[ k * p + rank ]; /* Below to the first bloc including the diagonal */
+ return k * p + rank ; /* Below to the first bloc including the diagonal */
else
- return perm[ tmp * a * p + rank ];
+ return tmp * a * p + rank ;
break;
case 1:
- tmp = arg->llvl->currpiv(arg->llvl, k, perm_m);
- return perm[ ( tmp == k / a ) ? k * p + rank : tmp * a * p + rank ];
+ tmp = arg->llvl->currpiv(arg->llvl, k, m);
+ return ( tmp == k / a ) ? k * p + rank : tmp * a * p + rank;
break;
case 2:
return m - p;
break;
case 3:
if ( arg->hlvl != NULL )
- return arg->hlvl->currpiv(arg->hlvl, k, perm_m);
+ return arg->hlvl->currpiv(arg->hlvl, k, m);
default:
return gmt;
}
@@ -1314,20 +1370,20 @@ static int hqr_currpiv(const libhqr_tree_t *qrtree, int k, int m)
tmp = lm / a;
/* tmpk = (k + p - 1 - m%p) / p / a; */
tmpk = k / (p * a);
- return perm[ ( tmp == tmpk ) ? k + (perm_m-k)%p : tmp * a * p + rank ];
+ return ( tmp == tmpk ) ? k + (m-k)%p : tmp * a * p + rank ;
break;
case 1:
- tmp = arg->llvl->currpiv(arg->llvl, k, perm_m);
+ tmp = arg->llvl->currpiv(arg->llvl, k, m);
/* tmpk = (k + p - 1 - m%p) / p / a; */
tmpk = k / (p * a);
- return perm[ ( tmp == tmpk ) ? k + (perm_m-k)%p : tmp * a * p + rank ];
+ return ( tmp == tmpk ) ? k + (m-k)%p : tmp * a * p + rank ;
break;
case 2:
- return perm[ perm_m - p];
+ return m - p;
break;
case 3:
if ( arg->hlvl != NULL )
- return perm[arg->hlvl->currpiv(arg->hlvl, k, perm_m)];
+ return arg->hlvl->currpiv(arg->hlvl, k, m);
default:
return gmt;
}
@@ -1357,31 +1413,23 @@ static int hqr_nextpiv(const libhqr_tree_t *qrtree, int k, int pivot, int start)
{
hqr_args_t *arg = (hqr_args_t*)(qrtree->args);
int tmp, ls, lp, nextp;
- int opivot, ostart; /* original values before permutation */
- int lpivot, rpivot, lstart, *perm;
+ int lpivot, rpivot, lstart;
int a = qrtree->a;
int p = qrtree->p;
int gmt = qrtree->mt;
- ostart = start;
- opivot = pivot;
- start = hqr_getinvperm( qrtree, k, ostart);
- pivot = hqr_getinvperm( qrtree, k, opivot);
-
lpivot = pivot / p; /* Local index in the distribution over p domains */
rpivot = pivot % p; /* Staring index in this distribution */
/* Local index in the distribution over p domains */
lstart = ( start == gmt ) ? arg->llvl->ldd * a : start / p;
- perm = arg->perm + (gmt+1) * k;
-
myassert( start > pivot && pivot >= k );
- myassert( start == gmt || opivot == hqr_currpiv( qrtree, k, ostart ) );
+ myassert( start == gmt || pivot == hqr_currpiv( qrtree, k, start ) );
/* TS level common to every case */
- ls = (start < gmt) ? hqr_gettype( qrtree, k, ostart ) : -1;
- lp = hqr_gettype( qrtree, k, opivot );
+ ls = (start < gmt) ? hqr_gettype( qrtree, k, start ) : -1;
+ lp = hqr_gettype( qrtree, k, pivot );
switch( ls )
{
@@ -1406,7 +1454,7 @@ static int hqr_nextpiv(const libhqr_tree_t *qrtree, int k, int pivot, int start)
if ( ( nextp < gmt ) &&
( nextp < pivot + a*p ) &&
( (nextp/p)%a != 0 ) )
- return perm[nextp];
+ return nextp;
start = gmt;
lstart = arg->llvl->ldd * a;
@@ -1424,7 +1472,7 @@ static int hqr_nextpiv(const libhqr_tree_t *qrtree, int k, int pivot, int start)
tmp = arg->llvl->nextpiv(arg->llvl, k, pivot, tmp);
if ( tmp != arg->llvl->ldd )
- return perm[tmp * a * p + rpivot];
+ return tmp * a * p + rpivot;
next_2:
/* no next of type 1, we reset start to search the next 2 */
@@ -1442,7 +1490,7 @@ static int hqr_nextpiv(const libhqr_tree_t *qrtree, int k, int pivot, int start)
(start == gmt) &&
(lpivot < k) &&
(pivot+p < gmt) ) {
- return perm[pivot+p];
+ return pivot+p;
}
/* no next of type 2, we reset start to search the next 3 */
@@ -1458,7 +1506,7 @@ static int hqr_nextpiv(const libhqr_tree_t *qrtree, int k, int pivot, int start)
if( arg->hlvl != NULL ) {
tmp = arg->hlvl->nextpiv( arg->hlvl, k, pivot, start );
if ( tmp != gmt )
- return perm[tmp];
+ return tmp;
}
default:
@@ -1489,28 +1537,21 @@ hqr_prevpiv(const libhqr_tree_t *qrtree, int k, int pivot, int start)
{
hqr_args_t *arg = (hqr_args_t*)(qrtree->args);
int tmp, ls, lp, nextp;
- int opivot, ostart; /* original values before permutation */
- int lpivot, rpivot, lstart, *perm;
+ int lpivot, rpivot, lstart;
int a = qrtree->a;
int p = qrtree->p;
int gmt = qrtree->mt;
- ostart = start;
- opivot = pivot;
- start = hqr_getinvperm( qrtree, k, ostart );
- pivot = hqr_getinvperm( qrtree, k, opivot );
-
lpivot = pivot / p; /* Local index in the distribution over p domains */
rpivot = pivot % p; /* Staring index in this distribution */
lstart = start / p; /* Local index in the distribution over p domains */
- perm = arg->perm + (gmt+1) * k;
myassert( start >= pivot && pivot >= k && start < gmt );
- myassert( start == pivot || opivot == hqr_currpiv( qrtree, k, ostart ) );
+ myassert( start == pivot || pivot == hqr_currpiv( qrtree, k, start ) );
/* T Slevel common to every case */
- ls = hqr_gettype( qrtree, k, ostart );
- lp = hqr_gettype( qrtree, k, opivot );
+ ls = hqr_gettype( qrtree, k, start );
+ lp = hqr_gettype( qrtree, k, pivot );
if ( lp == LIBHQR_KILLED_BY_TS )
return gmt;
@@ -1522,7 +1563,7 @@ hqr_prevpiv(const libhqr_tree_t *qrtree, int k, int pivot, int start)
if( arg->hlvl != NULL ) {
tmp = arg->hlvl->prevpiv( arg->hlvl, k, pivot, start );
if ( tmp != gmt )
- return perm[tmp];
+ return tmp;
}
start = pivot;
@@ -1534,7 +1575,7 @@ hqr_prevpiv(const libhqr_tree_t *qrtree, int k, int pivot, int start)
if ( ( start == pivot ) &&
(start+p < gmt ) )
- return perm[start+p];
+ return start+p;
if ( lp > LIBHQR_KILLED_BY_LOCALTREE )
return gmt;
@@ -1558,7 +1599,7 @@ hqr_prevpiv(const libhqr_tree_t *qrtree, int k, int pivot, int start)
tmp = arg->llvl->prevpiv(arg->llvl, k, pivot, tmp);
if ( tmp != arg->llvl->ldd )
- return perm[tmp * a * p + rpivot];
+ return tmp * a * p + rpivot;
start = pivot;
@@ -1566,7 +1607,7 @@ hqr_prevpiv(const libhqr_tree_t *qrtree, int k, int pivot, int start)
/* Search for predecessor in TS tree */
/* if ( ( start+p < gmt ) && */
/* ( (((start+p) / p) % a) != 0 ) ) */
- /* return perm[start + p]; */
+ /* return start + p; */
if ( start == pivot ) {
tmp = lpivot + a - 1 - lpivot%a;
@@ -1579,7 +1620,7 @@ hqr_prevpiv(const libhqr_tree_t *qrtree, int k, int pivot, int start)
}
assert(nextp < gmt);
if ( pivot < nextp )
- return perm[nextp];
+ return nextp;
default:
return gmt;
@@ -1665,21 +1706,20 @@ hqr_genperm( libhqr_tree_t *qrtree )
static int
hqr_getinvperm( const libhqr_tree_t *qrtree, int k, int m )
{
- hqr_args_t *arg = (hqr_args_t*)(qrtree->args);
int gmt = qrtree->mt + 1;
int a = qrtree->a;
int p = qrtree->p;
int pa = p * a;
int start = m / pa * pa;
int stop = libhqr_imin( start + pa, gmt ) - start;
- int *perm = arg->perm + gmt * k + start;
int i;
if (a == 1)
return m;
for ( i=m%p; i < stop; i+=p ) {
- if( perm[i] == m )
+ //if( perm[i] == m )
+ if( i == m )
return i+start;
}
@@ -1812,6 +1852,7 @@ libhqr_hqr_init( libhqr_tree_t *qrtree,
int a, int p,
int domino, int tsrr )
{
+ libhqr_tree_t qrtree_init;
double ratio = 0.0;
int low_mt, minMN;
hqr_args_t *arg;
@@ -1822,7 +1863,7 @@ libhqr_hqr_init( libhqr_tree_t *qrtree,
}
if ((trans != LIBHQR_QR) &&
(trans != LIBHQR_LQ)) {
- fprintf(stderr, "libhqr_hqr_ini, illegal value of trans");
+ fprintf(stderr, "libhqr_hqr_init, illegal value of trans");
return -2;
}
if (A == NULL) {
@@ -1851,94 +1892,116 @@ libhqr_hqr_init( libhqr_tree_t *qrtree,
}
}
- qrtree->getnbgeqrf = hqr_getnbgeqrf;
- qrtree->getm = hqr_getm;
- qrtree->geti = hqr_geti;
- qrtree->gettype = hqr_gettype;
- qrtree->currpiv = hqr_currpiv;
- qrtree->nextpiv = hqr_nextpiv;
- qrtree->prevpiv = hqr_prevpiv;
+ minMN = libhqr_imin(A->mt, A->nt);
- qrtree->mt = (trans == LIBHQR_QR) ? A->mt : A->nt;
- qrtree->nt = (trans == LIBHQR_QR) ? A->nt : A->mt;
+ /* Create a temporary qrtree structure based on the functions */
+ {
+ qrtree_init.domino = domino;
+ qrtree_init.getnbgeqrf = hqr_getnbgeqrf;
+ qrtree_init.getm = hqr_getm;
+ qrtree_init.geti = hqr_geti;
+ qrtree_init.gettype = hqr_gettype;
+ qrtree_init.currpiv = hqr_currpiv;
+ qrtree_init.nextpiv = hqr_nextpiv;
+ qrtree_init.prevpiv = hqr_prevpiv;
- a = libhqr_imin( a, qrtree->mt );
+ qrtree_init.mt = (trans == LIBHQR_QR) ? A->mt : A->nt;
+ qrtree_init.nt = minMN;
- qrtree->a = a;
- qrtree->p = p;
- qrtree->args = NULL;
+ a = libhqr_imin( a, qrtree_init.mt );
- arg = (hqr_args_t*) malloc( sizeof(hqr_args_t) );
- arg->domino = domino;
- arg->tsrr = tsrr;
- arg->perm = NULL;
+ qrtree_init.a = a;
+ qrtree_init.p = p;
+ qrtree_init.args = NULL;
- arg->llvl = (hqr_subpiv_t*) malloc( sizeof(hqr_subpiv_t) );
- arg->hlvl = NULL;
+ arg = (hqr_args_t*) malloc( sizeof(hqr_args_t) );
+ arg->domino = domino;
+ arg->tsrr = tsrr;
+ arg->perm = NULL;
- minMN = libhqr_imin(A->mt, A->nt);
- low_mt = (qrtree->mt + p * a - 1) / ( p * a );
+ arg->llvl = (hqr_subpiv_t*) malloc( sizeof(hqr_subpiv_t) );
+ arg->hlvl = NULL;
- arg->llvl->minMN = minMN;
- arg->llvl->ldd = low_mt;
- arg->llvl->a = a;
- arg->llvl->p = p;
- arg->llvl->domino = domino;
-
- switch( type_llvl ) {
- case LIBHQR_FLAT_TREE :
- hqr_low_flat_init(arg->llvl);
- break;
- case LIBHQR_FIBONACCI_TREE :
- hqr_low_fibonacci_init(arg->llvl, minMN);
- break;
- case LIBHQR_BINARY_TREE :
- hqr_low_binary_init(arg->llvl);
- break;
- case LIBHQR_GREEDY1P_TREE :
- hqr_low_greedy1p_init(arg->llvl, minMN);
- break;
- case LIBHQR_GREEDY_TREE :
- default:
- hqr_low_greedy_init(arg->llvl, minMN);
- }
-
- if ( p > 1 ) {
- arg->hlvl = (hqr_subpiv_t*) malloc( sizeof(hqr_subpiv_t) );
+ low_mt = (qrtree_init.mt + p * a - 1) / ( p * a );
arg->llvl->minMN = minMN;
- arg->hlvl->ldd = qrtree->mt;
- arg->hlvl->a = a;
- arg->hlvl->p = p;
- arg->hlvl->domino = domino;
+ arg->llvl->ldd = low_mt;
+ arg->llvl->a = a;
+ arg->llvl->p = p;
+ arg->llvl->domino = domino;
- switch( type_hlvl ) {
+ switch( type_llvl ) {
case LIBHQR_FLAT_TREE :
- hqr_high_flat_init(arg->hlvl);
- break;
- case LIBHQR_GREEDY_TREE :
- hqr_high_greedy_init(arg->hlvl, minMN);
+ hqr_low_flat_init(arg->llvl);
break;
- case LIBHQR_GREEDY1P_TREE :
- hqr_high_greedy1p_init(arg->hlvl);
+ case LIBHQR_FIBONACCI_TREE :
+ hqr_low_fibonacci_init(arg->llvl, minMN);
break;
case LIBHQR_BINARY_TREE :
- hqr_high_binary_init(arg->hlvl);
+ hqr_low_binary_init(arg->llvl);
break;
- case LIBHQR_FIBONACCI_TREE :
- hqr_high_fibonacci_init(arg->hlvl);
+ case LIBHQR_GREEDY1P_TREE :
+ hqr_low_greedy1p_init(arg->llvl, minMN);
break;
+ case LIBHQR_GREEDY_TREE :
default:
- if ( ratio >= 0.5 ) {
+ hqr_low_greedy_init(arg->llvl, minMN);
+ }
+
+ if ( p > 1 ) {
+ arg->hlvl = (hqr_subpiv_t*) malloc( sizeof(hqr_subpiv_t) );
+
+ arg->llvl->minMN = minMN;
+ arg->hlvl->ldd = qrtree_init.mt;
+ arg->hlvl->a = a;
+ arg->hlvl->p = p;
+ arg->hlvl->domino = domino;
+
+ switch( type_hlvl ) {
+ case LIBHQR_FLAT_TREE :
hqr_high_flat_init(arg->hlvl);
- } else {
+ break;
+ case LIBHQR_GREEDY_TREE :
+ hqr_high_greedy_init(arg->hlvl, minMN);
+ break;
+ case LIBHQR_GREEDY1P_TREE :
+ hqr_high_greedy1p_init(arg->hlvl);
+ break;
+ case LIBHQR_BINARY_TREE :
+ hqr_high_binary_init(arg->hlvl);
+ break;
+ case LIBHQR_FIBONACCI_TREE :
hqr_high_fibonacci_init(arg->hlvl);
+ break;
+ default:
+ if ( ratio >= 0.5 ) {
+ hqr_high_flat_init(arg->hlvl);
+ } else {
+ hqr_high_fibonacci_init(arg->hlvl);
+ }
}
}
+
+ qrtree_init.args = (void*)arg;
+ hqr_genperm( &qrtree_init );
}
- qrtree->args = (void*)arg;
- hqr_genperm( qrtree );
+ /* Initialize the final QR tree */
+ memcpy( qrtree, &qrtree_init, sizeof(libhqr_tree_t) );
+ qrtree->getnbgeqrf = hqr_getnbgeqrf;
+ qrtree->getm = hqr_getm;
+ qrtree->geti = hqr_geti;
+ qrtree->gettype = rdmtx_gettype;
+ qrtree->currpiv = rdmtx_currpiv;
+ qrtree->nextpiv = rdmtx_nextpiv;
+ qrtree->prevpiv = rdmtx_prevpiv;
+ qrtree->args = malloc( qrtree->mt * qrtree->nt * sizeof(libhqr_tileinfo_t) );
+
+ /* Initialize the matrix */
+ libhqr_matrix_init(qrtree, &qrtree_init);
+
+ /* Free the initial qrtree */
+ libhqr_hqr_finalize( &qrtree_init );
return 0;
}
@@ -2652,6 +2715,7 @@ libhqr_svd_init( libhqr_tree_t *qrtree,
{
int low_mt, minMN, a = -1;
hqr_args_t *arg;
+ libhqr_tree_t qrtree_init;
if (qrtree == NULL) {
fprintf(stderr,"libhqr_svd_init, illegal value of qrtree");
@@ -2670,73 +2734,123 @@ libhqr_svd_init( libhqr_tree_t *qrtree,
/* Compute parameters */
p = libhqr_imax( p, 1 );
- qrtree->getnbgeqrf = svd_getnbgeqrf;
- qrtree->getm = svd_getm;
- qrtree->geti = svd_geti;
- qrtree->gettype = svd_gettype;
- qrtree->currpiv = svd_currpiv;
- qrtree->nextpiv = svd_nextpiv;
- qrtree->prevpiv = svd_prevpiv;
+ /* Create a temporary qrtree structure based on the functions */
+ {
+ qrtree_init.getnbgeqrf = svd_getnbgeqrf;
+ qrtree_init.getm = svd_getm;
+ qrtree_init.geti = svd_geti;
+ qrtree_init.gettype = svd_gettype;
+ qrtree_init.currpiv = svd_currpiv;
+ qrtree_init.nextpiv = svd_nextpiv;
+ qrtree_init.prevpiv = svd_prevpiv;
- qrtree->mt = (trans == LIBHQR_QR) ? A->mt : A->nt;
- qrtree->nt = (trans == LIBHQR_QR) ? A->nt : A->mt;
+ qrtree_init.mt = (trans == LIBHQR_QR) ? A->mt : A->nt;
+ qrtree_init.nt = (trans == LIBHQR_QR) ? A->nt : A->mt;
- qrtree->a = a;
- qrtree->p = p;
- qrtree->args = NULL;
+ qrtree_init.a = a;
+ qrtree_init.p = p;
+ qrtree_init.args = NULL;
- arg = (hqr_args_t*) malloc( sizeof(hqr_args_t) );
- arg->domino = 0;
- arg->tsrr = 0;
- arg->perm = NULL;
+ arg = (hqr_args_t*) malloc( sizeof(hqr_args_t) );
+ arg->domino = 0;
+ arg->tsrr = 0;
+ arg->perm = NULL;
- arg->llvl = (hqr_subpiv_t*) malloc( sizeof(hqr_subpiv_t) );
- arg->hlvl = NULL;
+ arg->llvl = (hqr_subpiv_t*) malloc( sizeof(hqr_subpiv_t) );
+ arg->hlvl = NULL;
- minMN = libhqr_imin(A->mt, A->nt);
- low_mt = (qrtree->mt + p - 1) / ( p );
+ minMN = libhqr_imin(A->mt, A->nt);
+ low_mt = (qrtree_init.mt + p - 1) / ( p );
- arg->llvl->minMN = minMN;
- arg->llvl->ldd = low_mt;
- arg->llvl->a = a;
- arg->llvl->p = p;
- arg->llvl->domino = 0;
+ arg->llvl->minMN = minMN;
+ arg->llvl->ldd = low_mt;
+ arg->llvl->a = a;
+ arg->llvl->p = p;
+ arg->llvl->domino = 0;
- svd_low_adaptiv_init(arg->llvl, qrtree->mt, qrtree->nt,
- nbthread_per_node * (A->nodes / p), ratio );
+ svd_low_adaptiv_init(arg->llvl, qrtree_init.mt, qrtree_init.nt,
+ nbthread_per_node * (A->nodes / p), ratio );
- if ( p > 1 ) {
- arg->hlvl = (hqr_subpiv_t*) malloc( sizeof(hqr_subpiv_t) );
+ if ( p > 1 ) {
+ arg->hlvl = (hqr_subpiv_t*) malloc( sizeof(hqr_subpiv_t) );
- arg->llvl->minMN = minMN;
- arg->hlvl->ldd = qrtree->mt;
- arg->hlvl->a = a;
- arg->hlvl->p = p;
- arg->hlvl->domino = 0;
+ arg->llvl->minMN = minMN;
+ arg->hlvl->ldd = qrtree_init.mt;
+ arg->hlvl->a = a;
+ arg->hlvl->p = p;
+ arg->hlvl->domino = 0;
- switch( type_hlvl ) {
- case LIBHQR_FLAT_TREE :
- hqr_high_flat_init(arg->hlvl);
- break;
- case LIBHQR_GREEDY_TREE :
- hqr_high_greedy_init(arg->hlvl, minMN);
- break;
- case LIBHQR_GREEDY1P_TREE :
- hqr_high_greedy1p_init(arg->hlvl);
- break;
- case LIBHQR_BINARY_TREE :
- hqr_high_binary_init(arg->hlvl);
- break;
- case LIBHQR_FIBONACCI_TREE :
- hqr_high_fibonacci_init(arg->hlvl);
- break;
- default:
- hqr_high_fibonacci_init(arg->hlvl);
+ switch( type_hlvl ) {
+ case LIBHQR_FLAT_TREE :
+ hqr_high_flat_init(arg->hlvl);
+ break;
+ case LIBHQR_GREEDY_TREE :
+ hqr_high_greedy_init(arg->hlvl, minMN);
+ break;
+ case LIBHQR_GREEDY1P_TREE :
+ hqr_high_greedy1p_init(arg->hlvl);
+ break;
+ case LIBHQR_BINARY_TREE :
+ hqr_high_binary_init(arg->hlvl);
+ break;
+ case LIBHQR_FIBONACCI_TREE :
+ hqr_high_fibonacci_init(arg->hlvl);
+ break;
+ default:
+ hqr_high_fibonacci_init(arg->hlvl);
+ }
}
+
+ qrtree_init.args = (void*)arg;
+ hqr_genperm( &qrtree_init );
}
- qrtree->args = (void*)arg;
+ /* Initialize the final QR tree */
+ memcpy( qrtree, &qrtree_init, sizeof(libhqr_tree_t) );
+ qrtree->getnbgeqrf = svd_getnbgeqrf;
+ qrtree->getm = svd_getm;
+ qrtree->geti = svd_geti;
+ qrtree->gettype = rdmtx_gettype;
+ qrtree->currpiv = rdmtx_currpiv;
+ qrtree->nextpiv = rdmtx_nextpiv;
+ qrtree->prevpiv = rdmtx_prevpiv;
+ qrtree->args = malloc( qrtree->mt * qrtree->nt * sizeof(libhqr_tileinfo_t) );
+
+ /* Initialize the matrix */
+ libhqr_matrix_init(qrtree, &qrtree_init);
+
+ /* Free the initial qrtree */
+ libhqr_hqr_finalize( &qrtree_init );
return 0;
}
+void libhqr_matrix_init(libhqr_tree_t *qrtree, const libhqr_tree_t *qrtree_init){
+ int i, minMN, tile_id, p, k;
+ libhqr_tileinfo_t *arg = (libhqr_tileinfo_t*)(qrtree->args);
+ minMN = libhqr_imin(qrtree->nt, qrtree->mt);
+ for (k = 0; k < minMN; k++){
+ for (i = 0; i < qrtree->mt; i++){
+ tile_id = (k * qrtree->mt) + i;
+
+ arg[tile_id].type = qrtree_init->gettype(qrtree_init, k, i);
+ assert( (i < k) || (arg[tile_id].type >= 0) );
+ arg[tile_id].currpiv = qrtree_init->currpiv(qrtree_init, k, i);
+ p = arg[tile_id].currpiv;
+ arg[tile_id].first_nextpiv = qrtree_init->nextpiv(qrtree_init, k, i, qrtree->mt);
+ arg[tile_id].first_prevpiv = qrtree_init->prevpiv(qrtree_init, k, i, i);
+ arg[tile_id].nextpiv = qrtree_init->nextpiv(qrtree_init, k, p, i);
+ arg[tile_id].prevpiv = qrtree_init->prevpiv(qrtree_init, k, p, i);
+ assert( (arg[tile_id].nextpiv <= qrtree->mt && arg[tile_id].nextpiv >= k) || arg[tile_id].nextpiv==-1 );
+ assert( (arg[tile_id].prevpiv <= qrtree->mt && arg[tile_id].prevpiv >= k) || arg[tile_id].prevpiv==-1 );
+ assert( (arg[tile_id].first_nextpiv <= qrtree->mt && arg[tile_id].first_nextpiv >= k) || arg[tile_id].first_nextpiv==-1 );
+ assert( (arg[tile_id].first_prevpiv <= qrtree->mt && arg[tile_id].first_prevpiv >= k) || arg[tile_id].first_prevpiv==-1 );
+ }
+ }
+ qrtree->args = (void*)arg;
+}
+
+void libhqr_matrix_finalize(libhqr_tree_t *qrtree){
+ libhqr_tileinfo_t *arg = (libhqr_tileinfo_t*)(qrtree->args);
+ free(arg);
+}
diff --git a/src/libhqr_dbg.c b/src/libhqr_dbg.c
index 46d82ef0a7e92db340cc4ab08271813c3cbd97cd..cc333f41fbafe64b3061bd13356f7dfc53c8f959 100644
--- a/src/libhqr_dbg.c
+++ b/src/libhqr_dbg.c
@@ -19,13 +19,16 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
+#include <assert.h>
/* static int libhqr_qrtree_getinon0( const qr_piv_t *arg, */
/* const int k, int i, int mt ); */
#define ENDCHECK( test, ret ) \
- if ( !test ) \
- return ret;
+ if ( !test ) { \
+ assert( ret == 0 ); \
+ return ret; \
+ }
int libhqr_tree_check( libhqr_tiledesc_t *A, libhqr_tree_t *qrtree)
{
@@ -152,7 +155,7 @@ int libhqr_tree_check( libhqr_tiledesc_t *A, libhqr_tree_t *qrtree)
a, p, A->mt, A->nt,
m, k);
check = 0;
- return 3;
+ ENDCHECK( check, 3 );
}
}
}
@@ -195,7 +198,7 @@ int libhqr_tree_check( libhqr_tiledesc_t *A, libhqr_tree_t *qrtree)
a, p, A->mt, A->nt,
m, k);
check = 0;
- return 3;
+ ENDCHECK( check, 3 );
}
}
}
@@ -227,11 +230,11 @@ int libhqr_tree_check( libhqr_tiledesc_t *A, libhqr_tree_t *qrtree)
printf(" ----------------------------------------------------\n"
" - a = %d, p = %d, M = %d, N = %d\n"
" Check next/prev:\n"
- " next( m=%d, k=%d, start=%d ) => %d && prev( m=%d, k=%d, start=%d ) => %d\n ( %d != %d )",
+ " next( k=%d, m=%d, start=%d ) => %d && prev( k=%d, m=%d, start=%d ) => %d (instead of %d)\n",
a, p, A->mt, A->nt,
- m, k, start, next, m, k, next, prev, start, prev);
- check = 0;
- return 3;
+ k, m, start, next, k, m, next, prev, start);
+ check = 0;
+ ENDCHECK( check, 3 );
}
start = next;
} while ( start != A->mt );
diff --git a/src/libhqr_systolic.c b/src/libhqr_systolic.c
index 66153a3b0e66ccefdbc371c6aef8167cb8e87da9..d66d3bb9666b06bd142ad055f2e7f84d1b665248 100644
--- a/src/libhqr_systolic.c
+++ b/src/libhqr_systolic.c
@@ -19,6 +19,7 @@
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
+#include <string.h>
#define PRINT_PIVGEN 0
#ifdef PRINT_PIVGEN
@@ -352,6 +353,8 @@ libhqr_systolic_init( libhqr_tree_t *qrtree,
libhqr_typefacto_e trans, libhqr_tiledesc_t *A,
int p, int q )
{
+ libhqr_tree_t qrtree_init;
+
if (qrtree == NULL) {
fprintf(stderr, "libhqr_systolic_init, illegal value of qrtree");
return -1;
@@ -374,20 +377,39 @@ libhqr_systolic_init( libhqr_tree_t *qrtree,
return -5;
}
+ /* Create a temporary qrtree structure based on the functions */
+ {
+ qrtree_init.getnbgeqrf = systolic_getnbgeqrf;
+ qrtree_init.getm = systolic_getm;
+ qrtree_init.geti = systolic_geti;
+ qrtree_init.gettype = systolic_gettype;
+ qrtree_init.currpiv = systolic_currpiv;
+ qrtree_init.nextpiv = systolic_nextpiv;
+ qrtree_init.prevpiv = systolic_prevpiv;
+
+ qrtree_init.mt = (trans == LIBHQR_QR) ? A->mt : A->nt;
+ qrtree_init.nt = (trans == LIBHQR_QR) ? A->nt : A->mt;
+
+ qrtree_init.a = libhqr_imax( q, 1 );
+ qrtree_init.p = libhqr_imax( p, 1 );
+ qrtree_init.args = NULL;
+ }
+
+ memcpy( qrtree, &qrtree_init, sizeof(libhqr_tree_t) );
qrtree->getnbgeqrf = systolic_getnbgeqrf;
qrtree->getm = systolic_getm;
qrtree->geti = systolic_geti;
- qrtree->gettype = systolic_gettype;
- qrtree->currpiv = systolic_currpiv;
- qrtree->nextpiv = systolic_nextpiv;
- qrtree->prevpiv = systolic_prevpiv;
+ qrtree->gettype = rdmtx_gettype;
+ qrtree->currpiv = rdmtx_currpiv;
+ qrtree->nextpiv = rdmtx_nextpiv;
+ qrtree->prevpiv = rdmtx_prevpiv;
+ qrtree->args = malloc( qrtree->mt * qrtree->nt * sizeof(libhqr_tileinfo_t) );
- qrtree->mt = (trans == LIBHQR_QR) ? A->mt : A->nt;
- qrtree->nt = (trans == LIBHQR_QR) ? A->nt : A->mt;
+ /*Initialize the matrix */
+ libhqr_matrix_init(qrtree, &qrtree_init);
- qrtree->a = libhqr_imax( q, 1 );
- qrtree->p = libhqr_imax( p, 1 );
- qrtree->args = NULL;
+ /* Free the initial qrtree */
+ libhqr_systolic_finalize( &qrtree_init );
return 0;
}
diff --git a/testings/testing_pivgen.c b/testings/testing_pivgen.c
index 3fa5adde0a311b8dbcc183520723ca87cd2b784e..3b805735395defff14c665540cc4f458a6c0f216 100644
--- a/testings/testing_pivgen.c
+++ b/testings/testing_pivgen.c
@@ -76,8 +76,10 @@ int main(int argc, char ** argv)
fprintf(stderr, "-M %d -N %d --treel=%d --qr_a=%d --tsrr=%d FAILED(%d)\n",
allM[m], allN[n], alltreel[l], allA[a], r, rc);
ret |= 1;
+ libhqr_matrix_finalize( &qrtree );
+ return 0;
}
- libhqr_hqr_finalize( &qrtree );
+ libhqr_matrix_finalize( &qrtree );
done++;
printf("\r%6d / %6d", done, todo);
@@ -110,7 +112,7 @@ int main(int argc, char ** argv)
ret |= 1;
}
- libhqr_hqr_finalize( &qrtree );
+ libhqr_matrix_finalize( &qrtree );
done++;
printf("\r%6d / %6d", done, todo);
@@ -145,7 +147,7 @@ int main(int argc, char ** argv)
ret |= 1;
}
- libhqr_systolic_finalize( &qrtree );
+ libhqr_matrix_finalize( &qrtree );
done++;
printf("\r%6d / %6d", done, todo);
diff --git a/testings/testing_tileinit.c b/testings/testing_tileinit.c
new file mode 100644
index 0000000000000000000000000000000000000000..895aa63781f4c5a173e1c8b0f688f60a2b981bc5
--- /dev/null
+++ b/testings/testing_tileinit.c
@@ -0,0 +1,40 @@
+/**
+ *
+ * @file testting_tileinit.c
+ *
+ * Testing file for the functions stocking the informations of the tiles
+ *
+ * @copyright 2017 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
+ * Univ. Bordeaux. All rights reserved.
+ *
+ * @version 1.0.0
+ * @author Raphael Boucherie
+ * @author Mathieu Faverge
+ * @date 2017-04-18
+ *
+ */
+
+
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+#include <assert.h>
+#include "libhqr.h"
+
+int
+main(int argc, char ** argv)
+{
+ int rc;
+ libhqr_tree_t qrtree;
+ libhqr_tiledesc_t matrix;
+
+ memset( &qrtree, 0xdeadbeef, sizeof(libhqr_tree_t) );
+ matrix.nodes = 1;
+ matrix.p = 3;
+ matrix.mt = 1;
+ matrix.nt = 1;
+ libhqr_hqr_init( &qrtree, LIBHQR_QR, &matrix, 0, 0, 1, 3, 0, 0);
+ rc = libhqr_tree_check( &matrix, &qrtree );
+ libhqr_hqr_finalize( &qrtree );
+ return 1;
+}