Mercurial > hg > octave-lyh
annotate liboctave/dbleCHOL.cc @ 7482:29980c6b8604
don't check f77_exception_encountered
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Thu, 14 Feb 2008 21:57:50 -0500 |
| parents | a1dbe9d80eee |
| children | 40574114c514 |
| rev | line source |
|---|---|
| 457 | 1 /* |
| 2 | |
| 7017 | 3 Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 2007 |
| 4 John W. Eaton | |
| 457 | 5 |
| 6 This file is part of Octave. | |
| 7 | |
| 8 Octave is free software; you can redistribute it and/or modify it | |
| 9 under the terms of the GNU General Public License as published by the | |
| 7016 | 10 Free Software Foundation; either version 3 of the License, or (at your |
| 11 option) any later version. | |
| 457 | 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 | |
| 7016 | 19 along with Octave; see the file COPYING. If not, see |
| 20 <http://www.gnu.org/licenses/>. | |
| 457 | 21 |
| 22 */ | |
| 23 | |
| 24 #ifdef HAVE_CONFIG_H | |
| 1192 | 25 #include <config.h> |
| 457 | 26 #endif |
| 27 | |
| 6486 | 28 #include "dRowVector.h" |
| 457 | 29 #include "dbleCHOL.h" |
| 1847 | 30 #include "f77-fcn.h" |
| 457 | 31 #include "lo-error.h" |
| 32 | |
| 33 extern "C" | |
| 34 { | |
| 4552 | 35 F77_RET_T |
| 5275 | 36 F77_FUNC (dpotrf, DPOTRF) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, |
| 37 double*, const octave_idx_type&, octave_idx_type& | |
| 4552 | 38 F77_CHAR_ARG_LEN_DECL); |
| 5340 | 39 |
| 40 F77_RET_T | |
| 41 F77_FUNC (dpotri, DPOTRI) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, | |
| 42 double*, const octave_idx_type&, octave_idx_type& | |
| 43 F77_CHAR_ARG_LEN_DECL); | |
| 6486 | 44 |
| 45 F77_RET_T | |
| 46 F77_FUNC (dpocon, DPOCON) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, | |
| 47 double*, const octave_idx_type&, const double&, | |
| 48 double&, double*, octave_idx_type*, | |
| 49 octave_idx_type& F77_CHAR_ARG_LEN_DECL); | |
| 457 | 50 } |
| 51 | |
| 5275 | 52 octave_idx_type |
| 6486 | 53 CHOL::init (const Matrix& a, bool calc_cond) |
| 457 | 54 { |
| 5275 | 55 octave_idx_type a_nr = a.rows (); |
| 56 octave_idx_type a_nc = a.cols (); | |
| 1944 | 57 |
| 457 | 58 if (a_nr != a_nc) |
| 59 { | |
| 60 (*current_liboctave_error_handler) ("CHOL requires square matrix"); | |
| 61 return -1; | |
| 62 } | |
| 63 | |
| 5275 | 64 octave_idx_type n = a_nc; |
| 65 octave_idx_type info; | |
| 457 | 66 |
| 1944 | 67 chol_mat = a; |
| 68 double *h = chol_mat.fortran_vec (); | |
| 457 | 69 |
| 6486 | 70 // Calculate the norm of the matrix, for later use. |
| 71 double anorm = 0; | |
| 72 if (calc_cond) | |
| 73 anorm = chol_mat.abs().sum().row(static_cast<octave_idx_type>(0)).max(); | |
| 74 | |
| 4552 | 75 F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1), |
| 76 n, h, n, info | |
| 77 F77_CHAR_ARG_LEN (1))); | |
| 457 | 78 |
|
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
79 xrcond = 0.0; |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
80 if (info != 0) |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
81 info = -1; |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
82 else if (calc_cond) |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
83 { |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
84 octave_idx_type dpocon_info = 0; |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
85 |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
86 // Now calculate the condition number for non-singular matrix. |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
87 Array<double> z (3*n); |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
88 double *pz = z.fortran_vec (); |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
89 Array<octave_idx_type> iz (n); |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
90 octave_idx_type *piz = iz.fortran_vec (); |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
91 F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h, |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
92 n, anorm, xrcond, pz, piz, dpocon_info |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
93 F77_CHAR_ARG_LEN (1))); |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
94 |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
95 if (dpocon_info != 0) |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
96 info = -1; |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
97 } |
| 1944 | 98 else |
| 99 { | |
|
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
100 // If someone thinks of a more graceful way of doing this (or |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
101 // faster for that matter :-)), please let me know! |
| 457 | 102 |
|
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
103 if (n > 1) |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
104 for (octave_idx_type j = 0; j < a_nc; j++) |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
105 for (octave_idx_type i = j+1; i < a_nr; i++) |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
106 chol_mat.xelem (i, j) = 0.0; |
| 1944 | 107 } |
| 457 | 108 |
| 109 return info; | |
| 110 } | |
| 111 | |
| 5340 | 112 static Matrix |
| 113 chol2inv_internal (const Matrix& r) | |
| 114 { | |
| 115 Matrix retval; | |
| 116 | |
| 117 octave_idx_type r_nr = r.rows (); | |
| 118 octave_idx_type r_nc = r.cols (); | |
| 119 | |
| 120 if (r_nr == r_nc) | |
| 121 { | |
| 122 octave_idx_type n = r_nc; | |
| 6486 | 123 octave_idx_type info = 0; |
| 5340 | 124 |
| 125 Matrix tmp = r; | |
| 6486 | 126 double *v = tmp.fortran_vec(); |
| 5340 | 127 |
| 6486 | 128 if (info == 0) |
| 5340 | 129 { |
| 6486 | 130 F77_XFCN (dpotri, DPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n, |
| 131 v, n, info | |
| 132 F77_CHAR_ARG_LEN (1))); | |
| 5340 | 133 |
|
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
134 // If someone thinks of a more graceful way of doing this (or |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
135 // faster for that matter :-)), please let me know! |
| 5340 | 136 |
|
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
137 if (n > 1) |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
138 for (octave_idx_type j = 0; j < r_nc; j++) |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
139 for (octave_idx_type i = j+1; i < r_nr; i++) |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
140 tmp.xelem (i, j) = tmp.xelem (j, i); |
| 6486 | 141 |
|
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
142 retval = tmp; |
| 5340 | 143 } |
| 144 } | |
| 145 else | |
| 146 (*current_liboctave_error_handler) ("chol2inv requires square matrix"); | |
| 147 | |
| 148 return retval; | |
| 149 } | |
| 150 | |
| 151 // Compute the inverse of a matrix using the Cholesky factorization. | |
| 152 Matrix | |
| 153 CHOL::inverse (void) const | |
| 154 { | |
| 155 return chol2inv_internal (chol_mat); | |
| 156 } | |
| 157 | |
| 158 Matrix | |
| 159 chol2inv (const Matrix& r) | |
| 160 { | |
| 161 return chol2inv_internal (r); | |
| 162 } | |
| 163 | |
| 457 | 164 /* |
| 165 ;;; Local Variables: *** | |
| 166 ;;; mode: C++ *** | |
| 167 ;;; End: *** | |
| 168 */ |