Mercurial > matrix-functions
comparison mftoolbox/cosmsinm.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 [C,S,cost,s] = cosmsinm(A) | |
| 2 %COSMSINM Matrix cosine and sine by double angle algorithm. | |
| 3 % [C,S,COST,s] = COSMSINM(A) computes the cosine C and the sine S | |
| 4 % of the matrix A using the double angle algorithm with Pade | |
| 5 % approximation. The total number of matrix multiplications and | |
| 6 % linear systems solved is returned as COST and s specifies the | |
| 7 % amount of scaling (A is scaled to 2^(-s)*A). | |
| 8 | |
| 9 n = length(A); | |
| 10 theta = [0.00613443965526 0.11110098037055 0.43324697422211 0.98367255274344 ... | |
| 11 1.72463663220280 2.61357494421368 3.61521023400301 4.70271938553349 ... | |
| 12 5.85623410320942 7.06109053959248 0 9.58399601173102]; | |
| 13 | |
| 14 normA = norm(A,inf); | |
| 15 d = [2 4 6 8 10 12 14 16 20]; | |
| 16 for i = d(1:8) | |
| 17 if normA <= theta(i/2); | |
| 18 m = i; | |
| 19 cost = m/2+3; | |
| 20 [C,S] = cosmsinm_pade(A,m); | |
| 21 s = m; | |
| 22 return | |
| 23 end | |
| 24 end | |
| 25 | |
| 26 s = 0; | |
| 27 if normA > theta(20/2) | |
| 28 s = max( ceil( log2( normA/theta(20/2) ) ), 1 ); | |
| 29 A = A/(2^s); | |
| 30 normA = normA/(2^s); | |
| 31 end | |
| 32 | |
| 33 if normA > 2*theta(12/2) | |
| 34 m = 20; cost = 12; | |
| 35 elseif normA > theta(16/2) | |
| 36 A = A/2; m = 12; cost = 9; s = s+1; | |
| 37 else | |
| 38 m = 16; cost = 11; | |
| 39 end | |
| 40 | |
| 41 [C,S] = cosmsinm_pade(A,m); | |
| 42 for i = 1:s | |
| 43 S = 2*S*C; | |
| 44 C = 2*(C^2) - eye(n); | |
| 45 end | |
| 46 cost = cost + 2*s; |