GNU Octave  3.8.0 A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
SparseCmplxCHOL.cc
Go to the documentation of this file.
1 /*
2
3 Copyright (C) 2005-2013 David Bateman
5
6 This file is part of Octave.
7
8 Octave is free software; you can redistribute it and/or modify it
10 Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12
13 Octave is distributed in the hope that it will be useful, but WITHOUT
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
16 for more details.
17
18 You should have received a copy of the GNU General Public License
19 along with Octave; see the file COPYING. If not, see
21
22 */
23
24 #ifdef HAVE_CONFIG_H
25 #include <config.h>
26 #endif
27
28 #include "SparseCmplxCHOL.h"
29
30 // Instantiate the base CHOL class for the type we need
31 #define OCTAVE_CHOLMOD_TYPE CHOLMOD_COMPLEX
32 #include "sparse-base-chol.h"
33 #include "sparse-base-chol.cc"
35
36 // Compute the inverse of a matrix using the Cholesky factorization.
39 {
40  octave_idx_type r_nr = r.rows ();
41  octave_idx_type r_nc = r.cols ();
42  SparseComplexMatrix retval;
43
44  if (r_nr == r_nc)
45  {
46  MatrixType mattype (r);
47  int typ = mattype.type (false);
48  double rcond;
49  octave_idx_type info;
51
52  if (typ == MatrixType::Upper)
53  {
54  rinv = r.inverse (mattype, info, rcond, true, false);
55  retval = rinv.transpose () * rinv;
56  }
57  else if (typ == MatrixType::Lower)
58  {
59  rinv = r.transpose ().inverse (mattype, info, rcond, true, false);
60  retval = rinv.transpose () * rinv;
61  }
62  else
64  ("U must be a triangular matrix");
65  }
66  else
67  (*current_liboctave_error_handler) ("U must be a square matrix");
68
69  return retval;
70 }