#include <stdio.h>
#include <string.h>
#include <math.h>

double dotprod(int m,double u[m],double v[m]){
	double r=u[0]*v[0];
	for(int i=1;i<m;i++) r+=u[i]*v[i];
	return r;
}
double l2norm(int m,double u[m]){
	return sqrt(dotprod(m,u,u));
}
double gschmidt(int m,int n,double At[n][m],
	double Qt[n][m],double Rt[n][n]){
	for(int i=0;i<n;i++){
		double tq[m];
		memcpy(tq,At[i],m*sizeof(double));
		for(int j=0;j<i;j++){
			double r=dotprod(m,Qt[j],tq);
			for(int k=0;k<m;k++)
				tq[k]-=r*Qt[j][k];
		}
		double r=l2norm(m,tq);
		for(int k=0;k<m;k++)
			Qt[i][k]=tq[k]/r;
	}
	memset(Rt,0,n*n*sizeof(double));
	for(int i=0;i<n;i++)
		for(int j=0;j<n;j++)
			for(int k=0;k<m;k++)
				Rt[j][i]+=Qt[i][k]*At[j][k];
}

void matprint(int m,int n,double At[n][m]){
	for(int j=0;j<n;j++){
		for(int i=0;i<m;i++)
			printf("%g ",At[j][i]);
		printf("\n");
	}
}

int main(){
	double u[3]={1,2,3};
	double v[3]={2,3,1};
	double r=dotprod(3,u,v);
	printf("dotprod(u,v)=%g\n",r);
	printf("norm(u)=%g\n",l2norm(3,u));
	double At[2][3]={ {1,2,3}, {2,3,1} };
	double Qt[2][3];
	double Rt[2][2];
	gschmidt(3,2,At,Qt,Rt);
	printf("Qt=\n");
	matprint(3,2,Qt);
	printf("Rt=\n");
	matprint(2,2,Rt);
	return 0;
}