Compiler Explorer

Source code

#include <cstdio>
#include <cmath> // hypot, pow

// As of Intel oneAPI MKL version 2024.2, the direct call feature
// is only available for the following procedures:
//
//   BLAS: ?gemm, ?gemm3m, ?syrk, ?trsm, ?axpy, and ?dot
//   LAPACK: ?getrf, ?getrs, ?getri, ?potrf, and ?geqrf. (available for C applications only)
//
// More information can be found at:
//    https://www.intel.com/content/www/us/en/docs/onemkl/developer-guide-linux/2024-2/limitations-of-the-direct-call.html
//

// TODO:
//   - implement matrix padding, leading dimension should be a multiple of 8
//   - investigate use of aligned memory

template<int q>
inline void eval_basis(int n, 
                       const double * __restrict x, 
                       const double * __restrict y, 
                       double xc, 
                       double yc, 
                             double * __restrict b)
{
    #pragma omp simd
    for (int k = 0; k < n; k++) {
        double r = std::hypot(xc - x[k], yc - y[k]);
        b[k] = std::pow(r, q);
    }

// Evaluate polynomial part
    b[n+0]  = 1.0;
    b[n+1]  = xc;
    b[n+2]  = yc;
    b[n+3]  = (xc*xc);
    b[n+4]  = xc*yc;
    b[n+5]  = (yc*yc);
    b[n+6]  = (xc*xc) * xc;
    b[n+7]  = (xc*xc) * yc;
    b[n+8]  = xc * (yc*yc);
    b[n+9]  = yc * (yc*yc);
    b[n+10] = (xc*xc) * (xc*xc);
    b[n+11] = ((xc*xc) * xc) * yc;
    b[n+12] = (xc*xc) * (yc*yc);
    b[n+13] = xc * (yc *(yc*yc));
    b[n+14] = (yc*yc) * (yc*yc);

}

// Perform the transpose B = A^T
// 
// where A is of size m-by-n, and consequently B is of size n-by-m
// data layout is column-major with a leading dimension of lda or ldb,
// respectively
static void transpose(int m, int n, 
                      double * __restrict A, int lda, 
                      double * __restrict B, int ldb) 
{
#define ida(i,j) (i + j*lda)
#define idb(i,j) (i + j*ldb)

for (int j = 0; j < n; j++) {
        for (int i = 0; i < m; i++) {
            // A is read according to its natural order (column-major)
            B[idb(j,i)] = A[ida(i,j)];
        }
    }

#undef ida
#undef idb
}

// Zero an array A
// 
// where A is of size m-by-n
// data layout is column-major with a leading dimension of lda
static void zerofill(int m, int n, double *A, int lda) 
{
#define ida(i,j) (i + j*lda)
    for (int j = 0; j < n; j++) {
        for (int i = 0; i < m; i++) {
            A[ida(i,j)] = 0;
        }
    }
#undef ida
}

template<int q>
void AP(int n, const double *x, const double *y, double *M, int ld)
{
    const int np = 15; // Number of polynomial terms

#define idx(i,j) (i + ld*j)

//   v----------- set the left upper and lower blocks
    // [ A    * ]
    // [ P^T  * ]
    //
    // Loop over columns 
    for(int k = 0; k < n; k++) {
        eval_basis<q>(n,x,y,x[k],y[k],&M[idx(0,k)]);
    }

//    [ A    P ] <-- set the block P = (P^T)^T
    //    [ P^T  * ]
    transpose(np, n, &M[idx(n,0)], ld, &M[idx(0,n)], ld);

// Zero-block
    //    [ A    P ]
    //    [ P^T  0 ] <-- zero the bottom-right block
    zerofill(np, np, &M[idx(n,n)],ld);

#undef idx
}

#if DEBUG
// Helper functions for printing
extern "C" void disp(
        const int *n, 
        const int *ld, 
        const int *np, 
        const double *wrk );

extern "C" void bdisp(
        const int *n, 
        const int *ld, 
        const double *coeffs );
#endif

template<int q>
void rbf_poly4( int n, 
                const double *x,
                const double *y,
                int nc,
                const double *xc,
                const double *yc,
                double *coeffs,
                int ld,
                double *wrk,
                int *ipiv,
                int *ierr ) 
{
#define idx(i,j) (i + ld*j)
 
    const int np = 15;
    const int nt = n + np;

//
    // Evaluate matrix M
    //
    AP<q>(n,x,y,wrk,ld);

//
    // Evaluate right-hand sides b
    //
    for(int k = 0; k < nc; k++) {
        eval_basis<q>(n,x,y,xc[k],yc[k],&coeffs[idx(0,k)]);
    }

#undef idx
}

// 
// Export routines callable from C and Fortran
//
extern "C" {

void c_rbf3_poly4(const int *n, 
                const double *x,
                const double *y,
                const int *nc,
                const double *xc,
                const double *yc,
                double *coeffs,
                const int *ld,
                double *wrk,
                int *ipiv,
                int *ierr ) 
{
    rbf_poly4<3>(*n,x,y,*nc,xc,yc,coeffs,*ld,wrk,ipiv,ierr);
}

}#include <cstdio>
#include <cmath> // hypot, pow

// TODO:
//   - implement matrix padding, leading dimension should be a multiple of 8
//   - investigate use of aligned memory

}

for (int j = 0; j < n; j++) {
        for (int i = 0; i < m; i++) {
            // A is read according to its natural order (column-major)
            B[idb(j,i)] = A[ida(i,j)];
        }
    }

#undef ida
#undef idb
}

template<int q>
void AP(int n, const double *x, const double *y, double *M, int ld)
{
    const int np = 15; // Number of polynomial terms

#define idx(i,j) (i + ld*j)

//    [ A    P ] <-- set the block P = (P^T)^T
    //    [ P^T  * ]
    transpose(np, n, &M[idx(n,0)], ld, &M[idx(0,n)], ld);

// Zero-block
    //    [ A    P ]
    //    [ P^T  0 ] <-- zero the bottom-right block
    zerofill(np, np, &M[idx(n,n)],ld);

#undef idx
}

#if DEBUG
// Helper functions for printing
extern "C" void disp(
        const int *n, 
        const int *ld, 
        const int *np, 
        const double *wrk );

extern "C" void bdisp(
        const int *n, 
        const int *ld, 
        const double *coeffs );
#endif

//
    // Evaluate matrix M
    //
    AP<q>(n,x,y,wrk,ld);

//
    // Evaluate right-hand sides b
    //
    for(int k = 0; k < nc; k++) {
        eval_basis<q>(n,x,y,xc[k],yc[k],&coeffs[idx(0,k)]);
    }

#undef idx
}

// 
// Export routines callable from C and Fortran
//
extern "C" {

}