#include <stdio.h>
#include <stdlib.h>
#include <stddef.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
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;
}
}
}
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)
{
if (alpha==0 || m==0 || n==0) {
return;
}
if (incRowA>incColA) {
dger(n, m, alpha, y, incY, x, incX, A, incColA, incRowA);
return;
}
for (size_t j=0; j<n; ++j) {
for (size_t i=0; i<m; ++i) {
A[i*incRowA+j*incColA] += alpha*x[i*incX]*y[j*incY];
}
}
}
int
main()
{
size_t m = 3;
size_t n = 4;
ptrdiff_t incRowA = 1;
ptrdiff_t incColA = m;
double *A = malloc(m*n*sizeof(double));
if (!A) {
abort();
}
A[0*incRowA + 0*incColA] = 2;
A[0*incRowA + 1*incColA] = 1;
A[0*incRowA + 2*incColA] = 3;
A[0*incRowA + 3*incColA] = 1;
A[1*incRowA + 0*incColA] = 4;
A[1*incRowA + 1*incColA] = 5;
A[1*incRowA + 2*incColA] = 6;
A[1*incRowA + 3*incColA] = 2;
A[2*incRowA + 0*incColA] = 8;
A[2*incRowA + 1*incColA] = 8;
A[2*incRowA + 2*incColA] = 9;
A[2*incRowA + 3*incColA] = 3;
printf("A =\n");
printDGeMatrix(m, n, A, incRowA, incColA);
size_t k = m<n ? m : n;
for (size_t j=0; j<k; ++j) {
double alpha = 1/A[j*incRowA + j*incColA];
dscal(m-1-j, alpha, &A[(j+1)*incRowA + j*incColA], incRowA);
dger(m-1-j, n-1-j, -1,
&A[(j+1)*incRowA + j *incColA], incRowA,
&A[ j *incRowA + (j+1)*incColA], incColA,
&A[(j+1)*incRowA + (j+1)*incColA], incRowA, incColA);
}
printf("A =\n");
printDGeMatrix(m, n, A, incRowA, incColA);
free(A);
}