Mercurial > octave-nkf

comparison src/DLD-FUNCTIONS/chol.cc @ 13953:642e43164af6

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
fix behavior of chol (..., 'lower') to be compatible with matlab * chol.cc: transpose input matrix prior to factorization when chol (..., 'lower') is invoked so that only the lower triangular part is used.
author Carlo de Falco <kingcrimson@tiscali.it>
date Mon, 28 Nov 2011 12:39:39 +0100
parents 5fa482628bf6
children b0cdd60db5e5
comparison
equal deleted inserted replaced
13952:acaf33ccc04f 13953:642e43164af6
59 return octave_value (fact.chol_matrix ().transpose (), 59 return octave_value (fact.chol_matrix ().transpose (),
60 MatrixType (MatrixType::Lower)); 60 MatrixType (MatrixType::Lower));
61 } 61 }
62 62
63 DEFUN_DLD (chol, args, nargout, 63 DEFUN_DLD (chol, args, nargout,
64 "-*- texinfo -*-\n\ 64 "-*- texinfo -*-\n\
65 @deftypefn {Loadable Function} {@var{R} =} chol (@var{A})\n\ 65 @deftypefn {Loadable Function} {@var{R} =} chol (@var{A})\n\
66 @deftypefnx {Loadable Function} {[@var{R}, @var{p}] =} chol (@var{A})\n\ 66 @deftypefnx {Loadable Function} {[@var{R}, @var{p}] =} chol (@var{A})\n\
67 @deftypefnx {Loadable Function} {[@var{R}, @var{p}, @var{Q}] =} chol (@var{S})\n\ 67 @deftypefnx {Loadable Function} {[@var{R}, @var{p}, @var{Q}] =} chol (@var{S})\n\
68 @deftypefnx {Loadable Function} {[@var{R}, @var{p}, @var{Q}] =} chol (@var{S}, 'vector')\n\ 68 @deftypefnx {Loadable Function} {[@var{R}, @var{p}, @var{Q}] =} chol (@var{S}, 'vector')\n\
69 @deftypefnx {Loadable Function} {[@var{L}, @dots{}] =} chol (@dots{}, 'lower')\n\ 69 @deftypefnx {Loadable Function} {[@var{L}, @dots{}] =} chol (@dots{}, 'lower')\n\
70 @deftypefnx {Loadable Function} {[@var{L}, @dots{}] =} chol (@dots{}, 'upper')\n\
70 @cindex Cholesky factorization\n\ 71 @cindex Cholesky factorization\n\
71 Compute the Cholesky@tie{}factor, @var{R}, of the symmetric positive definite\n\ 72 Compute the Cholesky@tie{}factor, @var{R}, of the symmetric positive definite\n\
72 matrix @var{A}, where\n\ 73 matrix @var{A}, where\n\
73 @tex\n\ 74 @tex\n\
74 $ R^T R = A $.\n\ 75 $ R^T R = A $.\n\
125 @example\n\ 126 @example\n\
126 @var{L} * @var{L}' = @var{A}.\n\ 127 @var{L} * @var{L}' = @var{A}.\n\
127 @end example\n\ 128 @end example\n\
128 \n\ 129 \n\
129 @end ifnottex\n\ 130 @end ifnottex\n\
131 \n\
132 For full matrices, if the 'lower' flag is set only the lower triangular part of the matrix \
133 is used for the factorization, otherwise the upper triangular part is used.\n\
130 \n\ 134 \n\
131 In general the lower triangular factorization is significantly faster for\n\ 135 In general the lower triangular factorization is significantly faster for\n\
132 sparse matrices.\n\ 136 sparse matrices.\n\
133 @seealso{cholinv, chol2inv}\n\ 137 @seealso{cholinv, chol2inv}\n\
134 @end deftypefn") 138 @end deftypefn")
153 if (! error_state ) 157 if (! error_state )
154 { 158 {
155 if (tmp.compare ("vector") == 0) 159 if (tmp.compare ("vector") == 0)
156 vecout = true; 160 vecout = true;
157 else if (tmp.compare ("lower") == 0) 161 else if (tmp.compare ("lower") == 0)
162 // FIXME currently the option "lower" is handled by transposing the
163 // matrix, factorizing it with the lapack function DPOTRF ('U', ...)
164 // and finally transposing the factor. It would be more efficient to use
165 // DPOTRF ('L', ...) in this case.
158 LLt = true; 166 LLt = true;
159 else if (tmp.compare ("upper") == 0) 167 else if (tmp.compare ("upper") == 0)
160 LLt = false; 168 LLt = false;
161 else 169 else
162 error ("chol: unexpected second or third input"); 170 error ("chol: unexpected second or third input");
249 FloatMatrix m = arg.float_matrix_value (); 257 FloatMatrix m = arg.float_matrix_value ();
250 258
251 if (! error_state) 259 if (! error_state)
252 { 260 {
253 octave_idx_type info; 261 octave_idx_type info;
254 FloatCHOL fact (m, info); 262
263 FloatCHOL fact;
264 if (LLt)
265 fact = FloatCHOL (m.transpose (), info);
266 else
267 fact = FloatCHOL (m, info);
268
255 if (nargout == 2 || info == 0) 269 if (nargout == 2 || info == 0)
256 { 270 {
257 retval(1) = info; 271 retval(1) = info;
258 if (LLt) 272 if (LLt)
259 retval(0) = get_chol_l (fact); 273 retval(0) = get_chol_l (fact);
269 FloatComplexMatrix m = arg.float_complex_matrix_value (); 283 FloatComplexMatrix m = arg.float_complex_matrix_value ();
270 284
271 if (! error_state) 285 if (! error_state)
272 { 286 {
273 octave_idx_type info; 287 octave_idx_type info;
274 FloatComplexCHOL fact (m, info); 288
289 FloatComplexCHOL fact;
290 if (LLt)
291 fact = FloatComplexCHOL (m.transpose (), info);
292 else
293 fact = FloatComplexCHOL (m, info);
294
275 if (nargout == 2 || info == 0) 295 if (nargout == 2 || info == 0)
276 { 296 {
277 retval(1) = info; 297 retval(1) = info;
278 if (LLt) 298 if (LLt)
279 retval(0) = get_chol_l (fact); 299 retval(0) = get_chol_l (fact);
294 Matrix m = arg.matrix_value (); 314 Matrix m = arg.matrix_value ();
295 315
296 if (! error_state) 316 if (! error_state)
297 { 317 {
298 octave_idx_type info; 318 octave_idx_type info;
299 CHOL fact (m, info); 319
320 CHOL fact;
321 if (LLt)
322 fact = CHOL (m.transpose (), info);
323 else
324 fact = CHOL (m, info);
325
300 if (nargout == 2 || info == 0) 326 if (nargout == 2 || info == 0)
301 { 327 {
302 retval(1) = info; 328 retval(1) = info;
303 if (LLt) 329 if (LLt)
304 retval(0) = get_chol_l (fact); 330 retval(0) = get_chol_l (fact);
314 ComplexMatrix m = arg.complex_matrix_value (); 340 ComplexMatrix m = arg.complex_matrix_value ();
315 341
316 if (! error_state) 342 if (! error_state)
317 { 343 {
318 octave_idx_type info; 344 octave_idx_type info;
319 ComplexCHOL fact (m, info); 345
346 ComplexCHOL fact;
347 if (LLt)
348 fact = ComplexCHOL (m.transpose (), info);
349 else
350 fact = ComplexCHOL (m, info);
351
320 if (nargout == 2 || info == 0) 352 if (nargout == 2 || info == 0)
321 { 353 {
322 retval(1) = info; 354 retval(1) = info;
323 if (LLt) 355 if (LLt)
324 retval(0) = get_chol_l (fact); 356 retval(0) = get_chol_l (fact);
991 %! R1 = cholinsert(R,j,u2); A1 = R1'*R1; 1023 %! R1 = cholinsert(R,j,u2); A1 = R1'*R1;
992 %! 1024 %!
993 %! assert(norm(triu(R1)-R1,Inf) == 0) 1025 %! assert(norm(triu(R1)-R1,Inf) == 0)
994 %! assert(norm(A1(p,p) - single(Ac),Inf) < 2e1*eps('single')) 1026 %! assert(norm(A1(p,p) - single(Ac),Inf) < 2e1*eps('single'))
995 %! 1027 %!
1028
1029 %!test
1030 %! cu = chol (triu (A), 'upper');
1031 %! cl = chol (tril (A), 'lower');
1032 %! assert (cu, cl', eps)
1033 %!
1034 %!test
1035 %! cca = chol (Ac);
1036 %!
1037 %! ccal = chol (Ac, 'lower');
1038 %! ccal2 = chol (tril (Ac), 'lower');
1039 %!
1040 %! ccau = chol (Ac, 'upper');
1041 %! ccau2 = chol (triu (Ac), 'upper');
1042 %!
1043 %! assert (cca'*cca, Ac, eps)
1044 %! assert (ccau'*ccau, Ac, eps)
1045 %! assert (ccau2'*ccau2, Ac, eps)
1046 %!
1047 %! assert (cca, ccal', eps)
1048 %! assert (cca, ccau, eps)
1049 %! assert (cca, ccal2', eps)
1050 %! assert (cca, ccau2, eps)
1051 %!
1052 %!test
1053 %! cca = chol (single (Ac));
1054 %!
1055 %! ccal = chol (single (Ac), 'lower');
1056 %! ccal2 = chol (tril (single (Ac)), 'lower');
1057 %!
1058 %! ccau = chol (single (Ac), 'upper');
1059 %! ccau2 = chol (triu (single (Ac)), 'upper');
1060 %!
1061 %! assert (cca'*cca, single (Ac), eps ('single'))
1062 %! assert (ccau'*ccau, single (Ac), eps ('single'))
1063 %! assert (ccau2'*ccau2, single (Ac), eps ('single'))
1064 %!
1065 %! assert (cca, ccal', eps ('single'))
1066 %! assert (cca, ccau, eps ('single'))
1067 %! assert (cca, ccal2', eps ('single'))
1068 %! assert (cca, ccau2, eps ('single'))
1069
1070 %!test
1071 %! a = [12, 2, 3, 4;
1072 %! 2, 14, 5, 3;
1073 %! 3, 5, 16, 6;
1074 %! 4, 3, 6, 16];
1075 %!
1076 %! b = [0, 1, 2, 3;
1077 %! -1, 0, 1, 2;
1078 %! -2, -1, 0, 1;
1079 %! -3, -2, -1, 0];
1080 %!
1081 %! ca = a + i*b;
1082 %!
1083 %! cca = chol (ca);
1084 %!
1085 %! ccal = chol (ca, 'lower');
1086 %! ccal2 = chol (tril (ca), 'lower');
1087 %!
1088 %! ccau = chol (ca, 'upper');
1089 %! ccau2 = chol (triu (ca), 'upper');
1090 %!
1091 %! assert (cca'*cca, ca, 16*eps)
1092 %! assert (ccau'*ccau, ca, 16*eps)
1093 %! assert (ccau2'*ccau2, ca, 16*eps)
1094 %!
1095 %! assert (cca, ccal', 16*eps)
1096 %! assert (cca, ccau, 16*eps)
1097 %! assert (cca, ccal2', 16*eps)
1098 %! assert (cca, ccau2, 16*eps)
1099
996 */ 1100 */
997 1101
998 DEFUN_DLD (choldelete, args, , 1102 DEFUN_DLD (choldelete, args, ,
999 "-*- texinfo -*-\n\ 1103 "-*- texinfo -*-\n\
1000 @deftypefn {Loadable Function} {@var{R1} =} choldelete (@var{R}, @var{j})\n\ 1104 @deftypefn {Loadable Function} {@var{R1} =} choldelete (@var{R}, @var{j})\n\