线性求解器Ax=b的验证

前言

一般情况集成了一个线性求解器（Ax=b），我们需要验证其性能和精度，这时需要大量数据来做验证， 尤其是有不同性质的矩阵 $A$ 的数据，例如：稀疏性， 对称性， 正定性， 对角占优等。

Matrix Market

Matrix Market
Matrix Market 网站提供了友好的接口和数据， 方便我们验证求解器的精度和性能，尤其是对各种性质的矩阵 $A$ 都有很多数据， 包含不同维数。

Matlab IO

Read data

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 ( length(symm) == 0 )
   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 length(commentline) > 0 & commentline(1) == '%',
  commentline = fgets(mmfile);
end

% Read size information, then branch according to
% sparse or dense format

if ( strcmp(rep,'coordinate')) %  read matrix given in sparse 
                              %  coordinate matrix 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);
  
  if  ( strcmp(field,'real') )               % real valued entries:
  
    [T,count] = fscanf(mmfile,'%f',3);
    T = [T; fscanf(mmfile,'%f')];
    if ( size(T) ~= 3*entries )
       message = ...
       str2mat('Data file does not contain expected amount of data.',...
               'Check that number of data lines matches nonzero count.');
       disp(message);
       error('Invalid data.');
    end
    T = reshape(T,3,entries)';
    A = sparse(T(:,1), T(:,2), T(:,3), rows , cols);
  
  elseif   ( strcmp(field,'complex'))            % complex valued entries:
  
    T = fscanf(mmfile,'%f',4);
    T = [T; fscanf(mmfile,'%f')];
    if ( size(T) ~= 4*entries )
       message = ...
       str2mat('Data file does not contain expected amount of data.',...
               'Check that number of data lines matches nonzero count.');
       disp(message);
       error('Invalid data.');
    end
    T = reshape(T,4,entries)';
    A = sparse(T(:,1), T(:,2), T(:,3) + T(:,4)*sqrt(-1), rows , cols);
  
  elseif  ( strcmp(field,'pattern'))    % pattern matrix (no values given):
  
    T = fscanf(mmfile,'%f',2);
    T = [T; fscanf(mmfile,'%f')];
    if ( size(T) ~= 2*entries )
       message = ...
       str2mat('Data file does not contain expected amount of data.',...
               'Check that number of data lines matches nonzero count.');
       disp(message);
       error('Invalid data.');
    end
    T = reshape(T,2,entries)';
    A = sparse(T(:,1), T(:,2), ones(entries,1) , rows , cols);

  end

elseif ( strcmp(rep,'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;
  if  ( strcmp(field,'real') )               % real valued entries:
    A = fscanf(mmfile,'%f',1);
    A = [A; fscanf(mmfile,'%f')];
    if ( strcmp(symm,'symmetric') | strcmp(symm,'hermitian') | strcmp(symm,'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);
  elseif  ( strcmp(field,'complex'))         % complx valued entries:
    tmpr = fscanf(mmfile,'%f',1);
    tmpi = fscanf(mmfile,'%f',1);
    A  = tmpr+tmpi*i;
    for j=1:entries-1
      tmpr = fscanf(mmfile,'%f',1);
      tmpi = fscanf(mmfile,'%f',1);
      A  = [A; tmpr + tmpi*i];
    end
    if ( strcmp(symm,'symmetric') | strcmp(symm,'hermitian') | strcmp(symm,'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);
  elseif  ( strcmp(field,'pattern'))    % pattern (makes no sense for dense)
   disp('Matrix type:',field)
   error('Pattern matrix type invalid for array storage format.');
  else                                 % Unknown matrix type
   disp('Matrix type:',field)
   error('Invalid matrix type specification. Check header against MM documentation.');
  end
end

%
% If symmetric, skew-symmetric or Hermitian, duplicate lower
% triangular part and modify entries as appropriate:
%

if ( strcmp(symm,'symmetric') )
   A = A + A.' - diag(diag(A));
   entries = nnz(A);
elseif ( strcmp(symm,'hermitian') )
   A = A + A' - diag(diag(A));
   entries = nnz(A);
elseif ( strcmp(symm,'skew-symmetric') )
   A = A - A';
   entries = nnz(A);
end

fclose(mmfile);
% Done.

Write data

function [ err ] = mmwrite(filename,A,comment,field,precision)
%
% Function: mmwrite(filename,A,comment,field,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.
%
%                 field     -  '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, will check after FIND-ing A
  precision = 16;
elseif ( nargin == 2) 
  comment = '';
  % Check whether there is an imaginary part:
  mattype = 'real'; % placeholder, will check after FIND-ing A
  precision = 16;
end

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


[M,N] = size(A);

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

  [I,J,V] = find(A);
  if ( sum(abs(imag(nonzeros(V)))) > 0 )
    Vreal = 0; 
  else 
    Vreal = 1; 
  end

  if ( ~ strcmp(mattype,'pattern') & Vreal )
    mattype = 'real'; 
  elseif ( ~ strcmp(mattype,'pattern') )
    mattype = 'complex';
  end
%
% Determine symmetry:
%
  if ( M ~= N )
    symm = 'general';
    issymm = 0;
    NZ = length(V);
  else
    issymm = 1;
    NZ = length(V);
    for i=1:NZ
      if ( A(J(i),I(i)) ~= V(i) )
        issymm = 0;
        break;
      end
    end
    if ( issymm )
      symm = 'symmetric';
      ATEMP = tril(A);
      [I,J,V] = find(ATEMP);
      NZ = nnz(ATEMP);
    else
      isskew = 1;
      for i=1:NZ
        if ( A(J(i),I(i)) ~= - V(i) )
          isskew = 0;
          break;
        end
      end
      if ( isskew )
        symm = 'skew-symmetric';
        ATEMP = tril(A);
        [I,J,V] = find(ATEMP);
        NZ = nnz(ATEMP);
      elseif ( strcmp(mattype,'complex') )
        isherm = 1;
        for i=1:NZ
          if ( A(J(i),I(i)) ~= conj(V(i)) )
            isherm = 0;
            break;
          end
        end
        if ( isherm )
          symm = 'hermitian';
          ATEMP = tril(A);
          [I,J,V] = find(ATEMP);
          NZ = nnz(ATEMP);
        else 
          symm = 'general';
          NZ = nnz(A);
        end
      else
        symm = 'general';
        NZ = nnz(A);
      end
    end
  end

% Sparse coordinate format:

  rep = 'coordinate';


  fprintf(mmfile,'%%%%MatrixMarket matrix %s %s %s\n',rep,mattype,symm);
  [MC,NC] = 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);
  cplxformat = sprintf('%%d %%d %% .%dg %% .%dg\n',precision,precision);
  realformat = sprintf('%%d %%d %% .%dg\n',precision);
  if ( strcmp(mattype,'real') )
     for i=1:NZ
        fprintf(mmfile,realformat,I(i),J(i),V(i));
     end;
  elseif ( strcmp(mattype,'complex') )
  for i=1:NZ
     fprintf(mmfile,cplxformat,I(i),J(i),real(V(i)),imag(V(i)));
  end;
  elseif ( strcmp(mattype,'pattern') )
     for i=1:NZ
        fprintf(mmfile,'%d %d\n',I(i),J(i));
     end;
  else  
     err = -1;
     disp('Unsupported mattype:')
     mattype
  end;

%%%%%%%%%%%%%       This part for dense matrices      %%%%%%%%%%%%%%%%
else
  if ( sum(abs(imag(nonzeros(A)))) > 0 )
    Areal = 0; 
  else 
    Areal = 1; 
  end
  if ( ~strcmp(mattype,'pattern') & Areal )
    mattype = 'real';
  elseif ( ~strcmp(mattype,'pattern')  )
    mattype = 'complex';
  end
%
% Determine symmetry:
%
  if ( M ~= N )
    issymm = 0;
    symm = 'general';
  else
    issymm = 1;
    for j=1:N 
      for i=j+1:N
        if (A(i,j) ~= A(j,i) )
          issymm = 0;   
          break; 
        end
      end
      if ( ~ issymm ) break; end
    
    end
    if ( issymm )
      symm = 'symmetric';
    else
      isskew = 1;
      for j=1:N 
        for i=j+1:N
          if (A(i,j) ~= - A(j,i) )
            isskew = 0;   
            break; 
          end
        end
        if ( ~ isskew ) break; end
      end
      if ( isskew )
        symm = 'skew-symmetric';
      elseif ( strcmp(mattype,'complex') )
        isherm = 1;
        for j=1:N 
          for i=j+1:N
            if (A(i,j) ~= conj(A(j,i)) )
              isherm = 0;   
              break; 
            end
          end
          if ( ~ isherm ) break; end
        end
        if ( isherm )
          symm = 'hermitian';
        else 
          symm = 'general';
        end
      else
        symm = 'general';
      end
    end
  end

% Dense array format:

  rep = 'array';
  [MC,NC] = 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);
  cplxformat = sprintf('%% .%dg %% .%dg\n', precision,precision);
  realformat = sprintf('%% .%dg\n', precision);
  if ( ~ strcmp(symm,'general') )
     rowloop = 'j';
  else 
     rowloop = '1';
  end
  if ( strcmp(mattype,'real') )
     for j=1:N
       for i=eval(rowloop):M
          fprintf(mmfile,realformat,A(i,j));
       end
     end
  elseif ( strcmp(mattype,'complex') )
     for j=1:N
       for i=eval(rowloop):M
          fprintf(mmfile,cplxformat,real(A(i,j)),imag(A(i,j)));
       end
     end
  elseif ( strcmp(mattype,'pattern') )
     err = -2
     disp('Pattern type inconsistant with dense matrix')
  else
     err = -2
     disp('Unknown matrix type:')
     mattype
  end
end

fclose(mmfile);

测试

% 读取 bcsstk14.mtx 的矩阵信息存储于矩阵 A 中， rows cols 分别是行列
[A,rows,cols,entries,rep,field,symm] = mmread("bcsstk14.mtx");
% 将矩阵 A 写到 data.txt 文件中
mmwrite("data.txt", A, field)
[A2,rows,cols,entries,rep,field,symm] = mmread("data.txt");
error = A - A2;
error_norm = norm(error, 2)

C++ IO

在 网址 C++ IO 下载 mmio.h 和 mmio.c 文件

Read and write data

#include <iostream>
#include <fstream>
#include <cstdio>
#include <cstdlib>
#include <string.h>
#include <stdlib.h>
#include <string>
#include <map>
#include "mmio.h"
#include <stdio.h>

using namespace std;
 
int main(int argc, char *argv[])
{
#if 1
	// read data
	string fileName = "../data/bcsstk14.mtx"; 
    int ret_code;
    MM_typecode matcode;
    FILE *f;
    int M, N, nz;   
    int i, *I, *J;
    double *val;

	map<unsigned int, map<unsigned int, double> > A; // Storage in triplet mode， <row <col val> > ; A's id from zero

    if ((f = fopen(fileName.c_str(), "r")) == NULL) {
      cerr<<"read file error, please check."<<endl;
      exit(1);
    }

    if (mm_read_banner(f, &matcode) != 0){
        printf("Could not process Matrix Market banner.\n");
        exit(1);
    }

    /*  This is how one can screen matrix types if their application */
    /*  only supports a subset of the Matrix Market data types.      */
    if (mm_is_complex(matcode) && mm_is_matrix(matcode) && 
            mm_is_sparse(matcode) ) {
        printf("Sorry, this application does not support ");
        printf("Market Market type: [%s]\n", mm_typecode_to_str(matcode));
        exit(1);
    }

    /* find out size of sparse matrix .... */
    if ((ret_code = mm_read_mtx_crd_size(f, &M, &N, &nz)) !=0) {
        exit(1);
		}

    /* reseve memory for matrices */
    I = (int *) malloc(nz * sizeof(int));
    J = (int *) malloc(nz * sizeof(int));
    val = (double *) malloc(nz * sizeof(double));

    /* NOTE: when reading in doubles, ANSI C requires the use of the "l"  */
    /*   specifier as in "%lg", "%lf", "%le", otherwise errors will occur */
    /*  (ANSI C X3.159-1989, Sec. 4.9.6.2, p. 136 lines 13-15)            */
    for (i=0; i<nz; i++) {
        fscanf(f, "%d %d %lg\n", &I[i], &J[i], &val[i]);
        I[i]--;  /* adjust from 1-based to 0-based */
        J[i]--;
    }

    if (f != stdin)  {
		fclose(f);
	}
    /************************/
    /* now write out matrix */
    /************************/
    mm_write_banner(stdout, matcode);
    mm_write_mtx_crd_size(stdout, M, N, nz);
    for (i=0; i<nz; i++) {
        //fprintf(stdout, "%d %d %20.19g\n", I[i]+1, J[i]+1, val[i]);
		A[I[i]][J[i]] = val[i];
		if(mm_is_symmetric(matcode)) {
          _A[J[i]][I[i]] = val[i];
        }
	}

#endif

#if 0
		// write data
		    
	const int nz = 4;
	const int M = 10;
	const int N = 10;

	MM_typecode matcode;                        
	int I[nz] = { 0, 4, 2, 8 };
    int J[nz] = { 3, 8, 7, 5 };
    double val[nz] = {1.1, 2.2, 3.2, 4.4};
    int i;

    mm_initialize_typecode(&matcode);
    mm_set_matrix(&matcode);
    mm_set_coordinate(&matcode);
    mm_set_real(&matcode);

    mm_write_banner(stdout, matcode); 
    mm_write_mtx_crd_size(stdout, M, N, nz);

    /* NOTE: matrix market files use 1-based indices, i.e. first element
      of a vector has index 1, not 0.  */

    for (i=0; i<nz; i++) {
        fprintf(stdout, "%d %d %10.3g\n", I[i]+1, J[i]+1, val[i]);
	}
#endif
		return 0;
}

Download Matrix

在网址的主页面可以看到 Search by matrix properties，

在之后的红色框中选择所需要的矩阵性质 回车即可。

最后利用前面介绍的 IO 的代码， 即可把矩阵导入到自己的线性求解器中， 作为验证。

IO 代码下载

以上的所有代码包括一个测试的例子均可以到我的 GitHub 网页下载，其中包含 MakeFile文件。

参考网址

Matrix Market ： https://math.nist.gov/MatrixMarket/