add mmwrite, mmread, and test

This commit is contained in:
jan 2016-03-26 15:29:00 -07:00
parent e4eb232b1d
commit dcc0a1214a
3 changed files with 411 additions and 0 deletions

181
mmread.m Normal file
View File

@ -0,0 +1,181 @@
function [A,rows,cols,entries,rep,field,symm] = mmread(filename)
%
% function [A] = mmread(filename)
%
% function [A,rows,cols,entries,rep,field,symm] = mmread(filename)
%
% Reads the contents of the Matrix Market file 'filename'
% into the matrix 'A'. 'A' will be either sparse or full,
% depending on the Matrix Market format indicated by
% 'coordinate' (coordinate sparse storage), or
% 'array' (dense array storage). The data will be duplicated
% as appropriate if symmetry is indicated in the header.
%
% Optionally, size information about the matrix can be
% obtained by using the return values rows, cols, and
% entries, where entries is the number of nonzero entries
% in the final matrix. Type information can also be retrieved
% using the optional return values rep (representation), field,
% and symm (symmetry).
%
mmfile = fopen(filename,'r');
if ( mmfile == -1 )
disp(filename);
error('File not found');
end;
header = fgets(mmfile);
if (header == -1 )
error('Empty file.')
end
% NOTE: If using a version of Matlab for which strtok is not
% defined, substitute 'gettok' for 'strtok' in the
% following lines, and download gettok.m from the
% Matrix Market site.
[head0, header] = strtok(header); % see note above
[head1, header] = strtok(header);
[rep, header] = strtok(header);
[field, header] = strtok(header);
[symm, header] = strtok(header);
head1 = lower(head1);
rep = lower(rep);
field = lower(field);
symm = lower(symm);
if isempty(symm)
disp(['Not enough words in header line of file ',filename])
disp('Recognized format: ')
disp('%%MatrixMarket matrix representation field symmetry')
error('Check header line.')
end
if ~strcmp(head0,'%%MatrixMarket')
error('Not a valid MatrixMarket header.')
end
if ~strcmp(head1,'matrix')
disp(['This seems to be a MatrixMarket ',head1,' file.']);
disp('This function only knows how to read MatrixMarket matrix files.');
disp(' ');
error(' ');
end
% Read through comments, ignoring them
commentline = fgets(mmfile);
while ~isempty(commentline) && commentline(1) == '%',
commentline = fgets(mmfile);
end
% Read size information, then branch according to
% sparse or dense format
switch rep,
case 'coordinate', % read matrix given in sparse coordinate format
[sizeinfo,count] = sscanf(commentline, '%d%d%d');
while count == 0,
commentline = fgets(mmfile);
if commentline == -1,
error('End-of-file reached before size information was found.')
end
[sizeinfo,count] = sscanf(commentline, '%d%d%d');
if count > 0 && count ~= 3,
error('Invalid size specification line.')
end
end
rows = sizeinfo(1);
cols = sizeinfo(2);
entries = sizeinfo(3);
T = fscanf(mmfile, '%f');
switch field,
case 'real', % real valued entries:
check_size(T, 3 * entries);
T = reshape(T(:), 3, entries)';
A = sparse(T(:,1), T(:,2), T(:,3), rows , cols);
case 'complex', % complex valued entries:
check_size(T, 4 * entries);
T = reshape(T(:), 4, entries)';
A = sparse(T(:,1), T(:,2), T(:,3) + T(:,4)*1i, rows , cols);
case 'pattern', % pattern matrix (no values given):
check_size(T, 2 * entries);
T = reshape(T(:), 2, entries)';
A = sparse(T(:,1), T(:,2), ones(entries, 1), rows , cols);
end
case 'array', % read matrix given in dense array (column major) format
[sizeinfo, count] = sscanf(commentline, '%d%d');
while count == 0
commentline = fgets(mmfile);
if commentline == -1
error('End-of-file reached before size information was found.')
end
[sizeinfo, count] = sscanf(commentline, '%d%d');
if count > 0 && count ~= 2
error('Invalid size specification line.')
end
end
rows = sizeinfo(1);
cols = sizeinfo(2);
entries = rows*cols;
A = fscanf(mmfile, '%f');
switch field,
case 'real', % real valued entries:
% do nothing
case 'complex', % complex valued entries:
A = A(1:2:end-1) + 1i * A(2:2:end);
case 'pattern', % pattern (makes no sense for dense)
disp('Matrix type:',field)
error('Pattern matrix type invalid for array storage format.');
otherwise,
disp('Matrix type:',field)
error('Invalid matrix type specification. Check header against MM documentation.');
end
if any(strcmp(symm,{'symmetric','hermitian','skew-symmetric'}))
for j=1:cols-1,
currenti = j*rows;
A = [A(1:currenti);
zeros(j,1);
A(currenti+1:length(A))];
end
elseif ~strcmp(symm,'general'),
disp('Unrecognized symmetry')
disp(symm)
disp('Recognized choices:')
disp(' symmetric')
disp(' hermitian')
disp(' skew-symmetric')
disp(' general')
error('Check symmetry specification in header.');
end
A = reshape(A,rows,cols);
end
fclose(mmfile);
%
% If symmetric, skew-symmetric or Hermitian, duplicate lower
% triangular part and modify entries as appropriate:
%
switch symm,
case 'symmetric',
A = A + A.' - diag(diag(A));
entries = nnz(A);
case 'hermitian',
A = A + A' - diag(diag(A));
entries = nnz(A);
case 'skew-symmetric',
A = A - A.';
entries = nnz(A);
end
% Done.
end
function check_size(T, n)
if size(T) ~= n
disp('Data file does not contain expected amount of data.');
disp('Check that number of data lines matches nonzero count.');
error('Invalid data.');
end
end

60
mmtest.m Normal file
View File

@ -0,0 +1,60 @@
function mmtest()
tmpfile = [num2str(now()) '.mtx'];
fp_err = 1e-15;
function works = rwtest(m)
mmwrite(tmpfile, m);
m2 = mmread(tmpfile);
works = all(abs(m(:) - m2(:)) < fp_err*mean(abs(m(:))));
delete(tmpfile);
end
dense_r = rand(100) + 1i * rand(100);
sparse_r = sparse(dense_r * (abs(dense_r) > mean(abs(dense_r(:)))));
% full
assert(rwtest(dense_r));
assert(rwtest(sparse_r));
% symmetric
assert(rwtest( tril(dense_r) + tril(dense_r).' ));
assert(rwtest( tril(sparse_r) + tril(sparse_r).' ));
% skew-symmetric
assert(rwtest( tril(dense_r) - tril(dense_r).' ));
assert(rwtest( tril(sparse_r) - tril(sparse_r).' ));
% hermitian
assert(rwtest( tril(dense_r) + tril(dense_r)' ));
assert(rwtest( tril(sparse_r) + tril(sparse_r)' ));
dense_r = real(dense_r);
sparse_r = real(sparse_r);
% full
assert(rwtest(dense_r));
assert(rwtest(sparse_r));
% symmetric
assert(rwtest( tril(dense_r) + tril(dense_r).' ));
assert(rwtest( tril(sparse_r) + tril(sparse_r).' ));
% skew-symmetric
assert(rwtest( tril(dense_r) - tril(dense_r).' ));
assert(rwtest( tril(sparse_r) - tril(sparse_r).' ));
% hermitian
assert(rwtest( tril(dense_r) + tril(dense_r)' ));
assert(rwtest( tril(sparse_r) + tril(sparse_r)' ));
end

170
mmwrite.m Normal file
View File

@ -0,0 +1,170 @@
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);