QR_decomp.h

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

inline void QR_dec(double *A, double *Q, double *R, int rows, int cols) {
	// The function decomposes the input matrix A into the matrices Q and R: one
	// simmetric, one orthonormal and one upper triangular, by using the
	// Gram-Schmidt method. The input matrice A is defined as A[rows][cols], so
	// are the output matrices Q and R. This function is meant to be used in the
	// polar decomposition algorithm and has been tested with different sizes of
	// input matrices. Supposedly the algorithm works with any matrix, as long as
	// the columns vectors are independent.
	//
	// For tests for the standalone function please refer to the original github
	// repo this has been developed in.
	// https://github.com/DavidePatria/QR_decomposition_C/blob/main/README.md

	// As already mentioned in the README the matrices orders are: A mxn => Q mxn
	// , R nxn and rank(A) must be n The matrix A[m x n] = [A_00, A_01, ... A_0n;
	// ...... ; A_m0, ... , A_mn] can be accessed as a vector that has all its
	// rows consecutively written in a long vector, even if passed as a *A and
	// defined as A[m][n].
	//
	// If A(mxn) has m<n the function still returs R(mxn) and Q(nxn), but it is
	// enough to get the submatrices Q(mxm) and R(mxn) as a valid decomposition.
	// This is what also octave does.

	// vectors for internal coputations
	double T[rows];
	double S[rows];
	double norm;
	int i, ii, j, jj, k, kk;
	double r;

	for (i = 0; i < cols; i++) {
		printf("\n");

		// scrolling a column and copying it
		for (ii = 0; ii < rows; ii++) {
			Q[ii * cols + i] = A[ii * cols + i];
		}

		for (j = 0; j < i; j++) {

			// copying columns into auxiliary variables
			for (jj = 0; jj < rows; jj++) {
				T[jj] = Q[cols * jj + j];
				S[jj] = A[cols * jj + i];
			}

			// temporary storing T*K in r
			r = 0;
			for (k = 0; k < rows; k++) {
				r += T[k] * S[k];
			}

			// setting R[j][i] to r
			R[cols * j + i] = r;

			for (kk = 0; kk < rows; kk++) {
				// multiplying vector T by r
				T[kk] *= r;
				// subtract T[kk] from i-th column of Q
				Q[cols * kk + i] -= T[kk];
			}
		}

		// rezeroing norm at each cycle
		norm = 0;
		// norm of the i-th column
		for (k = 0; k < rows; k++) {
			// computing norm^2
			norm += Q[cols * k + i] * Q[cols * k + i];
		}
		norm = sqrt(norm);

		// assigning i-th element of R diagonal
		R[cols * i + i] = norm;

		for (k = 0; k < rows; k++) {
			Q[cols * k + i] /= R[cols * i + i];
		}
	}
}