#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");
}
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;
}
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()
{
size_t m = 5;
size_t n = 8;
size_t incRow = n;
size_t incCol = 1;
double *A = (double *) malloc(m*n*sizeof(double));
initMatrix(m, n, A, incRow, incCol);
printf("A=\n");
printMatrix(m, n, A, incRow, incCol);
printf("A(:,2) -> A(:, 3)\n");
dcopy(m, &A[2*incCol], incRow, &A[3*incCol], incRow);
printMatrix(m, n, A, incRow, incCol);
printf("A(3,:) <-> A(1,:)\n");
dswap(n, &A[3*incRow], incCol, &A[1*incRow], incCol);
printMatrix(m, n, A, incRow, incCol);
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);
printf("0.5*A(:,5) -> A(:,5)\n");
dscal(m, 0.5, &A[5*incCol], incRow);
printMatrix(m, n, A, incRow, incCol);
free(A);
}