Description

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

Contributor: Erwin Kalvelagen, Amsterdam Optimization

Small Model of Type : GAMS

Category : GAMS Test library

Main file : choles01.gms

$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;