Home Browse Files 
       1
       2
       3
       4
       5
       6
       7
       8
       9
      10
      11
      12
      13
      14
      15
      16
      17
      18
      19
      20
      21
      22
      23
      24
      25
      26
      27
      28
      29
      30
      31
      32
      33
      34
      35
      36
      37
      38
      39
      40
      41
      42
      43
      44
      45
      46
      47
      48
      49
      50
      51
      52
      53
      54
      55
      56
      57
      58
      59
      60
      61
      62
      63
      64
      65
      66
      67
      68
      69
      70
      71
      72
      73
      74
      75
      76
      77
      78
      79
      80
      81
      82
      83
      84
      85
      86
      87
      88
      89
      90
      91
      92
      93
      94
      95
      96
      97
      98
      99
     100
     101
     102
     103
     104
     105
     106
     107
     108
     109
     110
     111
     112
     113
     114
     115
     116
     117
     118
     119
     120
     121
     122
     123
     124
     125
     126
     127
     128
     129
     130
     131
     132
     133
     134
     135
     136
     137
     138
     139
     140
     141

#include <stddef.h>
#include <stdio.h>
#include <stdlib.h>

void
initMatrix(size_t m, size_t n,
     double *A,
     ptrdiff_t incRow, ptrdiff_t incCol)
{
    for (size_t i=0; i<m; ++i) {
        for (size_t j=0; j<n; ++j) {
            A[i*incRow+j*incCol] = i*n+j+1;
        }
    }
}

void
printMatrix(size_t m, size_t n,
            const double *A,
            ptrdiff_t incRow, ptrdiff_t incCol)
{
    for (size_t i=0; i<m; ++i) {
        for (size_t j=0; j<n; ++j) {
            printf("%10.3lf", A[i*incRow+j*incCol]);
        }
        printf("\n");
    }
    printf("\n");
}

//-- BLAS Level 1 --------------------------------------------------------------

void
dcopy(size_t n,
      const double *x, ptrdiff_t incX,
      double *y, ptrdiff_t incY)
{
    for (size_t i=0; i<n; ++i) {
        y[i*incY] = x[i*incX];
    }
}

void
dswap(size_t n,
      double *x, ptrdiff_t incX,
      double *y, ptrdiff_t incY)
{
    for (size_t i=0; i<n; ++i) {
        double tmp = x[i*incX];
        x[i*incX] = y[i*incY];
        y[i*incY] = tmp;
    }
}

void
daxpy(size_t n, double alpha,
      const double *x, ptrdiff_t incX,
      double *y, ptrdiff_t incY)
{
    if (alpha==0) {
        return;
    }
    for (size_t i=0; i<n; ++i) {
        y[i*incY] += alpha*x[i*incX];
    }
}

void
dscal(size_t n, double alpha,
      double *x, ptrdiff_t incX)
{
    if (alpha==1) {
        return;
    }
    if (alpha==0) {
        for (size_t i=0; i<n; ++i) {
            x[i*incX] = 0;
        }
    } else {
        for (size_t i=0; i<n; ++i) {
            x[i*incX] *= alpha;
        }
    }
}

double
ddot(size_t n,
     const double *x, ptrdiff_t incX,
     const double *y, ptrdiff_t incY)
{
    double alpha = 0;
    for (size_t i=0; i<n; ++i) {
        alpha += x[i*incX]*y[i*incY];
    }
    return alpha;
}

int
main()
{
    // Allocate matrix A with dimension m x n
    size_t m = 5;
    size_t n = 8;

    size_t incRow = n;
    size_t incCol = 1;

    double *A = (double *) malloc(m*n*sizeof(double));

    // Initialize matrix A
    initMatrix(m, n, A, incRow, incCol);

    // Print initial matrix A
    printf("A=\n");
    printMatrix(m, n, A, incRow, incCol);

    // Some simple tests for the BLAS Level 1 functions

    // 1) dcopy
    printf("A(:,2) -> A(:, 3)\n");
    dcopy(m, &A[2*incCol], incRow, &A[3*incCol], incRow);
    printMatrix(m, n, A, incRow, incCol);

    // 2) dswap
    printf("A(3,:) <-> A(1,:)\n");
    dswap(n, &A[3*incRow], incCol, &A[1*incRow], incCol);
    printMatrix(m, n, A, incRow, incCol);

    // 3) daxpy
    printf("2.5*A(0,:) + A(2,:) -> A(2,:)\n");
    daxpy(n, 2.5, &A[0*incRow], incCol, &A[2*incRow], incCol);
    printMatrix(m, n, A, incRow, incCol);

    // 4) dscal
    printf("0.5*A(:,5) -> A(:,5)\n");
    dscal(m, 0.5, &A[5*incCol], incRow);
    printMatrix(m, n, A, incRow, incCol);

    // Deallocate matrix A
    free(A);
}





































































































Home