LU decomposition experiment

The implementation of http://rosettacode.org/wiki/LU_decomposition for LU decomposition is the most stable way in my home experiments. Though I changed this pivot method to other pivot procedures, their results were not enhanced.
I wrote a code for LU decomposition using cv::Mat.

(addition) This code would make an error when elements on the diagonal in u would be 0 if the error can be avoid because this pivot does not consider the elements would be 0. To avoid it, check flag must be introduced.

void pivot(cv::Mat a, cv::Mat &p)
{
 int n = a.size().height;
 p = cv::Mat::zeros(a.size(), CV_64FC1);
 for (int i = 0; i < n; i++){
  p.at<double>(i, i) = 1;
 }

 for (int i = 0; i < n; i++){
  int maxj = i;
  for (int j = i; j < n; j++){
   if (fabs(a.at<double>(j, i)) > fabs(a.at<double>(maxj, i))) maxj = j;
  }
  //swap a_i and a_maxj
  if (i != maxj){
   for (int j = 0; j < n; j++){
    double temp = p.at<double>(i, j);
    p.at<double>(i, j) = p.at<double>(maxj, j);
    p.at<double>(maxj, j) = temp;
   }
  }
 }
}
  
void decomposition(cv::Mat a, cv::Mat &l, cv::Mat &u)
{
 cv::Mat p;
 pivot(a, p);
 cv::Mat ap = p*a;
 int n = ap.size().height;
 l = cv::Mat::zeros(ap.size(), CV_64FC1);
 u = cv::Mat::zeros(ap.size(), CV_64FC1);
 //init l
 for (int i = 0; i < n; i++){
  l.at<double>(i, i) = 1;
 }
 //init u
 for (int i = 0; i < n; i++){
  for (int j = 0; j < n; j++){
   double sum = 0;
   if (j <= i){
    for (int k = 0; k < j; k++){
     sum += l.at<double>(j, k)*u.at<double>(k, i);
    }
    u.at<double>(j, i) = ap.at<double>(j, i) - sum;
   }
   sum = 0;
   if (j >= i){
    for (int k = 0; k < i; k++){
     sum += l.at<double>(j, k)*u.at<double>(k, i);
    }
    l.at<double>(j, i) = (ap.at<double>(j, i) - sum) / u.at<double>(i, i);
   }
  }
 }
}

int main(int argc, _TCHAR* argv[])
{
 cv::Mat a = cv::Mat::zeros(cv::Size(4, 4), CV_64FC1);
 a.at<double>(0, 0) = 11; a.at<double>(0, 1) = 9; a.at<double>(0, 2) = 24; a.at<double>(0, 3) = 2;
 a.at<double>(1, 0) = 1; a.at<double>(1, 1) = 5; a.at<double>(1, 2) = 2; a.at<double>(1, 3) = 6;
 a.at<double>(2, 0) = 3; a.at<double>(2, 1) = 17; a.at<double>(2, 2) = 18; a.at<double>(2, 3) = 1;
 a.at<double>(3, 0) = 2; a.at<double>(3, 1) = 5; a.at<double>(3, 2) = 7; a.at<double>(3, 3) = 1;
 cv::Mat l, u;
 decomposition(a, l, u);
 std::cout << "a:\n" << a << "\nl:\n" << l << "\nu:\n" << u << "\nlu\n" << l*u << "\n\n";
}