coreblas: add ipiv_to_perm kernel

coreblas/compute/CMakeLists.txt

@@ -17,14 +17,14 @@
 #     Univ. of California Berkeley,
 #     Univ. of Colorado Denver.
 #
-# @version 1.2.0
+# @version 1.3.0
 #  @author Cedric Castagnede
 #  @author Emmanuel Agullo
 #  @author Mathieu Faverge
 #  @author Florent Pruvost
 #  @author Guillaume Sylvand
 #  @author Matthieu Kuhn
-#  @date 2022-02-22
+#  @date 2023-08-31
 #
 ###
 
@@ -132,9 +132,10 @@ precisions_rules_py(COREBLAS_SRCS_GENERATED "${ZSRC}"
                     PRECISIONS "${CHAMELEON_PRECISION}")
 
 set(COREBLAS_SRCS
-    global.c
-    ${COREBLAS_SRCS_GENERATED}
-    )
+  global.c
+  core_ipiv_to_perm.c
+  ${COREBLAS_SRCS_GENERATED}
+)
 
 # Force generation of sources
 # ---------------------------

coreblas/compute/core_ipiv_to_perm.c

+/**
+ *
+ * @file core_ipiv_to_perm.c
+ *
+ * @copyright 2023-2023 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
+ *                      Univ. Bordeaux. All rights reserved.
+ *
+ ***
+ *
+ * @brief Chameleon core_ipiv_to_perm CPU kernel
+ *
+ * @version 1.3.0
+ * @author Mathieu Faverge
+ * @date 2023-08-31
+ */
+#include "coreblas.h"
+
+/**
+ *******************************************************************************
+ *
+ * The idea here is to generate a permutation from the sequence of
+ * pivot.  To avoid storing one whole column at each step, we keep
+ * track of two vectors of nb elements, the first one contains the
+ * permutation of the first nb elements, and the second one contains
+ * the inverse permutation of those same elements.
+ *
+ * Lets have i the element to pivot with ip. ipiv[i] = ip;
+ * We set i_1 as such invp[ i_1  ] = i
+ *  and  ip_1 as such invp[ ip_1 ] = ip
+ *
+ * At each step we want to:
+ *   - swap perm[i] and perm[ip]
+ *   - set invp[i_1] to ip
+ *   - set invp[ip_1] to i
+ *
+ *******************************************************************************
+ *
+ * @param[in] m0
+ *          The base index for all values in ipiv, perm and invp. m0 >= 0.
+ *
+ * @param[in] m
+ *          The number of elements in perm and invp. m >= 0.
+ *
+ * @param[in] k
+ *          The number of elements in ipiv. k >= 0.
+ *
+ * @param[in] ipiv
+ *          The pivot array of size n. This is a (m0+1)-based indices array to follow
+ *          the Fortran standard.
+ *
+ * @param[out] perm
+ *          The permutation array of the destination row indices (m0-based) of the [1,n] set of rows.
+ *
+ * @param[out] invp
+ *          The permutation array of the origin row indices (m0-based) of the [1,n] set of rows.
+ *
+ */
+void CORE_ipiv_to_perm( int m0, int m, int k, int *ipiv, int *perm, int *invp )
+{
+    int i, j, ip;
+    int i_1, ip_1;
+
+    for(i=0; i < m; i++) {
+        perm[i] = i + m0;
+        invp[i] = i + m0;
+    }
+
+    for(i = 0; i < k; i++) {
+        ip = ipiv[i]-1;
+        assert( ip - m0 >= i );
+
+        if ( ip - m0 > i ) {
+
+            i_1 = perm[i];
+
+            if (ip-m0 < m) {
+                ip_1 = perm[ip-m0];
+                perm[ip-m0] = i_1;
+            } else {
+                ip_1 = ip;
+                for(j=0; j < m; j++) {
+                    if( invp[j] == ip ) {
+                        ip_1 = j + m0;
+                        break;
+                    }
+                }
+            }
+
+            perm[i] = ip_1;
+            i_1  -= m0;
+            ip_1 -= m0;
+
+            if (i_1  < m) invp[i_1 ] = ip;
+            if (ip_1 < m) invp[ip_1] = i + m0;
+        }
+    }
+}

coreblas/include/coreblas.h

@@ -11,14 +11,14 @@
  *
  * @brief Chameleon CPU kernels main header
  *
- * @version 1.2.0
+ * @version 1.3.0
  * @author Jakub Kurzak
  * @author Hatem Ltaief
  * @author Florent Pruvost
  * @author Guillaume Sylvand
  * @author Mathieu Faverge
  * @author Raphael Boucherie
- * @date 2022-02-22
+ * @date 2023-08-31
  *
  */
 #ifndef _coreblas_h_
@@ -87,6 +87,8 @@ void __coreblas_kernel_trace( const char *func, ... );
 
 #endif
 
+void CORE_ipiv_to_perm( int m0, int m, int k, int *ipiv, int *perm, int *invp );
+
 END_C_DECLS
 
 #endif /* _coreblas_h_ */