Weitere (einfache) BLAS Funktionen aufbohren

Analog zum Kopieren von Matrizen sollen folgende low-level BLAS Funktionen in hpc/ulmblas/ angepasst werden:

gescal zum Skalieren einer Matrix: \(X \leftarrow \alpha\,X\)
geaxpy zum Aufdatieren einer Matrix: \(Y \leftarrow \alpha\,X + Y\)

Dabei soll wieder berücksichtigt werden, ob die Matrix zeilen- oder spaltenweise im Speicher liegt. Natürlich sollte auch berücksichtigt werden, dass in bestimmten Spezialfällen keine Operation durchgeführt werden muss.

Beim Testen sollen wieder verschiedene Datentypen für A, B und alpha verwendet werden. Dabei kann folgende Vorlage benutzt werden:

session12/page01/hpc/tests/gescal_geaxpy.cc

#include <cstdio>
#include <hpc/ulmblas/geaxpy.h>
#include <hpc/ulmblas/gecopy.h>
#include <hpc/ulmblas/gescal.h>

template <typename T, typename Func>
void
geinit(long m, long n, T *A, long incRowA, long incColA, Func func)
{
    for (long i=0; i<m; ++i) {
        for (long j=0; j<n; ++j) {
            A[i*incRowA+j*incColA] = func(m, n, i, j);
        }
    }
}

void
print_value(float x)
{
    printf(" %6.1f", x);
}

void
print_value(double x)
{
    printf(" %6.1lf", x);
}

template <typename T>
void
geprint(long m, long n, const T *A, long incRowA, long incColA)
{
    for (long i=0; i<m; ++i) {
        for (long j=0; j<n; ++j) {
            print_value(A[i*incRowA+j*incColA]);
        }
        printf("\n");
    }
}

int
main()
{
    // Element type for A, B and alpha
    typedef float   TA;
    typedef double  TB;
    typedef double  Alpha;

    // Dimensions of A and B 
    long m       = 7;
    long n       = 8;

    // Storage order of A
    long incRowA = 1;
    long incColA = m;

    // Storage order of B
    long incRowB = n;
    long incColB = 1;

    // Allocate A and B
    TA *A = new TA[m*n];
    TB *B = new TB[m*n];

    // Value for alpha
    Alpha alpha = 1.5;
    printf("alpha =\n");
    geprint(1, 1, &alpha, 1, 1);

    // Initialize A
    auto fooInit = [](long m, long n, long i, long j) -> long
                   {
                       return m*j+i+1;
                   };
    geinit(m, n, A, incRowA, incColA, fooInit);

    // Print A
    printf("A = \n");
    geprint(m, n, A, incRowA, incColA);

    // Copy A to B
    hpc::ulmblas::gecopy(m, n, A, incRowA, incColA, B, incRowB, incColB);
    printf("B = \n");
    geprint(m, n, B, incRowB, incColB);

    //  B <- alpha*B
    printf("B <- alpha*B\n");
    printf("B = \n");
    hpc::ulmblas::gescal(m, n, alpha, B, incRowB, incColB);
    geprint(m, n, B, incRowB, incColB);

    //  B <- alpha*A + B
    printf("B <- alpha*A +B\n");
    printf("B = \n");
    hpc::ulmblas::geaxpy(m, n, alpha, A, incRowA, incColA, B, incRowB, incColB);
    geprint(m, n, B, incRowB, incColB);
}