$title Test cholesky utility (CHOLES01,SEQ=411)

$onText

Finds the cholesky decomposition A=LL' of a positive definite symmetric matrix
A through an external program

Contributor: Erwin Kalvelagen, Amsterdam Optimization

$offText

set i  / i1*i5 /;
alias (i,j,n);

table a(i,j) 'original matrix'
      i1     i2     i3     i4     i5
i1    64     48     24      8      8
i2    48     72     42     54     36
i3    24     42     89    107     95
i4     8     54    107    210    186
i5     8     36     95    186    187
;

scalar rc;
parameter L(i,j) 'cholesky factor';

execute_unload 'a.gdx', i, a;
executeTool.checkErrorLevel 'linalg.cholesky i a L -gdxin=a.gdx -gdxout=b.gdx';
execute_load 'b.gdx', L;

* Check if Cholesky factorization is correct
Parameter a_, adiff;
a_(i,j)    = sum(n, L(i,n)*L(j,n));
adiff(i,j) = round(a(i,j) - a_(i,j),1e-6);
option adiff:8:0:1;
abort$card(adiff) adiff;

option clear=L;
executeTool.checkErrorLevel 'linalg.cholesky i a L';

* Check if Cholesky factorization is correct
a_(i,j)    = sum(n, L(i,n)*L(j,n));
adiff(i,j) = round(a(i,j) - a_(i,j),1e-6);
option adiff:8:0:1;
abort$card(adiff) adiff;

*
* only lower triangular part of A is used
*
table a2(i,j) 'original matrix'
      i1     i2     i3     i4     i5
i1    64
i2    48     72
i3    24     42     89
i4     8     54    107    210
i5     8     36     95    186    187
;

parameter L2(i,j) 'cholesky factor';
execute_unload 'a.gdx', i, a2;
executeTool.checkErrorLevel 'linalg.cholesky i a2 L2 -gdxin=a.gdx -gdxout=b.gdx';
execute_load 'b.gdx', L2;

* Check if Cholesky factorization is correct
Parameter a2_, a2diff;
a2_(i,j)    = sum(n, L2(i,n)*L2(j,n));
a2diff(i,j)$(ord(j)<=ord(i)) = round(a2(i,j) - a2_(i,j),1e-6);
option a2diff:8:0:1;
abort$card(a2diff) a2diff;

option clear=L2;
executeTool.checkErrorLevel 'linalg.cholesky i a2 L2';

* Check if Cholesky factorization is correct
a2_(i,j)    = sum(n, L2(i,n)*L2(j,n));
a2diff(i,j)$(ord(j)<=ord(i)) = round(a2(i,j) - a2_(i,j),1e-6);
option a2diff:8:0:1;
abort$card(a2diff) a2diff;

*
* only upper triangular part of A is used
*
Parameter a3(i,j); a3(i,j) = a2(j,i)

parameter L3(i,j) 'cholesky factor';
execute_unload 'a.gdx', i, a3;
executeTool.checkErrorLevel 'linalg.cholesky i a3 L3 -gdxin=a.gdx -gdxout=b.gdx';
execute_load 'b.gdx', L3;

* Check if Cholesky factorization is correct
Parameter a3_, a3diff;
a3_(i,j)    = sum(n, L3(n,i)*L3(n,j));
a3diff(i,j)$(ord(j)>=ord(i)) = round(a3(i,j) - a3_(i,j),1e-6);
option a3diff:8:0:1;
abort$card(a3diff) a3,a3_,a3diff;

option clear=L3;
executeTool.checkErrorLevel 'linalg.cholesky i a3 L3';

* Check if Cholesky factorization is correct
a3_(i,j)    = sum(n, L3(n,i)*L3(n,j));
a3diff(i,j)$(ord(j)>=ord(i)) = round(a3(i,j) - a3_(i,j),1e-6);
option a3diff:8:0:1;
abort$card(a3diff) a3diff;