稀疏矩阵绘制

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:

下载网址：http://people.sc.fsu.edu/~burkardt/m_src/hb_to_mm/hb_to_mm.m

 

 下面附上绘制的结果图：

参考文献：

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

http://hi.baidu.com/guole205/blog/item/6be737fd51791f40d6887dce.html

注：

【问题】昨天这个程序似乎有些问题：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) )

即可。

【结果】

转自：http://blog.sina.com.cn/s/blog_4b25f6ad0100je58.html