Mercurial > octave-nkf
annotate liboctave/CmplxCHOL.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 "dMatrix.h" |
| 29 #include "dRowVector.h" | |
| 457 | 30 #include "CmplxCHOL.h" |
| 1847 | 31 #include "f77-fcn.h" |
| 457 | 32 #include "lo-error.h" |
| 33 | |
| 34 extern "C" | |
| 35 { | |
| 4552 | 36 F77_RET_T |
| 5340 | 37 F77_FUNC (zpotrf, ZPOTRF) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, |
| 38 Complex*, const octave_idx_type&, octave_idx_type& | |
| 39 F77_CHAR_ARG_LEN_DECL); | |
| 40 F77_RET_T | |
| 41 F77_FUNC (zpotri, ZPOTRI) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, | |
| 42 Complex*, const octave_idx_type&, octave_idx_type& | |
| 4552 | 43 F77_CHAR_ARG_LEN_DECL); |
| 6486 | 44 |
| 45 F77_RET_T | |
| 46 F77_FUNC (zpocon, ZPOCON) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, | |
| 47 Complex*, const octave_idx_type&, const double&, | |
| 48 double&, Complex*, double*, | |
| 49 octave_idx_type& F77_CHAR_ARG_LEN_DECL); | |
| 457 | 50 } |
| 51 | |
| 5275 | 52 octave_idx_type |
| 6486 | 53 ComplexCHOL::init (const ComplexMatrix& a, bool calc_cond) |
| 457 | 54 { |
| 5275 | 55 octave_idx_type a_nr = a.rows (); |
| 56 octave_idx_type a_nc = a.cols (); | |
| 457 | 57 |
| 1944 | 58 if (a_nr != a_nc) |
| 59 { | |
| 60 (*current_liboctave_error_handler) | |
| 61 ("ComplexCHOL requires square matrix"); | |
| 62 return -1; | |
| 63 } | |
| 457 | 64 |
| 5275 | 65 octave_idx_type n = a_nc; |
| 66 octave_idx_type info; | |
| 1944 | 67 |
| 68 chol_mat = a; | |
| 69 Complex *h = chol_mat.fortran_vec (); | |
| 457 | 70 |
| 6486 | 71 // Calculate the norm of the matrix, for later use. |
| 72 double anorm = 0; | |
| 73 if (calc_cond) | |
| 74 anorm = chol_mat.abs().sum().row(static_cast<octave_idx_type>(0)).max(); | |
| 75 | |
| 4552 | 76 F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1), 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 zpocon_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<Complex> z (2*n); |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
88 Complex *pz = z.fortran_vec (); |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
89 Array<double> rz (n); |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
90 double *prz = rz.fortran_vec (); |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
91 F77_XFCN (zpocon, ZPOCON, (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, prz, zpocon_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 (zpocon_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 } |
| 108 | |
| 109 return info; | |
| 457 | 110 } |
| 111 | |
| 5340 | 112 static ComplexMatrix |
| 113 chol2inv_internal (const ComplexMatrix& r) | |
| 114 { | |
| 115 ComplexMatrix 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; | |
| 123 octave_idx_type info; | |
| 124 | |
| 125 ComplexMatrix tmp = r; | |
| 126 | |
| 127 F77_XFCN (zpotri, ZPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n, | |
| 128 tmp.fortran_vec (), n, info | |
| 129 F77_CHAR_ARG_LEN (1))); | |
| 130 | |
|
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
131 // 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
|
132 // faster for that matter :-)), please let me know! |
| 5340 | 133 |
|
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
134 if (n > 1) |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
135 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
|
136 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
|
137 tmp.xelem (i, j) = std::conj (tmp.xelem (j, i)); |
| 5340 | 138 |
|
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
139 retval = tmp; |
| 5340 | 140 } |
| 141 else | |
| 142 (*current_liboctave_error_handler) ("chol2inv requires square matrix"); | |
| 143 | |
| 144 return retval; | |
| 145 } | |
| 146 | |
| 147 // Compute the inverse of a matrix using the Cholesky factorization. | |
| 148 ComplexMatrix | |
| 149 ComplexCHOL::inverse (void) const | |
| 150 { | |
| 151 return chol2inv_internal (chol_mat); | |
| 152 } | |
| 153 | |
| 154 ComplexMatrix | |
| 155 chol2inv (const ComplexMatrix& r) | |
| 156 { | |
| 157 return chol2inv_internal (r); | |
| 158 } | |
| 159 | |
| 457 | 160 /* |
| 161 ;;; Local Variables: *** | |
| 162 ;;; mode: C++ *** | |
| 163 ;;; End: *** | |
| 164 */ |