声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 2393|回复: 1

[其他相关] 请教有关“数值代数中Arnoldi 算法的块状化的比较”的问题

[复制链接]
发表于 2012-12-19 10:37 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
数值代数中Arnoldi 算法的块状化的比较;
%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%使用 double 数据类型%%
function [v,H] = blockArnoldi(A,m,r)
% % Input: A -- an n by n matrix
% m -- a positive integer
% r -- an Nxp block vector (v .ne. 0 assumed)
% % Output: V -- an n by m*p orthogonal matrix
% H -- a m*p by (m-1)*p upper Hessenberg matrix
if isnumeric(A)
  A = @(x)(A*x);
end
n = length(r);
p = size(r,2);
H (m*p,(m-1)*p)= 0;
v (n, m*p) = 0;
[v(:,1:p),H(1:p,1:p)] = sqr(r);
k = 0;
while k < m
     k = k +1;
     kp = k*p;
     tic;
     w = A(v(:,kp-p+1:kp));  
     toc,
for j =1:k
     jp = j*p;
     h = v(:,jp-p+1:jp)' * w;  
     w = w - v(:,jp-p+1:jp)*h;
     H(jp-p+1:jp,kp-p+1:kp) = h;   
end
    [v(:,kp+1:kp+p) ,H(kp+1:kp+p,kp-p+1:kp)] = sqr(w);  
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 使用 cell 数据类型

function [v,H] = blockArnoldi_cell(A,m,r)
% % Input: A -- an n by n matrix
% m -- a positive integer
% r -- an Nxp block vector (v .ne. 0 assumed)
% % Output: V -- an n by m*p orthogonal matrix
% H -- a m*p by (m-1)*p upper Hessenberg matrix
if isnumeric(A)
  A = @(x)(A*x);
end
%n = length(r);
p = size(r,2);
%H (m*p,(m-1)*p)= 0;
%v (n, m*p) = 0;
v = {}; h ={};
[v{1},h{1,1}] = sqr(r);
k = 0;
while k < m
     k = k + 1;
     w = A(v{k});   
for j =1:k
     h{j,k} = v{j}' * w;  
     w = w - v{j} * h{j,k};  
end
    [v{k+1},h{k+1,k}] = sqr(w);
end
H = zeros(p,p);
for j = 1: m-1
    for  k = j+2 : m+1
       h{k,j} = H;
    end
end
%  h
cell2mat(v);
H =cell2mat(h);
%%%%%%%%%%%%%%%%%%%%
实验结果: cell 数据类型的程序比double类型的程序快!
问题1,不是double类型数组访问是效率比cell数据类型的数组访问效率快很多吗? 为何这里的cell 数据类型的程序比double类型的程序快?
问题2, 如何把下列程序向量化?或者有更好的方法使得下列程序执行效率快一些?
for j = 1: m-1
    for  k = j+2 : m+1
       h{k,j} = H;
    end
end
%%%%%%%%%%%
在这里先谢谢各位了!!

回复
分享到:

使用道具 举报

 楼主| 发表于 2012-12-20 16:50 | 显示全部楼层
怎么没有人回帖啊.....................
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-17 09:53 , Processed in 0.087923 second(s), 18 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表