最近在做稀疏矩阵的相关工作，其中一项是绘制稀疏矩阵的结构图。经过一下午研究，将相关成果整理如下：

matlab源代码为：

1、文件 test_bcsstk01.m:

% 绘制稀疏矩阵 bcsstk01.rsa 和 bcsstm01.rsa 的结构图

clear

% 读入rsa形式的刚度矩阵

Dpi = 150;

hb_to_mmm('matrix\bcsstk01.rsa','matrix\bcsstk01');

load matrix\bcsstk01;

ni=bcsstk01(1,1); %稀疏矩阵的行

nj=bcsstk01(1,2); %稀疏矩阵的列

nz=bcsstk01(1,3); %稀疏矩阵的非零元总数

nl=nz+1; %总行数

KL=sparse( bcsstk01(2:nl,1), bcsstk01(2:nl,2), bcsstk01(2:nl,3),

ni, nj ); %转换为稀疏矩阵(刚度矩阵下三角)

[dr,va]=get_diag(bcsstk01(2:nl,1),bcsstk01(2:nl,2),bcsstk01(2:nl,3),nz); %提取对角线元素

KD=sparse( dr, dr, va, ni, nj

); %刚度矩阵的对角线元素

K = KL + KL' -

KD; %得到最后的刚度矩阵

clear KL;

clear KD;

clear dr;

clear va;

figure(1)

subplot(1, 2, 1)

spy(K)

title('矩阵 bcsstk01.rsa 的结构图');

% 读入rsa形式的质量矩阵

hb_to_mmm('matrix\bcsstm01.rsa','matrix\bcsstm01');

load matrix\bcsstm01;

ni=bcsstm01(1,1); %稀疏矩阵的行

nj=bcsstm01(1,2); %稀疏矩阵的列

nz=bcsstm01(1,3); %稀疏矩阵的非零元总数

nl=nz+1; %总行数%得到最后的刚度矩阵

ML=sparse( bcsstm01(2:nl,1), bcsstm01(2:nl,2), bcsstm01(2:nl,3),

ni, nj ); %转换为稀疏矩阵

[dr,va]=get_diag(bcsstm01(2:nl,1),bcsstm01(2:nl,2),bcsstm01(2:nl,3),nz); %提取对角线元素

MD=sparse( dr, dr, va, ni, nj

); %刚度矩阵的对角线元素

M = ML + ML' -

MD;

clear ML;

clear MD;

clear dr;

clear va;

subplot(1, 2, 2)

spy(M)

title('矩阵 bcsstm01.rsa 的结构图');

SDpi=['-r',num2str(Dpi)]; % 设定输出图像的 Dpi 值

print(1,'-dtiff', SDpi, 'bcsstk01_bcsstm01'

); % 将图像按规定 Dpi 和 文件名写入到文件

2、文件 get_diag.m:

function [dr,va]=get_diag(row,col,val,nz)

index=1;

sum =0;

for i=1:nz

if( row(i,1)

== col(i,1) )

dr(index,1) = row(i,1);

va(index,1) = val(i,1);

if( val(i,1) < 1e-20 )

sum = sum + 1;

end

index = index + 1;

end

end

if( sum > 0 )

warning('对角线元素上有 %d 个小于或等于0的元素.',sum)

end

end

3、文件 hb_to_mmm.m:

下面附上绘制的结果图：

参考文献：

【总结】Harwell-Boeing格式矩阵转换为稀疏矩阵的方法

注：

【问题】昨天这个程序似乎有些问题：hb_to_mmm.m文件得到的非零元个数似乎比实际的要少。

例如：bcsstk09.rsa这个文件，如果只考虑下三角部分，实际上共有9760个非零元，而若先用hb_to_mmm.m文件得到一个matlab格式的矩阵文件bcsstk09，再读入用spconvert函数转化成稀疏矩阵之后，虽然仍有9760个非零元，但其中许多非零元的数值却被设置为0了(直接查看bcsstk09.rsa文件，发现这些位置处的非零元数值为e-7/e-8量级)！这样，当再调用spy函数画出刚度矩阵下三角的非零元个数时，就只有9039个了。

【解决方案】修改hb_to_mmm.m文件的msm_to_mm_coordinate_real_general函数，将第1445行：

fprintf ( fid, ' %d %d %f\n', rows(k2), j, vals(k2) )

改为

fprintf ( fid, ' %d %d %e\n', rows(k2),

j, vals(k2) )

即可。

【结果】