Fortran:

parameter(nmax=100); dimension a(nmax,nmax)
n=2; a(1,1)=1; a(2,1)=2; a(1,2)=3; a(2,2)=4;
call my_routine(a,n,nmax)

C:

float **a; int n=2; a = (float **)malloc(n*sizeof(float *));
for ( i=0; i<n ; i++ ) a[i] = (float *)malloc(n*sizeof(float));
a[0][0]=1; a[0][1]=2; a[1][0]=3; a[1][1]=4;
my_function(a,n);

consider a matrix A which consists of n columns (a1,a2,...,an).
create a matrix Q = (q1,q2,...,qn) in the following way:
cycle over all columns of the matrix A : do i=1,n
take column qi as normalised ai : qi = ai/Ri,i ; Ri,i=√(ai·ai)
make the remaining columns orthogonal to qi :
	do j=i+1,n
		Ri,j = qi·aj
		aj = aj - qiRi,j
	end do
end do
now, apparently, QTA = R, or, A = QR