Gram-Schmidt for n x m matrix, n >= m

gram_schmidt.cpp

#include <cstdlib>
#include <iostream>
#include <math.h>

using namespace std;

/// matrix dimensions n >= m for my use case
const int n = 3;
const int m = 2;

/// static n x m matrix
double a[n][m] = {
    { 1.0, 0.0},
    {-1.0, 0.0},
    { 0.0, 1.0}
};

/// matrix Q and R
double r[n][m], q[n][m];

int main(int argc, char *argv[]) {
    int k, i, j;
    for (k=0; k<m; k++) {
      r[k][k]=0;
      for (i=0; i<n; i++) {
        r[k][k] = r[k][k] + a[i][k] * a[i][k]; 
      }
      
      r[k][k] = sqrt(r[k][k]); 
      
      for (i=0; i<n; i++) {
          q[i][k] = a[i][k]/r[k][k];
      }
     
      for(j=k+1; j<m; j++) {
        r[k][j]=0;
        for(i=0; i<n; i++) {
            r[k][j] += q[i][k] * a[i][j];
        }
        for (i=0; i<n; i++) { 
            a[i][j] = a[i][j] - r[k][j]*q[i][k];
        }
      }
    }
    return EXIT_SUCCESS;
}

Calculates Gram-Schmidt orthogonalization of a rectangular matrix n x m, with n > m. R matrix contains the orthogonal columns.