Mercurial > matrix-functions
comparison mftoolbox/polyvalm_ps.m @ 0:8f23314345f4 draft
Create local repository for matrix toolboxes. Step #0 done.
| author | Antonio Pino Robles <data.script93@gmail.com> |
|---|---|
| date | Wed, 06 May 2015 14:56:53 +0200 |
| parents | |
| children |
comparison
equal
deleted
inserted
replaced
| -1:000000000000 | 0:8f23314345f4 |
|---|---|
| 1 function [P,s,cost] = polyvalm_ps(c,A,s) | |
| 2 %POLYVALM_PS Evaluate polynomial at matrix argument by Paterson-Stockmeyer alg. | |
| 3 % [P,S,COST] = POLYVALM_PS(C,A,S) evaluates the polynomial whose | |
| 4 % coefficients are the vector C at the matrix A using the | |
| 5 % Paterson-Stockmeyer algorithm. If omitted, the integer parameter | |
| 6 % S is chosen automatically and its value is returned as an | |
| 7 % output argument. COST is the number of matrix multiplications used. | |
| 8 | |
| 9 m = length(c)-1; % Degree of poly. | |
| 10 c = c(end:-1:1); c = c(:); | |
| 11 n = length(A); | |
| 12 | |
| 13 if nargin < 3 | |
| 14 % Determine optimum parameter s. | |
| 15 s = ceil(sqrt(m)); | |
| 16 end | |
| 17 r = floor(m/s); | |
| 18 cost = s+r-(m==r*s)-1; | |
| 19 | |
| 20 % Apower{i+1} = A^i; | |
| 21 Apower = cell(s+1); | |
| 22 Apower{1} = eye(n); | |
| 23 for i=2:s+1 | |
| 24 Apower{i} = A*Apower{i-1}; | |
| 25 end | |
| 26 | |
| 27 B = cell(r+1); | |
| 28 for k=0:r-1 | |
| 29 temp = c(s*k+1)*eye(n); | |
| 30 for j=1:s-1 | |
| 31 temp = temp + c(s*k+j+1)*Apower{j+1}; | |
| 32 end | |
| 33 B{k+1} = temp; | |
| 34 end | |
| 35 B{r+1} = c(m+1)*Apower{m-s*r+1}; | |
| 36 for j=m-1:-1:s*r | |
| 37 if j == s*r | |
| 38 B{r+1} = B{r+1} + c(s*r+1)*eye(n); | |
| 39 else | |
| 40 B{r+1} = B{r+1} + c(j+1)*Apower{m-s*r-(m-j)+1}; | |
| 41 end | |
| 42 end | |
| 43 | |
| 44 As = Apower{s+1}; | |
| 45 P = zeros(n); | |
| 46 for k=r:-1:0 | |
| 47 P = P*As + B{k+1}; | |
| 48 end |