chol.cc

Go to the documentation of this file.
00001 /*
00002 
00003 Copyright (C) 1996-2012 John W. Eaton
00004 Copyright (C) 2008-2009 Jaroslav Hajek
00005 Copyright (C) 2008-2009 VZLU Prague
00006 
00007 This file is part of Octave.
00008 
00009 Octave is free software; you can redistribute it and/or modify it
00010 under the terms of the GNU General Public License as published by the
00011 Free Software Foundation; either version 3 of the License, or (at your
00012 option) any later version.
00013 
00014 Octave is distributed in the hope that it will be useful, but WITHOUT
00015 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00016 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00017 for more details.
00018 
00019 You should have received a copy of the GNU General Public License
00020 along with Octave; see the file COPYING.  If not, see
00021 <http://www.gnu.org/licenses/>.
00022 
00023 */
00024 
00025 
00026 #ifdef HAVE_CONFIG_H
00027 #include <config.h>
00028 #endif
00029 
00030 #include "CmplxCHOL.h"
00031 #include "dbleCHOL.h"
00032 #include "fCmplxCHOL.h"
00033 #include "floatCHOL.h"
00034 #include "SparseCmplxCHOL.h"
00035 #include "SparsedbleCHOL.h"
00036 #include "oct-spparms.h"
00037 #include "sparse-util.h"
00038 
00039 #include "ov-re-sparse.h"
00040 #include "ov-cx-sparse.h"
00041 #include "defun-dld.h"
00042 #include "error.h"
00043 #include "gripes.h"
00044 #include "oct-obj.h"
00045 #include "utils.h"
00046 
00047 template <class CHOLT>
00048 static octave_value
00049 get_chol_r (const CHOLT& fact)
00050 {
00051   return octave_value (fact.chol_matrix (),
00052                        MatrixType (MatrixType::Upper));
00053 }
00054 
00055 template <class CHOLT>
00056 static octave_value
00057 get_chol_l (const CHOLT& fact)
00058 {
00059   return octave_value (fact.chol_matrix ().transpose (),
00060                        MatrixType (MatrixType::Lower));
00061 }
00062 
00063 DEFUN_DLD (chol, args, nargout,
00064 "-*- texinfo -*-\n\
00065 @deftypefn  {Loadable Function} {@var{R} =} chol (@var{A})\n\
00066 @deftypefnx {Loadable Function} {[@var{R}, @var{p}] =} chol (@var{A})\n\
00067 @deftypefnx {Loadable Function} {[@var{R}, @var{p}, @var{Q}] =} chol (@var{S})\n\
00068 @deftypefnx {Loadable Function} {[@var{R}, @var{p}, @var{Q}] =} chol (@var{S}, 'vector')\n\
00069 @deftypefnx {Loadable Function} {[@var{L}, @dots{}] =} chol (@dots{}, 'lower')\n\
00070 @deftypefnx {Loadable Function} {[@var{L}, @dots{}] =} chol (@dots{}, 'upper')\n\
00071 @cindex Cholesky factorization\n\
00072 Compute the Cholesky@tie{}factor, @var{R}, of the symmetric positive definite\n\
00073 matrix @var{A}, where\n\
00074 @tex\n\
00075 $ R^T R = A $.\n\
00076 @end tex\n\
00077 @ifnottex\n\
00078 \n\
00079 @example\n\
00080 @var{R}' * @var{R} = @var{A}.\n\
00081 @end example\n\
00082 \n\
00083 @end ifnottex\n\
00084 \n\
00085 Called with one output argument @code{chol} fails if @var{A} or @var{S} is\n\
00086 not positive definite.  With two or more output arguments @var{p} flags\n\
00087 whether the matrix was positive definite and @code{chol} does not fail.  A\n\
00088 zero value indicated that the matrix was positive definite and the @var{R}\n\
00089 gives the factorization, and @var{p} will have a positive value otherwise.\n\
00090 \n\
00091 If called with 3 outputs then a sparsity preserving row/column permutation\n\
00092 is applied to @var{A} prior to the factorization.  That is @var{R}\n\
00093 is the factorization of @code{@var{A}(@var{Q},@var{Q})} such that\n\
00094 @tex\n\
00095 $ R^T R = Q^T A Q$.\n\
00096 @end tex\n\
00097 @ifnottex\n\
00098 \n\
00099 @example\n\
00100 @var{R}' * @var{R} = @var{Q}' * @var{A} * @var{Q}.\n\
00101 @end example\n\
00102 \n\
00103 @end ifnottex\n\
00104 \n\
00105 The sparsity preserving permutation is generally returned as a matrix.\n\
00106 However, given the flag 'vector', @var{Q} will be returned as a vector\n\
00107 such that\n\
00108 @tex\n\
00109 $ R^T R = A (Q, Q)$.\n\
00110 @end tex\n\
00111 @ifnottex\n\
00112 \n\
00113 @example\n\
00114 @var{R}' * @var{R} = @var{A}(@var{Q}, @var{Q}).\n\
00115 @end example\n\
00116 \n\
00117 @end ifnottex\n\
00118 \n\
00119 Called with either a sparse or full matrix and using the 'lower' flag,\n\
00120 @code{chol} returns the lower triangular factorization such that\n\
00121 @tex\n\
00122 $ L L^T = A $.\n\
00123 @end tex\n\
00124 @ifnottex\n\
00125 \n\
00126 @example\n\
00127 @var{L} * @var{L}' = @var{A}.\n\
00128 @end example\n\
00129 \n\
00130 @end ifnottex\n\
00131 \n\
00132 For full matrices, if the 'lower' flag is set only the lower triangular part\n\
00133 of the matrix is used for the factorization, otherwise the upper triangular\n\
00134 part is used.\n\
00135 \n\
00136 In general the lower triangular factorization is significantly faster for\n\
00137 sparse matrices.\n\
00138 @seealso{cholinv, chol2inv}\n\
00139 @end deftypefn")
00140 {
00141   octave_value_list retval;
00142   int nargin = args.length ();
00143   bool LLt = false;
00144   bool vecout = false;
00145 
00146   if (nargin < 1 || nargin > 3 || nargout > 3
00147       || (! args(0).is_sparse_type () && nargout > 2))
00148     {
00149       print_usage ();
00150       return retval;
00151     }
00152 
00153   int n = 1;
00154   while (n < nargin && ! error_state)
00155     {
00156       std::string tmp = args(n++).string_value ();
00157 
00158       if (! error_state )
00159         {
00160           if (tmp.compare ("vector") == 0)
00161             vecout = true;
00162           else if (tmp.compare ("lower") == 0)
00163             // FIXME currently the option "lower" is handled by transposing the
00164             //  matrix, factorizing it with the lapack function DPOTRF ('U', ...)
00165             //  and finally transposing the factor. It would be more efficient to use
00166             //  DPOTRF ('L', ...) in this case.
00167             LLt = true;
00168           else if (tmp.compare ("upper") == 0)
00169             LLt = false;
00170           else
00171             error ("chol: unexpected second or third input");
00172         }
00173       else
00174         error ("chol: expecting trailing string arguments");
00175     }
00176 
00177   if (! error_state)
00178     {
00179       octave_value arg = args(0);
00180 
00181       octave_idx_type nr = arg.rows ();
00182       octave_idx_type nc = arg.columns ();
00183       bool natural = (nargout != 3);
00184 
00185       int arg_is_empty = empty_arg ("chol", nr, nc);
00186 
00187       if (arg_is_empty < 0)
00188         return retval;
00189       if (arg_is_empty > 0)
00190         return octave_value (Matrix ());
00191 
00192       if (arg.is_sparse_type ())
00193         {
00194           if (arg.is_real_type ())
00195             {
00196               SparseMatrix m = arg.sparse_matrix_value ();
00197 
00198               if (! error_state)
00199                 {
00200                   octave_idx_type info;
00201                   SparseCHOL fact (m, info, natural);
00202                   if (nargout == 3)
00203                     {
00204                       if (vecout)
00205                         retval(2) = fact.perm ();
00206                       else
00207                         retval(2) = fact.Q();
00208                     }
00209 
00210                   if (nargout > 1 || info == 0)
00211                     {
00212                       retval(1) = fact.P();
00213                       if (LLt)
00214                         retval(0) = fact.L();
00215                       else
00216                         retval(0) = fact.R();
00217                     }
00218                   else
00219                     error ("chol: input matrix must be positive definite");
00220                 }
00221             }
00222           else if (arg.is_complex_type ())
00223             {
00224               SparseComplexMatrix m = arg.sparse_complex_matrix_value ();
00225 
00226               if (! error_state)
00227                 {
00228                   octave_idx_type info;
00229                   SparseComplexCHOL fact (m, info, natural);
00230 
00231                   if (nargout == 3)
00232                     {
00233                       if (vecout)
00234                         retval(2) = fact.perm ();
00235                       else
00236                         retval(2) = fact.Q();
00237                     }
00238 
00239                   if (nargout > 1 || info == 0)
00240                     {
00241                       retval(1) = fact.P();
00242                       if (LLt)
00243                         retval(0) = fact.L();
00244                       else
00245                         retval(0) = fact.R();
00246                     }
00247                   else
00248                     error ("chol: input matrix must be positive definite");
00249                 }
00250             }
00251           else
00252             gripe_wrong_type_arg ("chol", arg);
00253         }
00254       else if (arg.is_single_type ())
00255         {
00256           if (arg.is_real_type ())
00257             {
00258               FloatMatrix m = arg.float_matrix_value ();
00259 
00260               if (! error_state)
00261                 {
00262                   octave_idx_type info;
00263 
00264                   FloatCHOL fact;
00265                   if (LLt)
00266                     fact = FloatCHOL (m.transpose (), info);
00267                   else
00268                     fact = FloatCHOL (m, info);
00269 
00270                   if (nargout == 2 || info == 0)
00271                     {
00272                       retval(1) = info;
00273                       if (LLt)
00274                         retval(0) = get_chol_l (fact);
00275                       else
00276                         retval(0) = get_chol_r (fact);
00277                     }
00278                   else
00279                     error ("chol: input matrix must be positive definite");
00280                 }
00281             }
00282           else if (arg.is_complex_type ())
00283             {
00284               FloatComplexMatrix m = arg.float_complex_matrix_value ();
00285 
00286               if (! error_state)
00287                 {
00288                   octave_idx_type info;
00289 
00290                   FloatComplexCHOL fact;
00291                   if (LLt)
00292                     fact = FloatComplexCHOL (m.transpose (), info);
00293                   else
00294                     fact = FloatComplexCHOL (m, info);
00295 
00296                   if (nargout == 2 || info == 0)
00297                     {
00298                       retval(1) = info;
00299                       if (LLt)
00300                         retval(0) = get_chol_l (fact);
00301                       else
00302                         retval(0) = get_chol_r (fact);
00303                     }
00304                   else
00305                     error ("chol: input matrix must be positive definite");
00306                 }
00307             }
00308           else
00309             gripe_wrong_type_arg ("chol", arg);
00310         }
00311       else
00312         {
00313           if (arg.is_real_type ())
00314             {
00315               Matrix m = arg.matrix_value ();
00316 
00317               if (! error_state)
00318                 {
00319                   octave_idx_type info;
00320                   
00321                   CHOL fact;
00322                   if (LLt)
00323                      fact = CHOL (m.transpose (), info);
00324                   else
00325                     fact = CHOL (m, info);
00326 
00327                   if (nargout == 2 || info == 0)
00328                     {
00329                       retval(1) = info;
00330                       if (LLt)
00331                         retval(0) = get_chol_l (fact);
00332                       else
00333                         retval(0) = get_chol_r (fact);
00334                     }
00335                   else
00336                     error ("chol: input matrix must be positive definite");
00337                 }
00338             }
00339           else if (arg.is_complex_type ())
00340             {
00341               ComplexMatrix m = arg.complex_matrix_value ();
00342 
00343               if (! error_state)
00344                 {
00345                   octave_idx_type info;
00346                   
00347                   ComplexCHOL fact;
00348                   if (LLt)
00349                     fact = ComplexCHOL (m.transpose (), info);
00350                   else
00351                     fact = ComplexCHOL (m, info);
00352 
00353                   if (nargout == 2 || info == 0)
00354                     {
00355                       retval(1) = info;
00356                       if (LLt)
00357                         retval(0) = get_chol_l (fact);
00358                       else
00359                         retval(0) = get_chol_r (fact);
00360                     }
00361                   else
00362                     error ("chol: input matrix must be positive definite");
00363                 }
00364             }
00365           else
00366             gripe_wrong_type_arg ("chol", arg);
00367         }
00368     }
00369 
00370   return retval;
00371 }
00372 
00373 /*
00374 
00375 %!assert(chol ([2, 1; 1, 1]), [sqrt(2), 1/sqrt(2); 0, 1/sqrt(2)], sqrt (eps));
00376 %!assert(chol (single([2, 1; 1, 1])), single([sqrt(2), 1/sqrt(2); 0, 1/sqrt(2)]), sqrt (eps('single')));
00377 
00378 %!error chol ([1, 2; 3, 4]);
00379 %!error chol ([1, 2; 3, 4; 5, 6]);
00380 %!error <Invalid call to chol> chol ();
00381 %!error <unexpected second or third input> chol (1, 2);
00382 
00383  */
00384 
00385 DEFUN_DLD (cholinv, args, ,
00386   "-*- texinfo -*-\n\
00387 @deftypefn {Loadable Function} {} cholinv (@var{A})\n\
00388 Use the Cholesky@tie{}factorization to compute the inverse of the\n\
00389 symmetric positive definite matrix @var{A}.\n\
00390 @seealso{chol, chol2inv, inv}\n\
00391 @end deftypefn")
00392 {
00393   octave_value retval;
00394 
00395   int nargin = args.length ();
00396 
00397   if (nargin == 1)
00398     {
00399       octave_value arg = args(0);
00400 
00401       octave_idx_type nr = arg.rows ();
00402       octave_idx_type nc = arg.columns ();
00403 
00404       if (nr == 0 || nc == 0)
00405         retval = Matrix ();
00406       else
00407         {
00408           if (arg.is_sparse_type ())
00409             {
00410               if (arg.is_real_type ())
00411                 {
00412                   SparseMatrix m = arg.sparse_matrix_value ();
00413 
00414                   if (! error_state)
00415                     {
00416                       octave_idx_type info;
00417                       SparseCHOL chol (m, info);
00418                       if (info == 0)
00419                         retval = chol.inverse ();
00420                       else
00421                         error ("cholinv: A must be positive definite");
00422                     }
00423                 }
00424               else if (arg.is_complex_type ())
00425                 {
00426                   SparseComplexMatrix m = arg.sparse_complex_matrix_value ();
00427 
00428                   if (! error_state)
00429                     {
00430                       octave_idx_type info;
00431                       SparseComplexCHOL chol (m, info);
00432                       if (info == 0)
00433                         retval = chol.inverse ();
00434                       else
00435                         error ("cholinv: A must be positive definite");
00436                     }
00437                 }
00438               else
00439                 gripe_wrong_type_arg ("cholinv", arg);
00440             }
00441           else if (arg.is_single_type ())
00442             {
00443               if (arg.is_real_type ())
00444                 {
00445                   FloatMatrix m = arg.float_matrix_value ();
00446 
00447                   if (! error_state)
00448                     {
00449                       octave_idx_type info;
00450                       FloatCHOL chol (m, info);
00451                       if (info == 0)
00452                         retval = chol.inverse ();
00453                       else
00454                         error ("cholinv: A must be positive definite");
00455                     }
00456                 }
00457               else if (arg.is_complex_type ())
00458                 {
00459                   FloatComplexMatrix m = arg.float_complex_matrix_value ();
00460 
00461                   if (! error_state)
00462                     {
00463                       octave_idx_type info;
00464                       FloatComplexCHOL chol (m, info);
00465                       if (info == 0)
00466                         retval = chol.inverse ();
00467                       else
00468                         error ("cholinv: A must be positive definite");
00469                     }
00470                 }
00471               else
00472                 gripe_wrong_type_arg ("chol", arg);
00473             }
00474           else
00475             {
00476               if (arg.is_real_type ())
00477                 {
00478                   Matrix m = arg.matrix_value ();
00479 
00480                   if (! error_state)
00481                     {
00482                       octave_idx_type info;
00483                       CHOL chol (m, info);
00484                       if (info == 0)
00485                         retval = chol.inverse ();
00486                       else
00487                         error ("cholinv: A must be positive definite");
00488                     }
00489                 }
00490               else if (arg.is_complex_type ())
00491                 {
00492                   ComplexMatrix m = arg.complex_matrix_value ();
00493 
00494                   if (! error_state)
00495                     {
00496                       octave_idx_type info;
00497                       ComplexCHOL chol (m, info);
00498                       if (info == 0)
00499                         retval = chol.inverse ();
00500                       else
00501                         error ("cholinv: A must be positive definite");
00502                     }
00503                 }
00504               else
00505                 gripe_wrong_type_arg ("chol", arg);
00506             }
00507         }
00508     }
00509   else
00510     print_usage ();
00511 
00512   return retval;
00513 }
00514 
00515 /*
00516 
00517 %!shared A, Ainv
00518 %! A = [2,0.2;0.2,1];
00519 %! Ainv = inv(A);
00520 %!test
00521 %! Ainv1 = cholinv(A);
00522 %! assert (norm(Ainv-Ainv1),0,1e-10)
00523 %!testif HAVE_CHOLMOD
00524 %! Ainv2 = inv(sparse(A));
00525 %! assert (norm(Ainv-Ainv2),0,1e-10)
00526 %!testif HAVE_CHOLMOD
00527 %! Ainv3 = cholinv(sparse(A));
00528 %! assert (norm(Ainv-Ainv3),0,1e-10)
00529 
00530 */
00531 
00532 DEFUN_DLD (chol2inv, args, ,
00533   "-*- texinfo -*-\n\
00534 @deftypefn {Loadable Function} {} chol2inv (@var{U})\n\
00535 Invert a symmetric, positive definite square matrix from its Cholesky\n\
00536 decomposition, @var{U}.  Note that @var{U} should be an upper-triangular\n\
00537 matrix with positive diagonal elements.  @code{chol2inv (@var{U})}\n\
00538 provides @code{inv (@var{U}'*@var{U})} but it is much faster than\n\
00539 using @code{inv}.\n\
00540 @seealso{chol, cholinv, inv}\n\
00541 @end deftypefn")
00542 {
00543   octave_value retval;
00544 
00545   int nargin = args.length ();
00546 
00547   if (nargin == 1)
00548     {
00549       octave_value arg = args(0);
00550 
00551       octave_idx_type nr = arg.rows ();
00552       octave_idx_type nc = arg.columns ();
00553 
00554       if (nr == 0 || nc == 0)
00555         retval = Matrix ();
00556       else
00557         {
00558           if (arg.is_sparse_type ())
00559             {
00560               if (arg.is_real_type ())
00561                 {
00562                   SparseMatrix r = arg.sparse_matrix_value ();
00563 
00564                   if (! error_state)
00565                     retval = chol2inv (r);
00566                 }
00567               else if (arg.is_complex_type ())
00568                 {
00569                   SparseComplexMatrix r = arg.sparse_complex_matrix_value ();
00570 
00571                   if (! error_state)
00572                     retval = chol2inv (r);
00573                 }
00574               else
00575                 gripe_wrong_type_arg ("chol2inv", arg);
00576             }
00577           else if (arg.is_single_type ())
00578             {
00579               if (arg.is_real_type ())
00580                 {
00581                   FloatMatrix r = arg.float_matrix_value ();
00582 
00583                   if (! error_state)
00584                     retval = chol2inv (r);
00585                 }
00586               else if (arg.is_complex_type ())
00587                 {
00588                   FloatComplexMatrix r = arg.float_complex_matrix_value ();
00589 
00590                   if (! error_state)
00591                     retval = chol2inv (r);
00592                 }
00593               else
00594                 gripe_wrong_type_arg ("chol2inv", arg);
00595 
00596             }
00597           else
00598             {
00599               if (arg.is_real_type ())
00600                 {
00601                   Matrix r = arg.matrix_value ();
00602 
00603                   if (! error_state)
00604                     retval = chol2inv (r);
00605                 }
00606               else if (arg.is_complex_type ())
00607                 {
00608                   ComplexMatrix r = arg.complex_matrix_value ();
00609 
00610                   if (! error_state)
00611                     retval = chol2inv (r);
00612                 }
00613               else
00614                 gripe_wrong_type_arg ("chol2inv", arg);
00615             }
00616         }
00617     }
00618   else
00619     print_usage ();
00620 
00621   return retval;
00622 }
00623 
00624 DEFUN_DLD (cholupdate, args, nargout,
00625   "-*- texinfo -*-\n\
00626 @deftypefn {Loadable Function} {[@var{R1}, @var{info}] =} cholupdate (@var{R}, @var{u}, @var{op})\n\
00627 Update or downdate a Cholesky@tie{}factorization.  Given an upper triangular\n\
00628 matrix @var{R} and a column vector @var{u}, attempt to determine another\n\
00629 upper triangular matrix @var{R1} such that\n\
00630 @itemize @bullet\n\
00631 @item\n\
00632 @var{R1}'*@var{R1} = @var{R}'*@var{R} + @var{u}*@var{u}'\n\
00633 if @var{op} is \"+\"\n\
00634 \n\
00635 @item\n\
00636 @var{R1}'*@var{R1} = @var{R}'*@var{R} - @var{u}*@var{u}'\n\
00637 if @var{op} is \"-\"\n\
00638 @end itemize\n\
00639 \n\
00640 If @var{op} is \"-\", @var{info} is set to\n\
00641 @itemize\n\
00642 @item 0 if the downdate was successful,\n\
00643 \n\
00644 @item 1 if @var{R}'*@var{R} - @var{u}*@var{u}' is not positive definite,\n\
00645 \n\
00646 @item 2 if @var{R} is singular.\n\
00647 @end itemize\n\
00648 \n\
00649 If @var{info} is not present, an error message is printed in cases 1 and 2.\n\
00650 @seealso{chol, qrupdate}\n\
00651 @end deftypefn")
00652 {
00653   octave_idx_type nargin = args.length ();
00654 
00655   octave_value_list retval;
00656 
00657   if (nargin > 3 || nargin < 2)
00658     {
00659       print_usage ();
00660       return retval;
00661     }
00662 
00663   octave_value argr = args(0);
00664   octave_value argu = args(1);
00665 
00666   if (argr.is_numeric_type () && argu.is_numeric_type ()
00667       && (nargin < 3 || args(2).is_string ()))
00668     {
00669       octave_idx_type n = argr.rows ();
00670 
00671       std::string op = (nargin < 3) ? "+" : args(2).string_value ();
00672 
00673       bool down = op == "-";
00674 
00675       if (down || op == "+")
00676         if (argr.columns () == n && argu.rows () == n && argu.columns () == 1)
00677           {
00678             int err = 0;
00679             if (argr.is_single_type () || argu.is_single_type ())
00680               {
00681                 if (argr.is_real_type () && argu.is_real_type ())
00682                   {
00683                     // real case
00684                     FloatMatrix R = argr.float_matrix_value ();
00685                     FloatColumnVector u = argu.float_column_vector_value ();
00686 
00687                     FloatCHOL fact;
00688                     fact.set (R);
00689 
00690                     if (down)
00691                       err = fact.downdate (u);
00692                     else
00693                       fact.update (u);
00694 
00695                     retval(0) = get_chol_r (fact);
00696                   }
00697                 else
00698                   {
00699                     // complex case
00700                     FloatComplexMatrix R = argr.float_complex_matrix_value ();
00701                     FloatComplexColumnVector u = argu.float_complex_column_vector_value ();
00702 
00703                     FloatComplexCHOL fact;
00704                     fact.set (R);
00705 
00706                     if (down)
00707                       err = fact.downdate (u);
00708                     else
00709                       fact.update (u);
00710 
00711                     retval(0) = get_chol_r (fact);
00712                   }
00713               }
00714             else
00715               {
00716                 if (argr.is_real_type () && argu.is_real_type ())
00717                   {
00718                     // real case
00719                     Matrix R = argr.matrix_value ();
00720                     ColumnVector u = argu.column_vector_value ();
00721 
00722                     CHOL fact;
00723                     fact.set (R);
00724 
00725                     if (down)
00726                       err = fact.downdate (u);
00727                     else
00728                       fact.update (u);
00729 
00730                     retval(0) = get_chol_r (fact);
00731                   }
00732                 else
00733                   {
00734                     // complex case
00735                     ComplexMatrix R = argr.complex_matrix_value ();
00736                     ComplexColumnVector u = argu.complex_column_vector_value ();
00737 
00738                     ComplexCHOL fact;
00739                     fact.set (R);
00740 
00741                     if (down)
00742                       err = fact.downdate (u);
00743                     else
00744                       fact.update (u);
00745 
00746                     retval(0) = get_chol_r (fact);
00747                   }
00748               }
00749 
00750             if (nargout > 1)
00751               retval(1) = err;
00752             else if (err == 1)
00753               error ("cholupdate: downdate violates positiveness");
00754             else if (err == 2)
00755               error ("cholupdate: singular matrix");
00756           }
00757         else
00758           error ("cholupdate: dimension mismatch between R and U");
00759       else
00760         error ("cholupdate: OP must be \"+\" or \"-\"");
00761     }
00762   else
00763     print_usage ();
00764 
00765   return retval;
00766 }
00767 
00768 /*
00769 %!shared A, u, Ac, uc
00770 %! A = [  0.436997  -0.131721   0.124120  -0.061673 ;
00771 %!       -0.131721   0.738529   0.019851  -0.140295 ;
00772 %!        0.124120   0.019851   0.354879  -0.059472 ;
00773 %!       -0.061673  -0.140295  -0.059472   0.600939 ];
00774 %!
00775 %! u = [  0.98950 ;
00776 %!        0.39844 ;
00777 %!        0.63484 ;
00778 %!        0.13351 ];
00779 %! Ac = [  0.5585528 + 0.0000000i  -0.1662088 - 0.0315341i   0.0107873 + 0.0236411i  -0.0276775 - 0.0186073i ;
00780 %!        -0.1662088 + 0.0315341i   0.6760061 + 0.0000000i   0.0011452 - 0.0475528i   0.0145967 + 0.0247641i ;
00781 %!         0.0107873 - 0.0236411i   0.0011452 + 0.0475528i   0.6263149 - 0.0000000i  -0.1585837 - 0.0719763i ;
00782 %!        -0.0276775 + 0.0186073i   0.0145967 - 0.0247641i  -0.1585837 + 0.0719763i   0.6034234 - 0.0000000i ];
00783 %!
00784 %! uc = [ 0.54267 + 0.91519i ;
00785 %!        0.99647 + 0.43141i ;
00786 %!        0.83760 + 0.68977i ;
00787 %!        0.39160 + 0.90378i ];
00788 
00789 
00790 
00791 %!test
00792 %! R = chol(A);
00793 %!
00794 %! R1 = cholupdate(R,u);
00795 %!
00796 %! assert(norm(triu(R1)-R1,Inf) == 0)
00797 %! assert(norm(R1'*R1 - R'*R - u*u',Inf) < 1e1*eps)
00798 %!
00799 %! R1 = cholupdate(R1,u,"-");
00800 %!
00801 %! assert(norm(triu(R1)-R1,Inf) == 0)
00802 %! assert(norm(R1 - R,Inf) < 1e1*eps)
00803 %!
00804 %!test
00805 %! R = chol(Ac);
00806 %!
00807 %! R1 = cholupdate(R,uc);
00808 %!
00809 %! assert(norm(triu(R1)-R1,Inf) == 0)
00810 %! assert(norm(R1'*R1 - R'*R - uc*uc',Inf) < 1e1*eps)
00811 %!
00812 %! R1 = cholupdate(R1,uc,"-");
00813 %!
00814 %! assert(norm(triu(R1)-R1,Inf) == 0)
00815 %! assert(norm(R1 - R,Inf) < 1e1*eps)
00816 
00817 %!test
00818 %! R = chol(single(A));
00819 %!
00820 %! R1 = cholupdate(R,single(u));
00821 %!
00822 %! assert(norm(triu(R1)-R1,Inf) == 0)
00823 %! assert(norm(R1'*R1 - R'*R - single(u*u'),Inf) < 1e1*eps('single'))
00824 %!
00825 %! R1 = cholupdate(R1,single(u),"-");
00826 %!
00827 %! assert(norm(triu(R1)-R1,Inf) == 0)
00828 %! assert(norm(R1 - R,Inf) < 2e1*eps('single'))
00829 %!
00830 %!test
00831 %! R = chol(single(Ac));
00832 %!
00833 %! R1 = cholupdate(R,single(uc));
00834 %!
00835 %! assert(norm(triu(R1)-R1,Inf) == 0)
00836 %! assert(norm(R1'*R1 - R'*R - single(uc*uc'),Inf) < 1e1*eps('single'))
00837 %!
00838 %! R1 = cholupdate(R1,single(uc),"-");
00839 %!
00840 %! assert(norm(triu(R1)-R1,Inf) == 0)
00841 %! assert(norm(R1 - R,Inf) < 2e1*eps('single'))
00842 */
00843 
00844 DEFUN_DLD (cholinsert, args, nargout,
00845   "-*- texinfo -*-\n\
00846 @deftypefn  {Loadable Function} {@var{R1} =} cholinsert (@var{R}, @var{j}, @var{u})\n\
00847 @deftypefnx {Loadable Function} {[@var{R1}, @var{info}] =} cholinsert (@var{R}, @var{j}, @var{u})\n\
00848 Given a Cholesky@tie{}factorization of a real symmetric or complex Hermitian\n\
00849 positive definite matrix @w{@var{A} = @var{R}'*@var{R}}, @var{R}@tie{}upper\n\
00850 triangular, return the Cholesky@tie{}factorization of\n\
00851 @var{A1}, where @w{A1(p,p) = A}, @w{A1(:,j) = A1(j,:)' = u} and\n\
00852 @w{p = [1:j-1,j+1:n+1]}.  @w{u(j)} should be positive.\n\
00853 On return, @var{info} is set to\n\
00854 @itemize\n\
00855 @item 0 if the insertion was successful,\n\
00856 \n\
00857 @item 1 if @var{A1} is not positive definite,\n\
00858 \n\
00859 @item 2 if @var{R} is singular.\n\
00860 @end itemize\n\
00861 \n\
00862 If @var{info} is not present, an error message is printed in cases 1 and 2.\n\
00863 @seealso{chol, cholupdate, choldelete}\n\
00864 @end deftypefn")
00865 {
00866   octave_idx_type nargin = args.length ();
00867 
00868   octave_value_list retval;
00869 
00870   if (nargin != 3)
00871     {
00872       print_usage ();
00873       return retval;
00874     }
00875 
00876   octave_value argr = args(0);
00877   octave_value argj = args(1);
00878   octave_value argu = args(2);
00879 
00880   if (argr.is_numeric_type () && argu.is_numeric_type ()
00881       && argj.is_real_scalar ())
00882     {
00883       octave_idx_type n = argr.rows ();
00884       octave_idx_type j = argj.scalar_value ();
00885 
00886       if (argr.columns () == n && argu.rows () == n+1 && argu.columns () == 1)
00887         {
00888           if (j > 0 && j <= n+1)
00889             {
00890               int err = 0;
00891               if (argr.is_single_type () || argu.is_single_type ())
00892                 {
00893                   if (argr.is_real_type () && argu.is_real_type ())
00894                     {
00895                       // real case
00896                       FloatMatrix R = argr.float_matrix_value ();
00897                       FloatColumnVector u = argu.float_column_vector_value ();
00898 
00899                       FloatCHOL fact;
00900                       fact.set (R);
00901                       err = fact.insert_sym (u, j-1);
00902 
00903                       retval(0) = get_chol_r (fact);
00904                     }
00905                   else
00906                     {
00907                       // complex case
00908                       FloatComplexMatrix R = argr.float_complex_matrix_value ();
00909                       FloatComplexColumnVector u = argu.float_complex_column_vector_value ();
00910 
00911                       FloatComplexCHOL fact;
00912                       fact.set (R);
00913                       err = fact.insert_sym (u, j-1);
00914 
00915                       retval(0) = get_chol_r (fact);
00916                     }
00917                 }
00918               else
00919                 {
00920                   if (argr.is_real_type () && argu.is_real_type ())
00921                     {
00922                       // real case
00923                       Matrix R = argr.matrix_value ();
00924                       ColumnVector u = argu.column_vector_value ();
00925 
00926                       CHOL fact;
00927                       fact.set (R);
00928                       err = fact.insert_sym (u, j-1);
00929 
00930                       retval(0) = get_chol_r (fact);
00931                     }
00932                   else
00933                     {
00934                       // complex case
00935                       ComplexMatrix R = argr.complex_matrix_value ();
00936                       ComplexColumnVector u = argu.complex_column_vector_value ();
00937 
00938                       ComplexCHOL fact;
00939                       fact.set (R);
00940                       err = fact.insert_sym (u, j-1);
00941 
00942                       retval(0) = get_chol_r (fact);
00943                     }
00944                 }
00945 
00946               if (nargout > 1)
00947                 retval(1) = err;
00948               else if (err == 1)
00949                 error ("cholinsert: insertion violates positiveness");
00950               else if (err == 2)
00951                 error ("cholinsert: singular matrix");
00952               else if (err == 3)
00953                 error ("cholinsert: diagonal element must be real");
00954             }
00955           else
00956             error ("cholinsert: index J out of range");
00957         }
00958       else
00959         error ("cholinsert: dimension mismatch between R and U");
00960     }
00961   else
00962     print_usage ();
00963 
00964   return retval;
00965 }
00966 
00967 /*
00968 %!test
00969 %! u2 = [  0.35080 ;
00970 %!         0.63930 ;
00971 %!         3.31057 ;
00972 %!        -0.13825 ;
00973 %!         0.45266 ];
00974 %!
00975 %! R = chol(A);
00976 %!
00977 %! j = 3; p = [1:j-1, j+1:5];
00978 %! R1 = cholinsert(R,j,u2); A1 = R1'*R1;
00979 %!
00980 %! assert(norm(triu(R1)-R1,Inf) == 0)
00981 %! assert(norm(A1(p,p) - A,Inf) < 1e1*eps)
00982 %!
00983 %!test
00984 %! u2 = [  0.35080  + 0.04298i;
00985 %!         0.63930  + 0.23778i;
00986 %!         3.31057  + 0.00000i;
00987 %!        -0.13825  + 0.19879i;
00988 %!         0.45266  + 0.50020i];
00989 %!
00990 %! R = chol(Ac);
00991 %!
00992 %! j = 3; p = [1:j-1, j+1:5];
00993 %! R1 = cholinsert(R,j,u2); A1 = R1'*R1;
00994 %!
00995 %! assert(norm(triu(R1)-R1,Inf) == 0)
00996 %! assert(norm(A1(p,p) - Ac,Inf) < 1e1*eps)
00997 %!
00998 
00999 %!test
01000 %! u2 = single ([  0.35080 ;
01001 %!                 0.63930 ;
01002 %!                 3.31057 ;
01003 %!                -0.13825 ;
01004 %!                 0.45266 ]);
01005 %!
01006 %! R = chol(single(A));
01007 %!
01008 %! j = 3; p = [1:j-1, j+1:5];
01009 %! R1 = cholinsert(R,j,u2); A1 = R1'*R1;
01010 %!
01011 %! assert(norm(triu(R1)-R1,Inf) == 0)
01012 %! assert(norm(A1(p,p) - A,Inf) < 1e1*eps('single'))
01013 %!
01014 %!test
01015 %! u2 = single ([  0.35080  + 0.04298i;
01016 %!                 0.63930  + 0.23778i;
01017 %!                 3.31057  + 0.00000i;
01018 %!                -0.13825  + 0.19879i;
01019 %!                 0.45266  + 0.50020i]);
01020 %!
01021 %! R = chol(single(Ac));
01022 %!
01023 %! j = 3; p = [1:j-1, j+1:5];
01024 %! R1 = cholinsert(R,j,u2); A1 = R1'*R1;
01025 %!
01026 %! assert(norm(triu(R1)-R1,Inf) == 0)
01027 %! assert(norm(A1(p,p) - single(Ac),Inf) < 2e1*eps('single'))
01028 %!
01029 
01030 %!test
01031 %! cu = chol (triu (A), 'upper');
01032 %! cl = chol (tril (A), 'lower');
01033 %! assert (cu, cl', eps)
01034 %!
01035 %!test
01036 %! cca  = chol (Ac);
01037 %!
01038 %! ccal  = chol (Ac, 'lower');
01039 %! ccal2 = chol (tril (Ac), 'lower');
01040 %!
01041 %! ccau  = chol (Ac, 'upper');
01042 %! ccau2 = chol (triu (Ac), 'upper');
01043 %!
01044 %! assert (cca'*cca,     Ac, eps)
01045 %! assert (ccau'*ccau,   Ac, eps)
01046 %! assert (ccau2'*ccau2, Ac, eps)
01047 %!
01048 %! assert (cca, ccal',  eps)
01049 %! assert (cca, ccau,   eps)
01050 %! assert (cca, ccal2', eps)
01051 %! assert (cca, ccau2,  eps)
01052 %!
01053 %!test
01054 %! cca  = chol (single (Ac));
01055 %!
01056 %! ccal  = chol (single (Ac), 'lower');
01057 %! ccal2 = chol (tril (single (Ac)), 'lower');
01058 %!
01059 %! ccau  = chol (single (Ac), 'upper');
01060 %! ccau2 = chol (triu (single (Ac)), 'upper');
01061 %!
01062 %! assert (cca'*cca,     single (Ac), eps ('single'))
01063 %! assert (ccau'*ccau,   single (Ac), eps ('single'))
01064 %! assert (ccau2'*ccau2, single (Ac), eps ('single'))
01065 %!
01066 %! assert (cca, ccal',  eps ('single'))
01067 %! assert (cca, ccau,   eps ('single'))
01068 %! assert (cca, ccal2', eps ('single'))
01069 %! assert (cca, ccau2,  eps ('single'))
01070 
01071 %!test
01072 %! a = [12,  2,  3,  4;
01073 %!       2, 14,  5,  3;
01074 %!       3,  5, 16,  6;
01075 %!       4,  3,  6, 16];
01076 %!
01077 %! b = [0,  1,  2,  3;
01078 %!     -1,  0,  1,  2;
01079 %!     -2, -1,  0,  1;
01080 %!     -3, -2, -1,  0];
01081 %!
01082 %! ca = a + i*b;
01083 %!   
01084 %! cca  = chol (ca);
01085 %!
01086 %! ccal  = chol (ca, 'lower');
01087 %! ccal2 = chol (tril (ca), 'lower');
01088 %!
01089 %! ccau  = chol (ca, 'upper');
01090 %! ccau2 = chol (triu (ca), 'upper');
01091 %!
01092 %! assert (cca'*cca,     ca, 16*eps)
01093 %! assert (ccau'*ccau,   ca, 16*eps)
01094 %! assert (ccau2'*ccau2, ca, 16*eps)
01095 %!
01096 %! assert (cca, ccal',  16*eps)
01097 %! assert (cca, ccau,   16*eps)
01098 %! assert (cca, ccal2', 16*eps)
01099 %! assert (cca, ccau2,  16*eps)
01100 
01101 */
01102 
01103 DEFUN_DLD (choldelete, args, ,
01104   "-*- texinfo -*-\n\
01105 @deftypefn {Loadable Function} {@var{R1} =} choldelete (@var{R}, @var{j})\n\
01106 Given a Cholesky@tie{}factorization of a real symmetric or complex Hermitian\n\
01107 positive definite matrix @w{@var{A} = @var{R}'*@var{R}}, @var{R}@tie{}upper\n\
01108 triangular, return the Cholesky@tie{}factorization of @w{A(p,p)}, where\n\
01109 @w{p = [1:j-1,j+1:n+1]}.\n\
01110 @seealso{chol, cholupdate, cholinsert}\n\
01111 @end deftypefn")
01112 {
01113   octave_idx_type nargin = args.length ();
01114 
01115   octave_value_list retval;
01116 
01117   if (nargin != 2)
01118     {
01119       print_usage ();
01120       return retval;
01121     }
01122 
01123   octave_value argr = args(0);
01124   octave_value argj = args(1);
01125 
01126   if (argr.is_numeric_type () && argj.is_real_scalar ())
01127     {
01128       octave_idx_type n = argr.rows ();
01129       octave_idx_type j = argj.scalar_value ();
01130 
01131       if (argr.columns () == n)
01132         {
01133           if (j > 0 && j <= n)
01134             {
01135               if (argr.is_single_type ())
01136                 {
01137                   if (argr.is_real_type ())
01138                     {
01139                       // real case
01140                       FloatMatrix R = argr.float_matrix_value ();
01141 
01142                       FloatCHOL fact;
01143                       fact.set (R);
01144                       fact.delete_sym (j-1);
01145 
01146                       retval(0) = get_chol_r (fact);
01147                     }
01148                   else
01149                     {
01150                       // complex case
01151                       FloatComplexMatrix R = argr.float_complex_matrix_value ();
01152 
01153                       FloatComplexCHOL fact;
01154                       fact.set (R);
01155                       fact.delete_sym (j-1);
01156 
01157                       retval(0) = get_chol_r (fact);
01158                     }
01159                 }
01160               else
01161                 {
01162                   if (argr.is_real_type ())
01163                     {
01164                       // real case
01165                       Matrix R = argr.matrix_value ();
01166 
01167                       CHOL fact;
01168                       fact.set (R);
01169                       fact.delete_sym (j-1);
01170 
01171                       retval(0) = get_chol_r (fact);
01172                     }
01173                   else
01174                     {
01175                       // complex case
01176                       ComplexMatrix R = argr.complex_matrix_value ();
01177 
01178                       ComplexCHOL fact;
01179                       fact.set (R);
01180                       fact.delete_sym (j-1);
01181 
01182                       retval(0) = get_chol_r (fact);
01183                     }
01184                 }
01185             }
01186           else
01187             error ("choldelete: index J out of range");
01188         }
01189       else
01190         error ("choldelete: matrix R must be square");
01191     }
01192   else
01193     print_usage ();
01194 
01195   return retval;
01196 }
01197 
01198 /*
01199 %!test
01200 %! R = chol(A);
01201 %!
01202 %! j = 3; p = [1:j-1,j+1:4];
01203 %! R1 = choldelete(R,j);
01204 %!
01205 %! assert(norm(triu(R1)-R1,Inf) == 0)
01206 %! assert(norm(R1'*R1 - A(p,p),Inf) < 1e1*eps)
01207 %!
01208 %!test
01209 %! R = chol(Ac);
01210 %!
01211 %! j = 3; p = [1:j-1,j+1:4];
01212 %! R1 = choldelete(R,j);
01213 %!
01214 %! assert(norm(triu(R1)-R1,Inf) == 0)
01215 %! assert(norm(R1'*R1 - Ac(p,p),Inf) < 1e1*eps)
01216 
01217 %!test
01218 %! R = chol(single(A));
01219 %!
01220 %! j = 3; p = [1:j-1,j+1:4];
01221 %! R1 = choldelete(R,j);
01222 %!
01223 %! assert(norm(triu(R1)-R1,Inf) == 0)
01224 %! assert(norm(R1'*R1 - single(A(p,p)),Inf) < 1e1*eps('single'))
01225 %!
01226 %!test
01227 %! R = chol(single(Ac));
01228 %!
01229 %! j = 3; p = [1:j-1,j+1:4];
01230 %! R1 = choldelete(R,j);
01231 %!
01232 %! assert(norm(triu(R1)-R1,Inf) == 0)
01233 %! assert(norm(R1'*R1 - single(Ac(p,p)),Inf) < 1e1*eps('single'))
01234 */
01235 
01236 DEFUN_DLD (cholshift, args, ,
01237   "-*- texinfo -*-\n\
01238 @deftypefn {Loadable Function} {@var{R1} =} cholshift (@var{R}, @var{i}, @var{j})\n\
01239 Given a Cholesky@tie{}factorization of a real symmetric or complex Hermitian\n\
01240 positive definite matrix @w{@var{A} = @var{R}'*@var{R}}, @var{R}@tie{}upper\n\
01241 triangular, return the Cholesky@tie{}factorization of\n\
01242 @w{@var{A}(p,p)}, where @w{p} is the permutation @*\n\
01243 @code{p = [1:i-1, shift(i:j, 1), j+1:n]} if @w{@var{i} < @var{j}} @*\n\
01244  or @*\n\
01245 @code{p = [1:j-1, shift(j:i,-1), i+1:n]} if @w{@var{j} < @var{i}}.  @*\n\
01246 \n\
01247 @seealso{chol, cholinsert, choldelete}\n\
01248 @end deftypefn")
01249 {
01250   octave_idx_type nargin = args.length ();
01251 
01252   octave_value_list retval;
01253 
01254   if (nargin != 3)
01255     {
01256       print_usage ();
01257       return retval;
01258     }
01259 
01260   octave_value argr = args(0);
01261   octave_value argi = args(1);
01262   octave_value argj = args(2);
01263 
01264   if (argr.is_numeric_type () && argi.is_real_scalar () && argj.is_real_scalar ())
01265     {
01266       octave_idx_type n = argr.rows ();
01267       octave_idx_type i = argi.scalar_value ();
01268       octave_idx_type j = argj.scalar_value ();
01269 
01270       if (argr.columns () == n)
01271         {
01272           if (j > 0 && j <= n+1 && i > 0 && i <= n+1)
01273             {
01274 
01275               if (argr.is_single_type () && argi.is_single_type () &&
01276                   argj.is_single_type ())
01277                 {
01278                   if (argr.is_real_type ())
01279                     {
01280                       // real case
01281                       FloatMatrix R = argr.float_matrix_value ();
01282 
01283                       FloatCHOL fact;
01284                       fact.set (R);
01285                       fact.shift_sym (i-1, j-1);
01286 
01287                       retval(0) = get_chol_r (fact);
01288                     }
01289                   else
01290                     {
01291                       // complex case
01292                       FloatComplexMatrix R = argr.float_complex_matrix_value ();
01293 
01294                       FloatComplexCHOL fact;
01295                       fact.set (R);
01296                       fact.shift_sym (i-1, j-1);
01297 
01298                       retval(0) = get_chol_r (fact);
01299                     }
01300                 }
01301               else
01302                 {
01303                   if (argr.is_real_type ())
01304                     {
01305                       // real case
01306                       Matrix R = argr.matrix_value ();
01307 
01308                       CHOL fact;
01309                       fact.set (R);
01310                       fact.shift_sym (i-1, j-1);
01311 
01312                       retval(0) = get_chol_r (fact);
01313                     }
01314                   else
01315                     {
01316                       // complex case
01317                       ComplexMatrix R = argr.complex_matrix_value ();
01318 
01319                       ComplexCHOL fact;
01320                       fact.set (R);
01321                       fact.shift_sym (i-1, j-1);
01322 
01323                       retval(0) = get_chol_r (fact);
01324                     }
01325                 }
01326             }
01327           else
01328             error ("cholshift: index I or J is out of range");
01329         }
01330       else
01331         error ("cholshift: R must be a square matrix");
01332     }
01333   else
01334     print_usage ();
01335 
01336   return retval;
01337 }
01338 
01339 /*
01340 %!test
01341 %! R = chol(A);
01342 %!
01343 %! i = 1; j = 3; p = [1:i-1, shift(i:j,-1), j+1:4];
01344 %! R1 = cholshift(R,i,j);
01345 %!
01346 %! assert(norm(triu(R1)-R1,Inf) == 0)
01347 %! assert(norm(R1'*R1 - A(p,p),Inf) < 1e1*eps)
01348 %!
01349 %! j = 1; i = 3; p = [1:j-1, shift(j:i,+1), i+1:4];
01350 %! R1 = cholshift(R,i,j);
01351 %!
01352 %! assert(norm(triu(R1)-R1,Inf) == 0)
01353 %! assert(norm(R1'*R1 - A(p,p),Inf) < 1e1*eps)
01354 %!
01355 %!test
01356 %! R = chol(Ac);
01357 %!
01358 %! i = 1; j = 3; p = [1:i-1, shift(i:j,-1), j+1:4];
01359 %! R1 = cholshift(R,i,j);
01360 %!
01361 %! assert(norm(triu(R1)-R1,Inf) == 0)
01362 %! assert(norm(R1'*R1 - Ac(p,p),Inf) < 1e1*eps)
01363 %!
01364 %! j = 1; i = 3; p = [1:j-1, shift(j:i,+1), i+1:4];
01365 %! R1 = cholshift(R,i,j);
01366 %!
01367 %! assert(norm(triu(R1)-R1,Inf) == 0)
01368 %! assert(norm(R1'*R1 - Ac(p,p),Inf) < 1e1*eps)
01369 
01370 %!test
01371 %! R = chol(single(A));
01372 %!
01373 %! i = 1; j = 3; p = [1:i-1, shift(i:j,-1), j+1:4];
01374 %! R1 = cholshift(R,i,j);
01375 %!
01376 %! assert(norm(triu(R1)-R1,Inf) == 0)
01377 %! assert(norm(R1'*R1 - single(A(p,p)),Inf) < 1e1*eps('single'))
01378 %!
01379 %! j = 1; i = 3; p = [1:j-1, shift(j:i,+1), i+1:4];
01380 %! R1 = cholshift(R,i,j);
01381 %!
01382 %! assert(norm(triu(R1)-R1,Inf) == 0)
01383 %! assert(norm(R1'*R1 - single(A(p,p)),Inf) < 1e1*eps('single'))
01384 %!
01385 %!test
01386 %! R = chol(single(Ac));
01387 %!
01388 %! i = 1; j = 3; p = [1:i-1, shift(i:j,-1), j+1:4];
01389 %! R1 = cholshift(R,i,j);
01390 %!
01391 %! assert(norm(triu(R1)-R1,Inf) == 0)
01392 %! assert(norm(R1'*R1 - single(Ac(p,p)),Inf) < 1e1*eps('single'))
01393 %!
01394 %! j = 1; i = 3; p = [1:j-1, shift(j:i,+1), i+1:4];
01395 %! R1 = cholshift(R,i,j);
01396 %!
01397 %! assert(norm(triu(R1)-R1,Inf) == 0)
01398 %! assert(norm(R1'*R1 - single(Ac(p,p)),Inf) < 1e1*eps('single'))
01399 */
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines