#include <stdlib.h>
#include <stdio.h>
#include <stddef.h>
#include <math.h>
void
initGeMatrix(size_t m, size_t n,
double *A,
ptrdiff_t incRowA, ptrdiff_t incColA)
{
for (size_t i=0; i<m; ++i) {
for (size_t j=0; j<n; ++j) {
A[i*incRowA + j*incColA] = i*n + j +1;
}
}
}
void
printGeMatrix(size_t m, size_t n,
const double *A,
ptrdiff_t incRowA, ptrdiff_t incColA)
{
for (size_t i=0; i<m; ++i) {
for (size_t j=0; j<n; ++j) {
printf("%9.2lf ", A[i*incRowA + j*incColA]);
}
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
daxpy(size_t n, double alpha,
const double *x, ptrdiff_t incX,
double *y, ptrdiff_t incY)
{
for (size_t i=0; i<n; ++i) {
y[i*incY] += alpha*x[i*incX];
}
}
double
ddot(size_t n,
const double *x, ptrdiff_t incX,
const double *y, ptrdiff_t incY)
{
double result = 0;
for (size_t i=0; i<n; ++i) {
result += x[i*incX]*y[i*incY];
}
return result;
}
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] *= alpha;
}
} else {
for (size_t i=0; i<n; ++i) {
x[i*incX] = 0;
}
}
}
#ifndef COLMAJOR
#define COLMAJOR 1
#endif
int
main()
{
printf("COLMAJOR = %d\n", COLMAJOR);
size_t m = 5, n = 10;
ptrdiff_t incRowA = COLMAJOR ? 1 : n;
ptrdiff_t incColA = COLMAJOR ? m : 1;
double *A = malloc(m*n*sizeof(double));
if (!A) {
abort();
}
initGeMatrix(m, n, A, incRowA, incColA);
printf("A =\n");
printGeMatrix(m, n, A, incRowA, incColA);
printf("ddot(A(:,1), A(:,4)) = %lf\n",
ddot(m, &A[1*incColA], incRowA, &A[4*incColA], incRowA));
printf("compute: A(:,1) -= 0.5*A(:,4)\n");
daxpy(m, -0.5, &A[4*incColA], incRowA, &A[1*incColA], incRowA);
printf("A =\n");
printGeMatrix(m, n, A, incRowA, incColA);
printf("copy: A(1,:) = A(2,:)\n");
dcopy(n, &A[2*incRowA], incColA, &A[1*incRowA], incColA);
printf("A =\n");
printGeMatrix(m, n, A, incRowA, incColA);
printf("scal: A(1,:) = 1.5*A(1,:)\n");
dscal(n, 1.5, &A[1*incRowA], incColA);
printf("A =\n");
printGeMatrix(m, n, A, incRowA, incColA);
A[1*incRowA+0*incColA] = nan("");
A[1*incRowA+3*incColA] = nan("");
A[1*incRowA+9*incColA] = nan("");
printf("\n\nAfter inserting some NaNs into A:\n");
printf("A =\n");
printGeMatrix(m, n, A, incRowA, incColA);
printf("scal: A(1,:) = 1.5*A(1,:) <- will give wrong results\n");
dscal(n, 1.5, &A[1*incRowA], incColA);
printf("A =\n");
printGeMatrix(m, n, A, incRowA, incColA);
printf("scal: A(1,:) = 0*A(1,:)\n");
dscal(n, 0.0, &A[1*incRowA], incColA);
printf("A =\n");
printGeMatrix(m, n, A, incRowA, incColA);
free(A);
}