function [ err ] = mmwrite(filename, A, comment, mattype, precision)
%
% Function: mmwrite(filename,A,comment,mattype,precision)
%
%    Writes the sparse or dense matrix A to a Matrix Market (MM) 
%    formatted file.
%
% Required arguments: 
%
%                 filename  -  destination file
%
%                 A         -  sparse or full matrix
%
% Optional arguments: 
%
%                 comment   -  matrix of comments to prepend to
%                              the MM file.  To build a comment matrix,
%                              use str2mat. For example:
%
%                              comment = str2mat(' Comment 1' ,...
%                                                ' Comment 2',...
%                                                ' and so on.',...
%                                                ' to attach a date:',...
%                                                [' ',date]);
%                              If ommitted, a single line date stamp comment
%                              will be included.
%
%                 mattype   -  'real'
%                              'complex'
%                              'integer'
%                              'pattern'
%                              If ommitted, data will determine type.
%
%                 precision -  number of digits to display for real 
%                              or complex values
%                              If ommitted, full working precision is used.
%

if ( nargin == 5) 
  precision = 16;
elseif ( nargin == 4) 
  precision = 16;
elseif ( nargin == 3) 
  mattype = 'real'; % placeholder
  precision = 16;
elseif ( nargin == 2) 
  comment = '';
  % Check whether there is an imaginary part:
  mattype = 'real'; % placeholder
  precision = 16;
end

mmfile = fopen(filename,'w');
if ( mmfile == -1 )
 error('Cannot open file for output');
end;


[M, N] = size(A);
  
if ~strcmp(mattype,'pattern')
    if any(imag(A(:)))
        mattype = 'complex';
    else
        mattype = 'real';
    end
end

%
% Determine symmetry:
%
if issymmetric(A),
  symm = 'symmetric';
elseif issymmetric(A, 'skew'),
  symm = 'skew-symmetric';
elseif strcmp(mattype, 'complex') && ishermitian(A),
  symm = 'hermitian';
else
  symm = 'general';
end

%%%%%%%%%%%%%       This part for sparse matrices     %%%%%%%%%%%%%%%%
if issparse(A),

  if strcmp(symm, 'general')
      [I, J, V] = find(A);
      NZ = length(V);
  else
      ATEMP = tril(A);
      [I, J, V] = find(ATEMP);
      NZ = nnz(ATEMP);
  end

% Sparse coordinate format:

  rep = 'coordinate';

  fprintf(mmfile,'%%%%MatrixMarket matrix %s %s %s\n', rep, mattype, symm);
  [MC, ~] = size(comment);
  if ( MC == 0 )
    fprintf(mmfile,'%% Generated %s\n',[date]);
  else
    for i=1:MC,
      fprintf(mmfile,'%%%s\n',comment(i,:));
    end
  end
  fprintf(mmfile,'%d %d %d\n',M,N,NZ);
  
  switch mattype,
      case 'real',
          realformat = sprintf('%%d %%d %% .%dg\n',precision);
          fprintf(mmfile, realformat, [I, J, V].');
      case 'complex',
          cplxformat = sprintf('%%d %%d %% .%dg %% .%dg\n', precision, precision);
          fprintf(mmfile, cplxformat, [I, J, real(V), imag(V)].');
      case 'pattern',
          fprintf(mmfile, '%d %d\n', [I, J].');
      otherwise,
          err = -1;
          disp('Unsupported mattype:');
          disp(mattype);
  end

%%%%%%%%%%%%%       This part for dense matrices      %%%%%%%%%%%%%%%%
else

% Dense array format:

  rep = 'array';
  [MC, ~] = size(comment);
  fprintf(mmfile,'%%%%MatrixMarket matrix %s %s %s\n', rep, mattype, symm);
  for i=1:MC,
     fprintf(mmfile,'%%%s\n', comment(i,:));
  end;
  fprintf(mmfile,'%d %d\n',M,N);
  
  if ( ~ strcmp(symm,'general') )
     j_bool = 1;
  else 
     j_bool = 0;
  end
  rowloop = @(j) j_bool * j + (1-j_bool);
  
  switch mattype,
      case 'real',
         realformat = sprintf('%% .%dg\n', precision);
         
         for j=1:N
           for i=rowloop(j):M
              fprintf(mmfile,realformat,A(i,j));
           end
         end
      case 'complex',
         cplxformat = sprintf('%% .%dg %% .%dg\n', precision, precision);
         for j=1:N
           for i=rowloop(j):M
              fprintf(mmfile,cplxformat,real(A(i,j)),imag(A(i,j)));
           end
         end
      case 'pattern',
         err = -2;
         disp('Pattern type inconsistant with dense matrix')
      otherwise,
         err = -2;
         disp('Unknown matrix type:')
         disp(mattype);
  end
end

fclose(mmfile);