From dcc0a1214ac865f0dac58b245452904f3d7834cd Mon Sep 17 00:00:00 2001 From: jan Date: Sat, 26 Mar 2016 15:29:00 -0700 Subject: [PATCH] add mmwrite, mmread, and test --- mmread.m | 181 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ mmtest.m | 60 ++++++++++++++++++ mmwrite.m | 170 ++++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 411 insertions(+) create mode 100644 mmread.m create mode 100644 mmtest.m create mode 100644 mmwrite.m diff --git a/mmread.m b/mmread.m new file mode 100644 index 0000000..04f9b6d --- /dev/null +++ b/mmread.m @@ -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 \ No newline at end of file diff --git a/mmtest.m b/mmtest.m new file mode 100644 index 0000000..616f0f8 --- /dev/null +++ b/mmtest.m @@ -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 + + + + + + + + diff --git a/mmwrite.m b/mmwrite.m new file mode 100644 index 0000000..7569ba6 --- /dev/null +++ b/mmwrite.m @@ -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);