#include <stdio.h>
#include <stdlib.h>
#include <stddef.h>
#include <stdbool.h>
void
printDGeMatrix(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
printIGeMatrix(size_t m, size_t n,
const size_t *A,
ptrdiff_t incRowA, ptrdiff_t incColA)
{
for (size_t i=0; i<m; ++i) {
for (size_t j=0; j<n; ++j) {
printf("%6zu ", A[i*incRowA + j*incColA]);
}
printf("\n");
}
printf("\n");
}
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
dlaswp(size_t n, size_t k0, size_t k1,
const size_t *p, ptrdiff_t incP,
double *A, ptrdiff_t incRowA, ptrdiff_t incColA)
{
ptrdiff_t kInc = k0<k1 ? 1 : -1;
k1 += kInc;
for (size_t k=k0; k!=k1; k+=kInc) {
if (k!=p[k*incP]) {
dswap(n, &A[k*incRowA], incColA, &A[p[k*incP]], incColA);
}
}
}
#ifndef COLMAJOR
#define COLMAJOR 1
#endif
int
main()
{
printf("COLMAJOR = %d\n", COLMAJOR);
size_t m = 6;
size_t n = 2;
ptrdiff_t incRowA = COLMAJOR ? 1 : n;
ptrdiff_t incColA = COLMAJOR ? m : 1;
double *A = malloc(m*n*sizeof(double));
size_t k = 4;
ptrdiff_t incP = 1;
size_t *p = malloc(k*sizeof(size_t));
if (!A || !p) {
abort();
}
for (size_t j=0; j<n; ++j) {
for (size_t i=0; i<m; ++i) {
A[i*incRowA + j*incColA] = i*n + j + 1;
}
}
size_t k0 = 2;
size_t k1 = 3;
p[0*incP] = -1;
p[1*incP] = -1;
p[2*incP] = 5;
p[3*incP] = 4;
printf("A =\n");
printDGeMatrix(m, n, A, incRowA, incColA);
printf("p =\n");
printIGeMatrix(1, k, p, 0, incP);
printf("dlaswp with k0 = %zu, k1 = %zu\n", k0, k1);
dlaswp(n, k0, k1, p, incP, A, incRowA, incColA);
printf("A =\n");
printDGeMatrix(m, n, A, incRowA, incColA);
printf("dlaswp with k0 = %zu, k1 = %zu\n", k1, k0);
dlaswp(n, k1, k0, p, incP, A, incRowA, incColA);
printf("A =\n");
printDGeMatrix(m, n, A, incRowA, incColA);
free(A);
free(p);
}