clc;clear;
A = zeros(80,80);
ma1 = zeros(5,5);
ma1(1,1) = 5;
ma1(1,2) = 2;
ma1(5,5) = 5;
ma1(5,4) = 2;
ma2 = zeros(5,5);
for i = 1 : 5
ma2(i,i) = 1;
end
for i = 1 : 4
ma2(i,i+1) = 1;
end
for i = 2 : 4
ma1(i,i) = 5;
ma1(i,i-1) = 2;
ma1(i,i+1) = 2;
end
for i = 1 : 5 : 76
A(i:i+4,i:i+4) = ma1;
end
for i = 1: 5 : 36
A(i:i+4,i+40:i+44) = ma2;
end
for i = 41 : 5 : 76
A(i:i+4,i-40:i-36) = ma2;
end
b = zeros(80,1);
for i = 1 : 80
b(i) = i;
end
x1 = A\b;
norm(A*x1-b);
%x_2 = gmres(A,b);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%非0块的数量
num_block = 0;
rows = 80;
blockSize = 5;
row_num = rows / blockSize;
%遍历块矩阵
for i = 1 : 5 : 76
for j = 1 : 5 :76
flag = notzero(A(i:i+4,j:j+4),5);
if(flag == true)
num_block = num_block + 1;
end
end
end
colIndices = zeros(num_block,1);
row_offset = zeros(row_num+1,1);
row_offset(1) = 0;
row_offset(row_num+1) = num_block;
m_value = zeros(num_block * blockSize * blockSize,1);
k = 1;
num_block = 0;
%遍历块矩阵
for i = 1 : 5 : 76
for j = 1 : 5 :76
flag = notzero(A(i:i+4,j:j+4),5);
if(flag == true)
num_block = num_block + 1;
colIndices(num_block) = fix(j/5);
temp_block = A(i:i+4,j:j+4);
for p = 1:blockSize
for q = 1 : blockSize
m_value(k) = temp_block(p,q);
k = k+1;
end
end
end
end
row_offset(fix(i/5)+2) = num_block;
end
%%%%%%%%%%%%%%%%%%%%输出
fid = fopen('block_mat_csr.txt','wt');
fprintf(fid,'%.0f',row_num);%矩阵大小
fprintf(fid,'\n');
fprintf(fid,'%.0f',num_block);%非0矩阵个数
fprintf(fid,'\n');
for i = 1 : num_block
fprintf(fid,'%.0f',colIndices(i));
fprintf(fid,'\n');
end
for i = 1 : row_num+1
fprintf(fid,'%.0f',row_offset(i));
fprintf(fid,'\n');
end
fclose('all');
fid = fopen('block_mat_value.txt','wt');
for i = 1 : num_block*blockSize*blockSize
fprintf(fid,'%.6f',m_value(i));
fprintf(fid,'\n');
end
[matlab]生成块稀疏矩阵并转为csr输出
最新推荐文章于 2023-10-04 10:48:53 发布