Code:
. M = N + K
DO 90 J = N, 1, -1
D = 0.
DO 60 I = 1, M
F(I) = W(J, I) * V(I)
60 D = D + F(I) * W(J, I)
U(J, J) = D
IF ( J = 1 ) GO TO 90
IF ( D = 0. ) Go TO 80
λ = 1. / D
DO 75 I = 1, J - 1
β = 0.
DO 65 L = 1, M
65 β = β + W(I, L) * F(L)
β = λ * β
U(I, J) = β
DO 75 L = 1, M
75 W(I, L) = W(I, L) - β * W(J, L)
GO TO 90
80 DO 85 I = 1, J-1
85 U(I, J) = 0.
90 CONTINUE
EDIT: Here it is with 0-based arrays and a driver.
Code:
#include <stdio.h>
#define N 5 // max n value
#define K 2 // max k value
#define M (N + K)
void func( int n, int k,
double w[][M], double v[], double f[], double u[][N] )
{
int m = n + k;
for ( int j = n; j-- > 0; )
{
double d = 0.0;
for ( int i = 0; i < m; ++i )
{
f[i] = w[j][i] * v[i];
d += f[i] * w[j][i];
}
u[j][j] = d;
if ( j == 0 )
break;
if ( d == 0.0 )
for (int i = 0; i < j - 1; ++i)
u[i][j] = 0;
else
{
double lambda = 1.0 / d;
for ( int i = 0; i < j - 1; ++i )
{
double beta = 0.0;
for ( int x = 0; x < m; ++x )
beta += w[i][x] * f[x];
beta *= lambda;
u[i][j] = beta;
for ( int x = 0; x < m; ++x )
w[i][x] -= beta * w[j][x];
}
}
}
}
void prn( double *a, int size )
{
for ( int m = 0; m < size; ++m )
printf("%6.2f ", a[m]);
putchar('\n');
}
void prn2d( double *a, int rows, int cols )
{
for ( int r = 0; r < rows; ++r )
{
prn(a, cols);
a += cols;
}
}
int main()
{
double w[N][M] =
{
{ 1, 2, 3, 4, 5, 6, 7 },
{ 7, 6, 5, 4, 3, 2, 1 },
{ 1, 1, 1, 1, 1, 1, 1 },
{ 2, 2, 2, 2, 2, 2, 2 },
{ 3, 3, 3, 3, 3, 3, 3 }
};
double v[M] =
{ 1, 2, 3, 4, 5, 6, 7 };
double f[M] =
{ 7, 6, 5, 4, 3, 2, 1 };
double u[N][N] =
{
{ 1, 2, 3, 4, 5 },
{ 1, 2, 3, 4, 5 },
{ 1, 2, 3, 4, 5 },
{ 1, 2, 3, 4, 5 },
{ 1, 2, 3, 4, 5 }
};
func(N, K, w, v, f, u);
printf("w:\n");
prn2d((double*)w, N, M);
printf("u:\n");
prn2d((double*)u, N, N);
printf("v:\n");
prn(v, M);
printf("f:\n");
prn(f, M);
return 0;
}