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
     142
     143
     144
     145
     146
     147
     148
     149
     150
     151
     152
     153
     154
     155
     156
     157
     158
     159
     160
     161
     162
     163
     164
     165

#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;
}

//-- BLAS Level 2 --------------------------------------------------------------

void
dger(size_t m, size_t n, double alpha,
     const double *x, ptrdiff_t incX,
     const double *y, ptrdiff_t incY,
     double *A, ptrdiff_t incRowA, ptrdiff_t incColA)
{
    // ***Your code here***
}

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);

    // Matrix A refers to A_(1:m-1, 1:n-1)
    size_t    m       = m_ - 1;
    size_t    n       = n_ - 1;
    double    *A      = 0; /* ***Fix me *** */
    ptrdiff_t incRowA = 0; /* ***Fix me *** */
    ptrdiff_t incColA = 0; /* ***Fix me *** */

    // Vector x refers to A_(0, 1:n-1)
    double    *x      = 0; /* ***Fix me *** */
    ptrdiff_t incX    = 0; /* ***Fix me *** */

    // Vector y refers to A_(1:m-1, 0)
    double    *y      = 0; /* ***Fix me *** */
    ptrdiff_t incY    = 0; /* ***Fix me *** */

    // print A, x, y (for saving space: display vectors as row vectors)
    printf("A=\n");
    printMatrix(m, n, A, incRowA, incColA);
    printf("x=\n");
    printMatrix(1, n, x, 1, incX);
    printf("y=\n");
    printMatrix(1, m, y, 1, incY);

    // Some simple tests for the BLAS Level 2 functions

    // 1) dger
    printf("1.5*y * x^T + A -> A\n");
    // ***Your code here***
    printMatrix(m, n, A, incRow, incCol);

    // 2) dger
    printf("1.5*x * y^T + A^T -> A^T\n");
    // ***Your code here***
    printMatrix(m, n, A, incRow, incCol);

    // Deallocate matrix A_
    free(A_);
}





































































































Home