summaryrefslogtreecommitdiff
path: root/inst/misc/benchmarks/stk_benchmark_uppertriangular.m
blob: a9092a8ab29de8bfe6ac59352d2afc1dc2262053 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
% STK_BENCHMARK_UPPERTRIANGULAR times some computations with UT matrices

% Copyright Notice
%
%    Copyright (C) 2018 CentraleSupelec
%
%    Author:  Julien Bect  <julien.bect@centralesupelec.fr>

% Copying Permission Statement
%
%    This file is part of
%
%            STK: a Small (Matlab/Octave) Toolbox for Kriging
%               (https://github.com/stk-kriging/stk/)
%
%    STK is free software: you can redistribute it and/or modify it under
%    the terms of the GNU General Public License as published by the Free
%    Software Foundation,  either version 3  of the License, or  (at your
%    option) any later version.
%
%    STK is distributed  in the hope that it will  be useful, but WITHOUT
%    ANY WARRANTY;  without even the implied  warranty of MERCHANTABILITY
%    or FITNESS  FOR A  PARTICULAR PURPOSE.  See  the GNU  General Public
%    License for more details.
%
%    You should  have received a copy  of the GNU  General Public License
%    along with STK.  If not, see <http://www.gnu.org/licenses/>.


%% Prepare benchmark

n = 100;  NREP = 200;

model = stk_model (@stk_materncov32_iso);
model.param = [0 0];
x = (linspace (0, 1, n))';
K = stk_make_matcov (model, x);

% Assume that the Cholesky factorization of K is available
C = chol (K);


%% First test: K ^ (-1)

tic ();
for i = 1:NREP
    M1 = inv (K);
end
t1 = toc ();

tic
for i = 1:NREP
    B2 = inv (C');
    M2 = B2' * B2;
end
t2 = toc ();

tic
for i = 1:NREP
    B3 = linsolve (C, eye (n), struct ('UT', true, 'TRANSA', true));
    M3 = B3' * B3;
end
t3 = toc ();

t = [t1 t2 t3]'


%% Second test: W * K^(-1) * W'

W = randn (n);

tic ();
for i = 1:NREP
    M1 = W * (inv (K)) * (W');
end
t1 = toc ();

tic
for i = 1:NREP
    B2 = (C') \ (W');
    M2 = B2' * B2;
end
t2 = toc ();

tic
for i = 1:NREP
    B3 = linsolve (C, W', struct ('UT', true, 'TRANSA', true));
    M3 = B3' * B3;
end
t3 = toc ();

t = [t1 t2 t3]'