GramSchmidt.cpph

Go to the documentation of this file.

template<int Rows, int Cols, typename Scalar>
void tarch::la::modifiedGramSchmidt (
  Matrix<Rows,Cols,Scalar>   A,
  Matrix<Rows,Cols,Scalar>&  Q,
  Matrix<Cols,Cols,Scalar>&  R
) {
  for (int i=0; i < Rows; i++){
    Q(i,0) = A(i,0);
  }
 
  for (int k=0; k < Cols; k++){
    // Compute norm of k-th column vector of (modified) A
    R(k,k) = 0.0;
    for (int i=0; i < Rows; i++){
      R(k,k) += std::pow(A(i,k), 2);
    }
    R(k,k) = std::sqrt(R(k,k));
 
    // Normalize k-th column of matrix Q
    for (int i=0; i < Rows; i++){
      Q(i,k) = A(i,k) / R(k,k);
    }
 
    // Compute entries in R and next orthonormal vector
    for (int j=k+1; j < Cols; j++){
      // Compute entries of R from current Q and A
      for (int i=0; i < Rows; i++){
        R(k,j) += Q(i,k) * A(i,j);
      }
      // Subtract contributions from computed to open orthonormal vectors
      for (int i=0; i < Rows; i++){
        A(i,j) -= Q(i,k) * R(k,j);
      }
    }
  }
}