Matlab中Harwell-boeing格式稀疏矩阵的读取及绘图_模栖之鹰_新浪博客

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

 

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;                                                          %得到{zh1}的刚度矩阵

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;                %总行数%得到{zh1}的刚度矩阵

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) )

即可。

【结果】



郑重声明:资讯 【Matlab中Harwell-boeing格式稀疏矩阵的读取及绘图_模栖之鹰_新浪博客】由 发布,版权归原作者及其所在单位,其原创性以及文中陈述文字和内容未经(企业库qiyeku.com)证实,请读者仅作参考,并请自行核实相关内容。若本文有侵犯到您的版权, 请你提供相关证明及申请并与我们联系(qiyeku # qq.com)或【在线投诉】,我们审核后将会尽快处理。
—— 相关资讯 ——