Add low level trees functions inseparated files

src/low_binary.c

+/**
+ *
+ * @file low_binary.c
+ *
+ * @copyright 2010-2017 The University of Tennessee and The University
+ *                      of Tennessee Research Foundation.  All rights
+ *                      reserved.
+ * @copyright 2012-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-03-21
+ *
+ * Functions for low level binary tree
+ *
+ */
+#include "libhqr_internal.h"
+#include <math.h>
+
+static inline int
+hqr_low_binary_currpiv(const hqr_subpiv_t *arg, int k, int m)
+{
+    int k_a = arg->domino ? k / arg->a :  (k + arg->p - 1 - m%(arg->p)) / arg->p / arg->a;
+    int m_pa = (m / arg->p ) / arg->a;
+
+    int tmp1 = m_pa - k_a;
+    int tmp2 = 1;
+    (void)arg;
+
+    if ( tmp1 == 0)
+        return 0;
+    while( (tmp1 != 0 ) && (tmp1 % 2 == 0) ) {
+        tmp1 = tmp1 >> 1;
+        tmp2 = tmp2 << 1;
+    }
+    assert( m_pa - tmp2 >= k_a );
+    return m_pa - tmp2;
+}
+
+static inline int
+hqr_low_binary_nextpiv(const hqr_subpiv_t *arg, int k, int p, int start_pa)
+{
+    int k_a = arg->domino ? k / arg->a :  (k + arg->p - 1 - p%(arg->p)) / arg->p / arg->a;
+    int p_pa = (p / arg->p ) / arg->a;
+
+    int tmpp, bit;
+    myassert( (start_pa == arg->ldd) || (hqr_low_binary_currpiv( arg, k, start_pa*arg->a*arg->p ) == p_pa || !arg->domino) );
+
+    if ( start_pa <= p_pa )
+        return arg->ldd;
+
+    int offset = p_pa-k_a;
+    bit = 0;
+    if (start_pa != arg->ldd) {
+        while( ( (start_pa-k_a) & (1 << bit ) ) == 0 )
+            bit++;
+        bit++;
+    }
+
+    tmpp = offset | (1 << bit);
+    if ( ( tmpp != offset ) && ( tmpp+k_a < arg->ldd ) )
+        return tmpp + k_a;
+    else
+        return arg->ldd;
+}
+
+static inline int
+hqr_low_binary_prevpiv(const hqr_subpiv_t *arg, int k, int p, int start_pa)
+{
+    int k_a = arg->domino ? k / arg->a :  (k + arg->p - 1 - p%(arg->p)) / arg->p / arg->a;
+    int p_pa = (p / arg->p ) / arg->a;
+    int offset = p_pa - k_a;
+
+    myassert( start_pa >= p_pa && ( start_pa == p_pa || !arg->domino ||
+                                    hqr_low_binary_currpiv( arg, k, start_pa*arg->a*arg->p ) == p_pa ) );
+
+    if ( (start_pa == p_pa) && ( offset%2 == 0 ) ) {
+        int i, bit, tmp;
+        if ((p_pa - k_a) == 0)
+            bit = (int)( log( (double)(arg->ldd - k_a) ) / log( 2. ) );
+        else {
+            bit = 0;
+            while( (offset & (1 << bit )) == 0 )
+                bit++;
+        }
+        for( i=bit; i>-1; i--){
+            tmp = offset | (1 << i);
+            if ( ( offset != tmp ) && ( tmp+k_a < arg->ldd ) )
+                return tmp+k_a;
+        }
+        return arg->ldd;
+    }
+
+    if ( (start_pa - p_pa) > 1 )
+        return p_pa + ( (start_pa - p_pa) >> 1 );
+    else {
+        return arg->ldd;
+    }
+}
+
+void
+hqr_low_binary_init(hqr_subpiv_t *arg) {
+    arg->currpiv = hqr_low_binary_currpiv;
+    arg->nextpiv = hqr_low_binary_nextpiv;
+    arg->prevpiv = hqr_low_binary_prevpiv;
+    arg->ipiv = NULL;
+}

src/low_fibonacci.c

+/**
+ *
+ * @file low_fibonacci.c
+ *
+ * @copyright 2010-2017 The University of Tennessee and The University
+ *                      of Tennessee Research Foundation.  All rights
+ *                      reserved.
+ * @copyright 2012-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-03-21
+ *
+ * Functions for low level fibonacci tree
+ *
+ */
+#include "libhqr_internal.h"
+#include <stdlib.h>
+
+/* Return the pivot to use for the row m at step k */
+static inline int
+hqr_low_fibonacci_currpiv( const hqr_subpiv_t *qrpiv, int k, int m ) {
+    int k_a = qrpiv->domino ? k / qrpiv->a : (k + qrpiv->p - 1 - m%(qrpiv->p)) / qrpiv->p / qrpiv->a;
+    return (qrpiv->ipiv)[ k_a * (qrpiv->ldd) + ( (m / qrpiv->p) / qrpiv->a ) ];
+}
+
+/* Return the last row which has used the row m as a pivot in step k before the row start */
+static inline int
+hqr_low_fibonacci_prevpiv( const hqr_subpiv_t *qrpiv, int k, int p, int start_pa ) {
+    int i;
+    int k_a = qrpiv->domino ? k / qrpiv->a :  (k + qrpiv->p - 1 - p%(qrpiv->p)) / qrpiv->p / qrpiv->a;
+    int p_pa = (p / qrpiv->p ) / qrpiv->a;
+
+    for( i=start_pa+1; i<(qrpiv->ldd); i++ )
+        if ( (qrpiv->ipiv)[i +  k_a * (qrpiv->ldd)] == p_pa )
+            return i;
+    return i;
+}
+
+/* Return the next row which will use the row m as a pivot in step k after it has been used by row start */
+static inline int
+hqr_low_fibonacci_nextpiv( const hqr_subpiv_t *qrpiv, int k, int p, int start_pa ) {
+    int i;
+    int k_a = qrpiv->domino ? k / qrpiv->a :  (k + qrpiv->p - 1 - p%(qrpiv->p)) / qrpiv->p / qrpiv->a;
+    int p_pa = (p / qrpiv->p ) / qrpiv->a;
+
+    for( i=start_pa-1; i>k_a; i-- )
+        if ( (qrpiv->ipiv)[i + k_a * (qrpiv->ldd)] == p_pa )
+            return i;
+    return (qrpiv->ldd);
+}
+
+void
+hqr_low_fibonacci_init(hqr_subpiv_t *arg, int minMN) {
+    int *ipiv;
+    int mt;
+
+    arg->currpiv = hqr_low_fibonacci_currpiv;
+    arg->nextpiv = hqr_low_fibonacci_nextpiv;
+    arg->prevpiv = hqr_low_fibonacci_prevpiv;
+
+    mt = arg->ldd;
+
+    arg->ipiv = (int*)calloc( mt * minMN, sizeof(int) );
+    ipiv = arg->ipiv;
+
+    /*
+     * Fibonacci of order 1:
+     *    u_(n+1) = u_(n) + 1
+     */
+    {
+        int f1, k, m;
+
+        /* Fill in the first column */
+        f1 = 1;
+        for (m=1; m < mt; ) {
+            for (k=0; (k < f1) && (m < mt); k++, m++) {
+                ipiv[m] = m - f1;
+            }
+            f1++;
+        }
+
+        for( k=1; k<minMN; k++) {
+            for(m=k+1; m < mt; m++) {
+                ipiv[ k * mt + m ] = ipiv[ (k-1) * mt + m - 1 ] + 1;
+            }
+        }
+    }
+}

src/low_flat.c

+/**
+ *
+ * @file low_flat.c
+ *
+ * @copyright 2010-2017 The University of Tennessee and The University
+ *                      of Tennessee Research Foundation.  All rights
+ *                      reserved.
+ * @copyright 2012-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-03-21
+ *
+ * Functions for low level flat tree
+ *
+ */
+#include "libhqr_internal.h"
+
+static inline int
+hqr_low_flat_currpiv(const hqr_subpiv_t *arg, int k, int m)
+{
+    (void)m;
+    if ( arg->domino )
+        return k / arg->a;
+    else
+        return (k + arg->p - 1 - m%(arg->p)) / arg->p / arg->a ;
+}
+
+static inline int
+hqr_low_flat_nextpiv(const hqr_subpiv_t *arg, int k, int p, int start_pa)
+{
+    int k_a = arg->domino ? k / arg->a :  (k + arg->p - 1 - p%(arg->p)) / arg->p / arg->a;
+    int p_pa = (p / arg->p ) / arg->a;
+
+#ifdef FLAT_UP
+    if ( ( p_pa == k_a ) && (start_pa > k_a+1 ) )
+        return start_pa-1;
+    else
+#else /* FLAT_DOWN */
+        if ( start_pa <= p_pa )
+            return arg->ldd;
+
+    if ( p_pa == k_a && ( arg->ldd - k_a ) > 1 ) {
+        if ( start_pa == arg->ldd )
+            return p_pa+1;
+        else if ( start_pa < arg->ldd )
+            return start_pa+1;
+    }
+#endif
+    return arg->ldd;
+}
+
+static inline int
+hqr_low_flat_prevpiv(const hqr_subpiv_t *arg, int k, int p, int start_pa)
+{
+    int k_a = arg->domino ? k / arg->a :  (k + arg->p - 1 - p%(arg->p)) / arg->p / arg->a;
+    int p_pa = (p / arg->p ) / arg->a;
+
+#ifdef FLAT_UP
+    if ( p_pa == k_a && (start_pa+1 < arg->ldd) )
+        return start_pa+1;
+    else
+#else
+        if ( p_pa == k_a && ( arg->ldd - k_a ) > 1  ) {
+            if ( start_pa == p_pa )
+                return arg->ldd - 1;
+            else if ( start_pa > p_pa + 1 )
+                return start_pa-1;
+        }
+#endif
+    return arg->ldd;
+}
+
+void hqr_low_flat_init(hqr_subpiv_t *arg){
+    arg->currpiv = hqr_low_flat_currpiv;
+    arg->nextpiv = hqr_low_flat_nextpiv;
+    arg->prevpiv = hqr_low_flat_prevpiv;
+    arg->ipiv = NULL;
+}

src/low_greedy.c

+/**
+ *
+ * @file low_greedy.c
+ *
+ * @copyright 2010-2017 The University of Tennessee and The University
+ *                      of Tennessee Research Foundation.  All rights
+ *                      reserved.
+ * @copyright 2012-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-03-21
+ *
+ * Functions for low level greedy tree
+ *
+ */
+#include "libhqr_internal.h"
+#include <stdlib.h>
+#include <string.h>
+
+/* Return the pivot to use for the row m at step k */
+static inline int
+hqr_low_greedy_currpiv( const hqr_subpiv_t *qrpiv, int k, int m ) {
+    if (qrpiv->domino)
+        return (qrpiv->ipiv)[ k * (qrpiv->ldd) + ( (m / qrpiv->p) / qrpiv->a ) ];
+    else
+        return (qrpiv->ipiv)[ ( (m%qrpiv->p) * qrpiv->minMN + k ) * (qrpiv->ldd)
+                              + ( ( m  / qrpiv->p ) / qrpiv->a ) ];
+}
+
+/* Return the last row which has used the row m as a pivot in step k before the row start */
+static inline int
+hqr_low_greedy_prevpiv( const hqr_subpiv_t *qrpiv, int k, int p, int start_pa ) {
+    int i;
+    int p_pa = p / qrpiv->p / qrpiv->a;
+    int *ipiv = qrpiv->domino ? qrpiv->ipiv : qrpiv->ipiv + p%qrpiv->p * qrpiv->minMN *qrpiv->ldd;
+
+    for( i=start_pa+1; i<(qrpiv->ldd); i++ )
+        if ( ipiv[i +  k * (qrpiv->ldd)] == p_pa )
+            return i;
+    return i;
+ }
+
+/* Return the next row which will use the row m as a pivot in step k after it has been used by row start */
+static inline int
+hqr_low_greedy_nextpiv( const hqr_subpiv_t *qrpiv, int k, int p, int start_pa ) {
+    int i;
+    int pa = qrpiv->p * qrpiv->a;
+    int k_a = qrpiv->domino ? k / qrpiv->a :  (k + qrpiv->p - 1 - p%(qrpiv->p)) / qrpiv->p / qrpiv->a;
+    int p_pa = p / pa;
+    int *ipiv = qrpiv->domino ? qrpiv->ipiv : qrpiv->ipiv + p%qrpiv->p * qrpiv->minMN *qrpiv->ldd;
+
+    for( i=start_pa-1; i> k_a; i-- )
+        if ( ipiv[i + k * (qrpiv->ldd)] == p_pa )
+            return i;
+
+    return (qrpiv->ldd);
+}
+
+void
+hqr_low_greedy_init(hqr_subpiv_t *arg, int minMN) {
+    int *ipiv;
+    int mt, a, p, pa, domino;
+
+    arg->currpiv = hqr_low_greedy_currpiv;
+    arg->nextpiv = hqr_low_greedy_nextpiv;
+    arg->prevpiv = hqr_low_greedy_prevpiv;
+
+    mt = arg->ldd;
+    a = arg->a;
+    p = arg->p;
+    pa = p * a;
+    domino = arg->domino;
+
+    if ( domino )
+    {
+        int j, k, height, start, end, firstk = 0;
+        int *nT, *nZ;
+
+        arg->minMN =  libhqr_imin( minMN, mt*a );
+        minMN = arg->minMN;
+
+        arg->ipiv = (int*)calloc( mt * minMN, sizeof(int) );
+        ipiv = arg->ipiv;
+
+        nT = (int*)calloc(minMN, sizeof(int));
+        nZ = (int*)calloc(minMN, sizeof(int));
+
+        nT[0] = mt;
+
+        k = 0;
+        while ( (!( ( nT[minMN-1] == mt - ( (minMN - 1) / a ) ) &&
+                     ( nZ[minMN-1]+1 == nT[minMN-1] ) ) )
+                && ( firstk < minMN ) ) {
+            height = (nT[k] - nZ[k]) / 2;
+            if ( height == 0 ) {
+                while ( (firstk < minMN) &&
+                        ( nT[firstk] == mt - ( firstk / a ) ) &&
+                        ( nZ[firstk]+1 == nT[firstk] ) ) {
+                    if (  (( firstk % a) != a-1 )
+                          && ( firstk < minMN-1 ) )
+                        nT[firstk+1]++;
+                    firstk++;
+                }
+                k = firstk;
+                continue;
+            }
+
+            if (k < minMN-1) nT[k+1] += height;
+            start = mt - nZ[k] - 1;
+            end = start - height;
+            nZ[k] += height;
+
+            for( j=start; j > end; j-- ) {
+                ipiv[ k*mt + j ] = (j - height);
+            }
+
+            k++;
+            if (k > minMN-1) k = firstk;
+        }
+
+        free(nT);
+        free(nZ);
+    }
+    else
+    {
+        int j, k, myrank, height, start, end, firstk = 0;
+        int *nT2DO = (int*)malloc(minMN*sizeof(int));
+        int *nT = (int*)malloc(minMN*sizeof(int));
+        int *nZ = (int*)malloc(minMN*sizeof(int));
+
+        arg->ipiv = (int*)calloc( mt * minMN * p, sizeof(int) );
+        ipiv = arg->ipiv;
+
+        for ( myrank=0; myrank<p; myrank++ ) {
+            int lminMN = minMN;
+
+            memset( nT2DO, 0, minMN*sizeof(int));
+            memset( nT,    0, minMN*sizeof(int));
+            memset( nZ,    0, minMN*sizeof(int));
+
+            nT[0] = mt;
+
+            for(k=0; k<lminMN; k++) {
+                nT2DO[k] = libhqr_imax( mt - ((k + p - 1 - myrank) / pa), 0 );
+                if ( nT2DO[k] == 0 ) {
+                    lminMN = k;
+                    break;
+                }
+            }
+
+            k = 0; firstk = 0;
+            while ( (!( ( nT[lminMN-1] == nT2DO[lminMN-1] ) &&
+                        ( nZ[lminMN-1]+1 == nT[lminMN-1] ) ) )
+                    && ( firstk < lminMN ) ) {
+                height = (nT[k] - nZ[k]) / 2;
+                if ( height == 0 ) {
+                    while ( (firstk < lminMN) &&
+                            ( nT[firstk] == nT2DO[firstk] ) &&
+                            ( nZ[firstk]+1 == nT[firstk] ) ) {
+                        if (  ( firstk < lminMN-1 )  &&
+                              (( (firstk) % pa) != ((a-1)*p+myrank) ) )
+                            nT[firstk+1]++;
+                        firstk++;
+                    }
+                    k = firstk;
+                    continue;
+                }
+
+                if (k < lminMN-1) nT[k+1] += height;
+                start = mt - nZ[k] - 1;
+                end = start - height;
+                nZ[k] += height;
+
+                for( j=start; j > end; j-- ) {
+                    ipiv[ myrank*mt*minMN + k*mt + j ] = (j - height);
+                }
+
+                k++;
+                if (k > lminMN-1) k = firstk;
+            }
+        }
+
+        free(nT2DO);
+        free(nT);
+        free(nZ);
+    }
+
+#if 0
+    {
+        int m, k;
+        for(m=0; m<mt; m++) {
+            printf("%3d | ", m);
+            for (k=0; k<minMN; k++) {
+                printf( "%3d ", ipiv[k*mt + m] );
+            }
+            printf("\n");
+        }
+    }
+    if (!arg->domino) {
+        int m, k, myrank;
+        for ( myrank=1; myrank<p; myrank++ ) {
+            ipiv += mt*minMN;
+            printf("-------- rank %d ---------\n", myrank );
+            for(m=0; m<mt; m++) {
+                printf("%3d | ", m);
+                for (k=0; k<minMN; k++) {
+                  int k_a = (k + p - 1 - myrank) / p / a;
+                  if ( m >= k_a )
+                    printf( "%3d ", ipiv[k*mt + m] );
+                  else
+                    printf( "--- " );
+                }
+                printf("\n");
+            }
+        }
+    }
+#endif
+}

src/low_greedy1p.c

+/**
+ *
+ * @file low_greedy1p.c
+ *
+ * @copyright 2010-2017 The University of Tennessee and The University
+ *                      of Tennessee Research Foundation.  All rights
+ *                      reserved.
+ * @copyright 2012-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-03-21
+ *
+ * Functions for low level duplicated greedy tree
+ *
+ */
+#include "libhqr_internal.h"
+
+/* Return the pivot to use for the row m at step k */
+static inline int
+hqr_low_greedy1p_currpiv( const hqr_subpiv_t *qrpiv, int k, int m ) {
+    if (qrpiv->domino)
+        return (qrpiv->ipiv)[ k * (qrpiv->ldd) + ( (m / qrpiv->p) / qrpiv->a ) ];
+    else
+        return (qrpiv->ipiv)[ ( (m%qrpiv->p) * qrpiv->minMN + k ) * (qrpiv->ldd)
+                              + ( ( m  / qrpiv->p ) / qrpiv->a ) ];
+}
+
+/* Return the last row which has used the row m as a pivot in step k before the row start */
+static inline int
+hqr_low_greedy1p_prevpiv( const hqr_subpiv_t *qrpiv, int k, int p, int start_pa ) {
+    int i;
+    int p_pa = p / qrpiv->p / qrpiv->a;
+    int *ipiv = qrpiv->domino ? qrpiv->ipiv : qrpiv->ipiv + p%qrpiv->p * qrpiv->minMN *qrpiv->ldd;
+
+    for( i=start_pa+1; i<(qrpiv->ldd); i++ )
+        if ( ipiv[i +  k * (qrpiv->ldd)] == p_pa )
+            return i;
+    return i;
+}
+
+/* Return the next row which will use the row m as a pivot in step k after it has been used by row start */
+static inline int
+hqr_low_greedy1p_nextpiv( const hqr_subpiv_t *qrpiv, int k, int p, int start_pa ) {
+    int i;
+    int pa = qrpiv->p * qrpiv->a;
+    int k_a = qrpiv->domino ? k / qrpiv->a :  (k + qrpiv->p - 1 - p%(qrpiv->p)) / qrpiv->p / qrpiv->a;
+    int p_pa = p / pa;
+    int *ipiv = qrpiv->domino ? qrpiv->ipiv : qrpiv->ipiv + p%qrpiv->p * qrpiv->minMN *qrpiv->ldd;
+
+    for( i=start_pa-1; i> k_a; i-- )
+        if ( ipiv[i + k * (qrpiv->ldd)] == p_pa )
+            return i;
+
+    return (qrpiv->ldd);
+}
+
+void
+hqr_low_greedy1p_init(hqr_subpiv_t *arg, int minMN) {
+    int *ipiv;
+    int mt, a, p, pa, domino;
+    int j, k, height, start, end, nT, nZ;
+
+    arg->currpiv = hqr_low_greedy1p_currpiv;
+    arg->nextpiv = hqr_low_greedy1p_nextpiv;
+    arg->prevpiv = hqr_low_greedy1p_prevpiv;
+
+    mt = arg->ldd;
+    a = arg->a;
+    p = arg->p;
+    pa = p * a;
+    domino = arg->domino;
+
+    /* This section has not been coded yet, and will perform a classic greedy */
+    if ( domino )
+    {
+        arg->minMN =  libhqr_imin( minMN, mt*a );
+        minMN = arg->minMN;
+
+        arg->ipiv = (int*)malloc( mt * minMN * sizeof(int) );
+        ipiv = arg->ipiv;
+        memset(ipiv, 0, mt*minMN*sizeof(int));
+
+        /**
+         * Compute the local greedy tree of each column, on each node
+         */
+        for(k=0; k<minMN; k++) {
+            /* Number of tiles to factorized in this column on this rank */
+            nT = libhqr_imax( mt - (k / a), 0 );
+            /* Number of tiles already killed */
+            nZ = 0;
+
+            while( nZ < (nT-1) ) {
+                height = (nT - nZ) / 2;
+                start = mt - nZ - 1;
+                end = start - height;
+                nZ += height;
+
+                for( j=start; j > end; j-- ) {
+                    ipiv[ k*mt + j ] = (j - height);
+                }
+            }
+            assert( nZ+1 == nT );
+        }
+    }
+    else
+    {
+        int myrank;
+        end = 0;
+
+        arg->ipiv = (int*)malloc( mt * minMN * sizeof(int) * p );
+        ipiv = arg->ipiv;
+
+        memset( ipiv,  0, minMN*sizeof(int)*mt*p);
+
+        for ( myrank=0; myrank<p; myrank++ ) {
+
+            /**
+             * Compute the local greedy tree of each column, on each node
+             */
+            for(k=0; k<minMN; k++) {
+                /* Number of tiles to factorized in this column on this rank */
+                nT = libhqr_imax( mt - ((k + p - 1 - myrank) / pa), 0 );
+                /* Number of tiles already killed */
+                nZ = 0;
+
+                /* No more computations on this node */
+                if ( nT == 0 ) {
+                    break;
+                }
+
+                while( nZ < (nT-1) ) {
+                    height = (nT - nZ) / 2;
+                    start = mt - nZ - 1;
+                    end = start - height;
+                    nZ += height;
+
+                    for( j=start; j > end; j-- ) {
+                        ipiv[ myrank*mt*minMN + k*mt + j ] = (j - height);
+                    }
+                }
+                assert( nZ+1 == nT );
+            }
+        }
+    }
+
+#if 0
+    {
+        int m, k;
+        for(m=0; m<mt; m++) {
+            printf("%3d | ", m);
+            for (k=0; k<minMN; k++) {
+                printf( "%3d ", ipiv[k*mt + m] );
+            }
+            printf("\n");
+        }
+    }
+    if (!arg->domino) {
+        int m, k, myrank;
+        for ( myrank=1; myrank<p; myrank++ ) {
+            ipiv += mt*minMN;
+            printf("-------- rank %d ---------\n", myrank );
+            for(m=0; m<mt; m++) {
+                printf("%3d | ", m);
+                for (k=0; k<minMN; k++) {
+                    int k_a = (k + p - 1 - myrank) / p / a;
+                    if ( m >= k_a )
+                        printf( "%3d ", ipiv[k*mt + m] );
+                    else
+                        printf( "--- " );
+                }
+                printf("\n");
+            }
+        }
+    }
+#endif
+}