cpp-cheat/lapack/c.c

/*
How to use blas and lapack with the standard
interfaces provided by their respective projects, repectively through
`cblas.h` and `lapacke.h`
*/

#include <assert.h>
#include <math.h>
#include <stdbool.h>
#include <stdio.h>
#include <stdlib.h>

#include <cblas.h>
#include <lapacke.h>

/**
* assert two integers are equal
* if not, print them to stderr and assert false
*/
void assert_eqi( int i1, int i2 )
{
    if( i1 != i2 )
    {
        fprintf( stderr, "%d\n%d\n\n\n", i1, i2 );
        assert( false );
    }
}

/*
* assert two doubles are equal within err precision
* if not, print them to stderr
*/
void assert_eqd( double d1, double d2, double err )
{
    if( fabs( d1 - d2 ) > err )
    {
        fprintf( stderr, "%f\n%f\n\n\n", d1, d2 );
        assert( false );
    }
}

/*
* print an array of doubles to stderr
*/
void print_vecd( int n, double * v )
{
    int i;
    for(i=0; i<n; i++)
    {
        fprintf( stderr, "%.2f ", v[i] );
    }
    printf("\n");
}

/*
bool eq_vecd2( int n, double * v1, double * v2, double err)
{
    if( fabs( d1 - d2 ) > err )
        return false
    return true;
}

void assert_eqvd( int n, double * v1, double * v2, double err){
    if( eq_vecd( n, v1, v2, err ) != true ){
        print_vecd( n, v1 );
        print_vecd( n, v2 );
    }
}
*/

int main(void)
{
    int info, ipiv2[2];
    float err = 1e-6;
    float x2[2], b2[2], c2[2];
    float a2x2[2][2];

    //#cblas

        x2[0] =  1.0;
        x2[1] = -2.0;
        assert_eqd( cblas_snrm2(2, x2, 1), sqrt(5.0), err );

    //#lapacke

        //1 2 x = 5
        //3 4 y   11

        //x = 1
        //y = 1

        a2x2[0][0] = 1.0;
        a2x2[1][0] = 2.0;
        a2x2[0][1] = 3.0;
        a2x2[1][1] = 4.0;
        b2[0] = 5.;
        b2[1] = 11.;

        info = LAPACKE_sgesv(
            LAPACK_COL_MAJOR, //COL or ROW
            2,                //n
            1,                //nrhs
            &a2x2[0][0],
            2,                //lda
            &ipiv2[0],
            &b2[0],
            2                 //ldb
        );
        c2[0] = 1.0;
        c2[1] = 2.0;
        assert_eqi( info, 0 );
        assert_eqd( b2[0], c2[0], err );
        assert_eqd( b2[1], c2[1], err );

    return EXIT_SUCCESS;
}