#include <stdio.h>
#include <stdlib.h>
#include <stddef.h>
#include <math.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("%9zu ", 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];
}
}
}
size_t
idamax(size_t n, const double *x, ptrdiff_t incX)
{
double max = 0;
size_t I = 0;
for (size_t i=0; i<n; ++i) {
if (fabs(x[i*incX])>max) {
I = i;
max = fabs(x[i*incX]);
}
}
return I;
}
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;
}
}
int
main()
{
size_t m = 3;
size_t n = 4;
size_t min_mn = m<n ? m : n;
ptrdiff_t incRowA = 1;
ptrdiff_t incColA = m;
ptrdiff_t incP = 1;
double *A = malloc(m*n*sizeof(double));
size_t *p = malloc(min_mn*sizeof(size_t));
if (!A || !p) {
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) {
size_t pj = j + idamax(m-j, &A[j*incRowA + j*incColA], incRowA);
if (j!=pj) {
dswap(n, &A[j*incRowA], incColA, &A[pj*incRowA], incColA);
}
p[j*incP] = pj;
double alpha = A[j*incRowA + j*incColA];
if (alpha==0) {
return j;
} else {
alpha = 1/alpha;
}
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);
}
printDGeMatrix(m, n, A, incRowA, incColA);
printIGeMatrix(1, k, p, 0, incP);
free(A);
free(p);
}