function splitLU(filenameL, filenameU, LU)
% function splitLU(filenameL, filenameU, LU)
% split LU matrix to Lower and upper triangular parts and write
% em to disk in petsc's format
% 4.2.2003 Antti.Vanne@uku.fi
  
  MAT_FILE_COOKIE = 1211216;

  % open file for lower triangle part for writing, big-endian byte ordering
  fid = fopen(filenameL, 'w', 'b');
  
  disp('Saving L matrix');
  fwrite(fid, MAT_FILE_COOKIE, 'int');
  fwrite(fid, size(LU,1), 'int');
  fwrite(fid, size(LU,2), 'int');

  nonzeros = size(LU,1)^2/2+size(LU,1)/2;
  fwrite(fid, nonzeros, 'int');

  % nonzeros in each row
  nnzrow = 1:size(LU,2);
  fwrite(fid, nnzrow, 'int');

  colInds = zeros(nonzeros,1);
  inds = [0:size(LU, 2)-1];
  
  start = 1;
  for stride = 1:size(LU,2)
    colInds(start:start+stride-1) = inds(1:stride);
    start = start+stride;
    
  end
  
  fwrite(fid, colInds, 'int');
  
  fwrite(fid, 1.0, 'double');
  for kk = 2:size(LU, 2);
    fwrite(fid, LU(kk, 1:kk-1), 'double');
    fwrite(fid, 1.0, 'double');
  end
  
  fclose(fid);
  
  %% upper triangle
  
  fid = fopen(filenameU, 'w', 'b');
  
  disp('Saving U matrix');
  fwrite(fid, MAT_FILE_COOKIE, 'int');
  fwrite(fid, size(LU,1), 'int');
  fwrite(fid, size(LU,2), 'int');

  nonzeros = size(LU,1)^2/2+size(LU,1)/2;
  fwrite(fid, nonzeros, 'int');

  % nonzeros in each row
  nnzrow = size(LU,2):-1:1;
  fwrite(fid, nnzrow, 'int');

  colInds = zeros(nonzeros,1);
  inds = [0:size(LU, 2)-1];
  
  start = 1;
  kk = 1;
  for stride = size(LU,2):-1:1
    colInds(start:start+stride-1) = inds(kk:kk+stride-1);
    start = start+stride;
    kk = kk+1;
  end
  
  fwrite(fid, colInds, 'int');
  
  for kk = 1:size(LU, 2);
    fwrite(fid, LU(kk, kk:size(LU,2)), 'double');
  end
  
  fclose(fid);


% MATRIX
%   int    MAT_FILE_COOKIE
%   int    number of rows
%   int    number of columns
%   int    total number of nonzeros
%   int    *number nonzeros in each row
%   int    *column indices of all nonzeros (starting index is zero)
%   PetscScalar *values of all nonzeros